ref: 16b894430a5b71f4956ee305229bbbffd13605d0
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)
t0 = bigdup(&y1)
bigsub(t0, &y0)
umul(z1, t0)
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
}
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
}