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

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

このコミットは、Go言語の初期のbignumパッケージ、すなわち任意精度算術(多倍長整数)を扱うライブラリに対するものです。bignumパッケージは、通常のプリミティブ型では表現できない非常に大きな整数を扱うための基本的な算術演算(加算、減算、乗算、除算、剰余など)を提供します。このファイルbignum.goは、その中核となる実装を含んでいます。bignum_test.goは、これらの算術演算が正しく機能することを検証するためのテストコードです。

コミット

このコミットは、多倍長整数の除算および剰余演算のアルゴリズムを高速化し、コードベースをリファクタリングしてクリーンアップすることを目的としています。具体的には、内部的な演算単位を20ビットから30ビットに変更することで、除算/剰余演算の効率を向上させています。また、多数の新しいテストが追加され、ライブラリとしての完成度を高めています。

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

https://github.com/golang/go/commit/78b0013a07cc24557f3887a89cde283fe0c664ef

元コミット内容

commit 78b0013a07cc24557f3887a89cde283fe0c664ef
Author: Robert Griesemer <gri@golang.org>
Date:   Mon Nov 3 09:21:10 2008 -0800

    - changed general div/mod implementation to a faster algorithm
      (operates on 30bit values at a time instead of 20bit values)
    - refactored and cleaned up lots of code
    - more tests
    - close to check-in as complete library
    
    R=r
    OCL=18326
    CL=18326
---
 usr/gri/bignum/bignum.go      | 557 ++++++++++++++++++++++++------------------
 usr/gri/bignum/bignum_test.go |  44 +++-
 2 files changed, 367 insertions(+), 234 deletions(-)

diff --git a/usr/gri/bignum/bignum.go b/usr/gri/bignum/bignum.go
index 0821720048..0852f42262 100755
--- a/usr/gri/bignum/bignum.go
+++ b/usr/gri/bignum/bignum.go
@@ -38,22 +38,18 @@ const LogB = LogW - LogH;
 
 
 const (
-	L3 = LogB / 3;
-	B3 = 1 << L3;
-	M3 = B3 - 1;
-
-	L2 = L3 * 2;
+	L2 = LogB / 2;
  	B2 = 1 << L2;
  	M2 = B2 - 1;
 
-	L = L3 * 3;
+	L = L2 * 2;
  	B = 1 << L;
  	M = B - 1;
 )
 
 
 type (
-	Digit3 uint32;
+	Digit2 uint32;
  	Digit  uint64;
 )
 
@@ -69,26 +65,205 @@ func assert(p bool) {
 }
 
 
-func IsSmall(x Digit) bool {
-	return x < 1<<LogH;
+// ----------------------------------------------------------------------------
+// Raw operations
+
+func And1(z, x *[]Digit, y Digit) {
+	for i := len(x) - 1; i >= 0; i-- {
+		z[i] = x[i] & y;
+	}
 }
 
 
-func Split(x Digit) (Digit, Digit) {
-	return x>>L, x&M;
+func And(z, x, y *[]Digit) {
+	for i := len(x) - 1; i >= 0; i-- {
+		z[i] = x[i] & y[i];
+	}
 }
 
 
-export func Dump(x *[]Digit) {
-	print("[", len(x), "]");
+func Or1(z, x *[]Digit, y Digit) {
  	for i := len(x) - 1; i >= 0; i-- {
-		print(" ", x[i]);
+		z[i] = x[i] | y;
  	}
-	println();
 }
 
 
-export func Dump3(x *[]Digit3) {
+func Or(z, x, y *[]Digit) {
+	for i := len(x) - 1; i >= 0; i-- {
+		z[i] = x[i] | y[i];
+	}
+}
+
+
+func Xor1(z, x *[]Digit, y Digit) {
+	for i := len(x) - 1; i >= 0; i-- {
+		z[i] = x[i] ^ y;
+	}
+}
+
+
+func Xor(z, x, y *[]Digit) {
+	for i := len(x) - 1; i >= 0; i-- {
+		z[i] = x[i] ^ y[i];
+	}
+}
+
+
+func Add1(z, x *[]Digit, c Digit) Digit {
+	n := len(x);
+	for i := 0; i < n; i++ {
+		t := c + x[i];
+		c, z[i] = t>>L, t&M
+	}
+	return c;
+}
+
+
+func Add(z, x, y *[]Digit) Digit {
+	var c Digit;
+	n := len(x);
+	for i := 0; i < n; i++ {
+		t := c + x[i] + y[i];
+		c, z[i] = t>>L, t&M
+	}
+	return c;
+}
+
+
+func Sub1(z, x *[]Digit, c Digit) Digit {
+	n := len(x);
+	for i := 0; i < n; i++ {
+		t := c + x[i];
+		c, z[i] = Digit(int64(t)>>L), t&M;  // arithmetic shift!
+	}
+	return c;
+}
+
+
+func Sub(z, x, y *[]Digit) Digit {
+	var c Digit;
+	n := len(x);
+	for i := 0; i < n; i++ {
+		t := c + x[i] - y[i];
+		c, z[i] = Digit(int64(t)>>L), t&M;  // arithmetic shift!
+	}
+	return c;
+}
+
+
+// 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)),
+	// multiply the digits separately while avoiding overflow,
+	// and return the product as two separate digits.
+
+	const L0 = (L + 1)/2;
+	const L1 = L - L0;
+	const DL = L0 - L1;  // 0 or 1
+	const b  = 1<<L0;
+	const m  = b - 1;
+
+	// split x and y into sub-digits
+	// x = (x1*b + x0)
+	// y = (y1*b + y0)
+	x1, x0 := x>>L0, x&m;
+	y1, y0 := y>>L0, y&m;
+
+	// x*y = t2*b^2 + t1*b + t0
+	t0 := x0*y0;
+	t1 := x1*y0 + x0*y1;
+	t2 := x1*y1;
+
+	// compute the result digits but avoid overflow
+	// z = z1*B + z0 = x*y
+	z0 := (t1<<L0 + t0)&M;
+	z1 := t2<<DL + (t1 + t0>>L0)>>L1;
+	
+	return z1, z0;
+}
+
+
+func Mul(z, x, y *[]Digit) {
+	n := len(x);
+	m := len(y);
+	for j := 0; j < m; j++ {
+		d := y[j];
+		if d != 0 {
+			c := Digit(0);
+			for i := 0; i < n; i++ {
+				// z[i+j] += c + x[i]*d;
+				z1, z0 := Mul11(x[i], d);
+				t := c + z[i+j] + z0;
+				c, z[i+j] = t>>L, t&M;
+				c += z1;
+			}
+			z[n+j] = c;
+		}
+	}
+}
+
+
+func Mul1(z, x *[]Digit2, y Digit2) Digit2 {
+	n := len(x);
+	var c Digit;
+	f := Digit(y);
+	for i := 0; i < n; i++ {
+		t := c + Digit(x[i])*f;
+		c, z[i] = t>>L2, Digit2(t&M2);
+	}
+	return Digit2(c);
+}
+
+
+func Div1(z, x *[]Digit2, y Digit2) Digit2 {
+	n := len(x);
+	var c Digit;
+	d := Digit(y);
+	for i := n-1; i >= 0; i-- {
+		t := c*B2 + Digit(x[i]);
+		c, z[i] = t%d, Digit2(t/d);
+	}
+	return Digit2(c);
+}
+
+
+func Shl(z, x *[]Digit, s uint) Digit {
+	assert(s <= L);
+	n := len(x);
+	var c Digit;
+	for i := 0; i < n; i++ {
+		c, z[i] = x[i] >> (L-s), x[i] << s & M | c;
+	}
+	return c;
+}
+
+
+func Shr(z, x *[]Digit, s uint) Digit {
+	assert(s <= L);
+	n := len(x);
+	var c Digit;
+	for i := n - 1; i >= 0; i-- {
+		c, z[i] = x[i] << (L-s) & M, x[i] >> s | c;
+	}
+	return c;
+}
+
+
+// ----------------------------------------------------------------------------
+// Support
+
+func IsSmall(x Digit) bool {
+	return x < 1<<LogH;
+}
+
+
+func Split(x Digit) (Digit, Digit) {
+	return x>>L, x&M;
+}
+
+
+export func Dump(x *[]Digit) {
  	print("[", len(x), "]");
  	for i := len(x) - 1; i >= 0; i-- {
  		print(" ", x[i]);
@@ -140,7 +315,7 @@ func Normalize(x *Natural) *Natural {
 }
 
 
-func Normalize3(x *[]Digit3) *[]Digit3 {
+func Normalize2(x *[]Digit2) *[]Digit2 {
  	n := len(x);
  	for n > 0 && x[n - 1] == 0 { n-- }
  	if n < len(x) {
@@ -150,9 +325,10 @@ func Normalize3(x *[]Digit3) *[]Digit3 {
 }
 
 
-func (x *Natural) IsZero() bool {
-	return len(x) == 0;
-}
+// Predicates
+
+func (x *Natural) IsZero() bool { return len(x) == 0; }
+func (x *Natural) IsOdd() bool { return len(x) > 0 && x[0]&1 != 0; }
 
 
 func (x *Natural) Add(y *Natural) *Natural {
@@ -161,13 +337,10 @@ func (x *Natural) Add(y *Natural) *Natural {
  	if n < m {
  		return y.Add(x);
  	}
-	assert(n >= m);
-	z := new(Natural, n + 1);
-
-	c := Digit(0);
-	for i := 0; i < m; i++ { c, z[i] = Split(c + x[i] + y[i]); }
-	for i := m; i < n; i++ { c, z[i] = Split(c + x[i]); }
-	z[n] = c;
+	z := new(Natural, n + 1);
+	c := Add(z[0 : m], x[0 : m], y);
+	z[n] = Add1(z[m : n], x[m : n], c);
 
  	return Normalize(z);
 }
@@ -176,19 +349,15 @@ func (x *Natural) Sub(y *Natural) *Natural {
 func (x *Natural) Sub(y *Natural) *Natural {
  	n := len(x);
  	m := len(y);
-	assert(n >= m);
-	z := new(Natural, n);
-
-	c := Digit(0);
-	for i := 0; i < m; i++ {
-		t := c + x[i] - y[i];
-		c, z[i] = Digit(int64(t)>>L), t&M;  // arithmetic shift!
+	if n < m {
+		panic("underflow")
 	}
-	for i := m; i < n; i++ {
-		t := c + x[i];
-		c, z[i] = Digit(int64(t)>>L), t&M;  // arithmetic shift!
+
+	z := new(Natural, n);
+	c := Sub(z[0 : m], x[0 : m], y);
+	if Sub1(z[m : n], x[m : n], c) != 0 {
+		panic("underflow");
 	}
-	assert(c == 0);  // x.Sub(y) must be called with x >= y
  	return Normalize(z);
 }
 
@@ -207,56 +376,12 @@ func (x* Natural) MulAdd1(a, c Digit) *Natural {
 }
 
 
-// Returns c = x*y div B, z = x*y mod B.
-func Mul1(x, y Digit) (Digit, Digit) {
-	// Split x and y into 2 sub-digits each (in base sqrt(B)),
-	// multiply the digits separately while avoiding overflow,
-	// and return the product as two separate digits.
-
-	const L0 = (L + 1)/2;
-	const L1 = L - L0;
-	const DL = L0 - L1;  // 0 or 1
-	const b  = 1<<L0;
-	const m  = b - 1;
-
-	// split x and y into sub-digits
-	// x = (x1*b + x0)
-	// y = (y1*b + y0)
-	x1, x0 := x>>L0, x&m;
-	y1, y0 := y>>L0, y&m;
-
-	// x*y = t2*b^2 + t1*b + t0
-	t0 := x0*y0;
-	t1 := x1*y0 + x0*y1;
-	t2 := x1*y1;
-
-	// compute the result digits but avoid overflow
-	// z = z1*B + z0 = x*y
-	z0 := (t1<<L0 + t0)&M;
-	z1 := t2<<DL + (t1 + t0>>L0)>>L1;
-	
-	return z1, z0;
-}
-
-
 func (x *Natural) Mul(y *Natural) *Natural {
  	n := len(x);
  	m := len(y);
-	z := new(Natural, n + m);
-
-	for j := 0; j < m; j++ {
-		d := y[j];
-		if d != 0 {
-			c := Digit(0);
-			for i := 0; i < n; i++ {
-				// z[i+j] += c + x[i]*d;
-				z1, z0 := Mul1(x[i], d);
-				c, z[i+j] = Split(c + z[i+j] + z0);
-				c += z1;
-			}
-			z[n+j] = c;
-		}
-	}
+	z := new(Natural, n + m);
+	Mul(z, x, y);
 
  	return Normalize(z);
 }
@@ -294,42 +419,26 @@ func (x *Natural) Pow(n uint) *Natural {
 }
 
 
-func Shl1(x, c Digit, s uint) (Digit, Digit) {
-	assert(s <= LogB);
-	return x >> (LogB - s), x << s & M | c
-}
-
-
-func Shr1(x, c Digit, s uint) (Digit, Digit) {
-	assert(s <= LogB);
-	return x << (LogB - s) & M, x >> s | c
-}
-
-
 func (x *Natural) Shl(s uint) *Natural {
-	n := len(x);
-	tsi := int(s / LogB);
-	s = s % LogB;
-	z := new(Natural, n + si + 1);
+	n := uint(len(x));
+	m := n + s/L;
+	z := new(Natural, m+1);
  	
-	c := Digit(0);
-	for i := 0; i < n; i++ { c, z[i+si] = Shl1(x[i], c, s); }
-	z[n+si] = c;
+	z[m] = Shl(z[m-n : m], x, s%L);
  	
  	return Normalize(z);
 }
 
 
 func (x *Natural) Shr(s uint) *Natural {
-	n := len(x);
-	tsi := int(s / LogB);
-	if si >= n { si = n; }
-	s = s % LogB;
-	assert(si <= n);
-	z := new(Natural, n - si);
+	n := uint(len(x));
+	m := n - s/L;
+	if m > n {  // check for underflow
+		m = 0;
+	}
+	z := new(Natural, m);
  	
-	c := Digit(0);
-	for i := n - 1; i >= si; i-- { c, z[i-si] = Shr1(x[i], c, s); }
+	Shr(z, x[n-m : n], s%L);
  	
  	return Normalize(z);
 }
@@ -337,68 +446,36 @@ func (x *Natural) Shr(s uint) *Natural {
 
 // DivMod needs multi-precision division which is not available if Digit
 // is already using the largest uint size. Split base before division,\n-// and merge again after. Each Digit is split into 3 Digit3\'s.\n+// and merge again after. Each Digit is split into 2 Digit2\'s.\n 
-func SplitBase(x *Natural) *[]Digit3 {
-	// TODO Use Log() for better result - don\'t need Normalize3 at the end!\n+func Unpack(x *Natural) *[]Digit2 {
+	// TODO Use Log() for better result - don\'t need Normalize2 at the end!\n  	n := len(x);
-	z := new([]Digit3, n*3 + 1);  // add space for extra digit (used by DivMod)\n-	for i, j := 0, 0; i < n; i, j = i+1, j+3 {\n+	z := new([]Digit2, n*2 + 1);  // add space for extra digit (used by DivMod)\n+	for i := 0; i < n; i++ {\n  		t := x[i];
-		z[j+0] = Digit3(t >> (L3*0) & M3);\n-		z[j+1] = Digit3(t >> (L3*1) & M3);\n-		z[j+2] = Digit3(t >> (L3*2) & M3);\n+		z[i*2] = Digit2(t & M2);\n+		z[i*2 + 1] = Digit2(t >> L2 & M2);\n  	}\n-	return Normalize3(z);\n+	return Normalize2(z);\
 }
 
 
-func MergeBase(x *[]Digit3) *Natural {
-	i := len(x);
-	j := (i+2)/3;\n-	z := new(Natural, j);\n-
-	switch i%3 {\n-	case 1: z[j-1] = Digit(x[i-1]); i--; j--;\n-	case 2: z[j-1] = Digit(x[i-1])<<L3 | Digit(x[i-2]); i -= 2; j--;\n-	case 0:\n+func Pack(x *[]Digit2) *Natural {
+	n := (len(x) + 1) / 2;\n+	z := new(Natural, n);\n+	if len(x) & 1 == 1 {\n+		// handle odd len(x)\n+		n--;\n+		z[n] = Digit(x[n*2]);\n  	}\n-	
-	for i >= 3 {\n-		z[j-1] = ((Digit(x[i-1])<<L3) | Digit(x[i-2]))<<L3 | Digit(x[i-3]);\n-		i -= 3;\n-		j--;\n+	for i := 0; i < n; i++ {\n+		z[i] = Digit(x[i*2 + 1]) << L2 | Digit(x[i*2]);\n  	}\n-	assert(j == 0);\n-
  	return Normalize(z);
 }
 
 
-func Split3(x Digit) (Digit, Digit3) {
-	return uint64(int64(x)>>L3), Digit3(x&M3)\n}\n\n\nfunc Product(x *[]Digit3, y Digit) {
-	n := len(x);\n-	c := Digit(0);\n-	for i := 0; i < n; i++ { c, x[i] = Split3(c + Digit(x[i])*y) }\n-	assert(c == 0);\n}\n\n\nfunc Quotient(x *[]Digit3, y Digit) {
-	n := len(x);\n-	c := Digit(0);\n-	for i := n-1; i >= 0; i-- {
-		t := c*B3 + Digit(x[i]);
-		c, x[i] = t%y, Digit3(t/y);\n-	}\n-	assert(c == 0);\n}\n\n\n // Division and modulo computation - destroys x and y. Based on the\n // algorithms described in:\n //\n@@ -413,8 +490,8 @@ func Quotient(x *[]Digit3, y Digit) {\n // is described in 1), while 2) provides the background for a more\n // accurate initial guess of the trial digit.\n \n-func DivMod(x, y *[]Digit3) (*[]Digit3, *[]Digit3) {\n-	const b = B3;\n+func DivMod2(x, y *[]Digit2) (*[]Digit2, *[]Digit2) {\n+	const b = B2;\n  	\n  	n := len(x);\n  	m := len(y);\n@@ -424,14 +501,9 @@ func DivMod(x, y *[]Digit3) (*[]Digit3, *[]Digit3) {\n  	\n  	if m == 1 {\n  		// division by single digit\n-	\t\td := Digit(y[0]);
-	\t\tc := Digit(0);
-	\t\tfor i := n; i > 0; i-- {
-	\t\t\tt := c*b + Digit(x[i-1]);
-	\t\t\tc, x[i] = t%d, Digit3(t/d);
-	\t\t}
-	\t\tx[0] = Digit3(c);
-
+\t\t// result is shifted left by 1 in place!\n+\t\tx[0] = Div1(x[1 : n+1], x[0 : n], y[0]);\n+\t\t\n  	} else if m > n {\n  		// quotient = 0, remainder = x\n  		// TODO in this case we shouldn\'t even split base - FIX THIS\n@@ -444,24 +516,34 @@ func DivMod(x, y *[]Digit3) (*[]Digit3, *[]Digit3) {\n  		\n  		// normalize x and y\n  		f := b/(Digit(y[m-1]) + 1);\n-		Product(x, f);\n-		Product(y, f);\n+		Mul1(x, x, Digit2(f));\n+		Mul1(y, y, Digit2(f));\n  		assert(b/2 <= y[m-1] && y[m-1] < b);  // incorrect scaling\n  		\n-		d2 := Digit(y[m-1])*b + Digit(y[m-2]);\n+		y1, y2 := Digit(y[m-1]), Digit(y[m-2]);\n+		d2 := Digit(y1)*b + Digit(y2);\n  		for i := n-m; i >= 0; i-- {\n  		\tk := i+m;\n  		\t\n  		\t// compute trial digit\n-		\t\tr3 := (Digit(x[k])*b + Digit(x[k-1]))*b + Digit(x[k-2]);\n-		\t\tq := r3/d2;\n-		\t\tif q >= b { q = b-1 }\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\tif x0 != y1 {\n+\t\t\t\t\tq = (x0*b + x1)/y1;\n+\t\t\t\t} else {\n+\t\t\t\t\tq = b-1;\n+\t\t\t\t}\n+\t\t\t\tfor y2 * q > (x0*b + x1 - y1*q)*b + x2 {\n+\t\t\t\t\tq--\n+\t\t\t\t}\n+\t\t\t}\n  		\t\n  		\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\t\tc, x[i+j] = Digit(int64(t)>>L3), Digit3(t&M3);\n+		\t\t\t\tc, x[i+j] = Digit(int64(t)>>L2), Digit2(t&M2);\n  		\t\t}\n  		\t\t\n  		\t\t// correct if trial digit was too large\n@@ -469,18 +551,20 @@ func DivMod(x, y *[]Digit3) (*[]Digit3, *[]Digit3) {\n  		\t\t\t// add y\n  		\t\t\tc := Digit(0);\n  		\t\t\tfor j := 0; j < m; j++ {\n-		\t\t\t\tc, x[i+j] = Split3(c + Digit(x[i+j]) + Digit(y[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}\n  		\t\t\tassert(c + Digit(x[k]) == 0);\n  		\t\t\t// correct trial digit\n  		\t\t\tq--;\n  		\t\t}\n  		\t\t\n-		\t\tx[k] = Digit3(q);\n+		\t\tx[k] = Digit2(q);\n  		\t}\n  		\t\n  		\t// undo normalization for remainder\n-		\tQuotient(x[0 : m], f);\n+		\tc := Div1(x[0 : m], x[0 : m], Digit2(f));\n+		\tassert(c == 0);\n  		}\n  	}\n  \n  	return x[m : n+1], x[0 : m];\n@@ -488,14 +572,20 @@ func DivMod(x, y *[]Digit3) (*[]Digit3, *[]Digit3) {\n \n \n func (x *Natural) Div(y *Natural) *Natural {\n-	q, r := DivMod(SplitBase(x), SplitBase(y));\n-	return MergeBase(q);\n+	q, r := DivMod2(Unpack(x), Unpack(y));\n+	return Pack(q);\n }\n \n \n func (x *Natural) Mod(y *Natural) *Natural {\
-	q, r := DivMod(SplitBase(x), SplitBase(y));\n-	return MergeBase(r);\n+	q, r := DivMod2(Unpack(x), Unpack(y));\n+	return Pack(r);\n+}\n+\n+\n+func (x *Natural) DivMod(y *Natural) (*Natural, *Natural) {\n+	q, r := DivMod2(Unpack(x), Unpack(y));\n+	return Pack(q), Pack(r);\n }\n \n \n@@ -520,17 +610,17 @@ func (x *Natural) Cmp(y *Natural) int {\n }\n \n \n-func Log1(x Digit) int {\n+func Log2(x Digit) int {\n  	n := -1;\n  	for x != 0 { x = x >> 1; n++; }  // BUG >>= broken for uint64\n  	return n;\n }\n \n \n-func (x *Natural) Log() int {\
+func (x *Natural) Log2() int {\
  	n := len(x);\n  	if n > 0 {\n-	\t\tn = (n - 1)*L + Log1(x[n - 1]);\n+	\t\tn = (n - 1)*L + Log2(x[n - 1]);\n  	} else {\n  	\t\tn = -1;\n  	}\n@@ -544,11 +634,10 @@ func (x *Natural) And(y *Natural) *Natural {\
  	if n < m {\
  		return y.And(x);
  	}
-	assert(n >= m);
-	z := new(Natural, n);\
-
-	for i := 0; i < m; i++ { z[i] = x[i] & y[i]; }\
-	for i := m; i < n; i++ { z[i] = x[i]; }\
+	z := new(Natural, n);\
+	And(z[0 : m], x[0 : m], y);\
+	Or1(z[m : n], x[m : n], 0);\
 
  	return Normalize(z);
 }
@@ -560,11 +649,10 @@ func (x *Natural) Or(y *Natural) *Natural {\
  	if n < m {\
  		return y.Or(x);
  	}
-	assert(n >= m);
-	z := new(Natural, n);\
-
-	for i := 0; i < m; i++ { z[i] = x[i] | y[i]; }\
-	for i := m; i < n; i++ { z[i] = x[i]; }\
+	z := new(Natural, n);\
+	Or(z[0 : m], x[0 : m], y);\
+	Or1(z[m : n], x[m : n], 0);\
 
  	return Normalize(z);
 }
@@ -576,24 +664,15 @@ func (x *Natural) Xor(y *Natural) *Natural {\
  	if n < m {\
  		return y.Xor(x);
  	}
-	assert(n >= m);
-	z := new(Natural, n);\
-
-	for i := 0; i < m; i++ { z[i] = x[i] ^ y[i]; }\
-	for i := m; i < n; i++ { z[i] = x[i]; }\
+	z := new(Natural, n);\
+	Xor(z[0 : m], x[0 : m], y);\
+	Or1(z[m : n], x[m : n], 0);\
 
  	return Normalize(z);
 }
 
 
-func Copy(x *Natural) *Natural {
-	z := new(Natural, len(x));
-	//*z = *x;  // BUG assignment does\'t work yet
-	for i := len(x) - 1; i >= 0; i-- { z[i] = x[i]; }\
-	return z;
-}
-
-
 // Computes x = x div d (in place - the recv maybe modified) for "small" d\'s.\
 // Returns updated x and x mod d.\n func (x *Natural) DivMod1(d Digit) (*Natural, Digit) {\
  	n := len(x);\
@@ -601,36 +680,36 @@ func (x *Natural) DivMod1(d Digit) (*Natural, Digit) {\
  \tc := Digit(0);\
  	for i := len(x) - 1; i >= 0; i-- {\
-\t\t\tc = c<<L + x[i];\n-\t\t\tx[i] = c/d;\n-\t\t\tc %= d;\n+\t\t\tt := c<<L + x[i];\n+\t\t\tc, x[i] = t%d, t/d;\n  	}\
  \n  	return Normalize(x), c;\
  }\
  \n  \n-func (x *Natural) String(base Digit) string {\
+func (x *Natural) String(base uint) string {\
  	if x.IsZero() {\
  		return "0";
  	}\
  	\t
  	// allocate string\n-\t// TODO n is too small for bases < 10!!!\n-\tassert(base >= 10);  // for now\n-\t// approx. length: 1 char for 3 bits\n-\tn := x.Log()/3 + 10;  // +10 (round up) - what is the right number?\n+\tassert(2 <= base && base <= 16);\n+\tn := (x.Log2() + 1) / Log2(Digit(base)) + 1;  // TODO why the +1?\n  \ts := new([]byte, n);\
  \n  	// convert\n-\tconst hex = "0123456789abcdef";\n+\n+\t// don\'t destroy x, make a copy\n+\tt := new(Natural, len(x));\n+\tOr1(t, x, 0);  // copy x\n+\t\n  \ti := n;\
-\t\tx = Copy(x);  // don\'t destroy recv\n-\tfor !x.IsZero() {\n+\tfor !t.IsZero() {\
  \t\ti--;\
  \t\tvar d Digit;\
-\t\t\tx, d = x.DivMod1(base);\n-\t\t\ts[i] = hex[d];\n+\t\t\tt, d = t.DivMod1(Digit(base));\n+\t\t\ts[i] = "0123456789abcdef"[d];\n  \t};\
  \n  	return string(s[i : n]);
  }\
@@ -665,24 +744,24 @@ func (x *Natural) Gcd(y *Natural) *Natural {\
 }\n \n \n-func HexValue(ch byte) Digit {\
-\td := Digit(1 << LogH);\n+func HexValue(ch byte) uint {\
+\td := uint(1 << LogH);\
  \tswitch {\
-\tcase \'0\' <= ch && ch <= \'9\': d = Digit(ch - \'0\');\n-\tcase \'a\' <= ch && ch <= \'f\': d = Digit(ch - \'a\') + 10;\n-\tcase \'A\' <= ch && ch <= \'F\': d = Digit(ch - \'A\') + 10;\n+\tcase \'0\' <= ch && ch <= \'9\': d = uint(ch - \'0\');\n+\tcase \'a\' <= ch && ch <= \'f\': d = uint(ch - \'a\') + 10;\n+\tcase \'A\' <= ch && ch <= \'F\': d = uint(ch - \'A\') + 10;\
  \t}\
  \treturn d;\
  }\
  \n  \n  // TODO auto-detect base if base argument is 0\n-export func NatFromString(s string, base Digit) *Natural {\
+export func NatFromString(s string, base uint) *Natural {\
  \tx := NatZero;\
  \tfor i := 0; i < len(s); i++ {\
  \t\td := HexValue(s[i]);\
  \t\tif d < base {\
-\t\t\t\tx = x.MulAdd1(base, d);\n+\t\t\t\tx = x.MulAdd1(Digit(base), Digit(d));\n  \t\t} else {\
  \t\t\tbreak;\
  \t\t}\
  \t}\
@@ -776,19 +855,31 @@ func (x *Integer) Mul(y *Integer) *Integer {\
 \n \n func (x *Integer) Quo(y *Integer) *Integer {\
-\tpanic("UNIMPLEMENTED");\n-\treturn nil;\n+\t// x / y == x / y\n+\t// x / (-y) == -(x / y)\n+\t// (-x) / y == -(x / y)\n+\t// (-x) / (-y) == x / y\n+\treturn &Integer{x.sign != y.sign, x.mant.Div(y.mant)};\n }\n \n \n func (x *Integer) Rem(y *Integer) *Integer {\
-\tpanic("UNIMPLEMENTED");\n-\treturn nil;\n+\t// x % y == x % y\n+\t// x % (-y) == x % y\n+\t// (-x) % y == -(x % y)\n+\t// (-x) % (-y) == -(x % y)\n+\treturn &Integer{y.sign, x.mant.Mod(y.mant)};\n+}\n+\n+\n+func (x *Integer) QuoRem(y *Integer) (*Integer, *Integer) {\
+\tq, r := x.mant.DivMod(y.mant);\n+\treturn &Integer{x.sign != y.sign, q}, &Integer{y.sign, q};\n }\n \n \n func (x *Integer) Div(y *Integer) *Integer {\
-\tpanic("UNIMPLEMENTED");\n+\tq, r := x.mant.DivMod(y.mant);\n  \treturn nil;\
  }\n \n @@ -805,7 +896,7 @@ func (x *Integer) Cmp(y *Integer) int {\
 }\n \n \n-func (x *Integer) String(base Digit) string {\
+func (x *Integer) String(base uint) string {\
  	if x.mant.IsZero() {\
  		return "0";
  	}\
@@ -817,7 +908,7 @@ func (x *Integer) String(base Digit) string {\
 }\n \n \t\n-export func IntFromString(s string, base Digit) *Integer {\
+export func IntFromString(s string, base uint) *Integer {\
  \t// get sign, if any\n  \tsign := false;\
  \tif len(s) > 0 && (s[0] == \'-\' || s[0] == \'+\') {\ndiff --git a/usr/gri/bignum/bignum_test.go b/usr/gri/bignum/bignum_test.go
index ae81036b0f..10d8da7db7 100644
--- a/usr/gri/bignum/bignum_test.go
+++ b/usr/gri/bignum/bignum_test.go
@@ -38,14 +38,54 @@ func TEST_EQ(n uint, x, y *Big.Natural) {
 }\n \n \n+func TestLog2() {\
+\ttest_msg = "TestLog2A";\n+\tTEST(0, Big.Nat(1).Log2() == 0);\n+\tTEST(1, Big.Nat(2).Log2() == 1);\n+\tTEST(2, Big.Nat(3).Log2() == 1);\n+\tTEST(3, Big.Nat(4).Log2() == 2);\n+\t\n+\ttest_msg = "TestLog2B";\n+\tfor i := uint(0); i < 100; i++ {\n+\t\tTEST(i, Big.Nat(1).Shl(i).Log2() == int(i));\n+\t}\n+}\n+\n+\n func TestConv() {\
-\ttest_msg = "TestConv";\n+\ttest_msg = "TestConvA";\n  \tTEST(0, a.Cmp(Big.Nat(991)) == 0);\
  \tTEST(1, b.Cmp(Big.Fact(20)) == 0);\
  \tTEST(2, c.Cmp(Big.Fact(100)) == 0);\
  \tTEST(3, a.String(10) == sa);\
  \tTEST(4, b.String(10) == sb);\
  \tTEST(5, c.String(10) == sc);\
+\n+\ttest_msg = "TestConvB";\n+\tt := c.Mul(c);\n+\tfor base := uint(2); base <= 16; base++ {\n+\t\tTEST_EQ(base, Big.NatFromString(t.String(base), base), t);\n+\t}\n+}\n+\n+\n+func Sum(n uint, scale *Big.Natural) *Big.Natural {\
+\ts := Big.Nat(0);\n+\tfor ; n > 0; n-- {\n+\t\ts = s.Add(Big.Nat(uint64(n)).Mul(scale));\n+\t}\n+\treturn s;\n+}\n+\n+\n+func TestAdd() {\
+\ttest_msg = "TestAddA";\n+\n+\ttest_msg = "TestAddB";\n+\tfor i := uint(0); i < 100; i++ {\n+\t\tt := Big.Nat(uint64(i));\n+\t\tTEST_EQ(i, Sum(i, c), t.Mul(t).Add(t).Shr(1).Mul(c));\n+\t}\n }\n \n \n@@ -163,7 +203,9 @@ func TestPop() {\
 \n \n func main() {\
+\tTestLog2();\
  \tTestConv();\
+\tTestAdd();\
  \tTestShift();\
  \tTestMul();\
  \tTestDiv();

変更の背景

このコミットの主な背景は、多倍長整数の除算および剰余演算のパフォーマンス改善です。既存の実装では、内部的に20ビットの単位で演算を行っていましたが、これを30ビットの単位に変更することで、より効率的な計算が可能になります。これは、一度に処理できるデータの量が増えるため、ループ回数が減り、全体的な処理速度が向上するためです。

また、コードのリファクタリングとクリーンアップも行われています。これは、コードの可読性、保守性、および将来的な拡張性を向上させるための一般的なソフトウェア開発プラクティスです。新しいテストの追加は、変更が正しく機能することを確認し、将来の変更に対する回帰を防ぐための重要なステップです。最終的には、このbignumパッケージをGo言語の標準ライブラリとしてチェックインするための準備の一環として行われました。

前提知識の解説

任意精度算術(多倍長整数)

任意精度算術(Arbitrary-precision arithmetic)は、コンピュータの固定長ワードサイズに制限されずに、任意の桁数の整数を扱うための技術です。通常のintlong型では扱えない非常に大きな数(例えば、数千桁や数万桁の数)を計算する必要がある場合に用いられます。

多倍長整数は通常、数値の各「桁」(Digit)を配列やスライスに格納することで表現されます。この「桁」の基数(radix)は、コンピュータのワードサイズ(例えば32ビットや64ビット)に合わせて選ばれることが多く、これにより効率的な演算が可能になります。

除算・剰余アルゴリズム(KnuthのアルゴリズムD)

多倍長整数の除算は、通常の小学校で習う筆算の除算と似た原理で行われます。最も広く知られている効率的なアルゴリズムの一つに、ドナルド・クヌース(Donald Knuth)の著書『The Art of Computer Programming』で解説されている「アルゴリズムD」(Algorithm D)があります。

アルゴリズムDは、以下のようなステップで構成されます。

  1. 正規化(Normalization): 除数(y)の最上位桁が特定の範囲(例えば、基数の半分以上)になるように、被除数(x)と除数を同じ量だけ左シフトします。これにより、試行商(trial quotient)の推定精度が向上します。
  2. 試行商の推定(Trial Quotient Estimation): 被除数と除数の最上位の数桁を使って、現在の桁の商を推定します。この推定は、オーバーフローを避けるために慎重に行われます。
  3. 乗算と減算(Multiply and Subtract): 推定された試行商と除数を乗算し、その結果を被除数から減算します。
  4. 修正(Correction): 減算の結果が負になった場合(試行商が大きすぎた場合)、試行商をデクリメントし、除数を加算して修正します。
  5. 繰り返し: これを被除数のすべての桁に対して繰り返します。

このアルゴリズムの効率は、試行商の推定精度に大きく依存します。推定が正確であればあるほど、修正ステップの回数が減り、全体的なパフォーマンスが向上します。

基数(Radix)とビット幅の選択

多倍長整数演算において、内部的に数値を表現する「桁」の基数(またはビット幅)の選択は非常に重要です。

  • 小さいビット幅(例: 20ビット): 各桁が表現できる値の範囲が狭いため、配列の要素数が多くなり、ループの回数が増えます。これにより、演算のオーバーヘッドが大きくなる可能性があります。しかし、中間結果のオーバーフローを管理しやすくなるという利点もあります。
  • 大きいビット幅(例: 30ビット、60ビット): 各桁が表現できる値の範囲が広いため、配列の要素数が減り、ループの回数が少なくなります。これにより、演算の効率が向上し、特に乗算や除算のような複雑な演算で顕著なパフォーマンス改善が見られます。ただし、中間結果のオーバーフローをより注意深く管理する必要があります。

このコミットでは、uint64Digit型として使用しているため、最大で64ビットの値を扱うことができます。20ビットから30ビットへの変更は、uint64の半分(32ビット)に近い値を使用することで、CPUのワードサイズをより効率的に利用し、かつ中間演算でのオーバーフローを管理しやすいバランスの取れた選択と言えます。

技術的詳細

このコミットの技術的な核心は、多倍長整数の除算および剰余演算の基数を変更し、それに伴うアルゴリズムの最適化とコードのリファクタリングです。

1. 内部基数の変更 (L3 から L2 へ)

以前のバージョンでは、Digituint64)を3つのDigit3uint32)に分割して除算を行っていました。これは、LogB / 3で定義されるL3(約20ビット)を基数としていました。

// 変更前
const (
	L3 = LogB / 3; // 約20ビット
	B3 = 1 << L3;
	M3 = B3 - 1;

	L2 = L3 * 2; // 約40ビット
	B2 = 1 << L2;
	M2 = B2 - 1;

	L = L3 * 3; // 約60ビット
	B = 1 << L;
	M = B - 1;
)
type (
	Digit3 uint32;
	Digit  uint64;
)

このコミットでは、Digitを2つのDigit2uint32)に分割するように変更されました。これにより、基数がLogB / 2で定義されるL2(約30ビット)に変更されました。

// 変更後
const (
	L2 = LogB / 2; // 約30ビット
	B2 = 1 << L2;
	M2 = B2 - 1;

	L = L2 * 2; // 約60ビット
	B = 1 << L;
	M = B - 1;
)
type (
	Digit2 uint32; // 新しい内部演算単位
	Digit  uint64;
)

LogBLogW - LogHで定義されており、LogWDigitのビット幅(64)、LogHは不明ですが、おそらくオーバーヘッドを考慮した定数です。LDigitの有効ビット幅を表し、L2 * 2で約60ビットとなります。

この変更により、除算アルゴリズムが一度に処理する「桁」のビット幅が20ビットから30ビットに増加しました。これにより、同じ数値範囲を表現するために必要な内部的な桁数が減少し、ループの反復回数が減るため、特に除算のような計算量の多い操作でパフォーマンスが向上します。

2. 除算アルゴリズムの最適化 (DivMod から DivMod2 へ)

主要な除算関数がDivModからDivMod2にリネームされ、内部的なDigit3からDigit2への変更が反映されました。

  • UnpackPack関数の導入: 以前のSplitBaseMergeBaseに代わり、Natural型(多倍長整数)を[]Digit2スライスに展開(Unpack)し、計算後に再びNatural型にパック(Pack)する関数が導入されました。これにより、Natural型と内部演算用の[]Digit2スライス間の変換がより明確になりました。
  • 単一桁除算の最適化: DivMod2内のm == 1(除数が単一桁の場合)の処理が、新しいDiv1ヘルパー関数を使用するように変更されました。Div1は、単一桁除算を効率的に行うための新しい関数です。
  • KnuthのアルゴリズムDの適用: DivMod2の主要な部分では、KnuthのアルゴリズムDがより明示的に適用されています。特に、試行商qの推定ロジックが改善され、より正確な初期推定を行うことで、修正ステップ(q--)の回数を減らす試みがなされています。
    • 以前はr3/d2という単純な計算でしたが、新しいコードではy1y2(除数の上位2桁)とx0, x1, x2(被除数の上位3桁)を用いたKnuthの推奨する試行商推定ロジックが導入されています。これにより、試行商の精度が向上し、ループ内の修正回数が減少します。
  • Mul1関数の導入: Mul1[]Digit2スライスと単一のDigit2を乗算するための新しいヘルパー関数です。これは、正規化ステップでxyをスケーリングするために使用されます。

3. 新しい「Raw Operations」の追加

bignum.goには、Digitスライスに対する基本的なビット演算(And1, And, Or1, Or, Xor1, Xor)と、加算(Add1, Add)、減算(Sub1, Sub)、乗算(Mul11, Mul)、除算(Div1)、シフト(Shl, Shr)の「Raw Operations」が多数追加されました。

これらの関数は、多倍長整数の内部表現である[]Digitまたは[]Digit2スライスに対して直接操作を行う低レベルのヘルパー関数です。これにより、高レベルのNatural型メソッド(Add, Sub, Mul, Shl, Shrなど)が、これらの効率的な低レベル操作を呼び出すようにリファクタリングされ、コードのモジュール性と再利用性が向上しました。

特に注目すべきはMul11関数です。これは2つのDigitを乗算し、その結果を2つのDigit(上位と下位)として返す関数です。これは、多倍長乗算の基本的な構成要素であり、Karatsuba法のようなより高度な乗算アルゴリズムの基盤となります。

4. その他のリファクタリングとテストの追加

  • Natural型のAdd, Sub, Mul, Shl, Shrメソッドが、新しく追加された「Raw Operations」を使用するように書き換えられました。これにより、コードがより簡潔になり、重複が排除されました。
  • Integer型(符号付き多倍長整数)のQuo, Rem, QuoRem, Divメソッドが実装されました。以前はpanic("UNIMPLEMENTED")となっていましたが、Natural型のDivModを利用してこれらの演算が提供されるようになりました。
  • Log1Log2にリネームされ、Natural型のLogメソッドもLog2を使用するように変更されました。これは、おそらく対数の底が2であることをより明確にするためです。
  • String変換関数やNatFromStringIntFromString関数も、Digit型からuint型への変更が反映され、より汎用的に使用できるようになりました。
  • bignum_test.goには、TestLog2, TestConvB, TestAddなど、多数の新しいテストケースが追加されました。これにより、変更されたロジックが正しく機能すること、特に新しい基数での演算が期待通りに行われることが検証されます。

これらの変更は、bignumパッケージのパフォーマンスと堅牢性を大幅に向上させ、Go言語の標準ライブラリとしての採用に向けた重要なステップとなりました。

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

bignum.go

  1. 定数と型の変更:

    --- a/usr/gri/bignum/bignum.go
    +++ b/usr/gri/bignum/bignum.go
    @@ -38,22 +38,18 @@ const LogB = LogW - LogH;
     
     
     const (
    -	L3 = LogB / 3;
    -	B3 = 1 << L3;
    -	M3 = B3 - 1;
    -
    -	L2 = L3 * 2;
    +	L2 = LogB / 2;
      	B2 = 1 << L2;
      	M2 = B2 - 1;
      
    -	L = L3 * 3;
    +	L = L2 * 2;
      	B = 1 << L;
      	M = B - 1;
     )
     
     
     type (
    -	Digit3 uint32;
    +	Digit2 uint32;
      	Digit  uint64;
     )
    

    L3, B3, M3が削除され、L2, B2, M2の定義がLogB / 2に基づくように変更されました。Digit3型がDigit2型に置き換えられました。

  2. 新しいRaw Operationsの追加: And1, And, Or1, Or, Xor1, Xor, Add1, Add, Sub1, Sub, Mul11, Mul, Mul1, Div1, Shl, Shrといった多数の低レベル演算関数が追加されました。これらは[]Digitまたは[]Digit2スライスを直接操作します。

  3. Mul11関数の追加:

    // 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)),
    	// multiply the digits separately while avoiding overflow,
    	// and return the product as two separate digits.
    
    	const L0 = (L + 1)/2;
    	const L1 = L - L0;
    	const DL = L0 - L1;  // 0 or 1
    	const b  = 1<<L0;
    	const m  = b - 1;
    
    	// split x and y into sub-digits
    	// x = (x1*b + x0)
    	// y = (y1*b + y0)
    	x1, x0 := x>>L0, x&m;
    	y1, y0 := y>>L0, y&m;
    
    	// x*y = t2*b^2 + t1*b + t0
    	t0 := x0*y0;
    	t1 := x1*y0 + x0*y1;
    	t2 := x1*y1;
    
    	// compute the result digits but avoid overflow
    	// z = z1*B + z0 = x*y
    	z0 := (t1<<L0 + t0)&M;
    	z1 := t2<<DL + (t1 + t0>>L0)>>L1;
    	
    	return z1, z0;
    }
    

    2つのDigitを乗算し、結果を2つのDigit(上位と下位)として返す関数。多倍長乗算の基盤となります。

  4. UnpackPack関数の導入:

    --- a/usr/gri/bignum/bignum.go
    +++ b/usr/gri/bignum/bignum.go
    @@ -337,68 +446,36 @@ func (x *Natural) Shr(s uint) *Natural {
     
     // DivMod needs multi-precision division which is not available if Digit
     // is already using the largest uint size. Split base before division,\n-// and merge again after. Each Digit is split into 3 Digit3\'s.\n+// and merge again after. Each Digit is split into 2 Digit2\'s.\n 
    -func SplitBase(x *Natural) *[]Digit3 {
    -	// TODO Use Log() for better result - don\'t need Normalize3 at the end!\n+func Unpack(x *Natural) *[]Digit2 {
    +	// TODO Use Log() for better result - don\'t need Normalize2 at the end!\n      	n := len(x);
    -	z := new([]Digit3, n*3 + 1);  // add space for extra digit (used by DivMod)\n-	for i, j := 0, 0; i < n; i, j = i+1, j+3 {\n+	z := new([]Digit2, n*2 + 1);  // add space for extra digit (used by DivMod)\n+	for i := 0; i < n; i++ {\
      		t := x[i];
    -		z[j+0] = Digit3(t >> (L3*0) & M3);\n-		z[j+1] = Digit3(t >> (L3*1) & M3);\n-		z[j+2] = Digit3(t >> (L3*2) & M3);\n+		z[i*2] = Digit2(t & M2);\n+		z[i*2 + 1] = Digit2(t >> L2 & M2);\
      	}\n    -	return Normalize3(z);\n+	return Normalize2(z);\
     }
     
     
    -func MergeBase(x *[]Digit3) *Natural {
    -	i := len(x);
    -	j := (i+2)/3;\n-	z := new(Natural, j);\n-
    -	switch i%3 {\n-	case 1: z[j-1] = Digit(x[i-1]); i--; j--;\n-	case 2: z[j-1] = Digit(x[i-1])<<L3 | Digit(x[i-2]); i -= 2; j--;\n-	case 0:\n+func Pack(x *[]Digit2) *Natural {
    +	n := (len(x) + 1) / 2;\n+	z := new(Natural, n);\n+	if len(x) & 1 == 1 {\n+		// handle odd len(x)\n+		n--;\n+		z[n] = Digit(x[n*2]);\n      	}\n    -	
    -	for i >= 3 {\n    -		z[j-1] = ((Digit(x[i-1])<<L3) | Digit(x[i-2]))<<L3 | Digit(x[i-3]);\n    -		i -= 3;\n    -		j--;\n+	for i := 0; i < n; i++ {\
    +		z[i] = Digit(x[i*2 + 1]) << L2 | Digit(x[i*2]);\
      	}\n    -	assert(j == 0);\n    -
      	return Normalize(z);
     }
    

    Natural型と[]Digit2スライス間の変換を行うUnpackPack関数が追加されました。

  5. DivMod2関数の変更:

    --- a/usr/gri/bignum/bignum.go
    +++ b/usr/gri/bignum/bignum.go
    @@ -413,8 +490,8 @@ func Quotient(x *[]Digit3, y Digit) {\
     // is described in 1), while 2) provides the background for a more
     // accurate initial guess of the trial digit.\n 
    -func DivMod(x, y *[]Digit3) (*[]Digit3, *[]Digit3) {\
    -	const b = B3;\n+func DivMod2(x, y *[]Digit2) (*[]Digit2, *[]Digit2) {\
    +	const b = B2;\n      	\n      	n := len(x);\
      	m := len(y);\
    @@ -424,14 +501,9 @@ func DivMod(x, y *[]Digit3) (*[]Digit3, *[]Digit3) {\
      	\n      	if m == 1 {\
      		// division by single digit\n-	\t\td := Digit(y[0]);\n-	\t\tc := Digit(0);\n-	\t\tfor i := n; i > 0; i-- {\
    -	\t\t\tt := c*b + Digit(x[i-1]);\n-	\t\t\tc, x[i] = t%d, Digit3(t/d);\n-	\t\t}\n-	\t\tx[0] = Digit3(c);
    -
    +\t\t// result is shifted left by 1 in place!\n+\t\tx[0] = Div1(x[1 : n+1], x[0 : n], y[0]);\n+\t\t\n      	} else if m > n {\
      		// quotient = 0, remainder = x\n      		// TODO in this case we shouldn\'t even split base - FIX THIS\n    @@ -444,24 +516,34 @@ func DivMod(x, y *[]Digit3) (*[]Digit3, *[]Digit3) {\
      		\n      		// normalize x and y\n      		f := b/(Digit(y[m-1]) + 1);\n    -		Product(x, f);\n    -		Product(y, f);\n+		Mul1(x, x, Digit2(f));\n+		Mul1(y, y, Digit2(f));\n      		assert(b/2 <= y[m-1] && y[m-1] < b);  // incorrect scaling\n      		\n    -		d2 := Digit(y[m-1])*b + Digit(y[m-2]);\n+		y1, y2 := Digit(y[m-1]), Digit(y[m-2]);\n+		d2 := Digit(y1)*b + Digit(y2);\n      		for i := n-m; i >= 0; i-- {\
      		\tk := i+m;\
      		\t\n      		\t// compute trial digit\n    -		\t\tr3 := (Digit(x[k])*b + Digit(x[k-1]))*b + Digit(x[k-2]);\n    -		\t\tq := r3/d2;\n    -		\t\tif q >= b { q = b-1 }\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\tif x0 != y1 {\n+\t\t\t\t\tq = (x0*b + x1)/y1;\n+\t\t\t\t} else {\n+\t\t\t\t\tq = b-1;\n+\t\t\t\t}\n+\t\t\t\tfor y2 * q > (x0*b + x1 - y1*q)*b + x2 {\n+\t\t\t\t\tq--\n+\t\t\t\t}\n+\t\t\t}\n      		\t\n      		\t// subtract y*q\n      		\t\tc := Digit(0);\n      		\t\tfor j := 0; j < m; j++ {\
      		\t\t\tt := c + Digit(x[i+j]) - Digit(y[j])*q;  // arithmetic shift!\n    -		\t\t\t\tc, x[i+j] = Digit(int64(t)>>L3), Digit3(t&M3);\n+\t\t\t\t\tc, x[i+j] = Digit(int64(t)>>L2), Digit2(t&M2);\n      		\t\t}\n      		\t\t\n      		\t// correct if trial digit was too large\n    @@ -469,18 +551,20 @@ func DivMod(x, y *[]Digit3) (*[]Digit3, *[]Digit3) {\
      		\t\t\t// add y\n      		\t\t\tc := Digit(0);\n      		\t\t\tfor j := 0; j < m; j++ {\
    -		\t\t\t\tc, x[i+j] = Split3(c + Digit(x[i+j]) + Digit(y[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}\n      		\t\t\tassert(c + Digit(x[k]) == 0);\n      		\t\t\t// correct trial digit\n      		\t\t\tq--;\n      		\t\t}\n      		\t\t\n    -		\t\tx[k] = Digit3(q);\n+\t\t\t\tx[k] = Digit2(q);\n      		\t}\n      		\t\n      		\t// undo normalization for remainder\n    -		\tQuotient(x[0 : m], f);\n+\t\t\tc := Div1(x[0 : m], x[0 : m], Digit2(f));\n+\t\t\tassert(c == 0);\n      		}\n      	}\n      \n      	return x[m : n+1], x[0 : m];
    

    DivModDivMod2にリネームされ、内部的にDigit2を使用するように変更されました。単一桁除算のロジックがDiv1を使用するように簡略化され、KnuthのアルゴリズムDに基づく試行商の推定ロジックがより詳細に実装されました。

  6. Natural型のメソッドの変更: Add, Sub, Mul, Shl, Shrなどのメソッドが、新しく追加されたRaw Operationsを呼び出すように変更されました。

  7. Integer型の除算・剰余演算の実装: Quo, Rem, QuoRem, DivメソッドがNatural型のDivModを利用して実装されました。

bignum_test.go

  1. TestLog2の追加: Log2関数のテストが追加されました。

  2. TestConvBの追加: StringNatFromStringの相互変換のテストが追加され、異なる基数(2から16)での変換が正しく行われることを検証しています。

  3. TestAddの追加: 加算演算のテストが追加されました。

コアとなるコードの解説

定数と型の変更

  • const ( L2 = LogB / 2; B2 = 1 << L2; M2 = B2 - 1; ): これは、多倍長整数演算の内部的な基数を定義しています。LogBDigit型(uint64)の有効ビット幅に関連する定数です。L2は、Digitを2つの部分に分割する際の各部分のビット幅を約30ビットに設定しています。B2はその基数(2^L2)、M2はマスク(B2 - 1)です。これにより、以前の20ビット単位の演算から30ビット単位の演算へと移行し、一度に処理できるデータ量が増えることで、特に除算の効率が向上します。
  • type ( Digit2 uint32; Digit uint64; ): Digit2型が新しく導入され、内部的な演算単位としてuint32が使用されます。Digit型は引き続きuint64で、多倍長整数の各「桁」を表現します。

Raw Operations

新しく追加されたAnd1, And, Or1, Or, Xor1, Xor, Add1, Add, Sub1, Sub, Mul11, Mul, Mul1, Div1, Shl, Shrなどの関数は、多倍長整数の内部表現である[]Digitまたは[]Digit2スライスに対して直接作用する低レベルのプリミティブ演算です。

  • Mul11(x, y Digit) (Digit, Digit): この関数は、2つのDigituint64)を乗算し、その結果を2つのDigit(上位と下位)として返します。uint64 * uint64の結果は最大で128ビットになる可能性があるため、これを2つのuint64で表現します。内部では、DigitをさらにL0ビットとL1ビットの2つのサブ桁に分割し、それぞれのサブ桁を乗算して結果を合成するという「分割統治」的な手法を用いています。これは、多倍長乗算の基本的なビルディングブロックであり、より大きな数の乗算(Mul関数)で利用されます。
  • Mul1(z, x *[]Digit2, y Digit2) Digit2: []Digit2スライスで表される多倍長整数xと単一のDigit2``yを乗算し、結果をzに格納します。これは、除算の正規化ステップで被除数と除数をスケーリングするために使用されます。
  • Div1(z, x *[]Digit2, y Digit2) Digit2: []Digit2スライスで表される多倍長整数xを単一のDigit2``yで除算し、結果をzに格納します。これは、除数が単一桁の場合の除算を効率的に処理するために使用されます。

UnpackPack関数

  • Unpack(x *Natural) *[]Digit2: Natural型(多倍長整数)を、内部演算に適した[]Digit2スライスに変換します。Naturalの各Digituint64)を、2つのDigit2uint32)に分割してスライスに格納します。
  • Pack(x *[]Digit2) *Natural: []Digit2スライスで表された内部演算結果を、再びNatural型に変換します。[]Digit2の2つの要素を組み合わせて1つのDigituint64)を形成します。

これらの関数は、Natural型とDivMod2のような内部演算関数との間のインターフェースとして機能し、データ表現の抽象化とモジュール化を促進します。

DivMod2関数

  • func DivMod2(x, y *[]Digit2) (*[]Digit2, *[]Digit2): この関数は、[]Digit2スライスで表される2つの多倍長整数xyの除算と剰余を計算します。これは、KnuthのアルゴリズムDを実装したもので、多倍長除算の核心部分です。
    • 単一桁除算の最適化: m == 1(除数が単一桁)の場合、新しく導入されたDiv1関数を呼び出すことで、より効率的に処理します。
    • 正規化: 除数yの最上位桁が特定の範囲になるように、xyMul1関数でスケーリング(乗算)します。これにより、試行商の推定精度が向上します。
    • 試行商の推定: y1, y2(除数の上位2桁)とx0, x1, x2(被除数の上位3桁)を用いて、KnuthのアルゴリズムDに基づく試行商qを推定します。この推定は、ループ内で最も計算量の多い部分であり、その精度が全体のパフォーマンスに大きく影響します。
    • 乗算と減算: 推定された試行商qと除数yを乗算し、被除数xから減算します。
    • 修正: 減算の結果が負になった場合、試行商qをデクリメントし、除数yを加算して修正します。
    • 非正規化: 最後に、正規化のために行ったスケーリングを元に戻すために、剰余をDiv1関数で除算します。

Natural型の高レベルメソッド

Natural型のAdd, Sub, Mul, Shl, Shrなどのメソッドは、以前は直接的なループ処理を含んでいましたが、このコミットで新しく追加されたRaw Operations(Add, Sub, Mul, Shl, Shrなど)を呼び出すように変更されました。これにより、コードの重複が減り、可読性が向上し、低レベルの最適化がより容易になりました。

例えば、Natural.Mulは、以前は複雑なネストされたループを含んでいましたが、新しいMul Raw Operationを呼び出すように簡略化されました。

// 変更前 (簡略化)
func (x *Natural) Mul(y *Natural) *Natural {
	// ...
	for j := 0; j < m; j++ {
		d := y[j];
		if d != 0 {
			c := Digit(0);
			for i := 0; i < n; i++ {
				z1, z0 := Mul1(x[i], d); // 古いMul1
				c, z[i+j] = Split(c + z[i+j] + z0);
				c += z1;
			}
			z[n+j] = c;
		}
	}
	// ...
}

// 変更後
func (x *Natural) Mul(y *Natural) *Natural {
	n := len(x);
	m := len(y);
	z := new(Natural, n + m);
	Mul(z, x, y); // 新しいRaw OperationのMulを呼び出す
	return Normalize(z);
}

このように、このコミットは、多倍長整数演算のパフォーマンスを向上させるために、内部的な基数を変更し、KnuthのアルゴリズムDをより効率的に実装し、コードベース全体をリファクタリングしてモジュール化するという、包括的な技術的改善を行っています。

関連リンク

参考にした情報源リンク