Services & Resources / Wolfram Forums
-----
 /
MathGroup Archive
2004
*January
*February
*March
*April
*May
*June
*July
*August
*September
*October
*November
*December
*Archive Index
*Ask about this page
*Print this page
*Give us feedback
*Sign up for the Wolfram Insider

MathGroup Archive 2004

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: 2^20991011-1

  • To: mathgroup at smc.vnet.net
  • Subject: [mg45468] Re: [mg45441] 2^20991011-1
  • From: Daniel Lichtblau <danl at wolfram.com>
  • Date: Fri, 9 Jan 2004 05:20:35 -0500 (EST)
  • References: <200401072231.RAA10329@smc.vnet.net>
  • Sender: owner-wri-mathgroup at wolfram.com

Friedrich Laher wrote:
> 
> how can mathematica do that that fast ( 100 seconds on 2.4GHz Athalon )
> while the following needs almost 100 minutes
> #include<stdio.h>
> /*
> Zahl = 2^20 996 011 -  1 get decimal representation,
> but
> base 10^9 Zahl, 9 dec. digits codes binary in 32Bit.
> 
> 20 996 011*DekadischerLog(2) = 6 320 429, 100 ...
> Zahl
> has 6 320 431 dec. digits, so needs 702 271 unsigned long's ( of 32 bit ).
> 
> 20 996 011 bit's are 656 125 *32 + 11, 2^11-1 = 2074
> ---------------------
> 1 Stunde 38 Min. 3 Sek
> --------------------
> */
> unsigned long Zahl[702271], *DezZeiger;
> int main()
> {
> asm("
> #                       INITIALISIERUNG
> 
>         PUSHL %ebx
>         MOVD %esp,%mm3; LEA Zahl+4*702271,%esp
> 
>         MOVL  $656125,%ecx
>         MOVD    %ecx,%mm2; MOVL $1,%eax; MOVD %eax,%mm1
> 
>         CLC;CMC;  CLD
> 
>         LEA Zahl,%edi; MOVD %edi,%mm0
> 
>         MOV $2047,%edx; MOVL $1000000000,%ebx
> InitLoop:
>         SBB %eax,%eax; DIV %ebx; STOSL; LOOP InitLoop; JMP DezStore
> 
> #       HAUPTSCHLEIFE
> DivLoop:
>         LODSL; DIVL %ebx
> unShortened:
>         STOSL; LOOP DivLoop
> DezStore:
>         PUSHL %edx
>         MOVD  %mm2,%ecx; TESTL %ecx,%ecx; JE endCalc
> 
>         MOVD %mm0,%esi; MOVL %esi,%edi
> 
>         LODSL; SUBL %edx,%edx; DIVL %ebx
> 
>         TESTL %eax,%eax
>         JNE unShortened; PSUBD %mm1,%mm2; LOOP DivLoop; JMP DezStore
> endCalc:
>         MOVL %esp, DezZeiger; MOVD %mm3,%esp
>         POPL %ebx
> ");
> for(;DezZeiger != Zahl+702271;++DezZeiger)printf("%09u ",*DezZeiger);
> 
> printf("\n");
> }


For a number containing n bits, your base conversion code above will
have bit complexity (and hence speed) O(n^2). Divide-and-conquer methods
for base conversion, using fast multiplication/division of large
numbers, are asymptotically faster than that. Specifically, if M(n) is
the time needed to multiply a pair of n-bit numbers, then the complexity
of base conversion is also just O(M(n)).

of Mathematica; in version 5 I think you would get a timing under 30
seconds.

To see more specifically how the divide-and-conquer works, you can do a
crude implementation in Mathematica as below (you could also do a more
refined one at the expense of a bit more coding). I should mention that
I have not fully tested it for correctness but it seems to work as it
ought.

digits[n_,base_] /; base>n := {n}

digits[n_,base_] := With[{halfsize=Ceiling[Log[N[base],n]/2]},
  Join[PadRight[digits[#[[2]],base],halfsize],digits[#[[1]],base]]& [
    Internal`QuotientMod[n,base^halfsize]]
  ]

The runs below will give an indication of the asymptotic behavior.

In[73]:= Timing[dd1 = digits[2^209960-1,10^9];]
Out[73]= {0.61 Second, Null}

In[74]:= Timing[dd2 = digits[2^2099601-1,10^9];]
Out[74]= {7.34 Second, Null}

In[75]:= Timing[dd3 = digits[2^20996011-1,10^9];]
Out[75]= {103.69 Second, Null}


The build in code is faster though not tremendously so:

In[80]:= Timing[ee1 = Reverse[IntegerDigits[2^209960-1,10^9]];]
Out[80]= {0.04 Second, Null}

In[81]:= Timing[ee2 = Reverse[IntegerDigits[2^2099601-1,10^9]];]
Out[81]= {1.6 Second, Null}

In[82]:= Timing[ee3 = Reverse[IntegerDigits[2^20996011-1,10^9]];]
Out[82]= {41.33 Second, Null}

Note that digits is indeed correct at least for the inputs used in the
tests above.

In[83]:= {dd1,dd2,dd3} === {ee1,ee2,ee3}
Out[83]= True


Daniel Lichtblau
Wolfram Research


  • References:
    • 2^20991011-1
      • From: Friedrich Laher <mathefritz@schmieder-laher.de>
  • Prev by Date: Re: 2^20991011-1
  • Next by Date: RE: Re: what's wrong about my package?
  • Previous by thread: 2^20991011-1
  • Next by thread: Re: 2^20991011-1