git.distorted.org.uk Git - u/mdw/putty/blob - testdata/bignum.py

   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)