- steps towards implementation of div and mod

- string conversion in different bases
- tracked down a 6g bug, use work-around for now

R=r
OCL=17981
CL=17981

usr/gri/bignum/bignum.go

@@ -14,26 +14,52 @@ package Bignum
 
 
 // ----------------------------------------------------------------------------
-// Support
+// Representation
+
+type Word uint64
+const LogW = 32;
 
-type Word uint32
+const LogH = 4;  // bits for a hex digit (= "small" number)
+const H = 1 << LogH;
 
-const N = 4;
-const L = 28;  // = sizeof(Word) * 8
+const L = LogW - LogH;  // must be even (for Mul1)
 const B = 1 << L;
 const M = B - 1;
 
 
+// For division
+
+const (
+	L3 = L / 3;
+	B3 = 1 << L3;
+	M3 = B3 - 1;
+)
+
+
+type (
+	Word3 uint32;
+	Natural3 [] Word3;
+)
+
+
+// ----------------------------------------------------------------------------
+// Support
+
 // TODO replace this with a Go built-in assert
-func ASSERT(p bool) {
+func assert(p bool) {
 	if !p {
-		panic("ASSERT failed");
+		panic("assert failed");
 	}
 }
 
 
+func init() {
+	assert(L % 2 == 0);  // L must be even
+}
+
+
 func IsSmall(x Word) bool {
-	return x < 1 << N;
+	return x < H;
 }
 
 
@@ -42,96 +68,136 @@ func Update(x Word) (Word, Word) {
 }
 
 
+export func Dump(x *[]Word) {
+	print("[", len(x), "]");
+	for i := len(x) - 1; i >= 0; i-- {
+		print(" ", x[i]);
+	}
+	println();
+}
+
+
 // ----------------------------------------------------------------------------
-// Naturals
+// Natural numbers
 
 export type Natural []Word;
 export var NatZero *Natural = new(Natural, 0);
 
 
+export func NewNat(x Word) *Natural {
+	var z *Natural;
+	switch {
+	case x == 0:
+		z = NatZero;
+	case x < B:
+		z = new(Natural, 1);
+		z[0] = x;
+		return z;
+	default:
+		z = new(Natural, 2);
+		z[0], z[1] = Update(x);
+	}
+	return z;
+}
+
+
+func Normalize(x *Natural) *Natural {
+	i := len(x);
+	for i > 0 && x[i - 1] == 0 { i-- }
+	if i < len(x) {
+		x = x[0 : i];  // trim leading 0's
+	}
+	return x;
+}
+
+
+func Normalize3(x *Natural3) *Natural3 {
+	i := len(x);
+	for i > 0 && x[i - 1] == 0 { i-- }
+	if i < len(x) {
+		x = x[0 : i];  // trim leading 0's
+	}
+	return x;
+}
+
+
 func (x *Natural) IsZero() bool {
 	return len(x) == 0;
 }
 
 
-func (x *Natural) Add (y *Natural) *Natural {
-	xl := len(x);
-	yl := len(y);
-	if xl < yl {
+func (x *Natural) Add(y *Natural) *Natural {
+	n := len(x);
+	m := len(y);
+	if n < m {
 		return y.Add(x);
 	}
-	ASSERT(xl >= yl);
-	z := new(Natural, xl + 1);
+	assert(n >= m);
+	z := new(Natural, n + 1);
 
 	i := 0;
 	c := Word(0);
-	for i < yl { z[i], c = Update(x[i] + y[i] + c); i++; }
-	for i < xl { z[i], c = Update(x[i] + c); i++; }
-	if c != 0 { z[i] = c; i++; }
-	z = z[0 : i];
+	for i < m { z[i], c = Update(x[i] + y[i] + c); i++; }
+	for i < n { z[i], c = Update(x[i] + c); i++; }
+	z[i] = c;
 
-	return z;
+	return Normalize(z);
 }
 
 
-func (x *Natural) Sub (y *Natural) *Natural {
-	xl := len(x);
-	yl := len(y);
-	ASSERT(xl >= yl);
-	z := new(Natural, xl);
+func (x *Natural) Sub(y *Natural) *Natural {
+	n := len(x);
+	m := len(y);
+	assert(n >= m);
+	z := new(Natural, n);
 
 	i := 0;
 	c := Word(0);
-	for i < yl { z[i], c = Update(x[i] - y[i] + c); i++; }
-	for i < xl { z[i], c = Update(x[i] + c); i++; }
-	ASSERT(c == 0);  // usub(x, y) must be called with x >= y
-	for i > 0 && z[i - 1] == 0 { i--; }
-	z = z[0 : i];
+	for i < m { z[i], c = Update(x[i] - y[i] + c); i++; }
+	for i < n { z[i], c = Update(x[i] + c); i++; }
+	assert(c == 0);  // x.Sub(y) must be called with x >= y
 
-	return z;
+	return Normalize(z);
 }
 
 
 // Computes x = x*a + c (in place) for "small" a's.
-func (x* Natural) Mul1Add(a, c Word) *Natural {
-	ASSERT(IsSmall(a) && IsSmall(c));
-	if (x.IsZero() || a == 0) && c == 0 {
-		return NatZero;
+func (x* Natural) MulAdd1(a, c Word) *Natural {
+	assert(IsSmall(a-1) && IsSmall(c));
+	if x.IsZero() || a == 0 {
+		return NewNat(c);
 	}
-	xl := len(x);
+	n := len(x);
 
-	z := new(Natural, xl + 1);
-	i := 0;
-	for i < xl { z[i], c = Update(x[i] * a + c); i++; }
-	if c != 0 { z[i] = c; i++; }  
-	z = z[0 : i];
+	z := new(Natural, n + 1);
+	for i := 0; i < n; i++ { z[i], c = Update(x[i] * a + c); }
+	z[n] = c;
 
-	return z;
+	return Normalize(z);
 }
 
 
 // Returns z = (x * y) div B, c = (x * y) mod B.
-func Mul1(x, y Word) (z Word, c Word) {
-	const L2 = (L + 1) >> 1;
+func Mul1(x, y Word) (Word, Word) {
+	const L2 = (L + 1) / 2;  // TODO check if we can run with odd L
 	const B2 = 1 << L2;
 	const M2 = B2 - 1;
-
+	
 	x0 := x & M2;
 	x1 := x >> L2;
 
 	y0 := y & M2;
 	y1 := y >> L2;
 
-	z10 := x0*y0;
-	z21 := x1*y0 + x0*y1 + z10 >> L2;
-
-	cc := x1*y1 + z21 >> L2;  
-	zz := z21 & M2 << L2 | z10 & M2;
-	return zz, cc
+	z0 := x0*y0;
+	z1 := x1*y0 + x0*y1 + z0 >> L2;  z0 &= M2;
+	z2 := x1*y1 + z1 >> L2;  z1 &= M2;
+	
+	return z1 << L2 | z0, z2;
 }
 
 
-func (x *Natural) Mul (y *Natural) *Natural {
+func (x *Natural) Mul(y *Natural) *Natural {
 	if x.IsZero() || y.IsZero() {
 		return NatZero;
 	}
@@ -140,17 +206,16 @@ func (x *Natural) Mul (y *Natural) *Natural {
 	if xl < yl {
 		return y.Mul(x);  // for speed
 	}
-	ASSERT(xl >= yl && yl > 0);
+	assert(xl >= yl && yl > 0);
 
 	// initialize z
 	zl := xl + yl;
 	z := new(Natural, zl);
 
-	k := 0;
 	for j := 0; j < yl; j++ {
 		d := y[j];
 		if d != 0 {
-			k = j;
+			k := j;
 			c := Word(0);
 			for i := 0; i < xl; i++ {
 				// compute z[k] += x[i] * d + c;
@@ -162,31 +227,118 @@ func (x *Natural) Mul (y *Natural) *Natural {
 				c += t >> L;
 				k++;
 			}
-			if c != 0 {
-				z[k] = c;
-				k++;
-			}
+			z[k] = c;
 		}
 	}
-	z = z[0 : k];
 
-	return z;
+	return Normalize(z);
+}
+
+
+func Shl1(x Word, s int) (Word, Word) {
+	return 0, 0
+}
+
+
+func Shr1(x Word, s int) (Word, Word) {
+	return 0, 0
 }
 
 
-func (x *Natural) Div (y *Natural) *Natural {
+func (x *Natural) Shl(s int) *Natural {
+	panic("incomplete");
+	
+	if s == 0 {
+		return x;
+	}
+	
+	S := s/L;
+	s = s%L;
+	n := len(x) + S + 1;
+	z := new(Natural, n);
+	
+	c := Word(0);
+	for i := 0; i < n; i++ {
+		z[i + S], c = Shl1(x[i], s);
+	}
+	z[n + S] = c;
+	
+	return Normalize(z);
+}
+
+
+func (x *Natural) Shr(s uint) *Natural {
+	panic("incomplete");
+	
+	if s == 0 {
+		return x;
+	}
+	return nil
+}
+
+
+func SplitBase(x *Natural) *Natural3 {
+	xl := len(x);
+	z := new(Natural3, xl * 3);
+	for i, j := 0, 0; i < xl; i, j = i + 1, j + 3 {
+		t := x[i];
+		z[j] = Word3(t & M3);  t >>= L3;  j++;
+		z[j] = Word3(t & M3);  t >>= L3;  j++;
+		z[j] = Word3(t & M3);  t >>= L3;  j++;
+	}
+	return Normalize3(z);
+}
+
+
+func Scale(x *Natural, f Word) *Natural3 {
+	return nil;
+}
+
+
+func TrialDigit(r, d *Natural3, k, m int) Word {
+	km := k + m;
+	assert(2 <= m && m <= km);
+	r3 := (Word(r[km]) << L3 + Word(r[km - 1])) << L3 + Word(r[km - 2]);
+	d2 := Word(d[m - 1]) << L3 + Word(d[m - 2]);
+	qt := r3 / d2;
+	if qt >= B {
+		qt = B - 1;
+	}
+	return qt;
+}
+
+
+func DivMod(x, y *Natural) {
+	xl := len(x);
+	yl := len(y);
+	assert(2 <= yl && yl <= xl);  // use special-case algorithm otherwise
+	
+	f := B / (y[yl - 1] + 1);
+	r := Scale(x, f);
+	d := Scale(y, f);
+	n := len(r);
+	m := len(d);
+	
+	for k := n - m; k >= 0; k-- {
+		qt := TrialDigit(r, d, k, m);
+		
+	}
+}
+
+
+func (x *Natural) Div(y *Natural) *Natural {
 	panic("UNIMPLEMENTED");
 	return nil;
 }
 
 
-func (x *Natural) Mod (y *Natural) *Natural {
+func (x *Natural) Mod(y *Natural) *Natural {
 	panic("UNIMPLEMENTED");
 	return nil;
 }
 
 
-func (x *Natural) Cmp (y *Natural) int {
+func (x *Natural) Cmp(y *Natural) int {
 	xl := len(x);
 	yl := len(y);
 
@@ -202,176 +354,173 @@ func (x *Natural) Cmp (y *Natural) int {
 	case x[i] < y[i]: d = -1;
 	case x[i] > y[i]: d = 1;
 	}
+
 	return d;
 }
 
 
 func (x *Natural) Log() int {
-	xl := len(x);
-	if xl == 0 { return 0; }
-
-	n := (xl - 1) * L;
-	for t := x[xl - 1]; t != 0; t >>= 1 { n++ };
+	n := len(x);
+	if n == 0 { return 0; }
+	assert(n > 0);
+	
+	c := (n - 1) * L;
+	for t := x[n - 1]; t != 0; t >>= 1 { c++ };
 
-	return n;
+	return c;
 }
 
 
-func (x *Natural) And (y *Natural) *Natural {
-	xl := len(x);
-	yl := len(y);
-	if xl < yl {
+func (x *Natural) And(y *Natural) *Natural {
+	n := len(x);
+	m := len(y);
+	if n < m {
 		return y.And(x);
 	}
-	ASSERT(xl >= yl);
-	z := new(Natural, xl);
+	assert(n >= m);
+	z := new(Natural, n);
 
 	i := 0;
-	for i < yl { z[i] = x[i] & y[i]; i++; }
-	for i < xl { z[i] = x[i]; i++; }
-	for i > 0 && z[i - 1] == 0 { i--; }
-	z = z[0 : i];
+	for i < m { z[i] = x[i] & y[i]; i++; }
+	for i < n { z[i] = x[i]; i++; }
 
-	return z;
+	return Normalize(z);
 }
 
 
-func (x *Natural) Or (y *Natural) *Natural {
-	xl := len(x);
-	yl := len(y);
-	if xl < yl {
+func (x *Natural) Or(y *Natural) *Natural {
+	n := len(x);
+	m := len(y);
+	if n < m {
 		return y.Or(x);
 	}
-	ASSERT(xl >= yl);
-	z := new(Natural, xl);
+	assert(n >= m);
+	z := new(Natural, n);
 
 	i := 0;
-	for i < yl { z[i] = x[i] | y[i]; i++; }
-	for i < xl { z[i] = x[i]; i++; }
+	for i < m { z[i] = x[i] | y[i]; i++; }
+	for i < n { z[i] = x[i]; i++; }
 
-	return z;
+	return Normalize(z);
 }
 
 
-func (x *Natural) Xor (y *Natural) *Natural {
-	xl := len(x);
-	yl := len(y);
-	if xl < yl {
+func (x *Natural) Xor(y *Natural) *Natural {
+	n := len(x);
+	m := len(y);
+	if n < m {
 		return y.Xor(x);
 	}
-	ASSERT(xl >= yl);
-	z := new(Natural, xl);
+	assert(n >= m);
+	z := new(Natural, n);
 
 	i := 0;
-	for i < yl { z[i] = x[i] ^ y[i]; i++; }
-	for i < xl { z[i] = x[i]; i++; }
-	for i > 0 && z[i - 1] == 0 { i--; }
-	z = z[0 : i];
+	for i < m { z[i] = x[i] ^ y[i]; i++; }
+	for i < n { z[i] = x[i]; i++; }
 
-	return z;
+	return Normalize(z);
 }
 
 
-// Returns a copy of x with space for one extra digit (for Div/Mod use)
 func Copy(x *Natural) *Natural {
-	xl := len(x);
-
-	z := new(Natural, xl + 1);  // add space for one extra digit
-	for i := 0; i < xl; i++ { z[i] = x[i]; }
-	z = z[0 : xl];
-
+	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) for "small" d's. Returns updated x, x mod d.
-func (x *Natural) Mod1 (d Word) (*Natural, Word) {
-	ASSERT(IsSmall(d));
-	xl := len(x);
-	
+// Computes x = x div d (in place - the recv maybe modified) for "small" d's.
+// Returns updated x and x mod d.
+func (x *Natural) DivMod1(d Word) (*Natural, Word) {
+	assert(0 < d && IsSmall(d - 1));
+
 	c := Word(0);
-	for i := xl - 1; i >= 0; i-- {
-		c = c << L + x[i];
-		x[i], c = c / d, c %d;
-	}
-	if xl > 0 && x[xl - 1] == 0 {
-		x = x[0 : xl - 1];
+	for i := len(x) - 1; i >= 0; i-- {
+		var LL Word = L;  // BUG shift broken for const L
+		c = c << LL + x[i];
+		x[i] = c / d;
+		c %= d;
 	}
 
-	return x, c;
+	return Normalize(x), c;
 }
 
 
-func (x *Natural) String() string {
+func (x *Natural) String(base Word) string {
 	if x.IsZero() {
 		return "0";
 	}
 	
 	// allocate string
+	// TODO n is too small for bases < 10!!!
+	assert(base >= 10);  // for now
 	// approx. length: 1 char for 3 bits
 	n := x.Log()/3 + 1;  // +1 (round up)
 	s := new([]byte, n);
 
 	// convert
+	const hex = "0123456789abcdef";
 	i := n;
 	x = Copy(x);  // don't destroy recv
 	for !x.IsZero() {
 		i--;
 		var d Word;
-		x, d = x.Mod1(10);
-		s[i] = byte(d) + '0';
+		x, d = x.DivMod1(base);
+		s[i] = hex[d];
 	};
 
 	return string(s[i : n]);
 }
 
 
-export func NatFromWord(x Word) *Natural {
-	var z *Natural;
-	switch {
-	case x == 0:
-		z = NatZero;
-	case x < B:
-		z = new(Natural, 1);
-		z[0] = x;
-		return z;
-	default:
-		z = new(Natural, 2);
-		z[0], z[1] = Update(x);
-	}
-	return z;
-}
-
-
-// Support function for faster factorial computation.
 func MulRange(a, b Word) *Natural {
 	switch {
-	case a > b: return NatFromWord(1);
-	case a == b: return NatFromWord(a);
-	case a + 1 == b: return NatFromWord(a).Mul(NatFromWord(b));
+	case a > b: return NewNat(1);
+	case a == b: return NewNat(a);
+	case a + 1 == b: return NewNat(a).Mul(NewNat(b));
 	}
 	m := (a + b) >> 1;
-	ASSERT(a <= m && m < b);
+	assert(a <= m && m < b);
 	return MulRange(a, m).Mul(MulRange(m + 1, b));
 }
 
 
 export func Fact(n Word) *Natural {
-  return MulRange(2, n);
+	// Using MulRange() instead of the basic for-loop
+	// lead to faster factorial computation.
+	return MulRange(2, n);
 }
 
 
-export func NatFromString(s string) *Natural {
+func HexValue(ch byte) Word {
+	d := Word(H);
+	switch {
+	case '0' <= ch && ch <= '9': d = Word(ch - '0');
+	case 'a' <= ch && ch <= 'f': d = Word(ch - 'a') + 10;
+	case 'A' <= ch && ch <= 'F': d = Word(ch - 'A') + 10;
+	}
+	return d;
+}
+
+
+// TODO auto-detect base if base argument is 0
+export func NatFromString(s string, base Word) *Natural {
 	x := NatZero;
-	for i := 0; i < len(s) && '0' <= s[i] && s[i] <= '9'; i++ {
-		x = x.Mul1Add(10, Word(s[i] - '0'));
+	for i := 0; i < len(s); i++ {
+		d := HexValue(s[i]);
+		if d < base {
+			x = x.MulAdd1(base, d);
+		} else {
+			break;
+		}
 	}
 	return x;
 }
 
 
 // ----------------------------------------------------------------------------
-// Integers
+// Integer numbers
 
 export type Integer struct {
 	sign bool;
@@ -379,7 +528,7 @@ export type Integer struct {
 }
 
 
-func (x *Integer) Add (y *Integer) *Integer {
+func (x *Integer) Add(y *Integer) *Integer {
 	var z *Integer;
 	if x.sign == y.sign {
 		// x + y == x + y
@@ -401,7 +550,7 @@ func (x *Integer) Add (y *Integer) *Integer {
 }
 
 
-func (x *Integer) Sub (y *Integer) *Integer {
+func (x *Integer) Sub(y *Integer) *Integer {
 	var z *Integer;
 	if x.sign != y.sign {
 		// x - (-y) == x + y
@@ -423,7 +572,7 @@ func (x *Integer) Sub (y *Integer) *Integer {
 }
 
 
-func (x *Integer) Mul (y *Integer) *Integer {
+func (x *Integer) Mul(y *Integer) *Integer {
 	// x * y == x * y
 	// x * (-y) == -(x * y)
 	// (-x) * y == -(x * y)
@@ -432,50 +581,48 @@ func (x *Integer) Mul (y *Integer) *Integer {
 }
 
 
-func (x *Integer) Div (y *Integer) *Integer {
+func (x *Integer) Div(y *Integer) *Integer {
 	panic("UNIMPLEMENTED");
 	return nil;
 }
 
 
-func (x *Integer) Mod (y *Integer) *Integer {
+func (x *Integer) Mod(y *Integer) *Integer {
 	panic("UNIMPLEMENTED");
 	return nil;
 }
 
 
-func (x *Integer) Cmp (y *Integer) int {
+func (x *Integer) Cmp(y *Integer) int {
 	panic("UNIMPLEMENTED");
 	return 0;
 }
 
 
-func (x *Integer) String() string {
+func (x *Integer) String(base Word) string {
 	if x.mant.IsZero() {
 		return "0";
 	}
 	var s string;
 	if x.sign {
-		s = "-" + x.mant.String();
-	} else {
-		s = x.mant.String();
+		s = "-";
 	}
-	return s;
+	return s + x.mant.String(base);
 }
 
 	
-export func IntFromString(s string) *Integer {
+export func IntFromString(s string, base Word) *Integer {
 	// get sign, if any
 	sign := false;
 	if len(s) > 0 && (s[0] == '-' || s[0] == '+') {
 		sign = s[0] == '-';
 	}
-	return &Integer{sign, NatFromString(s[1 : len(s)])};
+	return &Integer{sign, NatFromString(s[1 : len(s)], base)};
 }
 
 
 // ----------------------------------------------------------------------------
-// Rationals
+// Rational numbers
 
 export type Rational struct {
 	a, b *Integer;  // a = numerator, b = denominator
@@ -488,33 +635,33 @@ func NewRat(a, b *Integer) *Rational {
 }
 
 
-func (x *Rational) Add (y *Rational) *Rational {
+func (x *Rational) Add(y *Rational) *Rational {
 	return NewRat((x.a.Mul(y.b)).Add(x.b.Mul(y.a)), x.b.Mul(y.b));
 }
 
 
-func (x *Rational) Sub (y *Rational) *Rational {
+func (x *Rational) Sub(y *Rational) *Rational {
 	return NewRat((x.a.Mul(y.b)).Sub(x.b.Mul(y.a)), x.b.Mul(y.b));
 }
 
 
-func (x *Rational) Mul (y *Rational) *Rational {
+func (x *Rational) Mul(y *Rational) *Rational {
 	return NewRat(x.a.Mul(y.a), x.b.Mul(y.b));
 }
 
 
-func (x *Rational) Div (y *Rational) *Rational {
+func (x *Rational) Div(y *Rational) *Rational {
 	return NewRat(x.a.Mul(y.b), x.b.Mul(y.a));
 }
 
 
-func (x *Rational) Mod (y *Rational) *Rational {
+func (x *Rational) Mod(y *Rational) *Rational {
 	panic("UNIMPLEMENTED");
 	return nil;
 }
 
 
-func (x *Rational) Cmp (y *Rational) int {
+func (x *Rational) Cmp(y *Rational) int {
 	panic("UNIMPLEMENTED");
 	return 0;
 }
@@ -527,7 +674,7 @@ export func RatFromString(s string) *Rational {
 
 
 // ----------------------------------------------------------------------------
-// Numbers
+// Scaled numbers
 
 export type Number struct {
 	mant *Rational;

usr/gri/bignum/bignum_test.go

@@ -14,9 +14,9 @@ const (
 
 
 var (
-	a = Bignum.NatFromString(sa);
-	b = Bignum.NatFromString(sb);
-	c = Bignum.NatFromString(sc);
+	a = Bignum.NatFromString(sa, 10);
+	b = Bignum.NatFromString(sb, 10);
+	c = Bignum.NatFromString(sc, 10);
 )
 
 
@@ -28,12 +28,12 @@ func TEST(msg string, b bool) {
 
 
 func TestConv() {
-	TEST("TC1", a.Cmp(Bignum.NatFromWord(991)) == 0);
+	TEST("TC1", a.Cmp(Bignum.NewNat(991)) == 0);
 	TEST("TC2", b.Cmp(Bignum.Fact(20)) == 0);
 	TEST("TC3", c.Cmp(Bignum.Fact(100)) == 0);
-	TEST("TC4", a.String() == sa);
-	TEST("TC5", b.String() == sb);
-	TEST("TC6", c.String() == sc);
+	TEST("TC4", a.String(10) == sa);
+	TEST("TC5", b.String(10) == sb);
+	TEST("TC6", c.String(10) == sc);
 }