ref: 3400b0428b9147012553c6f2280147d8dc99bbf1
dir: /lib/std/bigint.myr/
use "alloc.use" use "chartype.use" use "cmp.use" use "die.use" use "extremum.use" use "hasprefix.use" use "option.use" use "slcp.use" use "sldup.use" use "slfill.use" use "slpush.use" use "types.use" use "utf.use" use "errno.use" pkg std = type bigint = struct dig : uint32[:] /* little endian, no leading zeros. */ sign : int /* -1 for -ve, 0 for zero, 1 for +ve. */ ;; /* administrivia */ generic mkbigint : (v : @a::(numeric,integral) -> bigint#) const bigfree : (a : bigint# -> void) const bigdup : (a : bigint# -> bigint#) const bigassign : (d : bigint#, s : bigint# -> bigint#) const bigmove : (d : bigint#, s : bigint# -> bigint#) const bigparse : (s : byte[:] -> option(bigint#)) const bigclear : (a : bigint# -> bigint#) const bigbfmt : (b : byte[:], a : bigint#, base : int -> size) /* const bigtoint : (a : bigint# -> @a::(numeric,integral)) */ /* some useful predicates */ const bigiszero : (a : bigint# -> bool) const bigeq : (a : bigint#, b : bigint# -> bool) generic bigeqi : (a : bigint#, b : @a::(numeric,integral) -> bool) const bigcmp : (a : bigint#, b : bigint# -> order) /* bigint*bigint -> bigint ops */ const bigadd : (a : bigint#, b : bigint# -> bigint#) const bigsub : (a : bigint#, b : bigint# -> bigint#) const bigmul : (a : bigint#, b : bigint# -> bigint#) const bigdiv : (a : bigint#, b : bigint# -> bigint#) const bigmod : (a : bigint#, b : bigint# -> bigint#) const bigdivmod : (a : bigint#, b : bigint# -> (bigint#, bigint#)) const bigshl : (a : bigint#, b : bigint# -> bigint#) const bigshr : (a : bigint#, b : bigint# -> bigint#) const bigmodpow : (b : bigint#, e : bigint#, m : bigint# -> bigint#) /* const bigpow : (a : bigint#, b : bigint# -> bigint#) */ /* bigint*int -> bigint ops */ generic bigaddi : (a : bigint#, b : @a::(integral,numeric) -> bigint#) generic bigsubi : (a : bigint#, b : @a::(integral,numeric) -> bigint#) generic bigmuli : (a : bigint#, b : @a::(integral,numeric) -> bigint#) generic bigdivi : (a : bigint#, b : @a::(integral,numeric) -> bigint#) generic bigshli : (a : bigint#, b : @a::(integral,numeric) -> bigint#) generic bigshri : (a : bigint#, b : @a::(integral,numeric) -> bigint#) /* const bigpowi : (a : bigint#, b : uint64 -> bigint#) */ ;; const Base = 0x100000000ul generic mkbigint = {v : @a::(integral,numeric) var a var val a = zalloc() if v < 0 a.sign = -1 v = -v elif v > 0 a.sign = 1 ;; val = v castto(uint64) a.dig = slpush([][:], val castto(uint32)) if val > Base a.dig = slpush(a.dig, (val/Base) castto(uint32)) ;; -> trim(a) } const bigfree = {a slfree(a.dig) free(a) } const bigdup = {a -> bigassign(zalloc(), a) } const bigassign = {d, s slfree(d.dig) d# = s# d.dig = sldup(s.dig) -> d } const bigmove = {d, s slfree(d.dig) d# = s# s.dig = [][:] s.sign = 0 -> d } const bigclear = {v std.slfree(v.dig) v.sign = 0 v.dig = [][:] -> v } /* for now, just dump out something for debugging... */ const bigbfmt = {buf, x, base const digitchars = [ '0','1','2','3','4','5','6','7','8','9', 'a','b','c','d','e','f', 'g', 'h', 'i', 'j', 'k', 'l', 'm', 'n', 'o', 'p', 'q', 'r', 's', 't', 'u', 'v', 'w', 'x', 'y', 'z'] var v, val var n, i var tmp, rem, b if base < 0 || base > 36 die("invalid base in bigbfmt\n") ;; if bigiszero(x) n ;; if base == 0 b = mkbigint(10) else b = mkbigint(base) ;; n = 0 val = bigdup(x) /* generate the digits in reverse order */ while !bigiszero(val) (v, rem) = bigdivmod(val, b) if rem.dig.len > 0 n += encode(buf[n:], digitchars[rem.dig[0]]) else n += encode(buf[n:], '0') ;; bigfree(val) bigfree(rem) val = v ;; bigfree(val) bigfree(b) /* this is done last, so we get things right when we reverse the string */ if x.sign == 0 n += encode(buf[n:], '0') elif x.sign == -1 n += encode(buf[n:], '-') ;; /* we only generated ascii digits, so this works for reversing. */ for i = 0; i < n/2; i++ tmp = buf[i] buf[i] = buf[n - i - 1] buf[n - i - 1] = tmp ;; -> n } const bigparse = {str var c, val : int, base var v, b var a if hasprefix(str, "0x") || hasprefix(str, "0X") base = 16 elif hasprefix(str, "0o") || hasprefix(str, "0O") base = 8 elif hasprefix(str, "0b") || hasprefix(str, "0B") base = 2 else base = 10 ;; if base != 10 str = str[2:] ;; a = mkbigint(0) b = mkbigint(base) /* efficiency hack: to save allocations, just mutate v[0]. The value will always fit in one digit. */ v = mkbigint(1) while str.len != 0 (c, str) = striter(str) if c == '_' continue ;; val = charval(c, base) if val < 0 || val > base bigfree(a) bigfree(b) bigfree(v) -> `None ;; v.dig[0] = val castto(uint32) if val == 0 v.sign = 0 else v.sign = 1 ;; bigmul(a, b) bigadd(a, v) ;; -> `Some a } const bigiszero = {v -> v.dig.len == 0 } const bigeq = {a, b if a.sign != b.sign || a.dig.len != b.dig.len -> false ;; for var i = 0; i < a.dig.len; i++ if a.dig[i] != b.dig[i] -> false ;; ;; -> true } generic bigeqi = {a, b var v var dig : uint32[2] bigdigit(&v, b < 0, b castto(uint64), dig[:]) -> bigeq(a, &v) } const bigcmp = {a, b var da, db, sa, sb sa = a.sign castto(int64) sb = b.sign castto(int64) if sa < sb -> `Before elif sa > sb -> `After elif a.dig.len < b.dig.len -> signedorder(-sa) elif a.dig.len > b.dig.len -> signedorder(sa) else /* otherwise, the one with the first larger digit is bigger */ for var i = a.dig.len; i > 0; i-- da = a.dig[i - 1] castto(int64) db = b.dig[i - 1] castto(int64) if da != db -> signedorder(sa * (da - db)) ;; ;; ;; -> `Equal } const signedorder = {sign if sign < 0 -> `Before elif sign == 0 -> `Equal else -> `After ;; } /* a += b */ const bigadd = {a, b if a.sign == b.sign || a.sign == 0 a.sign = b.sign -> uadd(a, b) elif b.sign == 0 -> a else match bigcmp(a, b) | `Before: /* a is negative */ a.sign = b.sign -> usub(b, a) | `After: /* b is negative */ -> usub(a, b) | `Equal: die("Impossible. Equal vals with different sign.") ;; ;; } /* adds two unsigned values together. */ const uadd = {a, b var v, i var carry var n carry = 0 n = max(a.dig.len, b.dig.len) /* guaranteed to carry no more than one value */ a.dig = slzgrow(a.dig, n + 1) for i = 0; i < n; i++ v = (a.dig[i] castto(uint64)) + carry; if i < b.dig.len v += (b.dig[i] castto(uint64)) ;; if v >= Base carry = 1 else carry = 0 ;; a.dig[i] = v castto(uint32) ;; a.dig[i] += carry castto(uint32) -> trim(a) } /* a -= b */ const bigsub = {a, b /* 0 - x = -x */ if a.sign == 0 bigassign(a, b) a.sign = -b.sign -> a /* x - 0 = x */ elif b.sign == 0 -> a elif a.sign != b.sign -> uadd(a, b) else match bigcmp(a, b) | `Before: /* a is negative */ a.sign = b.sign -> usub(b, a) | `After: /* b is negative */ -> usub(a, b) | `Equal: -> bigclear(a) ;; ;; -> a } /* subtracts two unsigned values, where 'a' is strictly greater than 'b' */ const usub = {a, b var carry var v, i carry = 0 for i = 0; i < a.dig.len; i++ v = (a.dig[i] castto(int64)) - carry if i < b.dig.len v -= (b.dig[i] castto(int64)) ;; if v < 0 carry = 1 else carry = 0 ;; a.dig[i] = v castto(uint32) ;; -> trim(a) } /* a *= b */ const bigmul = {a, b var i, j var ai, bj, wij var carry, t var w if a.sign == 0 || b.sign == 0 a.sign = 0 slfree(a.dig) a.dig = [][:] -> a elif a.sign != b.sign a.sign = -1 else a.sign = 1 ;; w = slzalloc(a.dig.len + b.dig.len) for j = 0; j < b.dig.len; j++ carry = 0 for i = 0; i < a.dig.len; i++ ai = a.dig[i] castto(uint64) bj = b.dig[j] castto(uint64) wij = w[i+j] castto(uint64) t = ai * bj + wij + carry w[i + j] = (t castto(uint32)) carry = t >> 32 ;; w[i+j] = carry castto(uint32) ;; slfree(a.dig) a.dig = w -> trim(a) } const bigdiv = {a : bigint#, b : bigint# -> bigint# var q, r (q, r) = bigdivmod(a, b) bigfree(r) -> bigmove(a, q) } const bigmod = {a : bigint#, b : bigint# -> bigint# var q, r (q, r) = bigdivmod(a, b) bigfree(q) -> bigmove(a, r) } /* a /= b */ const bigdivmod = {a : bigint#, b : bigint# -> (bigint#, bigint#) /* Implements bigint division using Algorithm D from Knuth: Seminumerical algorithms, Section 4.3.1. */ var m : int64, n : int64 var qhat, rhat, carry, shift var x, y, z, w, p, t /* temporaries */ var b0, aj var u, v var i, j : int64 var q if bigiszero(b) die("divide by zero\n") ;; /* if b > a, we trucate to 0, with remainder 'a' */ if a.dig.len < b.dig.len -> (mkbigint(0), bigdup(a)) ;; q = zalloc() q.dig = slzalloc(max(a.dig.len, b.dig.len) + 1) if a.sign != b.sign q.sign = -1 else q.sign = 1 ;; /* handle single digit divisor separately: the knuth algorithm needs at least 2 digits. */ if b.dig.len == 1 carry = 0 b0 = (b.dig[0] castto(uint64)) for j = a.dig.len; j > 0; j-- aj = (a.dig[j - 1] castto(uint64)) q.dig[j - 1] = (((carry << 32) + aj)/b0) castto(uint32) carry = (carry << 32) + aj - (q.dig[j-1] castto(uint64))*b0 ;; q = trim(q) -> (q, trim(mkbigint(carry castto(int32)))) ;; u = bigdup(a) v = bigdup(b) m = u.dig.len n = v.dig.len shift = nlz(v.dig[n - 1]) bigshli(u, shift) bigshli(v, shift) u.dig = slzgrow(u.dig, u.dig.len + 1) for j = m - n; j >= 0; j-- /* load a few temps */ x = u.dig[j + n] castto(uint64) y = u.dig[j + n - 1] castto(uint64) z = v.dig[n - 1] castto(uint64) w = v.dig[n - 2] castto(uint64) t = u.dig[j + n - 2] castto(uint64) /* estimate qhat */ qhat = (x*Base + y)/z rhat = (x*Base + y) - qhat*z :divagain if qhat >= Base || (qhat * w) > (rhat*Base + t) qhat-- rhat += z if rhat < Base goto divagain ;; ;; /* multiply and subtract */ carry = 0 for i = 0; i < n; i++ p = qhat * (v.dig[i] castto(uint64)) t = (u.dig[i+j] castto(uint64)) - carry - (p % Base) u.dig[i+j] = t castto(uint32) carry = (((p castto(int64)) >> 32) - ((t castto(int64)) >> 32)) castto(uint64); ;; t = (u.dig[j + n] castto(uint64)) - carry u.dig[j + n] = t castto(uint32) q.dig[j] = qhat castto(uint32) /* adjust */ if t castto(int64) < 0 q.dig[j]-- carry = 0 for i = 0; i < n; i++ t = (u.dig[i+j] castto(uint64)) + (v.dig[i] castto(uint64)) + carry u.dig[i+j] = t castto(uint32) carry = t >> 32 ;; u.dig[j+n] = u.dig[j+n] + (carry castto(uint32)); ;; ;; /* undo the biasing for remainder */ u = bigshri(u, shift) -> (trim(q), trim(u)) } /* computes b^e % m */ const bigmodpow = {base, exp, mod var r, n r = mkbigint(1) n = 0 while !bigiszero(exp) if (exp.dig[0] & 1) != 0 bigmul(r, base) bigmod(r, mod) ;; bigshri(exp, 1) bigmul(base, base) bigmod(base, mod) ;; -> bigmove(base, r) } /* returns the number of leading zeros */ const nlz = {a : uint32 var n if a == 0 -> 32 ;; n = 0 if a <= 0x0000ffff n += 16 a <<= 16 ;; if a <= 0x00ffffff n += 8 a <<= 8 ;; if a <= 0x0fffffff n += 4 a <<= 4 ;; if a <= 0x3fffffff n += 2 a <<= 2 ;; if a <= 0x7fffffff n += 1 a <<= 1 ;; -> n } /* a <<= b */ const bigshl = {a, b match b.dig.len | 0: -> a | 1: -> bigshli(a, b.dig[0] castto(uint64)) | n: die("shift by way too much\n") ;; } /* a >>= b, unsigned */ const bigshr = {a, b match b.dig.len | 0: -> a | 1: -> bigshri(a, b.dig[0] castto(uint64)) | n: die("shift by way too much\n") ;; } /* a + b, b is integer. FIXME: acually make this a performace improvement */ generic bigaddi = {a, b var bigb : bigint var dig : uint32[2] bigdigit(&bigb, b < 0, b castto(uint64), dig[:]) bigadd(a, &bigb) -> a } generic bigsubi = {a, b : @a::(numeric,integral) var bigb : bigint var dig : uint32[2] bigdigit(&bigb, b < 0, b castto(uint64), dig[:]) bigsub(a, &bigb) -> a } generic bigmuli = {a, b var bigb : bigint var dig : uint32[2] bigdigit(&bigb, b < 0, b castto(uint64), dig[:]) bigmul(a, &bigb) -> a } generic bigdivi = {a, b var bigb : bigint var dig : uint32[2] bigdigit(&bigb, b < 0, b castto(uint64), dig[:]) bigdiv(a, &bigb) -> a } /* a << s, with integer arg. logical left shift. any other type would be illogical. */ generic bigshli = {a, s : @a::(numeric,integral) var off, shift var t, carry iassert(s >= 0, "shift amount must be positive") off = (s castto(uint64)) / 32 shift = (s castto(uint64)) % 32 /* zero shifted by anything is zero */ if a.sign == 0 -> a ;; a.dig = slzgrow(a.dig, 1 + a.dig.len + off castto(size)) /* blit over the base values */ for var i = a.dig.len; i > off; i-- a.dig[i - 1] = a.dig[i - 1 - off] ;; for var i = 0; i < off; i++ a.dig[i] = 0 ;; /* and shift over by the remainder */ carry = 0 for var i = 0; i < a.dig.len; i++ t = (a.dig[i] castto(uint64)) << shift a.dig[i] = (t | carry) castto(uint32) carry = t >> 32 ;; -> trim(a) } /* logical shift right, zero fills. sign remains untouched. */ generic bigshri = {a, s var off, shift var t, carry iassert(s >= 0, "shift amount must be positive") off = (s castto(uint64)) / 32 shift = (s castto(uint64)) % 32 /* blit over the base values */ for var i = 0; i < a.dig.len - off; i++ a.dig[i] = a.dig[i + off] ;; for var i = a.dig.len - off; i < a.dig.len; i++ a.dig[i] = 0 ;; /* and shift over by the remainder */ carry = 0 for var i = a.dig.len; i > 0; i-- t = (a.dig[i - 1] castto(uint64)) a.dig[i - 1] = (carry | (t >> shift)) castto(uint32) carry = t << (32 - shift) ;; -> trim(a) } /* creates a bigint on the stack; should not be modified. */ const bigdigit = {v, isneg : bool, val : uint64, dig v.sign = 1 if isneg val = -val v.sign = -1 ;; if val == 0 v.sign = 0 v.dig = [][:] elif val < Base v.dig = dig[:1] v.dig[0] = val castto(uint32) else v.dig = dig v.dig[0] = val castto(uint32) v.dig[1] = (val >> 32) castto(uint32) ;; } /* trims leading zeros */ const trim = {a var i for i = a.dig.len; i > 0; i-- if a.dig[i - 1] != 0 break ;; ;; a.dig = slgrow(a.dig, i) if i == 0 a.sign = 0 elif a.sign == 0 a.sign = 1 ;; -> a }