Karatsuba multiplication

src/bigints.nim

@@ -64,6 +64,7 @@ func initBigInt*(val: BigInt): BigInt =
 const
   zero = initBigInt(0)
   one = initBigInt(1)
+  karatsubaTreshold = 10
 
 func isZero(a: BigInt): bool {.inline.} =
   for i in countdown(a.limbs.high, 0):
@@ -388,7 +389,6 @@ template `-=`*(a: var BigInt, b: BigInt) =
     assert a == 3.initBigInt
   a = a - b
 
-
 func unsignedMultiplication(a: var BigInt, b, c: BigInt) {.inline.} =
   # always called with bl >= cl
   let
@@ -418,6 +418,26 @@ func unsignedMultiplication(a: var BigInt, b, c: BigInt) {.inline.} =
       inc pos
   normalize(a)
 
+func scalarMultiplication(a: var BigInt, b: BigInt, c: uint32) {.inline.} =
+  # always called with bl >= cl
+  let
+    bl = b.limbs.len
+  a.limbs.setLen(bl + 1)
+  var tmp = 0'u64
+
+  for i in 0 ..< bl:
+    tmp += uint64(b.limbs[i]) * uint64(c)
+    a.limbs[i] = uint32(tmp and uint32.high)
+    tmp = tmp shr 32
+
+  a.limbs[bl] = uint32(tmp)
+  normalize(a)
+
+# forward declaration for use in `multiplication`
+func karatsubaMultiplication(a: var BigInt, b, c: BigInt) {.inline.}
+func `shl`*(x: BigInt, y: Natural): BigInt
+func `shr`*(x: BigInt, y: Natural): BigInt
+
 func multiplication(a: var BigInt, b, c: BigInt) =
   # a = b * c
   if b.isZero or c.isZero:
@@ -428,11 +448,69 @@ func multiplication(a: var BigInt, b, c: BigInt) =
     cl = c.limbs.len
 
   if cl > bl:
-    unsignedMultiplication(a, c, b)
+    if bl <= karatsubaTreshold:
+      karatsubaMultiplication(a, c, b)
+    else:
+      unsignedMultiplication(a, c, b)
   else:
-    unsignedMultiplication(a, b, c)
+    if cl <= karatsubaTreshold:
+      karatsubaMultiplication(a, b, c)
+    else:
+      unsignedMultiplication(a, b, c)
   a.isNegative = b.isNegative xor c.isNegative
 
+func karatsubaMultiplication(a: var BigInt, b, c: BigInt) {.inline.} =
+  let
+    bl = b.limbs.len
+    cl = c.limbs.len
+  let n = max(bl, cl)
+  if bl == 1:
+    # base case : multiply the only limb with each limb of second term
+    scalarMultiplication(a, c, b.limbs[0])
+    return 
+  if cl == 1:
+    scalarMultiplication(a, b, c.limbs[0])
+    return
+  if bl < karatsubaTreshold:
+    if cl <= bl:
+      unsignedMultiplication(a, b, c)
+    else:
+      unsignedMultiplication(a, c, b)
+    return
+  if cl < karatsubaTreshold:
+    if bl <= cl:
+      unsignedMultiplication(a, c, b)
+    else:
+      unsignedMultiplication(a, b, c)
+    return
+  let k = n shr 1 # should it be ceil(n/2) ?
+  var
+    low_b, high_b, low_c, high_c: BigInt
+  # Decompose `b` and `c` in two parts of (almost) equal length
+  low_b.limbs = b.limbs[0 .. k-1]
+  high_b.limbs = b.limbs[k .. ^1]
+  low_c.limbs = c.limbs[0 .. k-1]
+  high_c.limbs = c.limbs[k .. ^1]
+
+  # subtractive version of Karatsuba's algorithm :
+  # limit carry handling in opposition to the additive version
+  var
+    lowProduct, highProduct, A3, A4, A5, middleTerm: BigInt = zero
+  karatsubaMultiplication(lowProduct, low_b, low_c)
+  karatsubaMultiplication(highProduct, high_b, high_c)
+  A3 = low_b - high_b # Additive variant of Karatsuba
+  A4 = low_c - high_c # would add them
+  if A4.limbs.len >= A3.limbs.len:
+    multiplication(A5, abs(A4), abs(A3))
+  else:
+    multiplication(A5, abs(A3), abs(A4))
+  middleTerm = lowProduct + highProduct + A5
+  a.limbs[0 .. k - 1] = lowProduct.limbs
+  # a += (middleTerm shr k) + (highProduct shr (2*k))
+  a.limbs[k .. 2*k-1] = middleTerm.limbs
+  a.limbs[2*k .. 3*k-1] = highProduct.limbs
+
+
 func `*`*(a, b: BigInt): BigInt =
   ## Multiplication for `BigInt`s.
   runnableExamples: