+++ /dev/null
-# Generate test cases for a bignum implementation.\r
-\r
-import sys\r
-\r
-# integer square roots\r
-def sqrt(n):\r
- d = long(n)\r
- a = 0L\r
- # b must start off as a power of 4 at least as large as n\r
- ndigits = len(hex(long(n)))\r
- b = 1L << (ndigits*4)\r
- while 1:\r
- a = a >> 1\r
- di = 2*a + b\r
- if di <= d:\r
- d = d - di\r
- a = a + b\r
- b = b >> 2\r
- if b == 0: break\r
- return a\r
-\r
-# continued fraction convergents of a rational\r
-def confrac(n, d):\r
- coeffs = [(1,0),(0,1)]\r
- while d != 0:\r
- i = n / d\r
- n, d = d, n % d\r
- coeffs.append((coeffs[-2][0]-i*coeffs[-1][0],\r
- coeffs[-2][1]-i*coeffs[-1][1]))\r
- return coeffs\r
-\r
-def findprod(target, dir = +1, ratio=(1,1)):\r
- # Return two numbers whose product is as close as we can get to\r
- # 'target', with any deviation having the sign of 'dir', and in\r
- # the same approximate ratio as 'ratio'.\r
-\r
- r = sqrt(target * ratio[0] * ratio[1])\r
- a = r / ratio[1]\r
- b = r / ratio[0]\r
- if a*b * dir < target * dir:\r
- a = a + 1\r
- b = b + 1\r
- assert a*b * dir >= target * dir\r
-\r
- best = (a,b,a*b)\r
-\r
- while 1:\r
- improved = 0\r
- a, b = best[:2]\r
-\r
- coeffs = confrac(a, b)\r
- for c in coeffs:\r
- # a*c[0]+b*c[1] is as close as we can get it to zero. So\r
- # if we replace a and b with a+c[1] and b+c[0], then that\r
- # will be added to our product, along with c[0]*c[1].\r
- da, db = c[1], c[0]\r
-\r
- # Flip signs as appropriate.\r
- if (a+da) * (b+db) * dir < target * dir:\r
- da, db = -da, -db\r
-\r
- # Multiply up. We want to get as close as we can to a\r
- # solution of the quadratic equation in n\r
- #\r
- # (a + n da) (b + n db) = target\r
- # => n^2 da db + n (b da + a db) + (a b - target) = 0\r
- A,B,C = da*db, b*da+a*db, a*b-target\r
- discrim = B^2-4*A*C\r
- if discrim > 0 and A != 0:\r
- root = sqrt(discrim)\r
- vals = []\r
- vals.append((-B + root) / (2*A))\r
- vals.append((-B - root) / (2*A))\r
- if root * root != discrim:\r
- root = root + 1\r
- vals.append((-B + root) / (2*A))\r
- vals.append((-B - root) / (2*A))\r
-\r
- for n in vals:\r
- ap = a + da*n\r
- bp = b + db*n\r
- pp = ap*bp\r
- if pp * dir >= target * dir and pp * dir < best[2]*dir:\r
- best = (ap, bp, pp)\r
- improved = 1\r
-\r
- if not improved:\r
- break\r
-\r
- return best\r
-\r
-def hexstr(n):\r
- s = hex(n)\r
- if s[:2] == "0x": s = s[2:]\r
- if s[-1:] == "L": s = s[:-1]\r
- return s\r
-\r
-# Tests of multiplication which exercise the propagation of the last\r
-# carry to the very top of the number.\r
-for i in range(1,4200):\r
- a, b, p = findprod((1<<i)+1, +1, (i, i*i+1))\r
- print "mul", hexstr(a), hexstr(b), hexstr(p)\r
- a, b, p = findprod((1<<i)+1, +1, (i, i+1))\r
- print "mul", hexstr(a), hexstr(b), hexstr(p)\r
-\r
-# Simple tests of modpow.\r
-for i in range(64, 4097, 63):\r
- modulus = sqrt(1<<(2*i-1)) | 1\r
- base = sqrt(3*modulus*modulus) % modulus\r
- expt = sqrt(modulus*modulus*2/5)\r
- print "pow", hexstr(base), hexstr(expt), hexstr(modulus), hexstr(pow(base, expt, modulus))\r
- if i <= 1024:\r
- # Test even moduli, which can't be done by Montgomery.\r
- modulus = modulus - 1\r
- print "pow", hexstr(base), hexstr(expt), hexstr(modulus), hexstr(pow(base, expt, modulus))\r
OSDN Git Service
