// - Natural unsigned integer numbers
// - Integer signed integer numbers
// - Rational rational numbers
-// - Number scaled rational numbers (contain exponent)
// ----------------------------------------------------------------------------
// Representation
+//
+// A natural number of the form
+//
+// x = x[n-1]*B^(n-1) + x[n-2]*B^(n-2) + ... + x[1]*B + x[0]
+//
+// with 0 <= x[i] < B and 0 <= i < n is stored in an array of length n,
+// with the digits x[i] as the array elements. 0 is represented as an
+// empty array (length == 0).
+//
+// A natural number is normalized if the array contains no leading 0 digits.
+// During arithmetic operations, denormalized values may occur which are
+// always normalized before returning the final result.
+//
+// The base B is chosen as large as possible on a given platform but there
+// are a few constraints besides the largest unsigned integer type available.
+// TODO describe the constraints.
-type Word uint64
-const LogW = 32;
+type Word uint64;
+const LogW = 64;
const LogH = 4; // bits for a hex digit (= "small" number)
const H = 1 << LogH;
-const L = LogW - LogH; // must be even (for Mul1)
-const B = 1 << L;
+const LogB = LogW - LogH;
+const L = LogB;
+const B = 1 << LogB;
const M = B - 1;
}
-func init() {
- assert(L % 2 == 0); // L must be even
-}
-
-
func IsSmall(x Word) bool {
return x < H;
}
-func Update(x Word) (Word, Word) {
- return x & M, x >> L;
+func Split(x Word) (Word, Word) {
+ return x>>L, x&M;
}
return z;
default:
z = new(Natural, 2);
- z[0], z[1] = Update(x);
+ z[1], z[0] = Split(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
+ n := len(x);
+ for n > 0 && x[n - 1] == 0 { n-- }
+ if n < len(x) {
+ x = x[0 : n]; // 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
+ n := len(x);
+ for n > 0 && x[n - 1] == 0 { n-- }
+ if n < len(x) {
+ x = x[0 : n]; // trim leading 0's
}
return x;
}
i := 0;
c := Word(0);
- for i < m { z[i], c = Update(x[i] + y[i] + c); i++; }
- for i < n { z[i], c = Update(x[i] + c); i++; }
+ for ; i < m; i++ { c, z[i] = Split(x[i] + y[i] + c); }
+ for ; i < n; i++ { c, z[i] = Split(x[i] + c); }
z[i] = c;
return Normalize(z);
i := 0;
c := Word(0);
- for i < m { z[i], c = Update(x[i] - y[i] + c); i++; }
- for i < n { z[i], c = Update(x[i] + c); i++; }
+ for ; i < m; i++ { c, z[i] = Split(x[i] - y[i] + c); }
+ for ; i < n; i++ { c, z[i] = Split(x[i] + c); }
assert(c == 0); // x.Sub(y) must be called with x >= y
return Normalize(z);
n := len(x);
z := new(Natural, n + 1);
- for i := 0; i < n; i++ { z[i], c = Update(x[i] * a + c); }
+ for i := 0; i < n; i++ { c, z[i] = Split(x[i]*a + c); }
z[n] = c;
return Normalize(z);
}
-// Returns z = (x * y) div B, c = (x * y) mod B.
+// Returns c = x*y div B, z = x*y mod B.
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;
+ // 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;
- x0 := x & M2;
- x1 := x >> L2;
-
- y0 := y & M2;
- y1 := y >> L2;
-
- z0 := x0*y0;
- z1 := x1*y0 + x0*y1 + z0 >> L2; z0 &= M2;
- z2 := x1*y1 + z1 >> L2; z1 &= M2;
-
- return z1 << L2 | z0, z2;
+ return z1, z0;
}
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);
+ n := len(x);
+ m := len(y);
+ z := new(Natural, n + m);
- for j := 0; j < yl; j++ {
+ for j := 0; j < m; 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++;
+ for i := 0; i < n; i++ {
+ // z[i+j] += x[i]*d + c;
+ z1, z0 := Mul1(x[i], d);
+ c, z[i+j] = Split(z[i+j] + z0 + c);
+ c += z1;
}
- z[k] = c;
+ z[n+j] = c;
}
}
}
-func Shl1(x Word, s int) (Word, Word) {
- return 0, 0
+// BUG use these until 6g shifts are working properly
+func shl(x Word, s uint) Word {
+ return x << s;
}
-func Shr1(x Word, s int) (Word, Word) {
- return 0, 0
+func shr(x Word, s uint) Word {
+ return x >> s;
}
-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);
+func Shl1(x, c Word, s uint) (Word, Word) {
+ assert(s <= LogB);
+ return shr(x, (LogB - s)), shl(x, s)&M | c
+}
+
+
+func (x *Natural) Shl(s uint) *Natural {
+ n := len(x);
+ si := int(s/LogB);
+ s = s%LogB;
+ z := new(Natural, n + si + 1);
+ i := 0;
c := Word(0);
- for i := 0; i < n; i++ {
- z[i + S], c = Shl1(x[i], s);
- }
- z[n + S] = c;
+ for ; i < n; i++ { c, z[i+si] = Shl1(x[i], c, s); }
+ z[i+si] = c;
return Normalize(z);
}
func (x *Natural) Shr(s uint) *Natural {
panic("incomplete");
-
- if s == 0 {
- return x;
- }
return nil
}
}
+func Log1(x Word) int {
+ n := -1;
+ for x != 0 { x >>= 1; n++; }
+ return n;
+}
+
+
func (x *Natural) Log() int {
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 c;
+ if n > 0 {
+ n = (n - 1)*L + Log1(x[n - 1]);
+ } else {
+ n = -1;
+ }
+ return n;
}
z := new(Natural, n);
i := 0;
- for i < m { z[i] = x[i] & y[i]; i++; }
- for i < n { z[i] = x[i]; i++; }
+ for ; i < m; i++ { z[i] = x[i] & y[i]; }
+ for ; i < n; i++ { z[i] = x[i]; }
return Normalize(z);
}
z := new(Natural, n);
i := 0;
- for i < m { z[i] = x[i] | y[i]; i++; }
- for i < n { z[i] = x[i]; i++; }
+ for ; i < m; i++ { z[i] = x[i] | y[i]; }
+ for ; i < n; i++ { z[i] = x[i]; }
return Normalize(z);
}
z := new(Natural, n);
i := 0;
- for i < m { z[i] = x[i] ^ y[i]; i++; }
- for i < n { z[i] = x[i]; i++; }
+ for ; i < m; i++ { z[i] = x[i] ^ y[i]; }
+ for ; i < n; i++ { z[i] = x[i]; }
return Normalize(z);
}
c := Word(0);
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 = c<<L + x[i];
+ x[i] = c/d;
c %= d;
}
// 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)
+ n := x.Log()/3 + 10; // +10 (round up) - what is the right number?
s := new([]byte, n);
// convert
case a == b: return NewNat(a);
case a + 1 == b: return NewNat(a).Mul(NewNat(b));
}
- m := (a + b) >> 1;
+ m := (a + b)>>1;
assert(a <= m && m < b);
return MulRange(a, m).Mul(MulRange(m + 1, b));
}
panic("UNIMPLEMENTED");
return nil;
}
-
-
-// ----------------------------------------------------------------------------
-// Scaled numbers
-
-export type Number struct {
- mant *Rational;
- exp Integer;
-}
)
-func TEST(msg string, b bool) {
+var test_msg string;
+func TEST(n int, b bool) {
if !b {
- panic("TEST failed: ", msg, "\n");
+ panic("TEST failed: ", test_msg, "(", n, ")\n");
}
}
func TestConv() {
- 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(10) == sa);
- TEST("TC5", b.String(10) == sb);
- TEST("TC6", c.String(10) == sc);
+ test_msg = "TestConv";
+ TEST(0, a.Cmp(Bignum.NewNat(991)) == 0);
+ TEST(1, b.Cmp(Bignum.Fact(20)) == 0);
+ TEST(2, c.Cmp(Bignum.Fact(100)) == 0);
+ TEST(3, a.String(10) == sa);
+ TEST(4, b.String(10) == sb);
+ TEST(5, c.String(10) == sc);
+}
+
+
+func TestShift() {
+ test_msg = "TestShiftA";
+ TEST(0, b.Shl(0).Cmp(b) == 0);
+ TEST(1, c.Shl(1).Cmp(c) > 0);
+
+ test_msg = "TestShiftB";
+ { const m = 3;
+ p := b;
+ f := Bignum.NewNat(1<<m);
+ for i := 0; i < 100; i++ {
+ TEST(i, b.Shl(uint(i*m)).Cmp(p) == 0);
+ p = p.Mul(f);
+ }
+ }
}
func main() {
TestConv();
+ TestShift();
print("PASSED\n");
}