Skip to content

G
golang
Commit 7cd11c1c authored by Robert Griesemer's avatar Robert Griesemer
Browse files
Options

- completed integer support (some logical functions missing) 

- completed rational support
- better documentation
- more tests
- cleanups

R=r
OCL=18438
CL=18438
parent e5d9a5c9
Show whitespace changes
Inline Side-by-side
Showing with 985 additions and 462 deletions
usr/gri/bignum/bignum.go
View file @ 7cd11c1c
...@@ -13,7 +13,7 @@ package Bignum ...@@ -13,7 +13,7 @@ package Bignum
// ---------------------------------------------------------------------------- // ----------------------------------------------------------------------------
// Representation // Internal representation
// //
// A natural number of the form // A natural number of the form
// //
...@@ -27,37 +27,56 @@ package Bignum ...@@ -27,37 +27,56 @@ package Bignum
// always normalized before returning the final result. The normalized // always normalized before returning the final result. The normalized
// representation of 0 is the empty array (length = 0). // 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 // 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 // are a few constraints besides the size of the largest unsigned integer
// type available. // type available:
// TODO describe the constraints. //
// 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 LogW = 64;
const LogH = 4; // bits for a hex digit (= "small" number) const LogH = 4; // bits for a hex digit (= "small" number)
const LogB = LogW - LogH; const LogB = LogW - LogH; // largest bit-width available
const ( const (
L2 = LogB / 2; // half-digits
B2 = 1 << L2; W2 = LogB / 2; // width
M2 = B2 - 1; B2 = 1 << W2; // base
M2 = B2 - 1; // mask
L = L2 * 2;
B = 1 << L; // full digits
M = B - 1; W = W2 * 2; // width
) B = 1 << W; // base
M = B - 1; // mask
type (
Digit2 uint32;
Digit uint64;
) )
// ---------------------------------------------------------------------------- // ----------------------------------------------------------------------------
// Support // Support functions
// TODO replace this with a Go built-in assert
func assert(p bool) { func assert(p bool) {
if !p { if !p {
panic("assert failed"); panic("assert failed");
...@@ -65,56 +84,36 @@ func assert(p bool) { ...@@ -65,56 +84,36 @@ func assert(p bool) {
} }
// ---------------------------------------------------------------------------- func IsSmall(x Digit) bool {
// Raw operations return x < 1<<LogH;
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 Xor1(z, x *[]Digit, y Digit) { export func Dump(x *[]Digit) {
print("[", len(x), "]");
for i := len(x) - 1; i >= 0; i-- { for i := len(x) - 1; i >= 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-- { // Raw operations on sequences of digits
z[i] = x[i] ^ y[i]; //
} // Naming conventions
} //
// c carry
// x, y operands
// z result
// n, m len(x), len(y)
func Add1(z, x *[]Digit, c Digit) Digit { func Add1(z, x *[]Digit, c Digit) Digit {
n := len(x); n := len(x);
for i := 0; i < n; i++ { for i := 0; i < n; i++ {
t := c + x[i]; t := c + x[i];
c, z[i] = t>>L, t&M c, z[i] = t>>W, t&M
} }
return c; return c;
} }
...@@ -125,7 +124,7 @@ func Add(z, x, y *[]Digit) Digit { ...@@ -125,7 +124,7 @@ func Add(z, x, y *[]Digit) Digit {
n := len(x); n := len(x);
for i := 0; i < n; i++ { for i := 0; i < n; i++ {
t := c + x[i] + y[i]; t := c + x[i] + y[i];
c, z[i] = t>>L, t&M c, z[i] = t>>W, t&M
} }
return c; return c;
} }
...@@ -135,7 +134,7 @@ func Sub1(z, x *[]Digit, c Digit) Digit { ...@@ -135,7 +134,7 @@ func Sub1(z, x *[]Digit, c Digit) Digit {
n := len(x); n := len(x);
for i := 0; i < n; i++ { for i := 0; i < n; i++ {
t := c + x[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; return c;
} }
...@@ -146,7 +145,7 @@ func Sub(z, x, y *[]Digit) Digit { ...@@ -146,7 +145,7 @@ func Sub(z, x, y *[]Digit) Digit {
n := len(x); n := len(x);
for i := 0; i < n; i++ { for i := 0; i < n; i++ {
t := c + x[i] - y[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; return c;
} }
...@@ -154,31 +153,33 @@ func Sub(z, x, y *[]Digit) Digit { ...@@ -154,31 +153,33 @@ func Sub(z, x, y *[]Digit) Digit {
// Returns c = x*y div B, z = x*y mod B. // Returns c = x*y div B, z = x*y mod B.
func Mul11(x, y Digit) (Digit, Digit) { 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, // multiply the digits separately while avoiding overflow,
// and return the product as two separate digits. // and return the product as two separate digits.
const L0 = (L + 1)/2; // This code also works for non-even bit widths W
const L1 = L - L0; // which is why there are separate constants below
const DL = L0 - L1; // 0 or 1 // for half-digits.
const b = 1<<L0; const W2 = (W + 1)/2;
const m = b - 1; const DW = W2*2 - W; // 0 or 1
const B2 = 1<<W2;
const M2 = B2 - 1;
// split x and y into sub-digits // split x and y into sub-digits
// x = (x1*b + x0) // x = (x1*B2 + x0)
// y = (y1*b + y0) // y = (y1*B2 + y0)
x1, x0 := x>>L0, x&m; x1, x0 := x>>W2, x&M2;
y1, y0 := y>>L0, y&m; 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; t0 := x0*y0;
t1 := x1*y0 + x0*y1; t1 := x1*y0 + x0*y1;
t2 := x1*y1; t2 := x1*y1;
// compute the result digits but avoid overflow // compute the result digits but avoid overflow
// z = z1*B + z0 = x*y // z = z1*B + z0 = x*y
z0 := (t1<<L0 + t0)&M; z0 := (t1<<W2 + t0)&M;
z1 := t2<<DL + (t1 + t0>>L0)>>L1; z1 := t2<<DW + (t1 + t0>>W2)>>(W-W2);
return z1, z0; return z1, z0;
} }
...@@ -195,7 +196,7 @@ func Mul(z, x, y *[]Digit) { ...@@ -195,7 +196,7 @@ func Mul(z, x, y *[]Digit) {
// z[i+j] += c + x[i]*d; // z[i+j] += c + x[i]*d;
z1, z0 := Mul11(x[i], d); z1, z0 := Mul11(x[i], d);
t := c + z[i+j] + z0; t := c + z[i+j] + z0;
c, z[i+j] = t>>L, t&M; c, z[i+j] = t>>W, t&M;
c += z1; c += z1;
} }
z[n+j] = c; z[n+j] = c;
...@@ -204,118 +205,113 @@ func Mul(z, x, y *[]Digit) { ...@@ -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); n := len(x);
var c Digit; var c Digit;
f := Digit(y);
for i := 0; i < n; i++ { for i := 0; i < n; i++ {
t := c + Digit(x[i])*f; c, z[i] = x[i] >> (W-s), x[i] << s & M | c;
c, z[i] = t>>L2, Digit2(t&M2);
} }
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); n := len(x);
var c Digit; var c Digit;
d := Digit(y); for i := n - 1; i >= 0; i-- {
for i := n-1; i >= 0; i-- { c, z[i] = x[i] << (W-s) & M, x[i] >> s | c;
t := c*B2 + Digit(x[i]);
c, z[i] = t%d, Digit2(t/d);
} }
return Digit2(c); return c;
} }
func Shl(z, x *[]Digit, s uint) Digit { func And1(z, x *[]Digit, y Digit) {
assert(s <= L); for i := len(x) - 1; i >= 0; i-- {
n := len(x); z[i] = x[i] & y;
var c Digit;
for i := 0; i < n; i++ {
c, z[i] = x[i] >> (L-s), x[i] << s & M | c;
} }
return c;
} }
func Shr(z, x *[]Digit, s uint) Digit { func And(z, x, y *[]Digit) {
assert(s <= L); for i := len(x) - 1; i >= 0; i-- {
n := len(x); z[i] = x[i] & y[i];
var c Digit;
for i := n - 1; i >= 0; i-- {
c, z[i] = x[i] << (L-s) & M, x[i] >> s | c;
} }
return c;
} }
// ---------------------------------------------------------------------------- func Or1(z, x *[]Digit, y Digit) {
// Support for i := len(x) - 1; i >= 0; i-- {
z[i] = x[i] | y;
}
}
func IsSmall(x Digit) bool {
return x < 1<<LogH; func Or(z, x, y *[]Digit) {
for i := len(x) - 1; i >= 0; i-- {
z[i] = x[i] | y[i];
}
} }
func Split(x Digit) (Digit, Digit) { func Xor1(z, x *[]Digit, y Digit) {
return x>>L, x&M; for i := len(x) - 1; i >= 0; i-- {
z[i] = x[i] ^ y;
}
} }
export func Dump(x *[]Digit) { func Xor(z, x, y *[]Digit) {
print("[", len(x), "]");
for i := len(x) - 1; i >= 0; i-- { for i := len(x) - 1; i >= 0; i-- {
print(" ", x[i]); z[i] = x[i] ^ y[i];
} }
println();
} }
// ---------------------------------------------------------------------------- // ----------------------------------------------------------------------------
// Natural numbers // 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 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; // Creation
switch {
case x == 0: export func Nat(x uint) *Natural {
z = NatZero; switch x {
case x < B: case 0: return NatZero;
z = new(Natural, 1); case 1: return NatOne;
z[0] = x; case 2: return NatTwo;
return z; case 10: return NatTen;
default:
z = new(Natural, 2);
z[1], z[0] = Split(x);
} }
return z; assert(Digit(x) < B);
return &Natural{Digit(x)};
} }
func Normalize(x *Natural) *Natural { // Predicates
n := len(x);
for n > 0 && x[n - 1] == 0 { n-- } func (x *Natural) IsOdd() bool {
if n < len(x) { return len(x) > 0 && x[0]&1 != 0;
x = x[0 : n]; // trim leading 0's
}
return x;
} }
func Normalize2(x *[]Digit2) *[]Digit2 { func (x *Natural) IsZero() bool {
return len(x) == 0;
}
// Operations
func Normalize(x *Natural) *Natural {
n := len(x); n := len(x);
for n > 0 && x[n - 1] == 0 { n-- } for n > 0 && x[n - 1] == 0 { n-- }
if n < len(x) { if n < len(x) {
...@@ -325,12 +321,6 @@ func Normalize2(x *[]Digit2) *[]Digit2 { ...@@ -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 { func (x *Natural) Add(y *Natural) *Natural {
n := len(x); n := len(x);
m := len(y); m := len(y);
...@@ -363,19 +353,6 @@ func (x *Natural) Sub(y *Natural) *Natural { ...@@ -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 { func (x *Natural) Mul(y *Natural) *Natural {
n := len(x); n := len(x);
m := len(y); m := len(y);
...@@ -387,77 +364,24 @@ func (x *Natural) Mul(y *Natural) *Natural { ...@@ -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 // DivMod needs multi-precision division which is not available if Digit
// is already using the largest uint size. Split base before division, // is already using the largest uint size. Instead, unpack each operand
// and merge again after. Each Digit is split into 2 Digit2's. // 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 { func Unpack(x *Natural) *[]Digit2 {
// TODO Use Log() for better result - don't need Normalize2 at the end!
n := len(x); n := len(x);
z := new([]Digit2, n*2 + 1); // add space for extra digit (used by DivMod) z := new([]Digit2, n*2 + 1); // add space for extra digit (used by DivMod)
for i := 0; i < n; i++ { for i := 0; i < n; i++ {
t := x[i]; t := x[i];
z[i*2] = Digit2(t & M2); 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 { ...@@ -470,34 +394,65 @@ func Pack(x *[]Digit2) *Natural {
z[n] = Digit(x[n*2]); z[n] = Digit(x[n*2]);
} }
for i := 0; i < n; i++ { 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); return Normalize(z);
} }
// Division and modulo computation - destroys x and y. Based on the func Mul1(z, x *[]Digit2, y Digit2) Digit2 {
// algorithms described in: 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 // 1) D. Knuth, "The Art of Computer Programming. Volume 2. Seminumerical
// Algorithms." Addison-Wesley, Reading, 1969. // Algorithms." Addison-Wesley, Reading, 1969.
// (Algorithm D, Sec. 4.3.1)
// //
// 2) P. Brinch Hansen, Multiple-length division revisited: A tour of the // 2) Henry S. Warren, Jr., "A Hacker's Delight". Addison-Wesley, 2003.
// (9-2 Multiword Division, p.140ff)
//
// 3) P. Brinch Hansen, Multiple-length division revisited: A tour of the
// minefield. "Software - Practice and Experience 24", (June 1994), // minefield. "Software - Practice and Experience 24", (June 1994),
// 579-601. John Wiley & Sons, Ltd. // 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); n := len(x);
m := len(y); m := len(y);
assert(m > 0); // division by zero if m == 0 {
assert(n+1 <= cap(x)); // space for one extra digit (should it be == ?) panic("division by zero");
}
assert(n+1 <= cap(x)); // space for one extra digit
x = x[0 : n + 1]; x = x[0 : n + 1];
assert(x[n] == 0);
if m == 1 { if m == 1 {
// division by single digit // division by single digit
...@@ -505,36 +460,39 @@ func DivMod2(x, y *[]Digit2) (*[]Digit2, *[]Digit2) { ...@@ -505,36 +460,39 @@ func DivMod2(x, y *[]Digit2) (*[]Digit2, *[]Digit2) {
x[0] = Div1(x[1 : n+1], x[0 : n], y[0]); x[0] = Div1(x[1 : n+1], x[0 : n], y[0]);
} else if m > n { } else if m > n {
// quotient = 0, remainder = x // y > x => quotient = 0, remainder = x
// TODO in this case we shouldn't even split base - FIX THIS // TODO in this case we shouldn't even unpack x and y
m = n; m = n;
} else { } else {
// general case // general case
assert(2 <= m && m <= n); assert(2 <= m && m <= n);
assert(x[n] == 0);
// normalize x and y // normalize x and y
f := b/(Digit(y[m-1]) + 1); // 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(x, x, Digit2(f));
Mul1(y, y, Digit2(f)); Mul1(y, y, Digit2(f));
assert(b/2 <= y[m-1] && y[m-1] < b); // incorrect scaling }
assert(B2/2 <= y[m-1] && y[m-1] < B2); // incorrect scaling
y1, y2 := Digit(y[m-1]), Digit(y[m-2]); y1, y2 := Digit(y[m-1]), Digit(y[m-2]);
d2 := Digit(y1)*b + Digit(y2); d2 := Digit(y1)<<W2 + Digit(y2);
for i := n-m; i >= 0; i-- { for i := n-m; i >= 0; i-- {
k := i+m; k := i+m;
// compute trial digit // compute trial digit (Knuth)
var q Digit; 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 { if x0 != y1 {
q = (x0*b + x1)/y1; q = (x0<<W2 + x1)/y1;
} else { } else {
q = b-1; q = B2 - 1;
} }
for y2 * q > (x0*b + x1 - y1*q)*b + x2 { for y2*q > (x0<<W2 + x1 - y1*q)<<W2 + x2 {
q-- q--
} }
} }
...@@ -542,8 +500,8 @@ func DivMod2(x, y *[]Digit2) (*[]Digit2, *[]Digit2) { ...@@ -542,8 +500,8 @@ func DivMod2(x, y *[]Digit2) (*[]Digit2, *[]Digit2) {
// subtract y*q // subtract y*q
c := Digit(0); c := Digit(0);
for j := 0; j < m; j++ { for j := 0; j < m; j++ {
t := c + Digit(x[i+j]) - Digit(y[j])*q; // arithmetic shift! t := c + Digit(x[i+j]) - Digit(y[j])*q;
c, x[i+j] = Digit(int64(t)>>L2), Digit2(t&M2); c, x[i+j] = Digit(int64(t)>>W2), Digit2(t&M2); // requires arithmetic shift!
} }
// correct if trial digit was too large // correct if trial digit was too large
...@@ -552,7 +510,7 @@ func DivMod2(x, y *[]Digit2) (*[]Digit2, *[]Digit2) { ...@@ -552,7 +510,7 @@ func DivMod2(x, y *[]Digit2) (*[]Digit2, *[]Digit2) {
c := Digit(0); c := Digit(0);
for j := 0; j < m; j++ { for j := 0; j < m; j++ {
t := c + Digit(x[i+j]) + Digit(y[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); assert(c + Digit(x[k]) == 0);
// correct trial digit // correct trial digit
...@@ -563,68 +521,56 @@ func DivMod2(x, y *[]Digit2) (*[]Digit2, *[]Digit2) { ...@@ -563,68 +521,56 @@ func DivMod2(x, y *[]Digit2) (*[]Digit2, *[]Digit2) {
} }
// undo normalization for remainder // undo normalization for remainder
if f != 1 {
c := Div1(x[0 : m], x[0 : m], Digit2(f)); c := Div1(x[0 : m], x[0 : m], Digit2(f));
assert(c == 0); assert(c == 0);
} }
}
return x[m : n+1], x[0 : m]; return x[m : n+1], x[0 : m];
} }
func (x *Natural) Div(y *Natural) *Natural { func (x *Natural) Div(y *Natural) *Natural {
q, r := DivMod2(Unpack(x), Unpack(y)); q, r := DivMod(Unpack(x), Unpack(y));
return Pack(q); return Pack(q);
} }
func (x *Natural) Mod(y *Natural) *Natural { func (x *Natural) Mod(y *Natural) *Natural {
q, r := DivMod2(Unpack(x), Unpack(y)); q, r := DivMod(Unpack(x), Unpack(y));
return Pack(r); return Pack(r);
} }
func (x *Natural) DivMod(y *Natural) (*Natural, *Natural) { 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); return Pack(q), Pack(r);
} }
func (x *Natural) Cmp(y *Natural) int { func (x *Natural) Shl(s uint) *Natural {
n := len(x); n := uint(len(x));
m := len(y); m := n + s/W;
z := new(Natural, m+1);
if n != m || n == 0 { z[m] = Shl(z[m-n : m], x, s%W);
return n - m;
}
i := n - 1; return Normalize(z);
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) Shr(s uint) *Natural {
n := uint(len(x));
m := n - s/W;
if m > n { // check for underflow
m = 0;
}
z := new(Natural, m);
func Log2(x Digit) int { Shr(z, x[n-m : n], s%W);
n := -1;
for x != 0 { x = x >> 1; n++; } // BUG >>= broken for uint64
return n;
}
func (x *Natural) Log2() int { return Normalize(z);
n := len(x);
if n > 0 {
n = (n - 1)*L + Log2(x[n - 1]);
} else {
n = -1;
}
return n;
} }
...@@ -673,14 +619,55 @@ func (x *Natural) Xor(y *Natural) *Natural { ...@@ -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. // 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)); assert(0 < d && IsSmall(d - 1));
c := Digit(0); c := Digit(0);
for i := len(x) - 1; i >= 0; i-- { for i := len(x) - 1; i >= 0; i-- {
t := c<<L + x[i]; t := c<<W + x[i];
c, x[i] = t%d, t/d; c, x[i] = t%d, t/d;
} }
...@@ -689,26 +676,25 @@ func (x *Natural) DivMod1(d Digit) (*Natural, Digit) { ...@@ -689,26 +676,25 @@ func (x *Natural) DivMod1(d Digit) (*Natural, Digit) {
func (x *Natural) String(base uint) string { func (x *Natural) String(base uint) string {
if x.IsZero() { if len(x) == 0 {
return "0"; return "0";
} }
// allocate string // allocate buffer for conversion
assert(2 <= base && base <= 16); assert(2 <= base && base <= 16);
n := (x.Log2() + 1) / Log2(Digit(base)) + 1; // TODO why the +1? n := (x.Log2() + 1) / Log2(Digit(base)) + 1; // +1: round up
s := new([]byte, n); s := new([]byte, n);
// convert // don't destroy x
// don't destroy x, make a copy
t := new(Natural, len(x)); t := new(Natural, len(x));
Or1(t, x, 0); // copy x Or1(t, x, 0); // copy
// convert
i := n; i := n;
for !t.IsZero() { for !t.IsZero() {
i--; i--;
var d Digit; var d Digit;
t, d = t.DivMod1(Digit(base)); t, d = DivMod1(t, Digit(base));
s[i] = "0123456789abcdef"[d]; s[i] = "0123456789abcdef"[d];
}; };
...@@ -716,7 +702,104 @@ func (x *Natural) String(base uint) string { ...@@ -716,7 +702,104 @@ func (x *Natural) String(base uint) string {
} }
export func MulRange(a, b Digit) *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;
}
return d;
}
// Computes x = x*d + c for "small" d's.
func MulAdd1(x *Natural, d, c Digit) *Natural {
assert(IsSmall(d-1) && IsSmall(c));
n := len(x);
z := new(Natural, n + 1);
for i := 0; i < n; i++ {
t := c + x[i]*d;
c, z[i] = t>>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 { switch {
case a > b: return Nat(1); case a > b: return Nat(1);
case a == b: return Nat(a); case a == b: return Nat(a);
...@@ -728,7 +811,7 @@ export func MulRange(a, b Digit) *Natural { ...@@ -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 // Using MulRange() instead of the basic for-loop
// lead to faster factorial computation. // lead to faster factorial computation.
return MulRange(2, n); return MulRange(2, n);
...@@ -744,60 +827,73 @@ func (x *Natural) Gcd(y *Natural) *Natural { ...@@ -744,60 +827,73 @@ func (x *Natural) Gcd(y *Natural) *Natural {
} }
func HexValue(ch byte) uint { // ----------------------------------------------------------------------------
d := uint(1 << LogH); // Integer numbers
switch { //
case '0' <= ch && ch <= '9': d = uint(ch - '0'); // Integers are normalized if the mantissa is normalized and the sign is
case 'a' <= ch && ch <= 'f': d = uint(ch - 'a') + 10; // false for mant == 0. Use MakeInt to create normalized Integers.
case 'A' <= ch && ch <= 'F': d = uint(ch - 'A') + 10;
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 Int(x int) *Integer {
export func NatFromString(s string, base uint) *Natural { sign := false;
x := NatZero; var ux uint;
for i := 0; i < len(s); i++ { if x < 0 {
d := HexValue(s[i]); sign = true;
if d < base { if -x == x {
x = x.MulAdd1(Digit(base), Digit(d)); // smallest negative integer
t := ^0;
ux = ^(uint(t) >> 1);
} else { } else {
break; ux = uint(-x);
} }
} else {
ux = uint(x);
} }
return x; return MakeInt(sign, Nat(ux));
} }
// ---------------------------------------------------------------------------- // Predicates
// Algorithms
export type T interface { func (x *Integer) IsOdd() bool {
IsZero() bool; return x.mant.IsOdd();
Mod(y T) bool;
} }
export func Gcd(x, y T) T {
// Euclidean algorithm. func (x *Integer) IsZero() bool {
for !y.IsZero() { return x.mant.IsZero();
x, y = y, x.Mod(y);
}
return x;
} }
// ---------------------------------------------------------------------------- func (x *Integer) IsNeg() bool {
// Integer numbers return x.sign && !x.mant.IsZero()
}
export type Integer struct {
sign bool; func (x *Integer) IsPos() bool {
mant *Natural; return !x.sign && !x.mant.IsZero()
} }
export func Int(x int64) *Integer { // Operations
return nil;
func (x *Integer) Neg() *Integer {
return MakeInt(!x.sign, x.mant);
} }
...@@ -806,14 +902,14 @@ func (x *Integer) Add(y *Integer) *Integer { ...@@ -806,14 +902,14 @@ func (x *Integer) Add(y *Integer) *Integer {
if x.sign == y.sign { if x.sign == y.sign {
// x + y == x + y // x + y == x + y
// (-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 { } else {
// x + (-y) == x - y == -(y - x) // x + (-y) == x - y == -(y - x)
// (-x) + y == y - x == -(x - y) // (-x) + y == y - x == -(x - y)
if x.mant.Cmp(y.mant) >= 0 { if x.mant.Cmp(y.mant) >= 0 {
z = &Integer{false, x.mant.Sub(y.mant)}; z = MakeInt(false, x.mant.Sub(y.mant));
} else { } else {
z = &Integer{true, y.mant.Sub(x.mant)}; z = MakeInt(true, y.mant.Sub(x.mant));
} }
} }
if x.sign { if x.sign {
...@@ -828,14 +924,14 @@ func (x *Integer) Sub(y *Integer) *Integer { ...@@ -828,14 +924,14 @@ func (x *Integer) Sub(y *Integer) *Integer {
if x.sign != y.sign { if x.sign != y.sign {
// x - (-y) == x + y // x - (-y) == x + y
// (-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 { } else {
// x - y == x - y == -(y - x) // x - y == x - y == -(y - x)
// (-x) - (-y) == y - x == -(x - y) // (-x) - (-y) == y - x == -(x - y)
if x.mant.Cmp(y.mant) >= 0 { if x.mant.Cmp(y.mant) >= 0 {
z = &Integer{false, x.mant.Sub(y.mant)}; z = MakeInt(false, x.mant.Sub(y.mant));
} else { } else {
z = &Integer{true, y.mant.Sub(x.mant)}; z = MakeInt(true, y.mant.Sub(x.mant));
} }
} }
if x.sign { if x.sign {
...@@ -850,16 +946,30 @@ func (x *Integer) Mul(y *Integer) *Integer { ...@@ -850,16 +946,30 @@ func (x *Integer) Mul(y *Integer) *Integer {
// x * (-y) == -(x * y) // x * (-y) == -(x * y)
// (-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)}; 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 { func (x *Integer) Quo(y *Integer) *Integer {
// x / y == x / y // x / y == x / y
// x / (-y) == -(x / y) // x / (-y) == -(x / y)
// (-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 { ...@@ -868,31 +978,116 @@ func (x *Integer) Rem(y *Integer) *Integer {
// x % (-y) == x % y // x % (-y) == x % y
// (-x) % y == -(x % y) // (-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) { func (x *Integer) QuoRem(y *Integer) (*Integer, *Integer) {
q, r := x.mant.DivMod(y.mant); 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 { func (x *Integer) Div(y *Integer) *Integer {
q, r := x.mant.DivMod(y.mant); q, r := x.QuoRem(y);
return nil; 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 { 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"); panic("UNIMPLEMENTED");
return nil; 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"); 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 { ...@@ -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 // get sign, if any
sign := false; sign := false;
if len(s) > 0 && (s[0] == '-' || s[0] == '+') { if len(s) > 0 && (s[0] == '-' || s[0] == '+') {
sign = 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 { ...@@ -922,56 +1127,119 @@ export func IntFromString(s string, base uint) *Integer {
// Rational numbers // Rational numbers
export type Rational struct { export type Rational struct {
a, b *Integer; // a = numerator, b = denominator a *Integer; // numerator
b *Natural; // denominator
} }
func (x *Rational) Normalize() *Rational { // Creation
f := x.a.mant.Gcd(x.b.mant);
x.a.mant = x.a.mant.Div(f); export func MakeRat(a *Integer, b *Natural) *Rational {
x.b.mant = x.b.mant.Div(f); f := a.mant.Gcd(b); // f > 0
return x; 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 { // Predicates
return (&Rational{a, b}).Normalize();
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 { 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 { 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 { 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 { func (x *Rational) Quo(y *Rational) *Rational {
return Rat(x.a.Mul(y.b), x.b.Mul(y.a)); 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 { func (x *Rational) Cmp(y *Rational) int {
panic("UNIMPLEMENTED"); return (x.a.MulNat(y.b)).Cmp(y.a.MulNat(x.b));
return nil;
} }
func (x *Rational) Cmp(y *Rational) int { func (x *Rational) String(base uint) string {
panic("UNIMPLEMENTED"); s := x.a.String(base);
return 0; if !x.IsInt() {
s += "/" + x.b.String(base);
}
return s;
} }
export func RatFromString(s string) *Rational { // Determines base (octal, decimal, hexadecimal) if base == 0.
panic("UNIMPLEMENTED"); export func RatFromString(s string, base uint, slen *int) *Rational {
return nil; // 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);
} }
usr/gri/bignum/bignum_test.go
View file @ 7cd11c1c
...@@ -9,28 +9,48 @@ import Big "bignum" ...@@ -9,28 +9,48 @@ import Big "bignum"
const ( const (
sa = "991"; sa = "991";
sb = "2432902008176640000"; // 20! sb = "2432902008176640000"; // 20!
sc = "93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000"; // 100! sc = "933262154439441526816992388562667004907159682643816214685929"
"638952175999932299156089414639761565182862536979208272237582"
"51185210916864000000000000000000000000"; // 100!
sp = "170141183460469231731687303715884105727"; // prime
) )
var ( var (
a = Big.NatFromString(sa, 10); nat_zero = Big.Nat(0);
b = Big.NatFromString(sb, 10); nat_one = Big.Nat(1);
c = Big.NatFromString(sc, 10); 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; var test_msg string;
func TEST(n uint, b bool) { func TEST(n uint, b bool) {
if !b { 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 { if x.Cmp(y) != 0 {
println("TEST failed:", test_msg, "(", n, ")\n"); println("TEST failed:", test_msg, "(", n, ")");
println("x =", x.String(10)); println("x =", x.String(10));
println("y =", y.String(10)); println("y =", y.String(10));
panic(); panic();
...@@ -38,180 +58,415 @@ func TEST_EQ(n uint, x, y *Big.Natural) { ...@@ -38,180 +58,415 @@ func TEST_EQ(n uint, x, y *Big.Natural) {
} }
func TestLog2() { func INT_EQ(n uint, x, y *Big.Integer) {
test_msg = "TestLog2A"; if x.Cmp(y) != 0 {
TEST(0, Big.Nat(1).Log2() == 0); println("TEST failed:", test_msg, "(", n, ")");
TEST(1, Big.Nat(2).Log2() == 1); println("x =", x.String(10));
TEST(2, Big.Nat(3).Log2() == 1); println("y =", y.String(10));
TEST(3, Big.Nat(4).Log2() == 2); panic();
}
}
test_msg = "TestLog2B";
for i := uint(0); i < 100; i++ { func RAT_EQ(n uint, x, y *Big.Rational) {
TEST(i, Big.Nat(1).Shl(i).Log2() == int(i)); if x.Cmp(y) != 0 {
println("TEST failed:", test_msg, "(", n, ")");
println("x =", x.String(10));
println("y =", y.String(10));
panic();
} }
} }
func TestConv() { func NatConv() {
test_msg = "TestConvA"; test_msg = "NatConvA";
TEST(0, a.Cmp(Big.Nat(991)) == 0); NAT_EQ(0, a, Big.Nat(991));
TEST(1, b.Cmp(Big.Fact(20)) == 0); NAT_EQ(1, b, Big.Fact(20));
TEST(2, c.Cmp(Big.Fact(100)) == 0); NAT_EQ(2, c, Big.Fact(100));
TEST(3, a.String(10) == sa); TEST(3, a.String(10) == sa);
TEST(4, b.String(10) == sb); TEST(4, b.String(10) == sb);
TEST(5, c.String(10) == sc); 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); t := c.Mul(c);
for base := uint(2); base <= 16; base++ { 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 { func Sum(n uint, scale *Big.Natural) *Big.Natural {
s := Big.Nat(0); s := nat_zero;
for ; n > 0; n-- { for ; n > 0; n-- {
s = s.Add(Big.Nat(uint64(n)).Mul(scale)); s = Add(s, Big.Nat(n).Mul(scale));
} }
return s; return s;
} }
func TestAdd() { func NatAdd() {
test_msg = "TestAddA"; test_msg = "NatAddA";
NAT_EQ(0, Add(nat_zero, nat_zero), nat_zero);
NAT_EQ(1, Add(nat_zero, c), c);
test_msg = "NatAddB";
for i := uint(0); i < 100; i++ {
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 = "TestAddB"; test_msg = "NatSubB";
for i := uint(0); i < 100; i++ { for i := uint(0); i < 100; i++ {
t := Big.Nat(uint64(i)); t := Sum(i, c);
TEST_EQ(i, Sum(i, c), t.Mul(t).Add(t).Shr(1).Mul(c)); for j := uint(0); j <= i; j++ {
t = t.Sub(Mul(Big.Nat(j), c));
}
NAT_EQ(i, t, nat_zero);
} }
} }
func TestShift() { func NatMul() {
test_msg = "TestShift1L"; 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 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(0, b.Shl(0).Cmp(b) == 0);
TEST(1, c.Shl(1).Cmp(c) > 0); TEST(1, c.Shl(1).Cmp(c) > 0);
test_msg = "TestShift1R"; test_msg = "NatShift1R";
TEST(0, b.Shr(0).Cmp(b) == 0); TEST(0, b.Shr(0).Cmp(b) == 0);
TEST(1, c.Shr(1).Cmp(c) < 0); TEST(1, c.Shr(1).Cmp(c) < 0);
test_msg = "TestShift2"; test_msg = "NatShift2";
for i := uint(0); i < 100; i++ { for i := uint(0); i < 100; i++ {
TEST(i, c.Shl(i).Shr(i).Cmp(c) == 0); TEST(i, c.Shl(i).Shr(i).Cmp(c) == 0);
} }
test_msg = "TestShift3L"; test_msg = "NatShift3L";
{ const m = 3; { const m = 3;
p := b; p := b;
f := Big.Nat(1<<m); f := Big.Nat(1<<m);
for i := uint(0); i < 100; i++ { for i := uint(0); i < 100; i++ {
TEST_EQ(i, b.Shl(i*m), p); NAT_EQ(i, b.Shl(i*m), p);
p = p.Mul(f); p = Mul(p, f);
} }
} }
test_msg = "TestShift3R"; test_msg = "NatShift3R";
{ p := c; { p := c;
for i := uint(0); c.Cmp(Big.NatZero) == 0; i++ { for i := uint(0); !p.IsZero(); i++ {
TEST_EQ(i, c.Shr(i), p); NAT_EQ(i, c.Shr(i), p);
p = p.Shr(1); p = p.Shr(1);
} }
} }
} }
func TestMul() { func IntShift() {
test_msg = "TestMulA"; test_msg = "IntShift1L";
TEST_EQ(0, b.Mul(Big.MulRange(0, 100)), Big.Nat(0)); TEST(0, ip.Shl(0).Cmp(ip) == 0);
TEST_EQ(0, b.Mul(Big.MulRange(21, 100)), c); TEST(1, ip.Shl(1).Cmp(ip) > 0);
test_msg = "TestMulB"; test_msg = "IntShift1R";
const n = 100; TEST(0, ip.Shr(0).Cmp(ip) == 0);
p := b.Mul(c).Shl(n); TEST(1, ip.Shr(1).Cmp(ip) < 0);
for i := uint(0); i < n; i++ {
TEST_EQ(i, b.Shl(i).Mul(c.Shl(n-i)), p);
}
}
test_msg = "IntShift2";
for i := uint(0); i < 100; i++ {
TEST(i, ip.Shl(i).Shr(i).Cmp(ip) == 0);
}
func TestDiv() { test_msg = "IntShift3L";
test_msg = "TestDivA"; { const m = 3;
TEST_EQ(0, c.Div(Big.Nat(1)), c); p := ip;
TEST_EQ(1, c.Div(Big.Nat(100)), Big.Fact(99)); f := Big.Int(1<<m);
TEST_EQ(2, b.Div(c), Big.Nat(0)); for i := uint(0); i < 100; i++ {
TEST_EQ(4, Big.Nat(1).Shl(100).Div(Big.Nat(1).Shl(90)), Big.Nat(1).Shl(10)); INT_EQ(i, ip.Shl(i*m), p);
TEST_EQ(5, c.Div(b), Big.MulRange(21, 100)); p = p.Mul(f);
}
}
test_msg = "TestDivB"; test_msg = "IntShift3R";
const n = 100; { p := ip;
p := Big.Fact(n); for i := uint(0); p.IsPos(); i++ {
for i := uint(0); i < n; i++ { INT_EQ(i, ip.Shr(i), p);
TEST_EQ(i, p.Div(Big.MulRange(1, uint64(i))), Big.MulRange(uint64(i+1), n)); p = p.Shr(1);
} }
}
test_msg = "IntShift4R";
//INT_EQ(0, Big.Int(-43).Shr(1), Big.Int(-43 >> 1));
//INT_EQ(1, ip.Neg().Shr(10), ip.Neg().Div(Big.Int(1).Shl(10)));
} }
func TestMod() { func NatCmp() {
test_msg = "TestModA"; test_msg = "NatCmp";
for i := uint(0); ; i++ { TEST(0, a.Cmp(a) == 0);
d := Big.Nat(1).Shl(i); TEST(1, a.Cmp(b) < 0);
if d.Cmp(c) < 0 { TEST(2, b.Cmp(a) > 0);
TEST_EQ(i, c.Add(d).Mod(c), d); TEST(3, a.Cmp(c) < 0);
} else { d := c.Add(b);
TEST_EQ(i, c.Add(d).Div(c), Big.Nat(2)); TEST(4, c.Cmp(d) < 0);
TEST_EQ(i, c.Add(d).Mod(c), d.Sub(c)); TEST(5, d.Cmp(c) > 0);
break; }
}
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() { func NatGcd() {
test_msg = "TestGcdA"; test_msg = "NatGcdA";
f := Big.Nat(99991); 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() { func NatPow() {
test_msg = "TestPowA"; test_msg = "NatPowA";
TEST_EQ(0, Big.Nat(2).Pow(0), Big.Nat(1)); NAT_EQ(0, nat_two.Pow(0), nat_one);
test_msg = "TestPowB"; test_msg = "NatPowB";
for i := uint(0); i < 100; i++ { 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() { func NatPop() {
test_msg = "TestPopA"; test_msg = "NatPopA";
TEST(0, Big.Nat(0).Pop() == 0); TEST(0, nat_zero.Pop() == 0);
TEST(1, Big.Nat(1).Pop() == 1); TEST(1, nat_one.Pop() == 1);
TEST(2, Big.Nat(10).Pop() == 2); TEST(2, Big.Nat(10).Pop() == 2);
TEST(3, Big.Nat(30).Pop() == 4); TEST(3, Big.Nat(30).Pop() == 4);
TEST(4, Big.Nat(0x1248f).Shl(33).Pop() == 8); TEST(4, Big.Nat(0x1248f).Shl(33).Pop() == 8);
test_msg = "TestPopB"; test_msg = "NatPopB";
for i := uint(0); i < 100; i++ { 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() { func main() {
TestLog2(); // Naturals
TestConv(); NatConv();
TestAdd(); NatAdd();
TestShift(); NatSub();
TestMul(); NatMul();
TestDiv(); NatDiv();
TestMod(); NatMod();
TestGcd(); NatShift();
TestPow(); NatCmp();
TestPop(); NatLog2();
NatGcd();
NatPow();
NatPop();
// Integers
// TODO add more tests
IntConv();
IntQuoRem();
IntDivMod();
IntShift();
// Rationals
// TODO add more tests
RatConv();
print("PASSED\n"); print("PASSED\n");
} }
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment