--- /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
