[インデックス 13461] ファイルの概要

このコミットは、Go言語の math/big パッケージにおける多倍長整数乗算の実装に関する重要な修正を含んでいます。具体的には、Mul 関数（乗算）がアンバランスな入力（桁数の大きく異なる2つの数）に対して二次的な空間計算量（O(n^2)）と線形的な再帰深度（O(n)）を持つ問題を修正し、より効率的なメモリ使用と再帰深度を実現しています。この修正は、特に大きな数値の乗算においてパフォーマンスと安定性を向上させます。

コミット

commit ac12131649391c6303f96514aee4424cb7a0b7d7
Author: Rémy Oudompheng <oudomphe@phare.normalesup.org>
Date:   Thu Jul 12 10:18:24 2012 -0700

    math/big: correct quadratic space complexity in Mul.
    
    The previous implementation used to have a O(n) recursion
    depth for unbalanced inputs. A test is added to check that a
    reasonable amount of bytes is allocated in this case.
    
    Fixes #3807.
    
    R=golang-dev, dsymonds, gri
    CC=golang-dev, remy
    https://golang.org/cl/6345075

GitHub上でのコミットページへのリンク

https://github.com/golang/go/commit/ac12131649391c6303f96514aee4424cb7a0b7d7

元コミット内容

math/big: correct quadratic space complexity in Mul.

以前の実装では、アンバランスな入力に対してO(n)の再帰深度を持っていました。このケースで妥当な量のバイトが割り当てられていることを確認するためのテストが追加されました。

Fixes #3807.

変更の背景

この変更の背景には、Go言語の math/big パッケージにおける多倍長整数乗算の効率性の問題がありました。特に、Karatsuba乗算アルゴリズムの実装において、乗算される2つの数値の桁数に大きな差がある（アンバランスな入力）場合に、非効率なメモリ割り当てと過度な再帰深度が発生していました。

具体的には、以下の問題が指摘されていました。

二次的な空間計算量 (Quadratic Space Complexity): アンバランスな入力に対して、乗算処理中に必要とされる一時的なメモリ領域が入力サイズの二乗に比例して増加してしまう問題。これは、非常に大きな数値を扱う際にメモリ不足やパフォーマンスの著しい低下を引き起こす可能性があります。
線形的な再帰深度 (O(n) Recursion Depth): Karatsubaアルゴリズムは再帰的に問題を分割して解決しますが、アンバランスな入力の場合、再帰の深さが入力の桁数に比例して深くなりすぎていました。これにより、スタックオーバーフローのリスクや、再帰呼び出しのオーバーヘッドによるパフォーマンス劣化が発生していました。

これらの問題は、Go Issue #3807として報告されており、このコミットはその問題を解決するために行われました。目的は、Mul 関数がより効率的にメモリを使用し、アンバランスな入力に対しても安定したパフォーマンスを発揮するようにすることです。

前提知識の解説

1. 多倍長整数 (Arbitrary-Precision Integers)

通常のプログラミング言語で扱える整数型（例: int32, int64）には、表現できる数値の範囲に上限があります。しかし、暗号学、科学計算、金融アプリケーションなど、非常に大きな整数を扱う必要がある場面があります。多倍長整数は、このような上限を持たず、メモリが許す限り任意の桁数の整数を表現・計算できるデータ型です。Go言語の math/big パッケージは、この多倍長整数（Int型）や自然数（nat型）を提供します。

2. Karatsuba乗算アルゴリズム

Karatsubaアルゴリズムは、大きな数の乗算を効率的に行うためのアルゴリズムです。通常の筆算による乗算（O(n^2)）よりも高速で、O(n^log2(3)) ≈ O(n^1.585) の計算量で実行されます。

基本的な考え方は以下の通りです。 2つのn桁の数 x と y を考えます。 x = x1 * B + x0 y = y1 * B + y0 ここで B は基数（通常は 10^(n/2) など）、x0, y0 は下位 n/2 桁、x1, y1 は上位 n/2 桁を表します。

通常の乗算 x * y は以下のように展開されます。 x * y = (x1 * B + x0) * (y1 * B + y0) = x1*y1*B^2 + (x1*y0 + x0*y1)*B + x0*y0

この式には4回の乗算（x1*y1, x1*y0, x0*y1, x0*y0）が含まれます。Karatsubaアルゴリズムは、以下の3回の乗算で同じ結果を得ます。 P0 = x0 * y0 P1 = x1 * y1 P2 = (x0 + x1) * (y0 + y1)

そして、 x1*y0 + x0*y1 = P2 - P0 - P1

したがって、 x * y = P1*B^2 + (P2 - P0 - P1)*B + P0

このアルゴリズムは再帰的に適用され、乗算の計算量を削減します。しかし、再帰の深さや、分割された部分のサイズが大きく異なる（アンバランスな入力）場合に、効率が低下したり、メモリ使用量が増加したりする問題が発生することがあります。

3. 空間計算量 (Space Complexity) と再帰深度 (Recursion Depth)

空間計算量: アルゴリズムが実行中に必要とするメモリの総量です。O(n^2) は入力サイズ n の二乗に比例してメモリが増加することを意味し、非常に非効率です。
再帰深度: 再帰関数が自身を呼び出す最大の回数です。再帰呼び出しごとにスタックフレームが積まれるため、再帰深度が深すぎるとスタックオーバーフローを引き起こす可能性があります。O(n) は入力サイズ n に比例して再帰が深くなることを意味します。

4. `nat` 型 (Natural Number)

Go言語の math/big パッケージでは、nat 型は符号なしの多倍長整数を表す内部的なスライス型です。これは、Int 型の内部表現として使用されます。_W は nat の各要素が保持するビット幅（通常は uint のサイズ、32ビットまたは64ビット）を表します。

技術的詳細

このコミットは、math/big パッケージの nat.go ファイルにある (z nat) mul(x, y nat) nat 関数、すなわち多倍長自然数の乗算処理を修正しています。主な変更点は、Karatsuba乗算アルゴリズムにおけるアンバランスな入力の処理方法です。

以前の実装では、x と y の長さが大きく異なる場合（例えば x が非常に長く、y が非常に短い場合）、Karatsuba分割の際に x1 や y1 の一部がゼロになるにもかかわらず、それらの項を明示的に追加する処理が非効率でした。特に、x1*y0 や x0*y1 のような中間項の計算と加算において、不必要な再帰呼び出しやメモリ割り当てが発生し、結果として二次的な空間計算量と線形的な再帰深度につながっていました。

新しい実装では、このアンバランスなケースをより効率的に処理するために、以下の変更が加えられています。

x1*y0 および x0*y1 項の加算ロジックの改善:
- 以前は t.mul(x1, y0.norm()) や t.mul(x0.norm(), y1) のように、x1 や y1 が短い場合でもKaratsubaの再帰呼び出しを伴う mul を使用していました。
- 新しいコードでは、y1 が短いことを利用し、x0*y1*b の項を addAt(z, t.mul(x0, y1), k) で直接加算しています。ここで t.mul(x0, y1) は、y1 が短い場合はKaratsubaではなく、より効率的な古典的な乗算（mulAddVWW など）にフォールバックする可能性があります。
- さらに重要なのは、xi*y0<<i および xi*y1*b<<(i+k) の項をループで処理するようになった点です。これは、x が k の倍数で分割されることを利用し、x の残りの部分 (x[k:]) を k 単位で繰り返し処理することで、アンバランスな乗算をより細かく、かつ効率的に分解しています。これにより、不必要な大きな一時配列の生成や深い再帰が回避されます。
メモリ割り当ての最適化:
- z = z.make(max(6*k, m+n)) で必要なメモリを確保した後、z[2*k:].clear() を呼び出して、Karatsubaの計算で使われない上位部分をクリアしています。これにより、以前のゴミデータが残ることを防ぎ、メモリの再利用を促進します。
- 全体として、中間結果の生成と加算のロジックが改善されたことで、一時的に必要となるメモリの量が大幅に削減され、二次的な空間計算量の問題が解消されています。
テストケースの追加:
- nat_test.go に TestMulUnbalanced という新しいテスト関数が追加されました。このテストは、非常に長い数 x と非常に短い数 y の乗算を行い、その際に割り当てられるメモリ量が入力サイズの合計に対して妥当な範囲に収まっているか（10倍以内）を検証します。これは、以前の二次的な空間計算量の問題を直接的に検出するためのものです。
- allocBytes ヘルパー関数も追加され、特定の関数呼び出しによって割り当てられたバイト数を測定できるようになりました。

これらの変更により、math/big.Mul はアンバランスな入力に対しても、より予測可能で効率的なメモリ使用と再帰深度を実現し、大規模な数値計算の安定性とパフォーマンスを向上させています。

コアとなるコードの変更箇所

`src/pkg/math/big/nat.go`

--- a/src/pkg/math/big/nat.go
+++ b/src/pkg/math/big/nat.go
@@ -342,7 +342,7 @@ func alias(x, y nat) bool {
 	return cap(x) > 0 && cap(y) > 0 && &x[0:cap(x)][cap(x)-1] == &y[0:cap(y)][cap(y)-1]
 }
 
-// addAt implements z += x*(1<<(_W*i)); z must be long enough.\n
+// addAt implements z += x<<(_W*i); z must be long enough.\n
 // (we don't use nat.add because we need z to stay the same\n
 // slice, and we don't need to normalize z after each addition)\n
 func addAt(z, x nat, i int) {
@@ -405,8 +405,8 @@ func (z nat) mul(x, y nat) nat {
 
 	// determine Karatsuba length k such that
 	//
-	//   x = x1*b + x0\n
-	//   y = y1*b + y0  (and k <= len(y), which implies k <= len(x))\n
+	//   x = xh*b + x0  (0 <= x0 < b)\n
+	//   y = yh*b + y0  (0 <= y0 < b)\n
 	//   b = 1<<(_W*k)  ("base" of digits xi, yi)\n
 	//
 	k := karatsubaLen(n)
@@ -417,27 +417,41 @@ func (z nat) mul(x, y nat) nat {
 	y0 := y[0:k]              // y0 is not normalized
 	z = z.make(max(6*k, m+n)) // enough space for karatsuba of x0*y0 and full result of x*y
 	karatsuba(z, x0, y0)
-	z = z[0 : m+n] // z has final length but may be incomplete, upper portion is garbage\n
-
-	// If x1 and/or y1 are not 0, add missing terms to z explicitly:\n
-	//\n
-	//     m+n       2*k       0\n
-	//   z = [   ...   | x0*y0 ]\n
-	//     +   [ x1*y1 ]\n
-	//     +   [ x1*y0 ]\n
-	//     +   [ x0*y1 ]\n
+	z = z[0 : m+n]  // z has final length but may be incomplete\n
+	z[2*k:].clear() // upper portion of z is garbage (and 2*k <= m+n since k <= n <= m)\n
+
+	// If xh != 0 or yh != 0, add the missing terms to z. For\n
+	// \n
+	//   xh = xi*b^i + ... + x2*b^2 + x1*b (0 <= xi < b) \n
+	//   yh =                         y1*b (0 <= y1 < b) \n
+	// \n
+	// the missing terms are \n
+	// \n
+	//   x0*y1*b and xi*y0*b^i, xi*y1*b^(i+1) for i > 0 \n
+	// \n
+	// since all the yi for i > 1 are 0 by choice of k: If any of them \n
+	// were > 0, then yh >= b^2 and thus y >= b^2. Then k' = k*2 would \n
+	// be a larger valid threshold contradicting the assumption about k. \n
 	//\n
 	if k < n || m != n {
-		x1 := x[k:] // x1 is normalized because x is\n
-		y1 := y[k:] // y1 is normalized because y is\n
 		var t nat
-		t = t.mul(x1, y1)\n
-		copy(z[2*k:], t)\n
-		z[2*k+len(t):].clear() // upper portion of z is garbage\n
-		t = t.mul(x1, y0.norm())\n
-		addAt(z, t, k)\n
-		t = t.mul(x0.norm(), y1)\n
-		addAt(z, t, k)\n
+
+		// add x0*y1*b\n
+		x0 := x0.norm()\n
+		y1 := y[k:] // y1 is normalized because y is\n
+		addAt(z, t.mul(x0, y1), k)\n
+
+		// add xi*y0<<i, xi*y1*b<<(i+k)\n
+		y0 := y0.norm()\n
+		for i := k; i < len(x); i += k {\n
+			xi := x[i:]\n
+			if len(xi) > k {\n
+				xi = xi[:k]\n
+			}\n
+			xi = xi.norm()\n
+			addAt(z, t.mul(xi, y0), i)\n
+			addAt(z, t.mul(xi, y1), i+k)\n
+		}\n
 	}\n
 
 	return z.norm()\n

`src/pkg/math/big/nat_test.go`

--- a/src/pkg/math/big/nat_test.go
+++ b/src/pkg/math/big/nat_test.go
@@ -7,6 +7,7 @@ package big
 import (
 	"io"
 	"math/rand"
+	"runtime"
 	"strings"
 	"testing"
 )
@@ -63,6 +64,36 @@ var prodNN = []argNN{
 	{nat{0, 0, 991 * 991}, nat{0, 991}, nat{0, 991}},
 	{nat{1 * 991, 2 * 991, 3 * 991, 4 * 991}, nat{1, 2, 3, 4}, nat{991}},
 	{nat{4, 11, 20, 30, 20, 11, 4}, nat{1, 2, 3, 4}, nat{4, 3, 2, 1}},
+	// 3^100 * 3^28 = 3^128
+	{
+		natFromString("11790184577738583171520872861412518665678211592275841109096961"),
+		natFromString("515377520732011331036461129765621272702107522001"),
+		natFromString("22876792454961"),
+	},
+	// z = 111....1 (70000 digits)
+	// x = 10^(99*700) + ... + 10^1400 + 10^700 + 1
+	// y = 111....1 (700 digits, larger than Karatsuba threshold on 32-bit and 64-bit)
+	{
+		natFromString(strings.Repeat("1", 70000)),
+		natFromString("1" + strings.Repeat(strings.Repeat("0", 699)+"1", 99)),
+		natFromString(strings.Repeat("1", 700)),
+	},
+	// z = 111....1 (20000 digits)
+	// x = 10^10000 + 1
+	// y = 111....1 (10000 digits)
+	{
+		natFromString(strings.Repeat("1", 20000)),
+		natFromString("1" + strings.Repeat("0", 9999) + "1"),
+		natFromString(strings.Repeat("1", 10000)),
+	},
+}
+
+func natFromString(s string) nat {
+	x, _, err := nat(nil).scan(strings.NewReader(s), 0)
+	if err != nil {
+		panic(err)
+	}
+	return x
 }
 
 func TestSet(t *testing.T) {
@@ -136,6 +167,31 @@ func TestMulRangeN(t *testing.T) {
 	}
 }
 
+// allocBytes returns the number of bytes allocated by invoking f. 
+func allocBytes(f func()) uint64 {
+	var stats runtime.MemStats
+	runtime.ReadMemStats(&stats)
+	t := stats.TotalAlloc
+	f()
+	runtime.ReadMemStats(&stats)
+	return stats.TotalAlloc - t
+}
+
+// TestMulUnbalanced tests that multiplying numbers of different lengths
+// does not cause deep recursion and in turn allocate too much memory.
+// test case for issue 3807
+func TestMulUnbalanced(t *testing.T) {
+	x := rndNat(50000)
+	y := rndNat(40)
+	allocSize := allocBytes(func() {
+		nat(nil).mul(x, y)
+	})
+	inputSize := uint64(len(x)+len(y)) * _S
+	if ratio := allocSize / uint64(inputSize); ratio > 10 {
+		t.Errorf("multiplication uses too much memory (%d > %d times the size of inputs)", allocSize, ratio)
+	}
+}
+
 var rnd = rand.New(rand.NewSource(0x43de683f473542af))\n
 var mulx = rndNat(1e4)\n
 var muly = rndNat(1e4)\n

コアとなるコードの解説

`src/pkg/math/big/nat.go` の変更点

mul 関数は、x と y という2つの nat 型の多倍長自然数を乗算し、結果を z に格納します。

コメントの変更:
- addAt 関数のコメントが z += x*(1<<(_W*i)) から z += x<<(_W*i) に修正されました。これは、1<<(_W*i) が b^i を意味し、x を i 桁分シフトして加算するという操作をより正確に表現しています。
- Karatsuba分割における x と y の表現が x1*b + x0 から xh*b + x0 に変更されました。これは、x の上位部分を xh (x high) と表現することで、より一般的なKaratsubaの表記に近づけています。
メモリのクリア:
- z = z[0 : m+n] の後に z[2*k:].clear() が追加されました。これは、Karatsubaの計算で x0*y0 の結果が z の下位 2*k 桁に格納された後、それより上位のメモリ領域に以前のゴミデータが残っている可能性があるため、それを明示的にゼロクリアしています。これにより、後続の加算処理がクリーンな状態で行われ、メモリの再利用が促進されます。
アンバランスな入力の処理ロジックの改善:
- if k < n || m != n のブロックが、アンバランスな入力（k が n より小さい、または m と n が異なる場合）を処理する部分です。
- 旧実装:
  - x1 := x[k:] と y1 := y[k:] で上位部分を取得し、t.mul(x1, y1) で x1*y1 を計算し、z の 2*k オフセットにコピーしていました。
  - その後、t.mul(x1, y0.norm()) と t.mul(x0.norm(), y1) で x1*y0 と x0*y1 を計算し、それぞれ addAt(z, t, k) で z に加算していました。この mul 呼び出しが、x1 や y1 が短い場合でもKaratsubaの再帰を深くする原因となっていました。
- 新実装:
  - x0*y1*b の項を addAt(z, t.mul(x0, y1), k) で加算します。ここで x0 は正規化されています。y1 は y[k:] で、y の上位部分です。
  - 最も重要な変更は、for i := k; i < len(x); i += k ループの導入です。このループは、x の残りの部分 (x[k:]) を k 単位で処理します。
    - xi := x[i:] で x の現在の k 桁のチャンクを取得します。必要に応じて len(xi) > k で長さを調整します。
    - xi = xi.norm() でチャンクを正規化します。
    - addAt(z, t.mul(xi, y0), i) で xi*y0 を計算し、i オフセットで z に加算します。
    - addAt(z, t.mul(xi, y1), i+k) で xi*y1 を計算し、i+k オフセットで z に加算します。
  - このループにより、x が y よりもはるかに長い場合でも、x を小さなチャンクに分割して y0 や y1 と乗算し、結果を適切に加算できるようになりました。これにより、不必要な大きな一時配列の生成や深い再帰が回避され、メモリ使用量と再帰深度が大幅に改善されます。

`src/pkg/math/big/nat_test.go` の変更点

新しいテストケースの追加:
- prodNN 変数に、非常に大きな数と比較的短い数の乗算を含む新しいテストケースが追加されました。これらは、アンバランスな入力に対する mul 関数の正確性を検証します。
- natFromString ヘルパー関数が追加され、文字列から nat 型の数値を簡単に生成できるようになりました。
メモリ割り当てテストの追加:
- allocBytes 関数が追加されました。これは、引数として渡された関数 f が実行中に割り当てたメモリの総バイト数を測定します。runtime.MemStats を使用して、関数実行前後のメモリ統計を比較することで実現されます。
- TestMulUnbalanced 関数が追加されました。
  - このテストは、長さ50000のランダムな数 x と長さ40のランダムな数 y を生成します。
  - allocBytes を使用して nat(nil).mul(x, y) の実行中に割り当てられるメモリ量を測定します。
  - 割り当てられたメモリ量 allocSize が、入力サイズ (len(x)+len(y)) * _S （_S は各ワードのバイト数）の10倍を超えた場合、テストは失敗します。これは、以前の二次的な空間計算量の問題を直接的に検出するための重要なテストです。

参考にした情報源リンク

Karatsuba algorithm - Wikipedia: https://en.wikipedia.org/wiki/Karatsuba_algorithm
Go math/big package documentation: https://pkg.go.dev/math/big
Big-O notation - Wikipedia: https://en.wikipedia.org/wiki/Big_O_notation
Space complexity - Wikipedia: https://en.wikipedia.org/wiki/Space_complexity
Recursion (computer science) - Wikipedia: https://en.wikipedia.org/wiki/Recursion_(computer_science)

comemo