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

[Julien Fischer]

Hi,

On Tue, 11 Feb 2014, Paul Bone wrote:

> 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.

Paul: I assume you won't look at the parallelisation stuff until after it is
committed?

The comments below apply mainly to whoever is handling this change from
now on.  (I'm happy to deal with the pull request and make the changes
below if that's the easiest thing to do.)

>  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

s/two integer/two positive integer/

>    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

In the possible improvmenents listed at the beginning of the implementation
section, the third entry, Karatsuba multiplication, should be marked as being
done.

> @@ -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

Sentences being with a capital letter and end with a full stop.

> +:- func karatsuba_threshold = int.
> +
> +karatsuba_threshold = 35.
> +
> +% use parallel execution if number of digits of split numbers is larger
> +% then this threshold

Ditto.  (And elsewhere).

> +:- 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))
> +    ;

Fix the indenation of the else branch here so that it aligns with
the then branch above.

> +      ( 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),

Avoid swapping between CamelCase and undscores in variables names there,
e.g. s/Ds1_upper/Ds1Upper/ etc.

> +          % 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))
> +          ),

The redundent parentheses around the then and else goals should be removed.

The diff looks fine otherwise, although the test case in
tests/hard_coded/integer_test also needs to be extended appropriately.

Thanks for implementing this!

Cheers,
Julien.