From: Dan Asimov <dasimov@earthlink.net>

Maybe I should be more specific.

Mathematica can do a loop of length 10000, say calculating

Floor[100*Mod[(3/2)^k,1]] [correction from later post]

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

Oh, you mean a million large, not 2^4096 large. Could it be that Mathematica's rational arithmetic is trying to reduce the fraction at every step in the exponentiation? Try maybe Floor[100*Mod[3^k,2^k]/2^k] ? For that matter, try 3^k

From: Bill Gosper <billgosper@gmail.com>

Oops, luckily, Julian was eavesdropping!

---------- Forwarded message ---------- From: Julian Ziegler Hunts <julianj.zh@gmail.com>

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.

You mean diverges from the precision Dan asked for, right? To get the *exact* answer you need about three times as many digits (and when q = 2 you would need to calculate in binary). If you want d digits of precision at the end, log10(p/q)*n + d digits is enough, right? (Seems to be what Julian is calculating.) Eg. for k = 1,000,000 you need 176,092 + d digits. Btw, we're talking about precision *mod 1* right? So .995 might come out .99 or .00 and that would be okay?

From: Dan Asimov <dasimov@earthlink.net>

Yes, Julian was right, in fact what I was calculating was

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

for long sequences of consecutive k.

That must be a lot faster to calculate in an incremental way, even in Mathematica, right? I don't have Mathematica but I wondered what might take it so long, so I timed something similar in Python. Python has big ints but only IEEE floats built in, so I took (p^n mod q^n) / q^n (This doesn't have the wraparound problem, it just takes more time and space.) def p_over_q_to_the_n_mod_1(p, q, n): pn = p ** n qn = q ** n m = pn % qn # Can't convert a number > 2**1023 to IEEE float, so... nbits = int(log(qn, 2)) if nbits > 1000: scale = 2 ** (nbits - 1000) m = m / scale qn = qn / scale return float(m) / qn $ ./p_over_q.py 3 2 1000000 pn = p^n... 0.974992990494 sec. qn = q^n... 0.0310521125793 sec. m = pn mod qn... 22.5215368271 sec. scale down... 0.12483215332 sec. r = float(m) / qn... 8.70227813721e-05 sec. Total: 23.6525011063 sec. 0.327352766151 The odd thing in Python's case is that it spent most of its time taking a number mod a power of two! Adding a special case for that... if q == 2: m = pn & (qn - 1) # mod using binary mask else: m = pn % qn $ ./p_over_q.py 3 2 1000000 pn = p^n... 0.976704120636 sec. qn = q^n... 0.0312979221344 sec. m = pn mod qn... 0.00104093551636 sec. scale down... 0.117168188095 sec. r = float(m) / qn... 8.20159912109e-05 sec. Total: 1.12629318237 sec. 0.327352766151 If that helps... --Steve