Keyboard shortcuts

Press ← or → to navigate between chapters

Press S or / to search in the book

Press ? to show this help

Press Esc to hide this help

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

このコミットは、Go言語のbignumパッケージにおける多倍長整数および有理数演算のサポートを大幅に強化し、既存の実装の改善、テストの拡充、およびコードのクリーンアップを行ったものです。特に、整数演算の完全なサポート（一部論理関数を除く）と有理数演算の完了に焦点が当てられています。

コミット

commit 7cd11c1c09729cc3ade1289014865ed26b22c354
Author: Robert Griesemer <gri@golang.org>
Date:   Tue Nov 4 09:33:15 2008 -0800

    - completed integer support (some logical functions missing)
    - completed rational support
    - better documentation
    - more tests
    - cleanups
    
    R=r
    OCL=18438
    CL=18438

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

https://github.com/golang/go/commit/7cd11c1c09729cc3ade1289014865ed26b22c354

元コミット内容

- completed integer support (some logical functions missing)
- completed rational support
- better documentation
- more tests
- cleanups

変更の背景

このコミットは、Go言語の初期開発段階におけるbignum（多倍長演算）パッケージの成熟化の一環として行われました。当時のGo言語には、標準で任意精度の数値演算を扱うための組み込み型やライブラリが不足していました。そのため、Robert Griesemer氏（Go言語の共同設計者の一人）によって、多倍長整数（Natural, Integer）および有理数（Rational）を扱うためのbignumパッケージが開発されていました。

このコミットの背景には、以下の目的があったと考えられます。

  1. 機能の拡充: 多倍長整数および有理数に対する基本的な算術演算（加算、減算、乗算、除算、剰余など）のサポートを完了させること。特に、整数演算の論理関数（AND, OR, XORなど）の一部が未実装であったり、有理数演算が不完全であったりした点を改善する必要がありました。
  2. 堅牢性の向上: 複雑な数値演算ライブラリにおいては、正確性と効率性が極めて重要です。既存のコードベースに対して、より多くのテストケースを追加し、エッジケースや潜在的なバグを特定・修正することで、ライブラリの堅牢性を高めることが求められました。
  3. コード品質の向上: 可読性、保守性、およびパフォーマンスを向上させるためのコードのクリーンアップとリファクタリング。これには、内部表現の最適化や、より明確な命名規則の適用などが含まれます。
  4. ドキュメンテーションの改善: 開発者や将来の利用者がライブラリを理解しやすくするために、内部実装の詳細やAPIの利用方法に関するドキュメンテーションを充実させること。

これらの変更は、Go言語がより広範な数値計算タスクに対応できるようになるための基盤を築く上で不可欠でした。

前提知識の解説

このコミットを理解するためには、以下の前提知識が役立ちます。

1. 多倍長整数演算 (Arbitrary-Precision Arithmetic)

通常のプログラミング言語で提供される整数型（例: int, long, int64など）は、表現できる数値の範囲が固定されています。しかし、暗号学、科学計算、金融アプリケーションなど、非常に大きな数値を扱う必要がある場合、これらの固定長整数では対応できません。多倍長整数演算は、メモリが許す限り任意の桁数の整数を表現し、それらに対して正確な算術演算を行う技術です。

  • 内部表現: 多倍長整数は通常、基数B（例: 2^32や2^64）を持つ「桁（digit）」の配列として表現されます。このコミットではDigit型（uint64）が使用されており、Bは1 << W（Wはビット幅）として定義されています。
  • 算術演算: 加算、減算、乗算、除算などの基本的な算術演算は、小学校で習う筆算と同様のアルゴリズムを、各桁に対して適用することで実現されます。ただし、桁上がり（carry）や桁借り（borrow）の処理、そして特に乗算や除算では、より複雑なアルゴリズム（例: Karatsuba法、Toom-Cook法、KnuthのAlgorithm Dなど）が用いられることがあります。

2. 有理数演算 (Rational Number Arithmetic)

有理数は、2つの整数の比a/b（b ≠ 0）として表される数です。浮動小数点数とは異なり、有理数演算は常に正確な結果を保証します。

  • 内部表現: 有理数は通常、分子（numerator）と分母（denominator）の2つの多倍長整数で構成されます。
  • 算術演算:
    • 加算: a/b + c/d = (ad + bc) / bd
    • 減算: a/b - c/d = (ad - bc) / bd
    • 乗算: a/b * c/d = ac / bd
    • 除算: (a/b) / (c/d) = ad / bc
  • 正規化: 有理数は、分子と分母の最大公約数（GCD）で割ることで正規化されます（例: 2/4は1/2に正規化）。これにより、一意な表現が保証され、計算が簡素化されます。

3. KnuthのAlgorithm D (多倍長除算アルゴリズム)

Donald Knuthの著書「The Art of Computer Programming, Volume 2: Seminumerical Algorithms」で詳述されているAlgorithm Dは、多倍長整数に対する効率的な除算アルゴリズムです。このアルゴリズムは、筆算の除算を多倍長数に拡張したもので、試行商（trial digit）の推定と修正を繰り返すことで商と剰余を求めます。

  • 正規化: 除算の前に、除数（分母）を正規化（最上位桁が特定の範囲に入るようにシフト）することで、試行商の推定精度を高めます。
  • 試行商の推定: 除数と被除数（分子）の最上位の数桁を用いて、試行商を推定します。この推定は、実際の商に近い値であることが重要です。
  • 乗算と減算: 推定された試行商と除数を乗算し、その結果を被除数から減算します。
  • 修正: 減算の結果が負になる場合、試行商が大きすぎたことを意味するため、試行商を減らし、被除数に除数を加算して修正します。

このコミットのDivMod関数（以前はDivMod2）は、このAlgorithm Dに基づいています。

4. Go言語の初期の型システムとuint

コミットの時期（2008年）はGo言語の初期段階であり、現在のGo言語とは異なる型定義や構文が見られます。例えば、uint型は現在のGo言語ではプラットフォーム依存の符号なし整数型ですが、このコミットではuintがuint64として扱われている箇所があります。また、panicの利用方法や、exportキーワード（現在のGoでは大文字で始まる識別子がエクスポートされる）など、初期のGo言語の特性が反映されています。

技術的詳細

このコミットにおける技術的な変更点は多岐にわたりますが、主要なものとして以下が挙げられます。

1. 内部表現の改善と定数の再定義

  • DigitとDigit2: 多倍長数の各桁を表す型としてDigit（uint64）とDigit2（uint32）が定義されています。Digit2は、特に多倍長除算において「半桁（half-digits）」として使用され、Digit型がプラットフォームの最大符号なし整数型である場合に、オーバーフローを避けて効率的な除算を行うために導入されています。
  • LogW, LogH, LogB:
    • LogW = 64: Digit型（uint64）のビット幅。
    • LogH = 4: 16進数1桁あたりのビット幅。
    • LogB = LogW - LogH: 多倍長数の基数Bのビット幅。これは、文字列との変換速度を向上させるために、10進数表現での除算・乗算が拡張精度演算なしで行えるようにするための「バッファ」として4ビットを確保していることを示唆しています。
  • W2, B2, M2: W2はLogBの半分（半桁のビット幅）、B2は半桁の基数、M2は半桁のマスクです。これらは、Mul11（単一桁の乗算）やDiv1（単一桁での除算）などの内部演算で、オーバーフローを避けるために使用されます。以前はL2, B2, M2という名前でしたが、W2に変更され、より意味が明確になりました。
  • W, B, M: WはW2 * 2（全桁のビット幅）、Bは全桁の基数、Mは全桁のマスクです。これらも以前はL, B, Mという名前でしたが、Wに変更され、半桁の定数との整合性が取られました。

これらの定数の変更は、多倍長演算の基数と桁のビット幅をより厳密に制御し、特に除算アルゴリズムの効率と正確性を向上させるためのものです。

2. Raw operationsの整理と追加

bignum.goのRaw operationsセクションは、多倍長数の内部配列（[]Digit）に対する低レベルな操作を定義しています。

  • 論理演算の移動: And1, And, Or1, Or, Xor1, Xorといった論理演算関数が、SupportセクションからRaw operationsセクションに移動されました。これにより、コードの構造がより論理的になりました。
  • シフト演算の修正: Shl（左シフト）とShr（右シフト）関数が修正され、シフト量sがW（全桁のビット幅）以下であることをassertで確認するようになりました。これにより、不正なシフト量に対する防御が強化されました。また、内部のビットシフト演算もLからWに、L-sからW-sに変更され、新しい定数定義に準拠しています。
  • Mul1とDiv1の追加: Digit2型を扱うMul1（単一桁での乗算）とDiv1（単一桁での除算）が追加されました。これらは、多倍長除算の内部で半桁を扱うために使用されます。

3. 多倍長除算アルゴリズムの改善 (DivMod)

このコミットの最も重要な変更点の一つは、多倍長除算関数DivMod（以前はDivMod2）の改善です。

  • UnpackとPackの正規化: Unpack関数はNatural型を[]Digit2に変換し、Pack関数は[]Digit2をNatural型に戻します。これらの関数は、変換後に結果を正規化（先行ゼロの除去）するようになりました。これにより、DivModの入力が常に正規化された状態であることが保証されます。
  • DivModのロジック改善:
    • ゼロ除算チェック: m == 0（除数がゼロ）の場合にpanicを発生させるようになりました。
    • 正規化係数fの計算: KnuthのAlgorithm Dでは、除数を正規化するために適切な係数fを乗算します。このコミットでは、fの計算ロジックが改善され、B2 / (Digit(y[m-1]) + 1)を使用するようになりました。また、f != 1の場合にのみ乗算を行うことで、不要な演算を避けています。
    • 試行商の推定: KnuthのAlgorithm Dに基づく試行商の推定ロジックがより明確に記述されています。x0, x1, x2といった被除数の上位桁と、y1（除数の最上位桁）を用いて試行商qを計算し、y2 * q > (x0*B2 + x1 - y1*q)*B2 + x2という条件でqを修正するループが含まれています。これは、試行商が大きすぎた場合に正確な値に調整するための重要なステップです。
    • 剰余の正規化解除: 除算の最後に、正規化のために乗算した係数fで剰余を割ることで、元のスケールに戻す処理が追加されました。

4. Natural（多倍長自然数）型の改善

  • Nat関数の最適化: Nat関数は、uint型の引数からNatural型を生成します。0, 1, 2, 10といった頻繁に使用される値に対しては、事前に定義されたグローバル変数（NatZero, NatOne, NatTwo, NatTen）を返すように最適化されました。これにより、これらの値の生成コストが削減されます。
  • Log2関数の修正: Log2関数は、Natural型のビット長を計算します。以前のバージョンではintを返していましたが、uintを返すように変更され、panic("Log2(0)")が追加され、ゼロの対数を計算しようとした場合にエラーを発生させるようになりました。
  • DivMod1関数の変更: DivMod1関数は、Natural型を単一桁dで除算し、商と剰余を返します。この関数は、レシーバメソッドから通常の関数に変更され、引数として*Naturalを受け取るようになりました。
  • String関数の改善: String関数は、Natural型を指定された基数で文字列に変換します。この関数は、DivMod1の変更に合わせて修正され、t.DivMod1(Digit(base))の代わりにDivMod1(t, Digit(base))を呼び出すようになりました。また、変換前にNatural型のコピーを作成することで、元の値を破壊しないように配慮されています。
  • NatFromString関数の追加: 文字列からNatural型を生成するNatFromString関数が追加されました。この関数は、基数（2進数、8進数、10進数、16進数）を自動的に判別する機能や、文字列のどの部分が数値として消費されたかを示すslen引数（ポインタ）をサポートしています。

5. Integer（多倍長整数）型の改善

  • MakeInt関数の追加: MakeInt関数は、符号と多倍長自然数からInteger型を生成します。ゼロのNaturalが渡された場合、符号をfalseに正規化するロジックが含まれています。
  • Int関数の改善: int型の引数からInteger型を生成するInt関数が改善されました。特に、最小負の整数（math.MinIntに相当）を正しく扱うためのロジックが追加されています。
  • 述語関数の追加: IsNeg(), IsPos()といった、整数の符号を判定する述語関数が追加されました。
  • 算術演算の修正: Add, Sub, Mul, Quo, Rem, QuoRemといった算術演算関数が、新しいMakeInt関数を使用するように修正されました。これにより、結果のIntegerが常に正規化された状態で生成されることが保証されます。
  • ユークリッド除算と剰余の導入: DivとMod関数が、C99スタイルの切り捨て除算（Quo, Rem）とは異なる、ユークリッド除算と剰余の定義に従うように変更されました。ユークリッド除算では、剰余は常に非負であり、0 <= m < |d|という条件を満たします。これは、数学的な性質がより明確であり、一部のアルゴリズムで好まれる定義です。
  • シフト演算の追加: Shl（左シフト）とShr（右シフト）関数がInteger型に追加されました。
  • 論理演算のプレースホルダー: And, Or, Xorといった論理演算関数が追加されましたが、これらはまだpanic("UNIMPLEMENTED")となっており、未実装であることが示されています。
  • Cmp関数の修正: Cmp関数は、2つのIntegerを比較します。符号を考慮した比較ロジックが改善されました。
  • IntFromString関数の改善: 文字列からInteger型を生成するIntFromString関数が改善され、符号の処理とslen引数のサポートが追加されました。

6. Rational（有理数）型の改善

  • MakeRat関数の追加: MakeRat関数は、分子*Integerと分母*NaturalからRational型を生成します。分子と分母の最大公約数（GCD）で割ることで、有理数を正規化するロジックが含まれています。
  • Rat関数の改善: int型の分子と分母からRational型を生成するRat関数が改善され、分母が負の場合に分子の符号を反転させるロジックが追加されました。
  • 述語関数の追加: IsZero(), IsNeg(), IsPos(), IsInt()といった、有理数の性質を判定する述語関数が追加されました。
  • 算術演算の修正: Add, Sub, Mul, Quoといった算術演算関数が、新しいMakeRat関数を使用するように修正されました。特に、Mul関数ではx.a.Mul(y.a)のようにIntegerの乗算を直接呼び出すようになりました。Quo関数は、DivからQuoに名前が変更され、分子と分母の符号を適切に処理するようになりました。
  • Cmp関数の修正: Cmp関数は、2つのRationalを比較します。x.a.MulNat(y.b)).Cmp(y.a.MulNat(x.b))というクロス乗算による比較ロジックが導入されました。
  • String関数の改善: String関数は、Rational型を指定された基数で文字列に変換します。分母が1の場合は整数として表示し、それ以外の場合は分子/分母の形式で表示するようになりました。
  • RatFromString関数の追加: 文字列からRational型を生成するRatFromString関数が追加されました。この関数は、分子/分母の形式の文字列を解析し、slen引数をサポートしています。

7. テストの拡充 (bignum_test.go)

  • テストヘルパー関数の追加: NAT_EQ, INT_EQ, RAT_EQといった、Natural, Integer, Rational型の比較を行うテストヘルパー関数が追加されました。これにより、テストコードの可読性と簡潔性が向上しました。
  • テストケースの追加と改善:
    • NatConv, IntConv, RatConv: 文字列と数値型間の変換に関するテストが大幅に拡充されました。特に、基数の自動判別やslen引数のテストが追加されています。
    • NatAdd, NatSub, NatMul, NatDiv, NatMod, NatShift: Natural型に対する加算、減算、乗算、除算、剰余、シフト演算のテストが追加・改善されました。特に、乗算と除算の逆演算関係のテストや、シフト演算の網羅的なテストが含まれています。
    • IntQuoRem, IntDivMod: Integer型に対するC99スタイルの切り捨て除算・剰余と、ユークリッド除算・剰余のテストが追加されました。これにより、異なる除算定義の正確性が検証されています。
    • IntShift: Integer型に対するシフト演算のテストが追加されました。
    • NatCmp, NatLog2, NatGcd, NatPow, NatPop: Natural型に対する比較、ビット長計算、最大公約数、べき乗、セットビット数カウントのテストが追加・改善されました。
  • テストメッセージの改善: TEST failed: ...という形式で、どのテストが失敗したかをより明確に表示するようになりました。

これらの変更は、bignumパッケージの機能性と堅牢性を大幅に向上させ、Go言語における高精度数値計算の基盤を強化するものでした。

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

このコミットにおけるコアとなるコードの変更箇所は、主にusr/gri/bignum/bignum.goファイルに集中しています。

  1. 定数定義の変更:

    --- a/usr/gri/bignum/bignum.go
+++ b/usr/gri/bignum/bignum.go
@@ -27,37 +27,56 @@ package Bignum
 // always normalized before returning the final result. The normalized
 // representation of 0 is the empty array (length = 0).
 //
+// The operations for all other numeric types are implemented on top of
+// the operations for natural numbers.
+//
 // The base B is chosen as large as possible on a given platform but there
 // are a few constraints besides the size of the largest unsigned integer
-// type available.
-// TODO describe the constraints.
+// type available:
+//
+// 1) To improve conversion speed between strings and numbers, the base B
+//    is chosen such that division and multiplication by 10 (for decimal
+//    string representation) can be done without using extended-precision
+//    arithmetic. This makes addition, subtraction, and conversion routines
+//    twice as fast. It requires a "buffer" of 4 bits per operand digit.
+//    That is, the size of B must be 4 bits smaller then the size of the
+//    type (Digit) in which these operations are performed. Having this
+//    buffer also allows for trivial (single-bit) carry computation in
+//    addition and subtraction (optimization suggested by Ken Thompson).
+//
+// 2) Long division requires extended-precision (2-digit) division per digit.
+//    Instead of sacrificing the largest base type for all other operations,
+//    for division the operands are unpacked into "half-digits", and the
+//    results are packed again. For faster unpacking/packing, the base size
+//    in bits must be even.
+
+type (
+	Digit  uint64;
+	Digit2 uint32;  // half-digits for division
+)
+

 const LogW = 64;
 const LogH = 4;  // bits for a hex digit (= "small" number)
-const LogB = LogW - LogH;
+const LogB = LogW - LogH;  // largest bit-width available
 
 
 const (
-	L2 = LogB / 2;
-	B2 = 1 << L2;
-	M2 = B2 - 1;
-
-	L = L2 * 2;
-	B = 1 << L;
-	M = B - 1;
-)
-
-
-type (
-	Digit2 uint32;
-	Digit  uint64;
+	// half-digits
+	W2 = LogB / 2;  // width
+	B2 = 1 << W2;   // base
+	M2 = B2 - 1;    // mask
+
+	// full digits
+	W = W2 * 2;     // width
+	B = 1 << W;     // base
+	M = B - 1;      // mask
 )

  2. Raw operationsの整理とMul11の修正:

    --- a/usr/gri/bignum/bignum.go
+++ b/usr/gri/bignum/bignum.go
@@ -154,31 +153,33 @@ func Sub(z, x, y *[]Digit) Digit {
 
 // Returns c = x*y div B, z = x*y mod B.
 func Mul11(x, y Digit) (Digit, Digit) {
-	// Split x and y into 2 sub-digits each (in base sqrt(B)),
+	// Split x and y into 2 sub-digits each,
 	// multiply the digits separately while avoiding overflow,\n
 	// and return the product as two separate digits.\n
 \n
-	const L0 = (L + 1)/2;\n
-	const L1 = L - L0;\n
-	const DL = L0 - L1;  // 0 or 1\n
-	const b  = 1<<L0;\n
-	const m  = b - 1;\n
+	// This code also works for non-even bit widths W\n
+	// which is why there are separate constants below\n
+	// for half-digits.\n
+	const W2 = (W + 1)/2;\n
+	const DW = W2*2 - W;  // 0 or 1\n
+	const B2  = 1<<W2;\n
+	const M2  = B2 - 1;\n
 \n     // split x and y into sub-digits\n    -	// x = (x1*b + x0)\n    -	// y = (y1*b + y0)\n    -	x1, x0 := x>>L0, x&m;\n    -	y1, y0 := y>>L0, y&m;\n    +	// x = (x1*B2 + x0)\n    +	// y = (y1*B2 + y0)\n    +	x1, x0 := x>>W2, x&M2;\n    +	y1, y0 := y>>W2, y&M2;\n     \n    -	// x*y = t2*b^2 + t1*b + t0\n    +	// x*y = t2*B2^2 + t1*B2 + t0\n     	t0 := x0*y0;\n     	t1 := x1*y0 + x0*y1;\n     	t2 := x1*y1;\n     \n     	// compute the result digits but avoid overflow\n     	// z = z1*B + z0 = x*y\n    -	z0 := (t1<<L0 + t0)&M;\n    -	z1 := t2<<DL + (t1 + t0>>L0)>>L1;\n    +	z0 := (t1<<W2 + t0)&M;\n    +	z1 := t2<<DW + (t1 + t0>>W2)>>(W-W2);\n     	\n     	return z1, z0;\n     }\n    ```

  3. DivMod関数の大幅な変更とMul1, Div1の追加:

    --- a/usr/gri/bignum/bignum.go
+++ b/usr/gri/bignum/bignum.go
@@ -470,34 +394,65 @@ func Pack(x *[]Digit2) *Natural {
 }
 
 
-// Division and modulo computation - destroys x and y. Based on the
-// algorithms described in:\n
+// func Mul1(z, x *[]Digit2, y Digit2) Digit2 {\n
+// \tn := len(x);\n
+// \tvar c Digit;\n
+// \tf := Digit(y);\n
+// \tfor i := 0; i < n; i++ {\n
+// \t\tt := c + Digit(x[i])*f;\n
+// \t\tc, z[i] = t>>W2, Digit2(t&M2);\n
+// \t}\n
+// \treturn Digit2(c);\n
+// }\n
+// \n
+// \n
+// func Div1(z, x *[]Digit2, y Digit2) Digit2 {\n
+// \tn := len(x);\n
+// \tvar c Digit;\n
+// \td := Digit(y);\n
+// \tfor i := n-1; i >= 0; i-- {\n
+// \t\tt := c*B2 + Digit(x[i]);\n
+// \t\tc, z[i] = t%d, Digit2(t/d);\n
+// \t}\n
+// \treturn Digit2(c);\n
+// }\n
+// \n
+// \n
+// DivMod returns q and r with x = y*q + r and 0 <= r < y.\n
+// x and y are destroyed in the process.\n
+//\n
+// The algorithm used here is based on 1). 2) describes the same algorithm\n
+// in C. A discussion and summary of the relevant theorems can be found in\n
+// 3). 3) also describes an easier way to obtain the trial digit - however\n
+// it relies on tripple-precision arithmetic which is why Knuth's method is\n
+// used here.\n
 //\n
 // 1) D. Knuth, "The Art of Computer Programming. Volume 2. Seminumerical\n
 //    Algorithms." Addison-Wesley, Reading, 1969.\n
+//    (Algorithm D, Sec. 4.3.1)\n
+//\n
+// 2) Henry S. Warren, Jr., "A Hacker's Delight". Addison-Wesley, 2003.\n
+//    (9-2 Multiword Division, p.140ff)\n
 //\n
-// 2) P. Brinch Hansen, Multiple-length division revisited: A tour of the\n
+// 3) P. Brinch Hansen, Multiple-length division revisited: A tour of the\n
 //    minefield. "Software - Practice and Experience 24", (June 1994),\n
 //    579-601. John Wiley & Sons, Ltd.\n
-//\n
-// Specifically, the inplace computation of quotient and remainder\n
-// is described in 1), while 2) provides the background for a more\n
-// accurate initial guess of the trial digit.\n
 \n    -func DivMod2(x, y *[]Digit2) (*[]Digit2, *[]Digit2) {\n    -	const b = B2;\n    -	\n    +func DivMod(x, y *[]Digit2) (*[]Digit2, *[]Digit2) {\n     	n := len(x);\n     	m := len(y);\n    -	assert(m > 0);  // division by zero\n    -	assert(n+1 <= cap(x));  // space for one extra digit (should it be == ?)\n    +	if m == 0 {\n    +		panic("division by zero");\n    +	}\n    +	assert(n+1 <= cap(x));  // space for one extra digit\n     	x = x[0 : n + 1];\n    +	assert(x[n] == 0);\n     	\n     	if m == 1 {\n     \t\t// division by single digit\n    -	\t\tx[0] = Div1(x[1 : n+1], x[0 : n], y[0]);\n    +	\t\tx[0] = Div1(x[1 : n+1], x[0 : n], y[0]);\n     \t\t\n     	} else if m > n {\n    -	\t\t// quotient = 0, remainder = x\n    -	\t\t// TODO in this case we shouldn't even split base - FIX THIS\n    +	\t\t// y > x => quotient = 0, remainder = x\n    +	\t\t// TODO in this case we shouldn't even unpack x and y\n     \t\tm = n;\n     \t\t\n     	} else {\n     \t\t// general case\n     \t\tassert(2 <= m && m <= n);\n    -	\t\tassert(x[n] == 0);\n     \t\t\n     \t\t// normalize x and y\n    -	\t\tf := b/(Digit(y[m-1]) + 1);\n    -	\t\tMul1(x, x, Digit2(f));\n    -	\t\tMul1(y, y, Digit2(f));\n    -	\t\tassert(b/2 <= y[m-1] && y[m-1] < b);  // incorrect scaling\n    +	\t\t// TODO Instead of multiplying, it would be sufficient to\n    +	\t\t//      shift y such that the normalization condition is\n    +	\t\t//      satisfied (as done in "Hacker's Delight").\n    +	\t\tf := B2 / (Digit(y[m-1]) + 1);\n    +	\t\tif f != 1 {\n    +	\t\t\tMul1(x, x, Digit2(f));\n    +	\t\t\tMul1(y, y, Digit2(f));\n    +	\t\t}\n    +	\t\tassert(B2/2 <= y[m-1] && y[m-1] < B2);  // incorrect scaling\n     \t\t\n     \t\ty1, y2 := Digit(y[m-1]), Digit(y[m-2]);\n    -	\t\td2 := Digit(y1)*b + Digit(y2);\n    +	\t\td2 := Digit(y1)<<W2 + Digit(y2);\n     \t\tfor i := n-m; i >= 0; i-- {\n     \t\t\tk := i+m;\n     \t\t\t\n    -	\t\t\t// compute trial digit\n    +	\t\t\t// compute trial digit (Knuth)\n     \t\t\tvar q Digit;\n    -	\t\t\t{\t// Knuth\n    -	\t\t\t\tx0, x1, x2 := Digit(x[k]), Digit(x[k-1]), Digit(x[k-2]);\n    +	\t\t\t{\tx0, x1, x2 := Digit(x[k]), Digit(x[k-1]), Digit(x[k-2]);\n     \t\t\t\tif x0 != y1 {\n    -	\t\t\t\t\tq = (x0*b + x1)/y1;\n    +	\t\t\t\t\tq = (x0<<W2 + x1)/y1;\n     \t\t\t\t} else {\n    -	\t\t\t\t\tq = b-1;\n    +	\t\t\t\t\tq = B2 - 1;\n     \t\t\t\t}\n    -	\t\t\t\tfor y2 * q > (x0*b + x1 - y1*q)*b + x2 {\n    +	\t\t\t\tfor y2*q > (x0<<W2 + x1 - y1*q)<<W2 + x2 {\n     \t\t\t\t\tq--\n     \t\t\t\t}\n     \t\t\t}\n    @@ -542,8 +500,8 @@ func DivMod2(x, y *[]Digit2) (*[]Digit2, *[]Digit2) {\n     \t\t// subtract y*q\n     \t\tc := Digit(0);\n     \t\tfor j := 0; j < m; j++ {\n    -	\t\t\tt := c + Digit(x[i+j]) - Digit(y[j])*q;  // arithmetic shift!\n    -	\t\t\tc, x[i+j] = Digit(int64(t)>>L2), Digit2(t&M2);\n    +	\t\t\tt := c + Digit(x[i+j]) - Digit(y[j])*q;\n    +	\t\t\tc, x[i+j] = Digit(int64(t)>>W2), Digit2(t&M2);  // requires arithmetic shift!\n     \t\t\t}\n     \t\t\t\n     \t\t\t// correct if trial digit was too large\n    @@ -552,7 +510,7 @@ func DivMod2(x, y *[]Digit2) (*[]Digit2, *[]Digit2) {\n     \t\t\tc := Digit(0);\n     \t\t\tfor j := 0; j < m; j++ {\n     \t\t\t\tt := c + Digit(x[i+j]) + Digit(y[j]);\n    -	\t\t\t\tc, x[i+j] = uint64(int64(t) >> L2), Digit2(t & M2)\n    +	\t\t\t\tc, x[i+j] = t >> W2, Digit2(t & M2)\n     \t\t\t}\n     \t\t\tassert(c + Digit(x[k]) == 0);\n     \t\t\t// correct trial digit\n    @@ -563,8 +521,10 @@ func DivMod2(x, y *[]Digit2) (*[]Digit2, *[]Digit2) {\n     \t\t}\n     \t\t\n     \t\t// undo normalization for remainder\n    -	\t\tc := Div1(x[0 : m], x[0 : m], Digit2(f));\n    -	\t\tassert(c == 0);\n    +	\t\tif f != 1 {\n    +	\t\t\tc := Div1(x[0 : m], x[0 : m], Digit2(f));\n    +	\t\t\tassert(c == 0);\n    +	\t\t}\n     \t}\n     \n     	return x[m : n+1], x[0 : m];\n    ```

  4. Natural型のNatFromStringとDivMod1の変更:

    --- a/usr/gri/bignum/bignum.go
+++ b/usr/gri/bignum/bignum.go
@@ -689,26 +676,25 @@ func (x *Natural) String(base uint) string {
 
 
 func (x *Natural) String(base uint) string {\n    -	if x.IsZero() {\n    +	if len(x) == 0 {\n     \t\treturn "0";\n     \t}\n     \t\n    -	// allocate string\n    +	// allocate buffer for conversion\n     \tassert(2 <= base && base <= 16);\n    -	n := (x.Log2() + 1) / Log2(Digit(base)) + 1;  // TODO why the +1?\n    +	n := (x.Log2() + 1) / Log2(Digit(base)) + 1;  // +1: round up\n     \ts := new([]byte, n);\n     \n    -	// convert\n    -\n    -	// don't destroy x, make a copy\n    +	// don't destroy x\n     \tt := new(Natural, len(x));\n    -	Or1(t, x, 0);  // copy x\n    +	Or1(t, x, 0);  // copy\n     \t\n    +	// convert\n     \ti := n;\n     \tfor !t.IsZero() {\n     \t\ti--;\n     \t\tvar d Digit;\n    -	\t\tt, d = t.DivMod1(Digit(base));\n    +	\t\tt, d = DivMod1(t, Digit(base));\n     \t\ts[i] = "0123456789abcdef"[d];\n     \t};\n     \n    @@ -716,7 +702,104 @@ func (x *Natural) String(base uint) string {
 }
 
 
-export func MulRange(a, b Digit) *Natural {\n    +func HexValue(ch byte) uint {\n    +	d := uint(1 << LogH);\n    +	switch {\n    +	case '0' <= ch && ch <= '9': d = uint(ch - '0');\n    +	case 'a' <= ch && ch <= 'f': d = uint(ch - 'a') + 10;\n    +	case 'A' <= ch && ch <= 'F': d = uint(ch - 'A') + 10;\n    +	}\n    +	return d;\n    +}\n    +\n    +\n    +// Computes x = x*d + c for "small" d's.\n    +func MulAdd1(x *Natural, d, c Digit) *Natural {\n    +	assert(IsSmall(d-1) && IsSmall(c));\n    +	n := len(x);\n    +	z := new(Natural, n + 1);\n    +\n    +	for i := 0; i < n; i++ {\n    +		t := c + x[i]*d;\n    +		c, z[i] = t>>W, t&M;\n    +	}\n    +	z[n] = c;\n    +\n    +	return Normalize(z);\n    +}\n    +\n    +\n    +// Determines base (octal, decimal, hexadecimal) if base == 0.\n    +export func NatFromString(s string, base uint, slen *int) *Natural {\n    +	// determine base if necessary\n    +	i, n := 0, len(s);\n    +	if base == 0 {\n    +		base = 10;\n    +		if n > 0 && s[0] == '0' {\n    +			if n > 1 && (s[1] == 'x' || s[1] == 'X') {\n    +				base, i = 16, 2;\n    +			} else {\n    +				base, i = 8, 1;\n    +			}\n    +		}\n    +	}\n    +	\n    +	// convert string\n    +	assert(2 <= base && base <= 16);\n    +	x := Nat(0);\n    +	for ; i < n; i++ {\n    +		d := HexValue(s[i]);\n    +		if d < base {\n    +			x = MulAdd1(x, Digit(base), Digit(d));\n    +		} else {\n    +			break;\n    +		}\n    +	}\n    +\n    +	// provide number of string bytes consumed if necessary\n    +	if slen != nil {\n    +		*slen = i;\n    +	}\n    +\n    +	return x;\n    +}\n    +\n    +\n    +// Natural number functions\n    +\n    +func Pop1(x Digit) uint {\n    +	n := uint(0);\n    +	for x != 0 {\n    +		x &= x-1;\n    +		n++;\n    +	}\n    +	return n;\n    +}\n    +\n    +\n    +func (x *Natural) Pop() uint {\n    +	n := uint(0);\n    +	for i := len(x) - 1; i >= 0; i-- {\n    +		n += Pop1(x[i]);\n    +	}\n    +	return n;\n    +}\n    +\n    +\n    +func (x *Natural) Pow(n uint) *Natural {\n    +	z := Nat(1);\n    +	for n > 0 {\n    +		// z * x^n == x^n0\n    +		if n&1 == 1 {\n    +			z = z.Mul(x);\n    +		}\n    +		x, n = x.Mul(x), n/2;\n    +	}\n    +	return z;\n    +}\n    +\n    +\n    +export func MulRange(a, b uint) *Natural {\n     	switch {\n     	case a > b: return Nat(1);\n     	case a == b: return Nat(a);\n    ```

  5. Integer型のIntFromStringとユークリッド除算の導入:

    --- a/usr/gri/bignum/bignum.go
+++ b/usr/gri/bignum/bignum.go
@@ -908,13 +1103,23 @@ func (x *Integer) String(base uint) string {
 }
 
 	\n    -export func IntFromString(s string, base uint) *Integer {\n    +// Determines base (octal, decimal, hexadecimal) if base == 0.\n    +export func IntFromString(s string, base uint, slen *int) *Integer {\n     	// get sign, if any\n     	sign := false;\n     	if len(s) > 0 && (s[0] == '-' || s[0] == '+') {\n     \t\tsign = s[0] == '-';\n    +	\ts = s[1 : len(s)];\n    +	}\n    +\n    +	z := MakeInt(sign, NatFromString(s, base, slen));\n    +\n    +	// correct slen if necessary\n    +	if slen != nil && sign {\n    +		*slen++;\n     \t}\n    -	return &Integer{sign, NatFromString(s[1 : len(s)], base)};\n    +\n    +	return z;\n     }

  6. Rational型のRatFromStringとMakeRatの追加:

    --- a/usr/gri/bignum/bignum.go
+++ b/usr/gri/bignum/bignum.go
@@ -922,56 +1127,119 @@ export func IntFromString(s string, base uint) *Integer {
 // Rational numbers
 
 export type Rational struct {\n    -	a, b *Integer;  // a = numerator, b = denominator\n    +	a *Integer;  // numerator\n    +	b *Natural;  // denominator\n     }\n     
 
-func (x *Rational) Normalize() *Rational {\n    -	f := x.a.mant.Gcd(x.b.mant);\n    -	x.a.mant = x.a.mant.Div(f);\n    -	x.b.mant = x.b.mant.Div(f);\n    -	return x;\n    +// Creation\n    +\n    +export func MakeRat(a *Integer, b *Natural) *Rational {\n    +	f := a.mant.Gcd(b);  // f > 0\n    +	if f.Cmp(Nat(1)) != 0 {\n    +		a = MakeInt(a.sign, a.mant.Div(f));\n    +		b = b.Div(f);\n    +	}\n    +	return &Rational{a, b};\n    +}\n    +\n    +\n    +export func Rat(a0 int, b0 int) *Rational {\n    +	a, b := Int(a0), Int(b0);\n    +	if b.sign {\n    +		a = a.Neg();\n    +	}\n    +	return MakeRat(a, b.mant);\n     }\n     
 
-func Rat(a, b *Integer) *Rational {\n    -	return (&Rational{a, b}).Normalize();\n    +// Predicates\n    +\n    +func (x *Rational) IsZero() bool {\n    +	return x.a.IsZero();\n    +}\n    +\n    +\n    +func (x *Rational) IsNeg() bool {\n    +	return x.a.IsNeg();\n    +}\n    +\n    +\n    +func (x *Rational) IsPos() bool {\n    +	return x.a.IsPos();\n    +}\n    +\n    +\n    +func (x *Rational) IsInt() bool {\n    +	return x.b.Cmp(Nat(1)) == 0;\n    +}\n    +\n    +\n    +// Operations\n    +\n    +func (x *Rational) Neg() *Rational {\n    +	return MakeRat(x.a.Neg(), x.b);\n     }\n     
 
 func (x *Rational) Add(y *Rational) *Rational {\n    -	return Rat((x.a.Mul(y.b)).Add(x.b.Mul(y.a)), x.b.Mul(y.b));\n    +	return MakeRat((x.a.MulNat(y.b)).Add(y.a.MulNat(x.b)), x.b.Mul(y.b));\n     }\n     
 
 func (x *Rational) Sub(y *Rational) *Rational {\n    -	return Rat((x.a.Mul(y.b)).Sub(x.b.Mul(y.a)), x.b.Mul(y.b));\n    +	return MakeRat((x.a.MulNat(y.b)).Sub(y.a.MulNat(x.b)), x.b.Mul(y.b));\n     }\n     
 
 func (x *Rational) Mul(y *Rational) *Rational {\n    -	return Rat(x.a.Mul(y.a), x.b.Mul(y.b));\n    +	return MakeRat(x.a.Mul(y.a), x.b.Mul(y.b));\n     }\n     
 
-func (x *Rational) Div(y *Rational) *Rational {\n    -	return Rat(x.a.Mul(y.b), x.b.Mul(y.a));\n    +func (x *Rational) Quo(y *Rational) *Rational {\n    +	a := x.a.MulNat(y.b);\n    +	b := y.a.MulNat(x.b);\n    +	if b.IsNeg() {\n    +		a = a.Neg();\n    +	}\n    +	return MakeRat(a, b.mant);\n     }\n     
 
-func (x *Rational) Mod(y *Rational) *Rational {\n    -	panic("UNIMPLEMENTED");\n    -	return nil;\n    +func (x *Rational) Cmp(y *Rational) int {\n    +	return (x.a.MulNat(y.b)).Cmp(y.a.MulNat(x.b));\n     }\n     
 
-func (x *Rational) Cmp(y *Rational) int {\n    -	panic("UNIMPLEMENTED");\n    -	return 0;\n    +func (x *Rational) String(base uint) string {\n    +	s := x.a.String(base);\n    +	if !x.IsInt() {\n    +		s += "/" + x.b.String(base);\n    +	}\n    +	return s;\n     }\n     
 
-export func RatFromString(s string) *Rational {\n    -	panic("UNIMPLEMENTED");\n    -	return nil;\n    +// Determines base (octal, decimal, hexadecimal) if base == 0.\n    +export func RatFromString(s string, base uint, slen *int) *Rational {\n    +	// read nominator\n    +	var alen, blen int;\n    +	a := IntFromString(s, base, &alen);\n    +	b := Nat(1);\n    +	\n    +	// read denominator, if any\n    +	if alen < len(s) && s[alen] == '/' {\n    +		alen++;\n    +		if alen < len(s) {\n    +			b = NatFromString(s[alen : len(s)], base, &blen);\n    +		}\n    +	}\n    +	\n    +	// provide number of string bytes consumed if necessary\n    +	if slen != nil {\n    +		*slen = alen + blen;\n    +	}\n    +\n    +	return MakeRat(a, b);\n     }

コアとなるコードの解説

1. 定数定義の変更

以前のL, L2といった定数名が、W, W2に変更されました。これは、Lが「Length」を意味するのに対し、Wが「Width」（ビット幅）を意味するため、より直感的で正確な命名です。特に、LogB = LogW - LogHという定義は、10進数変換の効率化のために4ビットのバッファを確保するという設計思想を明確に示しています。DigitとDigit2の型定義も、多倍長除算における半桁の概念を導入するために重要です。

2. Mul11関数の修正

Mul11は、単一のDigit同士の乗算を行い、その結果を2つのDigit（上位桁と下位桁）で返します。この関数は、多倍長乗算の基本的なビルディングブロックです。変更点では、内部で使用される定数名がL0, L1, b, mからW2, DW, B2, M2に変更され、新しい定数定義に準拠しています。これにより、半桁の概念がより明確に反映され、オーバーフローを避けるための計算が正確に行われます。

3. DivMod関数の大幅な変更とMul1, Div1の追加

DivMod関数は、多倍長除算の核心部分です。このコミットでは、KnuthのAlgorithm Dの実装が大幅に改善されました。

  • Mul1とDiv1: これらの関数は、Digit2（半桁）を扱うための単一桁乗算と除算です。DivMod内で、多倍長数を半桁の配列に「アンパック」してこれらの関数を呼び出すことで、Digit型が最大符号なし整数型である場合でも、オーバーフローを避けて正確な除算を実行できるようにしています。
  • DivModのロジック:
    • 正規化: 除数yを正規化するために、適切な係数fを計算し、xとyの両方に乗算します。これにより、試行商の推定精度が向上します。
    • 試行商の推定と修正: KnuthのAlgorithm Dの主要な部分であり、被除数と除数の上位桁から試行商qを推定し、そのqが大きすぎた場合に修正するループが含まれています。この修正ステップは、除算の正確性を保証するために不可欠です。
    • 剰余の正規化解除: 除算の最後に、正規化のために乗算した係数fで剰余を割ることで、元のスケールに戻します。

これらの変更により、DivMod関数はより堅牢で正確な多倍長除算を提供できるようになりました。

4. Natural型のNatFromStringとDivMod1の変更

  • NatFromString: 文字列からNatural型を生成するこの関数は、基数（10進数、16進数、8進数）の自動判別機能と、解析された文字列の長さを返すslen引数をサポートします。これは、ユーザー入力の解析や、より複雑な文字列処理において非常に便利です。
  • DivMod1: この関数は、Natural型を単一のDigitで除算します。以前はレシーバメソッドでしたが、通常の関数に変更され、引数として*Naturalを受け取るようになりました。これにより、関数の再利用性が向上し、String関数内でのコピー処理がより明確になりました。

5. Integer型のIntFromStringとユークリッド除算の導入

  • IntFromString: NatFromStringと同様に、符号の処理とslen引数のサポートが追加され、文字列からInteger型をより柔軟に解析できるようになりました。
  • ユークリッド除算: DivとMod関数が、C99スタイルの切り捨て除算とは異なるユークリッド除算の定義に従うように変更されました。ユークリッド除算は、剰余が常に非負であるという特性を持ち、数学的な性質がより一貫しています。これは、特に数論的なアルゴリズムを実装する際に重要となる場合があります。

6. Rational型のRatFromStringとMakeRatの追加

  • MakeRat: 分子と分母からRational型を生成し、自動的に正規化（最大公約数で割る）を行います。これにより、有理数が常に最も簡約された形で表現されることが保証されます。
  • RatFromString: 文字列からRational型を生成するこの関数は、分子/分母の形式の文字列を解析し、MakeRatを使用して正規化された有理数を返します。

これらのコアとなるコードの変更は、bignumパッケージの機能性、正確性、効率性を大幅に向上させ、Go言語における高精度数値計算の基盤をより強固なものにしました。

関連リンク

  • Go言語の公式ドキュメント (現在のmath/bigパッケージ): https://pkg.go.dev/math/big
  • Donald Knuth, "The Art of Computer Programming, Volume 2: Seminumerical Algorithms."

参考にした情報源リンク

  • コミットメッセージと変更されたファイルの内容 (usr/gri/bignum/bignum.go, usr/gri/bignum/bignum_test.go)
  • Go言語の歴史と初期開発に関する一般的な知識
  • 多倍長整数演算と有理数演算に関する一般的なアルゴリズム知識
  • KnuthのAlgorithm Dに関する情報 (特にコミット内のコメントで参照されている文献)
    • D. Knuth, "The Art of Computer Programming. Volume 2. Seminumerical Algorithms." Addison-Wesley, Reading, 1969. (Algorithm D, Sec. 4.3.1)
    • Henry S. Warren, Jr., "A Hacker's Delight". Addison-Wesley, 2003. (9-2 Multiword Division, p.140ff)
    • P. Brinch Hansen, Multiple-length division revisited: A tour of the minefield. "Software - Practice and Experience 24", (June 1994), 579-601. John Wiley & Sons, Ltd.
  • ユークリッド除算に関する情報
    • Raymond T. Boute, The Euclidian definition of the functions div and mod. "ACM Transactions on Programming Languages and Systems (TOPLAS)", 14(2):127-144, New York, NY, USA, 4/1992. ACM press.
  • Go言語の初期の構文や慣習に関する情報 (現在のGo言語との比較)