a, b, p = findprod((1<<i)+1, +1, (i, i+1))
print "mul", hexstr(a), hexstr(b), hexstr(p)
+# Bare tests of division/modulo.
+prefixes = [2**63, int(2**63.5), 2**64-1]
+for nsize in range(20, 200):
+ for dsize in range(20, 200):
+ for dprefix in prefixes:
+ d = sqrt(3<<(2*dsize)) + (dprefix<<dsize)
+ for nprefix in prefixes:
+ nbase = sqrt(3<<(2*nsize)) + (nprefix<<nsize)
+ for modulus in sorted({-1, 0, +1, d/2, nbase % d}):
+ n = nbase - (nbase % d) + modulus
+ if n < 0:
+ n += d
+ assert n >= 0
+ print "divmod", hexstr(n), hexstr(d), hexstr(n/d), hexstr(n%d)
+
+# Simple tests of modmul.
+for ai in range(20, 200, 60):
+ a = sqrt(3<<(2*ai-1))
+ for bi in range(20, 200, 60):
+ b = sqrt(5<<(2*bi-1))
+ for m in range(20, 600, 32):
+ m = sqrt(2**(m+1))
+ print "modmul", hexstr(a), hexstr(b), hexstr(m), hexstr((a*b) % m)
+
# Simple tests of modpow.
for i in range(64, 4097, 63):
modulus = sqrt(1<<(2*i-1)) | 1
# Test even moduli, which can't be done by Montgomery.
modulus = modulus - 1
print "pow", hexstr(base), hexstr(expt), hexstr(modulus), hexstr(pow(base, expt, modulus))
+ print "pow", hexstr(i), hexstr(expt), hexstr(modulus), hexstr(pow(i, expt, modulus))