+ print "mul", hexstr(a), hexstr(b), hexstr(p)
+
+# Simple tests of modpow.
+for i in range(64, 4097, 63):
+ modulus = sqrt(1<<(2*i-1)) | 1
+ base = sqrt(3*modulus*modulus) % modulus
+ expt = sqrt(modulus*modulus*2/5)
+ print "pow", hexstr(base), hexstr(expt), hexstr(modulus), hexstr(pow(base, expt, modulus))
+ if i <= 1024:
+ # 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))