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