diff options
Diffstat (limited to 'src/pkg/bignum/bignum.go')
| -rwxr-xr-x | src/pkg/bignum/bignum.go | 1464 |
1 files changed, 1464 insertions, 0 deletions
diff --git a/src/pkg/bignum/bignum.go b/src/pkg/bignum/bignum.go new file mode 100755 index 000000000..b9ea66587 --- /dev/null +++ b/src/pkg/bignum/bignum.go @@ -0,0 +1,1464 @@ +// Copyright 2009 The Go Authors. All rights reserved. +// Use of this source code is governed by a BSD-style +// license that can be found in the LICENSE file. + +// A package for arbitrary precision arithmethic. +// It implements the following numeric types: +// +// - Natural unsigned integers +// - Integer signed integers +// - Rational rational numbers +// +package bignum + +import "fmt" + + +// ---------------------------------------------------------------------------- +// Internal representation +// +// A natural number of the form +// +// x = x[n-1]*B^(n-1) + x[n-2]*B^(n-2) + ... + x[1]*B + x[0] +// +// with 0 <= x[i] < B and 0 <= i < n is stored in a slice of length n, +// with the digits x[i] as the slice elements. +// +// A natural number is normalized if the slice contains no leading 0 digits. +// During arithmetic operations, denormalized values may occur but are +// always normalized before returning the final result. The normalized +// representation of 0 is the empty slice (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: +// +// 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; + _LogH = 4; // bits for a hex digit (= small number) + _LogB = _LogW - _LogH; // largest bit-width available + + // 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 functions + +func assert(p bool) { + if !p { + panic("assert failed"); + } +} + + +func isSmall(x digit) bool { + return x < 1<<_LogH; +} + + +// For debugging. +func dump(x []digit) { + print("[", len(x), "]"); + for i := len(x) - 1; i >= 0; i-- { + print(" ", x[i]); + } + println(); +} + + +// ---------------------------------------------------------------------------- +// Natural numbers + +// Natural represents an unsigned integer value of arbitrary precision. +// +type Natural []digit; + +var ( + natZero Natural = Natural{}; + natOne Natural = Natural{1}; + natTwo Natural = Natural{2}; + natTen Natural = Natural{10}; +) + + +// Nat creates a small natural number with value x. +// Implementation restriction: At the moment, only values +// x < (1<<60) are supported. +// +func Nat(x uint) Natural { + switch x { + case 0: return natZero; + case 1: return natOne; + case 2: return natTwo; + case 10: return natTen; + } + assert(digit(x) < _B); + return Natural{digit(x)}; +} + + +// IsEven returns true iff x is divisible by 2. +// +func (x Natural) IsEven() bool { + return len(x) == 0 || x[0]&1 == 0; +} + + +// IsOdd returns true iff x is not divisible by 2. +// +func (x Natural) IsOdd() bool { + return len(x) > 0 && x[0]&1 != 0; +} + + +// IsZero returns true iff x == 0. +// +func (x Natural) IsZero() bool { + return len(x) == 0; +} + + +// Operations +// +// Naming conventions +// +// c carry +// x, y operands +// z result +// n, m len(x), len(y) + +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; +} + + +// Add returns the sum x + y. +// +func (x Natural) Add(y Natural) Natural { + n := len(x); + m := len(y); + if n < m { + return y.Add(x); + } + + c := digit(0); + z := make(Natural, n + 1); + i := 0; + for i < m { + t := c + x[i] + y[i]; + c, z[i] = t>>_W, t&_M; + i++; + } + for i < n { + t := c + x[i]; + c, z[i] = t>>_W, t&_M; + i++; + } + if c != 0 { + z[i] = c; + i++; + } + + return z[0 : i]; +} + + +// Sub returns the difference x - y for x >= y. +// If x < y, an underflow run-time error occurs (use Cmp to test if x >= y). +// +func (x Natural) Sub(y Natural) Natural { + n := len(x); + m := len(y); + if n < m { + panic("underflow") + } + + c := digit(0); + z := make(Natural, n); + i := 0; + for i < m { + t := c + x[i] - y[i]; + c, z[i] = digit(int64(t)>>_W), t&_M; // requires arithmetic shift! + i++; + } + for i < n { + t := c + x[i]; + c, z[i] = digit(int64(t)>>_W), t&_M; // requires arithmetic shift! + i++; + } + for i > 0 && z[i - 1] == 0 { // normalize + i--; + } + + return z[0 : i]; +} + + +// 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, + // multiply the digits separately while avoiding overflow, + // and return the product as two separate digits. + + // This code also works for non-even bit widths W + // which is why there are separate constants below + // for half-digits. + const W2 = (_W + 1)/2; + const DW = W2*2 - _W; // 0 or 1 + const B2 = 1<<W2; + const M2 = _B2 - 1; + + // split x and y into sub-digits + // x = (x1*B2 + x0) + // y = (y1*B2 + y0) + x1, x0 := x>>W2, x&M2; + y1, y0 := y>>W2, y&M2; + + // 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<<W2 + t0)&_M; + z1 := t2<<DW + (t1 + t0>>W2)>>(_W-W2); + + return z1, z0; +} + + +// Mul returns the product x * y. +// +func (x Natural) Mul(y Natural) Natural { + n := len(x); + m := len(y); + + z := make(Natural, n + m); + for j := 0; j < m; j++ { + d := y[j]; + if d != 0 { + c := digit(0); + for i := 0; i < n; i++ { + // z[i+j] += c + x[i]*d; + z1, z0 := mul11(x[i], d); + t := c + z[i+j] + z0; + c, z[i+j] = t>>_W, t&_M; + c += z1; + } + z[n+j] = c; + } + } + + return normalize(z); +} + + +// DivMod needs multi-precision division, which is not available if digit +// 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 { + n := len(x); + z := make([]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 >> _W2 & _M2); + } + + // normalize result + k := 2*n; + for k > 0 && z[k - 1] == 0 { k-- } + return z[0 : k]; // trim leading 0's +} + + +func pack(x []digit2) Natural { + n := (len(x) + 1) / 2; + z := make(Natural, n); + if len(x) & 1 == 1 { + // handle odd len(x) + n--; + z[n] = digit(x[n*2]); + } + for i := 0; i < n; i++ { + z[i] = digit(x[i*2 + 1]) << _W2 | digit(x[i*2]); + } + return normalize(z); +} + + +func mul1(z, x []digit2, y digit2) digit2 { + n := len(x); + c := digit(0); + 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); + c := digit(0); + 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., 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), +// 579-601. John Wiley & Sons, Ltd. + +func divmod(x, y []digit2) ([]digit2, []digit2) { + n := len(x); + m := len(y); + 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 + // result is shifted left by 1 in place! + x[0] = div1(x[1 : n+1], x[0 : n], y[0]); + + } else if m > n { + // 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); + + // normalize x and y + // 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)<<_W2 + digit(y2); + for i := n-m; i >= 0; i-- { + k := i+m; + + // compute trial digit (Knuth) + var q digit; + { x0, x1, x2 := digit(x[k]), digit(x[k-1]), digit(x[k-2]); + if x0 != y1 { + q = (x0<<_W2 + x1)/y1; + } else { + q = _B2 - 1; + } + for y2*q > (x0<<_W2 + x1 - y1*q)<<_W2 + x2 { + q-- + } + } + + // subtract y*q + c := digit(0); + for j := 0; j < m; j++ { + 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 + if c + digit(x[k]) != 0 { + // add y + c := digit(0); + for j := 0; j < m; j++ { + t := c + digit(x[i+j]) + digit(y[j]); + c, x[i+j] = t >> _W2, digit2(t & _M2) + } + assert(c + digit(x[k]) == 0); + // correct trial digit + q--; + } + + x[k] = digit2(q); + } + + // undo normalization for remainder + if f != 1 { + c := div1(x[0 : m], x[0 : m], digit2(f)); + assert(c == 0); + } + } + + return x[m : n+1], x[0 : m]; +} + + +// Div returns the quotient q = x / y for y > 0, +// with x = y*q + r and 0 <= r < y. +// If y == 0, a division-by-zero run-time error occurs. +// +func (x Natural) Div(y Natural) Natural { + q, r := divmod(unpack(x), unpack(y)); + return pack(q); +} + + +// Mod returns the modulus r of the division x / y for y > 0, +// with x = y*q + r and 0 <= r < y. +// If y == 0, a division-by-zero run-time error occurs. +// +func (x Natural) Mod(y Natural) Natural { + q, r := divmod(unpack(x), unpack(y)); + return pack(r); +} + + +// DivMod returns the pair (x.Div(y), x.Mod(y)) for y > 0. +// If y == 0, a division-by-zero run-time error occurs. +// +func (x Natural) DivMod(y Natural) (Natural, Natural) { + q, r := divmod(unpack(x), unpack(y)); + return pack(q), pack(r); +} + + +func shl(z, x []digit, s uint) digit { + assert(s <= _W); + n := len(x); + c := digit(0); + for i := 0; i < n; i++ { + c, z[i] = x[i] >> (_W-s), x[i] << s & _M | c; + } + return c; +} + + +// Shl implements ``shift left'' x << s. It returns x * 2^s. +// +func (x Natural) Shl(s uint) Natural { + n := uint(len(x)); + m := n + s/_W; + z := make(Natural, m+1); + + z[m] = shl(z[m-n : m], x, s%_W); + + return normalize(z); +} + + +func shr(z, x []digit, s uint) digit { + assert(s <= _W); + n := len(x); + c := digit(0); + for i := n - 1; i >= 0; i-- { + c, z[i] = x[i] << (_W-s) & _M, x[i] >> s | c; + } + return c; +} + + +// Shr implements ``shift right'' x >> s. It returns x / 2^s. +// +func (x Natural) Shr(s uint) Natural { + n := uint(len(x)); + m := n - s/_W; + if m > n { // check for underflow + m = 0; + } + z := make(Natural, m); + + shr(z, x[n-m : n], s%_W); + + return normalize(z); +} + + +// And returns the ``bitwise and'' x & y for the binary representation of x and y. +// +func (x Natural) And(y Natural) Natural { + n := len(x); + m := len(y); + if n < m { + return y.And(x); + } + + z := make(Natural, m); + for i := 0; i < m; i++ { + z[i] = x[i] & y[i]; + } + // upper bits are 0 + + return normalize(z); +} + + +func copy(z, x []digit) { + for i, e := range x { + z[i] = e + } +} + + +// Or returns the ``bitwise or'' x | y for the binary representation of x and y. +// +func (x Natural) Or(y Natural) Natural { + n := len(x); + m := len(y); + if n < m { + return y.Or(x); + } + + z := make(Natural, n); + for i := 0; i < m; i++ { + z[i] = x[i] | y[i]; + } + copy(z[m : n], x[m : n]); + + return z; +} + + +// Xor returns the ``bitwise exclusive or'' x ^ y for the binary representation of x and y. +// +func (x Natural) Xor(y Natural) Natural { + n := len(x); + m := len(y); + if n < m { + return y.Xor(x); + } + + z := make(Natural, n); + for i := 0; i < m; i++ { + z[i] = x[i] ^ y[i]; + } + copy(z[m : n], x[m : n]); + + return normalize(z); +} + + +// Cmp compares x and y. The result is an int value +// +// < 0 if x < y +// == 0 if x == y +// > 0 if x > y +// +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; +} + + +// Log2 computes the binary logarithm of x for x > 0. +// The result is the integer n for which 2^n <= x < 2^(n+1). +// If x == 0 a run-time error occurs. +// +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 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 + x[i]; + c, x[i] = t%d, t/d; + } + + return normalize(x), c; +} + + +// ToString converts x to a string for a given base, with 2 <= base <= 16. +// +func (x Natural) ToString(base uint) string { + if len(x) == 0 { + return "0"; + } + + // allocate buffer for conversion + assert(2 <= base && base <= 16); + n := (x.Log2() + 1) / log2(digit(base)) + 1; // +1: round up + s := make([]byte, n); + + // don't destroy x + t := make(Natural, len(x)); + copy(t, x); + + // convert + i := n; + for !t.IsZero() { + i--; + var d digit; + t, d = divmod1(t, digit(base)); + s[i] = "0123456789abcdef"[d]; + }; + + return string(s[i : n]); +} + + +// String converts x to its decimal string representation. +// x.String() is the same as x.ToString(10). +// +func (x Natural) String() string { + return x.ToString(10); +} + + +func fmtbase(c int) uint { + switch c { + case 'b': return 2; + case 'o': return 8; + case 'x': return 16; + } + return 10; +} + + +// Format is a support routine for fmt.Formatter. It accepts +// the formats 'b' (binary), 'o' (octal), and 'x' (hexadecimal). +// +func (x Natural) Format(h fmt.Formatter, c int) { + fmt.Fprintf(h, "%s", x.ToString(fmtbase(c))); +} + + +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 := make(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); +} + + +// NatFromString returns the natural number corresponding to the +// longest possible prefix of s representing a natural number in a +// given conversion base, the actual conversion base used, and the +// prefix length. +// +// If the base argument is 0, the string prefix determines the actual +// conversion base. A prefix of ``0x'' or ``0X'' selects base 16; the +// ``0'' prefix selects base 8. Otherwise the selected base is 10. +// +func NatFromString(s string, base uint) (Natural, uint, int) { + // 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; + } + } + + return x, base, i; +} + + +// Natural number functions + +func pop1(x digit) uint { + n := uint(0); + for x != 0 { + x &= x-1; + n++; + } + return n; +} + + +// Pop computes the ``population count'' of (the number of 1 bits in) x. +// +func (x Natural) Pop() uint { + n := uint(0); + for i := len(x) - 1; i >= 0; i-- { + n += pop1(x[i]); + } + return n; +} + + +// Pow computes x to the power of n. +// +func (xp Natural) Pow(n uint) Natural { + z := Nat(1); + x := xp; + 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; +} + + +// MulRange computes the product of all the unsigned integers +// in the range [a, b] inclusively. +// +func MulRange(a, b uint) Natural { + switch { + case a > b: return Nat(1); + case a == b: return Nat(a); + case a + 1 == b: return Nat(a).Mul(Nat(b)); + } + m := (a + b)>>1; + assert(a <= m && m < b); + return MulRange(a, m).Mul(MulRange(m + 1, b)); +} + + +// Fact computes the factorial of n (Fact(n) == MulRange(2, n)). +// +func Fact(n uint) Natural { + // Using MulRange() instead of the basic for-loop + // lead to faster factorial computation. + return MulRange(2, n); +} + + +// Binomial computes the binomial coefficient of (n, k). +// +func Binomial(n, k uint) Natural { + return MulRange(n-k+1, n).Div(MulRange(1, k)); +} + + +// Gcd computes the gcd of x and y. +// +func (x Natural) Gcd(y Natural) Natural { + // Euclidean algorithm. + a, b := x, y; + for !b.IsZero() { + a, b = b, a.Mod(b); + } + return a; +} + + +// ---------------------------------------------------------------------------- +// Integer numbers +// +// Integers are normalized if the mantissa is normalized and the sign is +// false for mant == 0. Use MakeInt to create normalized Integers. + +// Integer represents a signed integer value of arbitrary precision. +// +type Integer struct { + sign bool; + mant Natural; +} + + +// MakeInt makes an integer given a sign and a mantissa. +// The number is positive (>= 0) if sign is false or the +// mantissa is zero; it is negative otherwise. +// +func MakeInt(sign bool, mant Natural) *Integer { + if mant.IsZero() { + sign = false; // normalize + } + return &Integer{sign, mant}; +} + + +// Int creates a small integer with value x. +// Implementation restriction: At the moment, only values +// with an absolute value |x| < (1<<60) are supported. +// +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 { + ux = uint(-x); + } + } else { + ux = uint(x); + } + return MakeInt(sign, Nat(ux)); +} + + +// Predicates + +// IsEven returns true iff x is divisible by 2. +// +func (x *Integer) IsEven() bool { + return x.mant.IsEven(); +} + + +// IsOdd returns true iff x is not divisible by 2. +// +func (x *Integer) IsOdd() bool { + return x.mant.IsOdd(); +} + + +// IsZero returns true iff x == 0. +// +func (x *Integer) IsZero() bool { + return x.mant.IsZero(); +} + + +// IsNeg returns true iff x < 0. +// +func (x *Integer) IsNeg() bool { + return x.sign && !x.mant.IsZero() +} + + +// IsPos returns true iff x >= 0. +// +func (x *Integer) IsPos() bool { + return !x.sign && !x.mant.IsZero() +} + + +// Operations + +// Neg returns the negated value of x. +// +func (x *Integer) Neg() *Integer { + return MakeInt(!x.sign, x.mant); +} + + +// Add returns the sum x + y. +// +func (x *Integer) Add(y *Integer) *Integer { + var z *Integer; + if x.sign == y.sign { + // x + y == x + y + // (-x) + (-y) == -(x + y) + 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 = MakeInt(false, x.mant.Sub(y.mant)); + } else { + z = MakeInt(true, y.mant.Sub(x.mant)); + } + } + if x.sign { + z.sign = !z.sign; + } + return z; +} + + +// Sub returns the difference x - y. +// +func (x *Integer) Sub(y *Integer) *Integer { + var z *Integer; + if x.sign != y.sign { + // x - (-y) == x + y + // (-x) - y == -(x + y) + 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 = MakeInt(false, x.mant.Sub(y.mant)); + } else { + z = MakeInt(true, y.mant.Sub(x.mant)); + } + } + if x.sign { + z.sign = !z.sign; + } + return z; +} + + +// Mul returns the product x * y. +// +func (x *Integer) Mul(y *Integer) *Integer { + // x * y == x * y + // x * (-y) == -(x * y) + // (-x) * y == -(x * y) + // (-x) * (-y) == x * y + return MakeInt(x.sign != y.sign, x.mant.Mul(y.mant)); +} + + +// MulNat returns the product x * y, where y is a (unsigned) natural number. +// +func (x *Integer) MulNat(y Natural) *Integer { + // x * y == x * y + // (-x) * y == -(x * y) + return MakeInt(x.sign, x.mant.Mul(y)); +} + + +// Quo returns the quotient q = x / y for y != 0. +// If y == 0, a division-by-zero run-time error occurs. +// +// 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 MakeInt(x.sign != y.sign, x.mant.Div(y.mant)); +} + + +// Rem returns the remainder r of the division x / y for y != 0, +// with r = x - y*x.Quo(y). Unless r is zero, its sign corresponds +// to the sign of x. +// If y == 0, a division-by-zero run-time error occurs. +// +func (x *Integer) Rem(y *Integer) *Integer { + // x % y == x % y + // x % (-y) == x % y + // (-x) % y == -(x % y) + // (-x) % (-y) == -(x % y) + return MakeInt(x.sign, x.mant.Mod(y.mant)); +} + + +// QuoRem returns the pair (x.Quo(y), x.Rem(y)) for y != 0. +// If y == 0, a division-by-zero run-time error occurs. +// +func (x *Integer) QuoRem(y *Integer) (*Integer, *Integer) { + q, r := x.mant.DivMod(y.mant); + return MakeInt(x.sign != y.sign, q), MakeInt(x.sign, r); +} + + +// Div returns the quotient q = x / y for y != 0. +// If y == 0, a division-by-zero run-time error occurs. +// +// Div and Mod implement Euclidian division and modulus: +// +// q = x.Div(y) +// r = x.Mod(y) with: 0 <= r < |q| and: y = x*q + r +// +// (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.QuoRem(y); + if r.IsNeg() { + if y.IsPos() { + q = q.Sub(Int(1)); + } else { + q = q.Add(Int(1)); + } + } + return q; +} + + +// Mod returns the modulus r of the division x / y for y != 0, +// with r = x - y*x.Div(y). r is always positive. +// If y == 0, a division-by-zero run-time error occurs. +// +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; +} + + +// DivMod returns the pair (x.Div(y), x.Mod(y)). +// +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; +} + + +// Shl implements ``shift left'' x << s. It returns x * 2^s. +// +func (x *Integer) Shl(s uint) *Integer { + return MakeInt(x.sign, x.mant.Shl(s)); +} + + +// Shr implements ``shift right'' x >> s. It returns x / 2^s. +// Implementation restriction: Shl is not yet implemented for negative x. +// +func (x *Integer) Shr(s uint) *Integer { + z := MakeInt(x.sign, x.mant.Shr(s)); + if x.IsNeg() { + panic("UNIMPLEMENTED Integer.Shr of negative values"); + } + return z; +} + + +// And returns the ``bitwise and'' x & y for the binary representation of x and y. +// Implementation restriction: And is not implemented for negative x. +// +func (x *Integer) And(y *Integer) *Integer { + var z *Integer; + if !x.sign && !y.sign { + z = MakeInt(false, x.mant.And(y.mant)); + } else { + panic("UNIMPLEMENTED Integer.And of negative values"); + } + return z; +} + + +// Or returns the ``bitwise or'' x | y for the binary representation of x and y. +// Implementation restriction: Or is not implemented for negative x. +// +func (x *Integer) Or(y *Integer) *Integer { + var z *Integer; + if !x.sign && !y.sign { + z = MakeInt(false, x.mant.Or(y.mant)); + } else { + panic("UNIMPLEMENTED Integer.Or of negative values"); + } + return z; +} + + +// Xor returns the ``bitwise xor'' x | y for the binary representation of x and y. +// Implementation restriction: Xor is not implemented for negative integers. +// +func (x *Integer) Xor(y *Integer) *Integer { + var z *Integer; + if !x.sign && !y.sign { + z = MakeInt(false, x.mant.Xor(y.mant)); + } else { + panic("UNIMPLEMENTED Integer.Xor of negative values"); + } + return z; +} + + +// Cmp compares x and y. The result is an int value +// +// < 0 if x < y +// == 0 if x == y +// > 0 if x > y +// +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; +} + + +// ToString converts x to a string for a given base, with 2 <= base <= 16. +// +func (x *Integer) ToString(base uint) string { + if x.mant.IsZero() { + return "0"; + } + var s string; + if x.sign { + s = "-"; + } + return s + x.mant.ToString(base); +} + + +// String converts x to its decimal string representation. +// x.String() is the same as x.ToString(10). +// +func (x *Integer) String() string { + return x.ToString(10); +} + + +// Format is a support routine for fmt.Formatter. It accepts +// the formats 'b' (binary), 'o' (octal), and 'x' (hexadecimal). +// +func (x *Integer) Format(h fmt.Formatter, c int) { + fmt.Fprintf(h, "%s", x.ToString(fmtbase(c))); +} + + +// IntFromString returns the integer corresponding to the +// longest possible prefix of s representing an integer in a +// given conversion base, the actual conversion base used, and +// the prefix length. +// +// If the base argument is 0, the string prefix determines the actual +// conversion base. A prefix of ``0x'' or ``0X'' selects base 16; the +// ``0'' prefix selects base 8. Otherwise the selected base is 10. +// +func IntFromString(s string, base uint) (*Integer, uint, int) { + // skip sign, if any + i0 := 0; + if len(s) > 0 && (s[0] == '-' || s[0] == '+') { + i0 = 1; + } + + mant, base, slen := NatFromString(s[i0 : len(s)], base); + + return MakeInt(i0 > 0 && s[0] == '-', mant), base, i0 + slen; +} + + +// ---------------------------------------------------------------------------- +// Rational numbers + +// Rational represents a quotient a/b of arbitrary precision. +// +type Rational struct { + a *Integer; // numerator + b Natural; // denominator +} + + +// MakeRat makes a rational number given a numerator a and a denominator b. +// +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}; +} + + +// Rat creates a small rational number with value a0/b0. +// Implementation restriction: At the moment, only values a0, b0 +// with an absolute value |a0|, |b0| < (1<<60) are supported. +// +func Rat(a0 int, b0 int) *Rational { + a, b := Int(a0), Int(b0); + if b.sign { + a = a.Neg(); + } + return MakeRat(a, b.mant); +} + + +// Predicates + +// IsZero returns true iff x == 0. +// +func (x *Rational) IsZero() bool { + return x.a.IsZero(); +} + + +// IsNeg returns true iff x < 0. +// +func (x *Rational) IsNeg() bool { + return x.a.IsNeg(); +} + + +// IsPos returns true iff x > 0. +// +func (x *Rational) IsPos() bool { + return x.a.IsPos(); +} + + +// IsInt returns true iff x can be written with a denominator 1 +// in the form x == x'/1; i.e., if x is an integer value. +// +func (x *Rational) IsInt() bool { + return x.b.Cmp(Nat(1)) == 0; +} + + +// Operations + +// Neg returns the negated value of x. +// +func (x *Rational) Neg() *Rational { + return MakeRat(x.a.Neg(), x.b); +} + + +// Add returns the sum x + y. +// +func (x *Rational) Add(y *Rational) *Rational { + return MakeRat((x.a.MulNat(y.b)).Add(y.a.MulNat(x.b)), x.b.Mul(y.b)); +} + + +// Sub returns the difference x - y. +// +func (x *Rational) Sub(y *Rational) *Rational { + return MakeRat((x.a.MulNat(y.b)).Sub(y.a.MulNat(x.b)), x.b.Mul(y.b)); +} + + +// Mul returns the product x * y. +// +func (x *Rational) Mul(y *Rational) *Rational { + return MakeRat(x.a.Mul(y.a), x.b.Mul(y.b)); +} + + +// Quo returns the quotient x / y for y != 0. +// If y == 0, a division-by-zero run-time error occurs. +// +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); +} + + +// Cmp compares x and y. The result is an int value +// +// < 0 if x < y +// == 0 if x == y +// > 0 if x > y +// +func (x *Rational) Cmp(y *Rational) int { + return (x.a.MulNat(y.b)).Cmp(y.a.MulNat(x.b)); +} + + +// ToString converts x to a string for a given base, with 2 <= base <= 16. +// The string representation is of the form "n" if x is an integer; otherwise +// it is of form "n/d". +// +func (x *Rational) ToString(base uint) string { + s := x.a.ToString(base); + if !x.IsInt() { + s += "/" + x.b.ToString(base); + } + return s; +} + + +// String converts x to its decimal string representation. +// x.String() is the same as x.ToString(10). +// +func (x *Rational) String() string { + return x.ToString(10); +} + + +// Format is a support routine for fmt.Formatter. It accepts +// the formats 'b' (binary), 'o' (octal), and 'x' (hexadecimal). +// +func (x *Rational) Format(h fmt.Formatter, c int) { + fmt.Fprintf(h, "%s", x.ToString(fmtbase(c))); +} + + +// RatFromString returns the rational number corresponding to the +// longest possible prefix of s representing a rational number in a +// given conversion base, the actual conversion base used, and the +// prefix length. +// +// If the base argument is 0, the string prefix determines the actual +// conversion base. A prefix of ``0x'' or ``0X'' selects base 16; the +// ``0'' prefix selects base 8. Otherwise the selected base is 10. +// +func RatFromString(s string, base uint) (*Rational, uint, int) { + // read nominator + a, abase, alen := IntFromString(s, base); + b := Nat(1); + + // read denominator or fraction, if any + var blen int; + if alen < len(s) { + ch := s[alen]; + if ch == '/' { + alen++; + b, base, blen = NatFromString(s[alen : len(s)], base); + } else if ch == '.' { + alen++; + b, base, blen = NatFromString(s[alen : len(s)], abase); + assert(base == abase); + f := Nat(base).Pow(uint(blen)); + a = MakeInt(a.sign, a.mant.Mul(f).Add(b)); + b = f; + } + } + + return MakeRat(a, b), base, alen + blen; +} |