+ for (i = 0; i < 2 * mlen; i++)
+ a[i] = 0;
+ sfree(a);
+ for (i = 0; i < scratchlen; i++)
+ scratch[i] = 0;
+ sfree(scratch);
+ for (i = 0; i < 2 * mlen; i++)
+ b[i] = 0;
+ sfree(b);
+ for (i = 0; i < mlen; i++)
+ m[i] = 0;
+ sfree(m);
+ for (i = 0; i < mlen; i++)
+ n[i] = 0;
+ sfree(n);
+
+ freebn(base);
+
+ return result;
+}
+
+/*
+ * Compute (base ^ exp) % mod. Uses the Montgomery multiplication
+ * technique where possible, falling back to modpow_simple otherwise.
+ */
+Bignum modpow(Bignum base_in, Bignum exp, Bignum mod)
+{
+ BignumInt *a, *b, *x, *n, *mninv, *scratch;
+ int len, scratchlen, i, j;
+ Bignum base, base2, r, rn, inv, result;
+
+ /*
+ * The most significant word of mod needs to be non-zero. It
+ * should already be, but let's make sure.
+ */
+ assert(mod[mod[0]] != 0);
+
+ /*
+ * mod had better be odd, or we can't do Montgomery multiplication
+ * using a power of two at all.
+ */
+ if (!(mod[1] & 1))
+ return modpow_simple(base_in, exp, mod);
+
+ /*
+ * Make sure the base is smaller than the modulus, by reducing
+ * it modulo the modulus if not.
+ */
+ base = bigmod(base_in, mod);
+
+ /*
+ * Compute the inverse of n mod r, for monty_reduce. (In fact we
+ * want the inverse of _minus_ n mod r, but we'll sort that out
+ * below.)
+ */
+ len = mod[0];
+ r = bn_power_2(BIGNUM_INT_BITS * len);
+ inv = modinv(mod, r);
+
+ /*
+ * Multiply the base by r mod n, to get it into Montgomery
+ * representation.
+ */
+ base2 = modmul(base, r, mod);
+ freebn(base);
+ base = base2;
+
+ rn = bigmod(r, mod); /* r mod n, i.e. Montgomerified 1 */
+
+ freebn(r); /* won't need this any more */
+
+ /*
+ * Set up internal arrays of the right lengths, in big-endian
+ * format, containing the base, the modulus, and the modulus's
+ * inverse.
+ */
+ n = snewn(len, BignumInt);
+ for (j = 0; j < len; j++)
+ n[len - 1 - j] = mod[j + 1];
+
+ mninv = snewn(len, BignumInt);
+ for (j = 0; j < len; j++)
+ mninv[len - 1 - j] = (j < (int)inv[0] ? inv[j + 1] : 0);
+ freebn(inv); /* we don't need this copy of it any more */
+ /* Now negate mninv mod r, so it's the inverse of -n rather than +n. */
+ x = snewn(len, BignumInt);
+ for (j = 0; j < len; j++)
+ x[j] = 0;
+ internal_sub(x, mninv, mninv, len);
+
+ /* x = snewn(len, BignumInt); */ /* already done above */
+ for (j = 0; j < len; j++)
+ x[len - 1 - j] = (j < (int)base[0] ? base[j + 1] : 0);
+ freebn(base); /* we don't need this copy of it any more */
+
+ a = snewn(2*len, BignumInt);
+ b = snewn(2*len, BignumInt);
+ for (j = 0; j < len; j++)
+ a[2*len - 1 - j] = (j < (int)rn[0] ? rn[j + 1] : 0);
+ freebn(rn);
+
+ /* Scratch space for multiplies */
+ scratchlen = 3*len + mul_compute_scratch(len);
+ scratch = snewn(scratchlen, BignumInt);
+
+ /* Skip leading zero bits of exp. */
+ i = 0;
+ j = BIGNUM_INT_BITS-1;
+ while (i < (int)exp[0] && (exp[exp[0] - i] & (1 << j)) == 0) {
+ j--;
+ if (j < 0) {
+ i++;
+ j = BIGNUM_INT_BITS-1;
+ }
+ }
+
+ /* Main computation */
+ while (i < (int)exp[0]) {
+ while (j >= 0) {
+ internal_mul(a + len, a + len, b, len, scratch);
+ monty_reduce(b, n, mninv, scratch, len);
+ if ((exp[exp[0] - i] & (1 << j)) != 0) {
+ internal_mul(b + len, x, a, len, scratch);
+ monty_reduce(a, n, mninv, scratch, len);
+ } else {
+ BignumInt *t;
+ t = a;
+ a = b;
+ b = t;
+ }
+ j--;
+ }
+ i++;
+ j = BIGNUM_INT_BITS-1;
+ }
+
+ /*
+ * Final monty_reduce to get back from the adjusted Montgomery
+ * representation.
+ */
+ monty_reduce(a, n, mninv, scratch, len);
+
+ /* Copy result to buffer */
+ result = newbn(mod[0]);
+ for (i = 0; i < len; i++)
+ result[result[0] - i] = a[i + len];
+ while (result[0] > 1 && result[result[0]] == 0)
+ result[0]--;
+
+ /* Free temporary arrays */
+ for (i = 0; i < scratchlen; i++)
+ scratch[i] = 0;
+ sfree(scratch);
+ for (i = 0; i < 2 * len; i++)
+ a[i] = 0;
+ sfree(a);
+ for (i = 0; i < 2 * len; i++)
+ b[i] = 0;
+ sfree(b);
+ for (i = 0; i < len; i++)
+ mninv[i] = 0;
+ sfree(mninv);
+ for (i = 0; i < len; i++)
+ n[i] = 0;
+ sfree(n);
+ for (i = 0; i < len; i++)
+ x[i] = 0;
+ sfree(x);
+
+ return result;