Oops, luckily, Julian was eavesdropping! ---------- Forwarded message ---------- From: Julian Ziegler Hunts <julianj.zh@gmail.com> Date: Sat, Apr 23, 2011 at 8:40 PM Subject: Re: [math-fun] How can I accurately calculate (p/q)^n mod 1, for large n To: billgosper@gmail.com

Timing[Do[10^16*Floor[Mod[1.5^k, 1]], {k, 99999}]]

What? That's always 0. Mod[x,1] is ≥0 and <1, so Floor[Mod[x,1]]=0 for all x. The two methods should never differ, since even the numerics know that 0≤Mod[x,1]<1. Maybe you want Floor[10^16*Mod[1.5^k,1]], which diverges from the exact computation at k=21, or, with 10000 (decimal) digits of precision, at k=56758. Julian On Sat, Apr 23, 2011 at 7:32 PM, Bill Gosper <billgosper@gmail.com> wrote:

*From: *

At 05:20 PM 4/23/2011, Dan Asimov wrote:

...

Maybe I should be more specific.

...

Mathematica can do a loop of length 10000, say calculating

100*Floor[Mod[(3/2)^k,1]]

...

for 10001 <= k < 20000,

but when I try to do something like that for

...

a loop of size 100,000 or 1,000,000, it seems to take forever (or hang, I'm not sure which).

--Dan

It's computing the numerators and denominators, e.g., 2^100000, exactly.

In[622]:= Timing[Do[If[10^16*Floor[Mod[(3/2)^k, 1]] != 10^16*Floor[Mod[1.5^k, 1]],Print[k]], {k, 99999}]]

Out[622]= {59.0911, Null}

About 50 times faster:

In[623]:= Timing[Do[10^16*Floor[Mod[1.5^k, 1]], {k, 99999}]]

Out[623]= {1.21972, Null}

Or 60 this way: In[625]:= Block[{t = 10^16},Timing[Do[t*Floor[Mod[1.5^k, 1]], {k, 99999}]]]

Out[625]= {0.997455, Null}

At some point, the two methods will probably differ by 1 (in 10^16). --rwg (Boy, Neil must be busy, not to have beaten me to this one.)