# Fast Karatsuba Squaring, how to?

I want to improve timing of squaring for big numbers.

So I asked What is the fastest way to split a number into parts and What is the most efficient way to add zeros to end of a number with several answers posted.

Karatsuba multiplication must benefit from smaller size multiplication but mine suffer from this.

  FastSquare[a_] := (TotalLen = IntegerLength[a];
If[TotalLen == 1, Return[a*a]];
e = Quotient[TotalLen, 2];
{A, B} = QuotientRemainder[a, 10^e];
Result = A*A*10^(2*e) + 2*A*B*10^(e) + B*B;
Return[Result])


And timings are:

In[21]:= First[Timing[(2^200000000 + 1)^2]]

Out[21]= 2.281

In[20]:= First[Timing[FastSquare[2^200000000 + 1]]]

Out[20]= 15.86

(Debug) In[25]:= FastSquare[2^200000000 + 1]; // RuntimeToolsProfile

Calls   Time    Evaluation
1   15.766  FastSquare[2^200000000+1];
1   15.766  FastSquare[2^200000000+1]
1   15.75   FastSquare[a_]
1   15.75   TotalLen=IntegerLength[a];If[TotalLen==1,Return[a a]];e=Quotient[TotalLen,2];{A,B}=QuotientRemainder[a,10^e];Result=A A Power[<<2>>]+2 A B Power[<<2>>]+B B;Return[Result]
1   10.422  Result=A A 10^Times[<<2>>]+2 A B 10^e+B B
1   10.422  A A 10^(2 e)+2 A B 10^e+B B
1   5.328   {A,B}=QuotientRemainder[a,10^e]
1   5.328   QuotientRemainder[a,10^e]
1   5.203   A A 10^(2 e)
1   4.141   2 A B 10^e
1   1.344   10^(2 e)
1   1.031   B B
1   0.656   10^e
1   0.625   10^e
1   0.016   2^200000000+1


As you can see 5 seconds are spent on splitting the numbers (pre-processing). And each half size multiplication took double of original squaring time because of adding zeros.

Already I've shown here that it is quite possible to improve the speed of some of Mathematica functions by rewriting them.

Anyway to reduce the timing of FastSquare function to less than timing of ^2 ?

-

You are losing hugely due to a base 10 implementation. Integers are manipulated in base 2 in Mathematica. So you would want a base 2 variant to get any reasonable behavior.

Here is one way to code it. There might be tweaks that improve it.

fastSquare[a_] := Catch[Module[
{len = BitLength[a], len2, hi, lo, h2, l2, hl},
If[len < 100, Throw[a^2]];
len2 = BitShiftRight[len];
hi = BitShiftRight[a, len2];
lo = a - BitShiftLeft[hi, len2];
h2 = hi^2;
l2 = lo^2;
hl = hi*lo;
BitShiftLeft[h2, 2*len2] + BitShiftLeft[hl, len2 + 1] + l2
]]


(1) Might be I'm having a mental glitch here. I confess I do not see why we shift that middle term by len2+1 instead of len2. It seems to work on all examples I tried though. Feel free to enlighten me on what I am either seeing or doing wrong. (I mean above, not my generally misguided ways.)

(2) In theory one should do better by making this recursive. That would involve rewriting hl as ((hi+lo)^2-hi^2-lo^2)/2; the squares could be done recursively and the division is just a right shift.

(3) There is no reason to think any of this will beat native squaring because that can use FFT methods which should be faster still than Karatsuba.

(4) It turns out that item (3) is off base. I do not know why. Here is a test run.

ee = 2^200000000 + 1;

In[228]:= First[Timing[s1 = ee^2]]

Out[228]= 2.190000

In[229]:= First[Timing[s2 = fastSquare[ee]]]

Out[229]= 0.180000

In[230]:= s1 === s2

Out[230]= True


I find it strange that this can beat native squaring. I am guessing that some internal method switches could stand for better tuning. I'll ask around.

--- edit ---

It was shown to me that this is lightning fast due to the special simplicity, base 2, of the input. In general this will not beat native squaring for speed.

--- end edit ---

-

I've come up with the idea of expressing any number n as 2^Exp + B. The idea can be generalized to any other base, too.

FastSquare2[a_] := (
TotalLen = IntegerLength[a];
If[TotalLen == 1, Return[a*a]];
Exp = Ceiling[Log[10, (a)]*Log[10.]/Log[2.]];
If[2^Exp > a, Exp--];
B = a - 2^Exp;
Result = 2^(2 Exp) + 2*(2^Exp)*B + B^2;
Return[Result])


And timing for general number is quite acceptable while it computes any special base 2 number quite fast even faster than Daniel Lichtblau fastSquare function.

p = RandomPrime[10^1000];
e = p^100000;

In[56]:= IntegerLength[e]

Out[56]= 99980477

In[57]:= First[Timing[FastSquare2[e]]]

Out[57]= 4.

Out[58]= First[Timing[fastSquare[e]]] (* Daniel Lichtblau solution *)

Out[58]= 6.469

In[59]:= First[Timing[e^2]]

Out[59]= 3.766

In[60]:= First[Timing[FastSquare2[2^200000000 + 1]]]

Out[60]= 0.125
`
-