| // 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 |
| // |
| // This package has been designed for ease of use but the functions it provides |
| // are likely to be quite slow. It may be deprecated eventually. Use package |
| // big instead, if possible. |
| // |
| package bignum |
| |
| import ( |
| "fmt" |
| ) |
| |
| // TODO(gri) Complete the set of in-place operations. |
| |
| // ---------------------------------------------------------------------------- |
| // 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 // word width |
| 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. Keep around. |
| /* |
| func dump(x Natural) { |
| 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 |
| |
| |
| // Nat creates a small natural number with value x. |
| // |
| func Nat(x uint64) Natural { |
| if x == 0 { |
| return nil // len == 0 |
| } |
| |
| // single-digit values |
| // (note: cannot re-use preallocated values because |
| // the in-place operations may overwrite them) |
| if x < _B { |
| return Natural{digit(x)} |
| } |
| |
| // compute number of digits required to represent x |
| // (this is usually 1 or 2, but the algorithm works |
| // for any base) |
| n := 0 |
| for t := x; t > 0; t >>= _W { |
| n++ |
| } |
| |
| // split x into digits |
| z := make(Natural, n) |
| for i := 0; i < n; i++ { |
| z[i] = digit(x & _M) |
| x >>= _W |
| } |
| |
| return z |
| } |
| |
| |
| // Value returns the lowest 64bits of x. |
| // |
| func (x Natural) Value() uint64 { |
| // single-digit values |
| n := len(x) |
| switch n { |
| case 0: |
| return 0 |
| case 1: |
| return uint64(x[0]) |
| } |
| |
| // multi-digit values |
| // (this is usually 1 or 2, but the algorithm works |
| // for any base) |
| z := uint64(0) |
| s := uint(0) |
| for i := 0; i < n && s < 64; i++ { |
| z += uint64(x[i]) << s |
| s += _W |
| } |
| |
| return z |
| } |
| |
| |
| // Predicates |
| |
| // 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-- |
| } |
| return x[0:n] |
| } |
| |
| |
| // nalloc returns a Natural of n digits. If z is large |
| // enough, z is resized and returned. Otherwise, a new |
| // Natural is allocated. |
| // |
| func nalloc(z Natural, n int) Natural { |
| size := n |
| if size <= 0 { |
| size = 4 |
| } |
| if size <= cap(z) { |
| return z[0:n] |
| } |
| return make(Natural, n, size) |
| } |
| |
| |
| // Nadd sets *zp to the sum x + y. |
| // *zp may be x or y. |
| // |
| func Nadd(zp *Natural, x, y Natural) { |
| n := len(x) |
| m := len(y) |
| if n < m { |
| Nadd(zp, y, x) |
| return |
| } |
| |
| z := nalloc(*zp, n+1) |
| c := digit(0) |
| 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++ |
| } |
| *zp = z[0:i] |
| } |
| |
| |
| // Add returns the sum z = x + y. |
| // |
| func (x Natural) Add(y Natural) Natural { |
| var z Natural |
| Nadd(&z, x, y) |
| return z |
| } |
| |
| |
| // Nsub sets *zp to the difference x - y for x >= y. |
| // If x < y, an underflow run-time error occurs (use Cmp to test if x >= y). |
| // *zp may be x or y. |
| // |
| func Nsub(zp *Natural, x, y Natural) { |
| n := len(x) |
| m := len(y) |
| if n < m { |
| panic("underflow") |
| } |
| |
| z := nalloc(*zp, n) |
| c := digit(0) |
| 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++ |
| } |
| if int64(c) < 0 { |
| panic("underflow") |
| } |
| *zp = normalize(z) |
| } |
| |
| |
| // 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 { |
| var z Natural |
| Nsub(&z, x, y) |
| return z |
| } |
| |
| |
| // Returns z1 = (x*y + c) div B, z0 = (x*y + c) mod B. |
| // |
| func muladd11(x, y, c digit) (digit, digit) { |
| z1, z0 := MulAdd128(uint64(x), uint64(y), uint64(c)) |
| return digit(z1<<(64-logB) | z0>>logB), digit(z0 & _M) |
| } |
| |
| |
| func mul1(z, x Natural, y digit) (c digit) { |
| for i := 0; i < len(x); i++ { |
| c, z[i] = muladd11(x[i], y, c) |
| } |
| return |
| } |
| |
| |
| // Nscale sets *z to the scaled value (*z) * d. |
| // |
| func Nscale(z *Natural, d uint64) { |
| switch { |
| case d == 0: |
| *z = Nat(0) |
| return |
| case d == 1: |
| return |
| case d >= _B: |
| *z = z.Mul1(d) |
| return |
| } |
| |
| c := mul1(*z, *z, digit(d)) |
| |
| if c != 0 { |
| n := len(*z) |
| if n >= cap(*z) { |
| zz := make(Natural, n+1) |
| for i, d := range *z { |
| zz[i] = d |
| } |
| *z = zz |
| } else { |
| *z = (*z)[0 : n+1] |
| } |
| (*z)[n] = c |
| } |
| } |
| |
| |
| // 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) |
| } |
| |
| |
| // Mul1 returns the product x * d. |
| // |
| func (x Natural) Mul1(d uint64) Natural { |
| switch { |
| case d == 0: |
| return Nat(0) |
| case d == 1: |
| return x |
| case isSmall(digit(d)): |
| muladd1(x, digit(d), 0) |
| case d >= _B: |
| return x.Mul(Nat(d)) |
| } |
| |
| z := make(Natural, len(x)+1) |
| c := mul1(z, x, digit(d)) |
| z[len(x)] = c |
| return normalize(z) |
| } |
| |
| |
| // Mul returns the product x * y. |
| // |
| func (x Natural) Mul(y Natural) Natural { |
| n := len(x) |
| m := len(y) |
| if n < m { |
| return y.Mul(x) |
| } |
| |
| if m == 0 { |
| return Nat(0) |
| } |
| |
| if m == 1 && y[0] < _B { |
| return x.Mul1(uint64(y[0])) |
| } |
| |
| 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++ { |
| c, z[i+j] = muladd11(x[i], d, z[i+j]+c) |
| } |
| 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 mul21(z, x []digit2, y digit2) digit2 { |
| c := digit(0) |
| f := digit(y) |
| for i := 0; i < len(x); i++ { |
| t := c + digit(x[i])*f |
| c, z[i] = t>>_W2, digit2(t&_M2) |
| } |
| return digit2(c) |
| } |
| |
| |
| func div21(z, x []digit2, y digit2) digit2 { |
| c := digit(0) |
| d := digit(y) |
| for i := len(x) - 1; i >= 0; i-- { |
| t := c<<_W2 + 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] = div21(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 { |
| mul21(x, x, digit2(f)) |
| mul21(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]) |
| 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! |
| } |
| x[k] = digit2((c + digit(x[k])) & _M2) |
| |
| // correct if trial digit was too large |
| if 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) |
| } |
| x[k] = digit2((c + digit(x[k])) & _M2) |
| assert(x[k] == 0) |
| // correct trial digit |
| q-- |
| } |
| |
| x[k] = digit2(q) |
| } |
| |
| // undo normalization for remainder |
| if f != 1 { |
| c := div21(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, _ := 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 { |
| _, 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 Natural, 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 Natural, 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 2's-complement 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 Natural) { |
| for i, e := range x { |
| z[i] = e |
| } |
| } |
| |
| |
| // AndNot returns the ``bitwise clear'' x &^ y for the 2's-complement representation of x and y. |
| // |
| func (x Natural) AndNot(y Natural) Natural { |
| n := len(x) |
| m := len(y) |
| if n < m { |
| m = n |
| } |
| |
| 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) |
| } |
| |
| |
| // Or returns the ``bitwise or'' x | y for the 2's-complement 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 2's-complement 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 |
| } |
| |
| |
| // 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 log2(x uint64) 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(uint64(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(uint64(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.State, 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 |
| } |
| |
| |
| // 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. The syntax of natural numbers follows the syntax |
| // of unsigned integer literals in Go. |
| // |
| // 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(uint64(a)) |
| case a+1 == b: |
| return Nat(uint64(a)).Mul(Nat(uint64(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 |
| } |