www.usr.com/support/gpl/USR9107_release.1.4.tar.gz
[bcm963xx.git] / userapps / opensource / sshd / libtommath / bn_s_mp_exptmod.c
index fcd0e01..597e877 100755 (executable)
-/* LibTomMath, multiple-precision integer library -- Tom St Denis\r
- *\r
- * LibTomMath is library that provides for multiple-precision\r
- * integer arithmetic as well as number theoretic functionality.\r
- *\r
- * The library is designed directly after the MPI library by\r
- * Michael Fromberger but has been written from scratch with\r
- * additional optimizations in place.\r
- *\r
- * The library is free for all purposes without any express\r
- * guarantee it works.\r
- *\r
- * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org\r
- */\r
-#include <tommath.h>\r
-\r
-int\r
-s_mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y)\r
-{\r
-  mp_int  M[256], res, mu;\r
-  mp_digit buf;\r
-  int     err, bitbuf, bitcpy, bitcnt, mode, digidx, x, y, winsize;\r
-\r
-  /* find window size */\r
-  x = mp_count_bits (X);\r
-  if (x <= 7) {\r
-    winsize = 2;\r
-  } else if (x <= 36) {\r
-    winsize = 3;\r
-  } else if (x <= 140) {\r
-    winsize = 4;\r
-  } else if (x <= 450) {\r
-    winsize = 5;\r
-  } else if (x <= 1303) {\r
-    winsize = 6;\r
-  } else if (x <= 3529) {\r
-    winsize = 7;\r
-  } else {\r
-    winsize = 8;\r
-  }\r
-\r
-#ifdef MP_LOW_MEM\r
-    if (winsize > 5) {\r
-       winsize = 5;\r
-    }\r
-#endif\r
-\r
-  /* init M array */\r
-  for (x = 0; x < (1 << winsize); x++) {\r
-    if ((err = mp_init_size (&M[x], 1)) != MP_OKAY) {\r
-      for (y = 0; y < x; y++) {\r
-        mp_clear (&M[y]);\r
-      }\r
-      return err;\r
-    }\r
-  }\r
-\r
-  /* create mu, used for Barrett reduction */\r
-  if ((err = mp_init (&mu)) != MP_OKAY) {\r
-    goto __M;\r
-  }\r
-  if ((err = mp_reduce_setup (&mu, P)) != MP_OKAY) {\r
-    goto __MU;\r
-  }\r
-\r
-  /* create M table\r
-   *\r
-   * The M table contains powers of the base, \r
-   * e.g. M[x] = G**x mod P\r
-   *\r
-   * The first half of the table is not \r
-   * computed though accept for M[0] and M[1]\r
-   */\r
-  if ((err = mp_mod (G, P, &M[1])) != MP_OKAY) {\r
-    goto __MU;\r
-  }\r
-\r
-  /* compute the value at M[1<<(winsize-1)] by squaring \r
-   * M[1] (winsize-1) times \r
-   */\r
-  if ((err = mp_copy (&M[1], &M[1 << (winsize - 1)])) != MP_OKAY) {\r
-    goto __MU;\r
-  }\r
-\r
-  for (x = 0; x < (winsize - 1); x++) {\r
-    if ((err = mp_sqr (&M[1 << (winsize - 1)], \r
-                       &M[1 << (winsize - 1)])) != MP_OKAY) {\r
-      goto __MU;\r
-    }\r
-    if ((err = mp_reduce (&M[1 << (winsize - 1)], P, &mu)) != MP_OKAY) {\r
-      goto __MU;\r
-    }\r
-  }\r
-\r
-  /* create upper table */\r
-  for (x = (1 << (winsize - 1)) + 1; x < (1 << winsize); x++) {\r
-    if ((err = mp_mul (&M[x - 1], &M[1], &M[x])) != MP_OKAY) {\r
-      goto __MU;\r
-    }\r
-    if ((err = mp_reduce (&M[x], P, &mu)) != MP_OKAY) {\r
-      goto __MU;\r
-    }\r
-  }\r
-\r
-  /* setup result */\r
-  if ((err = mp_init (&res)) != MP_OKAY) {\r
-    goto __MU;\r
-  }\r
-  mp_set (&res, 1);\r
-\r
-  /* set initial mode and bit cnt */\r
-  mode   = 0;\r
-  bitcnt = 1;\r
-  buf    = 0;\r
-  digidx = X->used - 1;\r
-  bitcpy = 0;\r
-  bitbuf = 0;\r
-\r
-  for (;;) {\r
-    /* grab next digit as required */\r
-    if (--bitcnt == 0) {\r
-      if (digidx == -1) {\r
-        break;\r
-      }\r
-      buf = X->dp[digidx--];\r
-      bitcnt = (int) DIGIT_BIT;\r
-    }\r
-\r
-    /* grab the next msb from the exponent */\r
-    y = (buf >> (mp_digit)(DIGIT_BIT - 1)) & 1;\r
-    buf <<= (mp_digit)1;\r
-\r
-    /* if the bit is zero and mode == 0 then we ignore it\r
-     * These represent the leading zero bits before the first 1 bit\r
-     * in the exponent.  Technically this opt is not required but it\r
-     * does lower the # of trivial squaring/reductions used\r
-     */\r
-    if (mode == 0 && y == 0)\r
-      continue;\r
-\r
-    /* if the bit is zero and mode == 1 then we square */\r
-    if (mode == 1 && y == 0) {\r
-      if ((err = mp_sqr (&res, &res)) != MP_OKAY) {\r
-        goto __RES;\r
-      }\r
-      if ((err = mp_reduce (&res, P, &mu)) != MP_OKAY) {\r
-        goto __RES;\r
-      }\r
-      continue;\r
-    }\r
-\r
-    /* else we add it to the window */\r
-    bitbuf |= (y << (winsize - ++bitcpy));\r
-    mode = 2;\r
-\r
-    if (bitcpy == winsize) {\r
-      /* ok window is filled so square as required and multiply  */\r
-      /* square first */\r
-      for (x = 0; x < winsize; x++) {\r
-        if ((err = mp_sqr (&res, &res)) != MP_OKAY) {\r
-          goto __RES;\r
-        }\r
-        if ((err = mp_reduce (&res, P, &mu)) != MP_OKAY) {\r
-          goto __RES;\r
-        }\r
-      }\r
-\r
-      /* then multiply */\r
-      if ((err = mp_mul (&res, &M[bitbuf], &res)) != MP_OKAY) {\r
-        goto __MU;\r
-      }\r
-      if ((err = mp_reduce (&res, P, &mu)) != MP_OKAY) {\r
-        goto __MU;\r
-      }\r
-\r
-      /* empty window and reset */\r
-      bitcpy = 0;\r
-      bitbuf = 0;\r
-      mode = 1;\r
-    }\r
-  }\r
-\r
-  /* if bits remain then square/multiply */\r
-  if (mode == 2 && bitcpy > 0) {\r
-    /* square then multiply if the bit is set */\r
-    for (x = 0; x < bitcpy; x++) {\r
-      if ((err = mp_sqr (&res, &res)) != MP_OKAY) {\r
-        goto __RES;\r
-      }\r
-      if ((err = mp_reduce (&res, P, &mu)) != MP_OKAY) {\r
-        goto __RES;\r
-      }\r
-\r
-      bitbuf <<= 1;\r
-      if ((bitbuf & (1 << winsize)) != 0) {\r
-        /* then multiply */\r
-        if ((err = mp_mul (&res, &M[1], &res)) != MP_OKAY) {\r
-          goto __RES;\r
-        }\r
-        if ((err = mp_reduce (&res, P, &mu)) != MP_OKAY) {\r
-          goto __RES;\r
-        }\r
-      }\r
-    }\r
-  }\r
-\r
-  mp_exch (&res, Y);\r
-  err = MP_OKAY;\r
-__RES:mp_clear (&res);\r
-__MU:mp_clear (&mu);\r
-__M:\r
-  for (x = 0; x < (1 << winsize); x++) {\r
-    mp_clear (&M[x]);\r
-  }\r
-  return err;\r
-}\r
+#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, int redmode)
+{
+  mp_int  M[TAB_SIZE], res, mu;
+  mp_digit buf;
+  int     err, bitbuf, bitcpy, bitcnt, mode, digidx, x, y, winsize;
+  int (*redux)(mp_int*,mp_int*,mp_int*);
+
+  /* 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 (redmode == 0) {
+     if ((err = mp_reduce_setup (&mu, P)) != MP_OKAY) {
+        goto LBL_MU;
+     }
+     redux = mp_reduce;
+  } else {
+     if ((err = mp_reduce_2k_setup_l (P, &mu)) != MP_OKAY) {
+        goto LBL_MU;
+     }
+     redux = mp_reduce_2k_l;
+  }    
+
+  /* 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++) {
+    /* square it */
+    if ((err = mp_sqr (&M[1 << (winsize - 1)], 
+                       &M[1 << (winsize - 1)])) != MP_OKAY) {
+      goto LBL_MU;
+    }
+
+    /* reduce modulo P */
+    if ((err = redux (&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 = redux (&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 = redux (&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 = redux (&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 = redux (&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 = redux (&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 = redux (&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