[m-rev.] for review: Add integer multiplication based on Karatsuba's algorithm

[Paul Bone]

From: Matthias Güdemann <matthias.guedemann at googlemail.com>

This patch is contributed by Matthias Guedemann as pull request #13.

For review by Julien, he understands arbitrary precision arithmetic better
than I do and has offered to take a look.  I intend to take a look at the
parallelisation and possibly tweek it.

Thanks.

---

  library/integer.m:
    Adds Karatsuba based fast integer algorithm in the function
    `pos_mul_karatsuba' which replaces the call to `pos_mul_list' in `pos_mul'.

    It is a divide and conquer algorithm which uses O(n^1.585) basic operations
    instead of O(n^2). The function expects two integer variables, the first
    must be smaller than the second. It falls back to quadratic multiplication
    if the larger factor has less than `karatsuba_threshold' digits.

    If compiled in a parallel grade, the function uses parallelism for the
    divide and conquer part for numbers with more than
    `karatsuba_parallel_threshold' digits.

    The current values for `karatsuba_threshold' is 35,
    `karatsuba_parallel_threshold' is 10 times this value.
---
 library/integer.m | 52 ++++++++++++++++++++++++++++++++++++++++++++++++++--
 1 file changed, 50 insertions(+), 2 deletions(-)

diff --git a/library/integer.m b/library/integer.m
index 1b6bddb..500ab3a 100644
--- a/library/integer.m
+++ b/library/integer.m
@@ -748,11 +748,59 @@ diff_pairs_equal(Div, [X | Xs], [Y | Ys], [Mod | TailDs]) :-
 
 pos_mul(i(L1, Ds1), i(L2, Ds2)) =
     ( L1 < L2 ->
-        pos_mul_list(Ds1, integer.zero, i(L2, Ds2))
+        pos_mul_karatsuba(i(L1, Ds1), i(L2, Ds2))
     ;
-        pos_mul_list(Ds2, integer.zero, i(L1, Ds1))
+        pos_mul_karatsuba(i(L2, Ds2), i(L1, Ds1))
     ).
 
+% use quadratic multiplication for less than threshold digits
+:- func karatsuba_threshold = int.
+
+karatsuba_threshold = 35.
+
+% use parallel execution if number of digits of split numbers is larger
+% then this threshold
+:- func karatsuba_parallel_threshold = int.
+
+karatsuba_parallel_threshold = karatsuba_threshold * 10.
+
+% Karatsuba / Toom-2 multiplication in O(n^1.585)
+:- func pos_mul_karatsuba(integer, integer) = integer.
+
+pos_mul_karatsuba(i(L1, Ds1), i(L2, Ds2)) = Res :-
+    ( L1 < karatsuba_threshold ->
+        Res = pos_mul_list(Ds1, integer.zero, i(L2, Ds2))
+    ;
+      ( L2 < L1 ->
+          error("integer.pos_mul_karatsuba: second factor smaller")
+      ;
+          Middle = L2 div 2,
+          HiDigits = L2 - Middle,
+          HiDigitsSmall = max(0, L1 - Middle),
+          % split Ds1 in [L1 - Middle];[Middle] digits if L1 > Middle
+          % or leave as [L1] digits
+          list.split_upto(HiDigitsSmall, Ds1, Ds1_upper, Ds1_lower),
+          % Split Ds2 in [L2 - Middle; Middle] digits
+          list.split_upto(HiDigits, Ds2, Ds2_upper, Ds2_lower),
+          LoDs1 = i(min(L1, Middle), Ds1_lower),
+          LoDs2 = i(Middle, Ds2_lower),
+          HiDs1 = i(HiDigitsSmall, Ds1_upper),
+          HiDs2 = i(HiDigits, Ds2_upper),
+          ( Middle > karatsuba_parallel_threshold ->
+            ( Res0 = pos_mul(LoDs1, LoDs2) &
+              Res1 = pos_mul(LoDs1 + HiDs1, LoDs2 + HiDs2) &
+              Res2 = pos_mul(HiDs1, HiDs2))
+          ;
+            ( Res0 = pos_mul(LoDs1, LoDs2),
+              Res1 = pos_mul(LoDs1 + HiDs1, LoDs2 + HiDs2),
+              Res2 = pos_mul(HiDs1, HiDs2))
+          ),
+          Res = big_left_shift(Res2, 2*Middle*log2base) +
+                big_left_shift(Res1 - (Res2 + Res0), Middle*log2base) + Res0
+      )
+    ).
+
+
 :- func pos_mul_list(list(digit), integer, integer) = integer.
 
 pos_mul_list([], Carry, _Y) = Carry.
-- 
1.8.5.3