summaryrefslogtreecommitdiffstats
path: root/libtommath/bn_s_mp_exptmod.c
diff options
context:
space:
mode:
Diffstat (limited to 'libtommath/bn_s_mp_exptmod.c')
-rw-r--r--libtommath/bn_s_mp_exptmod.c236
1 files changed, 236 insertions, 0 deletions
diff --git a/libtommath/bn_s_mp_exptmod.c b/libtommath/bn_s_mp_exptmod.c
new file mode 100644
index 0000000..01a766f
--- /dev/null
+++ b/libtommath/bn_s_mp_exptmod.c
@@ -0,0 +1,236 @@
+#include <tommath.h>
+#ifdef BN_S_MP_EXPTMOD_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis
+ *
+ * LibTomMath is a library that provides multiple-precision
+ * integer arithmetic as well as number theoretic functionality.
+ *
+ * The library was designed directly after the MPI library by
+ * Michael Fromberger but has been written from scratch with
+ * additional optimizations in place.
+ *
+ * The library is free for all purposes without any express
+ * guarantee it works.
+ *
+ * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
+ */
+
+#ifdef MP_LOW_MEM
+ #define TAB_SIZE 32
+#else
+ #define TAB_SIZE 256
+#endif
+
+int s_mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y)
+{
+ mp_int M[TAB_SIZE], res, mu;
+ mp_digit buf;
+ int err, bitbuf, bitcpy, bitcnt, mode, digidx, x, y, winsize;
+
+ /* find window size */
+ x = mp_count_bits (X);
+ if (x <= 7) {
+ winsize = 2;
+ } else if (x <= 36) {
+ winsize = 3;
+ } else if (x <= 140) {
+ winsize = 4;
+ } else if (x <= 450) {
+ winsize = 5;
+ } else if (x <= 1303) {
+ winsize = 6;
+ } else if (x <= 3529) {
+ winsize = 7;
+ } else {
+ winsize = 8;
+ }
+
+#ifdef MP_LOW_MEM
+ if (winsize > 5) {
+ winsize = 5;
+ }
+#endif
+
+ /* init M array */
+ /* init first cell */
+ if ((err = mp_init(&M[1])) != MP_OKAY) {
+ return err;
+ }
+
+ /* now init the second half of the array */
+ for (x = 1<<(winsize-1); x < (1 << winsize); x++) {
+ if ((err = mp_init(&M[x])) != MP_OKAY) {
+ for (y = 1<<(winsize-1); y < x; y++) {
+ mp_clear (&M[y]);
+ }
+ mp_clear(&M[1]);
+ return err;
+ }
+ }
+
+ /* create mu, used for Barrett reduction */
+ if ((err = mp_init (&mu)) != MP_OKAY) {
+ goto LBL_M;
+ }
+ if ((err = mp_reduce_setup (&mu, P)) != MP_OKAY) {
+ goto LBL_MU;
+ }
+
+ /* create M table
+ *
+ * The M table contains powers of the base,
+ * e.g. M[x] = G**x mod P
+ *
+ * The first half of the table is not
+ * computed though accept for M[0] and M[1]
+ */
+ if ((err = mp_mod (G, P, &M[1])) != MP_OKAY) {
+ goto LBL_MU;
+ }
+
+ /* compute the value at M[1<<(winsize-1)] by squaring
+ * M[1] (winsize-1) times
+ */
+ if ((err = mp_copy (&M[1], &M[1 << (winsize - 1)])) != MP_OKAY) {
+ goto LBL_MU;
+ }
+
+ for (x = 0; x < (winsize - 1); x++) {
+ if ((err = mp_sqr (&M[1 << (winsize - 1)],
+ &M[1 << (winsize - 1)])) != MP_OKAY) {
+ goto LBL_MU;
+ }
+ if ((err = mp_reduce (&M[1 << (winsize - 1)], P, &mu)) != MP_OKAY) {
+ goto LBL_MU;
+ }
+ }
+
+ /* create upper table, that is M[x] = M[x-1] * M[1] (mod P)
+ * for x = (2**(winsize - 1) + 1) to (2**winsize - 1)
+ */
+ for (x = (1 << (winsize - 1)) + 1; x < (1 << winsize); x++) {
+ if ((err = mp_mul (&M[x - 1], &M[1], &M[x])) != MP_OKAY) {
+ goto LBL_MU;
+ }
+ if ((err = mp_reduce (&M[x], P, &mu)) != MP_OKAY) {
+ goto LBL_MU;
+ }
+ }
+
+ /* setup result */
+ if ((err = mp_init (&res)) != MP_OKAY) {
+ goto LBL_MU;
+ }
+ mp_set (&res, 1);
+
+ /* set initial mode and bit cnt */
+ mode = 0;
+ bitcnt = 1;
+ buf = 0;
+ digidx = X->used - 1;
+ bitcpy = 0;
+ bitbuf = 0;
+
+ for (;;) {
+ /* grab next digit as required */
+ if (--bitcnt == 0) {
+ /* if digidx == -1 we are out of digits */
+ if (digidx == -1) {
+ break;
+ }
+ /* read next digit and reset the bitcnt */
+ buf = X->dp[digidx--];
+ bitcnt = (int) DIGIT_BIT;
+ }
+
+ /* grab the next msb from the exponent */
+ y = (buf >> (mp_digit)(DIGIT_BIT - 1)) & 1;
+ buf <<= (mp_digit)1;
+
+ /* if the bit is zero and mode == 0 then we ignore it
+ * These represent the leading zero bits before the first 1 bit
+ * in the exponent. Technically this opt is not required but it
+ * does lower the # of trivial squaring/reductions used
+ */
+ if (mode == 0 && y == 0) {
+ continue;
+ }
+
+ /* if the bit is zero and mode == 1 then we square */
+ if (mode == 1 && y == 0) {
+ if ((err = mp_sqr (&res, &res)) != MP_OKAY) {
+ goto LBL_RES;
+ }
+ if ((err = mp_reduce (&res, P, &mu)) != MP_OKAY) {
+ goto LBL_RES;
+ }
+ continue;
+ }
+
+ /* else we add it to the window */
+ bitbuf |= (y << (winsize - ++bitcpy));
+ mode = 2;
+
+ if (bitcpy == winsize) {
+ /* ok window is filled so square as required and multiply */
+ /* square first */
+ for (x = 0; x < winsize; x++) {
+ if ((err = mp_sqr (&res, &res)) != MP_OKAY) {
+ goto LBL_RES;
+ }
+ if ((err = mp_reduce (&res, P, &mu)) != MP_OKAY) {
+ goto LBL_RES;
+ }
+ }
+
+ /* then multiply */
+ if ((err = mp_mul (&res, &M[bitbuf], &res)) != MP_OKAY) {
+ goto LBL_RES;
+ }
+ if ((err = mp_reduce (&res, P, &mu)) != MP_OKAY) {
+ goto LBL_RES;
+ }
+
+ /* empty window and reset */
+ bitcpy = 0;
+ bitbuf = 0;
+ mode = 1;
+ }
+ }
+
+ /* if bits remain then square/multiply */
+ if (mode == 2 && bitcpy > 0) {
+ /* square then multiply if the bit is set */
+ for (x = 0; x < bitcpy; x++) {
+ if ((err = mp_sqr (&res, &res)) != MP_OKAY) {
+ goto LBL_RES;
+ }
+ if ((err = mp_reduce (&res, P, &mu)) != MP_OKAY) {
+ goto LBL_RES;
+ }
+
+ bitbuf <<= 1;
+ if ((bitbuf & (1 << winsize)) != 0) {
+ /* then multiply */
+ if ((err = mp_mul (&res, &M[1], &res)) != MP_OKAY) {
+ goto LBL_RES;
+ }
+ if ((err = mp_reduce (&res, P, &mu)) != MP_OKAY) {
+ goto LBL_RES;
+ }
+ }
+ }
+ }
+
+ mp_exch (&res, Y);
+ err = MP_OKAY;
+LBL_RES:mp_clear (&res);
+LBL_MU:mp_clear (&mu);
+LBL_M:
+ mp_clear(&M[1]);
+ for (x = 1<<(winsize-1); x < (1 << winsize); x++) {
+ mp_clear (&M[x]);
+ }
+ return err;
+}
+#endif
generated by cgit v0.12 at 2024-11-10 04:48:28 (GMT)