* #35924 (zrtp): switch to libzrtpcpp
diff --git a/jni/libzrtp/sources/bnlib/lbn32.c b/jni/libzrtp/sources/bnlib/lbn32.c
new file mode 100644
index 0000000..73fedcb
--- /dev/null
+++ b/jni/libzrtp/sources/bnlib/lbn32.c
@@ -0,0 +1,4073 @@
+/*
+ * lbn32.c - Low-level bignum routines, 32-bit version.
+ *
+ * Copyright (c) 1995 Colin Plumb. All rights reserved.
+ * For licensing and other legal details, see the file legal.c.
+ *
+ * NOTE: the magic constants "32" and "64" appear in many places in this
+ * file, including inside identifiers. Because it is not possible to
+ * ask "#ifdef" of a macro expansion, it is not possible to use the
+ * preprocessor to conditionalize these properly. Thus, this file is
+ * intended to be edited with textual search and replace to produce
+ * alternate word size versions. Any reference to the number of bits
+ * in a word must be the string "32", and that string must not appear
+ * otherwise. Any reference to twice this number must appear as "64",
+ * which likewise must not appear otherwise. Is that clear?
+ *
+ * Remember, when doubling the bit size replace the larger number (64)
+ * first, then the smaller (32). When halving the bit size, do the
+ * opposite. Otherwise, things will get wierd. Also, be sure to replace
+ * every instance that appears. (:%s/foo/bar/g in vi)
+ *
+ * These routines work with a pointer to the least-significant end of
+ * an array of WORD32s. The BIG(x), LITTLE(y) and BIGLTTLE(x,y) macros
+ * defined in lbn.h (which expand to x on a big-edian machine and y on a
+ * little-endian machine) are used to conditionalize the code to work
+ * either way. If you have no assembly primitives, it doesn't matter.
+ * Note that on a big-endian machine, the least-significant-end pointer
+ * is ONE PAST THE END. The bytes are ptr[-1] through ptr[-len].
+ * On little-endian, they are ptr[0] through ptr[len-1]. This makes
+ * perfect sense if you consider pointers to point *between* bytes rather
+ * than at them.
+ *
+ * Because the array index values are unsigned integers, ptr[-i]
+ * may not work properly, since the index -i is evaluated as an unsigned,
+ * and if pointers are wider, zero-extension will produce a positive
+ * number rahter than the needed negative. The expression used in this
+ * code, *(ptr-i) will, however, work. (The array syntax is equivalent
+ * to *(ptr+-i), which is a pretty subtle difference.)
+ *
+ * Many of these routines will get very unhappy if fed zero-length inputs.
+ * They use assert() to enforce this. An higher layer of code must make
+ * sure that these aren't called with zero-length inputs.
+ *
+ * Any of these routines can be replaced with more efficient versions
+ * elsewhere, by just #defining their names. If one of the names
+ * is #defined, the C code is not compiled in and no declaration is
+ * made. Use the BNINCLUDE file to do that. Typically, you compile
+ * asm subroutines with the same name and just, e.g.
+ * #define lbnMulAdd1_32 lbnMulAdd1_32
+ *
+ * If you want to write asm routines, start with lbnMulAdd1_32().
+ * This is the workhorse of modular exponentiation. lbnMulN1_32() is
+ * also used a fair bit, although not as much and it's defined in terms
+ * of lbnMulAdd1_32 if that has a custom version. lbnMulSub1_32 and
+ * lbnDiv21_32 are used in the usual division and remainder finding.
+ * (Not the Montgomery reduction used in modular exponentiation, though.)
+ * Once you have lbnMulAdd1_32 defined, writing the other two should
+ * be pretty easy. (Just make sure you get the sign of the subtraction
+ * in lbnMulSub1_32 right - it's dest = dest - source * k.)
+ *
+ * The only definitions that absolutely need a double-word (BNWORD64)
+ * type are lbnMulAdd1_32 and lbnMulSub1_32; if those are provided,
+ * the rest follows. lbnDiv21_32, however, is a lot slower unless you
+ * have them, and lbnModQ_32 takes after it. That one is used quite a
+ * bit for prime sieving.
+ */
+
+#ifndef HAVE_CONFIG_H
+#define HAVE_CONFIG_H 0
+#endif
+#if HAVE_CONFIG_H
+#include <bnconfig.h>
+#endif
+
+/*
+ * Some compilers complain about #if FOO if FOO isn't defined,
+ * so do the ANSI-mandated thing explicitly...
+ */
+#ifndef NO_ASSERT_H
+#define NO_ASSERT_H 0
+#endif
+#ifndef NO_STRING_H
+#define NO_STRING_H 0
+#endif
+#ifndef HAVE_STRINGS_H
+#define HAVE_STRINGS_H 0
+#endif
+#ifndef NEED_MEMORY_H
+#define NEED_MEMORY_H 0
+#endif
+
+#if !NO_ASSERT_H
+#include <assert.h>
+#else
+#define assert(x) (void)0
+#endif
+
+#if !NO_STRING_H
+#include <string.h> /* For memcpy */
+#elif HAVE_STRINGS_H
+#include <strings.h>
+#endif
+#if NEED_MEMORY_H
+#include <memory.h>
+#endif
+
+#include "lbn.h"
+#include "lbn32.h"
+#include "lbnmem.h"
+
+#include "kludge.h"
+
+#ifndef BNWORD32
+#error 32-bit bignum library requires a 32-bit data type
+#endif
+
+/* If this is defined, include bnYield() calls */
+#if BNYIELD
+extern int (*bnYield)(void); /* From bn.c */
+#endif
+
+/*
+ * Most of the multiply (and Montgomery reduce) routines use an outer
+ * loop that iterates over one of the operands - a so-called operand
+ * scanning approach. One big advantage of this is that the assembly
+ * support routines are simpler. The loops can be rearranged to have
+ * an outer loop that iterates over the product, a so-called product
+ * scanning approach. This has the advantage of writing less data
+ * and doing fewer adds to memory, so is supposedly faster. Some
+ * code has been written using a product-scanning approach, but
+ * it appears to be slower, so it is turned off by default. Some
+ * experimentation would be appreciated.
+ *
+ * (The code is also annoying to get right and not very well commented,
+ * one of my pet peeves about math libraries. I'm sorry.)
+ */
+#ifndef PRODUCT_SCAN
+#define PRODUCT_SCAN 0
+#endif
+
+/*
+ * Copy an array of words. <Marvin mode on> Thrilling, isn't it? </Marvin>
+ * This is a good example of how the byte offsets and BIGLITTLE() macros work.
+ * Another alternative would have been
+ * memcpy(dest BIG(-len), src BIG(-len), len*sizeof(BNWORD32)), but I find that
+ * putting operators into conditional macros is confusing.
+ */
+#ifndef lbnCopy_32
+void
+lbnCopy_32(BNWORD32 *dest, BNWORD32 const *src, unsigned len)
+{
+ memcpy(BIGLITTLE(dest-len,dest), BIGLITTLE(src-len,src),
+ len * sizeof(*src));
+}
+#endif /* !lbnCopy_32 */
+
+/*
+ * Fill n words with zero. This does it manually rather than calling
+ * memset because it can assume alignment to make things faster while
+ * memset can't. Note how big-endian numbers are naturally addressed
+ * using predecrement, while little-endian is postincrement.
+ */
+#ifndef lbnZero_32
+void
+lbnZero_32(BNWORD32 *num, unsigned len)
+{
+ while (len--)
+ BIGLITTLE(*--num,*num++) = 0;
+}
+#endif /* !lbnZero_32 */
+
+/*
+ * Negate an array of words.
+ * Negation is subtraction from zero. Negating low-order words
+ * entails doing nothing until a non-zero word is hit. Once that
+ * is negated, a borrow is generated and never dies until the end
+ * of the number is hit. Negation with borrow, -x-1, is the same as ~x.
+ * Repeat that until the end of the number.
+ *
+ * Doesn't return borrow out because that's pretty useless - it's
+ * always set unless the input is 0, which is easy to notice in
+ * normalized form.
+ */
+#ifndef lbnNeg_32
+void
+lbnNeg_32(BNWORD32 *num, unsigned len)
+{
+ assert(len);
+
+ /* Skip low-order zero words */
+ while (BIGLITTLE(*--num,*num) == 0) {
+ if (!--len)
+ return;
+ LITTLE(num++;)
+ }
+ /* Negate the lowest-order non-zero word */
+ *num = -*num;
+ /* Complement all the higher-order words */
+ while (--len) {
+ BIGLITTLE(--num,++num);
+ *num = ~*num;
+ }
+}
+#endif /* !lbnNeg_32 */
+
+
+/*
+ * lbnAdd1_32: add the single-word "carry" to the given number.
+ * Used for minor increments and propagating the carry after
+ * adding in a shorter bignum.
+ *
+ * Technique: If we have a double-width word, presumably the compiler
+ * can add using its carry in inline code, so we just use a larger
+ * accumulator to compute the carry from the first addition.
+ * If not, it's more complex. After adding the first carry, which may
+ * be > 1, compare the sum and the carry. If the sum wraps (causing a
+ * carry out from the addition), the result will be less than each of the
+ * inputs, since the wrap subtracts a number (2^32) which is larger than
+ * the other input can possibly be. If the sum is >= the carry input,
+ * return success immediately.
+ * In either case, if there is a carry, enter a loop incrementing words
+ * until one does not wrap. Since we are adding 1 each time, the wrap
+ * will be to 0 and we can test for equality.
+ */
+#ifndef lbnAdd1_32 /* If defined, it's provided as an asm subroutine */
+#ifdef BNWORD64
+BNWORD32
+lbnAdd1_32(BNWORD32 *num, unsigned len, BNWORD32 carry)
+{
+ BNWORD64 t;
+ assert(len > 0); /* Alternative: if (!len) return carry */
+
+ t = (BNWORD64)BIGLITTLE(*--num,*num) + carry;
+ BIGLITTLE(*num,*num++) = (BNWORD32)t;
+ if ((t >> 32) == 0)
+ return 0;
+ while (--len) {
+ if (++BIGLITTLE(*--num,*num++) != 0)
+ return 0;
+ }
+ return 1;
+}
+#else /* no BNWORD64 */
+BNWORD32
+lbnAdd1_32(BNWORD32 *num, unsigned len, BNWORD32 carry)
+{
+ assert(len > 0); /* Alternative: if (!len) return carry */
+
+ if ((BIGLITTLE(*--num,*num++) += carry) >= carry)
+ return 0;
+ while (--len) {
+ if (++BIGLITTLE(*--num,*num++) != 0)
+ return 0;
+ }
+ return 1;
+}
+#endif
+#endif/* !lbnAdd1_32 */
+
+/*
+ * lbnSub1_32: subtract the single-word "borrow" from the given number.
+ * Used for minor decrements and propagating the borrow after
+ * subtracting a shorter bignum.
+ *
+ * Technique: Similar to the add, above. If there is a double-length type,
+ * use that to generate the first borrow.
+ * If not, after subtracting the first borrow, which may be > 1, compare
+ * the difference and the *negative* of the carry. If the subtract wraps
+ * (causing a borrow out from the subtraction), the result will be at least
+ * as large as -borrow. If the result < -borrow, then no borrow out has
+ * appeared and we may return immediately, except when borrow == 0. To
+ * deal with that case, use the identity that -x = ~x+1, and instead of
+ * comparing < -borrow, compare for <= ~borrow.
+ * Either way, if there is a borrow out, enter a loop decrementing words
+ * until a non-zero word is reached.
+ *
+ * Note the cast of ~borrow to (BNWORD32). If the size of an int is larger
+ * than BNWORD32, C rules say the number is expanded for the arithmetic, so
+ * the inversion will be done on an int and the value won't be quite what
+ * is expected.
+ */
+#ifndef lbnSub1_32 /* If defined, it's provided as an asm subroutine */
+#ifdef BNWORD64
+BNWORD32
+lbnSub1_32(BNWORD32 *num, unsigned len, BNWORD32 borrow)
+{
+ BNWORD64 t;
+ assert(len > 0); /* Alternative: if (!len) return borrow */
+
+ t = (BNWORD64)BIGLITTLE(*--num,*num) - borrow;
+ BIGLITTLE(*num,*num++) = (BNWORD32)t;
+ if ((t >> 32) == 0)
+ return 0;
+ while (--len) {
+ if ((BIGLITTLE(*--num,*num++))-- != 0)
+ return 0;
+ }
+ return 1;
+}
+#else /* no BNWORD64 */
+BNWORD32
+lbnSub1_32(BNWORD32 *num, unsigned len, BNWORD32 borrow)
+{
+ assert(len > 0); /* Alternative: if (!len) return borrow */
+
+ if ((BIGLITTLE(*--num,*num++) -= borrow) <= (BNWORD32)~borrow)
+ return 0;
+ while (--len) {
+ if ((BIGLITTLE(*--num,*num++))-- != 0)
+ return 0;
+ }
+ return 1;
+}
+#endif
+#endif /* !lbnSub1_32 */
+
+/*
+ * lbnAddN_32: add two bignums of the same length, returning the carry (0 or 1).
+ * One of the building blocks, along with lbnAdd1, of adding two bignums of
+ * differing lengths.
+ *
+ * Technique: Maintain a word of carry. If there is no double-width type,
+ * use the same technique as in lbnAdd1, above, to maintain the carry by
+ * comparing the inputs. Adding the carry sources is used as an OR operator;
+ * at most one of the two comparisons can possibly be true. The first can
+ * only be true if carry == 1 and x, the result, is 0. In that case the
+ * second can't possibly be true.
+ */
+#ifndef lbnAddN_32
+#ifdef BNWORD64
+BNWORD32
+lbnAddN_32(BNWORD32 *num1, BNWORD32 const *num2, unsigned len)
+{
+ BNWORD64 t;
+
+ assert(len > 0);
+
+ t = (BNWORD64)BIGLITTLE(*--num1,*num1) + BIGLITTLE(*--num2,*num2++);
+ BIGLITTLE(*num1,*num1++) = (BNWORD32)t;
+ while (--len) {
+ t = (BNWORD64)BIGLITTLE(*--num1,*num1) +
+ (BNWORD64)BIGLITTLE(*--num2,*num2++) + (t >> 32);
+ BIGLITTLE(*num1,*num1++) = (BNWORD32)t;
+ }
+
+ return (BNWORD32)(t>>32);
+}
+#else /* no BNWORD64 */
+BNWORD32
+lbnAddN_32(BNWORD32 *num1, BNWORD32 const *num2, unsigned len)
+{
+ BNWORD32 x, carry = 0;
+
+ assert(len > 0); /* Alternative: change loop to test at start */
+
+ do {
+ x = BIGLITTLE(*--num2,*num2++);
+ carry = (x += carry) < carry;
+ carry += (BIGLITTLE(*--num1,*num1++) += x) < x;
+ } while (--len);
+
+ return carry;
+}
+#endif
+#endif /* !lbnAddN_32 */
+
+/*
+ * lbnSubN_32: add two bignums of the same length, returning the carry (0 or 1).
+ * One of the building blocks, along with subn1, of subtracting two bignums of
+ * differing lengths.
+ *
+ * Technique: If no double-width type is availble, maintain a word of borrow.
+ * First, add the borrow to the subtrahend (did you have to learn all those
+ * awful words in elementary school, too?), and if it overflows, set the
+ * borrow again. Then subtract the modified subtrahend from the next word
+ * of input, using the same technique as in subn1, above.
+ * Adding the borrows is used as an OR operator; at most one of the two
+ * comparisons can possibly be true. The first can only be true if
+ * borrow == 1 and x, the result, is 0. In that case the second can't
+ * possibly be true.
+ *
+ * In the double-word case, (BNWORD32)-(t>>32) is subtracted, rather than
+ * adding t>>32, because the shift would need to sign-extend and that's
+ * not guaranteed to happen in ANSI C, even with signed types.
+ */
+#ifndef lbnSubN_32
+#ifdef BNWORD64
+BNWORD32
+lbnSubN_32(BNWORD32 *num1, BNWORD32 const *num2, unsigned len)
+{
+ BNWORD64 t;
+
+ assert(len > 0);
+
+ t = (BNWORD64)BIGLITTLE(*--num1,*num1) - BIGLITTLE(*--num2,*num2++);
+ BIGLITTLE(*num1,*num1++) = (BNWORD32)t;
+
+ while (--len) {
+ t = (BNWORD64)BIGLITTLE(*--num1,*num1) -
+ (BNWORD64)BIGLITTLE(*--num2,*num2++) - (BNWORD32)-(t >> 32);
+ BIGLITTLE(*num1,*num1++) = (BNWORD32)t;
+ }
+
+ return -(BNWORD32)(t>>32);
+}
+#else
+BNWORD32
+lbnSubN_32(BNWORD32 *num1, BNWORD32 const *num2, unsigned len)
+{
+ BNWORD32 x, borrow = 0;
+
+ assert(len > 0); /* Alternative: change loop to test at start */
+
+ do {
+ x = BIGLITTLE(*--num2,*num2++);
+ borrow = (x += borrow) < borrow;
+ borrow += (BIGLITTLE(*--num1,*num1++) -= x) > (BNWORD32)~x;
+ } while (--len);
+
+ return borrow;
+}
+#endif
+#endif /* !lbnSubN_32 */
+
+#ifndef lbnCmp_32
+/*
+ * lbnCmp_32: compare two bignums of equal length, returning the sign of
+ * num1 - num2. (-1, 0 or +1).
+ *
+ * Technique: Change the little-endian pointers to big-endian pointers
+ * and compare from the most-significant end until a difference if found.
+ * When it is, figure out the sign of the difference and return it.
+ */
+int
+lbnCmp_32(BNWORD32 const *num1, BNWORD32 const *num2, unsigned len)
+{
+ BIGLITTLE(num1 -= len, num1 += len);
+ BIGLITTLE(num2 -= len, num2 += len);
+
+ while (len--) {
+ if (BIGLITTLE(*num1++ != *num2++, *--num1 != *--num2)) {
+ if (BIGLITTLE(num1[-1] < num2[-1], *num1 < *num2))
+ return -1;
+ else
+ return 1;
+ }
+ }
+ return 0;
+}
+#endif /* !lbnCmp_32 */
+
+/*
+ * mul32_ppmmaa(ph,pl,x,y,a,b) is an optional routine that
+ * computes (ph,pl) = x * y + a + b. mul32_ppmma and mul32_ppmm
+ * are simpler versions. If you want to be lazy, all of these
+ * can be defined in terms of the others, so here we create any
+ * that have not been defined in terms of the ones that have been.
+ */
+
+/* Define ones with fewer a's in terms of ones with more a's */
+#if !defined(mul32_ppmma) && defined(mul32_ppmmaa)
+#define mul32_ppmma(ph,pl,x,y,a) mul32_ppmmaa(ph,pl,x,y,a,0)
+#endif
+
+#if !defined(mul32_ppmm) && defined(mul32_ppmma)
+#define mul32_ppmm(ph,pl,x,y) mul32_ppmma(ph,pl,x,y,0)
+#endif
+
+/*
+ * Use this definition to test the mul32_ppmm-based operations on machines
+ * that do not provide mul32_ppmm. Change the final "0" to a "1" to
+ * enable it.
+ */
+#if !defined(mul32_ppmm) && defined(BNWORD64) && 0 /* Debugging */
+#define mul32_ppmm(ph,pl,x,y) \
+ ({BNWORD64 _ = (BNWORD64)(x)*(y); (pl) = _; (ph) = _>>32;})
+#endif
+
+#if defined(mul32_ppmm) && !defined(mul32_ppmma)
+#define mul32_ppmma(ph,pl,x,y,a) \
+ (mul32_ppmm(ph,pl,x,y), (ph) += ((pl) += (a)) < (a))
+#endif
+
+#if defined(mul32_ppmma) && !defined(mul32_ppmmaa)
+#define mul32_ppmmaa(ph,pl,x,y,a,b) \
+ (mul32_ppmma(ph,pl,x,y,a), (ph) += ((pl) += (b)) < (b))
+#endif
+
+/*
+ * lbnMulN1_32: Multiply an n-word input by a 1-word input and store the
+ * n+1-word product. This uses either the mul32_ppmm and mul32_ppmma
+ * macros, or C multiplication with the BNWORD64 type. This uses mul32_ppmma
+ * if available, assuming you won't bother defining it unless you can do
+ * better than the normal multiplication.
+ */
+#ifndef lbnMulN1_32
+#ifdef lbnMulAdd1_32 /* If we have this asm primitive, use it. */
+void
+lbnMulN1_32(BNWORD32 *out, BNWORD32 const *in, unsigned len, BNWORD32 k)
+{
+ lbnZero_32(out, len);
+ BIGLITTLE(*(out-len-1),*(out+len)) = lbnMulAdd1_32(out, in, len, k);
+}
+#elif defined(mul32_ppmm)
+void
+lbnMulN1_32(BNWORD32 *out, BNWORD32 const *in, unsigned len, BNWORD32 k)
+{
+ BNWORD32 carry, carryin;
+
+ assert(len > 0);
+
+ BIG(--out;--in;);
+ mul32_ppmm(carry, *out, *in, k);
+ LITTLE(out++;in++;)
+
+ while (--len) {
+ BIG(--out;--in;)
+ carryin = carry;
+ mul32_ppmma(carry, *out, *in, k, carryin);
+ LITTLE(out++;in++;)
+ }
+ BIGLITTLE(*--out,*out) = carry;
+}
+#elif defined(BNWORD64)
+void
+lbnMulN1_32(BNWORD32 *out, BNWORD32 const *in, unsigned len, BNWORD32 k)
+{
+ BNWORD64 p;
+
+ assert(len > 0);
+
+ p = (BNWORD64)BIGLITTLE(*--in,*in++) * k;
+ BIGLITTLE(*--out,*out++) = (BNWORD32)p;
+
+ while (--len) {
+ p = (BNWORD64)BIGLITTLE(*--in,*in++) * k + (BNWORD32)(p >> 32);
+ BIGLITTLE(*--out,*out++) = (BNWORD32)p;
+ }
+ BIGLITTLE(*--out,*out) = (BNWORD32)(p >> 32);
+}
+#else
+#error No 32x32 -> 64 multiply available for 32-bit bignum package
+#endif
+#endif /* lbnMulN1_32 */
+
+/*
+ * lbnMulAdd1_32: Multiply an n-word input by a 1-word input and add the
+ * low n words of the product to the destination. *Returns the n+1st word
+ * of the product.* (That turns out to be more convenient than adding
+ * it into the destination and dealing with a possible unit carry out
+ * of *that*.) This uses either the mul32_ppmma and mul32_ppmmaa macros,
+ * or C multiplication with the BNWORD64 type.
+ *
+ * If you're going to write assembly primitives, this is the one to
+ * start with. It is by far the most commonly called function.
+ */
+#ifndef lbnMulAdd1_32
+#if defined(mul32_ppmm)
+BNWORD32
+lbnMulAdd1_32(BNWORD32 *out, BNWORD32 const *in, unsigned len, BNWORD32 k)
+{
+ BNWORD32 prod, carry, carryin;
+
+ assert(len > 0);
+
+ BIG(--out;--in;);
+ carryin = *out;
+ mul32_ppmma(carry, *out, *in, k, carryin);
+ LITTLE(out++;in++;)
+
+ while (--len) {
+ BIG(--out;--in;);
+ carryin = carry;
+ mul32_ppmmaa(carry, prod, *in, k, carryin, *out);
+ *out = prod;
+ LITTLE(out++;in++;)
+ }
+
+ return carry;
+}
+#elif defined(BNWORD64)
+BNWORD32
+lbnMulAdd1_32(BNWORD32 *out, BNWORD32 const *in, unsigned len, BNWORD32 k)
+{
+ BNWORD64 p;
+
+ assert(len > 0);
+
+ p = (BNWORD64)BIGLITTLE(*--in,*in++) * k + BIGLITTLE(*--out,*out);
+ BIGLITTLE(*out,*out++) = (BNWORD32)p;
+
+ while (--len) {
+ p = (BNWORD64)BIGLITTLE(*--in,*in++) * k +
+ (BNWORD32)(p >> 32) + BIGLITTLE(*--out,*out);
+ BIGLITTLE(*out,*out++) = (BNWORD32)p;
+ }
+
+ return (BNWORD32)(p >> 32);
+}
+#else
+#error No 32x32 -> 64 multiply available for 32-bit bignum package
+#endif
+#endif /* lbnMulAdd1_32 */
+
+/*
+ * lbnMulSub1_32: Multiply an n-word input by a 1-word input and subtract the
+ * n-word product from the destination. Returns the n+1st word of the product.
+ * This uses either the mul32_ppmm and mul32_ppmma macros, or
+ * C multiplication with the BNWORD64 type.
+ *
+ * This is rather uglier than adding, but fortunately it's only used in
+ * division which is not used too heavily.
+ */
+#ifndef lbnMulSub1_32
+#if defined(mul32_ppmm)
+BNWORD32
+lbnMulSub1_32(BNWORD32 *out, BNWORD32 const *in, unsigned len, BNWORD32 k)
+{
+ BNWORD32 prod, carry, carryin;
+
+ assert(len > 0);
+
+ BIG(--in;)
+ mul32_ppmm(carry, prod, *in, k);
+ LITTLE(in++;)
+ carry += (BIGLITTLE(*--out,*out++) -= prod) > (BNWORD32)~prod;
+
+ while (--len) {
+ BIG(--in;);
+ carryin = carry;
+ mul32_ppmma(carry, prod, *in, k, carryin);
+ LITTLE(in++;)
+ carry += (BIGLITTLE(*--out,*out++) -= prod) > (BNWORD32)~prod;
+ }
+
+ return carry;
+}
+#elif defined(BNWORD64)
+BNWORD32
+lbnMulSub1_32(BNWORD32 *out, BNWORD32 const *in, unsigned len, BNWORD32 k)
+{
+ BNWORD64 p;
+ BNWORD32 carry, t;
+
+ assert(len > 0);
+
+ p = (BNWORD64)BIGLITTLE(*--in,*in++) * k;
+ t = BIGLITTLE(*--out,*out);
+ carry = (BNWORD32)(p>>32) + ((BIGLITTLE(*out,*out++)=t-(BNWORD32)p) > t);
+
+ while (--len) {
+ p = (BNWORD64)BIGLITTLE(*--in,*in++) * k + carry;
+ t = BIGLITTLE(*--out,*out);
+ carry = (BNWORD32)(p>>32) +
+ ( (BIGLITTLE(*out,*out++)=t-(BNWORD32)p) > t );
+ }
+
+ return carry;
+}
+#else
+#error No 32x32 -> 64 multiply available for 32-bit bignum package
+#endif
+#endif /* !lbnMulSub1_32 */
+
+/*
+ * Shift n words left "shift" bits. 0 < shift < 32. Returns the
+ * carry, any bits shifted off the left-hand side (0 <= carry < 2^shift).
+ */
+#ifndef lbnLshift_32
+BNWORD32
+lbnLshift_32(BNWORD32 *num, unsigned len, unsigned shift)
+{
+ BNWORD32 x, carry;
+
+ assert(shift > 0);
+ assert(shift < 32);
+
+ carry = 0;
+ while (len--) {
+ BIG(--num;)
+ x = *num;
+ *num = (x<<shift) | carry;
+ LITTLE(num++;)
+ carry = x >> (32-shift);
+ }
+ return carry;
+}
+#endif /* !lbnLshift_32 */
+
+/*
+ * An optimized version of the above, for shifts of 1.
+ * Some machines can use add-with-carry tricks for this.
+ */
+#ifndef lbnDouble_32
+BNWORD32
+lbnDouble_32(BNWORD32 *num, unsigned len)
+{
+ BNWORD32 x, carry;
+
+ carry = 0;
+ while (len--) {
+ BIG(--num;)
+ x = *num;
+ *num = (x<<1) | carry;
+ LITTLE(num++;)
+ carry = x >> (32-1);
+ }
+ return carry;
+}
+#endif /* !lbnDouble_32 */
+
+/*
+ * Shift n words right "shift" bits. 0 < shift < 32. Returns the
+ * carry, any bits shifted off the right-hand side (0 <= carry < 2^shift).
+ */
+#ifndef lbnRshift_32
+BNWORD32
+lbnRshift_32(BNWORD32 *num, unsigned len, unsigned shift)
+{
+ BNWORD32 x, carry = 0;
+
+ assert(shift > 0);
+ assert(shift < 32);
+
+ BIGLITTLE(num -= len, num += len);
+
+ while (len--) {
+ LITTLE(--num;)
+ x = *num;
+ *num = (x>>shift) | carry;
+ BIG(num++;)
+ carry = x << (32-shift);
+ }
+ return carry >> (32-shift);
+}
+#endif /* !lbnRshift_32 */
+
+/*
+ * Multiply two numbers of the given lengths. prod and num2 may overlap,
+ * provided that the low len1 bits of prod are free. (This corresponds
+ * nicely to the place the result is returned from lbnMontReduce_32.)
+ *
+ * TODO: Use Karatsuba multiply. The overlap constraints may have
+ * to get rewhacked.
+ */
+#ifndef lbnMul_32
+void
+lbnMul_32(BNWORD32 *prod, BNWORD32 const *num1, unsigned len1,
+ BNWORD32 const *num2, unsigned len2)
+{
+ /* Special case of zero */
+ if (!len1 || !len2) {
+ lbnZero_32(prod, len1+len2);
+ return;
+ }
+
+ /* Multiply first word */
+ lbnMulN1_32(prod, num1, len1, BIGLITTLE(*--num2,*num2++));
+
+ /*
+ * Add in subsequent words, storing the most significant word,
+ * which is new each time.
+ */
+ while (--len2) {
+ BIGLITTLE(--prod,prod++);
+ BIGLITTLE(*(prod-len1-1),*(prod+len1)) =
+ lbnMulAdd1_32(prod, num1, len1, BIGLITTLE(*--num2,*num2++));
+ }
+}
+#endif /* !lbnMul_32 */
+
+/*
+ * lbnMulX_32 is a square multiply - both inputs are the same length.
+ * It's normally just a macro wrapper around the general multiply,
+ * but might be implementable in assembly more efficiently (such as
+ * when product scanning).
+ */
+#ifndef lbnMulX_32
+#if defined(BNWORD64) && PRODUCT_SCAN
+/*
+ * Test code to see whether product scanning is any faster. It seems
+ * to make the C code slower, so PRODUCT_SCAN is not defined.
+ */
+static void
+lbnMulX_32(BNWORD32 *prod, BNWORD32 const *num1, BNWORD32 const *num2,
+ unsigned len)
+{
+ BNWORD64 x, y;
+ BNWORD32 const *p1, *p2;
+ unsigned carry;
+ unsigned i, j;
+
+ /* Special case of zero */
+ if (!len)
+ return;
+
+ x = (BNWORD64)BIGLITTLE(num1[-1] * num2[-1], num1[0] * num2[0]);
+ BIGLITTLE(*--prod, *prod++) = (BNWORD32)x;
+ x >>= 32;
+
+ for (i = 1; i < len; i++) {
+ carry = 0;
+ p1 = num1;
+ p2 = BIGLITTLE(num2-i-1,num2+i+1);
+ for (j = 0; j <= i; j++) {
+ BIG(y = (BNWORD64)*--p1 * *p2++;)
+ LITTLE(y = (BNWORD64)*p1++ * *--p2;)
+ x += y;
+ carry += (x < y);
+ }
+ BIGLITTLE(*--prod,*prod++) = (BNWORD32)x;
+ x = (x >> 32) | (BNWORD64)carry << 32;
+ }
+ for (i = 1; i < len; i++) {
+ carry = 0;
+ p1 = BIGLITTLE(num1-i,num1+i);
+ p2 = BIGLITTLE(num2-len,num2+len);
+ for (j = i; j < len; j++) {
+ BIG(y = (BNWORD64)*--p1 * *p2++;)
+ LITTLE(y = (BNWORD64)*p1++ * *--p2;)
+ x += y;
+ carry += (x < y);
+ }
+ BIGLITTLE(*--prod,*prod++) = (BNWORD32)x;
+ x = (x >> 32) | (BNWORD64)carry << 32;
+ }
+
+ BIGLITTLE(*--prod,*prod) = (BNWORD32)x;
+}
+#else /* !defined(BNWORD64) || !PRODUCT_SCAN */
+/* Default trivial macro definition */
+#define lbnMulX_32(prod, num1, num2, len) lbnMul_32(prod, num1, len, num2, len)
+#endif /* !defined(BNWORD64) || !PRODUCT_SCAN */
+#endif /* !lbmMulX_32 */
+
+#if !defined(lbnMontMul_32) && defined(BNWORD64) && PRODUCT_SCAN
+/*
+ * Test code for product-scanning multiply. This seems to slow the C
+ * code down rather than speed it up.
+ * This does a multiply and Montgomery reduction together, using the
+ * same loops. The outer loop scans across the product, twice.
+ * The first pass computes the low half of the product and the
+ * Montgomery multipliers. These are stored in the product array,
+ * which contains no data as of yet. x and carry add up the columns
+ * and propagate carries forward.
+ *
+ * The second half multiplies the upper half, adding in the modulus
+ * times the Montgomery multipliers. The results of this multiply
+ * are stored.
+ */
+static void
+lbnMontMul_32(BNWORD32 *prod, BNWORD32 const *num1, BNWORD32 const *num2,
+ BNWORD32 const *mod, unsigned len, BNWORD32 inv)
+{
+ BNWORD64 x, y;
+ BNWORD32 const *p1, *p2, *pm;
+ BNWORD32 *pp;
+ BNWORD32 t;
+ unsigned carry;
+ unsigned i, j;
+
+ /* Special case of zero */
+ if (!len)
+ return;
+
+ /*
+ * This computes directly into the high half of prod, so just
+ * shift the pointer and consider prod only "len" elements long
+ * for the rest of the code.
+ */
+ BIGLITTLE(prod -= len, prod += len);
+
+ /* Pass 1 - compute Montgomery multipliers */
+ /* First iteration can have certain simplifications. */
+ x = (BNWORD64)BIGLITTLE(num1[-1] * num2[-1], num1[0] * num2[0]);
+ BIGLITTLE(prod[-1], prod[0]) = t = inv * (BNWORD32)x;
+ y = (BNWORD64)t * BIGLITTLE(mod[-1],mod[0]);
+ x += y;
+ /* Note: GCC 2.6.3 has a bug if you try to eliminate "carry" */
+ carry = (x < y);
+ assert((BNWORD32)x == 0);
+ x = x >> 32 | (BNWORD64)carry << 32;
+
+ for (i = 1; i < len; i++) {
+ carry = 0;
+ p1 = num1;
+ p2 = BIGLITTLE(num2-i-1,num2+i+1);
+ pp = prod;
+ pm = BIGLITTLE(mod-i-1,mod+i+1);
+ for (j = 0; j < i; j++) {
+ y = (BNWORD64)BIGLITTLE(*--p1 * *p2++, *p1++ * *--p2);
+ x += y;
+ carry += (x < y);
+ y = (BNWORD64)BIGLITTLE(*--pp * *pm++, *pp++ * *--pm);
+ x += y;
+ carry += (x < y);
+ }
+ y = (BNWORD64)BIGLITTLE(p1[-1] * p2[0], p1[0] * p2[-1]);
+ x += y;
+ carry += (x < y);
+ assert(BIGLITTLE(pp == prod-i, pp == prod+i));
+ BIGLITTLE(pp[-1], pp[0]) = t = inv * (BNWORD32)x;
+ assert(BIGLITTLE(pm == mod-1, pm == mod+1));
+ y = (BNWORD64)t * BIGLITTLE(pm[0],pm[-1]);
+ x += y;
+ carry += (x < y);
+ assert((BNWORD32)x == 0);
+ x = x >> 32 | (BNWORD64)carry << 32;
+ }
+
+ /* Pass 2 - compute reduced product and store */
+ for (i = 1; i < len; i++) {
+ carry = 0;
+ p1 = BIGLITTLE(num1-i,num1+i);
+ p2 = BIGLITTLE(num2-len,num2+len);
+ pm = BIGLITTLE(mod-i,mod+i);
+ pp = BIGLITTLE(prod-len,prod+len);
+ for (j = i; j < len; j++) {
+ y = (BNWORD64)BIGLITTLE(*--p1 * *p2++, *p1++ * *--p2);
+ x += y;
+ carry += (x < y);
+ y = (BNWORD64)BIGLITTLE(*--pm * *pp++, *pm++ * *--pp);
+ x += y;
+ carry += (x < y);
+ }
+ assert(BIGLITTLE(pm == mod-len, pm == mod+len));
+ assert(BIGLITTLE(pp == prod-i, pp == prod+i));
+ BIGLITTLE(pp[0],pp[-1]) = (BNWORD32)x;
+ x = (x >> 32) | (BNWORD64)carry << 32;
+ }
+
+ /* Last round of second half, simplified. */
+ BIGLITTLE(*(prod-len),*(prod+len-1)) = (BNWORD32)x;
+ carry = (x >> 32);
+
+ while (carry)
+ carry -= lbnSubN_32(prod, mod, len);
+ while (lbnCmp_32(prod, mod, len) >= 0)
+ (void)lbnSubN_32(prod, mod, len);
+}
+/* Suppress later definition */
+#define lbnMontMul_32 lbnMontMul_32
+#endif
+
+#if !defined(lbnSquare_32) && defined(BNWORD64) && PRODUCT_SCAN
+/*
+ * Trial code for product-scanning squaring. This seems to slow the C
+ * code down rather than speed it up.
+ */
+void
+lbnSquare_32(BNWORD32 *prod, BNWORD32 const *num, unsigned len)
+{
+ BNWORD64 x, y, z;
+ BNWORD32 const *p1, *p2;
+ unsigned carry;
+ unsigned i, j;
+
+ /* Special case of zero */
+ if (!len)
+ return;
+
+ /* Word 0 of product */
+ x = (BNWORD64)BIGLITTLE(num[-1] * num[-1], num[0] * num[0]);
+ BIGLITTLE(*--prod, *prod++) = (BNWORD32)x;
+ x >>= 32;
+
+ /* Words 1 through len-1 */
+ for (i = 1; i < len; i++) {
+ carry = 0;
+ y = 0;
+ p1 = num;
+ p2 = BIGLITTLE(num-i-1,num+i+1);
+ for (j = 0; j < (i+1)/2; j++) {
+ BIG(z = (BNWORD64)*--p1 * *p2++;)
+ LITTLE(z = (BNWORD64)*p1++ * *--p2;)
+ y += z;
+ carry += (y < z);
+ }
+ y += z = y;
+ carry += carry + (y < z);
+ if ((i & 1) == 0) {
+ assert(BIGLITTLE(--p1 == p2, p1 == --p2));
+ BIG(z = (BNWORD64)*p2 * *p2;)
+ LITTLE(z = (BNWORD64)*p1 * *p1;)
+ y += z;
+ carry += (y < z);
+ }
+ x += y;
+ carry += (x < y);
+ BIGLITTLE(*--prod,*prod++) = (BNWORD32)x;
+ x = (x >> 32) | (BNWORD64)carry << 32;
+ }
+ /* Words len through 2*len-2 */
+ for (i = 1; i < len; i++) {
+ carry = 0;
+ y = 0;
+ p1 = BIGLITTLE(num-i,num+i);
+ p2 = BIGLITTLE(num-len,num+len);
+ for (j = 0; j < (len-i)/2; j++) {
+ BIG(z = (BNWORD64)*--p1 * *p2++;)
+ LITTLE(z = (BNWORD64)*p1++ * *--p2;)
+ y += z;
+ carry += (y < z);
+ }
+ y += z = y;
+ carry += carry + (y < z);
+ if ((len-i) & 1) {
+ assert(BIGLITTLE(--p1 == p2, p1 == --p2));
+ BIG(z = (BNWORD64)*p2 * *p2;)
+ LITTLE(z = (BNWORD64)*p1 * *p1;)
+ y += z;
+ carry += (y < z);
+ }
+ x += y;
+ carry += (x < y);
+ BIGLITTLE(*--prod,*prod++) = (BNWORD32)x;
+ x = (x >> 32) | (BNWORD64)carry << 32;
+ }
+
+ /* Word 2*len-1 */
+ BIGLITTLE(*--prod,*prod) = (BNWORD32)x;
+}
+/* Suppress later definition */
+#define lbnSquare_32 lbnSquare_32
+#endif
+
+/*
+ * Square a number, using optimized squaring to reduce the number of
+ * primitive multiples that are executed. There may not be any
+ * overlap of the input and output.
+ *
+ * Technique: Consider the partial products in the multiplication
+ * of "abcde" by itself:
+ *
+ * a b c d e
+ * * a b c d e
+ * ==================
+ * ae be ce de ee
+ * ad bd cd dd de
+ * ac bc cc cd ce
+ * ab bb bc bd be
+ * aa ab ac ad ae
+ *
+ * Note that everything above the main diagonal:
+ * ae be ce de = (abcd) * e
+ * ad bd cd = (abc) * d
+ * ac bc = (ab) * c
+ * ab = (a) * b
+ *
+ * is a copy of everything below the main diagonal:
+ * de
+ * cd ce
+ * bc bd be
+ * ab ac ad ae
+ *
+ * Thus, the sum is 2 * (off the diagonal) + diagonal.
+ *
+ * This is accumulated beginning with the diagonal (which
+ * consist of the squares of the digits of the input), which is then
+ * divided by two, the off-diagonal added, and multiplied by two
+ * again. The low bit is simply a copy of the low bit of the
+ * input, so it doesn't need special care.
+ *
+ * TODO: Merge the shift by 1 with the squaring loop.
+ * TODO: Use Karatsuba. (a*W+b)^2 = a^2 * (W^2+W) + b^2 * (W+1) - (a-b)^2 * W.
+ */
+#ifndef lbnSquare_32
+void
+lbnSquare_32(BNWORD32 *prod, BNWORD32 const *num, unsigned len)
+{
+ BNWORD32 t;
+ BNWORD32 *prodx = prod; /* Working copy of the argument */
+ BNWORD32 const *numx = num; /* Working copy of the argument */
+ unsigned lenx = len; /* Working copy of the argument */
+
+ if (!len)
+ return;
+
+ /* First, store all the squares */
+ while (lenx--) {
+#ifdef mul32_ppmm
+ BNWORD32 ph, pl;
+ t = BIGLITTLE(*--numx,*numx++);
+ mul32_ppmm(ph,pl,t,t);
+ BIGLITTLE(*--prodx,*prodx++) = pl;
+ BIGLITTLE(*--prodx,*prodx++) = ph;
+#elif defined(BNWORD64) /* use BNWORD64 */
+ BNWORD64 p;
+ t = BIGLITTLE(*--numx,*numx++);
+ p = (BNWORD64)t * t;
+ BIGLITTLE(*--prodx,*prodx++) = (BNWORD32)p;
+ BIGLITTLE(*--prodx,*prodx++) = (BNWORD32)(p>>32);
+#else /* Use lbnMulN1_32 */
+ t = BIGLITTLE(numx[-1],*numx);
+ lbnMulN1_32(prodx, numx, 1, t);
+ BIGLITTLE(--numx,numx++);
+ BIGLITTLE(prodx -= 2, prodx += 2);
+#endif
+ }
+ /* Then, shift right 1 bit */
+ (void)lbnRshift_32(prod, 2*len, 1);
+
+ /* Then, add in the off-diagonal sums */
+ lenx = len;
+ numx = num;
+ prodx = prod;
+ while (--lenx) {
+ t = BIGLITTLE(*--numx,*numx++);
+ BIGLITTLE(--prodx,prodx++);
+ t = lbnMulAdd1_32(prodx, numx, lenx, t);
+ lbnAdd1_32(BIGLITTLE(prodx-lenx,prodx+lenx), lenx+1, t);
+ BIGLITTLE(--prodx,prodx++);
+ }
+
+ /* Shift it back up */
+ lbnDouble_32(prod, 2*len);
+
+ /* And set the low bit appropriately */
+ BIGLITTLE(prod[-1],prod[0]) |= BIGLITTLE(num[-1],num[0]) & 1;
+}
+#endif /* !lbnSquare_32 */
+
+/*
+ * lbnNorm_32 - given a number, return a modified length such that the
+ * most significant digit is non-zero. Zero-length input is okay.
+ */
+#ifndef lbnNorm_32
+unsigned
+lbnNorm_32(BNWORD32 const *num, unsigned len)
+{
+ BIGLITTLE(num -= len,num += len);
+ while (len && BIGLITTLE(*num++,*--num) == 0)
+ --len;
+ return len;
+}
+#endif /* lbnNorm_32 */
+
+/*
+ * lbnBits_32 - return the number of significant bits in the array.
+ * It starts by normalizing the array. Zero-length input is okay.
+ * Then assuming there's anything to it, it fetches the high word,
+ * generates a bit length by multiplying the word length by 32, and
+ * subtracts off 32/2, 32/4, 32/8, ... bits if the high bits are clear.
+ */
+#ifndef lbnBits_32
+unsigned
+lbnBits_32(BNWORD32 const *num, unsigned len)
+{
+ BNWORD32 t;
+ unsigned i;
+
+ len = lbnNorm_32(num, len);
+ if (len) {
+ t = BIGLITTLE(*(num-len),*(num+(len-1)));
+ assert(t);
+ len *= 32;
+ i = 32/2;
+ do {
+ if (t >> i)
+ t >>= i;
+ else
+ len -= i;
+ } while ((i /= 2) != 0);
+ }
+ return len;
+}
+#endif /* lbnBits_32 */
+
+/*
+ * If defined, use hand-rolled divide rather than compiler's native.
+ * If the machine doesn't do it in line, the manual code is probably
+ * faster, since it can assume normalization and the fact that the
+ * quotient will fit into 32 bits, which a general 64-bit divide
+ * in a compiler's run-time library can't do.
+ */
+#ifndef BN_SLOW_DIVIDE_64
+/* Assume that divisors of more than thirty-two bits are slow */
+#define BN_SLOW_DIVIDE_64 (64 > 0x20)
+#endif
+
+/*
+ * Return (nh<<32|nl) % d, and place the quotient digit into *q.
+ * It is guaranteed that nh < d, and that d is normalized (with its high
+ * bit set). If we have a double-width type, it's easy. If not, ooh,
+ * yuk!
+ */
+#ifndef lbnDiv21_32
+#if defined(BNWORD64) && !BN_SLOW_DIVIDE_64
+BNWORD32
+lbnDiv21_32(BNWORD32 *q, BNWORD32 nh, BNWORD32 nl, BNWORD32 d)
+{
+ BNWORD64 n = (BNWORD64)nh << 32 | nl;
+
+ /* Divisor must be normalized */
+ assert(d >> (32-1) == 1);
+
+ *q = n / d;
+ return n % d;
+}
+#else
+/*
+ * This is where it gets ugly.
+ *
+ * Do the division in two halves, using Algorithm D from section 4.3.1
+ * of Knuth. Note Theorem B from that section, that the quotient estimate
+ * is never more than the true quotient, and is never more than two
+ * too low.
+ *
+ * The mapping onto conventional long division is (everything a half word):
+ * _____________qh___ql_
+ * dh dl ) nh.h nh.l nl.h nl.l
+ * - (qh * d)
+ * -----------
+ * rrrr rrrr nl.l
+ * - (ql * d)
+ * -----------
+ * rrrr rrrr
+ *
+ * The implicit 3/2-digit d*qh and d*ql subtractors are computed this way:
+ * First, estimate a q digit so that nh/dh works. Subtracting qh*dh from
+ * the (nh.h nh.l) list leaves a 1/2-word remainder r. Then compute the
+ * low part of the subtractor, qh * dl. This also needs to be subtracted
+ * from (nh.h nh.l nl.h) to get the final remainder. So we take the
+ * remainder, which is (nh.h nh.l) - qh*dl, shift it and add in nl.h, and
+ * try to subtract qh * dl from that. Since the remainder is 1/2-word
+ * long, shifting and adding nl.h results in a single word r.
+ * It is possible that the remainder we're working with, r, is less than
+ * the product qh * dl, if we estimated qh too high. The estimation
+ * technique can produce a qh that is too large (never too small), leading
+ * to r which is too small. In that case, decrement the digit qh, add
+ * shifted dh to r (to correct for that error), and subtract dl from the
+ * product we're comparing r with. That's the "correct" way to do it, but
+ * just adding dl to r instead of subtracting it from the product is
+ * equivalent and a lot simpler. You just have to watch out for overflow.
+ *
+ * The process is repeated with (rrrr rrrr nl.l) for the low digit of the
+ * quotient ql.
+ *
+ * The various uses of 32/2 for shifts are because of the note about
+ * automatic editing of this file at the very top of the file.
+ */
+#define highhalf(x) ( (x) >> 32/2 )
+#define lowhalf(x) ( (x) & (((BNWORD32)1 << 32/2)-1) )
+BNWORD32
+lbnDiv21_32(BNWORD32 *q, BNWORD32 nh, BNWORD32 nl, BNWORD32 d)
+{
+ BNWORD32 dh = highhalf(d), dl = lowhalf(d);
+ BNWORD32 qh, ql, prod, r;
+
+ /* Divisor must be normalized */
+ assert((d >> (32-1)) == 1);
+
+ /* Do first half-word of division */
+ qh = nh / dh;
+ r = nh % dh;
+ prod = qh * dl;
+
+ /*
+ * Add next half-word of numerator to remainder and correct.
+ * qh may be up to two too large.
+ */
+ r = (r << (32/2)) | highhalf(nl);
+ if (r < prod) {
+ --qh; r += d;
+ if (r >= d && r < prod) {
+ --qh; r += d;
+ }
+ }
+ r -= prod;
+
+ /* Do second half-word of division */
+ ql = r / dh;
+ r = r % dh;
+ prod = ql * dl;
+
+ r = (r << (32/2)) | lowhalf(nl);
+ if (r < prod) {
+ --ql; r += d;
+ if (r >= d && r < prod) {
+ --ql; r += d;
+ }
+ }
+ r -= prod;
+
+ *q = (qh << (32/2)) | ql;
+
+ return r;
+}
+#endif
+#endif /* lbnDiv21_32 */
+
+
+/*
+ * In the division functions, the dividend and divisor are referred to
+ * as "n" and "d", which stand for "numerator" and "denominator".
+ *
+ * The quotient is (nlen-dlen+1) digits long. It may be overlapped with
+ * the high (nlen-dlen) words of the dividend, but one extra word is needed
+ * on top to hold the top word.
+ */
+
+/*
+ * Divide an n-word number by a 1-word number, storing the remainder
+ * and n-1 words of the n-word quotient. The high word is returned.
+ * It IS legal for rem to point to the same address as n, and for
+ * q to point one word higher.
+ *
+ * TODO: If BN_SLOW_DIVIDE_64, add a divnhalf_32 which uses 32-bit
+ * dividends if the divisor is half that long.
+ * TODO: Shift the dividend on the fly to avoid the last division and
+ * instead have a remainder that needs shifting.
+ * TODO: Use reciprocals rather than dividing.
+ */
+#ifndef lbnDiv1_32
+BNWORD32
+lbnDiv1_32(BNWORD32 *q, BNWORD32 *rem, BNWORD32 const *n, unsigned len,
+ BNWORD32 d)
+{
+ unsigned shift;
+ unsigned xlen;
+ BNWORD32 r;
+ BNWORD32 qhigh;
+
+ assert(len > 0);
+ assert(d);
+
+ if (len == 1) {
+ r = *n;
+ *rem = r%d;
+ return r/d;
+ }
+
+ shift = 0;
+ r = d;
+ xlen = 32/2;
+ do {
+ if (r >> xlen)
+ r >>= xlen;
+ else
+ shift += xlen;
+ } while ((xlen /= 2) != 0);
+ assert((d >> (32-1-shift)) == 1);
+ d <<= shift;
+
+ BIGLITTLE(q -= len-1,q += len-1);
+ BIGLITTLE(n -= len,n += len);
+
+ r = BIGLITTLE(*n++,*--n);
+ if (r < d) {
+ qhigh = 0;
+ } else {
+ qhigh = r/d;
+ r %= d;
+ }
+
+ xlen = len;
+ while (--xlen)
+ r = lbnDiv21_32(BIGLITTLE(q++,--q), r, BIGLITTLE(*n++,*--n), d);
+
+ /*
+ * Final correction for shift - shift the quotient up "shift"
+ * bits, and merge in the extra bits of quotient. Then reduce
+ * the final remainder mod the real d.
+ */
+ if (shift) {
+ d >>= shift;
+ qhigh = (qhigh << shift) | lbnLshift_32(q, len-1, shift);
+ BIGLITTLE(q[-1],*q) |= r/d;
+ r %= d;
+ }
+ *rem = r;
+
+ return qhigh;
+}
+#endif
+
+/*
+ * This function performs a "quick" modulus of a number with a divisor
+ * d which is guaranteed to be at most sixteen bits, i.e. less than 65536.
+ * This applies regardless of the word size the library is compiled with.
+ *
+ * This function is important to prime generation, for sieving.
+ */
+#ifndef lbnModQ_32
+/* If there's a custom lbnMod21_32, no normalization needed */
+#ifdef lbnMod21_32
+unsigned
+lbnModQ_32(BNWORD32 const *n, unsigned len, unsigned d)
+{
+ unsigned i, shift;
+ BNWORD32 r;
+
+ assert(len > 0);
+
+ BIGLITTLE(n -= len,n += len);
+
+ /* Try using a compare to avoid the first divide */
+ r = BIGLITTLE(*n++,*--n);
+ if (r >= d)
+ r %= d;
+ while (--len)
+ r = lbnMod21_32(r, BIGLITTLE(*n++,*--n), d);
+
+ return r;
+}
+#elif defined(BNWORD64) && !BN_SLOW_DIVIDE_64
+unsigned
+lbnModQ_32(BNWORD32 const *n, unsigned len, unsigned d)
+{
+ BNWORD32 r;
+
+ if (!--len)
+ return BIGLITTLE(n[-1],n[0]) % d;
+
+ BIGLITTLE(n -= len,n += len);
+ r = BIGLITTLE(n[-1],n[0]);
+
+ do {
+ r = (BNWORD32)((((BNWORD64)r<<32) | BIGLITTLE(*n++,*--n)) % d);
+ } while (--len);
+
+ return r;
+}
+#elif 32 >= 0x20
+/*
+ * If the single word size can hold 65535*65536, then this function
+ * is avilable.
+ */
+#ifndef highhalf
+#define highhalf(x) ( (x) >> 32/2 )
+#define lowhalf(x) ( (x) & ((1 << 32/2)-1) )
+#endif
+unsigned
+lbnModQ_32(BNWORD32 const *n, unsigned len, unsigned d)
+{
+ BNWORD32 r, x;
+
+ BIGLITTLE(n -= len,n += len);
+
+ r = BIGLITTLE(*n++,*--n);
+ while (--len) {
+ x = BIGLITTLE(*n++,*--n);
+ r = (r%d << 32/2) | highhalf(x);
+ r = (r%d << 32/2) | lowhalf(x);
+ }
+
+ return r%d;
+}
+#else
+/* Default case - use lbnDiv21_32 */
+unsigned
+lbnModQ_32(BNWORD32 const *n, unsigned len, unsigned d)
+{
+ unsigned i, shift;
+ BNWORD32 r;
+ BNWORD32 q;
+
+ assert(len > 0);
+
+ shift = 0;
+ r = d;
+ i = 32;
+ while (i /= 2) {
+ if (r >> i)
+ r >>= i;
+ else
+ shift += i;
+ }
+ assert(d >> (32-1-shift) == 1);
+ d <<= shift;
+
+ BIGLITTLE(n -= len,n += len);
+
+ r = BIGLITTLE(*n++,*--n);
+ if (r >= d)
+ r %= d;
+
+ while (--len)
+ r = lbnDiv21_32(&q, r, BIGLITTLE(*n++,*--n), d);
+
+ /*
+ * Final correction for shift - shift the quotient up "shift"
+ * bits, and merge in the extra bits of quotient. Then reduce
+ * the final remainder mod the real d.
+ */
+ if (shift)
+ r %= d >> shift;
+
+ return r;
+}
+#endif
+#endif /* lbnModQ_32 */
+
+/*
+ * Reduce n mod d and return the quotient. That is, find:
+ * q = n / d;
+ * n = n % d;
+ * d is altered during the execution of this subroutine by normalizing it.
+ * It must already have its most significant word non-zero; it is shifted
+ * so its most significant bit is non-zero.
+ *
+ * The quotient q is nlen-dlen+1 words long. To make it possible to
+ * overlap the quptient with the input (you can store it in the high dlen
+ * words), the high word of the quotient is *not* stored, but is returned.
+ * (If all you want is the remainder, you don't care about it, anyway.)
+ *
+ * This uses algorithm D from Knuth (4.3.1), except that we do binary
+ * (shift) normalization of the divisor. WARNING: This is hairy!
+ *
+ * This function is used for some modular reduction, but it is not used in
+ * the modular exponentiation loops; they use Montgomery form and the
+ * corresponding, more efficient, Montgomery reduction. This code
+ * is needed for the conversion to Montgomery form, however, so it
+ * has to be here and it might as well be reasonably efficient.
+ *
+ * The overall operation is as follows ("top" and "up" refer to the
+ * most significant end of the number; "bottom" and "down", the least):
+ *
+ * - Shift the divisor up until the most significant bit is set.
+ * - Shift the dividend up the same amount. This will produce the
+ * correct quotient, and the remainder can be recovered by shifting
+ * it back down the same number of bits. This may produce an overflow
+ * word, but the word is always strictly less than the most significant
+ * divisor word.
+ * - Estimate the first quotient digit qhat:
+ * - First take the top two words (one of which is the overflow) of the
+ * dividend and divide by the top word of the divisor:
+ * qhat = (nh,nm)/dh. This qhat is >= the correct quotient digit
+ * and, since dh is normalized, it is at most two over.
+ * - Second, correct by comparing the top three words. If
+ * (dh,dl) * qhat > (nh,nm,ml), decrease qhat and try again.
+ * The second iteration can be simpler because there can't be a third.
+ * The computation can be simplified by subtracting dh*qhat from
+ * both sides, suitably shifted. This reduces the left side to
+ * dl*qhat. On the right, (nh,nm)-dh*qhat is simply the
+ * remainder r from (nh,nm)%dh, so the right is (r,nl).
+ * This produces qhat that is almost always correct and at
+ * most (prob ~ 2/2^32) one too high.
+ * - Subtract qhat times the divisor (suitably shifted) from the dividend.
+ * If there is a borrow, qhat was wrong, so decrement it
+ * and add the divisor back in (once).
+ * - Store the final quotient digit qhat in the quotient array q.
+ *
+ * Repeat the quotient digit computation for successive digits of the
+ * quotient until the whole quotient has been computed. Then shift the
+ * divisor and the remainder down to correct for the normalization.
+ *
+ * TODO: Special case 2-word divisors.
+ * TODO: Use reciprocals rather than dividing.
+ */
+#ifndef divn_32
+BNWORD32
+lbnDiv_32(BNWORD32 *q, BNWORD32 *n, unsigned nlen, BNWORD32 *d, unsigned dlen)
+{
+ BNWORD32 nh,nm,nl; /* Top three words of the dividend */
+ BNWORD32 dh,dl; /* Top two words of the divisor */
+ BNWORD32 qhat; /* Extimate of quotient word */
+ BNWORD32 r; /* Remainder from quotient estimate division */
+ BNWORD32 qhigh; /* High word of quotient */
+ unsigned i; /* Temp */
+ unsigned shift; /* Bits shifted by normalization */
+ unsigned qlen = nlen-dlen; /* Size of quotient (less 1) */
+#ifdef mul32_ppmm
+ BNWORD32 t32;
+#elif defined(BNWORD64)
+ BNWORD64 t64;
+#else /* use lbnMulN1_32 */
+ BNWORD32 t2[2];
+#define t2high BIGLITTLE(t2[0],t2[1])
+#define t2low BIGLITTLE(t2[1],t2[0])
+#endif
+
+ assert(dlen);
+ assert(nlen >= dlen);
+
+ /*
+ * Special cases for short divisors. The general case uses the
+ * top top 2 digits of the divisor (d) to estimate a quotient digit,
+ * so it breaks if there are fewer digits available. Thus, we need
+ * special cases for a divisor of length 1. A divisor of length
+ * 2 can have a *lot* of administrivia overhead removed removed,
+ * so it's probably worth special-casing that case, too.
+ */
+ if (dlen == 1)
+ return lbnDiv1_32(q, BIGLITTLE(n-1,n), n, nlen,
+ BIGLITTLE(d[-1],d[0]));
+
+#if 0
+ /*
+ * @@@ This is not yet written... The general loop will do,
+ * albeit less efficiently
+ */
+ if (dlen == 2) {
+ /*
+ * divisor two digits long:
+ * use the 3/2 technique from Knuth, but we know
+ * it's exact.
+ */
+ dh = BIGLITTLE(d[-1],d[0]);
+ dl = BIGLITTLE(d[-2],d[1]);
+ shift = 0;
+ if ((sh & ((BNWORD32)1 << 32-1-shift)) == 0) {
+ do {
+ shift++;
+ } while (dh & (BNWORD32)1<<32-1-shift) == 0);
+ dh = dh << shift | dl >> (32-shift);
+ dl <<= shift;
+
+
+ }
+
+
+ for (shift = 0; (dh & (BNWORD32)1 << 32-1-shift)) == 0; shift++)
+ ;
+ if (shift) {
+ }
+ dh = dh << shift | dl >> (32-shift);
+ shift = 0;
+ while (dh
+ }
+#endif
+
+ dh = BIGLITTLE(*(d-dlen),*(d+(dlen-1)));
+ assert(dh);
+
+ /* Normalize the divisor */
+ shift = 0;
+ r = dh;
+ i = 32/2;
+ do {
+ if (r >> i)
+ r >>= i;
+ else
+ shift += i;
+ } while ((i /= 2) != 0);
+
+ nh = 0;
+ if (shift) {
+ lbnLshift_32(d, dlen, shift);
+ dh = BIGLITTLE(*(d-dlen),*(d+(dlen-1)));
+ nh = lbnLshift_32(n, nlen, shift);
+ }
+
+ /* Assert that dh is now normalized */
+ assert(dh >> (32-1));
+
+ /* Also get the second-most significant word of the divisor */
+ dl = BIGLITTLE(*(d-(dlen-1)),*(d+(dlen-2)));
+
+ /*
+ * Adjust pointers: n to point to least significant end of first
+ * first subtract, and q to one the most-significant end of the
+ * quotient array.
+ */
+ BIGLITTLE(n -= qlen,n += qlen);
+ BIGLITTLE(q -= qlen,q += qlen);
+
+ /* Fetch the most significant stored word of the dividend */
+ nm = BIGLITTLE(*(n-dlen),*(n+(dlen-1)));
+
+ /*
+ * Compute the first digit of the quotient, based on the
+ * first two words of the dividend (the most significant of which
+ * is the overflow word h).
+ */
+ if (nh) {
+ assert(nh < dh);
+ r = lbnDiv21_32(&qhat, nh, nm, dh);
+ } else if (nm >= dh) {
+ qhat = nm/dh;
+ r = nm % dh;
+ } else { /* Quotient is zero */
+ qhigh = 0;
+ goto divloop;
+ }
+
+ /* Now get the third most significant word of the dividend */
+ nl = BIGLITTLE(*(n-(dlen-1)),*(n+(dlen-2)));
+
+ /*
+ * Correct qhat, the estimate of quotient digit.
+ * qhat can only be high, and at most two words high,
+ * so the loop can be unrolled and abbreviated.
+ */
+#ifdef mul32_ppmm
+ mul32_ppmm(nm, t32, qhat, dl);
+ if (nm > r || (nm == r && t32 > nl)) {
+ /* Decrement qhat and adjust comparison parameters */
+ qhat--;
+ if ((r += dh) >= dh) {
+ nm -= (t32 < dl);
+ t32 -= dl;
+ if (nm > r || (nm == r && t32 > nl))
+ qhat--;
+ }
+ }
+#elif defined(BNWORD64)
+ t64 = (BNWORD64)qhat * dl;
+ if (t64 > ((BNWORD64)r << 32) + nl) {
+ /* Decrement qhat and adjust comparison parameters */
+ qhat--;
+ if ((r += dh) > dh) {
+ t64 -= dl;
+ if (t64 > ((BNWORD64)r << 32) + nl)
+ qhat--;
+ }
+ }
+#else /* Use lbnMulN1_32 */
+ lbnMulN1_32(BIGLITTLE(t2+2,t2), &dl, 1, qhat);
+ if (t2high > r || (t2high == r && t2low > nl)) {
+ /* Decrement qhat and adjust comparison parameters */
+ qhat--;
+ if ((r += dh) >= dh) {
+ t2high -= (t2low < dl);
+ t2low -= dl;
+ if (t2high > r || (t2high == r && t2low > nl))
+ qhat--;
+ }
+ }
+#endif
+
+ /* Do the multiply and subtract */
+ r = lbnMulSub1_32(n, d, dlen, qhat);
+ /* If there was a borrow, add back once. */
+ if (r > nh) { /* Borrow? */
+ (void)lbnAddN_32(n, d, dlen);
+ qhat--;
+ }
+
+ /* Remember the first quotient digit. */
+ qhigh = qhat;
+
+ /* Now, the main division loop: */
+divloop:
+ while (qlen--) {
+
+ /* Advance n */
+ nh = BIGLITTLE(*(n-dlen),*(n+(dlen-1)));
+ BIGLITTLE(++n,--n);
+ nm = BIGLITTLE(*(n-dlen),*(n+(dlen-1)));
+
+ if (nh == dh) {
+ qhat = ~(BNWORD32)0;
+ /* Optimized computation of r = (nh,nm) - qhat * dh */
+ r = nh + nm;
+ if (r < nh)
+ goto subtract;
+ } else {
+ assert(nh < dh);
+ r = lbnDiv21_32(&qhat, nh, nm, dh);
+ }
+
+ nl = BIGLITTLE(*(n-(dlen-1)),*(n+(dlen-2)));
+#ifdef mul32_ppmm
+ mul32_ppmm(nm, t32, qhat, dl);
+ if (nm > r || (nm == r && t32 > nl)) {
+ /* Decrement qhat and adjust comparison parameters */
+ qhat--;
+ if ((r += dh) >= dh) {
+ nm -= (t32 < dl);
+ t32 -= dl;
+ if (nm > r || (nm == r && t32 > nl))
+ qhat--;
+ }
+ }
+#elif defined(BNWORD64)
+ t64 = (BNWORD64)qhat * dl;
+ if (t64 > ((BNWORD64)r<<32) + nl) {
+ /* Decrement qhat and adjust comparison parameters */
+ qhat--;
+ if ((r += dh) >= dh) {
+ t64 -= dl;
+ if (t64 > ((BNWORD64)r << 32) + nl)
+ qhat--;
+ }
+ }
+#else /* Use lbnMulN1_32 */
+ lbnMulN1_32(BIGLITTLE(t2+2,t2), &dl, 1, qhat);
+ if (t2high > r || (t2high == r && t2low > nl)) {
+ /* Decrement qhat and adjust comparison parameters */
+ qhat--;
+ if ((r += dh) >= dh) {
+ t2high -= (t2low < dl);
+ t2low -= dl;
+ if (t2high > r || (t2high == r && t2low > nl))
+ qhat--;
+ }
+ }
+#endif
+
+ /*
+ * As a point of interest, note that it is not worth checking
+ * for qhat of 0 or 1 and installing special-case code. These
+ * occur with probability 2^-32, so spending 1 cycle to check
+ * for them is only worth it if we save more than 2^15 cycles,
+ * and a multiply-and-subtract for numbers in the 1024-bit
+ * range just doesn't take that long.
+ */
+subtract:
+ /*
+ * n points to the least significant end of the substring
+ * of n to be subtracted from. qhat is either exact or
+ * one too large. If the subtract gets a borrow, it was
+ * one too large and the divisor is added back in. It's
+ * a dlen+1 word add which is guaranteed to produce a
+ * carry out, so it can be done very simply.
+ */
+ r = lbnMulSub1_32(n, d, dlen, qhat);
+ if (r > nh) { /* Borrow? */
+ (void)lbnAddN_32(n, d, dlen);
+ qhat--;
+ }
+ /* Store the quotient digit */
+ BIGLITTLE(*q++,*--q) = qhat;
+ }
+ /* Tah dah! */
+
+ if (shift) {
+ lbnRshift_32(d, dlen, shift);
+ lbnRshift_32(n, dlen, shift);
+ }
+
+ return qhigh;
+}
+#endif
+
+/*
+ * Find the negative multiplicative inverse of x (x must be odd!) modulo 2^32.
+ *
+ * This just performs Newton's iteration until it gets the
+ * inverse. The initial estimate is always correct to 3 bits, and
+ * sometimes 4. The number of valid bits doubles each iteration.
+ * (To prove it, assume x * y == 1 (mod 2^n), and introduce a variable
+ * for the error mod 2^2n. x * y == 1 + k*2^n (mod 2^2n) and follow
+ * the iteration through.)
+ */
+#ifndef lbnMontInv1_32
+BNWORD32
+lbnMontInv1_32(BNWORD32 const x)
+{
+ BNWORD32 y = x, z;
+
+ assert(x & 1);
+
+ while ((z = x*y) != 1)
+ y *= 2 - z;
+ return -y;
+}
+#endif /* !lbnMontInv1_32 */
+
+#if defined(BNWORD64) && PRODUCT_SCAN
+/*
+ * Test code for product-scanning Montgomery reduction.
+ * This seems to slow the C code down rather than speed it up.
+ *
+ * The first loop computes the Montgomery multipliers, storing them over
+ * the low half of the number n.
+ *
+ * The second half multiplies the upper half, adding in the modulus
+ * times the Montgomery multipliers. The results of this multiply
+ * are stored.
+ */
+void
+lbnMontReduce_32(BNWORD32 *n, BNWORD32 const *mod, unsigned mlen, BNWORD32 inv)
+{
+ BNWORD64 x, y;
+ BNWORD32 const *pm;
+ BNWORD32 *pn;
+ BNWORD32 t;
+ unsigned carry;
+ unsigned i, j;
+
+ /* Special case of zero */
+ if (!mlen)
+ return;
+
+ /* Pass 1 - compute Montgomery multipliers */
+ /* First iteration can have certain simplifications. */
+ t = BIGLITTLE(n[-1],n[0]);
+ x = t;
+ t *= inv;
+ BIGLITTLE(n[-1], n[0]) = t;
+ x += (BNWORD64)t * BIGLITTLE(mod[-1],mod[0]); /* Can't overflow */
+ assert((BNWORD32)x == 0);
+ x = x >> 32;
+
+ for (i = 1; i < mlen; i++) {
+ carry = 0;
+ pn = n;
+ pm = BIGLITTLE(mod-i-1,mod+i+1);
+ for (j = 0; j < i; j++) {
+ y = (BNWORD64)BIGLITTLE(*--pn * *pm++, *pn++ * *--pm);
+ x += y;
+ carry += (x < y);
+ }
+ assert(BIGLITTLE(pn == n-i, pn == n+i));
+ y = t = BIGLITTLE(pn[-1], pn[0]);
+ x += y;
+ carry += (x < y);
+ BIGLITTLE(pn[-1], pn[0]) = t = inv * (BNWORD32)x;
+ assert(BIGLITTLE(pm == mod-1, pm == mod+1));
+ y = (BNWORD64)t * BIGLITTLE(pm[0],pm[-1]);
+ x += y;
+ carry += (x < y);
+ assert((BNWORD32)x == 0);
+ x = x >> 32 | (BNWORD64)carry << 32;
+ }
+
+ BIGLITTLE(n -= mlen, n += mlen);
+
+ /* Pass 2 - compute upper words and add to n */
+ for (i = 1; i < mlen; i++) {
+ carry = 0;
+ pm = BIGLITTLE(mod-i,mod+i);
+ pn = n;
+ for (j = i; j < mlen; j++) {
+ y = (BNWORD64)BIGLITTLE(*--pm * *pn++, *pm++ * *--pn);
+ x += y;
+ carry += (x < y);
+ }
+ assert(BIGLITTLE(pm == mod-mlen, pm == mod+mlen));
+ assert(BIGLITTLE(pn == n+mlen-i, pn == n-mlen+i));
+ y = t = BIGLITTLE(*(n-i),*(n+i-1));
+ x += y;
+ carry += (x < y);
+ BIGLITTLE(*(n-i),*(n+i-1)) = (BNWORD32)x;
+ x = (x >> 32) | (BNWORD64)carry << 32;
+ }
+
+ /* Last round of second half, simplified. */
+ t = BIGLITTLE(*(n-mlen),*(n+mlen-1));
+ x += t;
+ BIGLITTLE(*(n-mlen),*(n+mlen-1)) = (BNWORD32)x;
+ carry = (unsigned)(x >> 32);
+
+ while (carry)
+ carry -= lbnSubN_32(n, mod, mlen);
+ while (lbnCmp_32(n, mod, mlen) >= 0)
+ (void)lbnSubN_32(n, mod, mlen);
+}
+#define lbnMontReduce_32 lbnMontReduce_32
+#endif
+
+/*
+ * Montgomery reduce n, modulo mod. This reduces modulo mod and divides by
+ * 2^(32*mlen). Returns the result in the *top* mlen words of the argument n.
+ * This is ready for another multiplication using lbnMul_32.
+ *
+ * Montgomery representation is a very useful way to encode numbers when
+ * you're doing lots of modular reduction. What you do is pick a multiplier
+ * R which is relatively prime to the modulus and very easy to divide by.
+ * Since the modulus is odd, R is closen as a power of 2, so the division
+ * is a shift. In fact, it's a shift of an integral number of words,
+ * so the shift can be implicit - just drop the low-order words.
+ *
+ * Now, choose R *larger* than the modulus m, 2^(32*mlen). Then convert
+ * all numbers a, b, etc. to Montgomery form M(a), M(b), etc using the
+ * relationship M(a) = a*R mod m, M(b) = b*R mod m, etc. Note that:
+ * - The Montgomery form of a number depends on the modulus m.
+ * A fixed modulus m is assumed throughout this discussion.
+ * - Since R is relaitvely prime to m, multiplication by R is invertible;
+ * no information about the numbers is lost, they're just scrambled.
+ * - Adding (and subtracting) numbers in this form works just as usual.
+ * M(a+b) = (a+b)*R mod m = (a*R + b*R) mod m = (M(a) + M(b)) mod m
+ * - Multiplying numbers in this form produces a*b*R*R. The problem
+ * is to divide out the excess factor of R, modulo m as well as to
+ * reduce to the given length mlen. It turns out that this can be
+ * done *faster* than a normal divide, which is where the speedup
+ * in Montgomery division comes from.
+ *
+ * Normal reduction chooses a most-significant quotient digit q and then
+ * subtracts q*m from the number to be reduced. Choosing q is tricky
+ * and involved (just look at lbnDiv_32 to see!) and is usually
+ * imperfect, requiring a check for correction after the subtraction.
+ *
+ * Montgomery reduction *adds* a multiple of m to the *low-order* part
+ * of the number to be reduced. This multiple is chosen to make the
+ * low-order part of the number come out to zero. This can be done
+ * with no trickery or error using a precomputed inverse of the modulus.
+ * In this code, the "part" is one word, but any width can be used.
+ *
+ * Repeating this step sufficiently often results in a value which
+ * is a multiple of R (a power of two, remember) but is still (since
+ * the additions were to the low-order part and thus did not increase
+ * the value of the number being reduced very much) still not much
+ * larger than m*R. Then implicitly divide by R and subtract off
+ * m until the result is in the correct range.
+ *
+ * Since the low-order part being cancelled is less than R, the
+ * multiple of m added must have a multiplier which is at most R-1.
+ * Assuming that the input is at most m*R-1, the final number is
+ * at most m*(2*R-1)-1 = 2*m*R - m - 1, so subtracting m once from
+ * the high-order part, equivalent to subtracting m*R from the
+ * while number, produces a result which is at most m*R - m - 1,
+ * which divided by R is at most m-1.
+ *
+ * To convert *to* Montgomery form, you need a regular remainder
+ * routine, although you can just compute R*R (mod m) and do the
+ * conversion using Montgomery multiplication. To convert *from*
+ * Montgomery form, just Montgomery reduce the number to
+ * remove the extra factor of R.
+ *
+ * TODO: Change to a full inverse and use Karatsuba's multiplication
+ * rather than this word-at-a-time.
+ */
+#ifndef lbnMontReduce_32
+void
+lbnMontReduce_32(BNWORD32 *n, BNWORD32 const *mod, unsigned const mlen,
+ BNWORD32 inv)
+{
+ BNWORD32 t;
+ BNWORD32 c = 0;
+ unsigned len = mlen;
+
+ /* inv must be the negative inverse of mod's least significant word */
+ assert((BNWORD32)(inv * BIGLITTLE(mod[-1],mod[0])) == (BNWORD32)-1);
+
+ assert(len);
+
+ do {
+ t = lbnMulAdd1_32(n, mod, mlen, inv * BIGLITTLE(n[-1],n[0]));
+ c += lbnAdd1_32(BIGLITTLE(n-mlen,n+mlen), len, t);
+ BIGLITTLE(--n,++n);
+ } while (--len);
+
+ /*
+ * All that adding can cause an overflow past the modulus size,
+ * but it's unusual, and never by much, so a subtraction loop
+ * is the right way to deal with it.
+ * This subtraction happens infrequently - I've only ever seen it
+ * invoked once per reduction, and then just under 22.5% of the time.
+ */
+ while (c)
+ c -= lbnSubN_32(n, mod, mlen);
+ while (lbnCmp_32(n, mod, mlen) >= 0)
+ (void)lbnSubN_32(n, mod, mlen);
+}
+#endif /* !lbnMontReduce_32 */
+
+/*
+ * A couple of helpers that you might want to implement atomically
+ * in asm sometime.
+ */
+#ifndef lbnMontMul_32
+/*
+ * Multiply "num1" by "num2", modulo "mod", all of length "len", and
+ * place the result in the high half of "prod". "inv" is the inverse
+ * of the least-significant word of the modulus, modulo 2^32.
+ * This uses numbers in Montgomery form. Reduce using "len" and "inv".
+ *
+ * This is implemented as a macro to win on compilers that don't do
+ * inlining, since it's so trivial.
+ */
+#define lbnMontMul_32(prod, n1, n2, mod, len, inv) \
+ (lbnMulX_32(prod, n1, n2, len), lbnMontReduce_32(prod, mod, len, inv))
+#endif /* !lbnMontMul_32 */
+
+#ifndef lbnMontSquare_32
+/*
+ * Square "n", modulo "mod", both of length "len", and place the result
+ * in the high half of "prod". "inv" is the inverse of the least-significant
+ * word of the modulus, modulo 2^32.
+ * This uses numbers in Montgomery form. Reduce using "len" and "inv".
+ *
+ * This is implemented as a macro to win on compilers that don't do
+ * inlining, since it's so trivial.
+ */
+#define lbnMontSquare_32(prod, n, mod, len, inv) \
+ (lbnSquare_32(prod, n, len), lbnMontReduce_32(prod, mod, len, inv))
+
+#endif /* !lbnMontSquare_32 */
+
+/*
+ * Convert a number to Montgomery form - requires mlen + nlen words
+ * of memory in "n".
+ */
+void
+lbnToMont_32(BNWORD32 *n, unsigned nlen, BNWORD32 *mod, unsigned mlen)
+{
+ /* Move n up "mlen" words */
+ lbnCopy_32(BIGLITTLE(n-mlen,n+mlen), n, nlen);
+ lbnZero_32(n, mlen);
+ /* Do the division - dump the quotient in the high-order words */
+ (void)lbnDiv_32(BIGLITTLE(n-mlen,n+mlen), n, mlen+nlen, mod, mlen);
+}
+
+/*
+ * Convert from Montgomery form. Montgomery reduction is all that is
+ * needed.
+ */
+void
+lbnFromMont_32(BNWORD32 *n, BNWORD32 *mod, unsigned len)
+{
+ /* Zero the high words of n */
+ lbnZero_32(BIGLITTLE(n-len,n+len), len);
+ lbnMontReduce_32(n, mod, len, lbnMontInv1_32(mod[BIGLITTLE(-1,0)]));
+ /* Move n down len words */
+ lbnCopy_32(n, BIGLITTLE(n-len,n+len), len);
+}
+
+/*
+ * The windowed exponentiation algorithm, precomputes a table of odd
+ * powers of n up to 2^k. See the comment in bnExpMod_32 below for
+ * an explanation of how it actually works works.
+ *
+ * It takes 2^(k-1)-1 multiplies to compute the table, and (e-1)/(k+1)
+ * multiplies (on average) to perform the exponentiation. To minimize
+ * the sum, k must vary with e. The optimal window sizes vary with the
+ * exponent length. Here are some selected values and the boundary cases.
+ * (An underscore _ has been inserted into some of the numbers to ensure
+ * that magic strings like 32 do not appear in this table. It should be
+ * ignored.)
+ *
+ * At e = 1 bits, k=1 (0.000000) is best
+ * At e = 2 bits, k=1 (0.500000) is best
+ * At e = 4 bits, k=1 (1.500000) is best
+ * At e = 8 bits, k=2 (3.333333) < k=1 (3.500000)
+ * At e = 1_6 bits, k=2 (6.000000) is best
+ * At e = 26 bits, k=3 (9.250000) < k=2 (9.333333)
+ * At e = 3_2 bits, k=3 (10.750000) is best
+ * At e = 6_4 bits, k=3 (18.750000) is best
+ * At e = 82 bits, k=4 (23.200000) < k=3 (23.250000)
+ * At e = 128 bits, k=4 (3_2.400000) is best
+ * At e = 242 bits, k=5 (55.1_66667) < k=4 (55.200000)
+ * At e = 256 bits, k=5 (57.500000) is best
+ * At e = 512 bits, k=5 (100.1_66667) is best
+ * At e = 674 bits, k=6 (127.142857) < k=5 (127.1_66667)
+ * At e = 1024 bits, k=6 (177.142857) is best
+ * At e = 1794 bits, k=7 (287.125000) < k=6 (287.142857)
+ * At e = 2048 bits, k=7 (318.875000) is best
+ * At e = 4096 bits, k=7 (574.875000) is best
+ *
+ * The numbers in parentheses are the expected number of multiplications
+ * needed to do the computation. The normal russian-peasant modular
+ * exponentiation technique always uses (e-1)/2. For exponents as
+ * small as 192 bits (below the range of current factoring algorithms),
+ * half of the multiplies are eliminated, 45.2 as opposed to the naive
+ * 95.5. Counting the 191 squarings as 3/4 a multiply each (squaring
+ * proper is just over half of multiplying, but the Montgomery
+ * reduction in each case is also a multiply), that's 143.25
+ * multiplies, for totals of 188.45 vs. 238.75 - a 21% savings.
+ * For larger exponents (like 512 bits), it's 483.92 vs. 639.25, a
+ * 24.3% savings. It asymptotically approaches 25%.
+ *
+ * Um, actually there's a slightly more accurate way to count, which
+ * really is the average number of multiplies required, averaged
+ * uniformly over all 2^(e-1) e-bit numbers, from 2^(e-1) to (2^e)-1.
+ * It's based on the recurrence that for the last b bits, b <= k, at
+ * most one multiply is needed (and none at all 1/2^b of the time),
+ * while when b > k, the odds are 1/2 each way that the bit will be
+ * 0 (meaning no multiplies to reduce it to the b-1-bit case) and
+ * 1/2 that the bit will be 1, starting a k-bit window and requiring
+ * 1 multiply beyond the b-k-bit case. Since the most significant
+ * bit is always 1, a k-bit window always starts there, and that
+ * multiply is by 1, so it isn't a multiply at all. Thus, the
+ * number of multiplies is simply that needed for the last e-k bits.
+ * This recurrence produces:
+ *
+ * At e = 1 bits, k=1 (0.000000) is best
+ * At e = 2 bits, k=1 (0.500000) is best
+ * At e = 4 bits, k=1 (1.500000) is best
+ * At e = 6 bits, k=2 (2.437500) < k=1 (2.500000)
+ * At e = 8 bits, k=2 (3.109375) is best
+ * At e = 1_6 bits, k=2 (5.777771) is best
+ * At e = 24 bits, k=3 (8.437629) < k=2 (8.444444)
+ * At e = 3_2 bits, k=3 (10.437492) is best
+ * At e = 6_4 bits, k=3 (18.437500) is best
+ * At e = 81 bits, k=4 (22.6_40000) < k=3 (22.687500)
+ * At e = 128 bits, k=4 (3_2.040000) is best
+ * At e = 241 bits, k=5 (54.611111) < k=4 (54.6_40000)
+ * At e = 256 bits, k=5 (57.111111) is best
+ * At e = 512 bits, k=5 (99.777778) is best
+ * At e = 673 bits, k=6 (126.591837) < k=5 (126.611111)
+ * At e = 1024 bits, k=6 (176.734694) is best
+ * At e = 1793 bits, k=7 (286.578125) < k=6 (286.591837)
+ * At e = 2048 bits, k=7 (318.453125) is best
+ * At e = 4096 bits, k=7 (574.453125) is best
+ *
+ * This has the rollover points at 6, 24, 81, 241, 673 and 1793 instead
+ * of 8, 26, 82, 242, 674, and 1794. Not a very big difference.
+ * (The numbers past that are k=8 at 4609 and k=9 at 11521,
+ * vs. one more in each case for the approximation.)
+ *
+ * Given that exponents for which k>7 are useful are uncommon,
+ * a fixed size table for k <= 7 is used for simplicity.
+ *
+ * The basic number of squarings needed is e-1, although a k-bit
+ * window (for k > 1) can save, on average, k-2 of those, too.
+ * That savings currently isn't counted here. It would drive the
+ * crossover points slightly lower.
+ * (Actually, this win is also reduced in the DoubleExpMod case,
+ * meaning we'd have to split the tables. Except for that, the
+ * multiplies by powers of the two bases are independent, so
+ * the same logic applies to each as the single case.)
+ *
+ * Table entry i is the largest number of bits in an exponent to
+ * process with a window size of i+1. Entry 6 is the largest
+ * possible unsigned number, so the window will never be more
+ * than 7 bits, requiring 2^6 = 0x40 slots.
+ */
+#define BNEXPMOD_MAX_WINDOW 7
+static unsigned const bnExpModThreshTable[BNEXPMOD_MAX_WINDOW] = {
+ 5, 23, 80, 240, 672, 1792, (unsigned)-1
+/* 7, 25, 81, 241, 673, 1793, (unsigned)-1 ### The old approximations */
+};
+
+/*
+ * Perform modular exponentiation, as fast as possible! This uses
+ * Montgomery reduction, optimized squaring, and windowed exponentiation.
+ * The modulus "mod" MUST be odd!
+ *
+ * This returns 0 on success, -1 on out of memory.
+ *
+ * The window algorithm:
+ * The idea is to keep a running product of b1 = n^(high-order bits of exp),
+ * and then keep appending exponent bits to it. The following patterns
+ * apply to a 3-bit window (k = 3):
+ * To append 0: square
+ * To append 1: square, multiply by n^1
+ * To append 10: square, multiply by n^1, square
+ * To append 11: square, square, multiply by n^3
+ * To append 100: square, multiply by n^1, square, square
+ * To append 101: square, square, square, multiply by n^5
+ * To append 110: square, square, multiply by n^3, square
+ * To append 111: square, square, square, multiply by n^7
+ *
+ * Since each pattern involves only one multiply, the longer the pattern
+ * the better, except that a 0 (no multiplies) can be appended directly.
+ * We precompute a table of odd powers of n, up to 2^k, and can then
+ * multiply k bits of exponent at a time. Actually, assuming random
+ * exponents, there is on average one zero bit between needs to
+ * multiply (1/2 of the time there's none, 1/4 of the time there's 1,
+ * 1/8 of the time, there's 2, 1/32 of the time, there's 3, etc.), so
+ * you have to do one multiply per k+1 bits of exponent.
+ *
+ * The loop walks down the exponent, squaring the result buffer as
+ * it goes. There is a wbits+1 bit lookahead buffer, buf, that is
+ * filled with the upcoming exponent bits. (What is read after the
+ * end of the exponent is unimportant, but it is filled with zero here.)
+ * When the most-significant bit of this buffer becomes set, i.e.
+ * (buf & tblmask) != 0, we have to decide what pattern to multiply
+ * by, and when to do it. We decide, remember to do it in future
+ * after a suitable number of squarings have passed (e.g. a pattern
+ * of "100" in the buffer requires that we multiply by n^1 immediately;
+ * a pattern of "110" calls for multiplying by n^3 after one more
+ * squaring), clear the buffer, and continue.
+ *
+ * When we start, there is one more optimization: the result buffer
+ * is implcitly one, so squaring it or multiplying by it can be
+ * optimized away. Further, if we start with a pattern like "100"
+ * in the lookahead window, rather than placing n into the buffer
+ * and then starting to square it, we have already computed n^2
+ * to compute the odd-powers table, so we can place that into
+ * the buffer and save a squaring.
+ *
+ * This means that if you have a k-bit window, to compute n^z,
+ * where z is the high k bits of the exponent, 1/2 of the time
+ * it requires no squarings. 1/4 of the time, it requires 1
+ * squaring, ... 1/2^(k-1) of the time, it reqires k-2 squarings.
+ * And the remaining 1/2^(k-1) of the time, the top k bits are a
+ * 1 followed by k-1 0 bits, so it again only requires k-2
+ * squarings, not k-1. The average of these is 1. Add that
+ * to the one squaring we have to do to compute the table,
+ * and you'll see that a k-bit window saves k-2 squarings
+ * as well as reducing the multiplies. (It actually doesn't
+ * hurt in the case k = 1, either.)
+ *
+ * n must have mlen words allocated. Although fewer may be in use
+ * when n is passed in, all are in use on exit.
+ */
+int
+lbnExpMod_32(BNWORD32 *result, BNWORD32 const *n, unsigned nlen,
+ BNWORD32 const *e, unsigned elen, BNWORD32 *mod, unsigned mlen)
+{
+ BNWORD32 *table[1 << (BNEXPMOD_MAX_WINDOW-1)];
+ /* Table of odd powers of n */
+ unsigned ebits; /* Exponent bits */
+ unsigned wbits; /* Window size */
+ unsigned tblmask; /* Mask of exponentiation window */
+ BNWORD32 bitpos; /* Mask of current look-ahead bit */
+ unsigned buf; /* Buffer of exponent bits */
+ unsigned multpos; /* Where to do pending multiply */
+ BNWORD32 const *mult; /* What to multiply by */
+ unsigned i; /* Loop counter */
+ int isone; /* Flag: accum. is implicitly one */
+ BNWORD32 *a, *b; /* Working buffers/accumulators */
+ BNWORD32 *t; /* Pointer into the working buffers */
+ BNWORD32 inv; /* mod^-1 modulo 2^32 */
+ int y; /* bnYield() result */
+
+ assert(mlen);
+ assert(nlen <= mlen);
+
+ /* First, a couple of trivial cases. */
+ elen = lbnNorm_32(e, elen);
+ if (!elen) {
+ /* x ^ 0 == 1 */
+ lbnZero_32(result, mlen);
+ BIGLITTLE(result[-1],result[0]) = 1;
+ return 0;
+ }
+ ebits = lbnBits_32(e, elen);
+ if (ebits == 1) {
+ /* x ^ 1 == x */
+ if (n != result)
+ lbnCopy_32(result, n, nlen);
+ if (mlen > nlen)
+ lbnZero_32(BIGLITTLE(result-nlen,result+nlen),
+ mlen-nlen);
+ return 0;
+ }
+
+ /* Okay, now move the exponent pointer to the most-significant word */
+ e = BIGLITTLE(e-elen, e+elen-1);
+
+ /* Look up appropriate k-1 for the exponent - tblmask = 1<<(k-1) */
+ wbits = 0;
+ while (ebits > bnExpModThreshTable[wbits])
+ wbits++;
+
+ /* Allocate working storage: two product buffers and the tables. */
+ LBNALLOC(a, BNWORD32, 2*mlen);
+ if (!a)
+ return -1;
+ LBNALLOC(b, BNWORD32, 2*mlen);
+ if (!b) {
+ LBNFREE(a, 2*mlen);
+ return -1;
+ }
+
+ /* Convert to the appropriate table size: tblmask = 1<<(k-1) */
+ tblmask = 1u << wbits;
+
+ /* We have the result buffer available, so use it. */
+ table[0] = result;
+
+ /*
+ * Okay, we now have a minimal-sized table - expand it.
+ * This is allowed to fail! If so, scale back the table size
+ * and proceed.
+ */
+ for (i = 1; i < tblmask; i++) {
+ LBNALLOC(t, BNWORD32, mlen);
+ if (!t) /* Out of memory! Quit the loop. */
+ break;
+ table[i] = t;
+ }
+
+ /* If we stopped, with i < tblmask, shrink the tables appropriately */
+ while (tblmask > i) {
+ wbits--;
+ tblmask >>= 1;
+ }
+ /* Free up our overallocations */
+ while (--i > tblmask)
+ LBNFREE(table[i], mlen);
+
+ /* Okay, fill in the table */
+
+ /* Compute the necessary modular inverse */
+ inv = lbnMontInv1_32(mod[BIGLITTLE(-1,0)]); /* LSW of modulus */
+
+ /* Convert n to Montgomery form */
+
+ /* Move n up "mlen" words into a */
+ t = BIGLITTLE(a-mlen, a+mlen);
+ lbnCopy_32(t, n, nlen);
+ lbnZero_32(a, mlen);
+ /* Do the division - lose the quotient into the high-order words */
+ (void)lbnDiv_32(t, a, mlen+nlen, mod, mlen);
+ /* Copy into first table entry */
+ lbnCopy_32(table[0], a, mlen);
+
+ /* Square a into b */
+ lbnMontSquare_32(b, a, mod, mlen, inv);
+
+ /* Use high half of b to initialize the table */
+ t = BIGLITTLE(b-mlen, b+mlen);
+ for (i = 1; i < tblmask; i++) {
+ lbnMontMul_32(a, t, table[i-1], mod, mlen, inv);
+ lbnCopy_32(table[i], BIGLITTLE(a-mlen, a+mlen), mlen);
+#if BNYIELD
+ if (bnYield && (y = bnYield()) < 0)
+ goto yield;
+#endif
+ }
+
+ /* We might use b = n^2 later... */
+
+ /* Initialze the fetch pointer */
+ bitpos = (BNWORD32)1 << ((ebits-1) & (32-1)); /* Initialize mask */
+
+ /* This should point to the msbit of e */
+ assert((*e & bitpos) != 0);
+
+ /*
+ * Pre-load the window. Becuase the window size is
+ * never larger than the exponent size, there is no need to
+ * detect running off the end of e in here.
+ *
+ * The read-ahead is controlled by elen and the bitpos mask.
+ * Note that this is *ahead* of ebits, which tracks the
+ * most significant end of the window. The purpose of this
+ * initialization is to get the two wbits+1 bits apart,
+ * like they should be.
+ *
+ * Note that bitpos and e1len together keep track of the
+ * lookahead read pointer in the exponent that is used here.
+ */
+ buf = 0;
+ for (i = 0; i <= wbits; i++) {
+ buf = (buf << 1) | ((*e & bitpos) != 0);
+ bitpos >>= 1;
+ if (!bitpos) {
+ BIGLITTLE(e++,e--);
+ bitpos = (BNWORD32)1 << (32-1);
+ elen--;
+ }
+ }
+ assert(buf & tblmask);
+
+ /*
+ * Set the pending multiply positions to a location that will
+ * never be encountered, thus ensuring that nothing will happen
+ * until the need for a multiply appears and one is scheduled.
+ */
+ multpos = ebits; /* A NULL value */
+ mult = 0; /* Force a crash if we use these */
+
+ /*
+ * Okay, now begins the real work. The first step is
+ * slightly magic, so it's done outside the main loop,
+ * but it's very similar to what's inside.
+ */
+ ebits--; /* Start processing the first bit... */
+ isone = 1;
+
+ /*
+ * This is just like the multiply in the loop, except that
+ * - We know the msbit of buf is set, and
+ * - We have the extra value n^2 floating around.
+ * So, do the usual computation, and if the result is that
+ * the buffer should be multiplied by n^1 immediately
+ * (which we'd normally then square), we multiply it
+ * (which reduces to a copy, which reduces to setting a flag)
+ * by n^2 and skip the squaring. Thus, we do the
+ * multiply and the squaring in one step.
+ */
+ assert(buf & tblmask);
+ multpos = ebits - wbits;
+ while ((buf & 1) == 0) {
+ buf >>= 1;
+ multpos++;
+ }
+ /* Intermediates can wrap, but final must NOT */
+ assert(multpos <= ebits);
+ mult = table[buf>>1];
+ buf = 0;
+
+ /* Special case: use already-computed value sitting in buffer */
+ if (multpos == ebits)
+ isone = 0;
+
+ /*
+ * At this point, the buffer (which is the high half of b) holds
+ * either 1 (implicitly, as the "isone" flag is set), or n^2.
+ */
+
+ /*
+ * The main loop. The procedure is:
+ * - Advance the window
+ * - If the most-significant bit of the window is set,
+ * schedule a multiply for the appropriate time in the
+ * future (may be immediately)
+ * - Perform any pending multiples
+ * - Check for termination
+ * - Square the buffer
+ *
+ * At any given time, the acumulated product is held in
+ * the high half of b.
+ */
+ for (;;) {
+ ebits--;
+
+ /* Advance the window */
+ assert(buf < tblmask);
+ buf <<= 1;
+ /*
+ * This reads ahead of the current exponent position
+ * (controlled by ebits), so we have to be able to read
+ * past the lsb of the exponents without error.
+ */
+ if (elen) {
+ buf |= ((*e & bitpos) != 0);
+ bitpos >>= 1;
+ if (!bitpos) {
+ BIGLITTLE(e++,e--);
+ bitpos = (BNWORD32)1 << (32-1);
+ elen--;
+ }
+ }
+
+ /* Examine the window for pending multiplies */
+ if (buf & tblmask) {
+ multpos = ebits - wbits;
+ while ((buf & 1) == 0) {
+ buf >>= 1;
+ multpos++;
+ }
+ /* Intermediates can wrap, but final must NOT */
+ assert(multpos <= ebits);
+ mult = table[buf>>1];
+ buf = 0;
+ }
+
+ /* If we have a pending multiply, do it */
+ if (ebits == multpos) {
+ /* Multiply by the table entry remembered previously */
+ t = BIGLITTLE(b-mlen, b+mlen);
+ if (isone) {
+ /* Multiply by 1 is a trivial case */
+ lbnCopy_32(t, mult, mlen);
+ isone = 0;
+ } else {
+ lbnMontMul_32(a, t, mult, mod, mlen, inv);
+ /* Swap a and b */
+ t = a; a = b; b = t;
+ }
+ }
+
+ /* Are we done? */
+ if (!ebits)
+ break;
+
+ /* Square the input */
+ if (!isone) {
+ t = BIGLITTLE(b-mlen, b+mlen);
+ lbnMontSquare_32(a, t, mod, mlen, inv);
+ /* Swap a and b */
+ t = a; a = b; b = t;
+ }
+#if BNYIELD
+ if (bnYield && (y = bnYield()) < 0)
+ goto yield;
+#endif
+ } /* for (;;) */
+
+ assert(!isone);
+ assert(!buf);
+
+ /* DONE! */
+
+ /* Convert result out of Montgomery form */
+ t = BIGLITTLE(b-mlen, b+mlen);
+ lbnCopy_32(b, t, mlen);
+ lbnZero_32(t, mlen);
+ lbnMontReduce_32(b, mod, mlen, inv);
+ lbnCopy_32(result, t, mlen);
+ /*
+ * Clean up - free intermediate storage.
+ * Do NOT free table[0], which is the result
+ * buffer.
+ */
+ y = 0;
+#if BNYIELD
+yield:
+#endif
+ while (--tblmask)
+ LBNFREE(table[tblmask], mlen);
+ LBNFREE(b, 2*mlen);
+ LBNFREE(a, 2*mlen);
+
+ return y; /* Success */
+}
+
+/*
+ * Compute and return n1^e1 * n2^e2 mod "mod".
+ * result may be either input buffer, or something separate.
+ * It must be "mlen" words long.
+ *
+ * There is a current position in the exponents, which is kept in e1bits.
+ * (The exponents are swapped if necessary so e1 is the longer of the two.)
+ * At any given time, the value in the accumulator is
+ * n1^(e1>>e1bits) * n2^(e2>>e1bits) mod "mod".
+ * As e1bits is counted down, this is updated, by squaring it and doing
+ * any necessary multiplies.
+ * To decide on the necessary multiplies, two windows, each w1bits+1 bits
+ * wide, are maintained in buf1 and buf2, which read *ahead* of the
+ * e1bits position (with appropriate handling of the case when e1bits
+ * drops below w1bits+1). When the most-significant bit of either window
+ * becomes set, indicating that something needs to be multiplied by
+ * the accumulator or it will get out of sync, the window is examined
+ * to see which power of n1 or n2 to multiply by, and when (possibly
+ * later, if the power is greater than 1) the multiply should take
+ * place. Then the multiply and its location are remembered and the
+ * window is cleared.
+ *
+ * If we had every power of n1 in the table, the multiply would always
+ * be w1bits steps in the future. But we only keep the odd powers,
+ * so instead of waiting w1bits squarings and then multiplying
+ * by n1^k, we wait w1bits-k squarings and multiply by n1.
+ *
+ * Actually, w2bits can be less than w1bits, but the window is the same
+ * size, to make it easier to keep track of where we're reading. The
+ * appropriate number of low-order bits of the window are just ignored.
+ */
+int
+lbnDoubleExpMod_32(BNWORD32 *result,
+ BNWORD32 const *n1, unsigned n1len,
+ BNWORD32 const *e1, unsigned e1len,
+ BNWORD32 const *n2, unsigned n2len,
+ BNWORD32 const *e2, unsigned e2len,
+ BNWORD32 *mod, unsigned mlen)
+{
+ BNWORD32 *table1[1 << (BNEXPMOD_MAX_WINDOW-1)];
+ /* Table of odd powers of n1 */
+ BNWORD32 *table2[1 << (BNEXPMOD_MAX_WINDOW-1)];
+ /* Table of odd powers of n2 */
+ unsigned e1bits, e2bits; /* Exponent bits */
+ unsigned w1bits, w2bits; /* Window sizes */
+ unsigned tblmask; /* Mask of exponentiation window */
+ BNWORD32 bitpos; /* Mask of current look-ahead bit */
+ unsigned buf1, buf2; /* Buffer of exponent bits */
+ unsigned mult1pos, mult2pos; /* Where to do pending multiply */
+ BNWORD32 const *mult1, *mult2; /* What to multiply by */
+ unsigned i; /* Loop counter */
+ int isone; /* Flag: accum. is implicitly one */
+ BNWORD32 *a, *b; /* Working buffers/accumulators */
+ BNWORD32 *t; /* Pointer into the working buffers */
+ BNWORD32 inv; /* mod^-1 modulo 2^32 */
+ int y; /* bnYield() result */
+
+ assert(mlen);
+ assert(n1len <= mlen);
+ assert(n2len <= mlen);
+
+ /* First, a couple of trivial cases. */
+ e1len = lbnNorm_32(e1, e1len);
+ e2len = lbnNorm_32(e2, e2len);
+
+ /* Ensure that the first exponent is the longer */
+ e1bits = lbnBits_32(e1, e1len);
+ e2bits = lbnBits_32(e2, e2len);
+ if (e1bits < e2bits) {
+ i = e1len; e1len = e2len; e2len = i;
+ i = e1bits; e1bits = e2bits; e2bits = i;
+ t = (BNWORD32 *)n1; n1 = n2; n2 = t;
+ t = (BNWORD32 *)e1; e1 = e2; e2 = t;
+ }
+ assert(e1bits >= e2bits);
+
+ /* Handle a trivial case */
+ if (!e2len)
+ return lbnExpMod_32(result, n1, n1len, e1, e1len, mod, mlen);
+ assert(e2bits);
+
+ /* The code below fucks up if the exponents aren't at least 2 bits */
+ if (e1bits == 1) {
+ assert(e2bits == 1);
+
+ LBNALLOC(a, BNWORD32, n1len+n2len);
+ if (!a)
+ return -1;
+
+ lbnMul_32(a, n1, n1len, n2, n2len);
+ /* Do a direct modular reduction */
+ if (n1len + n2len >= mlen)
+ (void)lbnDiv_32(a+mlen, a, n1len+n2len, mod, mlen);
+ lbnCopy_32(result, a, mlen);
+ LBNFREE(a, n1len+n2len);
+ return 0;
+ }
+
+ /* Okay, now move the exponent pointers to the most-significant word */
+ e1 = BIGLITTLE(e1-e1len, e1+e1len-1);
+ e2 = BIGLITTLE(e2-e2len, e2+e2len-1);
+
+ /* Look up appropriate k-1 for the exponent - tblmask = 1<<(k-1) */
+ w1bits = 0;
+ while (e1bits > bnExpModThreshTable[w1bits])
+ w1bits++;
+ w2bits = 0;
+ while (e2bits > bnExpModThreshTable[w2bits])
+ w2bits++;
+
+ assert(w1bits >= w2bits);
+
+ /* Allocate working storage: two product buffers and the tables. */
+ LBNALLOC(a, BNWORD32, 2*mlen);
+ if (!a)
+ return -1;
+ LBNALLOC(b, BNWORD32, 2*mlen);
+ if (!b) {
+ LBNFREE(a, 2*mlen);
+ return -1;
+ }
+
+ /* Convert to the appropriate table size: tblmask = 1<<(k-1) */
+ tblmask = 1u << w1bits;
+ /* Use buf2 for its size, temporarily */
+ buf2 = 1u << w2bits;
+
+ LBNALLOC(t, BNWORD32, mlen);
+ if (!t) {
+ LBNFREE(b, 2*mlen);
+ LBNFREE(a, 2*mlen);
+ return -1;
+ }
+ table1[0] = t;
+ table2[0] = result;
+
+ /*
+ * Okay, we now have some minimal-sized tables - expand them.
+ * This is allowed to fail! If so, scale back the table sizes
+ * and proceed. We allocate both tables at the same time
+ * so if it fails partway through, they'll both be a reasonable
+ * size rather than one huge and one tiny.
+ * When i passes buf2 (the number of entries in the e2 window,
+ * which may be less than the number of entries in the e1 window),
+ * stop allocating e2 space.
+ */
+ for (i = 1; i < tblmask; i++) {
+ LBNALLOC(t, BNWORD32, mlen);
+ if (!t) /* Out of memory! Quit the loop. */
+ break;
+ table1[i] = t;
+ if (i < buf2) {
+ LBNALLOC(t, BNWORD32, mlen);
+ if (!t) {
+ LBNFREE(table1[i], mlen);
+ break;
+ }
+ table2[i] = t;
+ }
+ }
+
+ /* If we stopped, with i < tblmask, shrink the tables appropriately */
+ while (tblmask > i) {
+ w1bits--;
+ tblmask >>= 1;
+ }
+ /* Free up our overallocations */
+ while (--i > tblmask) {
+ if (i < buf2)
+ LBNFREE(table2[i], mlen);
+ LBNFREE(table1[i], mlen);
+ }
+ /* And shrink the second window too, if needed */
+ if (w2bits > w1bits) {
+ w2bits = w1bits;
+ buf2 = tblmask;
+ }
+
+ /*
+ * From now on, use the w2bits variable for the difference
+ * between w1bits and w2bits.
+ */
+ w2bits = w1bits-w2bits;
+
+ /* Okay, fill in the tables */
+
+ /* Compute the necessary modular inverse */
+ inv = lbnMontInv1_32(mod[BIGLITTLE(-1,0)]); /* LSW of modulus */
+
+ /* Convert n1 to Montgomery form */
+
+ /* Move n1 up "mlen" words into a */
+ t = BIGLITTLE(a-mlen, a+mlen);
+ lbnCopy_32(t, n1, n1len);
+ lbnZero_32(a, mlen);
+ /* Do the division - lose the quotient into the high-order words */
+ (void)lbnDiv_32(t, a, mlen+n1len, mod, mlen);
+ /* Copy into first table entry */
+ lbnCopy_32(table1[0], a, mlen);
+
+ /* Square a into b */
+ lbnMontSquare_32(b, a, mod, mlen, inv);
+
+ /* Use high half of b to initialize the first table */
+ t = BIGLITTLE(b-mlen, b+mlen);
+ for (i = 1; i < tblmask; i++) {
+ lbnMontMul_32(a, t, table1[i-1], mod, mlen, inv);
+ lbnCopy_32(table1[i], BIGLITTLE(a-mlen, a+mlen), mlen);
+#if BNYIELD
+ if (bnYield && (y = bnYield()) < 0)
+ goto yield;
+#endif
+ }
+
+ /* Convert n2 to Montgomery form */
+
+ t = BIGLITTLE(a-mlen, a+mlen);
+ /* Move n2 up "mlen" words into a */
+ lbnCopy_32(t, n2, n2len);
+ lbnZero_32(a, mlen);
+ /* Do the division - lose the quotient into the high-order words */
+ (void)lbnDiv_32(t, a, mlen+n2len, mod, mlen);
+ /* Copy into first table entry */
+ lbnCopy_32(table2[0], a, mlen);
+
+ /* Square it into a */
+ lbnMontSquare_32(a, table2[0], mod, mlen, inv);
+ /* Copy to b, low half */
+ lbnCopy_32(b, t, mlen);
+
+ /* Use b to initialize the second table */
+ for (i = 1; i < buf2; i++) {
+ lbnMontMul_32(a, b, table2[i-1], mod, mlen, inv);
+ lbnCopy_32(table2[i], t, mlen);
+#if BNYIELD
+ if (bnYield && (y = bnYield()) < 0)
+ goto yield;
+#endif
+ }
+
+ /*
+ * Okay, a recap: at this point, the low part of b holds
+ * n2^2, the high part holds n1^2, and the tables are
+ * initialized with the odd powers of n1 and n2 from 1
+ * through 2*tblmask-1 and 2*buf2-1.
+ *
+ * We might use those squares in b later, or we might not.
+ */
+
+ /* Initialze the fetch pointer */
+ bitpos = (BNWORD32)1 << ((e1bits-1) & (32-1)); /* Initialize mask */
+
+ /* This should point to the msbit of e1 */
+ assert((*e1 & bitpos) != 0);
+
+ /*
+ * Pre-load the windows. Becuase the window size is
+ * never larger than the exponent size, there is no need to
+ * detect running off the end of e1 in here.
+ *
+ * The read-ahead is controlled by e1len and the bitpos mask.
+ * Note that this is *ahead* of e1bits, which tracks the
+ * most significant end of the window. The purpose of this
+ * initialization is to get the two w1bits+1 bits apart,
+ * like they should be.
+ *
+ * Note that bitpos and e1len together keep track of the
+ * lookahead read pointer in the exponent that is used here.
+ * e2len is not decremented, it is only ever compared with
+ * e1len as *that* is decremented.
+ */
+ buf1 = buf2 = 0;
+ for (i = 0; i <= w1bits; i++) {
+ buf1 = (buf1 << 1) | ((*e1 & bitpos) != 0);
+ if (e1len <= e2len)
+ buf2 = (buf2 << 1) | ((*e2 & bitpos) != 0);
+ bitpos >>= 1;
+ if (!bitpos) {
+ BIGLITTLE(e1++,e1--);
+ if (e1len <= e2len)
+ BIGLITTLE(e2++,e2--);
+ bitpos = (BNWORD32)1 << (32-1);
+ e1len--;
+ }
+ }
+ assert(buf1 & tblmask);
+
+ /*
+ * Set the pending multiply positions to a location that will
+ * never be encountered, thus ensuring that nothing will happen
+ * until the need for a multiply appears and one is scheduled.
+ */
+ mult1pos = mult2pos = e1bits; /* A NULL value */
+ mult1 = mult2 = 0; /* Force a crash if we use these */
+
+ /*
+ * Okay, now begins the real work. The first step is
+ * slightly magic, so it's done outside the main loop,
+ * but it's very similar to what's inside.
+ */
+ isone = 1; /* Buffer is implicitly 1, so replace * by copy */
+ e1bits--; /* Start processing the first bit... */
+
+ /*
+ * This is just like the multiply in the loop, except that
+ * - We know the msbit of buf1 is set, and
+ * - We have the extra value n1^2 floating around.
+ * So, do the usual computation, and if the result is that
+ * the buffer should be multiplied by n1^1 immediately
+ * (which we'd normally then square), we multiply it
+ * (which reduces to a copy, which reduces to setting a flag)
+ * by n1^2 and skip the squaring. Thus, we do the
+ * multiply and the squaring in one step.
+ */
+ assert(buf1 & tblmask);
+ mult1pos = e1bits - w1bits;
+ while ((buf1 & 1) == 0) {
+ buf1 >>= 1;
+ mult1pos++;
+ }
+ /* Intermediates can wrap, but final must NOT */
+ assert(mult1pos <= e1bits);
+ mult1 = table1[buf1>>1];
+ buf1 = 0;
+
+ /* Special case: use already-computed value sitting in buffer */
+ if (mult1pos == e1bits)
+ isone = 0;
+
+ /*
+ * The first multiply by a power of n2. Similar, but
+ * we might not even want to schedule a multiply if e2 is
+ * shorter than e1, and the window might be shorter so
+ * we have to leave the low w2bits bits alone.
+ */
+ if (buf2 & tblmask) {
+ /* Remember low-order bits for later */
+ i = buf2 & ((1u << w2bits) - 1);
+ buf2 >>= w2bits;
+ mult2pos = e1bits - w1bits + w2bits;
+ while ((buf2 & 1) == 0) {
+ buf2 >>= 1;
+ mult2pos++;
+ }
+ assert(mult2pos <= e1bits);
+ mult2 = table2[buf2>>1];
+ buf2 = i;
+
+ if (mult2pos == e1bits) {
+ t = BIGLITTLE(b-mlen, b+mlen);
+ if (isone) {
+ lbnCopy_32(t, b, mlen); /* Copy low to high */
+ isone = 0;
+ } else {
+ lbnMontMul_32(a, t, b, mod, mlen, inv);
+ t = a; a = b; b = t;
+ }
+ }
+ }
+
+ /*
+ * At this point, the buffer (which is the high half of b)
+ * holds either 1 (implicitly, as the "isone" flag is set),
+ * n1^2, n2^2 or n1^2 * n2^2.
+ */
+
+ /*
+ * The main loop. The procedure is:
+ * - Advance the windows
+ * - If the most-significant bit of a window is set,
+ * schedule a multiply for the appropriate time in the
+ * future (may be immediately)
+ * - Perform any pending multiples
+ * - Check for termination
+ * - Square the buffers
+ *
+ * At any given time, the acumulated product is held in
+ * the high half of b.
+ */
+ for (;;) {
+ e1bits--;
+
+ /* Advance the windows */
+ assert(buf1 < tblmask);
+ buf1 <<= 1;
+ assert(buf2 < tblmask);
+ buf2 <<= 1;
+ /*
+ * This reads ahead of the current exponent position
+ * (controlled by e1bits), so we have to be able to read
+ * past the lsb of the exponents without error.
+ */
+ if (e1len) {
+ buf1 |= ((*e1 & bitpos) != 0);
+ if (e1len <= e2len)
+ buf2 |= ((*e2 & bitpos) != 0);
+ bitpos >>= 1;
+ if (!bitpos) {
+ BIGLITTLE(e1++,e1--);
+ if (e1len <= e2len)
+ BIGLITTLE(e2++,e2--);
+ bitpos = (BNWORD32)1 << (32-1);
+ e1len--;
+ }
+ }
+
+ /* Examine the first window for pending multiplies */
+ if (buf1 & tblmask) {
+ mult1pos = e1bits - w1bits;
+ while ((buf1 & 1) == 0) {
+ buf1 >>= 1;
+ mult1pos++;
+ }
+ /* Intermediates can wrap, but final must NOT */
+ assert(mult1pos <= e1bits);
+ mult1 = table1[buf1>>1];
+ buf1 = 0;
+ }
+
+ /*
+ * Examine the second window for pending multiplies.
+ * Window 2 can be smaller than window 1, but we
+ * keep the same number of bits in buf2, so we need
+ * to ignore any low-order bits in the buffer when
+ * computing what to multiply by, and recompute them
+ * later.
+ */
+ if (buf2 & tblmask) {
+ /* Remember low-order bits for later */
+ i = buf2 & ((1u << w2bits) - 1);
+ buf2 >>= w2bits;
+ mult2pos = e1bits - w1bits + w2bits;
+ while ((buf2 & 1) == 0) {
+ buf2 >>= 1;
+ mult2pos++;
+ }
+ assert(mult2pos <= e1bits);
+ mult2 = table2[buf2>>1];
+ buf2 = i;
+ }
+
+
+ /* If we have a pending multiply for e1, do it */
+ if (e1bits == mult1pos) {
+ /* Multiply by the table entry remembered previously */
+ t = BIGLITTLE(b-mlen, b+mlen);
+ if (isone) {
+ /* Multiply by 1 is a trivial case */
+ lbnCopy_32(t, mult1, mlen);
+ isone = 0;
+ } else {
+ lbnMontMul_32(a, t, mult1, mod, mlen, inv);
+ /* Swap a and b */
+ t = a; a = b; b = t;
+ }
+ }
+
+ /* If we have a pending multiply for e2, do it */
+ if (e1bits == mult2pos) {
+ /* Multiply by the table entry remembered previously */
+ t = BIGLITTLE(b-mlen, b+mlen);
+ if (isone) {
+ /* Multiply by 1 is a trivial case */
+ lbnCopy_32(t, mult2, mlen);
+ isone = 0;
+ } else {
+ lbnMontMul_32(a, t, mult2, mod, mlen, inv);
+ /* Swap a and b */
+ t = a; a = b; b = t;
+ }
+ }
+
+ /* Are we done? */
+ if (!e1bits)
+ break;
+
+ /* Square the buffer */
+ if (!isone) {
+ t = BIGLITTLE(b-mlen, b+mlen);
+ lbnMontSquare_32(a, t, mod, mlen, inv);
+ /* Swap a and b */
+ t = a; a = b; b = t;
+ }
+#if BNYIELD
+ if (bnYield && (y = bnYield()) < 0)
+ goto yield;
+#endif
+ } /* for (;;) */
+
+ assert(!isone);
+ assert(!buf1);
+ assert(!buf2);
+
+ /* DONE! */
+
+ /* Convert result out of Montgomery form */
+ t = BIGLITTLE(b-mlen, b+mlen);
+ lbnCopy_32(b, t, mlen);
+ lbnZero_32(t, mlen);
+ lbnMontReduce_32(b, mod, mlen, inv);
+ lbnCopy_32(result, t, mlen);
+
+ /* Clean up - free intermediate storage */
+ y = 0;
+#if BNYIELD
+yield:
+#endif
+ buf2 = tblmask >> w2bits;
+ while (--tblmask) {
+ if (tblmask < buf2)
+ LBNFREE(table2[tblmask], mlen);
+ LBNFREE(table1[tblmask], mlen);
+ }
+ t = table1[0];
+ LBNFREE(t, mlen);
+ LBNFREE(b, 2*mlen);
+ LBNFREE(a, 2*mlen);
+
+ return y; /* Success */
+}
+
+/*
+ * 2^exp (mod mod). This is an optimized version for use in Fermat
+ * tests. The input value of n is ignored; it is returned with
+ * "mlen" words valid.
+ */
+int
+lbnTwoExpMod_32(BNWORD32 *n, BNWORD32 const *exp, unsigned elen,
+ BNWORD32 *mod, unsigned mlen)
+{
+ unsigned e; /* Copy of high words of the exponent */
+ unsigned bits; /* Assorted counter of bits */
+ BNWORD32 const *bitptr;
+ BNWORD32 bitword, bitpos;
+ BNWORD32 *a, *b, *a1;
+ BNWORD32 inv;
+ int y; /* Result of bnYield() */
+
+ assert(mlen);
+
+ bitptr = BIGLITTLE(exp-elen, exp+elen-1);
+ bitword = *bitptr;
+ assert(bitword);
+
+ /* Clear n for future use. */
+ lbnZero_32(n, mlen);
+
+ bits = lbnBits_32(exp, elen);
+
+ /* First, a couple of trivial cases. */
+ if (bits <= 1) {
+ /* 2 ^ 0 == 1, 2 ^ 1 == 2 */
+ BIGLITTLE(n[-1],n[0]) = (BNWORD32)1<<elen;
+ return 0;
+ }
+
+ /* Set bitpos to the most significant bit */
+ bitpos = (BNWORD32)1 << ((bits-1) & (32-1));
+
+ /* Now, count the bits in the modulus. */
+ bits = lbnBits_32(mod, mlen);
+ assert(bits > 1); /* a 1-bit modulus is just stupid... */
+
+ /*
+ * We start with 1<<e, where "e" is as many high bits of the
+ * exponent as we can manage without going over the modulus.
+ * This first loop finds "e".
+ */
+ e = 1;
+ while (elen) {
+ /* Consume the first bit */
+ bitpos >>= 1;
+ if (!bitpos) {
+ if (!--elen)
+ break;
+ bitword = BIGLITTLE(*++bitptr,*--bitptr);
+ bitpos = (BNWORD32)1<<(32-1);
+ }
+ e = (e << 1) | ((bitpos & bitword) != 0);
+ if (e >= bits) { /* Overflow! Back out. */
+ e >>= 1;
+ break;
+ }
+ }
+ /*
+ * The bit in "bitpos" being examined by the bit buffer has NOT
+ * been consumed yet. This may be past the end of the exponent,
+ * in which case elen == 1.
+ */
+
+ /* Okay, now, set bit "e" in n. n is already zero. */
+ inv = (BNWORD32)1 << (e & (32-1));
+ e /= 32;
+ BIGLITTLE(n[-e-1],n[e]) = inv;
+ /*
+ * The effective length of n in words is now "e+1".
+ * This is used a little bit later.
+ */
+
+ if (!elen)
+ return 0; /* That was easy! */
+
+ /*
+ * We have now processed the first few bits. The next step
+ * is to convert this to Montgomery form for further squaring.
+ */
+
+ /* Allocate working storage: two product buffers */
+ LBNALLOC(a, BNWORD32, 2*mlen);
+ if (!a)
+ return -1;
+ LBNALLOC(b, BNWORD32, 2*mlen);
+ if (!b) {
+ LBNFREE(a, 2*mlen);
+ return -1;
+ }
+
+ /* Convert n to Montgomery form */
+ inv = BIGLITTLE(mod[-1],mod[0]); /* LSW of modulus */
+ assert(inv & 1); /* Modulus must be odd */
+ inv = lbnMontInv1_32(inv);
+ /* Move n (length e+1, remember?) up "mlen" words into b */
+ /* Note that we lie about a1 for a bit - it's pointing to b */
+ a1 = BIGLITTLE(b-mlen,b+mlen);
+ lbnCopy_32(a1, n, e+1);
+ lbnZero_32(b, mlen);
+ /* Do the division - dump the quotient into the high-order words */
+ (void)lbnDiv_32(a1, b, mlen+e+1, mod, mlen);
+ /*
+ * Now do the first squaring and modular reduction to put
+ * the number up in a1 where it belongs.
+ */
+ lbnMontSquare_32(a, b, mod, mlen, inv);
+ /* Fix up a1 to point to where it should go. */
+ a1 = BIGLITTLE(a-mlen,a+mlen);
+
+ /*
+ * Okay, now, a1 holds the number being accumulated, and
+ * b is a scratch register. Start working:
+ */
+ for (;;) {
+ /*
+ * Is the bit set? If so, double a1 as well.
+ * A modular doubling like this is very cheap.
+ */
+ if (bitpos & bitword) {
+ /*
+ * Double the number. If there was a carry out OR
+ * the result is greater than the modulus, subract
+ * the modulus.
+ */
+ if (lbnDouble_32(a1, mlen) ||
+ lbnCmp_32(a1, mod, mlen) > 0)
+ (void)lbnSubN_32(a1, mod, mlen);
+ }
+
+ /* Advance to the next exponent bit */
+ bitpos >>= 1;
+ if (!bitpos) {
+ if (!--elen)
+ break; /* Done! */
+ bitword = BIGLITTLE(*++bitptr,*--bitptr);
+ bitpos = (BNWORD32)1<<(32-1);
+ }
+
+ /*
+ * The elen/bitword/bitpos bit buffer is known to be
+ * non-empty, i.e. there is at least one more unconsumed bit.
+ * Thus, it's safe to square the number.
+ */
+ lbnMontSquare_32(b, a1, mod, mlen, inv);
+ /* Rename result (in b) back to a (a1, really). */
+ a1 = b; b = a; a = a1;
+ a1 = BIGLITTLE(a-mlen,a+mlen);
+#if BNYIELD
+ if (bnYield && (y = bnYield()) < 0)
+ goto yield;
+#endif
+ }
+
+ /* DONE! Just a little bit of cleanup... */
+
+ /*
+ * Convert result out of Montgomery form... this is
+ * just a Montgomery reduction.
+ */
+ lbnCopy_32(a, a1, mlen);
+ lbnZero_32(a1, mlen);
+ lbnMontReduce_32(a, mod, mlen, inv);
+ lbnCopy_32(n, a1, mlen);
+
+ /* Clean up - free intermediate storage */
+ y = 0;
+#if BNYIELD
+yield:
+#endif
+ LBNFREE(b, 2*mlen);
+ LBNFREE(a, 2*mlen);
+
+ return y; /* Success */
+}
+
+
+/*
+ * Returns a substring of the big-endian array of bytes representation
+ * of the bignum array based on two parameters, the least significant
+ * byte number (0 to start with the least significant byte) and the
+ * length. I.e. the number returned is a representation of
+ * (bn / 2^(8*lsbyte)) % 2 ^ (8*buflen).
+ *
+ * It is an error if the bignum is not at least buflen + lsbyte bytes
+ * long.
+ *
+ * This code assumes that the compiler has the minimal intelligence
+ * neded to optimize divides and modulo operations on an unsigned data
+ * type with a power of two.
+ */
+void
+lbnExtractBigBytes_32(BNWORD32 const *n, unsigned char *buf,
+ unsigned lsbyte, unsigned buflen)
+{
+ BNWORD32 t = 0; /* Needed to shut up uninitialized var warnings */
+ unsigned shift;
+
+ lsbyte += buflen;
+
+ shift = (8 * lsbyte) % 32;
+ lsbyte /= (32/8); /* Convert to word offset */
+ BIGLITTLE(n -= lsbyte, n += lsbyte);
+
+ if (shift)
+ t = BIGLITTLE(n[-1],n[0]);
+
+ while (buflen--) {
+ if (!shift) {
+ t = BIGLITTLE(*n++,*--n);
+ shift = 32;
+ }
+ shift -= 8;
+ *buf++ = (unsigned char)(t>>shift);
+ }
+}
+
+/*
+ * Merge a big-endian array of bytes into a bignum array.
+ * The array had better be big enough. This is
+ * equivalent to extracting the entire bignum into a
+ * large byte array, copying the input buffer into the
+ * middle of it, and converting back to a bignum.
+ *
+ * The buf is "len" bytes long, and its *last* byte is at
+ * position "lsbyte" from the end of the bignum.
+ *
+ * Note that this is a pain to get right. Fortunately, it's hardly
+ * critical for efficiency.
+ */
+void
+lbnInsertBigBytes_32(BNWORD32 *n, unsigned char const *buf,
+ unsigned lsbyte, unsigned buflen)
+{
+ BNWORD32 t = 0; /* Shut up uninitialized varibale warnings */
+
+ lsbyte += buflen;
+
+ BIGLITTLE(n -= lsbyte/(32/8), n += lsbyte/(32/8));
+
+ /* Load up leading odd bytes */
+ if (lsbyte % (32/8)) {
+ t = BIGLITTLE(*--n,*n++);
+ t >>= (lsbyte * 8) % 32;
+ }
+
+ /* The main loop - merge into t, storing at each word boundary. */
+ while (buflen--) {
+ t = (t << 8) | *buf++;
+ if ((--lsbyte % (32/8)) == 0)
+ BIGLITTLE(*n++,*--n) = t;
+ }
+
+ /* Merge odd bytes in t into last word */
+ lsbyte = (lsbyte * 8) % 32;
+ if (lsbyte) {
+ t <<= lsbyte;
+ t |= (((BNWORD32)1 << lsbyte) - 1) & BIGLITTLE(n[0],n[-1]);
+ BIGLITTLE(n[0],n[-1]) = t;
+ }
+
+ return;
+}
+
+/*
+ * Returns a substring of the little-endian array of bytes representation
+ * of the bignum array based on two parameters, the least significant
+ * byte number (0 to start with the least significant byte) and the
+ * length. I.e. the number returned is a representation of
+ * (bn / 2^(8*lsbyte)) % 2 ^ (8*buflen).
+ *
+ * It is an error if the bignum is not at least buflen + lsbyte bytes
+ * long.
+ *
+ * This code assumes that the compiler has the minimal intelligence
+ * neded to optimize divides and modulo operations on an unsigned data
+ * type with a power of two.
+ */
+void
+lbnExtractLittleBytes_32(BNWORD32 const *n, unsigned char *buf,
+ unsigned lsbyte, unsigned buflen)
+{
+ BNWORD32 t = 0; /* Needed to shut up uninitialized var warnings */
+
+ BIGLITTLE(n -= lsbyte/(32/8), n += lsbyte/(32/8));
+
+ if (lsbyte % (32/8)) {
+ t = BIGLITTLE(*--n,*n++);
+ t >>= (lsbyte % (32/8)) * 8 ;
+ }
+
+ while (buflen--) {
+ if ((lsbyte++ % (32/8)) == 0)
+ t = BIGLITTLE(*--n,*n++);
+ *buf++ = (unsigned char)t;
+ t >>= 8;
+ }
+}
+
+/*
+ * Merge a little-endian array of bytes into a bignum array.
+ * The array had better be big enough. This is
+ * equivalent to extracting the entire bignum into a
+ * large byte array, copying the input buffer into the
+ * middle of it, and converting back to a bignum.
+ *
+ * The buf is "len" bytes long, and its first byte is at
+ * position "lsbyte" from the end of the bignum.
+ *
+ * Note that this is a pain to get right. Fortunately, it's hardly
+ * critical for efficiency.
+ */
+void
+lbnInsertLittleBytes_32(BNWORD32 *n, unsigned char const *buf,
+ unsigned lsbyte, unsigned buflen)
+{
+ BNWORD32 t = 0; /* Shut up uninitialized varibale warnings */
+
+ /* Move to most-significant end */
+ lsbyte += buflen;
+ buf += buflen;
+
+ BIGLITTLE(n -= lsbyte/(32/8), n += lsbyte/(32/8));
+
+ /* Load up leading odd bytes */
+ if (lsbyte % (32/8)) {
+ t = BIGLITTLE(*--n,*n++);
+ t >>= (lsbyte * 8) % 32;
+ }
+
+ /* The main loop - merge into t, storing at each word boundary. */
+ while (buflen--) {
+ t = (t << 8) | *--buf;
+ if ((--lsbyte % (32/8)) == 0)
+ BIGLITTLE(*n++,*--n) = t;
+ }
+
+ /* Merge odd bytes in t into last word */
+ lsbyte = (lsbyte * 8) % 32;
+ if (lsbyte) {
+ t <<= lsbyte;
+ t |= (((BNWORD32)1 << lsbyte) - 1) & BIGLITTLE(n[0],n[-1]);
+ BIGLITTLE(n[0],n[-1]) = t;
+ }
+
+ return;
+}
+
+#ifdef DEADCODE /* This was a precursor to the more flexible lbnExtractBytes */
+/*
+ * Convert a big-endian array of bytes to a bignum.
+ * Returns the number of words in the bignum.
+ * Note the expression "32/8" for the number of bytes per word.
+ * This is so the word-size adjustment will work.
+ */
+unsigned
+lbnFromBytes_32(BNWORD32 *a, unsigned char const *b, unsigned blen)
+{
+ BNWORD32 t;
+ unsigned alen = (blen + (32/8-1))/(32/8);
+ BIGLITTLE(a -= alen, a += alen);
+
+ while (blen) {
+ t = 0;
+ do {
+ t = t << 8 | *b++;
+ } while (--blen & (32/8-1));
+ BIGLITTLE(*a++,*--a) = t;
+ }
+ return alen;
+}
+#endif
+
+/*
+ * Computes the GCD of a and b. Modifies both arguments; when it returns,
+ * one of them is the GCD and the other is trash. The return value
+ * indicates which: 0 for a, and 1 for b. The length of the retult is
+ * returned in rlen. Both inputs must have one extra word of precision.
+ * alen must be >= blen.
+ *
+ * TODO: use the binary algorithm (Knuth section 4.5.2, algorithm B).
+ * This is based on taking out common powers of 2, then repeatedly:
+ * gcd(2*u,v) = gcd(u,2*v) = gcd(u,v) - isolated powers of 2 can be deleted.
+ * gcd(u,v) = gcd(u-v,v) - the numbers can be easily reduced.
+ * It gets less reduction per step, but the steps are much faster than
+ * the division case.
+ */
+int
+lbnGcd_32(BNWORD32 *a, unsigned alen, BNWORD32 *b, unsigned blen,
+ unsigned *rlen)
+{
+#if BNYIELD
+ int y;
+#endif
+ assert(alen >= blen);
+
+ while (blen != 0) {
+ (void)lbnDiv_32(BIGLITTLE(a-blen,a+blen), a, alen, b, blen);
+ alen = lbnNorm_32(a, blen);
+ if (alen == 0) {
+ *rlen = blen;
+ return 1;
+ }
+ (void)lbnDiv_32(BIGLITTLE(b-alen,b+alen), b, blen, a, alen);
+ blen = lbnNorm_32(b, alen);
+#if BNYIELD
+ if (bnYield && (y = bnYield()) < 0)
+ return y;
+#endif
+ }
+ *rlen = alen;
+ return 0;
+}
+
+/*
+ * Invert "a" modulo "mod" using the extended Euclidean algorithm.
+ * Note that this only computes one of the cosequences, and uses the
+ * theorem that the signs flip every step and the absolute value of
+ * the cosequence values are always bounded by the modulus to avoid
+ * having to work with negative numbers.
+ * gcd(a,mod) had better equal 1. Returns 1 if the GCD is NOT 1.
+ * a must be one word longer than "mod". It is overwritten with the
+ * result.
+ * TODO: Use Richard Schroeppel's *much* faster algorithm.
+ */
+int
+lbnInv_32(BNWORD32 *a, unsigned alen, BNWORD32 const *mod, unsigned mlen)
+{
+ BNWORD32 *b; /* Hold a copy of mod during GCD reduction */
+ BNWORD32 *p; /* Temporary for products added to t0 and t1 */
+ BNWORD32 *t0, *t1; /* Inverse accumulators */
+ BNWORD32 cy;
+ unsigned blen, t0len, t1len, plen;
+ int y;
+
+ alen = lbnNorm_32(a, alen);
+ if (!alen)
+ return 1; /* No inverse */
+
+ mlen = lbnNorm_32(mod, mlen);
+
+ assert (alen <= mlen);
+
+ /* Inverse of 1 is 1 */
+ if (alen == 1 && BIGLITTLE(a[-1],a[0]) == 1) {
+ lbnZero_32(BIGLITTLE(a-alen,a+alen), mlen-alen);
+ return 0;
+ }
+
+ /* Allocate a pile of space */
+ LBNALLOC(b, BNWORD32, mlen+1);
+ if (b) {
+ /*
+ * Although products are guaranteed to always be less than the
+ * modulus, it can involve multiplying two 3-word numbers to
+ * get a 5-word result, requiring a 6th word to store a 0
+ * temporarily. Thus, mlen + 1.
+ */
+ LBNALLOC(p, BNWORD32, mlen+1);
+ if (p) {
+ LBNALLOC(t0, BNWORD32, mlen);
+ if (t0) {
+ LBNALLOC(t1, BNWORD32, mlen);
+ if (t1)
+ goto allocated;
+ LBNFREE(t0, mlen);
+ }
+ LBNFREE(p, mlen+1);
+ }
+ LBNFREE(b, mlen+1);
+ }
+ return -1;
+
+allocated:
+
+ /* Set t0 to 1 */
+ t0len = 1;
+ BIGLITTLE(t0[-1],t0[0]) = 1;
+
+ /* b = mod */
+ lbnCopy_32(b, mod, mlen);
+ /* blen = mlen (implicitly) */
+
+ /* t1 = b / a; b = b % a */
+ cy = lbnDiv_32(t1, b, mlen, a, alen);
+ *(BIGLITTLE(t1-(mlen-alen)-1,t1+(mlen-alen))) = cy;
+ t1len = lbnNorm_32(t1, mlen-alen+1);
+ blen = lbnNorm_32(b, alen);
+
+ /* while (b > 1) */
+ while (blen > 1 || BIGLITTLE(b[-1],b[0]) != (BNWORD32)1) {
+ /* q = a / b; a = a % b; */
+ if (alen < blen || (alen == blen && lbnCmp_32(a, a, alen) < 0))
+ assert(0);
+ cy = lbnDiv_32(BIGLITTLE(a-blen,a+blen), a, alen, b, blen);
+ *(BIGLITTLE(a-alen-1,a+alen)) = cy;
+ plen = lbnNorm_32(BIGLITTLE(a-blen,a+blen), alen-blen+1);
+ assert(plen);
+ alen = lbnNorm_32(a, blen);
+ if (!alen)
+ goto failure; /* GCD not 1 */
+
+ /* t0 += q * t1; */
+ assert(plen+t1len <= mlen+1);
+ lbnMul_32(p, BIGLITTLE(a-blen,a+blen), plen, t1, t1len);
+ plen = lbnNorm_32(p, plen + t1len);
+ assert(plen <= mlen);
+ if (plen > t0len) {
+ lbnZero_32(BIGLITTLE(t0-t0len,t0+t0len), plen-t0len);
+ t0len = plen;
+ }
+ cy = lbnAddN_32(t0, p, plen);
+ if (cy) {
+ if (t0len > plen) {
+ cy = lbnAdd1_32(BIGLITTLE(t0-plen,t0+plen),
+ t0len-plen, cy);
+ }
+ if (cy) {
+ BIGLITTLE(t0[-t0len-1],t0[t0len]) = cy;
+ t0len++;
+ }
+ }
+
+ /* if (a <= 1) return a ? t0 : FAIL; */
+ if (alen <= 1 && BIGLITTLE(a[-1],a[0]) == (BNWORD32)1) {
+ if (alen == 0)
+ goto failure; /* FAIL */
+ assert(t0len <= mlen);
+ lbnCopy_32(a, t0, t0len);
+ lbnZero_32(BIGLITTLE(a-t0len, a+t0len), mlen-t0len);
+ goto success;
+ }
+
+ /* q = b / a; b = b % a; */
+ if (blen < alen || (blen == alen && lbnCmp_32(b, a, alen) < 0))
+ assert(0);
+ cy = lbnDiv_32(BIGLITTLE(b-alen,b+alen), b, blen, a, alen);
+ *(BIGLITTLE(b-blen-1,b+blen)) = cy;
+ plen = lbnNorm_32(BIGLITTLE(b-alen,b+alen), blen-alen+1);
+ assert(plen);
+ blen = lbnNorm_32(b, alen);
+ if (!blen)
+ goto failure; /* GCD not 1 */
+
+ /* t1 += q * t0; */
+ assert(plen+t0len <= mlen+1);
+ lbnMul_32(p, BIGLITTLE(b-alen,b+alen), plen, t0, t0len);
+ plen = lbnNorm_32(p, plen + t0len);
+ assert(plen <= mlen);
+ if (plen > t1len) {
+ lbnZero_32(BIGLITTLE(t1-t1len,t1+t1len), plen-t1len);
+ t1len = plen;
+ }
+ cy = lbnAddN_32(t1, p, plen);
+ if (cy) {
+ if (t1len > plen) {
+ cy = lbnAdd1_32(BIGLITTLE(t1-plen,t0+plen),
+ t1len-plen, cy);
+ }
+ if (cy) {
+ BIGLITTLE(t1[-t1len-1],t1[t1len]) = cy;
+ t1len++;
+ }
+ }
+#if BNYIELD
+ if (bnYield && (y = bnYield() < 0))
+ goto yield;
+#endif
+ }
+
+ if (!blen)
+ goto failure; /* gcd(a, mod) != 1 -- FAIL */
+
+ /* return mod-t1 */
+ lbnCopy_32(a, mod, mlen);
+ assert(t1len <= mlen);
+ cy = lbnSubN_32(a, t1, t1len);
+ if (cy) {
+ assert(mlen > t1len);
+ cy = lbnSub1_32(BIGLITTLE(a-t1len, a+t1len), mlen-t1len, cy);
+ assert(!cy);
+ }
+
+success:
+ LBNFREE(t1, mlen);
+ LBNFREE(t0, mlen);
+ LBNFREE(p, mlen+1);
+ LBNFREE(b, mlen+1);
+
+ return 0;
+
+failure: /* GCD is not 1 - no inverse exists! */
+ y = 1;
+#if BNYIELD
+yield:
+#endif
+ LBNFREE(t1, mlen);
+ LBNFREE(t0, mlen);
+ LBNFREE(p, mlen+1);
+ LBNFREE(b, mlen+1);
+
+ return y;
+}
+
+/*
+ * Precompute powers of "a" mod "mod". Compute them every "bits"
+ * for "n" steps. This is sufficient to compute powers of g with
+ * exponents up to n*bits bits long, i.e. less than 2^(n*bits).
+ *
+ * This assumes that the caller has already initialized "array" to point
+ * to "n" buffers of size "mlen".
+ */
+int
+lbnBasePrecompBegin_32(BNWORD32 **array, unsigned n, unsigned bits,
+ BNWORD32 const *g, unsigned glen, BNWORD32 *mod, unsigned mlen)
+{
+ BNWORD32 *a, *b; /* Temporary double-width accumulators */
+ BNWORD32 *a1; /* Pointer to high half of a*/
+ BNWORD32 inv; /* Montgomery inverse of LSW of mod */
+ BNWORD32 *t;
+ unsigned i;
+
+ glen = lbnNorm_32(g, glen);
+ assert(glen);
+
+ assert (mlen == lbnNorm_32(mod, mlen));
+ assert (glen <= mlen);
+
+ /* Allocate two temporary buffers, and the array slots */
+ LBNALLOC(a, BNWORD32, mlen*2);
+ if (!a)
+ return -1;
+ LBNALLOC(b, BNWORD32, mlen*2);
+ if (!b) {
+ LBNFREE(a, 2*mlen);
+ return -1;
+ }
+
+ /* Okay, all ready */
+
+ /* Convert n to Montgomery form */
+ inv = BIGLITTLE(mod[-1],mod[0]); /* LSW of modulus */
+ assert(inv & 1); /* Modulus must be odd */
+ inv = lbnMontInv1_32(inv);
+ /* Move g up "mlen" words into a (clearing the low mlen words) */
+ a1 = BIGLITTLE(a-mlen,a+mlen);
+ lbnCopy_32(a1, g, glen);
+ lbnZero_32(a, mlen);
+
+ /* Do the division - dump the quotient into the high-order words */
+ (void)lbnDiv_32(a1, a, mlen+glen, mod, mlen);
+
+ /* Copy the first value into the array */
+ t = *array;
+ lbnCopy_32(t, a, mlen);
+ a1 = a; /* This first value is *not* shifted up */
+
+ /* Now compute the remaining n-1 array entries */
+ assert(bits);
+ assert(n);
+ while (--n) {
+ i = bits;
+ do {
+ /* Square a1 into b1 */
+ lbnMontSquare_32(b, a1, mod, mlen, inv);
+ t = b; b = a; a = t;
+ a1 = BIGLITTLE(a-mlen, a+mlen);
+ } while (--i);
+ t = *++array;
+ lbnCopy_32(t, a1, mlen);
+ }
+
+ /* Hooray, we're done. */
+ LBNFREE(b, 2*mlen);
+ LBNFREE(a, 2*mlen);
+ return 0;
+}
+
+/*
+ * result = base^exp (mod mod). "array" is a an array of pointers
+ * to procomputed powers of base, each 2^bits apart. (I.e. array[i]
+ * is base^(2^(i*bits))).
+ *
+ * The algorithm consists of:
+ * a = b = (powers of g to be raised to the power 2^bits-1)
+ * a *= b *= (powers of g to be raised to the power 2^bits-2)
+ * ...
+ * a *= b *= (powers of g to be raised to the power 1)
+ *
+ * All we do is walk the exponent 2^bits-1 times in groups of "bits" bits,
+ */
+int
+lbnBasePrecompExp_32(BNWORD32 *result, BNWORD32 const * const *array,
+ unsigned bits, BNWORD32 const *exp, unsigned elen,
+ BNWORD32 const *mod, unsigned mlen)
+{
+ BNWORD32 *a, *b, *c, *t;
+ BNWORD32 *a1, *b1;
+ int anull, bnull; /* Null flags: values are implicitly 1 */
+ unsigned i, j; /* Loop counters */
+ unsigned mask; /* Exponent bits to examime */
+ BNWORD32 const *eptr; /* Pointer into exp */
+ BNWORD32 buf, curbits, nextword; /* Bit-buffer varaibles */
+ BNWORD32 inv; /* Inverse of LSW of modulus */
+ unsigned ewords; /* Words of exponent left */
+ int bufbits; /* Number of valid bits */
+ int y = 0;
+
+ mlen = lbnNorm_32(mod, mlen);
+ assert (mlen);
+
+ elen = lbnNorm_32(exp, elen);
+ if (!elen) {
+ lbnZero_32(result, mlen);
+ BIGLITTLE(result[-1],result[0]) = 1;
+ return 0;
+ }
+ /*
+ * This could be precomputed, but it's so cheap, and it would require
+ * making the precomputation structure word-size dependent.
+ */
+ inv = lbnMontInv1_32(mod[BIGLITTLE(-1,0)]); /* LSW of modulus */
+
+ assert(elen);
+
+ /*
+ * Allocate three temporary buffers. The current numbers generally
+ * live in the upper halves of these buffers.
+ */
+ LBNALLOC(a, BNWORD32, mlen*2);
+ if (a) {
+ LBNALLOC(b, BNWORD32, mlen*2);
+ if (b) {
+ LBNALLOC(c, BNWORD32, mlen*2);
+ if (c)
+ goto allocated;
+ LBNFREE(b, 2*mlen);
+ }
+ LBNFREE(a, 2*mlen);
+ }
+ return -1;
+
+allocated:
+
+ anull = bnull = 1;
+
+ mask = (1u<<bits) - 1;
+ for (i = mask; i; --i) {
+ /* Set up bit buffer for walking the exponent */
+ eptr = exp;
+ buf = BIGLITTLE(*--eptr, *eptr++);
+ ewords = elen-1;
+ bufbits = 32;
+ for (j = 0; ewords || buf; j++) {
+ /* Shift down current buffer */
+ curbits = buf;
+ buf >>= bits;
+ /* If necessary, add next word */
+ bufbits -= bits;
+ if (bufbits < 0 && ewords > 0) {
+ nextword = BIGLITTLE(*--eptr, *eptr++);
+ ewords--;
+ curbits |= nextword << (bufbits+bits);
+ buf = nextword >> -bufbits;
+ bufbits += 32;
+ }
+ /* If appropriate, multiply b *= array[j] */
+ if ((curbits & mask) == i) {
+ BNWORD32 const *d = array[j];
+
+ b1 = BIGLITTLE(b-mlen-1,b+mlen);
+ if (bnull) {
+ lbnCopy_32(b1, d, mlen);
+ bnull = 0;
+ } else {
+ lbnMontMul_32(c, b1, d, mod, mlen, inv);
+ t = c; c = b; b = t;
+ }
+#if BNYIELD
+ if (bnYield && (y = bnYield() < 0))
+ goto yield;
+#endif
+ }
+ }
+
+ /* Multiply a *= b */
+ if (!bnull) {
+ a1 = BIGLITTLE(a-mlen-1,a+mlen);
+ b1 = BIGLITTLE(b-mlen-1,b+mlen);
+ if (anull) {
+ lbnCopy_32(a1, b1, mlen);
+ anull = 0;
+ } else {
+ lbnMontMul_32(c, a1, b1, mod, mlen, inv);
+ t = c; c = a; a = t;
+ }
+ }
+ }
+
+ assert(!anull); /* If it were, elen would have been 0 */
+
+ /* Convert out of Montgomery form and return */
+ a1 = BIGLITTLE(a-mlen-1,a+mlen);
+ lbnCopy_32(a, a1, mlen);
+ lbnZero_32(a1, mlen);
+ lbnMontReduce_32(a, mod, mlen, inv);
+ lbnCopy_32(result, a1, mlen);
+
+#if BNYIELD
+yield:
+#endif
+ LBNFREE(c, 2*mlen);
+ LBNFREE(b, 2*mlen);
+ LBNFREE(a, 2*mlen);
+
+ return y;
+}
+
+/*
+ * result = base1^exp1 *base2^exp2 (mod mod). "array1" and "array2" are
+ * arrays of pointers to procomputed powers of the corresponding bases,
+ * each 2^bits apart. (I.e. array1[i] is base1^(2^(i*bits))).
+ *
+ * Bits must be the same in both. (It could be made adjustable, but it's
+ * a bit of a pain. Just make them both equal to the larger one.)
+ *
+ * The algorithm consists of:
+ * a = b = (powers of base1 and base2 to be raised to the power 2^bits-1)
+ * a *= b *= (powers of base1 and base2 to be raised to the power 2^bits-2)
+ * ...
+ * a *= b *= (powers of base1 and base2 to be raised to the power 1)
+ *
+ * All we do is walk the exponent 2^bits-1 times in groups of "bits" bits,
+ */
+int
+lbnDoubleBasePrecompExp_32(BNWORD32 *result, unsigned bits,
+ BNWORD32 const * const *array1, BNWORD32 const *exp1, unsigned elen1,
+ BNWORD32 const * const *array2, BNWORD32 const *exp2,
+ unsigned elen2, BNWORD32 const *mod, unsigned mlen)
+{
+ BNWORD32 *a, *b, *c, *t;
+ BNWORD32 *a1, *b1;
+ int anull, bnull; /* Null flags: values are implicitly 1 */
+ unsigned i, j, k; /* Loop counters */
+ unsigned mask; /* Exponent bits to examime */
+ BNWORD32 const *eptr; /* Pointer into exp */
+ BNWORD32 buf, curbits, nextword; /* Bit-buffer varaibles */
+ BNWORD32 inv; /* Inverse of LSW of modulus */
+ unsigned ewords; /* Words of exponent left */
+ int bufbits; /* Number of valid bits */
+ int y = 0;
+ BNWORD32 const * const *array;
+
+ mlen = lbnNorm_32(mod, mlen);
+ assert (mlen);
+
+ elen1 = lbnNorm_32(exp1, elen1);
+ if (!elen1) {
+ return lbnBasePrecompExp_32(result, array2, bits, exp2, elen2,
+ mod, mlen);
+ }
+ elen2 = lbnNorm_32(exp2, elen2);
+ if (!elen2) {
+ return lbnBasePrecompExp_32(result, array1, bits, exp1, elen1,
+ mod, mlen);
+ }
+ /*
+ * This could be precomputed, but it's so cheap, and it would require
+ * making the precomputation structure word-size dependent.
+ */
+ inv = lbnMontInv1_32(mod[BIGLITTLE(-1,0)]); /* LSW of modulus */
+
+ assert(elen1);
+ assert(elen2);
+
+ /*
+ * Allocate three temporary buffers. The current numbers generally
+ * live in the upper halves of these buffers.
+ */
+ LBNALLOC(a, BNWORD32, mlen*2);
+ if (a) {
+ LBNALLOC(b, BNWORD32, mlen*2);
+ if (b) {
+ LBNALLOC(c, BNWORD32, mlen*2);
+ if (c)
+ goto allocated;
+ LBNFREE(b, 2*mlen);
+ }
+ LBNFREE(a, 2*mlen);
+ }
+ return -1;
+
+allocated:
+
+ anull = bnull = 1;
+
+ mask = (1u<<bits) - 1;
+ for (i = mask; i; --i) {
+ /* Walk each exponent in turn */
+ for (k = 0; k < 2; k++) {
+ /* Set up the exponent for walking */
+ array = k ? array2 : array1;
+ eptr = k ? exp2 : exp1;
+ ewords = (k ? elen2 : elen1) - 1;
+ /* Set up bit buffer for walking the exponent */
+ buf = BIGLITTLE(*--eptr, *eptr++);
+ bufbits = 32;
+ for (j = 0; ewords || buf; j++) {
+ /* Shift down current buffer */
+ curbits = buf;
+ buf >>= bits;
+ /* If necessary, add next word */
+ bufbits -= bits;
+ if (bufbits < 0 && ewords > 0) {
+ nextword = BIGLITTLE(*--eptr, *eptr++);
+ ewords--;
+ curbits |= nextword << (bufbits+bits);
+ buf = nextword >> -bufbits;
+ bufbits += 32;
+ }
+ /* If appropriate, multiply b *= array[j] */
+ if ((curbits & mask) == i) {
+ BNWORD32 const *d = array[j];
+
+ b1 = BIGLITTLE(b-mlen-1,b+mlen);
+ if (bnull) {
+ lbnCopy_32(b1, d, mlen);
+ bnull = 0;
+ } else {
+ lbnMontMul_32(c, b1, d, mod, mlen, inv);
+ t = c; c = b; b = t;
+ }
+#if BNYIELD
+ if (bnYield && (y = bnYield() < 0))
+ goto yield;
+#endif
+ }
+ }
+ }
+
+ /* Multiply a *= b */
+ if (!bnull) {
+ a1 = BIGLITTLE(a-mlen-1,a+mlen);
+ b1 = BIGLITTLE(b-mlen-1,b+mlen);
+ if (anull) {
+ lbnCopy_32(a1, b1, mlen);
+ anull = 0;
+ } else {
+ lbnMontMul_32(c, a1, b1, mod, mlen, inv);
+ t = c; c = a; a = t;
+ }
+ }
+ }
+
+ assert(!anull); /* If it were, elen would have been 0 */
+
+ /* Convert out of Montgomery form and return */
+ a1 = BIGLITTLE(a-mlen-1,a+mlen);
+ lbnCopy_32(a, a1, mlen);
+ lbnZero_32(a1, mlen);
+ lbnMontReduce_32(a, mod, mlen, inv);
+ lbnCopy_32(result, a1, mlen);
+
+#if BNYIELD
+yield:
+#endif
+ LBNFREE(c, 2*mlen);
+ LBNFREE(b, 2*mlen);
+ LBNFREE(a, 2*mlen);
+
+ return y;
+}