home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / crypl200.zip / BNLIB / LBN64.C < prev    next >
C/C++ Source or Header  |  1996-09-22  |  109KB  |  3,615 lines
  1. /*
  2.  * lbn64.c - Low-level bignum routines, 64-bit version.
  3.  *
  4.  * Copyright (c) 1995  Colin Plumb.  All rights reserved.
  5.  * For licensing and other legal details, see the file legal.c.
  6.  *
  7.  * NOTE: the magic constants "64" and "128" appear in many places in this
  8.  * file, including inside identifiers.  Because it is not possible to
  9.  * ask "#ifdef" of a macro expansion, it is not possible to use the
  10.  * preprocessor to conditionalize these properly.  Thus, this file is
  11.  * intended to be edited with textual search and replace to produce
  12.  * alternate word size versions.  Any reference to the number of bits
  13.  * in a word must be the string "64", and that string must not appear
  14.  * otherwise.  Any reference to twice this number must appear as "128",
  15.  * which likewise must not appear otherwise.  Is that clear?
  16.  *
  17.  * Remember, when doubling the bit size replace the larger number (128)
  18.  * first, then the smaller (64).  When halving the bit size, do the
  19.  * opposite.  Otherwise, things will get wierd.  Also, be sure to replace
  20.  * every instance that appears.  (:%s/foo/bar/g in vi)
  21.  *
  22.  * These routines work with a pointer to the least-significant end of
  23.  * an array of WORD64s.  The BIG(x), LITTLE(y) and BIGLTTLE(x,y) macros
  24.  * defined in lbn.h (which expand to x on a big-edian machine and y on a
  25.  * little-endian machine) are used to conditionalize the code to work
  26.  * either way.  If you have no assembly primitives, it doesn't matter.
  27.  * Note that on a big-endian machine, the least-significant-end pointer
  28.  * is ONE PAST THE END.  The bytes are ptr[-1] through ptr[-len].
  29.  * On little-endian, they are ptr[0] through ptr[len-1].  This makes
  30.  * perfect sense if you consider pointers to point *between* bytes rather
  31.  * than at them.
  32.  *
  33.  * Because the array index values are unsigned integers, ptr[-i]
  34.  * may not work properly, since the index -i is evaluated as an unsigned,
  35.  * and if pointers are wider, zero-extension will produce a positive
  36.  * number rahter than the needed negative.  The expression used in this
  37.  * code, *(ptr-i) will, however, work.  (The array syntax is equivalent
  38.  * to *(ptr+-i), which is a pretty subtle difference.)
  39.  *
  40.  * Many of these routines will get very unhappy if fed zero-length inputs.
  41.  * They use assert() to enforce this.  An higher layer of code must make
  42.  * sure that these aren't called with zero-length inputs.
  43.  *
  44.  * Any of these routines can be replaced with more efficient versions
  45.  * elsewhere, by just #defining their names.  If one of the names
  46.  * is #defined, the C code is not compiled in and no declaration is
  47.  * made.  Use the BNINCLUDE file to do that.  Typically, you compile
  48.  * asm subroutines with the same name and just, e.g.
  49.  * #define lbnMulAdd1_64 lbnMulAdd1_64
  50.  *
  51.  * If you want to write asm routines, start with lbnMulAdd1_64().
  52.  * This is the workhorse of modular exponentiation.  lbnMulN1_64() is
  53.  * also used a fair bit, although not as much and it's defined in terms
  54.  * of lbnMulAdd1_64 if that has a custom version.  lbnMulSub1_64 and
  55.  * lbnDiv21_64 are used in the usual division and remainder finding.
  56.  * (Not the Montgomery reduction used in modular exponentiation, though.)
  57.  * Once you have lbnMulAdd1_64 defined, writing the other two should
  58.  * be pretty easy.  (Just make sure you get the sign of the subtraction
  59.  * in lbnMulSub1_64 right - it's dest = dest - source * k.)
  60.  *
  61.  * The only definitions that absolutely need a double-word (BNWORD128)
  62.  * type are lbnMulAdd1_64 and lbnMulSub1_64; if those are provided,
  63.  * the rest follows.  lbnDiv21_64, however, is a lot slower unless you
  64.  * have them, and lbnModQ_64 takes after it.  That one is used quite a
  65.  * bit for prime sieving.
  66.  */
  67.  
  68. #ifndef HAVE_CONFIG_H
  69. #define HAVE_CONFIG_H 0
  70. #endif
  71. #if HAVE_CONFIG_H
  72. #if defined( INC_ALL ) || defined( INC_CHILD )
  73.   #include "config.h"
  74. #else
  75.   #include "bnlib/config.h"
  76. #endif /* Compiler-specific includes */
  77. #endif
  78.  
  79. /*
  80.  * Some compilers complain about #if FOO if FOO isn't defined,
  81.  * so do the ANSI-mandated thing explicitly...
  82.  */
  83. #ifndef NO_ASSERT_H
  84. #define NO_ASSERT_H 0
  85. #endif
  86. #ifndef NO_STRING_H
  87. #define NO_STRING_H 0
  88. #endif
  89. #ifndef HAVE_STRINGS_H
  90. #define HAVE_STRINGS_H 0
  91. #endif
  92. #ifndef NEED_MEMORY_H
  93. #define NEED_MEMORY_H 0
  94. #endif
  95.  
  96. #if !NO_ASSERT_H
  97. #include <assert.h>
  98. #else
  99. #define assert(x) (void)0
  100. #endif
  101.  
  102. #if !NO_STRING_H
  103. #include <string.h>    /* For memcpy */
  104. #elif HAVE_STRINGS_H
  105. #include <strings.h>
  106. #endif
  107. #if NEED_MEMORY_H
  108. #include <memory.h>
  109. #endif
  110.  
  111. #if defined( INC_ALL ) || defined( INC_CHILD )
  112.   #include "lbn.h"
  113.   #include "lbn64.h"
  114.   #include "lbnmem.h"
  115.   #include "legal.h"
  116.  
  117.   #include "kludge.h"
  118. #else
  119.   #include "bnlib/lbn.h"
  120.   #include "bnlib/lbn64.h"
  121.   #include "bnlib/lbnmem.h"
  122.   #include "bnlib/legal.h"
  123.  
  124.   #include "bnlib/kludge.h"
  125. #endif /* Compiler-specific includes */
  126.  
  127. #ifndef BNWORD64
  128. #error 64-bit bignum library requires a 64-bit data type
  129. #endif
  130.  
  131. /* Make sure the copyright notice gets included */
  132. volatile const char * volatile const lbnCopyright_64 = bnCopyright;
  133.  
  134. /*
  135.  * Most of the multiply (and Montgomery reduce) routines use an outer
  136.  * loop that iterates over one of the operands - a so-called operand
  137.  * scanning approach.  One big advantage of this is that the assembly
  138.  * support routines are simpler.  The loops can be rearranged to have
  139.  * an outer loop that iterates over the product, a so-called product
  140.  * scanning approach.  This has the advantage of writing less data
  141.  * and doing fewer adds to memory, so is supposedly faster.  Some
  142.  * code has been written using a product-scanning approach, but
  143.  * it appears to be slower, so it is turned off by default.  Some
  144.  * experimentation would be appreciated.
  145.  *
  146.  * (The code is also annoying to get right and not very well commented,
  147.  * one of my pet peeves about math libraries.  I'm sorry.)
  148.  */
  149. #ifndef PRODUCT_SCAN
  150. #define PRODUCT_SCAN 0
  151. #endif
  152.  
  153. /*
  154.  * Copy an array of words.  <Marvin mode on>  Thrilling, isn't it? </Marvin>
  155.  * This is a good example of how the byte offsets and BIGLITTLE() macros work.
  156.  * Another alternative would have been
  157.  * memcpy(dest BIG(-len), src BIG(-len), len*sizeof(BNWORD64)), but I find that
  158.  * putting operators into conditional macros is confusing.
  159.  */
  160. #ifndef lbnCopy_64
  161. void
  162. lbnCopy_64(BNWORD64 *dest, BNWORD64 const *src, unsigned len)
  163. {
  164.     memcpy(BIGLITTLE(dest-len,dest), BIGLITTLE(src-len,src),
  165.            len * sizeof(*src));
  166. }
  167. #endif /* !lbnCopy_64 */
  168.  
  169. /*
  170.  * Fill n words with zero.  This does it manually rather than calling
  171.  * memset because it can assume alignment to make things faster while
  172.  * memset can't.  Note how big-endian numbers are naturally addressed
  173.  * using predecrement, while little-endian is postincrement.
  174.  */
  175. #ifndef lbnZero_64
  176. void
  177. lbnZero_64(BNWORD64 *num, unsigned len)
  178. {
  179.     while (len--)
  180.         BIGLITTLE(*--num,*num++) = 0;
  181. }
  182. #endif /* !lbnZero_64 */
  183.  
  184. /*
  185.  * Negate an array of words.
  186.  * Negation is subtraction from zero.  Negating low-order words
  187.  * entails doing nothing until a non-zero word is hit.  Once that
  188.  * is negated, a borrow is generated and never dies until the end
  189.  * of the number is hit.  Negation with borrow, -x-1, is the same as ~x.
  190.  * Repeat that until the end of the number.
  191.  *
  192.  * Doesn't return borrow out because that's pretty useless - it's
  193.  * always set unless the input is 0, which is easy to notice in
  194.  * normalized form.
  195.  */
  196. #ifndef lbnNeg_64
  197. void
  198. lbnNeg_64(BNWORD64 *num, unsigned len)
  199. {
  200.     assert(len);
  201.  
  202.     /* Skip low-order zero words */
  203.     while (BIGLITTLE(*--num,*num) == 0) {
  204.         if (!--len)
  205.             return;
  206.         LITTLE(num++;)
  207.     }
  208.     /* Negate the lowest-order non-zero word */
  209.     *num = -*num;
  210.     /* Complement all the higher-order words */
  211.     while (--len) {
  212.         BIGLITTLE(--num,++num);
  213.         *num = ~*num;
  214.     }
  215. }
  216. #endif /* !lbnNeg_64 */
  217.  
  218.  
  219. /*
  220.  * lbnAdd1_64: add the single-word "carry" to the given number.
  221.  * Used for minor increments and propagating the carry after
  222.  * adding in a shorter bignum.
  223.  *
  224.  * Technique: If we have a double-width word, presumably the compiler
  225.  * can add using its carry in inline code, so we just use a larger
  226.  * accumulator to compute the carry from the first addition.
  227.  * If not, it's more complex.  After adding the first carry, which may
  228.  * be > 1, compare the sum and the carry.  If the sum wraps (causing a
  229.  * carry out from the addition), the result will be less than each of the
  230.  * inputs, since the wrap subtracts a number (2^64) which is larger than
  231.  * the other input can possibly be.  If the sum is >= the carry input,
  232.  * return success immediately.
  233.  * In either case, if there is a carry, enter a loop incrementing words
  234.  * until one does not wrap.  Since we are adding 1 each time, the wrap
  235.  * will be to 0 and we can test for equality.
  236.  */
  237. #ifndef lbnAdd1_64    /* If defined, it's provided as an asm subroutine */
  238. #ifdef BNWORD128
  239. BNWORD64
  240. lbnAdd1_64(BNWORD64 *num, unsigned len, BNWORD64 carry)
  241. {
  242.     BNWORD128 t;
  243.     assert(len > 0);    /* Alternative: if (!len) return carry */
  244.  
  245.     t = (BNWORD128)BIGLITTLE(*--num,*num) + carry;
  246.     BIGLITTLE(*num,*num++) = (BNWORD64)t;
  247.     if ((t >> 64) == 0)
  248.         return 0;
  249.     while (--len) {
  250.         if (++BIGLITTLE(*--num,*num++) != 0)
  251.             return 0;
  252.     }
  253.     return 1;
  254. }
  255. #else /* no BNWORD128 */
  256. BNWORD64
  257. lbnAdd1_64(BNWORD64 *num, unsigned len, BNWORD64 carry)
  258. {
  259.     assert(len > 0);    /* Alternative: if (!len) return carry */
  260.  
  261.     if ((BIGLITTLE(*--num,*num++) += carry) >= carry)
  262.         return 0;
  263.     while (--len) {
  264.         if (++BIGLITTLE(*--num,*num++) != 0)
  265.             return 0;
  266.     }
  267.     return 1;
  268. }
  269. #endif
  270. #endif/* !lbnAdd1_64 */
  271.  
  272. /*
  273.  * lbnSub1_64: subtract the single-word "borrow" from the given number.
  274.  * Used for minor decrements and propagating the borrow after
  275.  * subtracting a shorter bignum.
  276.  *
  277.  * Technique: Similar to the add, above.  If there is a double-length type,
  278.  * use that to generate the first borrow.
  279.  * If not, after subtracting the first borrow, which may be > 1, compare
  280.  * the difference and the *negative* of the carry.  If the subtract wraps
  281.  * (causing a borrow out from the subtraction), the result will be at least
  282.  * as large as -borrow.  If the result < -borrow, then no borrow out has
  283.  * appeared and we may return immediately, except when borrow == 0.  To
  284.  * deal with that case, use the identity that -x = ~x+1, and instead of
  285.  * comparing < -borrow, compare for <= ~borrow.
  286.  * Either way, if there is a borrow out, enter a loop decrementing words
  287.  * until a non-zero word is reached.
  288.  *
  289.  * Note the cast of ~borrow to (BNWORD64).  If the size of an int is larger
  290.  * than BNWORD64, C rules say the number is expanded for the arithmetic, so
  291.  * the inversion will be done on an int and the value won't be quite what
  292.  * is expected.
  293.  */
  294. #ifndef lbnSub1_64    /* If defined, it's provided as an asm subroutine */
  295. #ifdef BNWORD128
  296. BNWORD64
  297. lbnSub1_64(BNWORD64 *num, unsigned len, BNWORD64 borrow)
  298. {
  299.     BNWORD128 t;
  300.     assert(len > 0);    /* Alternative: if (!len) return borrow */
  301.  
  302.     t = (BNWORD128)BIGLITTLE(*--num,*num) - borrow;
  303.     BIGLITTLE(*num,*num++) = (BNWORD64)t;
  304.     if ((t >> 64) == 0)
  305.         return 0;
  306.     while (--len) {
  307.         if ((BIGLITTLE(*--num,*num++))-- != 0)
  308.             return 0;
  309.     }
  310.     return 1;
  311. }
  312. #else /* no BNWORD128 */
  313. BNWORD64
  314. lbnSub1_64(BNWORD64 *num, unsigned len, BNWORD64 borrow)
  315. {
  316.     assert(len > 0);    /* Alternative: if (!len) return borrow */
  317.  
  318.     if ((BIGLITTLE(*--num,*num++) -= borrow) <= (BNWORD64)~borrow)
  319.         return 0;
  320.     while (--len) {
  321.         if ((BIGLITTLE(*--num,*num++))-- != 0)
  322.             return 0;
  323.     }
  324.     return 1;
  325. }
  326. #endif
  327. #endif /* !lbnSub1_64 */
  328.  
  329. /*
  330.  * lbnAddN_64: add two bignums of the same length, returning the carry (0 or 1).
  331.  * One of the building blocks, along with lbnAdd1, of adding two bignums of
  332.  * differing lengths.
  333.  *
  334.  * Technique: Maintain a word of carry.  If there is no double-width type,
  335.  * use the same technique as in lbnAdd1, above, to maintain the carry by
  336.  * comparing the inputs.  Adding the carry sources is used as an OR operator;
  337.  * at most one of the two comparisons can possibly be true.  The first can
  338.  * only be true if carry == 1 and x, the result, is 0.  In that case the
  339.  * second can't possibly be true.
  340.  */
  341. #ifndef lbnAddN_64
  342. #ifdef BNWORD128
  343. BNWORD64
  344. lbnAddN_64(BNWORD64 *num1, BNWORD64 const *num2, unsigned len)
  345. {
  346.     BNWORD128 t;
  347.  
  348.     assert(len > 0);
  349.  
  350.     t = (BNWORD128)BIGLITTLE(*--num1,*num1) + BIGLITTLE(*--num2,*num2++);
  351.     BIGLITTLE(*num1,*num1++) = (BNWORD64)t;
  352.     while (--len) {
  353.         t = (BNWORD128)BIGLITTLE(*--num1,*num1) +
  354.             (BNWORD128)BIGLITTLE(*--num2,*num2++) + (t >> 64);
  355.         BIGLITTLE(*num1,*num1++) = (BNWORD64)t;
  356.     }
  357.  
  358.     return (BNWORD64)(t>>64);
  359. }
  360. #else /* no BNWORD128 */
  361. BNWORD64
  362. lbnAddN_64(BNWORD64 *num1, BNWORD64 const *num2, unsigned len)
  363. {
  364.     BNWORD64 x, carry = 0;
  365.  
  366.     assert(len > 0);    /* Alternative: change loop to test at start */
  367.  
  368.     do {
  369.         x = BIGLITTLE(*--num2,*num2++);
  370.         carry = (x += carry) < carry;
  371.         carry += (BIGLITTLE(*--num1,*num1++) += x) < x;
  372.     } while (--len);
  373.  
  374.     return carry;
  375. }
  376. #endif
  377. #endif /* !lbnAddN_64 */
  378.  
  379. /*
  380.  * lbnSubN_64: add two bignums of the same length, returning the carry (0 or 1).
  381.  * One of the building blocks, along with subn1, of subtracting two bignums of
  382.  * differing lengths.
  383.  *
  384.  * Technique: If no double-width type is availble, maintain a word of borrow.
  385.  * First, add the borrow to the subtrahend (did you have to learn all those
  386.  * awful words in elementary school, too?), and if it overflows, set the
  387.  * borrow again.  Then subtract the modified subtrahend from the next word
  388.  * of input, using the same technique as in subn1, above.
  389.  * Adding the borrows is used as an OR operator; at most one of the two
  390.  * comparisons can possibly be true.  The first can only be true if
  391.  * borrow == 1 and x, the result, is 0.  In that case the second can't
  392.  * possibly be true.
  393.  *
  394.  * In the double-word case, (BNWORD64)-(t>>64) is subtracted, rather than
  395.  * adding t>>64, because the shift would need to sign-extend and that's
  396.  * not guaranteed to happen in ANSI C, even with signed types.
  397.  */
  398. #ifndef lbnSubN_64
  399. #ifdef BNWORD128
  400. BNWORD64
  401. lbnSubN_64(BNWORD64 *num1, BNWORD64 const *num2, unsigned len)
  402. {
  403.     BNWORD128 t;
  404.  
  405.     assert(len > 0);
  406.  
  407.     t = (BNWORD128)BIGLITTLE(*--num1,*num1) - BIGLITTLE(*--num2,*num2++);
  408.     BIGLITTLE(*num1,*num1++) = (BNWORD64)t;
  409.  
  410.     while (--len) {
  411.         t = (BNWORD128)BIGLITTLE(*--num1,*num1) -
  412.             (BNWORD128)BIGLITTLE(*--num2,*num2++) - (BNWORD64)-(t >> 64);
  413.         BIGLITTLE(*num1,*num1++) = (BNWORD64)t;
  414.     }
  415.  
  416.     return -(BNWORD64)(t>>64);
  417. }
  418. #else
  419. BNWORD64
  420. lbnSubN_64(BNWORD64 *num1, BNWORD64 const *num2, unsigned len)
  421. {
  422.     BNWORD64 x, borrow = 0;
  423.  
  424.     assert(len > 0);    /* Alternative: change loop to test at start */
  425.  
  426.     do {
  427.         x = BIGLITTLE(*--num2,*num2++);
  428.         borrow = (x += borrow) < borrow;
  429.         borrow += (BIGLITTLE(*--num1,*num1++) -= x) > (BNWORD64)~x;
  430.     } while (--len);
  431.  
  432.     return borrow;
  433. }
  434. #endif
  435. #endif /* !lbnSubN_64 */
  436.  
  437. #ifndef lbnCmp_64
  438. /*
  439.  * lbnCmp_64: compare two bignums of equal length, returning the sign of
  440.  * num1 - num2. (-1, 0 or +1).
  441.  *
  442.  * Technique: Change the little-endian pointers to big-endian pointers
  443.  * and compare from the most-significant end until a difference if found.
  444.  * When it is, figure out the sign of the difference and return it.
  445.  */
  446. int
  447. lbnCmp_64(BNWORD64 const *num1, BNWORD64 const *num2, unsigned len)
  448. {
  449.     BIGLITTLE(num1 -= len, num1 += len);
  450.     BIGLITTLE(num2 -= len, num2 += len);
  451.  
  452.     while (len--) {
  453.         if (BIGLITTLE(*num1++ != *num2++, *--num1 != *--num2)) {
  454.             if (BIGLITTLE(num1[-1] < num2[-1], *num1 < *num2))
  455.                 return -1;
  456.             else
  457.                 return 1;
  458.         }
  459.     }
  460.     return 0;
  461. }
  462. #endif /* !lbnCmp_64 */
  463.  
  464. /*
  465.  * mul64_ppmmaa(ph,pl,x,y,a,b) is an optional routine that
  466.  * computes (ph,pl) = x * y + a + b.  mul64_ppmma and mul64_ppmm
  467.  * are simpler versions.  If you want to be lazy, all of these
  468.  * can be defined in terms of the others, so here we create any
  469.  * that have not been defined in terms of the ones that have been.
  470.  */
  471.  
  472. /* Define ones with fewer a's in terms of ones with more a's */
  473. #if !defined(mul64_ppmma) && defined(mul64_ppmmaa)
  474. #define mul64_ppmma(ph,pl,x,y,a) mul64_ppmmaa(ph,pl,x,y,a,0)
  475. #endif
  476.  
  477. #if !defined(mul64_ppmm) && defined(mul64_ppmma)
  478. #define mul64_ppmm(ph,pl,x,y) mul64_ppmma(ph,pl,x,y,0)
  479. #endif
  480.  
  481. /*
  482.  * Use this definition to test the mul64_ppmm-based operations on machines
  483.  * that do not provide mul64_ppmm.  Change the final "0" to a "1" to
  484.  * enable it.
  485.  */
  486. #if !defined(mul64_ppmm) && defined(BNWORD128) && 0    /* Debugging */
  487. #define mul64_ppmm(ph,pl,x,y) \
  488.     ({BNWORD128 _ = (BNWORD128)(x)*(y); (pl) = _; (ph) = _>>64;})
  489. #endif
  490.  
  491. #if defined(mul64_ppmm) && !defined(mul64_ppmma)
  492. #define mul64_ppmma(ph,pl,x,y,a) \
  493.     (mul64_ppmm(ph,pl,x,y), (ph) += ((pl) += (a)) < (a))
  494. #endif
  495.  
  496. #if defined(mul64_ppmma) && !defined(mul64_ppmmaa)
  497. #define mul64_ppmmaa(ph,pl,x,y,a,b) \
  498.     (mul64_ppmma(ph,pl,x,y,a), (ph) += ((pl) += (b)) < (b))
  499. #endif
  500.  
  501. /*
  502.  * lbnMulN1_64: Multiply an n-word input by a 1-word input and store the
  503.  * n+1-word product.  This uses either the mul64_ppmm and mul64_ppmma
  504.  * macros, or C multiplication with the BNWORD128 type.  This uses mul64_ppmma
  505.  * if available, assuming you won't bother defining it unless you can do
  506.  * better than the normal multiplication.
  507.  */
  508. #ifndef lbnMulN1_64
  509. #ifdef lbnMulAdd1_64    /* If we have this asm primitive, use it. */
  510. void
  511. lbnMulN1_64(BNWORD64 *out, BNWORD64 const *in, unsigned len, BNWORD64 k)
  512. {
  513.     lbnZero_64(out, len);
  514.     BIGLITTLE(*(out-len),*(out+len)) = lbnMulAdd1_64(out, in, len, k);
  515. }
  516. #elif defined(mul64_ppmm)
  517. void
  518. lbnMulN1_64(BNWORD64 *out, BNWORD64 const *in, unsigned len, BNWORD64 k)
  519. {
  520.     BNWORD64 prod, carry, carryin;
  521.  
  522.     assert(len > 0);
  523.  
  524.     BIG(--out;--in;);
  525.     mul64_ppmm(carry, *out, *in, k);
  526.     LITTLE(out++;in++;)
  527.  
  528.     while (--len) {
  529.         BIG(--out;--in;)
  530.         carryin = carry;
  531.         mul64_ppmma(carry, *out, *in, k, carryin);
  532.         LITTLE(out++;in++;)
  533.     }
  534.     BIGLITTLE(*--out,*out) = carry;
  535. }
  536. #elif defined(BNWORD128)
  537. void
  538. lbnMulN1_64(BNWORD64 *out, BNWORD64 const *in, unsigned len, BNWORD64 k)
  539. {
  540.     BNWORD128 p;
  541.  
  542.     assert(len > 0);
  543.  
  544.     p = (BNWORD128)BIGLITTLE(*--in,*in++) * k;
  545.     BIGLITTLE(*--out,*out++) = (BNWORD64)p;
  546.  
  547.     while (--len) {
  548.         p = (BNWORD128)BIGLITTLE(*--in,*in++) * k + (BNWORD64)(p >> 64);
  549.         BIGLITTLE(*--out,*out++) = (BNWORD64)p;
  550.     }
  551.     BIGLITTLE(*--out,*out) = (BNWORD64)(p >> 64);
  552. }
  553. #else
  554. #error No 64x64 -> 128 multiply available for 64-bit bignum package
  555. #endif
  556. #endif /* lbnMulN1_64 */
  557.  
  558. /*
  559.  * lbnMulAdd1_64: Multiply an n-word input by a 1-word input and add the
  560.  * low n words of the product to the destination.  *Returns the n+1st word
  561.  * of the product.*  (That turns out to be more convenient than adding
  562.  * it into the destination and dealing with a possible unit carry out
  563.  * of *that*.)  This uses either the mul64_ppmma and mul64_ppmmaa macros,
  564.  * or C multiplication with the BNWORD128 type.
  565.  *
  566.  * If you're going to write assembly primitives, this is the one to
  567.  * start with.  It is by far the most commonly called function.
  568.  */
  569. #ifndef lbnMulAdd1_64
  570. #if defined(mul64_ppmm)
  571. BNWORD64
  572. lbnMulAdd1_64(BNWORD64 *out, BNWORD64 const *in, unsigned len, BNWORD64 k)
  573. {
  574.     BNWORD64 prod, carry, carryin;
  575.  
  576.     assert(len > 0);
  577.  
  578.     BIG(--out;--in;);
  579.     carryin = *out;
  580.     mul64_ppmma(carry, *out, *in, k, carryin);
  581.     LITTLE(out++;in++;)
  582.  
  583.     while (--len) {
  584.         BIG(--out;--in;);
  585.         carryin = carry;
  586.         mul64_ppmmaa(carry, prod, *in, k, carryin, *out);
  587.         *out = prod;
  588.         LITTLE(out++;in++;)
  589.     }
  590.  
  591.     return carry;
  592. }
  593. #elif defined(BNWORD128)
  594. BNWORD64
  595. lbnMulAdd1_64(BNWORD64 *out, BNWORD64 const *in, unsigned len, BNWORD64 k)
  596. {
  597.     BNWORD128 p;
  598.  
  599.     assert(len > 0);
  600.  
  601.     p = (BNWORD128)BIGLITTLE(*--in,*in++) * k + BIGLITTLE(*--out,*out);
  602.     BIGLITTLE(*out,*out++) = (BNWORD64)p;
  603.  
  604.     while (--len) {
  605.         p = (BNWORD128)BIGLITTLE(*--in,*in++) * k +
  606.             (BNWORD64)(p >> 64) + BIGLITTLE(*--out,*out);
  607.         BIGLITTLE(*out,*out++) = (BNWORD64)p;
  608.     }
  609.  
  610.     return (BNWORD64)(p >> 64);
  611. }
  612. #else
  613. #error No 64x64 -> 128 multiply available for 64-bit bignum package
  614. #endif
  615. #endif /* lbnMulAdd1_64 */
  616.  
  617. /*
  618.  * lbnMulSub1_64: Multiply an n-word input by a 1-word input and subtract the
  619.  * n-word product from the destination.  Returns the n+1st word of the product.
  620.  * This uses either the mul64_ppmm and mul64_ppmma macros, or
  621.  * C multiplication with the BNWORD128 type.
  622.  *
  623.  * This is rather uglier than adding, but fortunately it's only used in
  624.  * division which is not used too heavily.
  625.  */
  626. #ifndef lbnMulN1_64
  627. #if defined(mul64_ppmm)
  628. BNWORD64
  629. lbnMulSub1_64(BNWORD64 *out, BNWORD64 const *in, unsigned len, BNWORD64 k)
  630. {
  631.     BNWORD64 prod, carry, carryin;
  632.  
  633.     assert(len > 0);
  634.  
  635.     BIG(--in;)
  636.     mul64_ppmm(carry, prod, *in, k);
  637.     LITTLE(in++;)
  638.     carry += (BIGLITTLE(*--out,*out++) -= prod) > (BNWORD64)~prod;
  639.  
  640.     while (--len) {
  641.         BIG(--in;);
  642.         carryin = carry;
  643.         mul64_ppmma(carry, prod, *in, k, carryin);
  644.         LITTLE(in++;)
  645.         carry += (BIGLITTLE(*--out,*out++) -= prod) > (BNWORD64)~prod;
  646.     }
  647.  
  648.     return carry;
  649. }
  650. #elif defined(BNWORD128)
  651. BNWORD64
  652. lbnMulSub1_64(BNWORD64 *out, BNWORD64 const *in, unsigned len, BNWORD64 k)
  653. {
  654.     BNWORD128 p;
  655.     BNWORD64 carry, t;
  656.  
  657.     assert(len > 0);
  658.  
  659.     p = (BNWORD128)BIGLITTLE(*--in,*in++) * k;
  660.     t = BIGLITTLE(*--out,*out);
  661.     carry = (BNWORD64)(p>>64) + ((BIGLITTLE(*out,*out++)=t-(BNWORD64)p) > t);
  662.  
  663.     while (--len) {
  664.         p = (BNWORD128)BIGLITTLE(*--in,*in++) * k + carry;
  665.         t = BIGLITTLE(*--out,*out);
  666.         carry = (BNWORD64)(p>>64) +
  667.             ( (BIGLITTLE(*out,*out++)=t-(BNWORD64)p) > t );
  668.     }
  669.  
  670.     return carry;
  671. }
  672. #else
  673. #error No 64x64 -> 128 multiply available for 64-bit bignum package
  674. #endif
  675. #endif /* !lbnMulSub1_64 */
  676.  
  677. /*
  678.  * Shift n words left "shift" bits.  0 < shift < 64.  Returns the
  679.  * carry, any bits shifted off the left-hand side (0 <= carry < 2^shift).
  680.  */
  681. #ifndef lbnLshift_64
  682. BNWORD64
  683. lbnLshift_64(BNWORD64 *num, unsigned len, unsigned shift)
  684. {
  685.     BNWORD64 x, carry;
  686.  
  687.     assert(shift > 0);
  688.     assert(shift < 64);
  689.  
  690.     carry = 0;
  691.     while (len--) {
  692.         BIG(--num;)
  693.         x = *num;
  694.         *num = (x<<shift) | carry;
  695.         LITTLE(num++;)
  696.         carry = x >> (64-shift);
  697.     }
  698.     return carry;
  699. }
  700. #endif /* !lbnLshift_64 */
  701.  
  702. /*
  703.  * An optimized version of the above, for shifts of 1.
  704.  * Some machines can use add-with-carry tricks for this.
  705.  */
  706. #ifndef lbnDouble_64
  707. BNWORD64
  708. lbnDouble_64(BNWORD64 *num, unsigned len)
  709. {
  710.     BNWORD64 x, carry;
  711.  
  712.     carry = 0;
  713.     while (len--) {
  714.         BIG(--num;)
  715.         x = *num;
  716.         *num = (x<<1) | carry;
  717.         LITTLE(num++;)
  718.         carry = x >> (64-1);
  719.     }
  720.     return carry;
  721. }
  722. #endif /* !lbnDouble_64 */
  723.  
  724. /*
  725.  * Shift n words right "shift" bits.  0 < shift < 64.  Returns the
  726.  * carry, any bits shifted off the right-hand side (0 <= carry < 2^shift).
  727.  */
  728. #ifndef lbnRshift_64
  729. BNWORD64
  730. lbnRshift_64(BNWORD64 *num, unsigned len, unsigned shift)
  731. {
  732.     BNWORD64 x, carry = 0;
  733.  
  734.     assert(shift > 0);
  735.     assert(shift < 64);
  736.  
  737.     BIGLITTLE(num -= len, num += len);
  738.  
  739.     while (len--) {
  740.         LITTLE(--num;)
  741.         x = *num;
  742.         *num = (x>>shift) | carry;
  743.         BIG(num++;)
  744.         carry = x << (64-shift);
  745.     }
  746.     return carry >> (64-shift);
  747. }
  748. #endif /* !lbnRshift_64 */
  749.  
  750. /*
  751.  * Multiply two numbers of the given lengths.  prod and num2 may overlap,
  752.  * provided that the low len1 bits of prod are free.  (This corresponds
  753.  * nicely to the place the result is returned from lbnMontReduce_64.)
  754.  *
  755.  * TODO: Use Karatsuba multiply.  The overlap constraints may have
  756.  * to get rewhacked.
  757.  */
  758. #ifndef lbnMul_64
  759. void
  760. lbnMul_64(BNWORD64 *prod, BNWORD64 const *num1, unsigned len1,
  761.                           BNWORD64 const *num2, unsigned len2)
  762. {
  763.     /* Special case of zero */
  764.     if (!len1 || !len2) {
  765.         lbnZero_64(prod, len1+len2);
  766.         return;
  767.     }
  768.  
  769.     /* Multiply first word */
  770.     lbnMulN1_64(prod, num1, len1, BIGLITTLE(*--num2,*num2++));
  771.  
  772.     /*
  773.      * Add in subsequent words, storing the most significant word,
  774.      * which is new each time.
  775.      */
  776.     while (--len2) {
  777.         BIGLITTLE(--prod,prod++);
  778.         BIGLITTLE(*(prod-len1-1),*(prod+len1)) =
  779.             lbnMulAdd1_64(prod, num1, len1, BIGLITTLE(*--num2,*num2++));
  780.     }
  781. }
  782. #endif /* !lbnMul_64 */
  783.  
  784. /*
  785.  * lbnMulX_64 is a square multiply - both inputs are the same length.
  786.  * It's normally just a macro wrapper around the general multiply,
  787.  * but might be implementable in assembly more efficiently (such as
  788.  * when product scanning).
  789.  */
  790. #ifndef lbnMulX_64
  791. #if defined(BNWORD128) && PRODUCT_SCAN
  792. /*
  793.  * Test code to see whether product scanning is any faster.  It seems
  794.  * to make the C code slower, so PRODUCT_SCAN is not defined.
  795.  */
  796. static void
  797. lbnMulX_64(BNWORD64 *prod, BNWORD64 const *num1, BNWORD64 const *num2,
  798.     unsigned len)
  799. {
  800.     BNWORD128 x, y;
  801.     BNWORD64 const *p1, *p2;
  802.     unsigned carry;
  803.     unsigned i, j;
  804.  
  805.     /* Special case of zero */
  806.     if (!len)
  807.         return;
  808.  
  809.     x = (BNWORD128)BIGLITTLE(num1[-1] * num2[-1], num1[0] * num2[0]);
  810.     BIGLITTLE(*--prod, *prod++) = (BNWORD64)x;
  811.     x >>= 64;
  812.  
  813.     for (i = 1; i < len; i++) {
  814.         carry = 0;
  815.         p1 = num1;
  816.         p2 = BIGLITTLE(num2-i-1,num2+i+1);
  817.         for (j = 0; j <= i; j++) {
  818.             BIG(y = (BNWORD128)*--p1 * *p2++;)
  819.             LITTLE(y = (BNWORD128)*p1++ * *--p2;)
  820.             x += y;
  821.             carry += (x < y);
  822.         }
  823.         BIGLITTLE(*--prod,*prod++) = (BNWORD64)x;
  824.         x = (x >> 64) | (BNWORD128)carry << 64;
  825.     }
  826.     for (i = 1; i < len; i++) {
  827.         carry = 0;
  828.         p1 = BIGLITTLE(num1-i,num1+i);
  829.         p2 = BIGLITTLE(num2-len,num2+len);
  830.         for (j = i; j < len; j++) {
  831.             BIG(y = (BNWORD128)*--p1 * *p2++;)
  832.             LITTLE(y = (BNWORD128)*p1++ * *--p2;)
  833.             x += y;
  834.             carry += (x < y);
  835.         }
  836.         BIGLITTLE(*--prod,*prod++) = (BNWORD64)x;
  837.         x = (x >> 64) | (BNWORD128)carry << 64;
  838.     }
  839.     
  840.     BIGLITTLE(*--prod,*prod) = (BNWORD64)x;
  841. }
  842. #else /* !defined(BNWORD128) || !PRODUCT_SCAN */
  843. /* Default trivial macro definition */
  844. #define lbnMulX_64(prod, num1, num2, len) lbnMul_64(prod, num1, len, num2, len)
  845. #endif /* !defined(BNWORD128) || !PRODUCT_SCAN */
  846. #endif /* !lbmMulX_64 */
  847.  
  848. #if !defined(lbnMontMul_64) && defined(BNWORD128) && PRODUCT_SCAN
  849. /*
  850.  * Test code for product-scanning multiply.  This seems to slow the C
  851.  * code down rather than speed it up.
  852.  * This does a multiply and Montgomery reduction together, using the
  853.  * same loops.  The outer loop scans across the product, twice.
  854.  * The first pass computes the low half of the product and the
  855.  * Montgomery multipliers.  These are stored in the product array,
  856.  * which contains no data as of yet.  x and carry add up the columns
  857.  * and propagate carries forward.
  858.  *
  859.  * The second half multiplies the upper half, adding in the modulus
  860.  * times the Montgomery multipliers.  The results of this multiply
  861.  * are stored.
  862.  */
  863. static void
  864. lbnMontMul_64(BNWORD64 *prod, BNWORD64 const *num1, BNWORD64 const *num2,
  865.     BNWORD64 const *mod, unsigned len, BNWORD64 inv)
  866. {
  867.     BNWORD128 x, y;
  868.     BNWORD64 const *p1, *p2, *pm;
  869.     BNWORD64 *pp;
  870.     BNWORD64 t;
  871.     unsigned carry;
  872.     unsigned i, j;
  873.  
  874.     /* Special case of zero */
  875.     if (!len)
  876.         return;
  877.  
  878.     /*
  879.      * This computes directly into the high half of prod, so just
  880.      * shift the pointer and consider prod only "len" elements long
  881.      * for the rest of the code.
  882.      */
  883.     BIGLITTLE(prod -= len, prod += len);
  884.  
  885.     /* Pass 1 - compute Montgomery multipliers */
  886.     /* First iteration can have certain simplifications. */
  887.     x = (BNWORD128)BIGLITTLE(num1[-1] * num2[-1], num1[0] * num2[0]);
  888.     BIGLITTLE(prod[-1], prod[0]) = t = inv * (BNWORD64)x;
  889.     y = (BNWORD128)t * BIGLITTLE(mod[-1],mod[0]);
  890.     x += y;
  891.     /* Note: GCC 2.6.3 has a bug if you try to eliminate "carry" */
  892.     carry = (x < y);
  893.     assert((BNWORD64)x == 0);
  894.     x = x >> 64 | (BNWORD128)carry << 64;
  895.  
  896.     for (i = 1; i < len; i++) {
  897.         carry = 0;
  898.         p1 = num1;
  899.         p2 = BIGLITTLE(num2-i-1,num2+i+1);
  900.         pp = prod;
  901.         pm = BIGLITTLE(mod-i-1,mod+i+1);
  902.         for (j = 0; j < i; j++) {
  903.             y = (BNWORD128)BIGLITTLE(*--p1 * *p2++, *p1++ * *--p2);
  904.             x += y;
  905.             carry += (x < y);
  906.             y = (BNWORD128)BIGLITTLE(*--pp * *pm++, *pp++ * *--pm);
  907.             x += y;
  908.             carry += (x < y);
  909.         }
  910.         y = (BNWORD128)BIGLITTLE(p1[-1] * p2[0], p1[0] * p2[-1]);
  911.         x += y;
  912.         carry += (x < y);
  913.         assert(BIGLITTLE(pp == prod-i, pp == prod+i));
  914.         BIGLITTLE(pp[-1], pp[0]) = t = inv * (BNWORD64)x;
  915.         assert(BIGLITTLE(pm == mod-1, pm == mod+1));
  916.         y = (BNWORD128)t * BIGLITTLE(pm[0],pm[-1]);
  917.         x += y;
  918.         carry += (x < y);
  919.         assert((BNWORD64)x == 0);
  920.         x = x >> 64 | (BNWORD128)carry << 64;
  921.     }
  922.  
  923.     /* Pass 2 - compute reduced product and store */
  924.     for (i = 1; i < len; i++) {
  925.         carry = 0;
  926.         p1 = BIGLITTLE(num1-i,num1+i);
  927.         p2 = BIGLITTLE(num2-len,num2+len);
  928.         pm = BIGLITTLE(mod-i,mod+i);
  929.         pp = BIGLITTLE(prod-len,prod+len);
  930.         for (j = i; j < len; j++) {
  931.             y = (BNWORD128)BIGLITTLE(*--p1 * *p2++, *p1++ * *--p2);
  932.             x += y;
  933.             carry += (x < y);
  934.             y = (BNWORD128)BIGLITTLE(*--pm * *pp++, *pm++ * *--pp);
  935.             x += y;
  936.             carry += (x < y);
  937.         }
  938.         assert(BIGLITTLE(pm == mod-len, pm == mod+len));
  939.         assert(BIGLITTLE(pp == prod-i, pp == prod+i));
  940.         BIGLITTLE(pp[0],pp[-1]) = (BNWORD64)x;
  941.         x = (x >> 64) | (BNWORD128)carry << 64;
  942.     }
  943.  
  944.     /* Last round of second half, simplified. */
  945.     BIGLITTLE(*(prod-len),*(prod+len-1)) = (BNWORD64)x;
  946.     carry = (x >> 64);
  947.  
  948.     while (carry)
  949.         carry -= lbnSubN_64(prod, mod, len);
  950.     while (lbnCmp_64(prod, mod, len) >= 0)
  951.         (void)lbnSubN_64(prod, mod, len);
  952. }
  953. /* Suppress later definition */
  954. #define lbnMontMul_64 lbnMontMul_64
  955. #endif
  956.  
  957. #if !defined(lbnSquare_64) && defined(BNWORD128) && PRODUCT_SCAN
  958. /*
  959.  * Trial code for product-scanning squaring.  This seems to slow the C
  960.  * code down rather than speed it up.
  961.  */
  962. void
  963. lbnSquare_64(BNWORD64 *prod, BNWORD64 const *num, unsigned len)
  964. {
  965.     BNWORD128 x, y, z;
  966.     BNWORD64 const *p1, *p2;
  967.     unsigned carry;
  968.     unsigned i, j;
  969.  
  970.     /* Special case of zero */
  971.     if (!len)
  972.         return;
  973.  
  974.     /* Word 0 of product */
  975.     x = (BNWORD128)BIGLITTLE(num[-1] * num[-1], num[0] * num[0]);
  976.     BIGLITTLE(*--prod, *prod++) = (BNWORD64)x;
  977.     x >>= 64;
  978.  
  979.     /* Words 1 through len-1 */
  980.     for (i = 1; i < len; i++) {
  981.         carry = 0;
  982.         y = 0;
  983.         p1 = num;
  984.         p2 = BIGLITTLE(num-i-1,num+i+1);
  985.         for (j = 0; j < (i+1)/2; j++) {
  986.             BIG(z = (BNWORD128)*--p1 * *p2++;)
  987.             LITTLE(z = (BNWORD128)*p1++ * *--p2;)
  988.             y += z;
  989.             carry += (y < z);
  990.         }
  991.         y += z = y;
  992.         carry += carry + (y < z);
  993.         if ((i & 1) == 0) {
  994.             assert(BIGLITTLE(--p1 == p2, p1 == --p2));
  995.             BIG(z = (BNWORD128)*p2 * *p2;)
  996.             LITTLE(z = (BNWORD128)*p1 * *p1;)
  997.             y += z;
  998.             carry += (y < z);
  999.         }
  1000.         x += y;
  1001.         carry += (x < y);
  1002.         BIGLITTLE(*--prod,*prod++) = (BNWORD64)x;
  1003.         x = (x >> 64) | (BNWORD128)carry << 64;
  1004.     }
  1005.     /* Words len through 2*len-2 */
  1006.     for (i = 1; i < len; i++) {
  1007.         carry = 0;
  1008.         y = 0;
  1009.         p1 = BIGLITTLE(num-i,num+i);
  1010.         p2 = BIGLITTLE(num-len,num+len);
  1011.         for (j = 0; j < (len-i)/2; j++) {
  1012.             BIG(z = (BNWORD128)*--p1 * *p2++;)
  1013.             LITTLE(z = (BNWORD128)*p1++ * *--p2;)
  1014.             y += z;
  1015.             carry += (y < z);
  1016.         }
  1017.         y += z = y;
  1018.         carry += carry + (y < z);
  1019.         if ((len-i) & 1) {
  1020.             assert(BIGLITTLE(--p1 == p2, p1 == --p2));
  1021.             BIG(z = (BNWORD128)*p2 * *p2;)
  1022.             LITTLE(z = (BNWORD128)*p1 * *p1;)
  1023.             y += z;
  1024.             carry += (y < z);
  1025.         }
  1026.         x += y;
  1027.         carry += (x < y);
  1028.         BIGLITTLE(*--prod,*prod++) = (BNWORD64)x;
  1029.         x = (x >> 64) | (BNWORD128)carry << 64;
  1030.     }
  1031.     
  1032.     /* Word 2*len-1 */
  1033.     BIGLITTLE(*--prod,*prod) = (BNWORD64)x;
  1034. }
  1035. /* Suppress later definition */
  1036. #define lbnSquare_64 lbnSquare_64
  1037. #endif
  1038.  
  1039. /*
  1040.  * Square a number, using optimized squaring to reduce the number of
  1041.  * primitive multiples that are executed.  There may not be any
  1042.  * overlap of the input and output.
  1043.  *
  1044.  * Technique: Consider the partial products in the multiplication
  1045.  * of "abcde" by itself:
  1046.  *
  1047.  *               a  b  c  d  e
  1048.  *            *  a  b  c  d  e
  1049.  *          ==================
  1050.  *              ae be ce de ee
  1051.  *           ad bd cd dd de
  1052.  *        ac bc cc cd ce
  1053.  *     ab bb bc bd be
  1054.  *  aa ab ac ad ae
  1055.  *
  1056.  * Note that everything above the main diagonal:
  1057.  *              ae be ce de = (abcd) * e
  1058.  *           ad bd cd       = (abc) * d
  1059.  *        ac bc             = (ab) * c
  1060.  *     ab                   = (a) * b
  1061.  *
  1062.  * is a copy of everything below the main diagonal:
  1063.  *                       de
  1064.  *                 cd ce
  1065.  *           bc bd be
  1066.  *     ab ac ad ae
  1067.  *
  1068.  * Thus, the sum is 2 * (off the diagonal) + diagonal.
  1069.  *
  1070.  * This is accumulated beginning with the diagonal (which
  1071.  * consist of the squares of the digits of the input), which is then
  1072.  * divided by two, the off-diagonal added, and multiplied by two
  1073.  * again.  The low bit is simply a copy of the low bit of the
  1074.  * input, so it doesn't need special care.
  1075.  *
  1076.  * TODO: Merge the shift by 1 with the squaring loop.
  1077.  * TODO: Use Karatsuba.  (a*W+b)^2 = a^2 * (W^2+W) + b^2 * (W+1) - (a-b)^2 * W.
  1078.  */
  1079. #ifndef lbnSquare_64
  1080. void
  1081. lbnSquare_64(BNWORD64 *prod, BNWORD64 const *num, unsigned len)
  1082. {
  1083.     BNWORD64 t;
  1084.     BNWORD64 *prodx = prod;        /* Working copy of the argument */
  1085.     BNWORD64 const *numx = num;    /* Working copy of the argument */
  1086.     unsigned lenx = len;        /* Working copy of the argument */
  1087.  
  1088.     if (!len)
  1089.         return;
  1090.  
  1091.     /* First, store all the squares */
  1092.     while (lenx--) {
  1093. #ifdef mul64_ppmm
  1094.         BNWORD64 ph, pl;
  1095.         t = BIGLITTLE(*--numx,*numx++);
  1096.         mul64_ppmm(ph,pl,t,t);
  1097.         BIGLITTLE(*--prodx,*prodx++) = pl;
  1098.         BIGLITTLE(*--prodx,*prodx++) = ph;
  1099. #elif defined(BNWORD128) /* use BNWORD128 */
  1100.         BNWORD128 p;
  1101.         t = BIGLITTLE(*--numx,*numx++);
  1102.         p = (BNWORD128)t * t;
  1103.         BIGLITTLE(*--prodx,*prodx++) = (BNWORD64)p;
  1104.         BIGLITTLE(*--prodx,*prodx++) = (BNWORD64)(p>>64);
  1105. #else    /* Use lbnMulN1_64 */
  1106.         t = BIGLITTLE(numx[-1],*numx);
  1107.         lbnMulN1_64(prodx, numx, 1, t);
  1108.         BIGLITTLE(--numx,numx++);
  1109.         BIGLITTLE(prodx -= 2, prodx += 2);
  1110. #endif
  1111.     }
  1112.     /* Then, shift right 1 bit */
  1113.     (void)lbnRshift_64(prod, 2*len, 1);
  1114.  
  1115.     /* Then, add in the off-diagonal sums */
  1116.     lenx = len;
  1117.     numx = num;
  1118.     prodx = prod;
  1119.     while (--lenx) {
  1120.         t = BIGLITTLE(*--numx,*numx++);
  1121.         BIGLITTLE(--prodx,prodx++);
  1122.         t = lbnMulAdd1_64(prodx, numx, lenx, t);
  1123.         lbnAdd1_64(BIGLITTLE(prodx-lenx,prodx+lenx), lenx+1, t);
  1124.         BIGLITTLE(--prodx,prodx++);
  1125.     }
  1126.  
  1127.     /* Shift it back up */
  1128.     lbnDouble_64(prod, 2*len);
  1129.  
  1130.     /* And set the low bit appropriately */
  1131.     BIGLITTLE(prod[-1],prod[0]) |= BIGLITTLE(num[-1],num[0]) & 1;
  1132. }
  1133. #endif /* !lbnSquare_64 */
  1134.  
  1135. /*
  1136.  * lbnNorm_64 - given a number, return a modified length such that the
  1137.  * most significant digit is non-zero.  Zero-length input is okay.
  1138.  */
  1139. #ifndef lbnNorm_64
  1140. unsigned
  1141. lbnNorm_64(BNWORD64 const *num, unsigned len)
  1142. {
  1143.     BIGLITTLE(num -= len,num += len);
  1144.     while (len && BIGLITTLE(*num++,*--num) == 0)
  1145.         --len;
  1146.     return len;
  1147. }
  1148. #endif /* lbnNorm_64 */
  1149.  
  1150. /*
  1151.  * lbnBits_64 - return the number of significant bits in the array.
  1152.  * It starts by normalizing the array.  Zero-length input is okay.
  1153.  * Then assuming there's anything to it, it fetches the high word,
  1154.  * generates a bit length by multiplying the word length by 64, and
  1155.  * subtracts off 64/2, 64/4, 64/8, ... bits if the high bits are clear.
  1156.  */
  1157. #ifndef lbnBits_64
  1158. unsigned
  1159. lbnBits_64(BNWORD64 const *num, unsigned len)
  1160. {
  1161.     BNWORD64 t;
  1162.     unsigned i;
  1163.  
  1164.     len = lbnNorm_64(num, len);
  1165.     if (len) {
  1166.         t = BIGLITTLE(*(num-len),*(num+(len-1)));
  1167.         assert(t);
  1168.         len *= 64;
  1169.         i = 64/2;
  1170.         do {
  1171.             if (t >> i)
  1172.                 t >>= i;
  1173.             else
  1174.                 len -= i;
  1175.         } while ((i /= 2) != 0);
  1176.     }
  1177.     return len;
  1178. }
  1179. #endif /* lbnBits_64 */
  1180.  
  1181. /*
  1182.  * If defined, use hand-rolled divide rather than compiler's native.
  1183.  * If the machine doesn't do it in line, the manual code is probably
  1184.  * faster, since it can assume normalization and the fact that the
  1185.  * quotient will fit into 64 bits, which a general 128-bit divide
  1186.  * in a compiler's run-time library can't do.
  1187.  */
  1188. #ifndef BN_SLOW_DIVIDE_128
  1189. #define BN_SLOW_DIVIDE_128 1
  1190. #endif
  1191.  
  1192. /*
  1193.  * Return (nh<<64|nl) % d, and place the quotient digit into *q.
  1194.  * It is guaranteed that nh < d, and that d is normalized (with its high
  1195.  * bit set).  If we have a double-width type, it's easy.  If not, ooh,
  1196.  * yuk!
  1197.  */
  1198. #ifndef lbnDiv21_64
  1199. #if defined(BNWORD128) && !BN_SLOW_DIVIDE_128
  1200. BNWORD64
  1201. lbnDiv21_64(BNWORD64 *q, BNWORD64 nh, BNWORD64 nl, BNWORD64 d)
  1202. {
  1203.     BNWORD128 n = (BNWORD128)nh << 64 | nl;
  1204.  
  1205.     /* Divisor must be normalized */
  1206.     assert(d >> (64-1) == 1);
  1207.  
  1208.     *q = n / d;
  1209.     return n % d;
  1210. }
  1211. #else
  1212. /*
  1213.  * This is where it gets ugly.
  1214.  *
  1215.  * Do the division in two halves, using Algorithm D from section 4.3.1
  1216.  * of Knuth.  Note Theorem B from that section, that the quotient estimate
  1217.  * is never more than the true quotient, and is never more than two
  1218.  * too low.
  1219.  *
  1220.  * The mapping onto conventional long division is (everything a half word):
  1221.  *        _____________qh___ql_
  1222.  * dh dl ) nh.h nh.l nl.h nl.l
  1223.  *             - (qh * d)
  1224.  *            -----------
  1225.  *              rrrr rrrr nl.l
  1226.  *                  - (ql * d)
  1227.  *                -----------
  1228.  *                  rrrr rrrr
  1229.  *
  1230.  * The implicit 3/2-digit d*qh and d*ql subtractors are computed this way:
  1231.  *   First, estimate a q digit so that nh/dh works.  Subtracting qh*dh from
  1232.  *   the (nh.h nh.l) list leaves a 1/2-word remainder r.  Then compute the
  1233.  *   low part of the subtractor, qh * dl.   This also needs to be subtracted
  1234.  *   from (nh.h nh.l nl.h) to get the final remainder.  So we take the
  1235.  *   remainder, which is (nh.h nh.l) - qh*dl, shift it and add in nl.h, and
  1236.  *   try to subtract qh * dl from that.  Since the remainder is 1/2-word
  1237.  *   long, shifting and adding nl.h results in a single word r.
  1238.  *   It is possible that the remainder we're working with, r, is less than
  1239.  *   the product qh * dl, if we estimated qh too high.  The estimation
  1240.  *   technique can produce a qh that is too large (never too small), leading
  1241.  *   to r which is too small.  In that case, decrement the digit qh, add
  1242.  *   shifted dh to r (to correct for that error), and subtract dl from the
  1243.  *   product we're comparing r with.  That's the "correct" way to do it, but
  1244.  *   just adding dl to r instead of subtracting it from the product is
  1245.  *   equivalent and a lot simpler.  You just have to watch out for overflow.
  1246.  *
  1247.  *   The process is repeated with (rrrr rrrr nl.l) for the low digit of the
  1248.  *   quotient ql.
  1249.  *
  1250.  * The various uses of 64/2 for shifts are because of the note about
  1251.  * automatic editing of this file at the very top of the file.
  1252.  */
  1253. #define highhalf(x) ( (x) >> 64/2 )
  1254. #define lowhalf(x) ( (x) & (((BNWORD64)1 << 64/2)-1) )
  1255. BNWORD64
  1256. lbnDiv21_64(BNWORD64 *q, BNWORD64 nh, BNWORD64 nl, BNWORD64 d)
  1257. {
  1258.     BNWORD64 dh = highhalf(d), dl = lowhalf(d);
  1259.     BNWORD64 qh, ql, prod, r;
  1260.  
  1261.     /* Divisor must be normalized */
  1262.     assert((d >> (64-1)) == 1);
  1263.  
  1264.     /* Do first half-word of division */
  1265.     qh = nh / dh;
  1266.     r = nh % dh;
  1267.     prod = qh * dl;
  1268.  
  1269.     /*
  1270.      * Add next half-word of numerator to remainder and correct.
  1271.      * qh may be up to two too large.
  1272.      */
  1273.     r = (r << (64/2)) | highhalf(nl);
  1274.     if (r < prod) {
  1275.         --qh; r += d;
  1276.         if (r >= d && r < prod) {
  1277.             --qh; r += d;
  1278.         }
  1279.     }
  1280.     r -= prod;
  1281.  
  1282.     /* Do second half-word of division */
  1283.     ql = r / dh;
  1284.     r = r % dh;
  1285.     prod = ql * dl;
  1286.  
  1287.     r = (r << (64/2)) | lowhalf(nl);
  1288.     if (r < prod) {
  1289.         --ql; r += d;
  1290.         if (r >= d && r < prod) {
  1291.             --ql; r += d;
  1292.         }
  1293.     }
  1294.     r -= prod;
  1295.  
  1296.     *q = (qh << (64/2)) | ql;
  1297.  
  1298.     return r;
  1299. }
  1300. #endif
  1301. #endif /* lbnDiv21_64 */
  1302.  
  1303.  
  1304. /*
  1305.  * In the division functions, the dividend and divisor are referred to
  1306.  * as "n" and "d", which stand for "numerator" and "denominator".
  1307.  *
  1308.  * The quotient is (nlen-dlen+1) digits long.  It may be overlapped with
  1309.  * the high (nlen-dlen) words of the dividend, but one extra word is needed
  1310.  * on top to hold the top word.
  1311.  */
  1312.  
  1313. /*
  1314.  * Divide an n-word number by a 1-word number, storing the remainder
  1315.  * and n-1 words of the n-word quotient.  The high word is returned.
  1316.  * It IS legal for rem to point to the same address as n, and for
  1317.  * q to point one word higher.
  1318.  *
  1319.  * TODO: If BN_SLOW_DIVIDE_128, add a divnhalf_64 which uses 64-bit
  1320.  *       dividends if the divisor is half that long.
  1321.  * TODO: Shift the dividend on the fly to avoid the last division and
  1322.  *       instead have a remainder that needs shifting.
  1323.  * TODO: Use reciprocals rather than dividing.
  1324.  */
  1325. #ifndef lbnDiv1_64
  1326. BNWORD64
  1327. lbnDiv1_64(BNWORD64 *q, BNWORD64 *rem, BNWORD64 const *n, unsigned len,
  1328.     BNWORD64 d)
  1329. {
  1330.     unsigned shift;
  1331.     unsigned xlen;
  1332.     BNWORD64 r;
  1333.     BNWORD64 qhigh;
  1334.  
  1335.     assert(len > 0);
  1336.     assert(d);
  1337.  
  1338.     if (len == 1) {
  1339.         r = *n;
  1340.         *rem = r%d;
  1341.         return r/d;
  1342.     }
  1343.  
  1344.     shift = 0;
  1345.     r = d;
  1346.     xlen = 64/2;
  1347.     do {
  1348.         if (r >> xlen)
  1349.             r >>= xlen;
  1350.         else
  1351.             shift += xlen;
  1352.     } while ((xlen /= 2) != 0);
  1353.     assert((d >> (64-1-shift)) == 1);
  1354.     d <<= shift;
  1355.  
  1356.     BIGLITTLE(q -= len-1,q += len-1);
  1357.     BIGLITTLE(n -= len,n += len);
  1358.  
  1359.     r = BIGLITTLE(*n++,*--n);
  1360.     if (r < d) {
  1361.         qhigh = 0;
  1362.     } else {
  1363.         qhigh = r/d;
  1364.         r %= d;
  1365.     }
  1366.  
  1367.     xlen = len;
  1368.     while (--xlen)
  1369.         r = lbnDiv21_64(BIGLITTLE(q++,--q), r, BIGLITTLE(*n++,*--n), d);
  1370.  
  1371.     /*
  1372.      * Final correction for shift - shift the quotient up "shift"
  1373.      * bits, and merge in the extra bits of quotient.  Then reduce
  1374.      * the final remainder mod the real d.
  1375.      */
  1376.     if (shift) {
  1377.         d >>= shift;
  1378.         qhigh = (qhigh << shift) | lbnLshift_64(q, len-1, shift);
  1379.         BIGLITTLE(q[-1],*q) |= r/d;
  1380.         r %= d;
  1381.     }
  1382.     *rem = r;
  1383.  
  1384.     return qhigh;
  1385. }
  1386. #endif
  1387.  
  1388. /*
  1389.  * This function performs a "quick" modulus of a number with a divisor
  1390.  * d which is guaranteed to be at most sixteen bits, i.e. less than 65536.
  1391.  * This applies regardless of the word size the library is compiled with.
  1392.  *
  1393.  * This function is important to prime generation, for sieving.
  1394.  */
  1395. #ifndef lbnModQ_64
  1396. /* If there's a custom lbnMod21_64, no normalization needed */
  1397. #ifdef lbnMod21_64
  1398. unsigned
  1399. lbnModQ_64(BNWORD64 const *n, unsigned len, unsigned d)
  1400. {
  1401.     unsigned i, shift;
  1402.     BNWORD64 r;
  1403.  
  1404.     assert(len > 0);
  1405.  
  1406.     BIGLITTLE(n -= len,n += len);
  1407.  
  1408.     /* Try using a compare to avoid the first divide */
  1409.     r = BIGLITTLE(*n++,*--n);
  1410.     if (r >= d)
  1411.         r %= d;
  1412.     while (--len)
  1413.         r = lbnMod21_64(r, BIGLITTLE(*n++,*--n), d);
  1414.  
  1415.     return r;
  1416. }
  1417. #elif defined(BNWORD128) && !BN_SLOW_DIVIDE_128
  1418. unsigned
  1419. lbnModQ_64(BNWORD64 const *n, unsigned len, unsigned d)
  1420. {
  1421.     BNWORD64 r;
  1422.  
  1423.     if (!--len)
  1424.         return BIGLITTLE(n[-1],n[0]) % d;
  1425.  
  1426.     BIGLITTLE(n -= len,n += len);
  1427.     r = BIGLITTLE(n[-1],n[0]);
  1428.  
  1429.     do {
  1430.         r = (BNWORD64)((((BNWORD128)r<<64) | BIGLITTLE(*n++,*--n)) % d);
  1431.     } while (--len);
  1432.  
  1433.     return r;
  1434. }
  1435. #elif 64 >= 0x20
  1436. /*
  1437.  * If the single word size can hold 65535*65536, then this function
  1438.  * is avilable.
  1439.  */
  1440. #ifndef highhalf
  1441. #define highhalf(x) ( (x) >> 64/2 )
  1442. #define lowhalf(x) ( (x) & ((1 << 64/2)-1) )
  1443. #endif
  1444. unsigned
  1445. lbnModQ_64(BNWORD64 const *n, unsigned len, unsigned d)
  1446. {
  1447.     BNWORD64 r, x;
  1448.  
  1449.     BIGLITTLE(n -= len,n += len);
  1450.  
  1451.     r = BIGLITTLE(*n++,*--n);
  1452.     while (--len) {
  1453.         x = BIGLITTLE(*n++,*--n);
  1454.         r = (r%d << 64/2) | highhalf(x);
  1455.         r = (r%d << 64/2) | lowhalf(x);
  1456.     }
  1457.  
  1458.     return r%d;
  1459. }
  1460. #else
  1461. /* Default case - use lbnDiv21_64 */
  1462. unsigned
  1463. lbnModQ_64(BNWORD64 const *n, unsigned len, unsigned d)
  1464. {
  1465.     unsigned i, shift;
  1466.     BNWORD64 r;
  1467.     BNWORD64 q;
  1468.  
  1469.     assert(len > 0);
  1470.  
  1471.     shift = 0;
  1472.     r = d;
  1473.     i = 64;
  1474.     while (i /= 2) {
  1475.         if (r >> i)
  1476.             r >>= i;
  1477.         else
  1478.             shift += i;
  1479.     }
  1480.     assert(d >> (64-1-shift) == 1);
  1481.     d <<= shift;
  1482.  
  1483.     BIGLITTLE(n -= len,n += len);
  1484.  
  1485.     r = BIGLITTLE(*n++,*--n);
  1486.     if (r >= d)
  1487.         r %= d;
  1488.  
  1489.     while (--len)
  1490.         r = lbnDiv21_64(&q, r, BIGLITTLE(*n++,*--n), d);
  1491.  
  1492.     /*
  1493.      * Final correction for shift - shift the quotient up "shift"
  1494.      * bits, and merge in the extra bits of quotient.  Then reduce
  1495.      * the final remainder mod the real d.
  1496.      */
  1497.     if (shift)
  1498.         r %= d >> shift;
  1499.  
  1500.     return r;
  1501. }
  1502. #endif
  1503. #endif /* lbnModQ_64 */
  1504.  
  1505. /*
  1506.  * Reduce n mod d and return the quotient.  That is, find:
  1507.  * q = n / d;
  1508.  * n = n % d;
  1509.  * d is altered during the execution of this subroutine by normalizing it.
  1510.  * It must already have its most significant word non-zero; it is shifted
  1511.  * so its most significant bit is non-zero.
  1512.  *
  1513.  * The quotient q is nlen-dlen+1 words long.  To make it possible to
  1514.  * overlap the quptient with the input (you can store it in the high dlen
  1515.  * words), the high word of the quotient is *not* stored, but is returned.
  1516.  * (If all you want is the remainder, you don't care about it, anyway.)
  1517.  *
  1518.  * This uses algorithm D from Knuth (4.3.1), except that we do binary
  1519.  * (shift) normalization of the divisor.  This is hairy!
  1520.  *
  1521.  * The overall operation is as follows ("top" and "up" refer to the
  1522.  * most significant end of the number; "bottom" and "down", the least):
  1523.  *
  1524.  * - Shift the divisor up until the most significant bit is set.
  1525.  * - Shift the dividend up the same amount.  This will produce the
  1526.  *   correct quotient, and the remainder can be recovered by shifting
  1527.  *   it back down the same number of bits.  This may produce an overflow
  1528.  *   word, but the word is always strictly less than the most significant
  1529.  *   divisor word.
  1530.  * - Estimate the first quotient digit qhat:
  1531.  *   - First take the top two words (one of which is the overflow) of the
  1532.  *     dividend and divide by the top word of the divisor:
  1533.  *     qhat = (nh,nm)/dh.  This qhat is >= the correct quotient digit
  1534.  *     and, since dh is normalized, it is at most two over.
  1535.  *   - Second, correct by comparing the top three words.  If
  1536.  *     (dh,dl) * qhat > (nh,nm,ml), decrease qhat and try again.
  1537.  *     The second iteration can be simpler because there can't be a third.
  1538.  *     The computation can be simplified by subtracting dh*qhat from
  1539.  *     both sides, suitably shifted.  This reduces the left side to
  1540.  *     dl*qhat.  On the right, (nh,nm)-dh*qhat is simply the
  1541.  *     remainder r from (nh,nm)%dh, so the right is (r,nl).
  1542.  *     This produces qhat that is almost always correct and at
  1543.  *     most (prob ~ 2/2^64) one too high.
  1544.  * - Subtract qhat times the divisor (suitably shifted) from the dividend.
  1545.  *   If there is a borrow, qhat was wrong, so decrement it
  1546.  *   and add the divisor back in (once).
  1547.  * - Store the final quotient digit qhat in the quotient array q.
  1548.  *
  1549.  * Repeat the quotient digit computation for successive digits of the
  1550.  * quotient until the whole quotient has been computed.  Then shift the
  1551.  * divisor and the remainder down to correct for the normalization.
  1552.  *
  1553.  * TODO: Special case 2-word divisors.
  1554.  * TODO: Use reciprocals rather than dividing.
  1555.  */
  1556. #ifndef divn_64
  1557. BNWORD64
  1558. lbnDiv_64(BNWORD64 *q, BNWORD64 *n, unsigned nlen, BNWORD64 *d, unsigned dlen)
  1559. {
  1560.     BNWORD64 nh,nm,nl;    /* Top three words of the dividend */
  1561.     BNWORD64 dh,dl;    /* Top two words of the divisor */
  1562.     BNWORD64 qhat;    /* Extimate of quotient word */
  1563.     BNWORD64 r;    /* Remainder from quotient estimate division */
  1564.     BNWORD64 qhigh;    /* High word of quotient */
  1565.     unsigned i;    /* Temp */
  1566.     unsigned shift;    /* Bits shifted by normalization */
  1567.     unsigned qlen = nlen-dlen; /* Size of quotient (less 1) */
  1568. #ifdef mul64_ppmm
  1569.     BNWORD64 t64;
  1570. #elif defined(BNWORD128)
  1571.     BNWORD128 t128;
  1572. #else /* use lbnMulN1_64 */
  1573.     BNWORD64 t2[2];
  1574. #define t2high BIGLITTLE(t2[0],t2[1])
  1575. #define t2low BIGLITTLE(t2[1],t2[0])
  1576. #endif
  1577.  
  1578.     assert(dlen);
  1579.     assert(nlen >= dlen);
  1580.  
  1581.     /*
  1582.      * Special cases for short divisors.  The general case uses the
  1583.      * top top 2 digits of the divisor (d) to estimate a quotient digit,
  1584.      * so it breaks if there are fewer digits available.  Thus, we need
  1585.      * special cases for a divisor of length 1.  A divisor of length
  1586.      * 2 can have a *lot* of administrivia overhead removed removed,
  1587.      * so it's probably worth special-casing that case, too.
  1588.      */
  1589.     if (dlen == 1)
  1590.         return lbnDiv1_64(q, BIGLITTLE(n-1,n), n, nlen,
  1591.                           BIGLITTLE(d[-1],d[0]));
  1592.  
  1593. #if 0
  1594.     /*
  1595.      * @@@ This is not yet written...  The general loop will do,
  1596.      * albeit less efficiently
  1597.      */
  1598.     if (dlen == 2) {
  1599.         /*
  1600.          * divisor two digits long:
  1601.          * use the 3/2 technique from Knuth, but we know
  1602.          * it's exact.
  1603.          */
  1604.         dh = BIGLITTLE(d[-1],d[0]);
  1605.         dl = BIGLITTLE(d[-2],d[1]);
  1606.         shift = 0;
  1607.         if ((sh & ((BNWORD64)1 << 64-1-shift)) == 0) {
  1608.             do {
  1609.                 shift++;
  1610.             } while (dh & (BNWORD64)1<<64-1-shift) == 0);
  1611.             dh = dh << shift | dl >> (64-shift);
  1612.             dl <<= shift;
  1613.  
  1614.  
  1615.         }
  1616.  
  1617.  
  1618.         for (shift = 0; (dh & (BNWORD64)1 << 64-1-shift)) == 0; shift++)
  1619.             ;
  1620.         if (shift) {
  1621.         }
  1622.         dh = dh << shift | dl >> (64-shift);
  1623.         shift = 0;
  1624.         while (dh
  1625.     }
  1626. #endif
  1627.  
  1628.     dh = BIGLITTLE(*(d-dlen),*(d+(dlen-1)));
  1629.     assert(dh);
  1630.  
  1631.     /* Normalize the divisor */
  1632.     shift = 0;
  1633.     r = dh;
  1634.     i = 64/2;
  1635.     do {
  1636.         if (r >> i)
  1637.             r >>= i;
  1638.         else
  1639.             shift += i;
  1640.     } while ((i /= 2) != 0);
  1641.  
  1642.     nh = 0;
  1643.     if (shift) {
  1644.         lbnLshift_64(d, dlen, shift);
  1645.         dh = BIGLITTLE(*(d-dlen),*(d+(dlen-1)));
  1646.         nh = lbnLshift_64(n, nlen, shift);
  1647.     }
  1648.  
  1649.     /* Assert that dh is now normalized */
  1650.     assert(dh >> (64-1));
  1651.  
  1652.     /* Also get the second-most significant word of the divisor */
  1653.     dl = BIGLITTLE(*(d-(dlen-1)),*(d+(dlen-2)));
  1654.  
  1655.     /*
  1656.      * Adjust pointers: n to point to least significant end of first
  1657.      * first subtract, and q to one the most-significant end of the
  1658.      * quotient array.
  1659.      */
  1660.     BIGLITTLE(n -= qlen,n += qlen);
  1661.     BIGLITTLE(q -= qlen,q += qlen);
  1662.  
  1663.     /* Fetch the most significant stored word of the dividend */
  1664.     nm = BIGLITTLE(*(n-dlen),*(n+(dlen-1)));
  1665.  
  1666.     /*
  1667.      * Compute the first digit of the quotient, based on the
  1668.      * first two words of the dividend (the most significant of which
  1669.      * is the overflow word h).
  1670.      */
  1671.     if (nh) {
  1672.         assert(nh < dh);
  1673.         r = lbnDiv21_64(&qhat, nh, nm, dh);
  1674.     } else if (nm >= dh) {
  1675.         qhat = nm/dh;
  1676.         r = nm % dh;
  1677.     } else {    /* Quotient is zero */
  1678.         qhigh = 0;
  1679.         goto divloop;
  1680.     }
  1681.  
  1682.     /* Now get the third most significant word of the dividend */
  1683.     nl = BIGLITTLE(*(n-(dlen-1)),*(n+(dlen-2)));
  1684.  
  1685.     /*
  1686.      * Correct qhat, the estimate of quotient digit.
  1687.      * qhat can only be high, and at most two words high,
  1688.      * so the loop can be unrolled and abbreviated.
  1689.      */
  1690. #ifdef mul64_ppmm
  1691.     mul64_ppmm(nm, t64, qhat, dl);
  1692.     if (nm > r || (nm == r && t64 > nl)) {
  1693.         /* Decrement qhat and adjust comparison parameters */
  1694.         qhat--;
  1695.         if ((r += dh) >= dh) {
  1696.             nm -= (t64 < dl);
  1697.             t64 -= dl;
  1698.             if (nm > r || (nm == r && t64 > nl))
  1699.                 qhat--;
  1700.         }
  1701.     }
  1702. #elif defined(BNWORD128)
  1703.     t128 = (BNWORD128)qhat * dl;
  1704.     if (t128 > ((BNWORD128)r << 64) + nl) {
  1705.         /* Decrement qhat and adjust comparison parameters */
  1706.         qhat--;
  1707.         if ((r += dh) > dh) {
  1708.             t128 -= dl;
  1709.             if (t128 > ((BNWORD128)r << 64) + nl)
  1710.                 qhat--;
  1711.         }
  1712.     }
  1713. #else /* Use lbnMulN1_64 */
  1714.     lbnMulN1_64(BIGLITTLE(t2+2,t2), &dl, 1, qhat);
  1715.     if (t2high > r || (t2high == r && t2low > nl)) {
  1716.         /* Decrement qhat and adjust comparison parameters */
  1717.         qhat--;
  1718.         if ((r += dh) >= dh) {
  1719.             t2high -= (t2low < dl);
  1720.             t2low -= dl;
  1721.             if (t2high > r || (t2high == r && t2low > nl))
  1722.                 qhat--;
  1723.         }
  1724.     }
  1725. #endif
  1726.  
  1727.     /* Do the multiply and subtract */
  1728.     r = lbnMulSub1_64(n, d, dlen, qhat);
  1729.     /* If there was a borrow, add back once. */
  1730.     if (r > nh) {    /* Borrow? */
  1731.         (void)lbnAddN_64(n, d, dlen);
  1732.         qhat--;
  1733.     }
  1734.  
  1735.     /* Remember the first quotient digit. */
  1736.     qhigh = qhat;
  1737.  
  1738.     /* Now, the main division loop: */
  1739. divloop:
  1740.     while (qlen--) {
  1741.  
  1742.         /* Advance n */
  1743.         nh = BIGLITTLE(*(n-dlen),*(n+(dlen-1)));
  1744.         BIGLITTLE(++n,--n);
  1745.         nm = BIGLITTLE(*(n-dlen),*(n+(dlen-1)));
  1746.  
  1747.         if (nh == dh) {
  1748.             qhat = ~(BNWORD64)0;
  1749.             /* Optimized computation of r = (nh,nm) - qhat * dh */
  1750.             r = nh + nm;
  1751.             if (r < nh)
  1752.                 goto subtract;
  1753.         } else {
  1754.             assert(nh < dh);
  1755.             r = lbnDiv21_64(&qhat, nh, nm, dh);
  1756.         }
  1757.  
  1758.         nl = BIGLITTLE(*(n-(dlen-1)),*(n+(dlen-2)));
  1759. #ifdef mul64_ppmm
  1760.         mul64_ppmm(nm, t64, qhat, dl);
  1761.         if (nm > r || (nm == r && t64 > nl)) {
  1762.             /* Decrement qhat and adjust comparison parameters */
  1763.             qhat--;
  1764.             if ((r += dh) >= dh) {
  1765.                 nm -= (t64 < dl);
  1766.                 t64 -= dl;
  1767.                 if (nm > r || (nm == r && t64 > nl))
  1768.                     qhat--;
  1769.             }
  1770.         }
  1771. #elif defined(BNWORD128)
  1772.         t128 = (BNWORD128)qhat * dl;
  1773.         if (t128 > ((BNWORD128)r<<64) + nl) {
  1774.             /* Decrement qhat and adjust comparison parameters */
  1775.             qhat--;
  1776.             if ((r += dh) >= dh) {
  1777.                 t128 -= dl;
  1778.                 if (t128 > ((BNWORD128)r << 64) + nl)
  1779.                     qhat--;
  1780.             }
  1781.         }
  1782. #else /* Use lbnMulN1_64 */
  1783.         lbnMulN1_64(BIGLITTLE(t2+2,t2), &dl, 1, qhat);
  1784.         if (t2high > r || (t2high == r && t2low > nl)) {
  1785.             /* Decrement qhat and adjust comparison parameters */
  1786.             qhat--;
  1787.             if ((r += dh) >= dh) {
  1788.                 t2high -= (t2low < dl);
  1789.                 t2low -= dl;
  1790.                 if (t2high > r || (t2high == r && t2low > nl))
  1791.                     qhat--;
  1792.             }
  1793.         }
  1794. #endif
  1795.  
  1796.         /*
  1797.          * As a point of interest, note that it is not worth checking
  1798.          * for qhat of 0 or 1 and installing special-case code.  These
  1799.          * occur with probability 2^-64, so spending 1 cycle to check
  1800.          * for them is only worth it if we save more than 2^15 cycles,
  1801.          * and a multiply-and-subtract for numbers in the 1024-bit
  1802.          * range just doesn't take that long.
  1803.          */
  1804. subtract:
  1805.         /*
  1806.          * n points to the least significant end of the substring
  1807.          * of n to be subtracted from.  qhat is either exact or
  1808.          * one too large.  If the subtract gets a borrow, it was
  1809.          * one too large and the divisor is added back in.  It's
  1810.          * a dlen+1 word add which is guaranteed to produce a
  1811.          * carry out, so it can be done very simply.
  1812.          */
  1813.         r = lbnMulSub1_64(n, d, dlen, qhat);
  1814.         if (r > nh) {    /* Borrow? */
  1815.             (void)lbnAddN_64(n, d, dlen);
  1816.             qhat--;
  1817.         }
  1818.         /* Store the quotient digit */
  1819.         BIGLITTLE(*q++,*--q) = qhat;
  1820.     }
  1821.     /* Tah dah! */
  1822.  
  1823.     if (shift) {
  1824.         lbnRshift_64(d, dlen, shift);
  1825.         lbnRshift_64(n, dlen, shift);
  1826.     }
  1827.  
  1828.     return qhigh;
  1829. }
  1830. #endif
  1831.  
  1832. /*
  1833.  * Find the negative multiplicative inverse of x (x must be odd!) modulo 2^64.
  1834.  *
  1835.  * This just performs Newton's iteration until it gets the
  1836.  * inverse.  The initial estimate is always correct to 3 bits, and
  1837.  * sometimes 4.  The number of valid bits doubles each iteration.
  1838.  * (To prove it, assume x * y == 1 (mod 2^n), and introduce a variable
  1839.  * for the error mod 2^2n.  x * y == 1 + k*2^n (mod 2^2n) and follow
  1840.  * the iteration through.)
  1841.  */
  1842. #ifndef lbnMontInv1_64
  1843. BNWORD64
  1844. lbnMontInv1_64(BNWORD64 const x)
  1845. {
  1846.         BNWORD64 y = x, z;
  1847.  
  1848.     assert(x & 1);
  1849.  
  1850.         while ((z = x*y) != 1)
  1851.                 y *= 2 - z;
  1852.         return -y;
  1853. }
  1854. #endif /* !lbnMontInv1_64 */
  1855.  
  1856. #if defined(BNWORD128) && PRODUCT_SCAN
  1857. /*
  1858.  * Test code for product-scanning Montgomery reduction.
  1859.  * This seems to slow the C code down rather than speed it up.
  1860.  *
  1861.  * The first loop computes the Montgomery multipliers, storing them over
  1862.  * the low half of the number n.
  1863.  *
  1864.  * The second half multiplies the upper half, adding in the modulus
  1865.  * times the Montgomery multipliers.  The results of this multiply
  1866.  * are stored.
  1867.  */
  1868. void
  1869. lbnMontReduce_64(BNWORD64 *n, BNWORD64 const *mod, unsigned mlen, BNWORD64 inv)
  1870. {
  1871.     BNWORD128 x, y;
  1872.     BNWORD64 const *pm;
  1873.     BNWORD64 *pn;
  1874.     BNWORD64 t;
  1875.     unsigned carry;
  1876.     unsigned i, j;
  1877.  
  1878.     /* Special case of zero */
  1879.     if (!mlen)
  1880.         return;
  1881.  
  1882.     /* Pass 1 - compute Montgomery multipliers */
  1883.     /* First iteration can have certain simplifications. */
  1884.     t = BIGLITTLE(n[-1],n[0]);
  1885.     x = t;
  1886.     t *= inv;
  1887.     BIGLITTLE(n[-1], n[0]) = t;
  1888.     x += (BNWORD128)t * BIGLITTLE(mod[-1],mod[0]); /* Can't overflow */
  1889.     assert((BNWORD64)x == 0);
  1890.     x = x >> 64;
  1891.  
  1892.     for (i = 1; i < mlen; i++) {
  1893.         carry = 0;
  1894.         pn = n;
  1895.         pm = BIGLITTLE(mod-i-1,mod+i+1);
  1896.         for (j = 0; j < i; j++) {
  1897.             y = (BNWORD128)BIGLITTLE(*--pn * *pm++, *pn++ * *--pm);
  1898.             x += y;
  1899.             carry += (x < y);
  1900.         }
  1901.         assert(BIGLITTLE(pn == n-i, pn == n+i));
  1902.         y = t = BIGLITTLE(pn[-1], pn[0]);
  1903.         x += y;
  1904.         carry += (x < y);
  1905.         BIGLITTLE(pn[-1], pn[0]) = t = inv * (BNWORD64)x;
  1906.         assert(BIGLITTLE(pm == mod-1, pm == mod+1));
  1907.         y = (BNWORD128)t * BIGLITTLE(pm[0],pm[-1]);
  1908.         x += y;
  1909.         carry += (x < y);
  1910.         assert((BNWORD64)x == 0);
  1911.         x = x >> 64 | (BNWORD128)carry << 64;
  1912.     }
  1913.  
  1914.     BIGLITTLE(n -= mlen, n += mlen);
  1915.  
  1916.     /* Pass 2 - compute upper words and add to n */
  1917.     for (i = 1; i < mlen; i++) {
  1918.         carry = 0;
  1919.         pm = BIGLITTLE(mod-i,mod+i);
  1920.         pn = n;
  1921.         for (j = i; j < mlen; j++) {
  1922.             y = (BNWORD128)BIGLITTLE(*--pm * *pn++, *pm++ * *--pn);
  1923.             x += y;
  1924.             carry += (x < y);
  1925.         }
  1926.         assert(BIGLITTLE(pm == mod-mlen, pm == mod+mlen));
  1927.         assert(BIGLITTLE(pn == n+mlen-i, pn == n-mlen+i));
  1928.         y = t = BIGLITTLE(*(n-i),*(n+i-1));
  1929.         x += y;
  1930.         carry += (x < y);
  1931.         BIGLITTLE(*(n-i),*(n+i-1)) = (BNWORD64)x;
  1932.         x = (x >> 64) | (BNWORD128)carry << 64;
  1933.     }
  1934.  
  1935.     /* Last round of second half, simplified. */
  1936.     t = BIGLITTLE(*(n-mlen),*(n+mlen-1));
  1937.     x += t;
  1938.     BIGLITTLE(*(n-mlen),*(n+mlen-1)) = (BNWORD64)x;
  1939.     carry = (unsigned)(x >> 64);
  1940.  
  1941.     while (carry)
  1942.         carry -= lbnSubN_64(n, mod, mlen);
  1943.     while (lbnCmp_64(n, mod, mlen) >= 0)
  1944.         (void)lbnSubN_64(n, mod, mlen);
  1945. }
  1946. #define lbnMontReduce_64 lbnMontReduce_64
  1947. #endif
  1948.  
  1949. /*
  1950.  * Montgomery reduce n, modulo mod.  This reduces modulo mod and divides by
  1951.  * 2^(64*mlen).  Returns the result in the *top* mlen words of the argument n.
  1952.  * This is ready for another multiplication using lbnMul_64.
  1953.  *
  1954.  * Montgomery representation is a very useful way to encode numbers when
  1955.  * you're doing lots of modular reduction.  What you do is pick a multiplier
  1956.  * R which is relatively prime to the modulus and very easy to divide by.
  1957.  * Since the modulus is odd, R is closen as a power of 2, so the division
  1958.  * is a shift.  In fact, it's a shift of an integral number of words,
  1959.  * so the shift can be implicit - just drop the low-order words.
  1960.  *
  1961.  * Now, choose R *larger* than the modulus m, 2^(64*mlen).  Then convert
  1962.  * all numbers a, b, etc. to Montgomery form M(a), M(b), etc using the
  1963.  * relationship M(a) = a*R mod m, M(b) = b*R mod m, etc.  Note that:
  1964.  * - The Montgomery form of a number depends on the modulus m.
  1965.  *   A fixed modulus m is assumed throughout this discussion.
  1966.  * - Since R is relaitvely prime to m, multiplication by R is invertible;
  1967.  *   no information about the numbers is lost, they're just scrambled.
  1968.  * - Adding (and subtracting) numbers in this form works just as usual.
  1969.  *   M(a+b) = (a+b)*R mod m = (a*R + b*R) mod m = (M(a) + M(b)) mod m
  1970.  * - Multiplying numbers in this form produces a*b*R*R.  The problem
  1971.  *   is to divide out the excess factor of R, modulo m as well as to
  1972.  *   reduce to the given length mlen.  It turns out that this can be
  1973.  *   done *faster* than a normal divide, which is where the speedup
  1974.  *   in Montgomery division comes from.
  1975.  *
  1976.  * Normal reduction chooses a most-significant quotient digit q and then
  1977.  * subtracts q*m from the number to be reduced.  Choosing q is tricky
  1978.  * and involved (just look at lbnDiv_64 to see!) and is usually
  1979.  * imperfect, requiring a check for correction after the subtraction.
  1980.  *
  1981.  * Montgomery reduction *adds* a multiple of m to the *low-order* part
  1982.  * of the number to be reduced.  This multiple is chosen to make the
  1983.  * low-order part of the number come out to zero.  This can be done
  1984.  * with no trickery or error using a precomputed inverse of the modulus.
  1985.  * In this code, the "part" is one word, but any width can be used.
  1986.  *
  1987.  * Repeating this step sufficiently often results in a value which
  1988.  * is a multiple of R (a power of two, remember) but is still (since
  1989.  * the additions were to the low-order part and thus did not increase
  1990.  * the value of the number being reduced very much) still not much
  1991.  * larger than m*R.  Then implicitly divide by R and subtract off
  1992.  * m until the result is in the correct range.
  1993.  *
  1994.  * Since the low-order part being cancelled is less than R, the
  1995.  * multiple of m added must have a multiplier which is at most R-1.
  1996.  * Assuming that the input is at most m*R-1, the final number is
  1997.  * at most m*(2*R-1)-1 = 2*m*R - m - 1, so subtracting m once from
  1998.  * the high-order part, equivalent to subtracting m*R from the
  1999.  * while number, produces a result which is at most m*R - m - 1,
  2000.  * which divided by R is at most m-1.
  2001.  *
  2002.  * To convert *to* Montgomery form, you need a regular remainder
  2003.  * routine, although you can just compute R*R (mod m) and do the
  2004.  * conversion using Montgomery multiplication.  To convert *from*
  2005.  * Montgomery form, just Montgomery reduce the number to
  2006.  * remove the extra factor of R.
  2007.  *
  2008.  * TODO: Change to a full inverse and use Karatsuba's multiplication
  2009.  * rather than this word-at-a-time.
  2010.  */
  2011. #ifndef lbnMontReduce_64
  2012. void
  2013. lbnMontReduce_64(BNWORD64 *n, BNWORD64 const *mod, unsigned const mlen,
  2014.                 BNWORD64 inv)
  2015. {
  2016.     BNWORD64 t;
  2017.     BNWORD64 c = 0;
  2018.     unsigned len = mlen;
  2019.  
  2020.     /* inv must be the negative inverse of mod's least significant word */
  2021.     assert((BNWORD64)(inv * BIGLITTLE(mod[-1],mod[0])) == (BNWORD64)-1);
  2022.  
  2023.     assert(len);
  2024.  
  2025.     do {
  2026.         t = lbnMulAdd1_64(n, mod, mlen, inv * BIGLITTLE(n[-1],n[0]));
  2027.         c += lbnAdd1_64(BIGLITTLE(n-mlen,n+mlen), len, t);
  2028.         BIGLITTLE(--n,++n);
  2029.     } while (--len);
  2030.  
  2031.     /*
  2032.      * All that adding can cause an overflow past the modulus size,
  2033.      * but it's unusual, and never by much, so a subtraction loop
  2034.      * is the right way to deal with it.
  2035.      * This subtraction happens infrequently - I've only ever seen it
  2036.      * invoked once per reduction, and then just under 22.5% of the time.
  2037.      */
  2038.     while (c)
  2039.         c -= lbnSubN_64(n, mod, mlen);
  2040.     while (lbnCmp_64(n, mod, mlen) >= 0)
  2041.         (void)lbnSubN_64(n, mod, mlen);
  2042. }
  2043. #endif /* !lbnMontReduce_64 */
  2044.  
  2045. /*
  2046.  * A couple of helpers that you might want to implement atomically
  2047.  * in asm sometime.
  2048.  */
  2049. #ifndef lbnMontMul_64
  2050. /*
  2051.  * Multiply "num1" by "num2", modulo "mod", all of length "len", and
  2052.  * place the result in the high half of "prod".  "inv" is the inverse
  2053.  * of the least-significant word of the modulus, modulo 2^64.
  2054.  * This uses numbers in Montgomery form.  Reduce using "len" and "inv".
  2055.  *
  2056.  * This is implemented as a macro to win on compilers that don't do
  2057.  * inlining, since it's so trivial.
  2058.  */
  2059. #define lbnMontMul_64(prod, n1, n2, mod, len, inv) \
  2060.     (lbnMulX_64(prod, n1, n2, len), lbnMontReduce_64(prod, mod, len, inv))
  2061. #endif /* !lbnMontMul_64 */
  2062.  
  2063. #ifndef lbnMontSquare_64
  2064. /*
  2065.  * Square "num", modulo "mod", both of length "len", and place the result
  2066.  * in the high half of "prod".  "inv" is the inverse of the least-significant
  2067.  * word of the modulus, modulo 2^64.
  2068.  * This uses numbers in Montgomery form.  Reduce using "len" and "inv".
  2069.  *
  2070.  * This is implemented as a macro to win on compilers that don't do
  2071.  * inlining, since it's so trivial.
  2072.  */
  2073. #define lbnMontSquare_64(prod, n, mod, len, inv) \
  2074.     (lbnSquare_64(prod, n, len), lbnMontReduce_64(prod, mod, len, inv))
  2075.     
  2076. #endif /* !lbnMontSquare_64 */
  2077.  
  2078. /*
  2079.  * Convert a number to Montgomery form - requires mlen + nlen words
  2080.  * of memory in "n".
  2081.  */
  2082. void
  2083. lbnToMont_64(BNWORD64 *n, unsigned nlen, BNWORD64 *mod, unsigned mlen)
  2084. {
  2085.     /* Move n up "mlen" words */
  2086.     lbnCopy_64(BIGLITTLE(n-mlen,n+mlen), n, nlen);
  2087.     lbnZero_64(n, mlen);
  2088.     /* Do the division - dump the quotient in the high-order words */
  2089.     (void)lbnDiv_64(BIGLITTLE(n-mlen,n+mlen), n, mlen+nlen, mod, mlen);
  2090. }
  2091.  
  2092. /*
  2093.  * Convert from Montgomery form.  Montgomery reduction is all that is
  2094.  * needed.
  2095.  */
  2096. void
  2097. lbnFromMont_64(BNWORD64 *n, BNWORD64 *mod, unsigned len)
  2098. {
  2099.     /* Zero the high words of n */
  2100.     lbnZero_64(BIGLITTLE(n-len,n+len), len);
  2101.     lbnMontReduce_64(n, mod, len, lbnMontInv1_64(BIGLITTLE(mod[-1],mod[0])));
  2102.     /* Move n down len words */
  2103.     lbnCopy_64(n, BIGLITTLE(n-len,n+len), len);
  2104. }
  2105.  
  2106. /*
  2107.  * The windowed exponentiation algorithm, precomputes a table of odd
  2108.  * powers of n up to 2^k.  It takes 2^(k-1)-1 multiplies to compute
  2109.  * the table, and (e-1)/(k+1) multiplies (on average) to perform the
  2110.  * exponentiation.  To minimize the sum, k must vary with e.
  2111.  * The optimal window sizes vary with the exponent length.  Here are
  2112.  * some selected values and the boundary cases.
  2113.  * (An underscore _ has been inserted into some of the numbers to ensure
  2114.  * that magic strings like 64 do not appear in this table.  It should be
  2115.  * ignored.)
  2116.  *
  2117.  * At e =    1 bits, k=1   (0.000000) is best.
  2118.  * At e =    2 bits, k=1   (0.500000) is best.
  2119.  * At e =    4 bits, k=1   (1.500000) is best.
  2120.  * At e =    8 bits, k=2   (3.333333) < k=1   (3.500000)
  2121.  * At e =  1_6 bits, k=2   (6.000000) is best.
  2122.  * At e =   26 bits, k=3   (9.250000) < k=2   (9.333333)
  2123.  * At e =  3_2 bits, k=3  (10.750000) is best.
  2124.  * At e =  6_4 bits, k=3  (18.750000) is best.
  2125.  * At e =   82 bits, k=4  (23.200000) < k=3  (23.250000)
  2126.  * At e =  128 bits, k=4 (3_2.400000) is best.
  2127.  * At e =  242 bits, k=5  (55.1_66667) < k=4 (55.200000)
  2128.  * At e =  256 bits, k=5  (57.500000) is best.
  2129.  * At e =  512 bits, k=5 (100.1_66667) is best.
  2130.  * At e =  674 bits, k=6 (127.142857) < k=5 (127.1_66667)
  2131.  * At e = 1024 bits, k=6 (177.142857) is best.
  2132.  * At e = 1794 bits, k=7 (287.125000) < k=6 (287.142857)
  2133.  * At e = 2048 bits, k=7 (318.875000) is best.
  2134.  * At e = 4096 bits, k=7 (574.875000) is best.
  2135.  *
  2136.  * The numbers in parentheses are the expected number of multiplications
  2137.  * needed to do the computation.  The normal russian-peasant modular
  2138.  * exponentiation technique always uses (e-1)/2.  For exponents as
  2139.  * small as 192 bits (below the range of current factoring algorithms),
  2140.  * half of the multiplies are eliminated, 45.2 as opposed to the naive
  2141.  * 95.5.  Counting the 191 squarings as 3/4 a multiply each (squaring
  2142.  * proper is just over half of multiplying, but the Montgomery
  2143.  * reduction in each case is also a multiply), that's 143.25
  2144.  * multiplies, for totals of 188.45 vs. 238.75 - a 21% savings.
  2145.  * For larger exponents (like 512 bits), it's 483.92 vs. 639.25, a
  2146.  * 24.3% savings.  It asymptotically approaches 25%.
  2147.  *
  2148.  * Given that exponents for which k>7 are useful are uncommon,
  2149.  * a fixed size table for k <= 7 is used for simplicity.
  2150.  * k = 8 is uzeful at 4610 bits, k = 9 at 11522 bits.
  2151.  *
  2152.  * The basic number of squarings needed is e-1, although a k-bit
  2153.  * window (for k > 1) can save, on average, k-2 of those, too.
  2154.  * That savings currently isn't counted here.  It would drive the
  2155.  * crossover points slightly lower.
  2156.  * (Actually, this win is also reduced in the DoubleExpMod case,
  2157.  * meaning we'd have to split the tables.  Except for that, the
  2158.  * multiplies by powers of the two bases are independent, so
  2159.  * the same logic applies to each as the single case.)
  2160.  *
  2161.  * Table entry i is the largest number of bits in an exponent to
  2162.  * process with a window size of i+1.  So the window never goes above 7
  2163.  * bits, requiring 2^(7-1) = 0x40 precomputed multiples.
  2164.  */
  2165. #define BNEXPMOD_MAX_WINDOW    7
  2166. static unsigned const bnExpModThreshTable[BNEXPMOD_MAX_WINDOW] = {
  2167.     7, 25, 81, 241, 673, 1793, (unsigned)-1
  2168. };
  2169.  
  2170. /*
  2171.  * Perform modular exponentiation, as fast as possible!  This uses
  2172.  * Montgomery reduction, optimized squaring, and windowed exponentiation.
  2173.  * The modulus "mod" MUST be odd!
  2174.  *
  2175.  * This returns 0 on success, -1 on out of memory.
  2176.  *
  2177.  * The window algorithm:
  2178.  * The idea is to keep a running product of b1 = n^(high-order bits of exp),
  2179.  * and then keep appending exponent bits to it.  The following patterns
  2180.  * apply to a 3-bit window (k = 3):
  2181.  * To append   0: square
  2182.  * To append   1: square, multiply by n^1
  2183.  * To append  10: square, multiply by n^1, square
  2184.  * To append  11: square, square, multiply by n^3
  2185.  * To append 100: square, multiply by n^1, square, square
  2186.  * To append 101: square, square, square, multiply by n^5
  2187.  * To append 110: square, square, multiply by n^3, square
  2188.  * To append 111: square, square, square, multiply by n^7
  2189.  *
  2190.  * Since each pattern involves only one multiply, the longer the pattern
  2191.  * the better, except that a 0 (no multiplies) can be appended directly.
  2192.  * We precompute a table of odd powers of n, up to 2^k, and can then
  2193.  * multiply k bits of exponent at a time.  Actually, assuming random
  2194.  * exponents, there is on average one zero bit between needs to
  2195.  * multiply (1/2 of the time there's none, 1/4 of the time there's 1,
  2196.  * 1/8 of the time, there's 2, 1/64 of the time, there's 3, etc.), so
  2197.  * you have to do one multiply per k+1 bits of exponent.
  2198.  *
  2199.  * The loop walks down the exponent, squaring the result buffer as
  2200.  * it goes.  There is a wbits+1 bit lookahead buffer, buf, that is
  2201.  * filled with the upcoming exponent bits.  (What is read after the
  2202.  * end of the exponent is unimportant, but it is filled with zero here.)
  2203.  * When the most-significant bit of this buffer becomes set, i.e.
  2204.  * (buf & tblmask) != 0, we have to decide what pattern to multiply
  2205.  * by, and when to do it.  We decide, remember to do it in future
  2206.  * after a suitable number of squarings have passed (e.g. a pattern
  2207.  * of "100" in the buffer requires that we multiply by n^1 immediately;
  2208.  * a pattern of "110" calls for multiplying by n^3 after one more
  2209.  * squaring), clear the buffer, and continue.
  2210.  *
  2211.  * When we start, there is one more optimization: the result buffer
  2212.  * is implcitly one, so squaring it or multiplying by it can be
  2213.  * optimized away.  Further, if we start with a pattern like "100"
  2214.  * in the lookahead window, rather than placing n into the buffer
  2215.  * and then starting to square it, we have already computed n^2
  2216.  * to compute the odd-powers table, so we can place that into
  2217.  * the buffer and save a squaring.
  2218.  *
  2219.  * This means that if you have a k-bit window, to compute n^z,
  2220.  * where z is the high k bits of the exponent, 1/2 of the time
  2221.  * it requires no squarings.  1/4 of the time, it requires 1
  2222.  * squaring, ... 1/2^(k-1) of the time, it reqires k-2 squarings.
  2223.  * And the remaining 1/2^(k-1) of the time, the top k bits are a
  2224.  * 1 followed by k-1 0 bits, so it again only requires k-2
  2225.  * squarings, not k-1.  The average of these is 1.  Add that
  2226.  * to the one squaring we have to do to compute the table,
  2227.  * and you'll see that a k-bit window saves k-2 squarings
  2228.  * as well as reducing the multiplies.  (It actually doesn't
  2229.  * hurt in the case k = 1, either.)
  2230.  *
  2231.  * n must have mlen words allocated.  Although fewer may be in use
  2232.  * when n is passed in, all are in use on exit.
  2233.  */
  2234. int
  2235. lbnExpMod_64(BNWORD64 *result, BNWORD64 const *n, unsigned nlen,
  2236.     BNWORD64 const *e, unsigned elen, BNWORD64 *mod, unsigned mlen)
  2237. {
  2238.     BNWORD64 *table[1 << (BNEXPMOD_MAX_WINDOW-1)];
  2239.                 /* Table of odd powers of n */
  2240.     unsigned ebits;        /* Exponent bits */
  2241.     unsigned wbits;        /* Window size */
  2242.     unsigned tblmask;    /* Mask of exponentiation window */
  2243.     BNWORD64 bitpos;    /* Mask of current look-ahead bit */
  2244.     unsigned buf;        /* Buffer of exponent bits */
  2245.     unsigned multpos;    /* Where to do pending multiply */
  2246.     BNWORD64 const *mult;    /* What to multiply by */
  2247.     unsigned i;        /* Loop counter */
  2248.     int isone;        /* Flag: accum. is implicitly one */
  2249.     BNWORD64 *a, *b;    /* Working buffers/accumulators */
  2250.     BNWORD64 *t;        /* Pointer into the working buffers */
  2251.     BNWORD64 inv;        /* mod^-1 modulo 2^64 */
  2252.  
  2253.     assert(mlen);
  2254.     assert(nlen <= mlen);
  2255.  
  2256.     /* First, a couple of trivial cases. */
  2257.     elen = lbnNorm_64(e, elen);
  2258.     if (!elen) {
  2259.         /* x ^ 0 == 1 */
  2260.         lbnZero_64(result, mlen);
  2261.         BIGLITTLE(result[-1],result[0]) = 1;
  2262.         return 0;
  2263.     }
  2264.     ebits = lbnBits_64(e, elen);
  2265.     if (ebits == 1) {
  2266.         /* x ^ 1 == x */
  2267.         if (n != result)
  2268.             lbnCopy_64(result, n, nlen);
  2269.         if (mlen > nlen)
  2270.             lbnZero_64(BIGLITTLE(result-nlen,result+nlen),
  2271.                        mlen-nlen);
  2272.         return 0;
  2273.     }
  2274.  
  2275.     /* Okay, now move the exponent pointer to the most-significant word */
  2276.     e = BIGLITTLE(e-elen, e+elen-1);
  2277.  
  2278.     /* Look up appropriate k-1 for the exponent - tblmask = 1<<(k-1) */
  2279.     wbits = 0;
  2280.     while (ebits > bnExpModThreshTable[wbits])
  2281.         wbits++;
  2282.  
  2283.     /* Allocate working storage: two product buffers and the tables. */
  2284.     LBNALLOC(a, 2*mlen);
  2285.     if (!a)
  2286.         return -1;
  2287.     LBNALLOC(b, 2*mlen);
  2288.     if (!b) {
  2289.         LBNFREE(a, 2*mlen);
  2290.         return -1;
  2291.     }
  2292.  
  2293.     /* Convert to the appropriate table size: tblmask = 1<<(k-1) */
  2294.     tblmask = 1u << wbits;
  2295.  
  2296.     LBNALLOC(t, mlen);
  2297.     if (!t) {
  2298.         LBNFREE(b, 2*mlen);
  2299.         LBNFREE(a, 2*mlen);
  2300.         return -1;
  2301.     }
  2302.     /* We have the result buffer available, so use it. */
  2303.     table[0] = result;
  2304.  
  2305.     /*
  2306.      * Okay, we now have a minimal-sized table - expand it.
  2307.      * This is allowed to fail!  If so, scale back the table size
  2308.      * and proceed.
  2309.      */
  2310.     for (i = 1; i < tblmask; i++) {
  2311.         LBNALLOC(t, mlen);
  2312.         if (!t)    /* Out of memory!  Quit the loop. */
  2313.             break;
  2314.         table[i] = t;
  2315.     }
  2316.  
  2317.     /* If we stopped, with i < tblmask, shrink the tables appropriately */
  2318.     while (tblmask > i) {
  2319.         wbits--;
  2320.         tblmask >>= 1;
  2321.     }
  2322.     /* Free up our overallocations */
  2323.     while (--i > tblmask)
  2324.         LBNFREE(table[i], mlen);
  2325.  
  2326.     /* Okay, fill in the table */
  2327.  
  2328.     /* Compute the necessary modular inverse */
  2329.     inv = lbnMontInv1_64(mod[BIGLITTLE(-1,0)]);    /* LSW of modulus */
  2330.  
  2331.     /* Convert n to Montgomery form */
  2332.  
  2333.     /* Move n up "mlen" words into a */
  2334.     t = BIGLITTLE(a-mlen, a+mlen);
  2335.     lbnCopy_64(t, n, nlen);
  2336.     lbnZero_64(a, mlen);
  2337.     /* Do the division - lose the quotient into the high-order words */
  2338.     (void)lbnDiv_64(t, a, mlen+nlen, mod, mlen);
  2339.     /* Copy into first table entry */
  2340.     lbnCopy_64(table[0], a, mlen);
  2341.  
  2342.     /* Square a into b */
  2343.     lbnMontSquare_64(b, a, mod, mlen, inv);
  2344.  
  2345.     /* Use high half of b to initialize the table */
  2346.     t = BIGLITTLE(b-mlen, b+mlen);
  2347.     for (i = 1; i < tblmask; i++) {
  2348.         lbnMontMul_64(a, t, table[i-1], mod, mlen, inv);
  2349.         lbnCopy_64(table[i], BIGLITTLE(a-mlen, a+mlen), mlen);
  2350.     }
  2351.  
  2352.     /* We might use b = n^2 later... */
  2353.  
  2354.     /* Initialze the fetch pointer */
  2355.     bitpos = (BNWORD64)1 << ((ebits-1) & (64-1));    /* Initialize mask */
  2356.  
  2357.     /* This should point to the msbit of e */
  2358.     assert((*e & bitpos) != 0);
  2359.  
  2360.     /*
  2361.      * Pre-load the window.  Becuase the window size is
  2362.      * never larger than the exponent size, there is no need to
  2363.      * detect running off the end of e in here.
  2364.      *
  2365.      * The read-ahead is controlled by elen and the bitpos mask.
  2366.      * Note that this is *ahead* of ebits, which tracks the
  2367.      * most significant end of the window.  The purpose of this
  2368.      * initialization is to get the two wbits+1 bits apart,
  2369.      * like they should be.
  2370.      *
  2371.      * Note that bitpos and e1len together keep track of the
  2372.      * lookahead read pointer in the exponent that is used here.
  2373.      */
  2374.     buf = 0;
  2375.     for (i = 0; i <= wbits; i++) {
  2376.         buf = (buf << 1) | ((*e & bitpos) != 0);
  2377.         bitpos >>= 1;
  2378.         if (!bitpos) {
  2379.             BIGLITTLE(e++,e--);
  2380.             bitpos = (BNWORD64)1 << (64-1);
  2381.             elen--;
  2382.         }
  2383.     }
  2384.     assert(buf & tblmask);
  2385.  
  2386.     /*
  2387.      * Set the pending multiply positions to a location that will
  2388.      * never be encountered, thus ensuring that nothing will happen
  2389.      * until the need for a multiply appears and one is scheduled.
  2390.      */
  2391.     multpos = ebits;    /* A NULL value */
  2392.     mult = 0;    /* Force a crash if we use these */
  2393.  
  2394.     /*
  2395.      * Okay, now begins the real work.  The first step is
  2396.      * slightly magic, so it's done outside the main loop,
  2397.      * but it's very similar to what's inside.
  2398.      */
  2399.     ebits--;    /* Start processing the first bit... */
  2400.     isone = 1;
  2401.  
  2402.     /*
  2403.      * This is just like the multiply in the loop, except that
  2404.      * - We know the msbit of buf is set, and
  2405.      * - We have the extra value n^2 floating around.
  2406.      * So, do the usual computation, and if the result is that
  2407.      * the buffer should be multiplied by n^1 immediately
  2408.      * (which we'd normally then square), we multiply it
  2409.      * (which reduces to a copy, which reduces to setting a flag)
  2410.      * by n^2 and skip the squaring.  Thus, we do the
  2411.      * multiply and the squaring in one step.
  2412.      */
  2413.     assert(buf & tblmask);
  2414.     multpos = ebits - wbits;
  2415.     while ((buf & 1) == 0) {
  2416.         buf >>= 1;
  2417.         multpos++;
  2418.     }
  2419.     /* Intermediates can wrap, but final must NOT */
  2420.     assert(multpos <= ebits);
  2421.     mult = table[buf>>1];
  2422.     buf = 0;
  2423.  
  2424.     /* Special case: use already-computed value sitting in buffer */
  2425.     if (multpos == ebits)
  2426.         isone = 0;
  2427.  
  2428.     /*
  2429.      * At this point, the buffer (which is the high half of b) holds
  2430.      * either 1 (implicitly, as the "isone" flag is set), or n^2.
  2431.      */
  2432.  
  2433.     /*
  2434.      * The main loop.  The procedure is:
  2435.      * - Advance the window
  2436.      * - If the most-significant bit of the window is set,
  2437.      *   schedule a multiply for the appropriate time in the
  2438.      *   future (may be immediately)
  2439.      * - Perform any pending multiples
  2440.      * - Check for termination
  2441.      * - Square the buffer
  2442.      *
  2443.      * At any given time, the acumulated product is held in
  2444.      * the high half of b.
  2445.      */
  2446.     for (;;) {
  2447.         ebits--;
  2448.  
  2449.         /* Advance the window */
  2450.         assert(buf < tblmask);
  2451.         buf <<= 1;
  2452.         /*
  2453.          * This reads ahead of the current exponent position
  2454.          * (controlled by ebits), so we have to be able to read
  2455.          * past the lsb of the exponents without error.
  2456.          */
  2457.         if (elen) {
  2458.             buf |= ((*e & bitpos) != 0);
  2459.             bitpos >>= 1;
  2460.             if (!bitpos) {
  2461.                 BIGLITTLE(e++,e--);
  2462.                 bitpos = (BNWORD64)1 << (64-1);
  2463.                 elen--;
  2464.             }
  2465.         }
  2466.  
  2467.         /* Examine the window for pending multiplies */
  2468.         if (buf & tblmask) {
  2469.             multpos = ebits - wbits;
  2470.             while ((buf & 1) == 0) {
  2471.                 buf >>= 1;
  2472.                 multpos++;
  2473.             }
  2474.             /* Intermediates can wrap, but final must NOT */
  2475.             assert(multpos <= ebits);
  2476.             mult = table[buf>>1];
  2477.             buf = 0;
  2478.         }
  2479.  
  2480.         /* If we have a pending multiply, do it */
  2481.         if (ebits == multpos) {
  2482.             /* Multiply by the table entry remembered previously */
  2483.             t = BIGLITTLE(b-mlen, b+mlen);
  2484.             if (isone) {
  2485.                 /* Multiply by 1 is a trivial case */
  2486.                 lbnCopy_64(t, mult, mlen);
  2487.                 isone = 0;
  2488.             } else {
  2489.                 lbnMontMul_64(a, t, mult, mod, mlen, inv);
  2490.                 /* Swap a and b */
  2491.                 t = a; a = b; b = t;
  2492.             }
  2493.         }
  2494.  
  2495.         /* Are we done? */
  2496.         if (!ebits)
  2497.             break;
  2498.  
  2499.         /* Square the input */
  2500.         if (!isone) {
  2501.             t = BIGLITTLE(b-mlen, b+mlen);
  2502.             lbnMontSquare_64(a, t, mod, mlen, inv);
  2503.             /* Swap a and b */
  2504.             t = a; a = b; b = t;
  2505.         }
  2506.     } /* for (;;) */
  2507.  
  2508.     assert(!isone);
  2509.     assert(!buf);
  2510.  
  2511.     /* DONE! */
  2512.  
  2513.     /* Convert result out of Montgomery form */
  2514.     t = BIGLITTLE(b-mlen, b+mlen);
  2515.     lbnCopy_64(b, t, mlen);
  2516.     lbnZero_64(t, mlen);
  2517.     lbnMontReduce_64(b, mod, mlen, inv);
  2518.     lbnCopy_64(result, t, mlen);
  2519.     /*
  2520.      * Clean up - free intermediate storage.
  2521.      * Do NOT free table[0], which is the result
  2522.      * buffer.
  2523.      */
  2524.     while (--tblmask)
  2525.         LBNFREE(table[tblmask], mlen);
  2526.     LBNFREE(b, 2*mlen);
  2527.     LBNFREE(a, 2*mlen);
  2528.  
  2529.     return 0;    /* Success */
  2530. }
  2531.  
  2532. /*
  2533.  * Compute and return n1^e1 * n2^e2 mod "mod".
  2534.  * result may be either input buffer, or something separate.
  2535.  * It must be "mlen" words long.
  2536.  *
  2537.  * There is a current position in the exponents, which is kept in e1bits.
  2538.  * (The exponents are swapped if necessary so e1 is the longer of the two.)
  2539.  * At any given time, the value in the accumulator is
  2540.  * n1^(e1>>e1bits) * n2^(e2>>e1bits) mod "mod".
  2541.  * As e1bits is counted down, this is updated, by squaring it and doing
  2542.  * any necessary multiplies.
  2543.  * To decide on the necessary multiplies, two windows, each w1bits+1 bits
  2544.  * wide, are maintained in buf1 and buf2, which read *ahead* of the
  2545.  * e1bits position (with appropriate handling of the case when e1bits
  2546.  * drops below w1bits+1).  When the most-significant bit of either window
  2547.  * becomes set, indicating that something needs to be multiplied by
  2548.  * the accumulator or it will get out of sync, the window is examined
  2549.  * to see which power of n1 or n2 to multiply by, and when (possibly
  2550.  * later, if the power is greater than 1) the multiply should take
  2551.  * place.  Then the multiply and its location are remembered and the
  2552.  * window is cleared.
  2553.  *
  2554.  * If we had every power of n1 in the table, the multiply would always
  2555.  * be w1bits steps in the future.  But we only keep the odd powers,
  2556.  * so instead of waiting w1bits squarings and then multiplying
  2557.  * by n1^k, we wait w1bits-k squarings and multiply by n1.
  2558.  *
  2559.  * Actually, w2bits can be less than w1bits, but the window is the same
  2560.  * size, to make it easier to keep track of where we're reading.  The
  2561.  * appropriate number of low-order bits of the window are just ignored.
  2562.  */
  2563. int
  2564. lbnDoubleExpMod_64(BNWORD64 *result,
  2565.                    BNWORD64 const *n1, unsigned n1len,
  2566.                    BNWORD64 const *e1, unsigned e1len,
  2567.                    BNWORD64 const *n2, unsigned n2len,
  2568.                    BNWORD64 const *e2, unsigned e2len,
  2569.                    BNWORD64 *mod, unsigned mlen)
  2570. {
  2571.     BNWORD64 *table1[1 << (BNEXPMOD_MAX_WINDOW-1)];
  2572.                     /* Table of odd powers of n1 */
  2573.     BNWORD64 *table2[1 << (BNEXPMOD_MAX_WINDOW-1)];
  2574.                     /* Table of odd powers of n2 */
  2575.     unsigned e1bits, e2bits;    /* Exponent bits */
  2576.     unsigned w1bits, w2bits;    /* Window sizes */
  2577.     unsigned tblmask;        /* Mask of exponentiation window */
  2578.     BNWORD64 bitpos;        /* Mask of current look-ahead bit */
  2579.     unsigned buf1, buf2;        /* Buffer of exponent bits */
  2580.     unsigned mult1pos, mult2pos;    /* Where to do pending multiply */
  2581.     BNWORD64 const *mult1, *mult2;    /* What to multiply by */
  2582.     unsigned i;            /* Loop counter */
  2583.     int isone;            /* Flag: accum. is implicitly one */
  2584.     BNWORD64 *a, *b;        /* Working buffers/accumulators */
  2585.     BNWORD64 *t;            /* Pointer into the working buffers */
  2586.     BNWORD64 inv;            /* mod^-1 modulo 2^64 */
  2587.  
  2588.     assert(mlen);
  2589.     assert(n1len <= mlen);
  2590.     assert(n2len <= mlen);
  2591.  
  2592.     /* First, a couple of trivial cases. */
  2593.     e1len = lbnNorm_64(e1, e1len);
  2594.     e2len = lbnNorm_64(e2, e2len);
  2595.  
  2596.     /* Ensure that the first exponent is the longer */
  2597.     e1bits = lbnBits_64(e1, e1len);
  2598.     e2bits = lbnBits_64(e2, e2len);
  2599.     if (e1bits < e2bits) {
  2600.         i = e1len; e1len = e2len; e2len = i;
  2601.         i = e1bits; e1bits = e2bits; e2bits = i;
  2602.         t = (BNWORD64 *)n1; n1 = n2; n2 = t;
  2603.         t = (BNWORD64 *)e1; e1 = e2; e2 = t;
  2604.     }
  2605.     assert(e1bits >= e2bits);
  2606.  
  2607.     /* Handle a trivial case */
  2608.     if (!e2len)
  2609.         return lbnExpMod_64(result, n1, n1len, e1, e1len, mod, mlen);
  2610.     assert(e2bits);
  2611.  
  2612.     /* The code below fucks up if the exponents aren't at least 2 bits */
  2613.     if (e1bits == 1) {
  2614.         assert(e2bits == 1);
  2615.  
  2616.         LBNALLOC(a, n1len+n2len);
  2617.         if (!a)
  2618.             return -1;
  2619.  
  2620.         lbnMul_64(a, n1, n1len, n2, n2len);
  2621.         /* Do a direct modular reduction */
  2622.         if (n1len + n2len >= mlen)
  2623.             (void)lbnDiv_64(a+mlen, a, n1len+n2len, mod, mlen);
  2624.         lbnCopy_64(result, a, mlen);
  2625.         LBNFREE(a, n1len+n2len);
  2626.         return 0;
  2627.     }
  2628.  
  2629.     /* Okay, now move the exponent pointers to the most-significant word */
  2630.     e1 = BIGLITTLE(e1-e1len, e1+e1len-1);
  2631.     e2 = BIGLITTLE(e2-e2len, e2+e2len-1);
  2632.  
  2633.     /* Look up appropriate k-1 for the exponent - tblmask = 1<<(k-1) */
  2634.     w1bits = 0;
  2635.     while (e1bits > bnExpModThreshTable[w1bits])
  2636.         w1bits++;
  2637.     w2bits = 0;
  2638.     while (e2bits > bnExpModThreshTable[w2bits])
  2639.         w2bits++;
  2640.  
  2641.     assert(w1bits >= w2bits);
  2642.  
  2643.     /* Allocate working storage: two product buffers and the tables. */
  2644.     LBNALLOC(a, 2*mlen);
  2645.     if (!a)
  2646.         return -1;
  2647.     LBNALLOC(b, 2*mlen);
  2648.     if (!b) {
  2649.         LBNFREE(a, 2*mlen);
  2650.         return -1;
  2651.     }
  2652.  
  2653.     /* Convert to the appropriate table size: tblmask = 1<<(k-1) */
  2654.     tblmask = 1u << w1bits;
  2655.     /* Use buf2 for its size, temporarily */
  2656.     buf2 = 1u << w2bits;
  2657.  
  2658.     LBNALLOC(t, mlen);
  2659.     if (!t) {
  2660.         LBNFREE(b, 2*mlen);
  2661.         LBNFREE(a, 2*mlen);
  2662.         return -1;
  2663.     }
  2664.     table1[0] = t;
  2665.     table2[0] = result;
  2666.  
  2667.     /*
  2668.      * Okay, we now have some minimal-sized tables - expand them.
  2669.      * This is allowed to fail!  If so, scale back the table sizes
  2670.      * and proceed.  We allocate both tables at the same time
  2671.      * so if it fails partway through, they'll both be a reasonable
  2672.      * size rather than one huge and one tiny.
  2673.      * When i passes buf2 (the number of entries in the e2 window,
  2674.      * which may be less than the number of entries in the e1 window),
  2675.      * stop allocating e2 space.
  2676.      */
  2677.     for (i = 1; i < tblmask; i++) {
  2678.         LBNALLOC(t, mlen);
  2679.         if (!t)    /* Out of memory!  Quit the loop. */
  2680.             break;
  2681.         table1[i] = t;
  2682.         if (i < buf2) {
  2683.             LBNALLOC(t, mlen);
  2684.             if (!t) {
  2685.                 LBNFREE(table1[i], mlen);
  2686.                 break;
  2687.             }
  2688.             table2[i] = t;
  2689.         }
  2690.     }
  2691.  
  2692.     /* If we stopped, with i < tblmask, shrink the tables appropriately */
  2693.     while (tblmask > i) {
  2694.         w1bits--;
  2695.         tblmask >>= 1;
  2696.     }
  2697.     /* Free up our overallocations */
  2698.     while (--i > tblmask) {
  2699.         if (i < buf2)
  2700.             LBNFREE(table2[i], mlen);
  2701.         LBNFREE(table1[i], mlen);
  2702.     }
  2703.     /* And shrink the second window too, if needed */
  2704.     if (w2bits > w1bits) {
  2705.         w2bits = w1bits;
  2706.         buf2 = tblmask;
  2707.     }
  2708.  
  2709.     /*
  2710.      * From now on, use the w2bits variable for the difference
  2711.      * between w1bits and w2bits.
  2712.      */
  2713.     w2bits = w1bits-w2bits;
  2714.  
  2715.     /* Okay, fill in the tables */
  2716.  
  2717.     /* Compute the necessary modular inverse */
  2718.     inv = lbnMontInv1_64(mod[BIGLITTLE(-1,0)]);    /* LSW of modulus */
  2719.  
  2720.     /* Convert n1 to Montgomery form */
  2721.  
  2722.     /* Move n1 up "mlen" words into a */
  2723.     t = BIGLITTLE(a-mlen, a+mlen);
  2724.     lbnCopy_64(t, n1, n1len);
  2725.     lbnZero_64(a, mlen);
  2726.     /* Do the division - lose the quotient into the high-order words */
  2727.     (void)lbnDiv_64(t, a, mlen+n1len, mod, mlen);
  2728.     /* Copy into first table entry */
  2729.     lbnCopy_64(table1[0], a, mlen);
  2730.  
  2731.     /* Square a into b */
  2732.     lbnMontSquare_64(b, a, mod, mlen, inv);
  2733.  
  2734.     /* Use high half of b to initialize the first table */
  2735.     t = BIGLITTLE(b-mlen, b+mlen);
  2736.     for (i = 1; i < tblmask; i++) {
  2737.         lbnMontMul_64(a, t, table1[i-1], mod, mlen, inv);
  2738.         lbnCopy_64(table1[i], BIGLITTLE(a-mlen, a+mlen), mlen);
  2739.     }
  2740.  
  2741.     /* Convert n2 to Montgomery form */
  2742.  
  2743.     t = BIGLITTLE(a-mlen, a+mlen);
  2744.     /* Move n2 up "mlen" words into a */
  2745.     lbnCopy_64(t, n2, n2len);
  2746.     lbnZero_64(a, mlen);
  2747.     /* Do the division - lose the quotient into the high-order words */
  2748.     (void)lbnDiv_64(t, a, mlen+n2len, mod, mlen);
  2749.     /* Copy into first table entry */
  2750.     lbnCopy_64(table2[0], a, mlen);
  2751.  
  2752.     /* Square it into a */
  2753.     lbnMontSquare_64(a, table2[0], mod, mlen, inv);
  2754.     /* Copy to b, low half */
  2755.     lbnCopy_64(b, t, mlen);
  2756.  
  2757.     /* Use b to initialize the second table */
  2758.     for (i = 1; i < buf2; i++) {
  2759.         lbnMontMul_64(a, b, table2[i-1], mod, mlen, inv);
  2760.         lbnCopy_64(table2[i], t, mlen);
  2761.     }
  2762.  
  2763.     /*
  2764.      * Okay, a recap: at this point, the low part of b holds
  2765.      * n2^2, the high part holds n1^2, and the tables are
  2766.      * initialized with the odd powers of n1 and n2 from 1
  2767.      * through 2*tblmask-1 and 2*buf2-1.
  2768.      *
  2769.      * We might use those squares in b later, or we might not.
  2770.      */
  2771.  
  2772.     /* Initialze the fetch pointer */
  2773.     bitpos = (BNWORD64)1 << ((e1bits-1) & (64-1));    /* Initialize mask */
  2774.  
  2775.     /* This should point to the msbit of e1 */
  2776.     assert((*e1 & bitpos) != 0);
  2777.  
  2778.     /*
  2779.      * Pre-load the windows.  Becuase the window size is
  2780.      * never larger than the exponent size, there is no need to
  2781.      * detect running off the end of e1 in here.
  2782.      *
  2783.      * The read-ahead is controlled by e1len and the bitpos mask.
  2784.      * Note that this is *ahead* of e1bits, which tracks the
  2785.      * most significant end of the window.  The purpose of this
  2786.      * initialization is to get the two w1bits+1 bits apart,
  2787.      * like they should be.
  2788.      *
  2789.      * Note that bitpos and e1len together keep track of the
  2790.      * lookahead read pointer in the exponent that is used here.
  2791.      * e2len is not decremented, it is only ever compared with
  2792.      * e1len as *that* is decremented.
  2793.      */
  2794.     buf1 = buf2 = 0;
  2795.     for (i = 0; i <= w1bits; i++) {
  2796.         buf1 = (buf1 << 1) | ((*e1 & bitpos) != 0);
  2797.         if (e1len <= e2len)
  2798.             buf2 = (buf2 << 1) | ((*e2 & bitpos) != 0);
  2799.         bitpos >>= 1;
  2800.         if (!bitpos) {
  2801.             BIGLITTLE(e1++,e1--);
  2802.             if (e1len <= e2len)
  2803.                 BIGLITTLE(e2++,e2--);
  2804.             bitpos = (BNWORD64)1 << (64-1);
  2805.             e1len--;
  2806.         }
  2807.     }
  2808.     assert(buf1 & tblmask);
  2809.  
  2810.     /*
  2811.      * Set the pending multiply positions to a location that will
  2812.      * never be encountered, thus ensuring that nothing will happen
  2813.      * until the need for a multiply appears and one is scheduled.
  2814.      */
  2815.     mult1pos = mult2pos = e1bits;    /* A NULL value */
  2816.     mult1 = mult2 = 0;    /* Force a crash if we use these */
  2817.  
  2818.     /*
  2819.      * Okay, now begins the real work.  The first step is
  2820.      * slightly magic, so it's done outside the main loop,
  2821.      * but it's very similar to what's inside.
  2822.      */
  2823.     isone = 1;    /* Buffer is implicitly 1, so replace * by copy */
  2824.     e1bits--;    /* Start processing the first bit... */
  2825.  
  2826.     /*
  2827.      * This is just like the multiply in the loop, except that
  2828.      * - We know the msbit of buf1 is set, and
  2829.      * - We have the extra value n1^2 floating around.
  2830.      * So, do the usual computation, and if the result is that
  2831.      * the buffer should be multiplied by n1^1 immediately
  2832.      * (which we'd normally then square), we multiply it
  2833.      * (which reduces to a copy, which reduces to setting a flag)
  2834.      * by n1^2 and skip the squaring.  Thus, we do the
  2835.      * multiply and the squaring in one step.
  2836.      */
  2837.     assert(buf1 & tblmask);
  2838.     mult1pos = e1bits - w1bits;
  2839.     while ((buf1 & 1) == 0) {
  2840.         buf1 >>= 1;
  2841.         mult1pos++;
  2842.     }
  2843.     /* Intermediates can wrap, but final must NOT */
  2844.     assert(mult1pos <= e1bits);
  2845.     mult1 = table1[buf1>>1];
  2846.     buf1 = 0;
  2847.  
  2848.     /* Special case: use already-computed value sitting in buffer */
  2849.     if (mult1pos == e1bits)
  2850.         isone = 0;
  2851.  
  2852.     /*
  2853.      * The first multiply by a power of n2.  Similar, but
  2854.      * we might not even want to schedule a multiply if e2 is
  2855.      * shorter than e1, and the window might be shorter so
  2856.      * we have to leave the low w2bits bits alone.
  2857.      */
  2858.     if (buf2 & tblmask) {
  2859.         /* Remember low-order bits for later */
  2860.         i = buf2 & ((1u << w2bits) - 1);
  2861.         buf2 >>= w2bits;
  2862.         mult2pos = e1bits - w1bits + w2bits;
  2863.         while ((buf2 & 1) == 0) {
  2864.             buf2 >>= 1;
  2865.             mult2pos++;
  2866.         }
  2867.         assert(mult2pos <= e1bits);
  2868.         mult2 = table2[buf2>>1];
  2869.         buf2 = i;
  2870.  
  2871.         if (mult2pos == e1bits) {
  2872.             t = BIGLITTLE(b-mlen, b+mlen);
  2873.             if (isone) {
  2874.                 lbnCopy_64(t, b, mlen);    /* Copy low to high */
  2875.                 isone = 0;
  2876.             } else {
  2877.                 lbnMontMul_64(a, t, b, mod, mlen, inv);
  2878.                 t = a; a = b; b = t;
  2879.             }
  2880.         }
  2881.     }
  2882.  
  2883.     /*
  2884.      * At this point, the buffer (which is the high half of b)
  2885.      * holds either 1 (implicitly, as the "isone" flag is set),
  2886.      * n1^2, n2^2 or n1^2 * n2^2.
  2887.      */
  2888.  
  2889.     /*
  2890.      * The main loop.  The procedure is:
  2891.      * - Advance the windows
  2892.      * - If the most-significant bit of a window is set,
  2893.      *   schedule a multiply for the appropriate time in the
  2894.      *   future (may be immediately)
  2895.      * - Perform any pending multiples
  2896.      * - Check for termination
  2897.      * - Square the buffers
  2898.      *
  2899.      * At any given time, the acumulated product is held in
  2900.      * the high half of b.
  2901.      */
  2902.     for (;;) {
  2903.         e1bits--;
  2904.  
  2905.         /* Advance the windows */
  2906.         assert(buf1 < tblmask);
  2907.         buf1 <<= 1;
  2908.         assert(buf2 < tblmask);
  2909.         buf2 <<= 1;
  2910.         /*
  2911.          * This reads ahead of the current exponent position
  2912.          * (controlled by e1bits), so we have to be able to read
  2913.          * past the lsb of the exponents without error.
  2914.          */
  2915.         if (e1len) {
  2916.             buf1 |= ((*e1 & bitpos) != 0);
  2917.             if (e1len <= e2len)
  2918.                 buf2 |= ((*e2 & bitpos) != 0);
  2919.             bitpos >>= 1;
  2920.             if (!bitpos) {
  2921.                 BIGLITTLE(e1++,e1--);
  2922.                 if (e1len <= e2len)
  2923.                     BIGLITTLE(e2++,e2--);
  2924.                 bitpos = (BNWORD64)1 << (64-1);
  2925.                 e1len--;
  2926.             }
  2927.         }
  2928.  
  2929.         /* Examine the first window for pending multiplies */
  2930.         if (buf1 & tblmask) {
  2931.             mult1pos = e1bits - w1bits;
  2932.             while ((buf1 & 1) == 0) {
  2933.                 buf1 >>= 1;
  2934.                 mult1pos++;
  2935.             }
  2936.             /* Intermediates can wrap, but final must NOT */
  2937.             assert(mult1pos <= e1bits);
  2938.             mult1 = table1[buf1>>1];
  2939.             buf1 = 0;
  2940.         }
  2941.  
  2942.         /*
  2943.          * Examine the second window for pending multiplies.
  2944.          * Window 2 can be smaller than window 1, but we
  2945.          * keep the same number of bits in buf2, so we need
  2946.          * to ignore any low-order bits in the buffer when
  2947.          * computing what to multiply by, and recompute them
  2948.          * later.
  2949.          */
  2950.         if (buf2 & tblmask) {
  2951.             /* Remember low-order bits for later */
  2952.             i = buf2 & ((1u << w2bits) - 1);
  2953.             buf2 >>= w2bits;
  2954.             mult2pos = e1bits - w1bits + w2bits;
  2955.             while ((buf2 & 1) == 0) {
  2956.                 buf2 >>= 1;
  2957.                 mult2pos++;
  2958.             }
  2959.             assert(mult2pos <= e1bits);
  2960.             mult2 = table2[buf2>>1];
  2961.             buf2 = i;
  2962.         }
  2963.  
  2964.  
  2965.         /* If we have a pending multiply for e1, do it */
  2966.         if (e1bits == mult1pos) {
  2967.             /* Multiply by the table entry remembered previously */
  2968.             t = BIGLITTLE(b-mlen, b+mlen);
  2969.             if (isone) {
  2970.                 /* Multiply by 1 is a trivial case */
  2971.                 lbnCopy_64(t, mult1, mlen);
  2972.                 isone = 0;
  2973.             } else {
  2974.                 lbnMontMul_64(a, t, mult1, mod, mlen, inv);
  2975.                 /* Swap a and b */
  2976.                 t = a; a = b; b = t;
  2977.             }
  2978.         }
  2979.  
  2980.         /* If we have a pending multiply for e2, do it */
  2981.         if (e1bits == mult2pos) {
  2982.             /* Multiply by the table entry remembered previously */
  2983.             t = BIGLITTLE(b-mlen, b+mlen);
  2984.             if (isone) {
  2985.                 /* Multiply by 1 is a trivial case */
  2986.                 lbnCopy_64(t, mult2, mlen);
  2987.                 isone = 0;
  2988.             } else {
  2989.                 lbnMontMul_64(a, t, mult2, mod, mlen, inv);
  2990.                 /* Swap a and b */
  2991.                 t = a; a = b; b = t;
  2992.             }
  2993.         }
  2994.  
  2995.         /* Are we done? */
  2996.         if (!e1bits)
  2997.             break;
  2998.  
  2999.         /* Square the buffer */
  3000.         if (!isone) {
  3001.             t = BIGLITTLE(b-mlen, b+mlen);
  3002.             lbnMontSquare_64(a, t, mod, mlen, inv);
  3003.             /* Swap a and b */
  3004.             t = a; a = b; b = t;
  3005.         }
  3006.     } /* for (;;) */
  3007.  
  3008.     assert(!isone);
  3009.     assert(!buf1);
  3010.     assert(!buf2);
  3011.  
  3012.     /* DONE! */
  3013.  
  3014.     /* Convert result out of Montgomery form */
  3015.     t = BIGLITTLE(b-mlen, b+mlen);
  3016.     lbnCopy_64(b, t, mlen);
  3017.     lbnZero_64(t, mlen);
  3018.     lbnMontReduce_64(b, mod, mlen, inv);
  3019.     lbnCopy_64(result, t, mlen);
  3020.  
  3021.     /* Clean up - free intermediate storage */
  3022.     buf2 = tblmask >> w2bits;
  3023.     while (--tblmask) {
  3024.         if (tblmask < buf2)
  3025.             LBNFREE(table2[tblmask], mlen);
  3026.         LBNFREE(table1[tblmask], mlen);
  3027.     }
  3028.     t = table1[0];
  3029.     LBNFREE(t, mlen);
  3030.     LBNFREE(b, 2*mlen);
  3031.     LBNFREE(a, 2*mlen);
  3032.  
  3033.     return 0;    /* Success */
  3034. }
  3035.  
  3036. /*
  3037.  * 2^exp (mod mod).  This is an optimized version for use in Fermat
  3038.  * tests.  The input value of n is ignored; it is returned with
  3039.  * "mlen" words valid.
  3040.  */
  3041. int
  3042. lbnTwoExpMod_64(BNWORD64 *n, BNWORD64 const *exp, unsigned elen,
  3043.     BNWORD64 *mod, unsigned mlen)
  3044. {
  3045.     unsigned e;    /* Copy of high words of the exponent */
  3046.     unsigned bits;    /* Assorted counter of bits */
  3047.     BNWORD64 const *bitptr;
  3048.     BNWORD64 bitword, bitpos;
  3049.     BNWORD64 *a, *b, *a1;
  3050.     BNWORD64 inv;
  3051.  
  3052.     assert(mlen);
  3053.  
  3054.     bitptr = BIGLITTLE(exp-elen, exp+elen-1);
  3055.     bitword = *bitptr;
  3056.     assert(bitword);
  3057.  
  3058.     /* Clear n for future use. */
  3059.     lbnZero_64(n, mlen);
  3060.  
  3061.     bits = lbnBits_64(exp, elen);
  3062.     
  3063.     /* First, a couple of trivial cases. */
  3064.     if (bits <= 1) {
  3065.         /* 2 ^ 0 == 1,  2 ^ 1 == 2 */
  3066.         BIGLITTLE(n[-1],n[0]) = (BNWORD64)1<<elen;
  3067.         return 0;
  3068.     }
  3069.  
  3070.     /* Set bitpos to the most significant bit */
  3071.     bitpos = (BNWORD64)1 << ((bits-1) & (64-1));
  3072.  
  3073.     /* Now, count the bits in the modulus. */
  3074.     bits = lbnBits_64(mod, mlen);
  3075.     assert(bits > 1);    /* a 1-bit modulus is just stupid... */
  3076.  
  3077.     /*
  3078.      * We start with 1<<e, where "e" is as many high bits of the
  3079.      * exponent as we can manage without going over the modulus.
  3080.      * This first loop finds "e".
  3081.      */
  3082.     e = 1;
  3083.     while (elen) {
  3084.         /* Consume the first bit */
  3085.         bitpos >>= 1;
  3086.         if (!bitpos) {
  3087.             if (!--elen)
  3088.                 break;
  3089.             bitword = BIGLITTLE(*++bitptr,*--bitptr);
  3090.             bitpos = (BNWORD64)1<<(64-1);
  3091.         }
  3092.         e = (e << 1) | ((bitpos & bitword) != 0);
  3093.         if (e >= bits) {    /* Overflow!  Back out. */
  3094.             e >>= 1;
  3095.             break;
  3096.         }
  3097.     }
  3098.     /*
  3099.      * The bit in "bitpos" being examined by the bit buffer has NOT
  3100.      * been consumed yet.  This may be past the end of the exponent,
  3101.      * in which case elen == 1.
  3102.      */
  3103.  
  3104.     /* Okay, now, set bit "e" in n.  n is already zero. */
  3105.     inv = (BNWORD64)1 << (e & (64-1));
  3106.     e /= 64;
  3107.     BIGLITTLE(n[-e-1],n[e]) = inv;
  3108.     /*
  3109.      * The effective length of n in words is now "e+1".
  3110.      * This is used a little bit later.
  3111.      */
  3112.  
  3113.     if (!elen)
  3114.         return 0;    /* That was easy! */
  3115.  
  3116.     /*
  3117.      * We have now processed the first few bits.  The next step
  3118.      * is to convert this to Montgomery form for further squaring.
  3119.      */
  3120.  
  3121.     /* Allocate working storage: two product buffers */
  3122.     LBNALLOC(a, 2*mlen);
  3123.     if (!a)
  3124.         return -1;
  3125.     LBNALLOC(b, 2*mlen);
  3126.     if (!b) {
  3127.         LBNFREE(a, 2*mlen);
  3128.         return -1;
  3129.     }
  3130.  
  3131.     /* Convert n to Montgomery form */
  3132.     inv = BIGLITTLE(mod[-1],mod[0]);    /* LSW of modulus */
  3133.     assert(inv & 1);    /* Modulus must be odd */
  3134.     inv = lbnMontInv1_64(inv);
  3135.     /* Move n (length e+1, remember?) up "mlen" words into b */
  3136.     /* Note that we lie about a1 for a bit - it's pointing to b */
  3137.     a1 = BIGLITTLE(b-mlen,b+mlen);
  3138.     lbnCopy_64(a1, n, e+1);
  3139.     lbnZero_64(b, mlen);
  3140.     /* Do the division - dump the quotient into the high-order words */
  3141.     (void)lbnDiv_64(a1, b, mlen+e+1, mod, mlen);
  3142.     /*
  3143.      * Now do the first squaring and modular reduction to put
  3144.      * the number up in a1 where it belongs.
  3145.      */
  3146.     lbnMontSquare_64(a, b, mod, mlen, inv);
  3147.     /* Fix up a1 to point to where it should go. */
  3148.     a1 = BIGLITTLE(a-mlen,a+mlen);
  3149.  
  3150.     /*
  3151.      * Okay, now, a1 holds the number being accumulated, and
  3152.      * b is a scratch register.  Start working:
  3153.      */
  3154.     for (;;) {
  3155.         /*
  3156.          * Is the bit set?  If so, double a1 as well.
  3157.          * A modular doubling like this is very cheap.
  3158.          */
  3159.         if (bitpos & bitword) {
  3160.             /*
  3161.              * Double the number.  If there was a carry out OR
  3162.              * the result is greater than the modulus, subract
  3163.              * the modulus.
  3164.              */
  3165.             if (lbnDouble_64(a1, mlen) ||
  3166.                 lbnCmp_64(a1, mod, mlen) > 0)
  3167.                 (void)lbnSubN_64(a1, mod, mlen);
  3168.         }
  3169.  
  3170.         /* Advance to the next exponent bit */
  3171.         bitpos >>= 1;
  3172.         if (!bitpos) {
  3173.             if (!--elen)
  3174.                 break;    /* Done! */
  3175.             bitword = BIGLITTLE(*++bitptr,*--bitptr);
  3176.             bitpos = (BNWORD64)1<<(64-1);
  3177.         }
  3178.  
  3179.         /*
  3180.          * The elen/bitword/bitpos bit buffer is known to be
  3181.          * non-empty, i.e. there is at least one more unconsumed bit.
  3182.          * Thus, it's safe to square the number.
  3183.          */
  3184.         lbnMontSquare_64(b, a1, mod, mlen, inv);
  3185.         /* Rename result (in b) back to a (a1, really). */
  3186.         a1 = b; b = a; a = a1;
  3187.         a1 = BIGLITTLE(a-mlen,a+mlen);
  3188.     }
  3189.  
  3190.     /* DONE!  Just a little bit of cleanup... */
  3191.  
  3192.     /*
  3193.      * Convert result out of Montgomery form... this is
  3194.      * just a Montgomery reduction.
  3195.      */
  3196.     lbnCopy_64(a, a1, mlen);
  3197.     lbnZero_64(a1, mlen);
  3198.     lbnMontReduce_64(a, mod, mlen, inv);
  3199.     lbnCopy_64(n, a1, mlen);
  3200.  
  3201.     /* Clean up - free intermediate storage */
  3202.     LBNFREE(b, 2*mlen);
  3203.     LBNFREE(a, 2*mlen);
  3204.  
  3205.     return 0;    /* Success */
  3206. }
  3207.  
  3208.  
  3209. /*
  3210.  * Returns a substring of the big-endian array of bytes representation
  3211.  * of the bignum array based on two parameters, the least significant
  3212.  * byte number (0 to start with the least significant byte) and the
  3213.  * length.  I.e. the number returned is a representation of
  3214.  * (bn / 2^(8*lsbyte)) % 2 ^ (8*buflen).
  3215.  *
  3216.  * It is an error if the bignum is not at least buflen + lsbyte bytes
  3217.  * long.
  3218.  *
  3219.  * This code assumes that the compiler has the minimal intelligence
  3220.  * neded to optimize divides and modulo operations on an unsigned data
  3221.  * type with a power of two.
  3222.  */
  3223. void
  3224. lbnExtractBigBytes_64(BNWORD64 const *n, unsigned char *buf,
  3225.     unsigned lsbyte, unsigned buflen)
  3226. {
  3227.     BNWORD64 t = 0;    /* Needed to shut up uninitialized var warnings */
  3228.     unsigned shift;
  3229.  
  3230.     lsbyte += buflen;
  3231.  
  3232.     shift = (8 * lsbyte) % 64;
  3233.     lsbyte /= (64/8);    /* Convert to word offset */
  3234.     BIGLITTLE(n -= lsbyte, n += lsbyte);
  3235.  
  3236.     if (shift)
  3237.         t = BIGLITTLE(n[-1],n[0]);
  3238.  
  3239.     while (buflen--) {
  3240.         if (!shift) {
  3241.             t = BIGLITTLE(*n++,*--n);
  3242.             shift = 64;
  3243.         }
  3244.         shift -= 8;
  3245.         *buf++ = (unsigned char)(t>>shift);
  3246.     }
  3247. }
  3248.  
  3249. /*
  3250.  * Merge a big-endian array of bytes into a bignum array.
  3251.  * The array had better be big enough.  This is
  3252.  * equivalent to extracting the entire bignum into a
  3253.  * large byte array, copying the input buffer into the
  3254.  * middle of it, and converting back to a bignum.
  3255.  *
  3256.  * The buf is "len" bytes long, and its *last* byte is at
  3257.  * position "lsbyte" from the end of the bignum.
  3258.  *
  3259.  * Note that this is a pain to get right.  Fortunately, it's hardly
  3260.  * critical for efficiency.
  3261.  */
  3262. void
  3263. lbnInsertBigBytes_64(BNWORD64 *n, unsigned char const *buf,
  3264.                   unsigned lsbyte,  unsigned buflen)
  3265. {
  3266.     BNWORD64 t = 0;    /* Shut up uninitialized varibale warnings */
  3267.  
  3268.     lsbyte += buflen;
  3269.  
  3270.     BIGLITTLE(n -= lsbyte/(64/8), n += lsbyte/(64/8));
  3271.  
  3272.     /* Load up leading odd bytes */
  3273.     if (lsbyte % (64/8)) {
  3274.         t = BIGLITTLE(*--n,*n++);
  3275.         t >>= (lsbyte * 8) % 64;
  3276.     }
  3277.  
  3278.     /* The main loop - merge into t, storing at each word boundary. */
  3279.     while (buflen--) {
  3280.         t = (t << 8) | *buf++;
  3281.         if ((--lsbyte % (64/8)) == 0)
  3282.             BIGLITTLE(*n++,*--n) = t;
  3283.     }
  3284.  
  3285.     /* Merge odd bytes in t into last word */
  3286.     lsbyte = (lsbyte * 8) % 64;
  3287.     if (lsbyte) {
  3288.         t <<= lsbyte;
  3289.         t |= (((BNWORD64)1 << lsbyte) - 1) & BIGLITTLE(n[0],n[-1]);
  3290.         BIGLITTLE(n[0],n[-1]) = t;
  3291.     }
  3292.  
  3293.     return;
  3294. }
  3295.  
  3296. /*
  3297.  * Returns a substring of the little-endian array of bytes representation
  3298.  * of the bignum array based on two parameters, the least significant
  3299.  * byte number (0 to start with the least significant byte) and the
  3300.  * length.  I.e. the number returned is a representation of
  3301.  * (bn / 2^(8*lsbyte)) % 2 ^ (8*buflen).
  3302.  *
  3303.  * It is an error if the bignum is not at least buflen + lsbyte bytes
  3304.  * long.
  3305.  *
  3306.  * This code assumes that the compiler has the minimal intelligence
  3307.  * neded to optimize divides and modulo operations on an unsigned data
  3308.  * type with a power of two.
  3309.  */
  3310. void
  3311. lbnExtractLittleBytes_64(BNWORD64 const *n, unsigned char *buf,
  3312.     unsigned lsbyte, unsigned buflen)
  3313. {
  3314.     BNWORD64 t = 0;    /* Needed to shut up uninitialized var warnings */
  3315.  
  3316.     BIGLITTLE(n -= lsbyte/(64/8), n += lsbyte/(64/8));
  3317.  
  3318.     if (lsbyte % (64/8)) {
  3319.         t = BIGLITTLE(*--n,*n++);
  3320.         t >>= (lsbyte % (64/8)) * 8 ;
  3321.     }
  3322.  
  3323.     while (buflen--) {
  3324.         if ((lsbyte++ % (64/8)) == 0)
  3325.             t = BIGLITTLE(*--n,*n++);
  3326.         *buf++ = (unsigned char)t;
  3327.         t >>= 8;
  3328.     }
  3329. }
  3330.  
  3331. /*
  3332.  * Merge a little-endian array of bytes into a bignum array.
  3333.  * The array had better be big enough.  This is
  3334.  * equivalent to extracting the entire bignum into a
  3335.  * large byte array, copying the input buffer into the
  3336.  * middle of it, and converting back to a bignum.
  3337.  *
  3338.  * The buf is "len" bytes long, and its first byte is at
  3339.  * position "lsbyte" from the end of the bignum.
  3340.  *
  3341.  * Note that this is a pain to get right.  Fortunately, it's hardly
  3342.  * critical for efficiency.
  3343.  */
  3344. void
  3345. lbnInsertLittleBytes_64(BNWORD64 *n, unsigned char const *buf,
  3346.                   unsigned lsbyte,  unsigned buflen)
  3347. {
  3348.     BNWORD64 t = 0;    /* Shut up uninitialized varibale warnings */
  3349.  
  3350.     /* Move to most-significant end */
  3351.     lsbyte += buflen;
  3352.     buf += buflen;
  3353.  
  3354.     BIGLITTLE(n -= lsbyte/(64/8), n += lsbyte/(64/8));
  3355.  
  3356.     /* Load up leading odd bytes */
  3357.     if (lsbyte % (64/8)) {
  3358.         t = BIGLITTLE(*--n,*n++);
  3359.         t >>= (lsbyte * 8) % 64;
  3360.     }
  3361.  
  3362.     /* The main loop - merge into t, storing at each word boundary. */
  3363.     while (buflen--) {
  3364.         t = (t << 8) | *--buf;
  3365.         if ((--lsbyte % (64/8)) == 0)
  3366.             BIGLITTLE(*n++,*--n) = t;
  3367.     }
  3368.  
  3369.     /* Merge odd bytes in t into last word */
  3370.     lsbyte = (lsbyte * 8) % 64;
  3371.     if (lsbyte) {
  3372.         t <<= lsbyte;
  3373.         t |= (((BNWORD64)1 << lsbyte) - 1) & BIGLITTLE(n[0],n[-1]);
  3374.         BIGLITTLE(n[0],n[-1]) = t;
  3375.     }
  3376.  
  3377.     return;
  3378. }
  3379.  
  3380. #ifdef DEADCODE    /* This was a precursor to the more flexible lbnExtractBytes */
  3381. /*
  3382.  * Convert a big-endian array of bytes to a bignum.
  3383.  * Returns the number of words in the bignum.
  3384.  * Note the expression "64/8" for the number of bytes per word.
  3385.  * This is so the word-size adjustment will work.
  3386.  */
  3387. unsigned
  3388. lbnFromBytes_64(BNWORD64 *a, unsigned char const *b, unsigned blen)
  3389. {
  3390.     BNWORD64 t;
  3391.     unsigned alen = (blen + (64/8-1))/(64/8);
  3392.     BIGLITTLE(a -= alen, a += alen);
  3393.  
  3394.     while (blen) {
  3395.         t = 0;
  3396.         do {
  3397.             t = t << 8 | *b++;
  3398.         } while (--blen & (64/8-1));
  3399.         BIGLITTLE(*a++,*--a) = t;
  3400.     }
  3401.     return alen;
  3402. }
  3403. #endif
  3404.  
  3405. /*
  3406.  * Computes the GCD of a and b.  Modifies both arguments;
  3407.  * when it returns, one of them is the GCD and the other is trash.
  3408.  * The return value is the length of the GCD, with the sign telling
  3409.  * whether it is in a (+ve) or b (-ve).  Both inputs must have
  3410.  * one extra word of precision.  alen must be >= blen.
  3411.  *
  3412.  * TODO: use the binary algorithm (Knuth section 4.5.2, algorithm B).
  3413.  * This is based on taking out common powers of 2, then repeatedly:
  3414.  * gcd(2*u,v) = gcd(u,2*v) = gcd(u,v) - isolated powers of 2 can be deleted.
  3415.  * gcd(u,v) = gcd(u-v,v) - the numbers can be easily reduced.
  3416.  * It gets less reduction per step, but the steps are much faster than
  3417.  * the division case.
  3418.  */
  3419. int
  3420. lbnGcd_64(BNWORD64 *a, unsigned alen, BNWORD64 *b, unsigned blen)
  3421. {
  3422.     assert(alen >= blen);
  3423.  
  3424.     while (blen != 0) {
  3425.         (void)lbnDiv_64(BIGLITTLE(a-blen,a+blen), a, alen, b, blen);
  3426.         alen = lbnNorm_64(a, blen);
  3427.         if (alen == 0)
  3428.             return -(int)blen;
  3429.         (void)lbnDiv_64(BIGLITTLE(b-alen,b+alen), b, blen, a, alen);
  3430.         blen = lbnNorm_64(b, alen);
  3431.     }
  3432.     return alen;
  3433. }
  3434.  
  3435. /*
  3436.  * Invert "a" modulo "mod" using the extended Euclidean algorithm.
  3437.  * Note that this only computes one of the cosequences, and uses the
  3438.  * theorem that the signs flip every step and the absolute value of
  3439.  * the cosequence values are always bounded by the modulus to avoid
  3440.  * having to work with negative numbers.
  3441.  * gcd(a,mod) had better equal 1.  Returns 1 if the GCD is NOT 1.
  3442.  * a must be one word longer than "mod".  It is overwritten with the
  3443.  * result.
  3444.  * TODO: Use Richard Schroeppel's *much* faster algorithm.
  3445.  */
  3446. int
  3447. lbnInv_64(BNWORD64 *a, unsigned alen, BNWORD64 const *mod, unsigned mlen)
  3448. {
  3449.     BNWORD64 *b;    /* Hold a copy of mod during GCD reduction */
  3450.     BNWORD64 *p;    /* Temporary for products added to t0 and t1 */
  3451.     BNWORD64 *t0, *t1;    /* Inverse accumulators */
  3452.     BNWORD64 cy;
  3453.     unsigned blen, t0len, t1len, plen;
  3454.  
  3455.     alen = lbnNorm_64(a, alen);
  3456.     if (!alen)
  3457.         return 1;    /* No inverse */
  3458.  
  3459.     mlen = lbnNorm_64(mod, mlen);
  3460.  
  3461.     assert (alen <= mlen);
  3462.  
  3463.     /* Inverse of 1 is 1 */
  3464.     if (alen == 1 && BIGLITTLE(a[-1],a[0]) == 1) {
  3465.         lbnZero_64(BIGLITTLE(a-alen,a+alen), mlen-alen);
  3466.         return 0;
  3467.     }
  3468.  
  3469.     /* Allocate a pile of space */
  3470.     LBNALLOC(b, mlen+1);
  3471.     if (b) {
  3472.         /*
  3473.          * Although products are guaranteed to always be less than the
  3474.          * modulus, it can involve multiplying two 3-word numbers to
  3475.          * get a 5-word result, requiring a 6th word to store a 0
  3476.          * temporarily.  Thus, mlen + 1.
  3477.          */
  3478.         LBNALLOC(p, mlen+1);
  3479.         if (p) {
  3480.             LBNALLOC(t0, mlen);
  3481.             if (t0) {
  3482.                 LBNALLOC(t1, mlen);
  3483.                 if (t1)
  3484.                         goto allocated;
  3485.                 LBNFREE(t1, mlen);
  3486.             }
  3487.             LBNFREE(p, mlen+1);
  3488.         }
  3489.         LBNFREE(b, mlen+1);
  3490.     }
  3491.     return -1;
  3492.  
  3493. allocated:
  3494.  
  3495.     /* Set t0 to 1 */
  3496.     t0len = 1;
  3497.     BIGLITTLE(t0[-1],t0[0]) = 1;
  3498.     
  3499.     /* b = mod */
  3500.     lbnCopy_64(b, mod, mlen);
  3501.     /* blen = mlen (implicitly) */
  3502.     
  3503.     /* t1 = b / a; b = b % a */
  3504.     cy = lbnDiv_64(t1, b, mlen, a, alen);
  3505.     *(BIGLITTLE(t1-(mlen-alen)-1,t1+(mlen-alen))) = cy;
  3506.     t1len = lbnNorm_64(t1, mlen-alen+1);
  3507.     blen = lbnNorm_64(b, alen);
  3508.  
  3509.     /* while (b > 1) */
  3510.     while (blen > 1 || BIGLITTLE(b[-1],b[0]) != (BNWORD64)1) {
  3511.         /* q = a / b; a = a % b; */
  3512.         if (alen < blen || (alen == blen && lbnCmp_64(a, a, alen) < 0))
  3513.             assert(0);
  3514.         cy = lbnDiv_64(BIGLITTLE(a-blen,a+blen), a, alen, b, blen);
  3515.         *(BIGLITTLE(a-alen-1,a+alen)) = cy;
  3516.         plen = lbnNorm_64(BIGLITTLE(a-blen,a+blen), alen-blen+1);
  3517.         assert(plen);
  3518.         alen = lbnNorm_64(a, blen);
  3519.         if (!alen)
  3520.             goto failure;    /* GCD not 1 */
  3521.  
  3522.         /* t0 += q * t1; */
  3523.         assert(plen+t1len <= mlen+1);
  3524.         lbnMul_64(p, BIGLITTLE(a-blen,a+blen), plen, t1, t1len);
  3525.         plen = lbnNorm_64(p, plen + t1len);
  3526.         assert(plen <= mlen);
  3527.         if (plen > t0len) {
  3528.             lbnZero_64(BIGLITTLE(t0-t0len,t0+t0len), plen-t0len);
  3529.             t0len = plen;
  3530.         }
  3531.         cy = lbnAddN_64(t0, p, plen);
  3532.         if (cy) {
  3533.             if (t0len > plen) {
  3534.                 cy = lbnAdd1_64(BIGLITTLE(t0-plen,t0+plen),
  3535.                         t0len-plen, cy);
  3536.             }
  3537.             if (cy) {
  3538.                 BIGLITTLE(t0[-t0len-1],t0[t0len]) = cy;
  3539.                 t0len++;
  3540.             }
  3541.         }
  3542.  
  3543.         /* if (a <= 1) return a ? t0 : FAIL; */
  3544.         if (alen <= 1 && BIGLITTLE(a[-1],a[0]) == (BNWORD64)1) {
  3545.             if (alen == 0)
  3546.                 goto failure;    /* FAIL */
  3547.             assert(t0len <= mlen);
  3548.             lbnCopy_64(a, t0, t0len);
  3549.             lbnZero_64(BIGLITTLE(a-t0len, a+t0len), mlen-t0len);
  3550.             goto success;
  3551.         }
  3552.  
  3553.         /* q = b / a; b = b % a; */
  3554.         if (blen < alen || (blen == alen && lbnCmp_64(b, a, alen) < 0))
  3555.             assert(0);
  3556.         cy = lbnDiv_64(BIGLITTLE(b-alen,b+alen), b, blen, a, alen);
  3557.         *(BIGLITTLE(b-blen-1,b+blen)) = cy;
  3558.         plen = lbnNorm_64(BIGLITTLE(b-alen,b+alen), blen-alen+1);
  3559.         assert(plen);
  3560.         blen = lbnNorm_64(b, alen);
  3561.         if (!blen)
  3562.             goto failure;    /* GCD not 1 */
  3563.  
  3564.         /* t1 += q * t0; */
  3565.         assert(plen+t0len <= mlen+1);
  3566.         lbnMul_64(p, BIGLITTLE(b-alen,b+alen), plen, t0, t0len);
  3567.         plen = lbnNorm_64(p, plen + t0len);
  3568.         assert(plen <= mlen);
  3569.         if (plen > t1len) {
  3570.             lbnZero_64(BIGLITTLE(t1-t1len,t1+t1len), plen-t1len);
  3571.             t1len = plen;
  3572.         }
  3573.         cy = lbnAddN_64(t1, p, plen);
  3574.         if (cy) {
  3575.             if (t1len > plen) {
  3576.                 cy = lbnAdd1_64(BIGLITTLE(t1-plen,t0+plen),
  3577.                         t1len-plen, cy);
  3578.             }
  3579.             if (cy) {
  3580.                 BIGLITTLE(t1[-t1len-1],t1[t1len]) = cy;
  3581.                 t1len++;
  3582.             }
  3583.         }
  3584.     }
  3585.  
  3586.     if (!blen)
  3587.         goto failure;    /* gcd(a, mod) != 1 -- FAIL */
  3588.  
  3589.     /* return mod-t1 */
  3590.     lbnCopy_64(a, mod, mlen);
  3591.     assert(t1len <= mlen);
  3592.     cy = lbnSubN_64(a, t1, t1len);
  3593.     if (cy) {
  3594.         assert(mlen > t1len);
  3595.         cy = lbnSub1_64(BIGLITTLE(a-t1len, a+t1len), mlen-t1len, cy);
  3596.         assert(!cy);
  3597.     }
  3598.  
  3599. success:
  3600.     LBNFREE(t1, mlen);
  3601.     LBNFREE(t0, mlen);
  3602.     LBNFREE(p, mlen+1);
  3603.     LBNFREE(b, mlen+1);
  3604.     
  3605.     return 0;
  3606.  
  3607. failure:
  3608.     LBNFREE(t1, mlen);
  3609.     LBNFREE(t0, mlen);
  3610.     LBNFREE(p, mlen+1);
  3611.     LBNFREE(b, mlen+1);
  3612.     
  3613.     return 1;
  3614. }
  3615.