ref: 0c678f7b8a8ea747663cc0e6cc18c90fdb1aa470
parent: f922a2c6dfb9a8043f82a19946a06b85469b3be1
author: Ori Bernstein <[email protected]>
date: Wed Jan 22 18:09:27 EST 2014
Implement bigint division.
--- a/libstd/bigint.myr
+++ b/libstd/bigint.myr
@@ -6,6 +6,7 @@
use "option.use"
use "slcp.use"
use "sldup.use"
+use "slfill.use"
use "slpush.use"
use "types.use"
use "utf.use"
@@ -12,7 +13,7 @@
pkg std =
type bigint = struct
- dig : uint32[:] /* little endian, no leading zeros. */
+ dig : uint32[:] /* little endian, no leading zeros. */
sign : int /* -1 for -ve, 0 for zero, 1 for +ve. */
;;
@@ -36,7 +37,6 @@
const bigdivmod : (a : bigint#, b : bigint# -> [bigint#, bigint#])
const bigshl : (a : bigint#, b : bigint# -> bigint#)
const bigshr : (a : bigint#, b : bigint# -> bigint#)
- const bigshra : (a : bigint#, b : bigint# -> bigint#)
/* bigint*int -> bigint ops */
const bigaddi : (a : bigint#, b : int64 -> bigint#)
@@ -45,7 +45,6 @@
const bigdivi : (a : bigint#, b : int64 -> bigint#)
const bigshli : (a : bigint#, b : uint64 -> bigint#)
const bigshri : (a : bigint#, b : uint64 -> bigint#)
- const bigshrai : (a : bigint#, b : uint64 -> bigint#)
;;
const mkbigint = {v@@ -316,19 +315,25 @@
/* a /= b */
const bigdivmod = {a : bigint#, b : bigint# -> [bigint#, bigint#]- /*var u, v /* normalized a and b */*/
- var q/*, qhat /* quotient and estimated quotient */*/
/*
- var rhat /* remainder */
- var s, i, j, t
- var p
- */
- var j
- var carry, b0, aj
+ Implements bigint division using Algorithm D from
+ Knuth: Seminumerical algorithms, Section 4.3.1.
+ */
+ var qhat, rhat, carry, shift
+ var x, y, z, w, p, t /* temporaries */
+ var b0, aj
+ var u, v
+ var m : int64, n : int64
+ 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))
@@ -338,7 +343,7 @@
q.sign = 1
;;
- /* handle single digit divisor separately for performance */
+ /* 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))
@@ -347,14 +352,99 @@
q.dig[j - 1] = (((carry << 32) + aj)/b0) castto(uint32)
carry = (carry << 32) + aj - (q.dig[j-1] castto(uint64))*b0
;;
- for v in q.dig
- ;;
-> (trim(q), trim(mkbigint(carry castto(int32))))
;;
- die("big bigint division not implemented\n")- -> (trim(a), mkbigint(carry castto(int32)))
+
+ u = bigdup(a)
+ v = bigdup(b)
+ q = zalloc()
+ q.dig = slalloc(max(u.dig.len, v.dig.len))
+ m = u.dig.len
+ n = v.dig.len
+
+ shift = nlz(v.dig[n - 1])
+ bigshli(u, shift)
+ bigshli(v, shift)
+ for j = m - n; j >= 0; j--
+ /* load a few temps */
+ x = (u.dig[j + n] castto(uint64)) << 32
+ y = u.dig[j + n - 1] castto(uint64)
+ z = v.dig[n - 1] castto(uint64)
+ w = v.dig[n - 2] castto(uint64)
+
+ /* estimate qhat */
+ qhat = (x + y)/z
+ rhat = (x + y) - (qhat * z)
+:divagain
+ if qhat > 0xfffffffful || qhat * w > ((rhat << 32) + w)
+ qhat--
+ rhat += z
+ if rhat <= ~0
+ 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 & (0xffffffff castto(uint64)))
+ u.dig[i+j] = t castto(uint32)
+ carry = (p >> 32) - (t >> 32);
+ ;;
+ t = x - carry
+ u.dig[j + n] = t castto(uint32)
+ q.dig[j] = qhat castto(uint32)
+ /* adjust */
+ if x < carry
+ 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 */
+ bigshri(u, shift)
+ -> (trim(q), trim(u))
}
+/* 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, bmatch b.dig.len
@@ -373,30 +463,10 @@
;;
}
-/* a >>= b, sign extending */
-const bigshra = {a, b- match b.dig.len
- | 0: -> a
- | 1: -> bigshrai(a, b.dig[0] castto(uint64))
- | n: die("shift by way too much\n")- ;;
-}
-
-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
- ;;
- -> a
-}
-
+/*
+ a << s, with integer arg.
+ logical left shift. any other type would be illogical.
+ */
const bigshli = {a, svar off, shift
var t, carry
@@ -423,19 +493,8 @@
-> trim(a)
}
+/* logical shift right, zero fills */
const bigshri = {a, s- -> bigshrfill(a, s, 0)
-}
-
-const bigshrai = {a, s- if a.sign == -1
- -> bigshrfill(a, s, ~0)
- else
- -> bigshrfill(a, s, 0)
- ;;
-}
-
-const bigshrfill = {a, s, fillvar off, shift
var t, carry
var i
@@ -448,7 +507,7 @@
a.dig[i] = a.dig[i + off]
;;
for i = a.dig.len; i < a.dig.len + off; i++
- a.dig[i] = fill
+ a.dig[i] = 0
;;
/* and shift over by the remainder */
carry = 0
@@ -459,3 +518,20 @@
;;
-> trim(a)
}
+
+/* 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
+ ;;
+ -> a
+}
+