From: Robert Griesemer Date: Tue, 4 Nov 2008 17:33:15 +0000 (-0800) Subject: - completed integer support (some logical functions missing) X-Git-Tag: weekly.2009-11-06~2807 X-Git-Url: http://www.git.cypherpunks.su/?a=commitdiff_plain;h=7cd11c1c09729cc3ade1289014865ed26b22c354;p=gostls13.git - completed integer support (some logical functions missing) - completed rational support - better documentation - more tests - cleanups R=r OCL=18438 CL=18438 --- diff --git a/usr/gri/bignum/bignum.go b/usr/gri/bignum/bignum.go index 0852f42262..47232e67b9 100755 --- a/usr/gri/bignum/bignum.go +++ b/usr/gri/bignum/bignum.go @@ -13,7 +13,7 @@ package Bignum // ---------------------------------------------------------------------------- -// Representation +// Internal representation // // A natural number of the form // @@ -27,37 +27,56 @@ package Bignum // always normalized before returning the final result. The normalized // representation of 0 is the empty array (length = 0). // +// The operations for all other numeric types are implemented on top of +// the operations for natural numbers. +// // The base B is chosen as large as possible on a given platform but there // are a few constraints besides the size of the largest unsigned integer -// type available. -// TODO describe the constraints. +// type available: +// +// 1) To improve conversion speed between strings and numbers, the base B +// is chosen such that division and multiplication by 10 (for decimal +// string representation) can be done without using extended-precision +// arithmetic. This makes addition, subtraction, and conversion routines +// twice as fast. It requires a "buffer" of 4 bits per operand digit. +// That is, the size of B must be 4 bits smaller then the size of the +// type (Digit) in which these operations are performed. Having this +// buffer also allows for trivial (single-bit) carry computation in +// addition and subtraction (optimization suggested by Ken Thompson). +// +// 2) Long division requires extended-precision (2-digit) division per digit. +// Instead of sacrificing the largest base type for all other operations, +// for division the operands are unpacked into "half-digits", and the +// results are packed again. For faster unpacking/packing, the base size +// in bits must be even. + +type ( + Digit uint64; + Digit2 uint32; // half-digits for division +) + const LogW = 64; const LogH = 4; // bits for a hex digit (= "small" number) -const LogB = LogW - LogH; +const LogB = LogW - LogH; // largest bit-width available const ( - L2 = LogB / 2; - B2 = 1 << L2; - M2 = B2 - 1; - - L = L2 * 2; - B = 1 << L; - M = B - 1; -) - - -type ( - Digit2 uint32; - Digit uint64; + // half-digits + W2 = LogB / 2; // width + B2 = 1 << W2; // base + M2 = B2 - 1; // mask + + // full digits + W = W2 * 2; // width + B = 1 << W; // base + M = B - 1; // mask ) // ---------------------------------------------------------------------------- -// Support +// Support functions -// TODO replace this with a Go built-in assert func assert(p bool) { if !p { panic("assert failed"); @@ -65,56 +84,36 @@ func assert(p bool) { } -// ---------------------------------------------------------------------------- -// Raw operations - -func And1(z, x *[]Digit, y Digit) { - for i := len(x) - 1; i >= 0; i-- { - z[i] = x[i] & y; - } -} - - -func And(z, x, y *[]Digit) { - for i := len(x) - 1; i >= 0; i-- { - z[i] = x[i] & y[i]; - } -} - - -func Or1(z, x *[]Digit, y Digit) { - for i := len(x) - 1; i >= 0; i-- { - z[i] = x[i] | y; - } -} - - -func Or(z, x, y *[]Digit) { - for i := len(x) - 1; i >= 0; i-- { - z[i] = x[i] | y[i]; - } +func IsSmall(x Digit) bool { + return x < 1<= 0; i-- { - z[i] = x[i] ^ y; + print(" ", x[i]); } + println(); } -func Xor(z, x, y *[]Digit) { - for i := len(x) - 1; i >= 0; i-- { - z[i] = x[i] ^ y[i]; - } -} +// ---------------------------------------------------------------------------- +// Raw operations on sequences of digits +// +// Naming conventions +// +// c carry +// x, y operands +// z result +// n, m len(x), len(y) 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 + c, z[i] = t>>W, t&M } return c; } @@ -125,7 +124,7 @@ func Add(z, x, y *[]Digit) Digit { n := len(x); for i := 0; i < n; i++ { t := c + x[i] + y[i]; - c, z[i] = t>>L, t&M + c, z[i] = t>>W, t&M } return c; } @@ -135,7 +134,7 @@ 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! + c, z[i] = Digit(int64(t)>>W), t&M; // requires arithmetic shift! } return c; } @@ -146,7 +145,7 @@ func Sub(z, x, y *[]Digit) 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! + c, z[i] = Digit(int64(t)>>W), t&M; // requires arithmetic shift! } return c; } @@ -154,31 +153,33 @@ func Sub(z, x, y *[]Digit) Digit { // 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)), + // Split x and y into 2 sub-digits each, // 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, x&m; - y1, y0 := y>>L0, y&m; + // x = (x1*B2 + x0) + // y = (y1*B2 + y0) + x1, x0 := x>>W2, x&M2; + y1, y0 := y>>W2, y&M2; - // x*y = t2*b^2 + t1*b + t0 + // x*y = t2*B2^2 + t1*B2 + 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)>>L1; + z0 := (t1<>W2)>>(W-W2); return z1, z0; } @@ -195,7 +196,7 @@ func Mul(z, x, y *[]Digit) { // 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, z[i+j] = t>>W, t&M; c += z1; } z[n+j] = c; @@ -204,118 +205,113 @@ func Mul(z, x, y *[]Digit) { } -func Mul1(z, x *[]Digit2, y Digit2) Digit2 { +func Shl(z, x *[]Digit, s uint) Digit { + assert(s <= W); 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); + c, z[i] = x[i] >> (W-s), x[i] << s & M | c; } - return Digit2(c); + return c; } -func Div1(z, x *[]Digit2, y Digit2) Digit2 { +func Shr(z, x *[]Digit, s uint) Digit { + assert(s <= W); 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); + for i := n - 1; i >= 0; i-- { + c, z[i] = x[i] << (W-s) & M, x[i] >> s | c; } - return Digit2(c); + return 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; +func And1(z, x *[]Digit, y Digit) { + for i := len(x) - 1; i >= 0; i-- { + z[i] = x[i] & y; } - 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; +func And(z, x, y *[]Digit) { + for i := len(x) - 1; i >= 0; i-- { + z[i] = x[i] & y[i]; } - return c; } -// ---------------------------------------------------------------------------- -// Support +func Or1(z, x *[]Digit, y Digit) { + for i := len(x) - 1; i >= 0; i-- { + z[i] = x[i] | y; + } +} -func IsSmall(x Digit) bool { - return x < 1<= 0; i-- { + z[i] = x[i] | y[i]; + } } -func Split(x Digit) (Digit, Digit) { - return x>>L, x&M; +func Xor1(z, x *[]Digit, y Digit) { + for i := len(x) - 1; i >= 0; i-- { + z[i] = x[i] ^ y; + } } -export func Dump(x *[]Digit) { - print("[", len(x), "]"); +func Xor(z, x, y *[]Digit) { for i := len(x) - 1; i >= 0; i-- { - print(" ", x[i]); + z[i] = x[i] ^ y[i]; } - println(); } // ---------------------------------------------------------------------------- // Natural numbers -// -// Naming conventions -// -// B, b bases -// c carry -// x, y operands -// z result -// n, m n = len(x), m = len(y) export type Natural []Digit; -export var NatZero *Natural = new(Natural, 0); + +var ( + NatZero *Natural = &Natural{}; + NatOne *Natural = &Natural{1}; + NatTwo *Natural = &Natural{2}; + NatTen *Natural = &Natural{10}; +) -export func Nat(x Digit) *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[1], z[0] = Split(x); +// Creation + +export func Nat(x uint) *Natural { + switch x { + case 0: return NatZero; + case 1: return NatOne; + case 2: return NatTwo; + case 10: return NatTen; } - return z; + assert(Digit(x) < B); + return &Natural{Digit(x)}; } -func Normalize(x *Natural) *Natural { - 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; +// Predicates + +func (x *Natural) IsOdd() bool { + return len(x) > 0 && x[0]&1 != 0; +} + + +func (x *Natural) IsZero() bool { + return len(x) == 0; } -func Normalize2(x *[]Digit2) *[]Digit2 { +// Operations + +func Normalize(x *Natural) *Natural { n := len(x); for n > 0 && x[n - 1] == 0 { n-- } if n < len(x) { @@ -325,12 +321,6 @@ func Normalize2(x *[]Digit2) *[]Digit2 { } -// 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 { n := len(x); m := len(y); @@ -363,19 +353,6 @@ func (x *Natural) Sub(y *Natural) *Natural { } -// Computes x = x*a + c (in place) for "small" a's. -func (x* Natural) MulAdd1(a, c Digit) *Natural { - assert(IsSmall(a-1) && IsSmall(c)); - n := len(x); - z := new(Natural, n + 1); - - for i := 0; i < n; i++ { c, z[i] = Split(c + x[i]*a); } - z[n] = c; - - return Normalize(z); -} - - func (x *Natural) Mul(y *Natural) *Natural { n := len(x); m := len(y); @@ -387,77 +364,24 @@ func (x *Natural) Mul(y *Natural) *Natural { } -func Pop1(x Digit) uint { - n := uint(0); - for x != 0 { - x &= x-1; - n++; - } - return n; -} - - -func (x *Natural) Pop() uint { - n := uint(0); - for i := len(x) - 1; i >= 0; i-- { - n += Pop1(x[i]); - } - return n; -} - - -func (x *Natural) Pow(n uint) *Natural { - z := Nat(1); - for n > 0 { - // z * x^n == x^n0 - if n&1 == 1 { - z = z.Mul(x); - } - x, n = x.Mul(x), n/2; - } - return z; -} - - -func (x *Natural) Shl(s uint) *Natural { - n := uint(len(x)); - m := n + s/L; - z := new(Natural, m+1); - - z[m] = Shl(z[m-n : m], x, s%L); - - return Normalize(z); -} - - -func (x *Natural) Shr(s uint) *Natural { - n := uint(len(x)); - m := n - s/L; - if m > n { // check for underflow - m = 0; - } - z := new(Natural, m); - - Shr(z, x[n-m : n], s%L); - - return Normalize(z); -} - - // DivMod needs multi-precision division which is not available if Digit -// is already using the largest uint size. Split base before division, -// and merge again after. Each Digit is split into 2 Digit2's. +// is already using the largest uint size. Instead, unpack each operand +// into operands with twice as many digits of half the size (Digit2), do +// DivMod, and then pack the results again. func Unpack(x *Natural) *[]Digit2 { - // TODO Use Log() for better result - don't need Normalize2 at the end! n := len(x); z := new([]Digit2, n*2 + 1); // add space for extra digit (used by DivMod) for i := 0; i < n; i++ { t := x[i]; z[i*2] = Digit2(t & M2); - z[i*2 + 1] = Digit2(t >> L2 & M2); + z[i*2 + 1] = Digit2(t >> W2 & M2); } - return Normalize2(z); + + // normalize result + k := 2*n; + for k > 0 && z[k - 1] == 0 { k-- } + return z[0 : k]; // trim leading 0's } @@ -470,34 +394,65 @@ func Pack(x *[]Digit2) *Natural { z[n] = Digit(x[n*2]); } for i := 0; i < n; i++ { - z[i] = Digit(x[i*2 + 1]) << L2 | Digit(x[i*2]); + z[i] = Digit(x[i*2 + 1]) << W2 | Digit(x[i*2]); } return Normalize(z); } -// Division and modulo computation - destroys x and y. Based on the -// algorithms described in: +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>>W2, 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); +} + + +// DivMod returns q and r with x = y*q + r and 0 <= r < y. +// x and y are destroyed in the process. +// +// The algorithm used here is based on 1). 2) describes the same algorithm +// in C. A discussion and summary of the relevant theorems can be found in +// 3). 3) also describes an easier way to obtain the trial digit - however +// it relies on tripple-precision arithmetic which is why Knuth's method is +// used here. // // 1) D. Knuth, "The Art of Computer Programming. Volume 2. Seminumerical // Algorithms." Addison-Wesley, Reading, 1969. +// (Algorithm D, Sec. 4.3.1) +// +// 2) Henry S. Warren, Jr., "A Hacker's Delight". Addison-Wesley, 2003. +// (9-2 Multiword Division, p.140ff) // -// 2) P. Brinch Hansen, Multiple-length division revisited: A tour of the +// 3) P. Brinch Hansen, Multiple-length division revisited: A tour of the // minefield. "Software - Practice and Experience 24", (June 1994), // 579-601. John Wiley & Sons, Ltd. -// -// Specifically, the inplace computation of quotient and remainder -// is described in 1), while 2) provides the background for a more -// accurate initial guess of the trial digit. -func DivMod2(x, y *[]Digit2) (*[]Digit2, *[]Digit2) { - const b = B2; - +func DivMod(x, y *[]Digit2) (*[]Digit2, *[]Digit2) { n := len(x); m := len(y); - assert(m > 0); // division by zero - assert(n+1 <= cap(x)); // space for one extra digit (should it be == ?) + if m == 0 { + panic("division by zero"); + } + assert(n+1 <= cap(x)); // space for one extra digit x = x[0 : n + 1]; + assert(x[n] == 0); if m == 1 { // division by single digit @@ -505,36 +460,39 @@ func DivMod2(x, y *[]Digit2) (*[]Digit2, *[]Digit2) { x[0] = Div1(x[1 : n+1], x[0 : n], y[0]); } else if m > n { - // quotient = 0, remainder = x - // TODO in this case we shouldn't even split base - FIX THIS + // y > x => quotient = 0, remainder = x + // TODO in this case we shouldn't even unpack x and y m = n; } else { // general case assert(2 <= m && m <= n); - assert(x[n] == 0); // normalize x and y - f := b/(Digit(y[m-1]) + 1); - Mul1(x, x, Digit2(f)); - Mul1(y, y, Digit2(f)); - assert(b/2 <= y[m-1] && y[m-1] < b); // incorrect scaling + // TODO Instead of multiplying, it would be sufficient to + // shift y such that the normalization condition is + // satisfied (as done in "Hacker's Delight"). + f := B2 / (Digit(y[m-1]) + 1); + if f != 1 { + Mul1(x, x, Digit2(f)); + Mul1(y, y, Digit2(f)); + } + assert(B2/2 <= y[m-1] && y[m-1] < B2); // incorrect scaling y1, y2 := Digit(y[m-1]), Digit(y[m-2]); - d2 := Digit(y1)*b + Digit(y2); + d2 := Digit(y1)<= 0; i-- { k := i+m; - // compute trial digit + // compute trial digit (Knuth) var q Digit; - { // Knuth - x0, x1, x2 := Digit(x[k]), Digit(x[k-1]), Digit(x[k-2]); + { x0, x1, x2 := Digit(x[k]), Digit(x[k-1]), Digit(x[k-2]); if x0 != y1 { - q = (x0*b + x1)/y1; + q = (x0< (x0*b + x1 - y1*q)*b + x2 { + for y2*q > (x0<>L2), Digit2(t&M2); + t := c + Digit(x[i+j]) - Digit(y[j])*q; + c, x[i+j] = Digit(int64(t)>>W2), Digit2(t&M2); // requires arithmetic shift! } // correct if trial digit was too large @@ -552,7 +510,7 @@ func DivMod2(x, y *[]Digit2) (*[]Digit2, *[]Digit2) { c := Digit(0); for j := 0; j < m; j++ { t := c + Digit(x[i+j]) + Digit(y[j]); - c, x[i+j] = uint64(int64(t) >> L2), Digit2(t & M2) + c, x[i+j] = t >> W2, Digit2(t & M2) } assert(c + Digit(x[k]) == 0); // correct trial digit @@ -563,8 +521,10 @@ func DivMod2(x, y *[]Digit2) (*[]Digit2, *[]Digit2) { } // undo normalization for remainder - c := Div1(x[0 : m], x[0 : m], Digit2(f)); - assert(c == 0); + if f != 1 { + c := Div1(x[0 : m], x[0 : m], Digit2(f)); + assert(c == 0); + } } return x[m : n+1], x[0 : m]; @@ -572,59 +532,45 @@ func DivMod2(x, y *[]Digit2) (*[]Digit2, *[]Digit2) { func (x *Natural) Div(y *Natural) *Natural { - q, r := DivMod2(Unpack(x), Unpack(y)); + q, r := DivMod(Unpack(x), Unpack(y)); return Pack(q); } func (x *Natural) Mod(y *Natural) *Natural { - q, r := DivMod2(Unpack(x), Unpack(y)); + q, r := DivMod(Unpack(x), Unpack(y)); return Pack(r); } func (x *Natural) DivMod(y *Natural) (*Natural, *Natural) { - q, r := DivMod2(Unpack(x), Unpack(y)); + q, r := DivMod(Unpack(x), Unpack(y)); return Pack(q), Pack(r); } -func (x *Natural) Cmp(y *Natural) int { - n := len(x); - m := len(y); +func (x *Natural) Shl(s uint) *Natural { + n := uint(len(x)); + m := n + s/W; + z := new(Natural, m+1); + + z[m] = Shl(z[m-n : m], x, s%W); + + return Normalize(z); +} - if n != m || n == 0 { - return n - m; - } - i := n - 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 Log2(x Digit) int { - n := -1; - for x != 0 { x = x >> 1; n++; } // BUG >>= broken for uint64 - return n; -} - - -func (x *Natural) Log2() int { - n := len(x); - if n > 0 { - n = (n - 1)*L + Log2(x[n - 1]); - } else { - n = -1; +func (x *Natural) Shr(s uint) *Natural { + n := uint(len(x)); + m := n - s/W; + if m > n { // check for underflow + m = 0; } - return n; + z := new(Natural, m); + + Shr(z, x[n-m : n], s%W); + + return Normalize(z); } @@ -673,14 +619,55 @@ func (x *Natural) Xor(y *Natural) *Natural { } -// Computes x = x div d (in place - the recv maybe modified) for "small" d's. +func (x *Natural) Cmp(y *Natural) int { + n := len(x); + m := len(y); + + if n != m || n == 0 { + return n - m; + } + + i := n - 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 Log2(x Digit) uint { + assert(x > 0); + n := uint(0); + for x > 0 { + x >>= 1; + n++; + } + return n - 1; +} + + +func (x *Natural) Log2() uint { + n := len(x); + if n > 0 { + return (uint(n) - 1)*W + Log2(x[n - 1]); + } + panic("Log2(0)"); +} + + +// Computes x = x div d in place (modifies x) for "small" d's. // Returns updated x and x mod d. -func (x *Natural) DivMod1(d Digit) (*Natural, Digit) { +func DivMod1(x *Natural, d Digit) (*Natural, Digit) { assert(0 < d && IsSmall(d - 1)); c := Digit(0); for i := len(x) - 1; i >= 0; i-- { - t := c<>W, t&M; + } + z[n] = c; + + return Normalize(z); +} + + +// Determines base (octal, decimal, hexadecimal) if base == 0. +export func NatFromString(s string, base uint, slen *int) *Natural { + // determine base if necessary + i, n := 0, len(s); + if base == 0 { + base = 10; + if n > 0 && s[0] == '0' { + if n > 1 && (s[1] == 'x' || s[1] == 'X') { + base, i = 16, 2; + } else { + base, i = 8, 1; + } + } + } + + // convert string + assert(2 <= base && base <= 16); + x := Nat(0); + for ; i < n; i++ { + d := HexValue(s[i]); + if d < base { + x = MulAdd1(x, Digit(base), Digit(d)); + } else { + break; + } + } + + // provide number of string bytes consumed if necessary + if slen != nil { + *slen = i; + } + + return x; +} + + +// Natural number functions + +func Pop1(x Digit) uint { + n := uint(0); + for x != 0 { + x &= x-1; + n++; + } + return n; +} + + +func (x *Natural) Pop() uint { + n := uint(0); + for i := len(x) - 1; i >= 0; i-- { + n += Pop1(x[i]); + } + return n; +} + + +func (x *Natural) Pow(n uint) *Natural { + z := Nat(1); + for n > 0 { + // z * x^n == x^n0 + if n&1 == 1 { + z = z.Mul(x); + } + x, n = x.Mul(x), n/2; + } + return z; +} + + +export func MulRange(a, b uint) *Natural { switch { case a > b: return Nat(1); case a == b: return Nat(a); @@ -728,7 +811,7 @@ export func MulRange(a, b Digit) *Natural { } -export func Fact(n Digit) *Natural { +export func Fact(n uint) *Natural { // Using MulRange() instead of the basic for-loop // lead to faster factorial computation. return MulRange(2, n); @@ -744,60 +827,73 @@ func (x *Natural) Gcd(y *Natural) *Natural { } -func HexValue(ch byte) uint { - d := uint(1 << LogH); - switch { - case '0' <= ch && ch <= '9': d = uint(ch - '0'); - case 'a' <= ch && ch <= 'f': d = uint(ch - 'a') + 10; - case 'A' <= ch && ch <= 'F': d = uint(ch - 'A') + 10; +// ---------------------------------------------------------------------------- +// Integer numbers +// +// Integers are normalized if the mantissa is normalized and the sign is +// false for mant == 0. Use MakeInt to create normalized Integers. + +export type Integer struct { + sign bool; + mant *Natural; +} + + +// Creation + +export func MakeInt(sign bool, mant *Natural) *Integer { + if mant.IsZero() { + sign = false; // normalize } - return d; + return &Integer{sign, mant}; } -// TODO auto-detect base if base argument is 0 -export func NatFromString(s string, base uint) *Natural { - x := NatZero; - for i := 0; i < len(s); i++ { - d := HexValue(s[i]); - if d < base { - x = x.MulAdd1(Digit(base), Digit(d)); +export func Int(x int) *Integer { + sign := false; + var ux uint; + if x < 0 { + sign = true; + if -x == x { + // smallest negative integer + t := ^0; + ux = ^(uint(t) >> 1); } else { - break; + ux = uint(-x); } + } else { + ux = uint(x); } - return x; + return MakeInt(sign, Nat(ux)); } -// ---------------------------------------------------------------------------- -// Algorithms +// Predicates -export type T interface { - IsZero() bool; - Mod(y T) bool; +func (x *Integer) IsOdd() bool { + return x.mant.IsOdd(); } -export func Gcd(x, y T) T { - // Euclidean algorithm. - for !y.IsZero() { - x, y = y, x.Mod(y); - } - return x; + +func (x *Integer) IsZero() bool { + return x.mant.IsZero(); } -// ---------------------------------------------------------------------------- -// Integer numbers +func (x *Integer) IsNeg() bool { + return x.sign && !x.mant.IsZero() +} -export type Integer struct { - sign bool; - mant *Natural; + +func (x *Integer) IsPos() bool { + return !x.sign && !x.mant.IsZero() } -export func Int(x int64) *Integer { - return nil; +// Operations + +func (x *Integer) Neg() *Integer { + return MakeInt(!x.sign, x.mant); } @@ -806,14 +902,14 @@ func (x *Integer) Add(y *Integer) *Integer { if x.sign == y.sign { // x + y == x + y // (-x) + (-y) == -(x + y) - z = &Integer{x.sign, x.mant.Add(y.mant)}; + z = MakeInt(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)}; + z = MakeInt(false, x.mant.Sub(y.mant)); } else { - z = &Integer{true, y.mant.Sub(x.mant)}; + z = MakeInt(true, y.mant.Sub(x.mant)); } } if x.sign { @@ -828,14 +924,14 @@ func (x *Integer) Sub(y *Integer) *Integer { if x.sign != y.sign { // x - (-y) == x + y // (-x) - y == -(x + y) - z = &Integer{x.sign, x.mant.Add(y.mant)}; + z = MakeInt(false, 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)}; + z = MakeInt(false, x.mant.Sub(y.mant)); } else { - z = &Integer{true, y.mant.Sub(x.mant)}; + z = MakeInt(true, y.mant.Sub(x.mant)); } } if x.sign { @@ -850,16 +946,30 @@ func (x *Integer) Mul(y *Integer) *Integer { // x * (-y) == -(x * y) // (-x) * y == -(x * y) // (-x) * (-y) == x * y - return &Integer{x.sign != y.sign, x.mant.Mul(y.mant)}; + return MakeInt(x.sign != y.sign, x.mant.Mul(y.mant)); } +func (x *Integer) MulNat(y *Natural) *Integer { + // x * y == x * y + // (-x) * y == -(x * y) + return MakeInt(x.sign, x.mant.Mul(y)); +} + + +// Quo and Rem implement T-division and modulus (like C99): +// +// q = x.Quo(y) = trunc(x/y) (truncation towards zero) +// r = x.Rem(y) = x - y*q +// +// ( Daan Leijen, "Division and Modulus for Computer Scientists". ) + func (x *Integer) Quo(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.Div(y.mant)}; + return MakeInt(x.sign != y.sign, x.mant.Div(y.mant)); } @@ -868,31 +978,116 @@ func (x *Integer) Rem(y *Integer) *Integer { // x % (-y) == x % y // (-x) % y == -(x % y) // (-x) % (-y) == -(x % y) - return &Integer{y.sign, x.mant.Mod(y.mant)}; + return MakeInt(x.sign, x.mant.Mod(y.mant)); } func (x *Integer) QuoRem(y *Integer) (*Integer, *Integer) { q, r := x.mant.DivMod(y.mant); - return &Integer{x.sign != y.sign, q}, &Integer{y.sign, q}; + return MakeInt(x.sign != y.sign, q), MakeInt(x.sign, r); } +// Div and Mod implement Euclidian division and modulus: +// +// d = x.Div(y) +// m = x.Mod(y) with: 0 <= m < |d| and: y = x*d + m +// +// ( Raymond T. Boute, The Euclidian definition of the functions +// div and mod. "ACM Transactions on Programming Languages and +// Systems (TOPLAS)", 14(2):127-144, New York, NY, USA, 4/1992. +// ACM press. ) + + func (x *Integer) Div(y *Integer) *Integer { - q, r := x.mant.DivMod(y.mant); - return nil; + q, r := x.QuoRem(y); + if r.IsNeg() { + if y.IsPos() { + q = q.Sub(Int(1)); + } else { + q = q.Add(Int(1)); + } + } + return q; } func (x *Integer) Mod(y *Integer) *Integer { + r := x.Rem(y); + if r.IsNeg() { + if y.IsPos() { + r = r.Add(y); + } else { + r = r.Sub(y); + } + } + return r; +} + + +func (x *Integer) DivMod(y *Integer) (*Integer, *Integer) { + q, r := x.QuoRem(y); + if r.IsNeg() { + if y.IsPos() { + q = q.Sub(Int(1)); + r = r.Add(y); + } else { + q = q.Add(Int(1)); + r = r.Sub(y); + } + } + return q, r; +} + + +func (x *Integer) Shl(s uint) *Integer { + return MakeInt(x.sign, x.mant.Shl(s)); +} + + +func (x *Integer) Shr(s uint) *Integer { + z := MakeInt(x.sign, x.mant.Shr(s)); + if x.IsNeg() { + panic("UNIMPLEMENTED"); + } + return z; +} + + +func (x *Integer) And(y *Integer) *Integer { panic("UNIMPLEMENTED"); return nil; } -func (x *Integer) Cmp(y *Integer) int { +func (x *Integer) Or(y *Integer) *Integer { + panic("UNIMPLEMENTED"); + return nil; +} + + +func (x *Integer) Xor(y *Integer) *Integer { panic("UNIMPLEMENTED"); - return 0; + return nil; +} + + +func (x *Integer) Cmp(y *Integer) int { + // x cmp y == x cmp y + // x cmp (-y) == x + // (-x) cmp y == y + // (-x) cmp (-y) == -(x cmp y) + var r int; + switch { + case x.sign == y.sign: + r = x.mant.Cmp(y.mant); + if x.sign { + r = -r; + } + case x.sign: r = -1; + case y.sign: r = 1; + } + return r; } @@ -908,13 +1103,23 @@ func (x *Integer) String(base uint) string { } -export func IntFromString(s string, base uint) *Integer { +// Determines base (octal, decimal, hexadecimal) if base == 0. +export func IntFromString(s string, base uint, slen *int) *Integer { // get sign, if any sign := false; if len(s) > 0 && (s[0] == '-' || s[0] == '+') { sign = s[0] == '-'; + s = s[1 : len(s)]; + } + + z := MakeInt(sign, NatFromString(s, base, slen)); + + // correct slen if necessary + if slen != nil && sign { + *slen++; } - return &Integer{sign, NatFromString(s[1 : len(s)], base)}; + + return z; } @@ -922,56 +1127,119 @@ export func IntFromString(s string, base uint) *Integer { // Rational numbers export type Rational struct { - a, b *Integer; // a = numerator, b = denominator + a *Integer; // numerator + b *Natural; // denominator } -func (x *Rational) Normalize() *Rational { - f := x.a.mant.Gcd(x.b.mant); - x.a.mant = x.a.mant.Div(f); - x.b.mant = x.b.mant.Div(f); - return x; +// Creation + +export func MakeRat(a *Integer, b *Natural) *Rational { + f := a.mant.Gcd(b); // f > 0 + if f.Cmp(Nat(1)) != 0 { + a = MakeInt(a.sign, a.mant.Div(f)); + b = b.Div(f); + } + return &Rational{a, b}; +} + + +export func Rat(a0 int, b0 int) *Rational { + a, b := Int(a0), Int(b0); + if b.sign { + a = a.Neg(); + } + return MakeRat(a, b.mant); } -func Rat(a, b *Integer) *Rational { - return (&Rational{a, b}).Normalize(); +// Predicates + +func (x *Rational) IsZero() bool { + return x.a.IsZero(); +} + + +func (x *Rational) IsNeg() bool { + return x.a.IsNeg(); +} + + +func (x *Rational) IsPos() bool { + return x.a.IsPos(); +} + + +func (x *Rational) IsInt() bool { + return x.b.Cmp(Nat(1)) == 0; +} + + +// Operations + +func (x *Rational) Neg() *Rational { + return MakeRat(x.a.Neg(), x.b); } func (x *Rational) Add(y *Rational) *Rational { - return Rat((x.a.Mul(y.b)).Add(x.b.Mul(y.a)), x.b.Mul(y.b)); + return MakeRat((x.a.MulNat(y.b)).Add(y.a.MulNat(x.b)), x.b.Mul(y.b)); } func (x *Rational) Sub(y *Rational) *Rational { - return Rat((x.a.Mul(y.b)).Sub(x.b.Mul(y.a)), x.b.Mul(y.b)); + return MakeRat((x.a.MulNat(y.b)).Sub(y.a.MulNat(x.b)), x.b.Mul(y.b)); } func (x *Rational) Mul(y *Rational) *Rational { - return Rat(x.a.Mul(y.a), x.b.Mul(y.b)); + return MakeRat(x.a.Mul(y.a), x.b.Mul(y.b)); } -func (x *Rational) Div(y *Rational) *Rational { - return Rat(x.a.Mul(y.b), x.b.Mul(y.a)); +func (x *Rational) Quo(y *Rational) *Rational { + a := x.a.MulNat(y.b); + b := y.a.MulNat(x.b); + if b.IsNeg() { + a = a.Neg(); + } + return MakeRat(a, b.mant); } -func (x *Rational) Mod(y *Rational) *Rational { - panic("UNIMPLEMENTED"); - return nil; +func (x *Rational) Cmp(y *Rational) int { + return (x.a.MulNat(y.b)).Cmp(y.a.MulNat(x.b)); } -func (x *Rational) Cmp(y *Rational) int { - panic("UNIMPLEMENTED"); - return 0; +func (x *Rational) String(base uint) string { + s := x.a.String(base); + if !x.IsInt() { + s += "/" + x.b.String(base); + } + return s; } -export func RatFromString(s string) *Rational { - panic("UNIMPLEMENTED"); - return nil; +// Determines base (octal, decimal, hexadecimal) if base == 0. +export func RatFromString(s string, base uint, slen *int) *Rational { + // read nominator + var alen, blen int; + a := IntFromString(s, base, &alen); + b := Nat(1); + + // read denominator, if any + if alen < len(s) && s[alen] == '/' { + alen++; + if alen < len(s) { + b = NatFromString(s[alen : len(s)], base, &blen); + } + } + + // provide number of string bytes consumed if necessary + if slen != nil { + *slen = alen + blen; + } + + return MakeRat(a, b); } diff --git a/usr/gri/bignum/bignum_test.go b/usr/gri/bignum/bignum_test.go index 10d8da7db7..eb7d3beeec 100644 --- a/usr/gri/bignum/bignum_test.go +++ b/usr/gri/bignum/bignum_test.go @@ -9,28 +9,48 @@ import Big "bignum" const ( sa = "991"; sb = "2432902008176640000"; // 20! - sc = "93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000"; // 100! + sc = "933262154439441526816992388562667004907159682643816214685929" + "638952175999932299156089414639761565182862536979208272237582" + "51185210916864000000000000000000000000"; // 100! + sp = "170141183460469231731687303715884105727"; // prime ) var ( - a = Big.NatFromString(sa, 10); - b = Big.NatFromString(sb, 10); - c = Big.NatFromString(sc, 10); + nat_zero = Big.Nat(0); + nat_one = Big.Nat(1); + nat_two = Big.Nat(2); + + a = Big.NatFromString(sa, 10, nil); + b = Big.NatFromString(sb, 10, nil); + c = Big.NatFromString(sc, 10, nil); + p = Big.NatFromString(sp, 10, nil); + + int_zero = Big.Int(0); + int_one = Big.Int(1); + int_two = Big.Int(2); + + ip = Big.IntFromString(sp, 10, nil); + + rat_zero = Big.Rat(0, 1); + rat_half = Big.Rat(1, 2); + rat_one = Big.Rat(1, 1); + rat_two = Big.Rat(2, 1); ) var test_msg string; func TEST(n uint, b bool) { if !b { - panic("TEST failed: ", test_msg, "(", n, ")\n"); + println("TEST failed: ", test_msg, "(", n, ")"); + panic(); } } -func TEST_EQ(n uint, x, y *Big.Natural) { +func NAT_EQ(n uint, x, y *Big.Natural) { if x.Cmp(y) != 0 { - println("TEST failed:", test_msg, "(", n, ")\n"); + println("TEST failed:", test_msg, "(", n, ")"); println("x =", x.String(10)); println("y =", y.String(10)); panic(); @@ -38,180 +58,415 @@ func TEST_EQ(n uint, x, y *Big.Natural) { } -func TestLog2() { - test_msg = "TestLog2A"; - TEST(0, Big.Nat(1).Log2() == 0); - TEST(1, Big.Nat(2).Log2() == 1); - TEST(2, Big.Nat(3).Log2() == 1); - TEST(3, Big.Nat(4).Log2() == 2); - - test_msg = "TestLog2B"; - for i := uint(0); i < 100; i++ { - TEST(i, Big.Nat(1).Shl(i).Log2() == int(i)); +func INT_EQ(n uint, x, y *Big.Integer) { + if x.Cmp(y) != 0 { + println("TEST failed:", test_msg, "(", n, ")"); + println("x =", x.String(10)); + println("y =", y.String(10)); + panic(); } } -func TestConv() { - test_msg = "TestConvA"; - TEST(0, a.Cmp(Big.Nat(991)) == 0); - TEST(1, b.Cmp(Big.Fact(20)) == 0); - TEST(2, c.Cmp(Big.Fact(100)) == 0); +func RAT_EQ(n uint, x, y *Big.Rational) { + if x.Cmp(y) != 0 { + println("TEST failed:", test_msg, "(", n, ")"); + println("x =", x.String(10)); + println("y =", y.String(10)); + panic(); + } +} + + +func NatConv() { + test_msg = "NatConvA"; + NAT_EQ(0, a, Big.Nat(991)); + NAT_EQ(1, b, Big.Fact(20)); + NAT_EQ(2, c, Big.Fact(100)); TEST(3, a.String(10) == sa); TEST(4, b.String(10) == sb); TEST(5, c.String(10) == sc); - test_msg = "TestConvB"; + test_msg = "NatConvB"; + var slen int; + NAT_EQ(0, Big.NatFromString("0", 0, nil), nat_zero); + NAT_EQ(1, Big.NatFromString("123", 0, nil), Big.Nat(123)); + NAT_EQ(2, Big.NatFromString("077", 0, nil), Big.Nat(7*8 + 7)); + NAT_EQ(3, Big.NatFromString("0x1f", 0, nil), Big.Nat(1*16 + 15)); + NAT_EQ(4, Big.NatFromString("0x1fg", 0, &slen), Big.Nat(1*16 + 15)); + TEST(4, slen == 4); + + test_msg = "NatConvC"; t := c.Mul(c); for base := uint(2); base <= 16; base++ { - TEST_EQ(base, Big.NatFromString(t.String(base), base), t); + NAT_EQ(base, Big.NatFromString(t.String(base), base, nil), t); + } +} + + +func IntConv() { + test_msg = "IntConv"; + var slen int; + INT_EQ(0, Big.IntFromString("0", 0, nil), int_zero); + INT_EQ(1, Big.IntFromString("-0", 0, nil), int_zero); + INT_EQ(2, Big.IntFromString("123", 0, nil), Big.Int(123)); + INT_EQ(3, Big.IntFromString("-123", 0, nil), Big.Int(-123)); + INT_EQ(4, Big.IntFromString("077", 0, nil), Big.Int(7*8 + 7)); + INT_EQ(5, Big.IntFromString("-077", 0, nil), Big.Int(-(7*8 + 7))); + INT_EQ(6, Big.IntFromString("0x1f", 0, nil), Big.Int(1*16 + 15)); + INT_EQ(7, Big.IntFromString("-0x1f", 0, nil), Big.Int(-(1*16 + 15))); + INT_EQ(8, Big.IntFromString("0x1fg", 0, &slen), Big.Int(1*16 + 15)); + INT_EQ(9, Big.IntFromString("-0x1fg", 0, &slen), Big.Int(-(1*16 + 15))); + TEST(10, slen == 5); +} + + +func RatConv() { + test_msg = "RatConv"; + var slen int; + RAT_EQ(0, Big.RatFromString("0", 0, nil), rat_zero); + RAT_EQ(1, Big.RatFromString("0/", 0, nil), rat_zero); + RAT_EQ(2, Big.RatFromString("0/1", 0, nil), rat_zero); + RAT_EQ(3, Big.RatFromString("010/8", 0, nil), rat_one); + RAT_EQ(4, Big.RatFromString("20/0xa", 0, &slen), rat_two); + TEST(5, slen == 6); +} + + +func Add(x, y *Big.Natural) *Big.Natural { + z1 := x.Add(y); + z2 := y.Add(x); + if z1.Cmp(z2) != 0 { + println("addition not symmetric"); + println("x =", x.String(10)); + println("y =", y.String(10)); + panic(); } + return z1; } func Sum(n uint, scale *Big.Natural) *Big.Natural { - s := Big.Nat(0); + s := nat_zero; for ; n > 0; n-- { - s = s.Add(Big.Nat(uint64(n)).Mul(scale)); + s = Add(s, Big.Nat(n).Mul(scale)); } return s; } -func TestAdd() { - test_msg = "TestAddA"; +func NatAdd() { + test_msg = "NatAddA"; + NAT_EQ(0, Add(nat_zero, nat_zero), nat_zero); + NAT_EQ(1, Add(nat_zero, c), c); - test_msg = "TestAddB"; + test_msg = "NatAddB"; for i := uint(0); i < 100; i++ { - t := Big.Nat(uint64(i)); - TEST_EQ(i, Sum(i, c), t.Mul(t).Add(t).Shr(1).Mul(c)); + t := Big.Nat(i); + NAT_EQ(i, Sum(i, c), t.Mul(t).Add(t).Shr(1).Mul(c)); + } +} + + +func Mul(x, y *Big.Natural) *Big.Natural { + z1 := x.Mul(y); + z2 := y.Mul(x); + if z1.Cmp(z2) != 0 { + println("multiplication not symmetric"); + println("x =", x.String(10)); + println("y =", y.String(10)); + panic(); + } + if !x.IsZero() && z1.Div(x).Cmp(y) != 0 { + println("multiplication/division not inverse (A)"); + println("x =", x.String(10)); + println("y =", y.String(10)); + panic(); + } + if !y.IsZero() && z1.Div(y).Cmp(x) != 0 { + println("multiplication/division not inverse (B)"); + println("x =", x.String(10)); + println("y =", y.String(10)); + panic(); + } + return z1; +} + + +func NatSub() { + test_msg = "NatSubA"; + NAT_EQ(0, nat_zero.Sub(nat_zero), nat_zero); + NAT_EQ(1, c.Sub(nat_zero), c); + + test_msg = "NatSubB"; + for i := uint(0); i < 100; i++ { + t := Sum(i, c); + for j := uint(0); j <= i; j++ { + t = t.Sub(Mul(Big.Nat(j), c)); + } + NAT_EQ(i, t, nat_zero); + } +} + + +func NatMul() { + test_msg = "NatMulA"; + NAT_EQ(0, Mul(c, nat_zero), nat_zero); + NAT_EQ(1, Mul(c, nat_one), c); + + test_msg = "NatMulB"; + NAT_EQ(0, b.Mul(Big.MulRange(0, 100)), nat_zero); + NAT_EQ(1, b.Mul(Big.MulRange(21, 100)), c); + + test_msg = "NatMulC"; + const n = 100; + p := b.Mul(c).Shl(n); + for i := uint(0); i < n; i++ { + NAT_EQ(i, Mul(b.Shl(i), c.Shl(n-i)), p); + } +} + + +func NatDiv() { + test_msg = "NatDivA"; + NAT_EQ(0, c.Div(nat_one), c); + NAT_EQ(1, c.Div(Big.Nat(100)), Big.Fact(99)); + NAT_EQ(2, b.Div(c), nat_zero); + NAT_EQ(4, nat_one.Shl(100).Div(nat_one.Shl(90)), nat_one.Shl(10)); + NAT_EQ(5, c.Div(b), Big.MulRange(21, 100)); + + test_msg = "NatDivB"; + const n = 100; + p := Big.Fact(n); + for i := uint(0); i < n; i++ { + NAT_EQ(i, p.Div(Big.MulRange(1, i)), Big.MulRange(i+1, n)); + } +} + + +func IntQuoRem() { + test_msg = "IntQuoRem"; + type T struct { x, y, q, r int }; + a := []T{ + T{+8, +3, +2, +2}, + T{+8, -3, -2, +2}, + T{-8, +3, -2, -2}, + T{-8, -3, +2, -2}, + T{+1, +2, 0, +1}, + T{+1, -2, 0, +1}, + T{-1, +2, 0, -1}, + T{-1, -2, 0, -1}, + }; + for i := uint(0); i < len(a); i++ { + e := &a[i]; + x, y := Big.Int(e.x).Mul(ip), Big.Int(e.y).Mul(ip); + q, r := Big.Int(e.q), Big.Int(e.r).Mul(ip); + qq, rr := x.QuoRem(y); + INT_EQ(4*i+0, x.Quo(y), q); + INT_EQ(4*i+1, x.Rem(y), r); + INT_EQ(4*i+2, qq, q); + INT_EQ(4*i+3, rr, r); } } -func TestShift() { - test_msg = "TestShift1L"; +func IntDivMod() { + test_msg = "IntDivMod"; + type T struct { x, y, q, r int }; + a := []T{ + T{+8, +3, +2, +2}, + T{+8, -3, -2, +2}, + T{-8, +3, -3, +1}, + T{-8, -3, +3, +1}, + T{+1, +2, 0, +1}, + T{+1, -2, 0, +1}, + T{-1, +2, -1, +1}, + T{-1, -2, +1, +1}, + }; + for i := uint(0); i < len(a); i++ { + e := &a[i]; + x, y := Big.Int(e.x).Mul(ip), Big.Int(e.y).Mul(ip); + q, r := Big.Int(e.q), Big.Int(e.r).Mul(ip); + qq, rr := x.DivMod(y); + INT_EQ(4*i+0, x.Div(y), q); + INT_EQ(4*i+1, x.Mod(y), r); + INT_EQ(4*i+2, qq, q); + INT_EQ(4*i+3, rr, r); + } +} + + +func NatMod() { + test_msg = "NatModA"; + for i := uint(0); ; i++ { + d := nat_one.Shl(i); + if d.Cmp(c) < 0 { + NAT_EQ(i, c.Add(d).Mod(c), d); + } else { + NAT_EQ(i, c.Add(d).Div(c), nat_two); + NAT_EQ(i, c.Add(d).Mod(c), d.Sub(c)); + break; + } + } +} + + +func NatShift() { + test_msg = "NatShift1L"; TEST(0, b.Shl(0).Cmp(b) == 0); TEST(1, c.Shl(1).Cmp(c) > 0); - test_msg = "TestShift1R"; + test_msg = "NatShift1R"; TEST(0, b.Shr(0).Cmp(b) == 0); TEST(1, c.Shr(1).Cmp(c) < 0); - test_msg = "TestShift2"; + test_msg = "NatShift2"; for i := uint(0); i < 100; i++ { TEST(i, c.Shl(i).Shr(i).Cmp(c) == 0); } - test_msg = "TestShift3L"; + test_msg = "NatShift3L"; { const m = 3; p := b; f := Big.Nat(1< 0); - test_msg = "TestMulB"; - const n = 100; - p := b.Mul(c).Shl(n); - for i := uint(0); i < n; i++ { - TEST_EQ(i, b.Shl(i).Mul(c.Shl(n-i)), p); + test_msg = "IntShift1R"; + TEST(0, ip.Shr(0).Cmp(ip) == 0); + TEST(1, ip.Shr(1).Cmp(ip) < 0); + + test_msg = "IntShift2"; + for i := uint(0); i < 100; i++ { + TEST(i, ip.Shl(i).Shr(i).Cmp(ip) == 0); } -} + test_msg = "IntShift3L"; + { const m = 3; + p := ip; + f := Big.Int(1<> 1)); + //INT_EQ(1, ip.Neg().Shr(10), ip.Neg().Div(Big.Int(1).Shl(10))); } -func TestMod() { - test_msg = "TestModA"; - for i := uint(0); ; i++ { - d := Big.Nat(1).Shl(i); - if d.Cmp(c) < 0 { - TEST_EQ(i, c.Add(d).Mod(c), d); - } else { - TEST_EQ(i, c.Add(d).Div(c), Big.Nat(2)); - TEST_EQ(i, c.Add(d).Mod(c), d.Sub(c)); - break; - } +func NatCmp() { + test_msg = "NatCmp"; + TEST(0, a.Cmp(a) == 0); + TEST(1, a.Cmp(b) < 0); + TEST(2, b.Cmp(a) > 0); + TEST(3, a.Cmp(c) < 0); + d := c.Add(b); + TEST(4, c.Cmp(d) < 0); + TEST(5, d.Cmp(c) > 0); +} + + +func NatLog2() { + test_msg = "NatLog2A"; + TEST(0, nat_one.Log2() == 0); + TEST(1, nat_two.Log2() == 1); + TEST(2, Big.Nat(3).Log2() == 1); + TEST(3, Big.Nat(4).Log2() == 2); + + test_msg = "NatLog2B"; + for i := uint(0); i < 100; i++ { + TEST(i, nat_one.Shl(i).Log2() == i); } } -func TestGcd() { - test_msg = "TestGcdA"; +func NatGcd() { + test_msg = "NatGcdA"; f := Big.Nat(99991); - TEST_EQ(0, b.Mul(f).Gcd(c.Mul(f)), Big.MulRange(1, 20).Mul(f)); + NAT_EQ(0, b.Mul(f).Gcd(c.Mul(f)), Big.MulRange(1, 20).Mul(f)); } -func TestPow() { - test_msg = "TestPowA"; - TEST_EQ(0, Big.Nat(2).Pow(0), Big.Nat(1)); +func NatPow() { + test_msg = "NatPowA"; + NAT_EQ(0, nat_two.Pow(0), nat_one); - test_msg = "TestPowB"; + test_msg = "NatPowB"; for i := uint(0); i < 100; i++ { - TEST_EQ(i, Big.Nat(2).Pow(i), Big.Nat(1).Shl(i)); + NAT_EQ(i, nat_two.Pow(i), nat_one.Shl(i)); } } -func TestPop() { - test_msg = "TestPopA"; - TEST(0, Big.Nat(0).Pop() == 0); - TEST(1, Big.Nat(1).Pop() == 1); +func NatPop() { + test_msg = "NatPopA"; + TEST(0, nat_zero.Pop() == 0); + TEST(1, nat_one.Pop() == 1); TEST(2, Big.Nat(10).Pop() == 2); TEST(3, Big.Nat(30).Pop() == 4); TEST(4, Big.Nat(0x1248f).Shl(33).Pop() == 8); - test_msg = "TestPopB"; + test_msg = "NatPopB"; for i := uint(0); i < 100; i++ { - TEST(i, Big.Nat(1).Shl(i).Sub(Big.Nat(1)).Pop() == i); + TEST(i, nat_one.Shl(i).Sub(nat_one).Pop() == i); } } func main() { - TestLog2(); - TestConv(); - TestAdd(); - TestShift(); - TestMul(); - TestDiv(); - TestMod(); - TestGcd(); - TestPow(); - TestPop(); + // Naturals + NatConv(); + NatAdd(); + NatSub(); + NatMul(); + NatDiv(); + NatMod(); + NatShift(); + NatCmp(); + NatLog2(); + NatGcd(); + NatPow(); + NatPop(); + + // Integers + // TODO add more tests + IntConv(); + IntQuoRem(); + IntDivMod(); + IntShift(); + + // Rationals + // TODO add more tests + RatConv(); + print("PASSED\n"); }