ref: 74d91a0021de012908cfdb35fb61a1473a376130
dir: /lib/std/bigint.myr/
use "alloc" use "chartype" use "cmp" use "die" use "errno" use "extremum" use "hasprefix" use "memops" use "option" use "slcp" use "sldup" use "slfill" use "slpush" use "striter" use "types" use "utf" 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 -> bigint#) :: numeric,integral @a const bigfrombytes : (isneg : bool, v : byte[:] -> 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 bigsteal : (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 bigiseven : (a : bigint# -> bool) const bigcmp : (a : bigint#, b : bigint# -> order) const bigcmpabs : (a : bigint#, b : bigint# -> order) generic bigcmpi : (a : bigint#, b : @a -> order) :: numeric,integral @a /* shorthand for comparisons */ const bigeq : (a : bigint#, b : bigint# -> bool) const biglt : (a : bigint#, b : bigint# -> bool) const bigle : (a : bigint#, b : bigint# -> bool) const biggt : (a : bigint#, b : bigint# -> bool) const bigge : (a : bigint#, b : bigint# -> bool) generic bigeqi : (a : bigint#, b : @a -> bool) :: numeric,integral @a generic biglti : (a : bigint#, b : @a -> bool) :: numeric,integral @a generic biglei : (a : bigint#, b : @a -> bool) :: numeric,integral @a generic biggti : (a : bigint#, b : @a -> bool) :: numeric,integral @a generic biggei : (a : bigint#, b : @a -> bool) :: numeric,integral @a /* 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 bigand : (a : bigint#, b : bigint# -> bigint#) const bigor : (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 -> bigint#) :: integral,numeric @a generic bigsubi : (a : bigint#, b : @a -> bigint#) :: integral,numeric @a generic bigmuli : (a : bigint#, b : @a -> bigint#) :: integral,numeric @a generic bigdivi : (a : bigint#, b : @a -> bigint#) :: integral,numeric @a generic bigmodi : (a : bigint#, b : @a -> bigint#) :: integral,numeric @a generic bigshli : (a : bigint#, b : @a -> bigint#) :: integral,numeric @a generic bigshri : (a : bigint#, b : @a -> bigint#) :: integral,numeric @a generic bigandi : (a : bigint#, b : @a -> bigint#) :: integral,numeric @a generic bigori : (a : bigint#, b : @a -> bigint#) :: integral,numeric @a //const bigpowi : (a : bigint#, b : uint64 -> bigint#) //const bigmodpowi : (b : bigint#, e : bigint#, m : bigint# -> bigint#) /* information about bigints */ const bigbitcount : (a : bigint# -> size) ;; /* put for debugging */ extern const put : (fmt : byte[:], args : ... -> size) const Base = 0x100000000ul const Kmin = 64 generic mkbigint = {v : @a :: integral,numeric @a var a var val a = zalloc() if v < 0 a.sign = -1 v = -v elif v > 0 a.sign = 1 ;; val = (v : uint64) slpush(&a.dig, (val : uint32)) if val >= Base slpush(&a.dig, (val/Base : uint32)) ;; -> trim(a) } const bigfrombytes = {isneg, v var i, off, a, last a = mkbigint(0) if isneg a.sign = -1 else a.sign = 1 ;; for i = 0; i + 4 <= v.len; i += 4 slpush(&a.dig, \ (v[i + 0] << 0 : uint32) | \ (v[i + 1] << 8 : uint32) | \ (v[i + 2] << 16 : uint32) | \ (v[i + 3] << 24 : uint32)) ;; last = 0 for i; i < v.len; i++ off = i & 0x3 last |= (v[off] : uint32) << (8 *off) ;; slpush(&a.dig, last) -> 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 bigsteal = {d, s bigmove(d, s); bigfree(s) -> d } const bigclear = {v 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") ;; n = 0 if base == 0 b = mkbigint(10) else b = mkbigint(base) ;; 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 val : int, base var v, b var a var s = 1 if hasprefix(str, "-") s = -1 str = str[1:] ;; 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) for c : bychar(str) if c == '_' continue ;; val = charval(c, base) if val < 0 || val > base bigfree(a) bigfree(b) bigfree(v) -> `None ;; v.dig[0] = (val : uint32) if val == 0 v.sign = 0 else v.sign = s ;; bigmul(a, b) bigadd(a, v) ;; bigfree(b) bigfree(v) -> `Some a } const bigiszero = {v -> v.dig.len == 0 } const bigiseven = {v -> v.dig.len == 0 || v.dig[0] & 1 == 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 } const biglt = {a, b match bigcmp(a, b) | `Before: -> true | _: -> false ;; } const bigle = {a, b match bigcmp(a, b) | `Before: -> true | `Equal: -> true | _: -> false ;; } const biggt = {a, b match bigcmp(a, b) | `After: -> true | _: -> false ;; } const bigge = {a, b match bigcmp(a, b) | `After: -> true | `Equal: -> true | _: -> false ;; } generic bigeqi = {a, b var v var dig : uint32[2] bigdigit(&v, b < 0, (b : uint64), dig[:]) -> bigeq(a, &v) } generic biglti = {a, b match bigcmpi(a, b) | `Before: -> true | _: -> false ;; } generic biglei = {a, b match bigcmpi(a, b) | `Before: -> true | `Equal: -> true | _: -> false ;; } generic biggti = {a, b match bigcmpi(a, b) | `After: -> true | _: -> false ;; } generic biggei = {a, b match bigcmpi(a, b) | `After: -> true | `Equal: -> true | _: -> false ;; } generic bigcmpi = {a, b var v var dig : uint32[2] bigdigit(&v, b < 0, (b : uint64), dig[:]) -> bigcmp(a, &v) } const bigcmp = {a, b var da, db, sa, sb sa = (a.sign : int64) sb = (b.sign : 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] : int64) db = (b.dig[i - 1] : int64) if da != db -> signedorder(sa * (da - db)) ;; ;; ;; -> `Equal } const bigcmpabs = {a, b if a.dig.len < b.dig.len -> `Before elif a.dig.len > b.dig.len -> `After else for var i = a.dig.len; i > 0; i-- var da = (a.dig[i - 1] : int64) var db = (b.dig[i - 1] : int64) if da != db -> signedorder(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 bigcmpabs(a, b) | `Before: /* (1) + (-2) or (-1) + (2) */ a.sign = b.sign var r = bigdup(b) usub(r, a) -> bigsteal(a, r) | `After: /* (2) + (-1) or (-2) + (1) */ -> usub(a, b) | `Equal: -> bigclear(a) ;; ;; } /* 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 */ slzgrow(&a.dig, n + 1) for i = 0; i < n; i++ v = (a.dig[i] : uint64) + carry; if i < b.dig.len v += (b.dig[i] : uint64) ;; if v >= Base carry = 1 else carry = 0 ;; a.dig[i] = (v : uint32) ;; a.dig[i] += (carry : 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 bigcmpabs(a, b) | `Before: /* (-1) - (-2) or (1) - (2) */ var r = bigdup(b) usub(r, a) r.sign = -1 * b.sign -> bigsteal(a, r) | `After: /* (-2) - (-1) or (2) - (1) */ -> 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] : int64) - carry if i < b.dig.len v -= (b.dig[i] : int64) ;; if v < 0 carry = 1 else carry = 0 ;; a.dig[i] = (v : uint32) ;; -> trim(a) } /* a *= b */ const bigmul = {a, b var s if a.sign == 0 || b.sign == 0 a.sign = 0 slfree(a.dig) a.dig = [][:] -> a elif a.sign != b.sign s = -1 else s = 1 ;; umul(a, b) a.sign = s -> trim(a) } const umul = {a, b var r if a.dig.len < Kmin || b.dig.len < Kmin smallmul(a, b) else r = mkbigint(0) kmul(r, a, b) bigmove(a, r) ;; } const kmul = {r, a, b var x0, x1, y0, y1, n var z0, z1, z2, t0 if a.dig.len < b.dig.len t0 = a a = b b = t0 ;; n = min(a.dig.len / 2, b.dig.len - 1) x0 = [.sign=1, .dig=a.dig[:n]] x1 = [.sign=1, .dig=a.dig[n:]] y0 = [.sign=1, .dig=b.dig[:n]] y1 = [.sign=1, .dig=b.dig[n:]] z0 = bigdup(&x0) trim(z0) umul(z0, &y0) z2 = bigdup(&x1) trim(z2) umul(z2, &y1) z1 = bigdup(&x0) trim(z1) bigsub(z1, &x1) var z1sign = z1.sign t0 = bigdup(&y1) bigsub(t0, &y0) umul(z1, t0) z1.sign = z1sign * t0.sign bigadd(z1, z0) bigadd(z1, z2) bigshli(z1, 32*n) bigshli(z2, 32*2*n) bigclear(r) bigadd(r, z0) bigadd(r, z1) bigadd(r, z2) bigfree(z0) bigfree(z1) bigfree(z2) bigfree(t0) } const smallmul = {a, b var i, j var ai, bj, wij var carry, t var w 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] : uint64) bj = (b.dig[j] : uint64) wij = (w[i+j] : uint64) t = ai * bj + wij + carry w[i+j] = (t : uint32) carry = t >> 32 ;; w[i + j] = (carry : 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) -> bigsteal(a, q) } const bigmod = {a : bigint#, b : bigint# -> bigint# var q, r (q, r) = bigdivmod(a, b) bigfree(q) -> bigsteal(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 pt, tt 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] : uint64)) for j = a.dig.len; j > 0; j-- aj = ((a.dig[j - 1] : uint64)) q.dig[j - 1] = ((((carry << 32) + aj)/b0) : uint32) carry = (carry << 32) + aj - (q.dig[j-1] : uint64)*b0 ;; q = trim(q) -> (q, trim(mkbigint(carry))) ;; u = bigdup(a) v = bigdup(b) m = u.dig.len n = v.dig.len /* normalize */ shift = nlz(v.dig[n - 1]) bigshli(u, shift) bigshli(v, shift) slzgrow(&u.dig, u.dig.len + 1) /* Since we're little endian, we iterate backwards from Knuth */ for j = m - n; j >= 0; j-- /* load a few temps for less casting */ x = (u.dig[j + n] : uint64) y = (u.dig[j + n - 1] : uint64) z = (v.dig[n - 1] : uint64) w = (v.dig[n - 2] : uint64) t = (u.dig[j + n - 2] : uint64) /* calculate 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] : uint64)) t = (u.dig[i+j] : uint64) - carry - (p % Base) u.dig[i+j] = (t : uint32) tt = (t : int64) >> 32 pt = (p >> 32) carry = ((pt : int64) - (tt : int64) : uint64) ;; t = (u.dig[j + n] : uint64) - carry u.dig[j + n] = (t : uint32) q.dig[j] = (qhat : uint32) /* adjust */ if (t : int64) < 0 q.dig[j]-- carry = 0 for i = 0; i < n; i++ t = (u.dig[i+j] : uint64) + (v.dig[i] : uint64) + carry u.dig[i+j] = (t : uint32) carry = t >> 32 ;; u.dig[j+n] = u.dig[j+n] + (carry : uint32) ;; ;; /* undo the biasing for remainder */ bigshri(u, shift) bigfree(v) -> (trim(q), trim(u)) } const bigand = {a, b if a.dig.len > b.dig.len slzgrow(&b.dig, a.dig.len) ;; if b.dig.len > a.dig.len slzgrow(&a.dig, b.dig.len) ;; for var i = 0; i < a.dig.len; i++ a.dig[i] &= b.dig[i] ;; -> trim(a) } const bigor = {a, b slzgrow(&a.dig, max(a.dig.len, b.dig.len)) for var i = 0; i < min(a.dig.len, b.dig.len); i++ a.dig[i] |= b.dig[i] ;; -> trim(a) } /* 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) ;; -> bigsteal(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] : 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]: 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 : uint64), dig[:]) bigadd(a, &bigb) -> a } generic bigsubi = {a, b : @a :: numeric,integral @a var bigb : bigint var dig : uint32[2] bigdigit(&bigb, b < 0, (b : uint64), dig[:]) bigsub(a, &bigb) -> a } generic bigmuli = {a, b var bigb : bigint var dig : uint32[2] bigdigit(&bigb, b < 0, (b : uint64), dig[:]) bigmul(a, &bigb) -> a } generic bigdivi = {a, b var bigb : bigint var dig : uint32[2] bigdigit(&bigb, b < 0, (b : uint64), dig[:]) bigdiv(a, &bigb) -> a } generic bigmodi = {a, b var bigb : bigint var dig : uint32[2] bigdigit(&bigb, b < 0, (b : uint64), dig[:]) bigmod(a, &bigb) -> a } /* a << s, with integer arg. logical left shift. any other type would be illogical. */ generic bigshli = {a, s : @a :: numeric,integral @a var off, shift var t, carry iassert(s >= 0, "shift amount must be positive") off = (s : uint64) / 32 shift = (s : uint64) % 32 /* zero shifted by anything is zero */ if a.sign == 0 -> a ;; slzgrow(&a.dig, (1 + a.dig.len + off : 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] : uint64) << shift a.dig[i] = (t | carry: 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 : uint64) / 32 shift = (s : uint64) % 32 if off > a.dig.len a.dig = a.dig[:0] else /* blit over the base values */ for var i = 0; i < a.dig.len - off; i++ a.dig[i] = a.dig[i + off] ;; a.dig = a.dig[:a.dig.len - off] /* and shift over by the remainder */ carry = 0 for var i = a.dig.len; i > 0; i-- t = ((a.dig[i - 1] : uint64)) a.dig[i - 1] = (carry | (t >> shift): uint32) carry = t << (32 - shift) ;; ;; -> trim(a) } generic bigandi = {a, b var v var dig : uint32[2] bigdigit(&v, b < 0, (b : uint64), dig[:]) -> bigand(a, &v) } generic bigori = {a, b var v var dig : uint32[2] bigdigit(&v, b < 0, (b : uint64), dig[:]) -> bigor(a, &v) } /* 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 : uint32) else v.dig = dig v.dig[0] = (val : uint32) v.dig[1] = ((val >> 32) : 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 = a.dig[:i] if i == 0 a.sign = 0 ;; -> a } const bigbitcount = {a var top, len, mask len = 32*a.dig.len if len > 0 top = a.dig[a.dig.len - 1] mask = 1 << 31 while top & mask == 0 len-- mask >>= 1 ;; ;; -> len }