9
$\begingroup$

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]; // RuntimeTools`Profile 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 ?

Improve this question
$\endgroup$
2
  • $\begingroup$ Try RPF method by Krishil Sheth , it's faster than karatsuba for squaring numbers upto 100x faster $\endgroup$ Commented May 6 at 18:52
  • 2
    $\begingroup$ @KRISHILSHETH RPF is not especially useful if there is no available description or reference implementation. $\endgroup$ Commented May 7 at 1:33

2 Answers 2

Reset to default
15
$\begingroup$

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

Improve this answer
$\endgroup$
3
$\begingroup$

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
Improve this answer
$\endgroup$

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.