1 # Generate test cases for a bignum implementation.
9 # b must start off as a power of 4 at least as large as n
10 ndigits = len(hex(long(n)))
22 # continued fraction convergents of a rational
24 coeffs = [(1,0),(0,1)]
28 coeffs.append((coeffs[-2][0]-i*coeffs[-1][0],
29 coeffs[-2][1]-i*coeffs[-1][1]))
32 def findprod(target, dir = +1, ratio=(1,1)):
33 # Return two numbers whose product is as close as we can get to
34 # 'target', with any deviation having the sign of 'dir', and in
35 # the same approximate ratio as 'ratio'.
37 r = sqrt(target * ratio[0] * ratio[1])
40 if a*b * dir < target * dir:
43 assert a*b * dir >= target * dir
51 coeffs = confrac(a, b)
53 # a*c[0]+b*c[1] is as close as we can get it to zero. So
54 # if we replace a and b with a+c[1] and b+c[0], then that
55 # will be added to our product, along with c[0]*c[1].
58 # Flip signs as appropriate.
59 if (a+da) * (b+db) * dir < target * dir:
62 # Multiply up. We want to get as close as we can to a
63 # solution of the quadratic equation in n
65 # (a + n da) (b + n db) = target
66 # => n^2 da db + n (b da + a db) + (a b - target) = 0
67 A,B,C = da*db, b*da+a*db, a*b-target
69 if discrim > 0 and A != 0:
72 vals.append((-B + root) / (2*A))
73 vals.append((-B - root) / (2*A))
74 if root * root != discrim:
76 vals.append((-B + root) / (2*A))
77 vals.append((-B - root) / (2*A))
83 if pp * dir >= target * dir and pp * dir < best[2]*dir:
94 if s[:2] == "0x": s = s[2:]
95 if s[-1:] == "L": s = s[:-1]
98 # Tests of multiplication which exercise the propagation of the last
99 # carry to the very top of the number.
100 for i in range(1,4200):
101 a, b, p = findprod((1<<i)+1, +1, (i, i*i+1))
102 print "mul", hexstr(a), hexstr(b), hexstr(p)
103 a, b, p = findprod((1<<i)+1, +1, (i, i+1))
104 print "mul", hexstr(a), hexstr(b), hexstr(p)
106 # Simple tests of modmul.
107 for ai in range(20, 200, 60):
108 a = sqrt(3<<(2*ai-1))
109 for bi in range(20, 200, 60):
110 b = sqrt(5<<(2*bi-1))
111 for m in range(20, 600, 32):
113 print "modmul", hexstr(a), hexstr(b), hexstr(m), hexstr((a*b) % m)
115 # Simple tests of modpow.
116 for i in range(64, 4097, 63):
117 modulus = sqrt(1<<(2*i-1)) | 1
118 base = sqrt(3*modulus*modulus) % modulus
119 expt = sqrt(modulus*modulus*2/5)
120 print "pow", hexstr(base), hexstr(expt), hexstr(modulus), hexstr(pow(base, expt, modulus))
122 # Test even moduli, which can't be done by Montgomery.
123 modulus = modulus - 1
124 print "pow", hexstr(base), hexstr(expt), hexstr(modulus), hexstr(pow(base, expt, modulus))
125 print "pow", hexstr(i), hexstr(expt), hexstr(modulus), hexstr(pow(i, expt, modulus))