From: Robert Griesemer Date: Wed, 29 Oct 2008 01:42:26 +0000 (-0700) Subject: - steps towards implementation of div and mod X-Git-Tag: weekly.2009-11-06~2876 X-Git-Url: http://www.git.cypherpunks.su/?a=commitdiff_plain;h=379b5a39213ab8c995a1726c9e461df31aa72aac;p=gostls13.git - 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 --- diff --git a/usr/gri/bignum/bignum.go b/usr/gri/bignum/bignum.go index 598035c14f..3e0c940f8e 100755 --- a/usr/gri/bignum/bignum.go +++ b/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; diff --git a/usr/gri/bignum/bignum_test.go b/usr/gri/bignum/bignum_test.go index fe3b1e4e85..e5e33a02d2 100644 --- a/usr/gri/bignum/bignum_test.go +++ b/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); }