--- /dev/null
+// Copyright 2009 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package Bignum
+
+// A package for arbitrary precision arithmethic.
+// It implements the following numeric types:
+//
+// - Natural unsigned integer numbers
+// - Integer signed integer numbers
+// - Rational rational numbers
+// - Number scaled rational numbers (contain exponent)
+
+
+// ----------------------------------------------------------------------------
+// Support
+
+type Word uint32
+
+const N = 4;
+const L = 28; // = sizeof(Word) * 8
+const M = 1 << L - 1;
+
+
+// TODO replace this with a Go built-in assert
+func ASSERT(p bool) {
+ if !p {
+ panic("ASSERT failed");
+ }
+}
+
+
+func Update(x Word) (Word, Word) {
+ return x & M, x >> L;
+}
+
+
+// ----------------------------------------------------------------------------
+// Naturals
+
+export type Natural []Word;
+export var NatZero *Natural = new(Natural, 0);
+
+
+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 {
+ return y.Add(x);
+ }
+ ASSERT(xl >= yl);
+ z := new(Natural, xl + 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];
+
+ return z;
+}
+
+
+func (x *Natural) Sub (y *Natural) *Natural {
+ xl := len(x);
+ yl := len(y);
+ ASSERT(xl >= yl);
+ z := new(Natural, xl);
+
+ 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];
+
+ return z;
+}
+
+
+// Computes x = x*a + c (in place) for "small" a's.
+func (x* Natural) Mul1Add(a, c Word) *Natural {
+ ASSERT(0 <= a && a < 1 << N);
+ ASSERT(0 <= c && c < 1 << N);
+ if (x.IsZero() || a == 0) && c == 0 {
+ return NatZero;
+ }
+ xl := 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];
+
+ return 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;
+ 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
+}
+
+
+func (x *Natural) Mul (y *Natural) *Natural {
+ if x.IsZero() || y.IsZero() {
+ return NatZero;
+ }
+ xl := len(x);
+ yl := len(y);
+ if xl < yl {
+ return y.Mul(x); // for speed
+ }
+ 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;
+ c := Word(0);
+ for i := 0; i < xl; i++ {
+ // compute z[k] += x[i] * d + c;
+ t := z[k] + c;
+ var z1 Word;
+ z1, c = Mul1(x[i], d);
+ t += z1;
+ z[k] = t & M;
+ c += t >> L;
+ k++;
+ }
+ if c != 0 {
+ z[k] = Word(c);
+ k++;
+ }
+ }
+ }
+ z = z[0 : k];
+
+ return z;
+}
+
+
+func (x *Natural) Div (y *Natural) *Natural {
+ panic("UNIMPLEMENTED");
+ return nil;
+}
+
+
+func (x *Natural) Mod (y *Natural) *Natural {
+ panic("UNIMPLEMENTED");
+ return nil;
+}
+
+
+func (x *Natural) Cmp (y *Natural) int {
+ xl := len(x);
+ yl := len(y);
+
+ if xl != yl || xl == 0 {
+ return xl - yl;
+ }
+
+ i := xl - 1;
+ for i > 0 && x[i] == y[i] { i--; }
+
+ d := 0;
+ switch {
+ 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++ };
+
+ return n;
+}
+
+
+func (x *Natural) And (y *Natural) *Natural {
+ xl := len(x);
+ yl := len(y);
+ if xl < yl {
+ return y.And(x);
+ }
+ ASSERT(xl >= yl);
+ z := new(Natural, xl);
+
+ 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];
+
+ return z;
+}
+
+
+func (x *Natural) Or (y *Natural) *Natural {
+ xl := len(x);
+ yl := len(y);
+ if xl < yl {
+ return y.And(x);
+ }
+ ASSERT(xl >= yl);
+ z := new(Natural, xl);
+
+ i := 0;
+ for i < yl { z[i] = x[i] | y[i]; i++; }
+ for i < xl { z[i] = x[i]; i++; }
+
+ return z;
+}
+
+
+func (x *Natural) Xor (y *Natural) *Natural {
+ xl := len(x);
+ yl := len(y);
+ if xl < yl {
+ return y.And(x);
+ }
+ ASSERT(xl >= yl);
+ z := new(Natural, xl);
+
+ 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];
+
+ return 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];
+
+ return z;
+}
+
+
+// Computes x = x div d (in place) for "small" d's. Returns x mod d.
+func (x *Natural) Mod1 (d Word) (*Natural, Word) {
+ ASSERT(0 < d && d < (1 << N));
+ xl := len(x);
+ c := Word(0);
+
+ i := xl;
+ for i > 0 {
+ i--;
+ c = c << L + x[i];
+
+ q := c / d;
+ x[i] = q;
+
+ //x[i] = c / d; // BUG
+
+ c = c % d;
+ }
+ if xl > 0 && x[xl - 1] == 0 {
+ x = x[0 : xl - 1];
+ if xl - 1 == 0 && len(x) != 0 {
+ panic();
+ }
+ }
+
+ return x, c;
+}
+
+
+func (x *Natural) String() string {
+ if x.IsZero() {
+ return "0";
+ }
+
+ // allocate string
+ // approx. length: 1 char for 3 bits
+ n := x.Log()/3 + 1; // +1 (round up)
+ s := new([]byte, n);
+
+ // convert
+ 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';
+ };
+
+ return string(s[i : n]);
+}
+
+
+export func NatFromWord(x Word) *Natural {
+ var z *Natural;
+ switch {
+ case x == 0:
+ z = NatZero;
+ case x < 2 << L:
+ 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));
+ }
+ m := (a + b) >> 1;
+ ASSERT(a <= m && m < b);
+ return MulRange(a, m).Mul(MulRange(m + 1, b));
+}
+
+
+export func Fact(n Word) *Natural {
+ return MulRange(2, n);
+}
+
+
+export func NatFromString(s string) *Natural {
+ x := NatZero;
+ for i := 0; i < len(s) && '0' <= s[i] && s[i] <= '9'; i++ {
+ x = x.Mul1Add(10, Word(s[i] - '0'));
+ }
+ return x;
+}
+
+
+// ----------------------------------------------------------------------------
+// Integers
+
+export type Integer struct {
+ sign bool;
+ mant *Natural;
+}
+
+
+func (x *Integer) Add (y *Integer) *Integer {
+ var z *Integer;
+ if x.sign == y.sign {
+ // x + y == x + y
+ // (-x) + (-y) == -(x + y)
+ z = &Integer{x.sign, x.mant.Add(y.mant)};
+ } else {
+ // x + (-y) == x - y == -(y - x)
+ // (-x) + y == y - x == -(x - y)
+ if x.mant.Cmp(y.mant) >= 0 {
+ z = &Integer{false, x.mant.Sub(y.mant)};
+ } else {
+ z = &Integer{true, y.mant.Sub(x.mant)};
+ }
+ }
+ if x.sign {
+ z.sign = !z.sign;
+ }
+ return z;
+}
+
+
+func (x *Integer) Sub (y *Integer) *Integer {
+ var z *Integer;
+ if x.sign != y.sign {
+ // x - (-y) == x + y
+ // (-x) - y == -(x + y)
+ z = &Integer{x.sign, x.mant.Add(y.mant)};
+ } else {
+ // x - y == x - y == -(y - x)
+ // (-x) - (-y) == y - x == -(x - y)
+ if x.mant.Cmp(y.mant) >= 0 {
+ z = &Integer{false, x.mant.Sub(y.mant)};
+ } else {
+ z = &Integer{true, y.mant.Sub(x.mant)};
+ }
+ }
+ if x.sign {
+ z.sign = !z.sign;
+ }
+ return z;
+}
+
+
+func (x *Integer) Mul (y *Integer) *Integer {
+ // x * y == x * y
+ // x * (-y) == -(x * y)
+ // (-x) * y == -(x * y)
+ // (-x) * (-y) == x * y
+ return &Integer{x.sign != y.sign, x.mant.Mul(y.mant)};
+}
+
+
+func (x *Integer) Div (y *Integer) *Integer {
+ panic("UNIMPLEMENTED");
+ return nil;
+}
+
+
+func (x *Integer) Mod (y *Integer) *Integer {
+ panic("UNIMPLEMENTED");
+ return nil;
+}
+
+
+func (x *Integer) Cmp (y *Integer) int {
+ panic("UNIMPLEMENTED");
+ return 0;
+}
+
+
+export func IntFromString(s string) *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)])};
+}
+
+
+// ----------------------------------------------------------------------------
+// Rationals
+
+export type Rational struct {
+ a, b *Integer; // a = numerator, b = denominator
+}
+
+
+func NewRat(a, b *Integer) *Rational {
+ // TODO normalize the rational
+ return &Rational{a, b};
+}
+
+
+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 {
+ 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 {
+ return NewRat(x.a.Mul(y.a), x.b.Mul(y.b));
+}
+
+
+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 {
+ panic("UNIMPLEMENTED");
+ return nil;
+}
+
+
+func (x *Rational) Cmp (y *Rational) int {
+ panic("UNIMPLEMENTED");
+ return 0;
+}
+
+
+export func RatFromString(s string) *Rational {
+ panic("UNIMPLEMENTED");
+ return nil;
+}
+
+
+// ----------------------------------------------------------------------------
+// Numbers
+
+export type Number struct {
+ mant *Rational;
+ exp Integer;
+}