开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > 加密解密 > CA认证 > 查看源码
bntest16.c
资源名称：certlib.tar.gz　[点击查看]
上传用户：zbbssh
上传日期：2007-01-08
资源大小：196k
文件大小：21k
源码类别：

CA认证

开发平台：

C/C++

bntest16.c：源码内容
  1. /*
  2.  * Test driver for low-level bignum library (16-bit version).
  3.  * This access the low-level library directly.  It is NOT an
  4.  * example of how to program with the library normally!  By
  5.  * accessing the library at a low level, it is possible to
  6.  * exercise the smallest components and thus localize bugs
  7.  * more accurately.  This is especially useful when writing
  8.  * assembly-language primitives.
  9.  *
  10.  * This also does timing tests on modular exponentiation.
  11.  * Modular exponentiation is so computationally expensive that
  12.  * the fact that this code omits one level of interface glue
  13.  * has no perceptible effect on the results.
  14.  */
  15. #if HAVE_CONFIG_H
  16. #include "config.h"
  17. #endif
  19. #include <stdio.h>
  20. #if !NO_STDLIB_H
  21. #include <stdlib.h> /* For strtol */
  22. #endif
  23. #if !NO_STRING_H
  24. #include <string.h>
  25. #elif HAVE_STRINGS_H
  26. #include <strings.h>
  27. #endif
  28. #if NEED_MEMORY_H
  29. #include <memory.h>
  30. #endif
  32. #include "cputime.h"
  33. #include "lbn16.h"
  35. #include "kludge.h"
  37. /* Work with up to 2048-bit numbers */
  38. #define MAXBITS 2048
  39. #define SIZE (MAXBITS/16 + 1)
  41. /* Additive congruential random number generator, x[i] = x[i-24] + x[i-55] */
  42. static BNWORD16 randp[55];
  43. static BNWORD16 *randp1 = randp, *randp2 = randp+24;
  45. static BNWORD16
  46. rand16(void)
  47. {
  48.         if (++randp2 == randp+55) {
  49.                 randp2 = randp;
  50.                 randp1++;
  51.         } else if (++randp1 == randp+55) {
  52.                 randp1 = randp;
  53.         }
  55.         return  *randp1 += *randp2;
  56. }
  58. /*
  59.  * CRC-3_2: x^3_2+x^26+x^23+x^22+x^1_6+x^12+x^11+x^10+x^8+x^7+x^5+x^4+x^2+x+1
  60.  *
  61.  * The additive congruential RNG is seeded with a single integer,
  62.  * which is shuffled with a CRC polynomial to generate the initial
  63.  * table values.  The Polynomial is the same size as the words being
  64.  * used.
  65.  *
  66.  * Thus, in the various versions of this library, we actually use this
  67.  * polynomial as-is, this polynomial mod x^17, and this polynomial with
  68.  * the leading coefficient deleted and replaced with x^6_4.  As-is,
  69.  * it's irreducible, so it has a long period.  Modulo x^17, it factors as
  70.  * (x^4+x^3+x^2+x+1) * (x^12+x^11+x^8+x^7+x^6+x^5+x^4+x^3+1),
  71.  * which still has a large enough period (4095) for the use it's put to.
  72.  * With the leading coefficient moved up, it factors as
  73.  * (x^50+x^49+x^48+x^47+x^46+x^43+x^41+x^40+x^38+x^37+x^36+x^35+x^34+x^33+
  74.  *  x^31+x^30+x^29+x^28+x^27+x^25+x^23+x^18+x^1_6+x^15+x^14+x^13+x^11+x^9+
  75.  *  x^8+x^7+x^6+x^5+x^3+x^2+1)*(x^11+x^10+x^9+x^5+x^4+x^3+1)*(x^3+x+1),
  76.  * which definitely has a long enough period to serve for initialization.
  77.  * 
  78.  * The effort put into this PRNG is kind of unwarranted given the trivial
  79.  * use it's being put to, but oh, well.  It does have the nice advantage
  80.  * of producing numbers that are portable between platforms, so if there's
  81.  * a problem with one platform, you can compare all the intermediate
  82.  * results with another platform.
  83.  */
  84. #define POLY (BNWORD16)0x04c11db7
  86. static void
  87. srand16(BNWORD16 seed)
  88. {
  89.         int i, j;
  91.         for (i = 0; i < 55; i++) {
  92.                 for (j = 0; j < 16; j++)
  93.                         if (seed >> (16-1))
  94.                                 seed = (seed << 1) ^ POLY;
  95.                         else
  96.                                 seed <<= 1;
  97.                 randp[i] = seed;
  98.         }
  99.         for (i = 0; i < 3*55; i ++)
  100.                 rand16();
  101. }
  103. static void
  104. randnum(BNWORD16 *num, unsigned len)
  105. {
  106.         while (len--)
  107.                 BIGLITTLE(*--num,*num++) = rand16();
  108. }
  110. static void
  111. bnprint16(BNWORD16 const *num, unsigned len)
  112. {
  113.         BIGLITTLE(num -= len, num += len);
  115.         while (len--)
  116.                 printf("%0*lX", 16/4, (unsigned long)BIGLITTLE(*num++,*--num));
  117. }
  119. static void
  120. bnput16(char const *prompt, BNWORD16 const *num, unsigned len)
  121. {
  122.         fputs(prompt, stdout);
  123.         bnprint16(num, len);
  124.         putchar('n');
  125. }
  127. /*
  128.  * One of our tests uses a known prime.  The following selections were
  129.  * taken from the tables at the end of Hans Reisel's "Prime Numbers and
  130.  * Computer Methods for Factorization", second edition - an excellent book.
  131.  * (ISBN 0-8176-3743-5 ISBN 3-7643-3743-5)
  132.  */
  133. #if 0
  134. /* P31=1839605 17620282 38179967 87333633 from the factors of 3^256+2^256 */
  135. static unsigned char const prime[] = {
  136.         0x17,0x38,0x15,0xBC,0x8B,0xBB,0xE9,0xEF,0x01,0xA9,0xFD,0x3A,0x01
  137. };
  138. #elif 0
  139. /* P48=40554942 04557502 46193993 36199835 4279613_2 73199617 from the same */
  140. static unsigned char const prime[] = {
  141.         0x47,0x09,0x77,0x07,0xCF,0xFD,0xE1,0x54,0x3E,0x24,
  142.         0xF7,0xF1,0x7A,0x3E,0x91,0x51,0xCC,0xC7,0xD4,0x01
  143. };
  144. #elif 0
  145. /*
  146.  * P75 = 450 55287640 97906895 47687014 5808213_2
  147.  *  05219565 99525911 39967964 66003_258 91979521
  148.  * from the factors of 4^128+3+128
  149.  * (The "026" and "062" are to prevent a Bad String from appearing here.)
  150.  */
  151. static unsigned char const prime[] = {
  152. 0xFF,0x00,0xFF,0x00,0xFF,0x01,0x06,0x4F,0xF8,0xED,
  153. 0xA3,0x37,0x23,0x2A,0x04,0xEA,0xF9,0x5F,0x30,0x4C,
  154. 0xAE,0xCD, 026,0x4E, 062,0x10,0x04,0x7D,0x0D,0x79,
  155. 0x01
  156. };
  157. #else
  158. /*
  159.  * P75 = 664 85659796 45277755 9123_2190 67300940
  160.  *  51844953 78793489 59444670 35675855 57440257
  161.  * from the factors of 5^128+4^128
  162.  * (The "026" is to prevent a Bad String from appearing here.)
  163.  */
  164. static unsigned char const prime[] = {
  165. 0x01,0x78,0x4B,0xA5,0xD3,0x30,0x03,0xEB,0x73,0xE6,
  166. 0x0F,0x4E,0x31,0x7D,0xBC,0xE2,0xA0,0xD4, 026,0x3F,
  167. 0x3C,0xEA,0x1B,0x44,0xAD,0x39,0xE7,0xE5,0xAD,0x19,
  168. 0x67,0x01
  169. };
  170. #endif
  172. static int
  173. usage(char const *name)
  174. {
  175.         fprintf(stderr, "Usage: %s [modbits [expbits [expbits2]]n"
  176. "With no arguments, just runs test suite.  If modbits is given, runsn"
  177. "quick validation test, then runs timing tests of modular exponentiation.n"
  178. "If expbits is given, it is used as an exponent size, otherwise it defaultsn"
  179. "to the same as modbits.  If expbits2 is given it is used as the secondn"
  180. "exponent size in the double-exponentiation tests, otherwise it defaultsn"
  181. "to the same as expbits.  All are limited to %u bits.n",
  182.                 name, (unsigned)MAXBITS);
  183.         return 1;
  184. }
  186. int
  187. main(int argc, char **argv)
  188. {
  189. unsigned i, j, k, l, m;
  190. int z;
  191. BNWORD16 t, carry, borrow;
  192. BNWORD16 a[SIZE], b[SIZE], c[SIZE], d[SIZE];
  193. BNWORD16 e[SIZE], f[SIZE];
  194. long modbits = 0, expbits = 0, expbits2 = 0;
  195. char *p;
  196. #define A BIGLITTLE((a+SIZE),a)
  197. #define B BIGLITTLE((b+SIZE),b)
  198. #define C BIGLITTLE((c+SIZE),c)
  199. #define D BIGLITTLE((d+SIZE),d)
  200. #define E BIGLITTLE((e+SIZE),e)
  201. #define F BIGLITTLE((f+SIZE),f)
  202.         static unsigned const smallprimes[] = {
  203.                 2, 3, 5, 7, 11, 13, 17, 19, 23, 27, 29, 31, 37, 41, 43
  204.         };
  206.         srand16(1);
  208.         puts(BIGLITTLE("Big-endian machine","Little-endian machine"));
  210. /*
  211.  * Parse the command line.  The do{} loop is a dummy, to have
  212.  * something to break out of when done with the command line.
  213.  */
  214. do {
  215. if (argc <= 1)
  216. break;
  217.                 modbits = strtol(argv[1], &p, 0);
  218.                 if (modbits < 3 || modbits > MAXBITS || *p) {
  219.                         fprintf(stderr, "Illegal modbits: "%s"n", argv[1]);
  220.                         return usage(argv[0]);
  221.                 }
  222. if (argc <= 2) {
  223. expbits2 = expbits = modbits;
  224. break;
  225. }
  226. expbits = strtol(argv[2], &p, 0);
  227. if (expbits < 3 || expbits > MAXBITS || *p) {
  228. fprintf(stderr, "Illegal expbits: "%s"n", argv[2]);
  229. return usage(argv[0]);
  230. }
  231. if (argc <= 3) {
  232. expbits2 = expbits;
  233. break;
  234. }
  235. expbits2 = strtol(argv[3], &p, 0);
  236. if (expbits2 < 3 || expbits2 > MAXBITS || *p) {
  237. fprintf(stderr, "Illegal expbits2: "%s"n", argv[3]);
  238. return usage(argv[0]);
  239. }
  240. if (argc <= 4)
  241. break;
  242. fprintf(stderr, "Unexpected fourth argument "%s"n",
  243.         argv[4]);
  244. return usage(argv[0]);
  245. } while (0);
  247.         /* B is a nice not-so-little prime */
  248.         lbnInsertBigBytes_16(B, prime, 0, sizeof(prime));
  249.         ((unsigned char *)c)[0] = 0;
  250.         lbnInsertBigBytes_16(B, (unsigned char *)c, sizeof(prime), 1);
  251.         lbnExtractBigBytes_16(B, (unsigned char *)c, 0, sizeof(prime)+1);
  252.         i = (sizeof(prime)-1)/(16/8)+1;        /* Size of array in words */
  253.         if (((unsigned char *)c)[0] ||
  254.             memcmp(prime, (unsigned char *)c+1, sizeof(prime)) != 0)
  255.         {
  256.                 printf("Input != output!:n   ");
  257.                 for (k = 0; k < sizeof(prime); k++)
  258.                         printf("%02X ", prime[k]);
  259.                 putchar('n');
  260.                 for (k = 0; k < sizeof(prime)+1; k++)
  261.                         printf("%02X ", ((unsigned char *)c)[k]);
  262.                 putchar('n');
  263.                 bnput16("p = ", B, i);
  265.         }
  267.         /* Timing test code - only if requested on the command line */
  268. if (modbits) {
  269. #if CLOCK_AVAIL
  270. timetype start, stop;
  271. unsigned long cursec, expsec, twoexpsec, dblexpsec;
  272. unsigned curms, expms, twoexpms, dblexpms;
  274. expsec = twoexpsec = dblexpsec = 0;
  275. expms = twoexpms = dblexpms = 0;
  276. #endif
  278. lbnCopy_16(C,B,i);
  279. lbnSub1_16(C,i,1);        /* C is exponent: p-1 */
  281. puts("Testing modexp with a known prime.  "
  282.      "All results should be 1.");
  283. bnput16("p   = ", B, i);
  284. bnput16("p-1 = ", C, i);
  285. lbnTwoExpMod_16(A, C, i, B, i);
  286. bnput16("2^(p-1) mod p = ", A, i);
  287. for (j = 0; j < 10; j++) {
  288. randnum(A,i);
  289. (void)lbnDiv_16(D,A,i,B,i);
  291. bnput16("a = ", A, i);
  292. lbnExpMod_16(A, A, i, C, i, B, i);
  293. bnput16("a^(p-1) mod p = ", A, i);
  294. }
  296. printf("n"
  297.        "Timing exponentiations modulo a %d-bit modulus, i.e.n"
  298.        "2^<%d> mod <%d> bits, <%d>^<%d> mod <%d> bits andn"
  299.        "<%d>^<%d> * <%d>^<%d> mod <%d> bitsn",
  300.        (int)modbits, (int)expbits, (int)modbits,
  301.        (int)modbits, (int)expbits, (int)modbits,
  302.        (int)modbits, (int)expbits, (int)modbits, (int)expbits2,
  303.        (int)modbits);
  305. i = ((int)modbits-1)/16+1;
  306. k = ((int)expbits-1)/16+1;
  307. l = ((int)expbits2-1)/16+1;
  308. for (j = 0; j < 25; j++) {
  309. randnum(A,i);        /* Base */
  310. randnum(B,k);        /* Exponent */
  311. randnum(C,i);        /* Modulus */
  312. randnum(D,i);        /* Base2 */
  313. randnum(E,l);        /* Exponent */
  314. /* Clip bases and mod to appropriate number of bits */
  315. t = ((BNWORD16)2<<((modbits-1)%16)) - 1;
  316. *(BIGLITTLE(A-i,A+i-1)) &= t;
  317. *(BIGLITTLE(C-i,C+i-1)) &= t;
  318. *(BIGLITTLE(D-i,D+i-1)) &= t;
  319. /* Make modulus large (msbit set) and odd (lsbit set) */
  320. *(BIGLITTLE(C-i,C+i-1)) |= (t >> 1) + 1;
  321. BIGLITTLE(C[-1],C[0]) |= 1;
  323. /* Clip exponent to appropriate number of bits */
  324. t = ((BNWORD16)2<<((expbits-1)%16)) - 1;
  325. *(BIGLITTLE(B-k,B+k-1)) &= t;
  326. /* Make exponent large (msbit set) */
  327. *(BIGLITTLE(B-k,B+k-1)) |= (t >> 1) + 1;
  328. /* The same for exponent 2 */
  329. t = ((BNWORD16)2<<((expbits2-1)%16)) - 1;
  330. *(BIGLITTLE(E-l,E+l-1)) &= t;
  331. *(BIGLITTLE(E-l,E+l-1)) |= (t >> 1) + 1;
  333. m = lbnBits_16(A, i);
  334. if (m > (unsigned)modbits) {
  335. bnput16("a = ", a, i);
  336. printf("%u bits, should be <= %dn",
  337.        m, (int)modbits);
  338. }
  339. m = lbnBits_16(B, k);
  340. if (m != (unsigned)expbits) {
  341. bnput16("b = ", b, i);
  342. printf("%u bits, should be %dn",
  343.        m, (int)expbits);
  344. }
  345. m = lbnBits_16(C, i);
  346. if (m != (unsigned)modbits) {
  347. bnput16("c = ", c, k);
  348. printf("%u bits, should be %dn",
  349.        m, (int)modbits);
  350. }
  351. m = lbnBits_16(D, i);
  352. if (m > (unsigned)modbits) {
  353. bnput16("d = ", d, i);
  354. printf("%u bits, should be <= %dn",
  355.        m, (int)modbits);
  356. }
  357. m = lbnBits_16(E, l);
  358. if (m != (unsigned)expbits2) {
  359. bnput16("e = ", e, i);
  360. printf("%u bits, should be %dn",
  361.        m, (int)expbits2);
  362. }
  363. #if CLOCK_AVAIL
  364. gettime(&start);
  365. #endif
  366. lbnTwoExpMod_16(A, B, k, C, i);
  367. #if CLOCK_AVAIL
  368. gettime(&stop);
  370. subtime(stop, start);
  371. twoexpsec += cursec = sec(stop);
  372. twoexpms += curms = msec(stop);
  374. printf("2^<%d>:%3lu.%03u   ",
  375.        (int)expbits, cursec, curms);
  376. #else
  377. printf("2^<%d>    ", (int)expbits);
  378. #endif
  379. fflush(stdout);
  381. #if CLOCK_AVAIL
  382. gettime(&start);
  383. #endif
  384. lbnExpMod_16(A, A, i, B, k, C, i);
  385. #if CLOCK_AVAIL
  386. gettime(&stop);
  388. subtime(stop, start);
  389. expsec += cursec = sec(stop);
  390. expms += curms = msec(stop);
  391. printf("<%d>^<%d>:%3lu.%03u   ",
  392.        (int)modbits, (int)expbits, cursec, curms);
  393. #else
  394. printf("<%d>^<%d>    ", (int)modbits, (int)expbits);
  395. #endif
  396. fflush(stdout);
  398. #if CLOCK_AVAIL
  399. gettime(&start);
  400. #endif
  401. lbnDoubleExpMod_16(D, A, i, B, k, D, i, E, l, C, i);
  402. #if CLOCK_AVAIL
  403. gettime(&stop);
  405. subtime(stop, start);
  406. dblexpsec += cursec = sec(stop);
  407. dblexpms += curms = msec(stop);
  408. printf("<%d>^<%d>*<%d>^<%d>:%3lu.%03un",
  409.        (int)modbits, (int)expbits,
  410.        (int)modbits, (int)expbits2,
  411.        cursec, curms);
  412. #else
  413. printf("<%d>^<%d>*<%d>^<%d>n",
  414.        (int)modbits, (int)expbits,
  415.        (int)modbits, (int)expbits2);
  416. #endif
  417. }
  418. #if CLOCK_AVAIL
  419. putchar('n');
  420. twoexpms += 1000 * (unsigned)(twoexpsec % j);
  421. while (twoexpms > 1000*j) {
  422. twoexpsec += j;
  423. twoexpms -= 1000*j;
  424. }
  425. printf("2^<%d> mod <%d> bits AVERAGE: %lu.%03u sn",
  426.        (int)expbits, (int)modbits,
  427.        twoexpsec/j, twoexpms/j);
  429. expms += 1000 * (unsigned)(expsec % j);
  430. while (expms > 1000*j) {
  431. expsec += j;
  432. expms -= 1000*j;
  433. }
  434. printf("<%d>^<%d> mod <%d> bits AVERAGE: %lu.%03u sn",
  435.        (int)modbits, (int)expbits, (int)modbits,
  436.        expsec/j, expms/j);
  438. dblexpms += 1000 * (unsigned)(dblexpsec % j);
  439. while (dblexpms > 1000*j) {
  440. dblexpsec += j;
  441. dblexpms -= 1000*j;
  442. }
  443. printf("<%d>^<%d> * <%d>^<%d> mod <%d> bits AVERAGE:"
  444.        " %lu.%03u sn",
  445.        (int)modbits, (int)expbits, (int)modbits, (int)expbits2,
  446.        (int)modbits, dblexpsec/j, dblexpms/j);
  447. #endif
  448. putchar('n');
  449. }
  451.         puts("Beginning 1000 interations of sanity checking.n"
  452.              "Any output indicates a bug.  No output is very strongn"
  453.              "evidence that all the important low-level bignum routinesn"
  454.              "are working properly.n");
  456. /*
  457.  * If you change this loop to have an iteration 0, all results
  458.  * are primted on that iteration.  Useful to see what's going
  459.  * on in case of major wierdness, but it produces a *lot* of
  460.  * output.
  461.  */
  462.         for (j = 1; j <= 1000; j++) {
  463. /* Dop the tests for lots of different number sizes. */
  464. for (i = 1; i <= SIZE/2; i++) {
  465. /* Make a random number i words long */
  466. do {
  467. randnum(A,i);
  468. } while (lbnNorm_16(A,i) < i);
  470. /* Checl lbnCmp - does a == a? */
  471. if (lbnCmp_16(A,A,i) || !j) {
  472. bnput16("a = ", A, i);
  473. printf("(a <=> a) = %dn", lbnCmp_16(A,A,i));
  474. }
  476. memcpy(c, a, sizeof(a));
  478. /* Check that the difference, after copy, is good. */
  479. if (lbnCmp_16(A,C,i) || !j) {
  480. bnput16("a = ", A, i);
  481. bnput16("c = ", C, i);
  482. printf("(a <=> c) = %dn", lbnCmp_16(A,C,i));
  483. }
  485. /* Generate a non-zero random t */
  486. do {
  487. t = rand16();
  488. } while (!t);
  490. /*
  491.  * Add t to A.  Check that:
  492.  * - lbnCmp works in both directions, and
  493.  * - A + t is greater than A.  If there was a carry,
  494.  *   the result, less the carry, should be *less*
  495.  *   than A.
  496.  */
  497. carry = lbnAdd1_16(A,i,t);
  498. if (lbnCmp_16(A,C,i) + lbnCmp_16(C,A,i) != 0 ||
  499.     lbnCmp_16(A,C,i) != (carry ? -1 : 1) || !j)
  500. {
  501. bnput16("c       = ", C, i);
  502. printf("t = %lXn", (unsigned long)t);
  503. bnput16("a = c+t = ", A, i);
  504. printf("carry = %lXn", (unsigned long)carry);
  505. printf("(a <=> c) = %dn", lbnCmp_16(A,C,i));
  506. printf("(c <=> a) = %dn", lbnCmp_16(C,A,i));
  507. }
  509. /* Subtract t again */
  510. memcpy(d, a, sizeof(a));
  511. borrow = lbnSub1_16(A,i,t);
  513. if (carry != borrow || lbnCmp_16(A,C,i) || !j) {
  514. bnput16("a = ", C, i);
  515. printf("t = %lXn", (unsigned long)t);
  516. lbnAdd1_16(A,i,t);
  517. bnput16("a += t = ", A, i);
  518. printf("Carry = %lXn", (unsigned long)carry);
  519. lbnSub1_16(A,i,t);
  520. bnput16("a -= t = ", A, i);
  521. printf("Borrow = %lXn", (unsigned long)borrow);
  522. printf("(a <=> c) = %dn", lbnCmp_16(A,C,i));
  523. }
  525. /* Generate a random B */
  526. do {
  527. randnum(B,i);
  528. } while (lbnNorm_16(B,i) < i);
  530. carry = lbnAddN_16(A,B,i);
  531. memcpy(d, a, sizeof(a));
  532. borrow = lbnSubN_16(A,B,i);
  534. if (carry != borrow || lbnCmp_16(A,C,i) || !j) {
  535. bnput16("a = ", C, i);
  536. bnput16("b = ", B, i);
  537. bnput16("a += b = ", D, i);
  538. printf("Carry = %lXn", (unsigned long)carry);
  539. bnput16("a -= b = ", A, i);
  540. printf("Borrow = %lXn", (unsigned long)borrow);
  541. printf("(a <=> c) = %dn", lbnCmp_16(A,C,i));
  542. }
  544. /* D = B * t */
  545. lbnMulN1_16(D, B, i, t);
  546. memcpy(e, d, sizeof(a));
  547. carry = *(BIGLITTLE(D-i-1,D+i));
  548. /* D = A + B * t */
  549. carry = lbnAddN_16(D,A,i);
  551. /* For comparison with the result of lbnMulAdd1_16 */
  552. borrow = carry + *(BIGLITTLE(D-i-1,D+i));
  554. *(BIGLITTLE(A-i-1,A+i)) = 0;
  555. carry = lbnMulAdd1_16(A, B, i, t);
  557. /* Did MulAdd get the same answer as mul then add? */
  558. if (carry != borrow || lbnCmp_16(A, D, i) || !j) {
  559. bnput16("a = ", C, i);
  560. bnput16("b = ", B, i);
  561. printf("t = %lXn", (unsigned long)t);
  562. bnput16("e = b * t = ", E, i+1);
  563. bnput16("a + e = ", D, i);
  564. printf("carry = %lXn", (unsigned long)borrow);
  565. bnput16("a + b * t = ", A, i);
  566. printf("carry = %lXn", (unsigned long)carry);
  567. }
  569. memcpy(d, a, sizeof(a));
  570. borrow = lbnMulSub1_16(A, B, i, t);
  572. /* Did MulSub perform the inverse of MulAdd */
  573. if (carry != borrow || lbnCmp_16(A,C,i) || !j) {
  574. bnput16("a = ", C, i);
  575. bnput16("b = ", B, i);
  576. bnput16("a += b*t = ", D, i+1);
  577. printf("Carry = %lXn", (unsigned long)carry);
  578. bnput16("a -= b*t = ", A, i+1);
  579. printf("Borrow = %lXn", (unsigned long)borrow);
  580. printf("(a <=> c) = %dn", lbnCmp_16(A,C,i));
  581. bnput16("b*t = ", E, i+1);
  582. }
  584. /* At this point we're done with t, so it's scratch */
  586. /* Does Mul work both ways symmetrically */
  587. lbnMul_16(C,A,i,B,i);
  588. lbnMul_16(D,B,i,A,i);
  589. if (lbnCmp_16(C,D,i+i) || !j) {
  590. bnput16("a = ", A, i);
  591. bnput16("b = ", B, i);
  592. bnput16("a * b = ", C, i+i);
  593. bnput16("b * a = ", D, i+i);
  594. printf("(a*b <=> b*a) = %dn",
  595.        lbnCmp_16(C,D,i+i));
  596. }
  598. /* Generate an F less than A and B */
  599. do {
  600. randnum(F,i);
  601. } while (lbnCmp_16(F,A,i) >= 0 ||
  602.          lbnCmp_16(F,B,i) >= 0);
  604. /* Add F to D */
  605. lbnAdd1_16(BIGLITTLE(D-i,D+i), i, lbnAddN_16(D, F, i));
  606. memcpy(c, d, sizeof(d));
  608. /*
  609.  * Divide by A and check that quotient and remainder
  610.  * match
  611.  */
  612. t = lbnDiv_16(E,C,i+i,A,i);
  613. if (t || lbnCmp_16(E,B,i) || lbnCmp_16(C, F, i) || !j) {
  614. bnput16("a = ", A, i);
  615. bnput16("b = ", B, i);
  616. bnput16("f = ", F, i);
  617. bnput16("a * b + f = ", D, i+i);
  618. printf("qhigh = %lXn", (unsigned long)t);
  619. bnput16("(a*b+f) / a = ", E, i);
  620. bnput16("(a*b+f) % a = ", C, i);
  621. }
  623. memcpy(c, d, sizeof(d));
  625. /* Divide by B and check similarly */
  626. t = lbnDiv_16(E,C,i+i,B,i);
  627. if (lbnCmp_16(E,A,i) || lbnCmp_16(C, F, i) || !j) {
  628. bnput16("a = ", A, i);
  629. bnput16("b = ", B, i);
  630. bnput16("f = ", F, i);
  631. bnput16("a * b + f = ", D, i+i);
  632. printf("qhigh = %lXn", (unsigned long)t);
  633. bnput16("(a*b+f) / b = ", E, i);
  634. bnput16("(a*b+f) % b = ", C, i);
  635. }
  637. /* Checl that multiply and square do the same thing */
  638. lbnMul_16(C,A,i,A,i);
  639. lbnSquare_16(D,A,i);
  640. if (lbnCmp_16(C,D,i+i) || !j) {
  641. bnput16("a*a = ", C, i+i);
  642. bnput16("a^2 = ", D, i+i);
  643. printf("(a * a == a^2) = %dn",
  644.        lbnCmp_16(C,D,i+i));
  645. }
  647. /* Compute a GCD */
  648. lbnCopy_16(C,A,i);
  649. lbnCopy_16(D,B,i);
  650. z = lbnGcd_16(C,i,D,i);
  652. /* Approximate check that the GCD came out right */
  653. for (l = 0;
  654.      l < sizeof(smallprimes)/sizeof(*smallprimes);
  655.      l++)
  656. {
  657. t = lbnModQ_16(z > 0 ? C : D,
  658.                (unsigned)(z > 0 ? z : -z),
  659.        smallprimes[l]);
  660. carry = lbnModQ_16(A, i, smallprimes[l]);
  661. borrow = lbnModQ_16(B, i, smallprimes[l]);
  662. if (!t != (!carry && !borrow)) {
  663. bnput16("a = ", A, i);
  664. printf("a mod %u = %lun",
  665.        smallprimes[l],
  666.        (unsigned long)carry);
  667. bnput16("b = ", B, i);
  668. printf("b mod %u = %lun",
  669.        smallprimes[l],
  670.        (unsigned long)borrow);
  671. bnput16("gcd(a,b) = ", z > 0 ? C : D,
  672. (unsigned)(z > 0 ? z : -z));
  673. printf("gcd(a,b) mod %u = %lun",
  674.        smallprimes[l],
  675.        (unsigned long)t);
  676. }
  677. }
  679. /*
  680.  * Do some Montgomery operations
  681.  * Start with A > B, and also place a copy of B
  682.  * into C.
  683.  */
  684. if (lbnCmp_16(A, B, i) < 0) {
  685. memcpy(c, a, sizeof(c));
  686. memcpy(a, b, sizeof(a));
  687. memcpy(b, c, sizeof(b));
  688. } else {
  689. memcpy(c, b, sizeof(c));
  690. }
  692. /* Convert to and from */
  693. BIGLITTLE(A[-1],A[0]) |= 1;
  694. lbnToMont_16(B, i, A, i);
  695. lbnFromMont_16(B, A, i);
  696. if (lbnCmp_16(B, C, i)) {
  697. memcpy(b, c, sizeof(c));
  698. bnput16("mod = ", A, i);
  699. bnput16("input = ", B, i);
  700. lbnToMont_16(B, i, A, i);
  701. bnput16("mont = ", B, i);
  702. lbnFromMont_16(B, A, i);
  703. bnput16("output = ", B, i);
  704. }
  705. /* E = B^5 (mod A), the simple way */
  706. lbnSquare_16(E, B, i);
  707. (void)lbnDiv_16(BIGLITTLE(E-i,E+i),E,i+i,A,i);
  708. lbnSquare_16(D, E, i);
  709. (void)lbnDiv_16(BIGLITTLE(D-i,D+i),D,i+i,A,i);
  710. lbnMul_16(E, D, i, B, i);
  711. (void)lbnDiv_16(BIGLITTLE(E-i,E+i),E,i+i,A,i);
  713. /* D = B^5, using ExpMod */
  714. BIGLITTLE(F[-1],F[0]) = 5;
  715. lbnExpMod_16(D, B, i, F, 1, A, i);
  716. if (lbnCmp_16(D, E, i)  || !j) {
  717. bnput16("mod = ", A, i);
  718. bnput16("input = ", B, i);
  719. bnput16("input^5 = ", E, i);
  720. bnput16("input^5 = ", D, i);
  721. printf("a>b (x <=> y) = %dn",
  722.        lbnCmp_16(D,E,i));
  723. }
  724. /* TODO: Test lbnTwoExpMod, lbnDoubleExpMod */
  725. } /* for (i) */
  726. printf("r%d ", j);
  727. fflush(stdout);
  728.         } /* for (j) */
  729.         printf("%d iterations of up to %d 16-bit words completed.n",
  730.        j-1, i-1);
  731.         return 0;
  732. }