## File: bignum.py

package info (click to toggle)
putty 0.62-9+deb7u3
 `123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125` ``````# Generate test cases for a bignum implementation. import sys # integer square roots def sqrt(n): d = long(n) a = 0L # b must start off as a power of 4 at least as large as n ndigits = len(hex(long(n))) b = 1L << (ndigits*4) while 1: a = a >> 1 di = 2*a + b if di <= d: d = d - di a = a + b b = b >> 2 if b == 0: break return a # continued fraction convergents of a rational def confrac(n, d): coeffs = [(1,0),(0,1)] while d != 0: i = n / d n, d = d, n % d coeffs.append((coeffs[-2]-i*coeffs[-1], coeffs[-2]-i*coeffs[-1])) return coeffs def findprod(target, dir = +1, ratio=(1,1)): # Return two numbers whose product is as close as we can get to # 'target', with any deviation having the sign of 'dir', and in # the same approximate ratio as 'ratio'. r = sqrt(target * ratio * ratio) a = r / ratio b = r / ratio if a*b * dir < target * dir: a = a + 1 b = b + 1 assert a*b * dir >= target * dir best = (a,b,a*b) while 1: improved = 0 a, b = best[:2] coeffs = confrac(a, b) for c in coeffs: # a*c+b*c is as close as we can get it to zero. So # if we replace a and b with a+c and b+c, then that # will be added to our product, along with c*c. da, db = c, c # Flip signs as appropriate. if (a+da) * (b+db) * dir < target * dir: da, db = -da, -db # Multiply up. We want to get as close as we can to a # solution of the quadratic equation in n # # (a + n da) (b + n db) = target # => n^2 da db + n (b da + a db) + (a b - target) = 0 A,B,C = da*db, b*da+a*db, a*b-target discrim = B^2-4*A*C if discrim > 0 and A != 0: root = sqrt(discrim) vals = [] vals.append((-B + root) / (2*A)) vals.append((-B - root) / (2*A)) if root * root != discrim: root = root + 1 vals.append((-B + root) / (2*A)) vals.append((-B - root) / (2*A)) for n in vals: ap = a + da*n bp = b + db*n pp = ap*bp if pp * dir >= target * dir and pp * dir < best*dir: best = (ap, bp, pp) improved = 1 if not improved: break return best def hexstr(n): s = hex(n) if s[:2] == "0x": s = s[2:] if s[-1:] == "L": s = s[:-1] return s # Tests of multiplication which exercise the propagation of the last # carry to the very top of the number. for i in range(1,4200): a, b, p = findprod((1<