summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorOri Bernstein <ori@eigenstate.org>2018-03-17 01:26:17 -0700
committerOri Bernstein <ori@eigenstate.org>2018-03-17 01:29:58 -0700
commitf127b8250ecbd5f44622f9bf620af933f7102e78 (patch)
tree0c88b717a00971c7533b7912f4c6f9b5b4afef41
parentdc9fbed2061a1c5b6de5ce1627aee10c189d123f (diff)
downloadmc-f127b8250ecbd5f44622f9bf620af933f7102e78.tar.gz
Implement karatsuba multiplication.
This is a nice optimization. In the most optimal case, it makes this stupid benchmark about 10 times faster: a = std.mkbigint(17193) b = std.mkbigint(1091327) for var i = 0; i < 12; i++ mul(a, b) mul(b, a) ;; It also makes a decent change (about 10%) on `mbld bench` for bigfactorial, in spite of us mostly not falling into the karatsuba case.
-rw-r--r--bench/bigfactorial.myr1
-rw-r--r--lib/std/bigint.myr92
2 files changed, 82 insertions, 11 deletions
diff --git a/bench/bigfactorial.myr b/bench/bigfactorial.myr
index cbaab6c..a893867 100644
--- a/bench/bigfactorial.myr
+++ b/bench/bigfactorial.myr
@@ -7,6 +7,7 @@ const main = {
[.name="bigfactorial-100", .fn={ctx; bigfact(100)}],
[.name="bigfactorial-1000", .fn={ctx; bigfact(1000)}],
[.name="bigfactorial-10000", .fn={ctx; bigfact(10000)}],
+ [.name="bigfactorial-20000", .fn={ctx; bigfact(20000)}],
][:])
}
diff --git a/lib/std/bigint.myr b/lib/std/bigint.myr
index 78a5422..29286c4 100644
--- a/lib/std/bigint.myr
+++ b/lib/std/bigint.myr
@@ -92,6 +92,7 @@ pkg std =
extern const put : (fmt : byte[:], args : ... -> size)
const Base = 0x100000000ul
+const Kmin = 64
generic mkbigint = {v : @a :: integral,numeric @a
var a
@@ -124,7 +125,7 @@ const bigfrombytes = {isneg, v
;;
for i = 0; i + 4 <= v.len; i += 4
- std.slpush(&a.dig, \
+ slpush(&a.dig, \
(v[i + 0] << 0 : uint32) | \
(v[i + 1] << 8 : uint32) | \
(v[i + 2] << 16 : uint32) | \
@@ -135,7 +136,7 @@ const bigfrombytes = {isneg, v
off = i & 0x3
last |= (v[off] : uint32) << (8 *off)
;;
- std.slpush(&a.dig, last)
+ slpush(&a.dig, last)
-> trim(a)
}
@@ -170,7 +171,7 @@ const bigsteal = {d, s
}
const bigclear = {v
- std.slfree(v.dig)
+ slfree(v.dig)
v.sign = 0
v.dig = [][:]
-> v
@@ -259,7 +260,7 @@ const bigparse = {str
fit in one digit.
*/
v = mkbigint(1)
- for c : std.bychar(str)
+ for c : bychar(str)
if c == '_'
continue
;;
@@ -515,10 +516,7 @@ const usub = {a, b
/* a *= b */
const bigmul = {a, b
- var i, j
- var ai, bj, wij
- var carry, t
- var w
+ var s
if a.sign == 0 || b.sign == 0
a.sign = 0
@@ -526,10 +524,83 @@ const bigmul = {a, b
a.dig = [][:]
-> a
elif a.sign != b.sign
- a.sign = -1
+ s = -1
else
- a.sign = 1
+ s = 1
+ ;;
+
+ umul(a, b)
+
+ a.sign = s
+ -> trim(a)
+}
+
+const umul = {a, b
+ var r
+
+ if a.dig.len < Kmin || b.dig.len < 2
+ 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++
@@ -546,7 +617,6 @@ const bigmul = {a, b
;;
slfree(a.dig)
a.dig = w
- -> trim(a)
}
const bigdiv = {a : bigint#, b : bigint# -> bigint#