Fast remainder

[Albert Chan]

(05-12-2021 03:10 PM)Albert Chan Wrote:  If n = 128 = 2^7, then log2(d) < 7

F ≥ N + log2(d)
max(N) = 32 - 7 = 25

If we combined 2 fastmod into 1, we have:

x = ((d + x[i-1])*n + x[i])*n < 2*n^3 = 2^(1 + 3*log2(n))

N ≤ 1+3*7 = 22 ≤ 25, thus we are safe.

Instead of using general rule, we can get more closely to the fastmod limit.

M = 2^F - 1
c = floor(M/d) + 1 = (M+1)/d + (d-1-M%d)/d = 2^F/d + cerr

fastmod(x, d, c) = floor((x*c & M) * d / (M+1))

For fastmod to work, x must not be too big to cause error ≥ 1, where floor cannot correct.

x * cerr * d / (M+1) < 1
x * (d-1 - M%d) ≤ M

max(x) = floor(M / (d-1 - M%d))

Note that max(x) goes infinite when RHS denominator goes 0, when d = pow-of-2.
To confirm the equation, lets redo previous post examples, by direct calculation.

>>> M = (1<<32) - 1 # F = 32
>>> maxx = lambda d: int(M // (d-1 - M%d))
>>> print hex(maxx(17)), hex(maxx(29))

0xfffffff  0x13b13b13

Now, we are ready to get closer fastmod limit.

Code:

for n in range(128, 256):
    try: limit = maxx(n-1)
    except ZeroDivisionError: continue
    if 2*n*n*n > limit: print 'max base =',n-1; break

max base = 217

Actually, looping it all is unnecessary.
We maximize denominator (assume M%d=0), replace integer divide with float divide.

2*n^3 = M / (n-2)
n ≈ ⁴√(M/2) ≈ 215.3

Consider effect of n-2 ≈ n, this is below solution for actual n.
Start from n=216, we hit the sweet spot in 3 tries (n=218 failed max(x) test)

---

This is the patch to code (post #2)
For n=60, timing goes from 14.5 sec, down to 12.3 sec.

Code:

33c33
<         d = fastmod((d + x[i-1])*n, k, c);
---
>         d = (d + x[i-1])*n;
79,80c79,80
<     if (n-1 < 181) {puzzle(n); return 0;}
<     return printf("Usage %s n, 0 < n <= 181\n", argv[0]);
---
>     if (n-1 < 217) {puzzle(n); return 0;}
>     return printf("Usage %s n, 0 < n <= 217\n", argv[0]);