f3c29e34 |
1 | # Generate test cases for a bignum implementation. |
2 | |
3 | import sys |
4 | import mathlib |
5 | |
6 | def findprod(target, dir = +1, ratio=(1,1)): |
7 | # Return two numbers whose product is as close as we can get to |
8 | # 'target', with any deviation having the sign of 'dir', and in |
9 | # the same approximate ratio as 'ratio'. |
10 | |
11 | r = mathlib.sqrt(target * ratio[0] * ratio[1]) |
12 | a = r / ratio[1] |
13 | b = r / ratio[0] |
14 | if a*b * dir < target * dir: |
15 | a = a + 1 |
16 | b = b + 1 |
17 | assert a*b * dir >= target * dir |
18 | |
19 | best = (a,b,a*b) |
20 | |
21 | while 1: |
22 | improved = 0 |
23 | a, b = best[:2] |
24 | |
25 | terms = mathlib.confracr(a, b, output=None) |
26 | coeffs = [(1,0),(0,1)] |
27 | for t in terms: |
28 | coeffs.append((coeffs[-2][0]-t*coeffs[-1][0], |
29 | coeffs[-2][1]-t*coeffs[-1][1])) |
30 | for c in coeffs: |
31 | # a*c[0]+b*c[1] is as close as we can get it to zero. So |
32 | # if we replace a and b with a+c[1] and b+c[0], then that |
33 | # will be added to our product, along with c[0]*c[1]. |
34 | da, db = c[1], c[0] |
35 | |
36 | # Flip signs as appropriate. |
37 | if (a+da) * (b+db) * dir < target * dir: |
38 | da, db = -da, -db |
39 | |
40 | # Multiply up. We want to get as close as we can to a |
41 | # solution of the quadratic equation in n |
42 | # |
43 | # (a + n da) (b + n db) = target |
44 | # => n^2 da db + n (b da + a db) + (a b - target) = 0 |
45 | A,B,C = da*db, b*da+a*db, a*b-target |
46 | discrim = B^2-4*A*C |
47 | if discrim > 0 and A != 0: |
48 | root = mathlib.sqrt(discrim) |
49 | vals = [] |
50 | vals.append((-B + root) / (2*A)) |
51 | vals.append((-B - root) / (2*A)) |
52 | if root * root != discrim: |
53 | root = root + 1 |
54 | vals.append((-B + root) / (2*A)) |
55 | vals.append((-B - root) / (2*A)) |
56 | |
57 | for n in vals: |
58 | ap = a + da*n |
59 | bp = b + db*n |
60 | pp = ap*bp |
61 | if pp * dir >= target * dir and pp * dir < best[2]*dir: |
62 | best = (ap, bp, pp) |
63 | improved = 1 |
64 | |
65 | if not improved: |
66 | break |
67 | |
68 | return best |
69 | |
70 | def hexstr(n): |
71 | s = hex(n) |
72 | if s[:2] == "0x": s = s[2:] |
73 | if s[-1:] == "L": s = s[:-1] |
74 | return s |
75 | |
76 | # Tests of multiplication which exercise the propagation of the last |
77 | # carry to the very top of the number. |
78 | for i in range(1,4200): |
79 | a, b, p = findprod((1<<i)+1, +1, (i, i*i+1)) |
80 | print hexstr(a), hexstr(b), hexstr(p) |
81 | a, b, p = findprod((1<<i)+1, +1, (i, i+1)) |
82 | print hexstr(a), hexstr(b), hexstr(p) |