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 modpow.
107 for i
in range(64, 4097, 63):
108 modulus
= sqrt(1<<(2*i
-1)) |
1
109 base
= sqrt(3*modulus
*modulus
) % modulus
110 expt
= sqrt(modulus
*modulus
*2/5)
111 print "pow", hexstr(base
), hexstr(expt
), hexstr(modulus
), hexstr(pow(base
, expt
, modulus
))
113 # Test even moduli, which can't be done by Montgomery.
114 modulus
= modulus
- 1
115 print "pow", hexstr(base
), hexstr(expt
), hexstr(modulus
), hexstr(pow(base
, expt
, modulus
))