Browse files

added libtommath-0.12

  • Loading branch information...
1 parent 33c5019 commit 57354e11ac7611384bc4e6e8dd23eff281ce2a30 Tom St Denis committed with sjaeckel Feb 28, 2003
Showing with 5,878 additions and 4,610 deletions.
  1. +0 −2 b.bat
  2. +0 −3,908 bn.c
  3. BIN bn.pdf
  4. +56 −13 bn.tex
  5. +168 −0 bn_fast_mp_invmod.c
  6. +116 −0 bn_fast_mp_montgomery_reduce.c
  7. +113 −0 bn_fast_s_mp_mul_digs.c
  8. +86 −0 bn_fast_s_mp_mul_high_digs.c
  9. +112 −0 bn_fast_s_mp_sqr.c
  10. +31 −0 bn_mp_2expt.c
  11. +27 −0 bn_mp_abs.c
  12. +56 −0 bn_mp_add.c
  13. +33 −0 bn_mp_add_d.c
  14. +36 −0 bn_mp_addmod.c
  15. +51 −0 bn_mp_and.c
  16. +26 −0 bn_mp_clamp.c
  17. +33 −0 bn_mp_clear.c
  18. +30 −0 bn_mp_cmp.c
  19. +37 −0 bn_mp_cmp_d.c
  20. +40 −0 bn_mp_cmp_mag.c
  21. +48 −0 bn_mp_copy.c
  22. +35 −0 bn_mp_count_bits.c
  23. +200 −0 bn_mp_div.c
  24. +38 −0 bn_mp_div_2.c
  25. +78 −0 bn_mp_div_2d.c
  26. +44 −0 bn_mp_div_d.c
  27. +25 −0 bn_mp_exch.c
  28. +49 −0 bn_mp_expt_d.c
  29. +213 −0 bn_mp_exptmod.c
  30. +229 −0 bn_mp_exptmod_fast.c
  31. +118 −0 bn_mp_gcd.c
  32. +40 −0 bn_mp_grow.c
  33. +35 −0 bn_mp_init.c
  34. +28 −0 bn_mp_init_copy.c
  35. +33 −0 bn_mp_init_size.c
  36. +203 −0 bn_mp_invmod.c
  37. +114 −0 bn_mp_jacobi.c
  38. +142 −0 bn_mp_karatsuba_mul.c
  39. +106 −0 bn_mp_karatsuba_sqr.c
  40. +42 −0 bn_mp_lcm.c
  41. +46 −0 bn_mp_lshd.c
  42. +43 −0 bn_mp_mod.c
  43. +51 −0 bn_mp_mod_2d.c
  44. +47 −0 bn_mp_mod_d.c
  45. +80 −0 bn_mp_montgomery_reduce.c
  46. +53 −0 bn_mp_montgomery_setup.c
  47. +30 −0 bn_mp_mul.c
  48. +44 −0 bn_mp_mul_2.c
  49. +57 −0 bn_mp_mul_2d.c
  50. +46 −0 bn_mp_mul_d.c
  51. +36 −0 bn_mp_mulmod.c
  52. +115 −0 bn_mp_n_root.c
  53. +27 −0 bn_mp_neg.c
  54. +45 −0 bn_mp_or.c
  55. +48 −0 bn_mp_rand.c
  56. +28 −0 bn_mp_read_signed_bin.c
  57. +39 −0 bn_mp_read_unsigned_bin.c
  58. +95 −0 bn_mp_reduce.c
  59. +45 −0 bn_mp_rshd.c
  60. +24 −0 bn_mp_set.c
  61. +45 −0 bn_mp_set_int.c
  62. +28 −0 bn_mp_shrink.c
  63. +22 −0 bn_mp_signed_bin_size.c
  64. +29 −0 bn_mp_sqr.c
  65. +36 −0 bn_mp_sqrmod.c
  66. +58 −0 bn_mp_sub.c
  67. +33 −0 bn_mp_sub_d.c
  68. +36 −0 bn_mp_submod.c
  69. +28 −0 bn_mp_to_signed_bin.c
  70. +43 −0 bn_mp_to_unsigned_bin.c
  71. +23 −0 bn_mp_unsigned_bin_size.c
  72. +45 −0 bn_mp_xor.c
  73. +24 −0 bn_mp_zero.c
  74. +139 −0 bn_radix.c
  75. +33 −0 bn_reverse.c
  76. +91 −0 bn_s_mp_add.c
  77. +83 −0 bn_s_mp_mul_digs.c
  78. +77 −0 bn_s_mp_mul_high_digs.c
  79. +89 −0 bn_s_mp_sqr.c
  80. +74 −0 bn_s_mp_sub.c
  81. +19 −0 bncore.c
  82. +10 −0 changes.txt
  83. +19 −29 { → demo}/demo.c
  84. +15 −1 etc/makefile
  85. +116 −112 etc/mersenne.c
  86. +316 −250 etc/pprime.c
  87. +96 −0 etc/tune.c
  88. 0 { → etclib}/poly.c
  89. 0 { → etclib}/poly.h
  90. BIN ltmtest.exe
  91. +52 −12 makefile
  92. BIN mpitest.exe
  93. +242 −242 mtest/mtest.c
  94. +0 −34 timer.asm
  95. +39 −0 timings.txt
  96. +48 −7 bn.h → tommath.h
View
2 b.bat
@@ -1,2 +0,0 @@
-nasm -f elf timer.asm
-gcc -Wall -W -O3 -fomit-frame-pointer -funroll-loops -DTIMER_X86 demo.c bn.c timer.o -o ltmdemo
View
3,908 bn.c
0 additions, 3,908 deletions not shown because the diff is too large. Please use a local Git client to view these changes.
View
BIN bn.pdf
Binary file not shown.
View
69 bn.tex
@@ -1,7 +1,7 @@
\documentclass{article}
\begin{document}
-\title{LibTomMath v0.11 \\ A Free Multiple Precision Integer Library}
+\title{LibTomMath v0.12 \\ A Free Multiple Precision Integer Library}
\author{Tom St Denis \\ tomstdenis@iahu.ca}
\maketitle
\newpage
@@ -35,23 +35,24 @@ \section{Introduction}
\section{Building Against LibTomMath}
-Building against LibTomMath is very simple because there is only one source file. Simply add ``bn.c'' to your project and
-copy both ``bn.c'' and ``bn.h'' into your project directory. There is no configuration nor building required before hand.
+As of v0.12 LibTomMath is not a simple single source file project like MPI. LibTomMath retains the exact same API as MPI
+but is implemented differently. To build LibTomMath you will need a copy of GNU cc and GNU make. Both are free so if you
+don't have a copy don't whine to me about it.
-If you are porting an MPI application to LibTomMath the first step will be to remove all references to MPI and replace them
-with references to LibTomMath. For example, substitute
+To build the library type
\begin{verbatim}
-#include "mpi.h"
+make
\end{verbatim}
-with
+This will build the library file libtommath.a. If you want to build the library and also install it (in /usr/bin and /usr/include) then
+type
\begin{verbatim}
-#include "bn.h"
+make install
\end{verbatim}
-Remove ``mpi.c'' from your project and replace it with ``bn.c''.
+Now within your application include ``tommath.h'' and link against libtommath.a to get MPI-like functionality.
\section{Programming with LibTomMath}
@@ -183,6 +184,28 @@ \subsection{Digit Manipulations}
/* c = a mod 2^d */
int mp_mod_2d(mp_int *a, int b, mp_int *c);
+
+/* computes a = 2^b */
+int mp_2expt(mp_int *a, int b);
+
+/* makes a pseudo-random int of a given size */
+int mp_rand(mp_int *a, int digits);
+
+\end{verbatim}
+
+\subsection{Binary Operations}
+
+\begin{verbatim}
+
+/* c = a XOR b */
+int mp_xor(mp_int *a, mp_int *b, mp_int *c);
+
+/* c = a OR b */
+int mp_or(mp_int *a, mp_int *b, mp_int *c);
+
+/* c = a AND b */
+int mp_and(mp_int *a, mp_int *b, mp_int *c);
+
\end{verbatim}
\subsection{Basic Arithmetic}
@@ -366,6 +389,27 @@ \subsubsection{mp\_mod\_2d(mp\_int *a, int b, mp\_int *c)}
This function requires $O(N)$ additional digits of memory and $O(2 \cdot N)$ time.
+\subsubsection{mp\_2expt(mp\_int *a, int b)}
+Computes $a = 2^b$ by first setting $a$ to zero then OR'ing the correct bit to get the right value.
+
+\subsubsection{mp\_rand(mp\_int *a, int digits)}
+Computes a pseudo-random (\textit{via rand()}) integer that is always ``$digits$'' digits in length. Not for
+cryptographic use.
+
+\subsection{Binary Arithmetic}
+\subsubsection{mp\_xor(mp\_int *a, mp\_int *b, mp\_int *c)}
+Computes $c = a \oplus b$, pseudo-extends with zeroes whichever of $a$ or $b$ is shorter such that the length
+of $c$ is the maximum length of the two inputs.
+
+\subsubsection{mp\_or(mp\_int *a, mp\_int *b, mp\_int *c)}
+Computes $c = a \lor b$, pseudo-extends with zeroes whichever of $a$ or $b$ is shorter such that the length
+of $c$ is the maximum length of the two inputs.
+
+\subsubsection{mp\_and(mp\_int *a, mp\_int *b, mp\_int *c)}
+Computes $c = a \land b$, pseudo-extends with zeroes whichever of $a$ or $b$ is shorter such that the length
+of $c$ is the maximum length of the two inputs.
+
+
\subsection{Basic Arithmetic}
\subsubsection{mp\_cmp(mp\_int *a, mp\_int *b)}
@@ -622,14 +666,13 @@ \subsection{Exponentiation Algorithms}
124 multiplications. The MPI method would have 512 squarings and 512 multiplications. Randomly every $2k$ bits another
multiplication is saved via the sliding-window technique on top of the savings the $k$-ary method provides.
-Both LibTomMath and MPI use Barrett reduction instead of division to reduce the numbers modulo the modulus given.
+Both LibTomMath and MPI use Barrett reduction instead of division to reduce the numbers modulo the modulus given.
However, LibTomMath can take advantage of the fact that the multiplications required within the Barrett reduction
do not have to give full precision. As a result the reduction step is much faster and just as accurate. The LibTomMath code
will automatically determine at run-time (e.g. when its called) whether the faster multiplier can be used. The
faster multipliers have also been optimized into the two variants (baseline and comba baseline).
-As a result of all these changes exponentiation in LibTomMath is much faster than compared to MPI.
-
-
+LibTomMath also has a variant of the exptmod function that uses Montgomery reductions instead of Barrett reductions
+which is faser. As a result of all these changes exponentiation in LibTomMath is much faster than compared to MPI.
\end{document}
View
168 bn_fast_mp_invmod.c
@@ -0,0 +1,168 @@
+/* LibTomMath, multiple-precision integer library -- Tom St Denis
+ *
+ * LibTomMath is library that provides for multiple-precision
+ * integer arithmetic as well as number theoretic functionality.
+ *
+ * The library is 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://libtommath.iahu.ca
+ */
+#include <tommath.h>
+
+/* computes the modular inverse via binary extended euclidean algorithm, that is c = 1/a mod b */
+int
+fast_mp_invmod (mp_int * a, mp_int * b, mp_int * c)
+{
+ mp_int x, y, u, v, B, D;
+ int res, neg;
+
+
+ if ((res = mp_init (&x)) != MP_OKAY) {
+ goto __ERR;
+ }
+
+ if ((res = mp_init (&y)) != MP_OKAY) {
+ goto __X;
+ }
+
+ if ((res = mp_init (&u)) != MP_OKAY) {
+ goto __Y;
+ }
+
+ if ((res = mp_init (&v)) != MP_OKAY) {
+ goto __U;
+ }
+
+ if ((res = mp_init (&B)) != MP_OKAY) {
+ goto __V;
+ }
+
+ if ((res = mp_init (&D)) != MP_OKAY) {
+ goto __B;
+ }
+
+ /* x == modulus, y == value to invert */
+ if ((res = mp_copy (b, &x)) != MP_OKAY) {
+ goto __D;
+ }
+ if ((res = mp_copy (a, &y)) != MP_OKAY) {
+ goto __D;
+ }
+
+ if ((res = mp_abs (&y, &y)) != MP_OKAY) {
+ goto __D;
+ }
+
+ /* 2. [modified] if x,y are both even then return an error! */
+ if (mp_iseven (&x) == 1 && mp_iseven (&y) == 1) {
+ res = MP_VAL;
+ goto __D;
+ }
+
+ /* 3. u=x, v=y, A=1, B=0, C=0,D=1 */
+ if ((res = mp_copy (&x, &u)) != MP_OKAY) {
+ goto __D;
+ }
+ if ((res = mp_copy (&y, &v)) != MP_OKAY) {
+ goto __D;
+ }
+ mp_set (&D, 1);
+
+
+top:
+ /* 4. while u is even do */
+ while (mp_iseven (&u) == 1) {
+ /* 4.1 u = u/2 */
+ if ((res = mp_div_2 (&u, &u)) != MP_OKAY) {
+ goto __D;
+ }
+ /* 4.2 if A or B is odd then */
+ if (mp_iseven (&B) == 0) {
+ if ((res = mp_sub (&B, &x, &B)) != MP_OKAY) {
+ goto __D;
+ }
+ }
+ /* A = A/2, B = B/2 */
+ if ((res = mp_div_2 (&B, &B)) != MP_OKAY) {
+ goto __D;
+ }
+ }
+
+
+ /* 5. while v is even do */
+ while (mp_iseven (&v) == 1) {
+ /* 5.1 v = v/2 */
+ if ((res = mp_div_2 (&v, &v)) != MP_OKAY) {
+ goto __D;
+ }
+ /* 5.2 if C,D are even then */
+ if (mp_iseven (&D) == 0) {
+ /* C = (C+y)/2, D = (D-x)/2 */
+ if ((res = mp_sub (&D, &x, &D)) != MP_OKAY) {
+ goto __D;
+ }
+ }
+ /* C = C/2, D = D/2 */
+ if ((res = mp_div_2 (&D, &D)) != MP_OKAY) {
+ goto __D;
+ }
+ }
+
+ /* 6. if u >= v then */
+ if (mp_cmp (&u, &v) != MP_LT) {
+ /* u = u - v, A = A - C, B = B - D */
+ if ((res = mp_sub (&u, &v, &u)) != MP_OKAY) {
+ goto __D;
+ }
+
+ if ((res = mp_sub (&B, &D, &B)) != MP_OKAY) {
+ goto __D;
+ }
+ } else {
+ /* v - v - u, C = C - A, D = D - B */
+ if ((res = mp_sub (&v, &u, &v)) != MP_OKAY) {
+ goto __D;
+ }
+
+ if ((res = mp_sub (&D, &B, &D)) != MP_OKAY) {
+ goto __D;
+ }
+ }
+
+ /* if not zero goto step 4 */
+ if (mp_iszero (&u) == 0)
+ goto top;
+
+ /* now a = C, b = D, gcd == g*v */
+
+ /* if v != 1 then there is no inverse */
+ if (mp_cmp_d (&v, 1) != MP_EQ) {
+ res = MP_VAL;
+ goto __D;
+ }
+
+ /* b is now the inverse */
+ neg = a->sign;
+ while (D.sign == MP_NEG) {
+ if ((res = mp_add (&D, b, &D)) != MP_OKAY) {
+ goto __D;
+ }
+ }
+ mp_exch (&D, c);
+ c->sign = neg;
+ res = MP_OKAY;
+
+__D:mp_clear (&D);
+__B:mp_clear (&B);
+__V:mp_clear (&v);
+__U:mp_clear (&u);
+__Y:mp_clear (&y);
+__X:mp_clear (&x);
+__ERR:
+ return res;
+}
View
116 bn_fast_mp_montgomery_reduce.c
@@ -0,0 +1,116 @@
+/* LibTomMath, multiple-precision integer library -- Tom St Denis
+ *
+ * LibTomMath is library that provides for multiple-precision
+ * integer arithmetic as well as number theoretic functionality.
+ *
+ * The library is 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://libtommath.iahu.ca
+ */
+#include <tommath.h>
+
+/* computes xR^-1 == x (mod N) via Montgomery Reduction (comba) */
+int
+fast_mp_montgomery_reduce (mp_int * a, mp_int * m, mp_digit mp)
+{
+ int ix, res, olduse;
+ mp_word W[512];
+
+ /* get old used count */
+ olduse = a->used;
+
+ /* grow a as required */
+ if (a->alloc < m->used + 1) {
+ if ((res = mp_grow (a, m->used + 1)) != MP_OKAY) {
+ return res;
+ }
+ }
+
+ /* copy the digits of a */
+ for (ix = 0; ix < a->used; ix++) {
+ W[ix] = a->dp[ix];
+ }
+
+ /* zero the high words */
+ for (; ix < m->used * 2 + 1; ix++) {
+ W[ix] = 0;
+ }
+
+ for (ix = 0; ix < m->used; ix++) {
+ /* ui = ai * m' mod b
+ *
+ * We avoid a double precision multiplication (which isn't required)
+ * by casting the value down to a mp_digit. Note this requires that W[ix-1] have
+ * the carry cleared (see after the inner loop)
+ */
+ register mp_digit ui;
+ ui = (((mp_digit) (W[ix] & MP_MASK)) * mp) & MP_MASK;
+
+ /* a = a + ui * m * b^i
+ *
+ * This is computed in place and on the fly. The multiplication
+ * by b^i is handled by offseting which columns the results
+ * are added to.
+ *
+ * Note the comba method normally doesn't handle carries in the inner loop
+ * In this case we fix the carry from the previous column since the Montgomery
+ * reduction requires digits of the result (so far) [see above] to work. This is
+ * handled by fixing up one carry after the inner loop. The carry fixups are done
+ * in order so after these loops the first m->used words of W[] have the carries
+ * fixed
+ */
+ {
+ register int iy;
+ register mp_digit *tmpx;
+ register mp_word *_W;
+
+ /* aliases */
+ tmpx = m->dp;
+ _W = W + ix;
+
+ /* inner loop */
+ for (iy = 0; iy < m->used; iy++) {
+ *_W++ += ((mp_word) ui) * ((mp_word) * tmpx++);
+ }
+ }
+
+ /* now fix carry for next digit, W[ix+1] */
+ W[ix + 1] += W[ix] >> ((mp_word) DIGIT_BIT);
+ }
+
+ /* nox fix rest of carries */
+ for (++ix; ix <= m->used * 2 + 1; ix++) {
+ W[ix] += (W[ix - 1] >> ((mp_word) DIGIT_BIT));
+ }
+
+ /* copy out, A = A/b^n
+ *
+ * The result is A/b^n but instead of converting from an array of mp_word
+ * to mp_digit than calling mp_rshd we just copy them in the right
+ * order
+ */
+ for (ix = 0; ix < m->used + 1; ix++) {
+ a->dp[ix] = W[ix + m->used] & ((mp_word) MP_MASK);
+ }
+
+ /* set the max used */
+ a->used = m->used + 1;
+
+ /* zero oldused digits, if the input a was larger than
+ * m->used+1 we'll have to clear the digits */
+ for (; ix < olduse; ix++) {
+ a->dp[ix] = 0;
+ }
+ mp_clamp (a);
+
+ /* if A >= m then A = A - m */
+ if (mp_cmp_mag (a, m) != MP_LT) {
+ return s_mp_sub (a, m, a);
+ }
+ return MP_OKAY;
+}
View
113 bn_fast_s_mp_mul_digs.c
@@ -0,0 +1,113 @@
+/* LibTomMath, multiple-precision integer library -- Tom St Denis
+ *
+ * LibTomMath is library that provides for multiple-precision
+ * integer arithmetic as well as number theoretic functionality.
+ *
+ * The library is 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://libtommath.iahu.ca
+ */
+#include <tommath.h>
+
+/* Fast (comba) multiplier
+ *
+ * This is the fast column-array [comba] multiplier. It is designed to compute
+ * the columns of the product first then handle the carries afterwards. This
+ * has the effect of making the nested loops that compute the columns very
+ * simple and schedulable on super-scalar processors.
+ *
+ */
+int
+fast_s_mp_mul_digs (mp_int * a, mp_int * b, mp_int * c, int digs)
+{
+ int olduse, res, pa, ix;
+ mp_word W[512];
+
+ if (c->alloc < digs) {
+ if ((res = mp_grow (c, digs)) != MP_OKAY) {
+ return res;
+ }
+ }
+
+ /* clear temp buf (the columns) */
+ memset (W, 0, sizeof (mp_word) * digs);
+
+ /* calculate the columns */
+ pa = a->used;
+ for (ix = 0; ix < pa; ix++) {
+
+ /* this multiplier has been modified to allow you to control how many digits
+ * of output are produced. So at most we want to make upto "digs" digits
+ * of output
+ */
+
+
+ /* this adds products to distinct columns (at ix+iy) of W
+ * note that each step through the loop is not dependent on
+ * the previous which means the compiler can easily unroll
+ * the loop without scheduling problems
+ */
+ {
+ register mp_digit tmpx, *tmpy;
+ register mp_word *_W;
+ register int iy, pb;
+
+ /* alias for the the word on the left e.g. A[ix] * A[iy] */
+ tmpx = a->dp[ix];
+
+ /* alias for the right side */
+ tmpy = b->dp;
+
+ /* alias for the columns, each step through the loop adds a new
+ term to each column
+ */
+ _W = W + ix;
+
+ /* the number of digits is limited by their placement. E.g.
+ we avoid multiplying digits that will end up above the # of
+ digits of precision requested
+ */
+ pb = MIN (b->used, digs - ix);
+
+ for (iy = 0; iy < pb; iy++) {
+ *_W++ += ((mp_word) tmpx) * ((mp_word) * tmpy++);
+ }
+ }
+
+ }
+
+ /* setup dest */
+ olduse = c->used;
+ c->used = digs;
+
+
+ /* At this point W[] contains the sums of each column. To get the
+ * correct result we must take the extra bits from each column and
+ * carry them down
+ *
+ * Note that while this adds extra code to the multiplier it saves time
+ * since the carry propagation is removed from the above nested loop.
+ * This has the effect of reducing the work from N*(N+N*c)==N^2 + c*N^2 to
+ * N^2 + N*c where c is the cost of the shifting. On very small numbers
+ * this is slower but on most cryptographic size numbers it is faster.
+ */
+
+ for (ix = 1; ix < digs; ix++) {
+ W[ix] += (W[ix - 1] >> ((mp_word) DIGIT_BIT));
+ c->dp[ix - 1] = (mp_digit) (W[ix - 1] & ((mp_word) MP_MASK));
+ }
+ c->dp[digs - 1] = (mp_digit) (W[digs - 1] & ((mp_word) MP_MASK));
+
+ /* clear unused */
+ for (; ix < olduse; ix++) {
+ c->dp[ix] = 0;
+ }
+
+ mp_clamp (c);
+ return MP_OKAY;
+}
View
86 bn_fast_s_mp_mul_high_digs.c
@@ -0,0 +1,86 @@
+/* LibTomMath, multiple-precision integer library -- Tom St Denis
+ *
+ * LibTomMath is library that provides for multiple-precision
+ * integer arithmetic as well as number theoretic functionality.
+ *
+ * The library is 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://libtommath.iahu.ca
+ */
+#include <tommath.h>
+
+/* this is a modified version of fast_s_mp_mul_digs that only produces
+ * output digits *above* digs. See the comments for fast_s_mp_mul_digs
+ * to see how it works.
+ *
+ * This is used in the Barrett reduction since for one of the multiplications
+ * only the higher digits were needed. This essentially halves the work.
+ */
+int
+fast_s_mp_mul_high_digs (mp_int * a, mp_int * b, mp_int * c, int digs)
+{
+ int oldused, newused, res, pa, pb, ix;
+ mp_word W[512];
+
+
+ newused = a->used + b->used + 1;
+ if (c->alloc < newused) {
+ if ((res = mp_grow (c, newused)) != MP_OKAY) {
+ return res;
+ }
+ }
+
+ /* like the other comba method we compute the columns first */
+ pa = a->used;
+ pb = b->used;
+ memset (&W[digs], 0, (pa + pb + 1 - digs) * sizeof (mp_word));
+ for (ix = 0; ix < pa; ix++) {
+ {
+ register mp_digit tmpx, *tmpy;
+ register int iy;
+ register mp_word *_W;
+
+ /* work todo, that is we only calculate digits that are at "digs" or above */
+ iy = digs - ix;
+
+ /* copy of word on the left of A[ix] * B[iy] */
+ tmpx = a->dp[ix];
+
+ /* alias for right side */
+ tmpy = b->dp + iy;
+
+ /* alias for the columns of output. Offset to be equal to or above the
+ * smallest digit place requested
+ */
+ _W = &(W[digs]);
+
+ /* compute column products for digits above the minimum */
+ for (; iy < pb; iy++) {
+ *_W++ += ((mp_word) tmpx) * ((mp_word) * tmpy++);
+ }
+ }
+ }
+
+ /* setup dest */
+ oldused = c->used;
+ c->used = newused;
+
+ /* now convert the array W downto what we need */
+ for (ix = digs + 1; ix < newused; ix++) {
+ W[ix] += (W[ix - 1] >> ((mp_word) DIGIT_BIT));
+ c->dp[ix - 1] = (mp_digit) (W[ix - 1] & ((mp_word) MP_MASK));
+ }
+ c->dp[(pa + pb + 1) - 1] =
+ (mp_digit) (W[(pa + pb + 1) - 1] & ((mp_word) MP_MASK));
+
+ for (; ix < oldused; ix++) {
+ c->dp[ix] = 0;
+ }
+ mp_clamp (c);
+ return MP_OKAY;
+}
View
112 bn_fast_s_mp_sqr.c
@@ -0,0 +1,112 @@
+/* LibTomMath, multiple-precision integer library -- Tom St Denis
+ *
+ * LibTomMath is library that provides for multiple-precision
+ * integer arithmetic as well as number theoretic functionality.
+ *
+ * The library is 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://libtommath.iahu.ca
+ */
+#include <tommath.h>
+
+/* fast squaring
+ *
+ * This is the comba method where the columns of the product are computed first
+ * then the carries are computed. This has the effect of making a very simple
+ * inner loop that is executed the most
+ *
+ * W2 represents the outer products and W the inner.
+ *
+ * A further optimizations is made because the inner products are of the form
+ * "A * B * 2". The *2 part does not need to be computed until the end which is
+ * good because 64-bit shifts are slow!
+ *
+ *
+ */
+int
+fast_s_mp_sqr (mp_int * a, mp_int * b)
+{
+ int olduse, newused, res, ix, pa;
+ mp_word W2[512], W[512];
+
+ pa = a->used;
+ newused = pa + pa + 1;
+ if (b->alloc < newused) {
+ if ((res = mp_grow (b, newused)) != MP_OKAY) {
+ return res;
+ }
+ }
+
+ /* zero temp buffer (columns)
+ * Note that there are two buffers. Since squaring requires
+ * a outter and inner product and the inner product requires
+ * computing a product and doubling it (a relatively expensive
+ * op to perform n^2 times if you don't have to) the inner and
+ * outer products are computed in different buffers. This way
+ * the inner product can be doubled using n doublings instead of
+ * n^2
+ */
+ memset (W, 0, newused * sizeof (mp_word));
+ memset (W2, 0, newused * sizeof (mp_word));
+
+ /* This computes the inner product. To simplify the inner N^2 loop
+ * the multiplication by two is done afterwards in the N loop.
+ */
+ for (ix = 0; ix < pa; ix++) {
+ /* compute the outer product */
+ W2[ix + ix] += ((mp_word) a->dp[ix]) * ((mp_word) a->dp[ix]);
+
+ {
+ register mp_digit tmpx, *tmpy;
+ register mp_word *_W;
+ register int iy;
+
+ /* copy of left side */
+ tmpx = a->dp[ix];
+
+ /* alias for right side */
+ tmpy = a->dp + (ix + 1);
+
+ /* the column to store the result in */
+ _W = W + (ix + ix + 1);
+
+ /* inner products */
+ for (iy = ix + 1; iy < pa; iy++) {
+ *_W++ += ((mp_word) tmpx) * ((mp_word) * tmpy++);
+ }
+ }
+ }
+
+ /* setup dest */
+ olduse = b->used;
+ b->used = newused;
+
+ /* double first value, since the inner products are half of what they should be */
+ W[0] += W[0] + W2[0];
+
+ /* now compute digits */
+ for (ix = 1; ix < newused; ix++) {
+ /* double/add next digit */
+ W[ix] += W[ix] + W2[ix];
+
+ W[ix] = W[ix] + (W[ix - 1] >> ((mp_word) DIGIT_BIT));
+ b->dp[ix - 1] = (mp_digit) (W[ix - 1] & ((mp_word) MP_MASK));
+ }
+ b->dp[(newused) - 1] = (mp_digit) (W[(newused) - 1] & ((mp_word) MP_MASK));
+
+ /* clear high */
+ for (; ix < olduse; ix++) {
+ b->dp[ix] = 0;
+ }
+
+ /* fix the sign (since we no longer make a fresh temp) */
+ b->sign = MP_ZPOS;
+
+ mp_clamp (b);
+ return MP_OKAY;
+}
View
31 bn_mp_2expt.c
@@ -0,0 +1,31 @@
+/* LibTomMath, multiple-precision integer library -- Tom St Denis
+ *
+ * LibTomMath is library that provides for multiple-precision
+ * integer arithmetic as well as number theoretic functionality.
+ *
+ * The library is 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://libtommath.iahu.ca
+ */
+#include <tommath.h>
+
+/* computes a = 2^b */
+int
+mp_2expt (mp_int * a, int b)
+{
+ int res;
+
+ mp_zero (a);
+ if ((res = mp_grow (a, b / DIGIT_BIT + 1)) != MP_OKAY) {
+ return res;
+ }
+ a->used = b / DIGIT_BIT + 1;
+ a->dp[b / DIGIT_BIT] = 1 << (b % DIGIT_BIT);
+
+ return MP_OKAY;
+}
View
27 bn_mp_abs.c
@@ -0,0 +1,27 @@
+/* LibTomMath, multiple-precision integer library -- Tom St Denis
+ *
+ * LibTomMath is library that provides for multiple-precision
+ * integer arithmetic as well as number theoretic functionality.
+ *
+ * The library is 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://libtommath.iahu.ca
+ */
+#include <tommath.h>
+
+/* b = |a| */
+int
+mp_abs (mp_int * a, mp_int * b)
+{
+ int res;
+ if ((res = mp_copy (a, b)) != MP_OKAY) {
+ return res;
+ }
+ b->sign = MP_ZPOS;
+ return MP_OKAY;
+}
View
56 bn_mp_add.c
@@ -0,0 +1,56 @@
+/* LibTomMath, multiple-precision integer library -- Tom St Denis
+ *
+ * LibTomMath is library that provides for multiple-precision
+ * integer arithmetic as well as number theoretic functionality.
+ *
+ * The library is 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://libtommath.iahu.ca
+ */
+#include <tommath.h>
+
+/* high level addition (handles signs) */
+int
+mp_add (mp_int * a, mp_int * b, mp_int * c)
+{
+ int sa, sb, res;
+
+
+ sa = a->sign;
+ sb = b->sign;
+
+ /* handle four cases */
+ if (sa == MP_ZPOS && sb == MP_ZPOS) {
+ /* both positive */
+ res = s_mp_add (a, b, c);
+ c->sign = MP_ZPOS;
+ } else if (sa == MP_ZPOS && sb == MP_NEG) {
+ /* a + -b == a - b, but if b>a then we do it as -(b-a) */
+ if (mp_cmp_mag (a, b) == MP_LT) {
+ res = s_mp_sub (b, a, c);
+ c->sign = MP_NEG;
+ } else {
+ res = s_mp_sub (a, b, c);
+ c->sign = MP_ZPOS;
+ }
+ } else if (sa == MP_NEG && sb == MP_ZPOS) {
+ /* -a + b == b - a, but if a>b then we do it as -(a-b) */
+ if (mp_cmp_mag (a, b) == MP_GT) {
+ res = s_mp_sub (a, b, c);
+ c->sign = MP_NEG;
+ } else {
+ res = s_mp_sub (b, a, c);
+ c->sign = MP_ZPOS;
+ }
+ } else {
+ /* -a + -b == -(a + b) */
+ res = s_mp_add (a, b, c);
+ c->sign = MP_NEG;
+ }
+ return res;
+}
View
33 bn_mp_add_d.c
@@ -0,0 +1,33 @@
+/* LibTomMath, multiple-precision integer library -- Tom St Denis
+ *
+ * LibTomMath is library that provides for multiple-precision
+ * integer arithmetic as well as number theoretic functionality.
+ *
+ * The library is 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://libtommath.iahu.ca
+ */
+#include <tommath.h>
+
+/* single digit addition */
+int
+mp_add_d (mp_int * a, mp_digit b, mp_int * c)
+{
+ mp_int t;
+ int res;
+
+
+ if ((res = mp_init (&t)) != MP_OKAY) {
+ return res;
+ }
+ mp_set (&t, b);
+ res = mp_add (a, &t, c);
+
+ mp_clear (&t);
+ return res;
+}
View
36 bn_mp_addmod.c
@@ -0,0 +1,36 @@
+/* LibTomMath, multiple-precision integer library -- Tom St Denis
+ *
+ * LibTomMath is library that provides for multiple-precision
+ * integer arithmetic as well as number theoretic functionality.
+ *
+ * The library is 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://libtommath.iahu.ca
+ */
+#include <tommath.h>
+
+/* d = a + b (mod c) */
+int
+mp_addmod (mp_int * a, mp_int * b, mp_int * c, mp_int * d)
+{
+ int res;
+ mp_int t;
+
+
+ if ((res = mp_init (&t)) != MP_OKAY) {
+ return res;
+ }
+
+ if ((res = mp_add (a, b, &t)) != MP_OKAY) {
+ mp_clear (&t);
+ return res;
+ }
+ res = mp_mod (&t, c, d);
+ mp_clear (&t);
+ return res;
+}
View
51 bn_mp_and.c
@@ -0,0 +1,51 @@
+/* LibTomMath, multiple-precision integer library -- Tom St Denis
+ *
+ * LibTomMath is library that provides for multiple-precision
+ * integer arithmetic as well as number theoretic functionality.
+ *
+ * The library is 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://libtommath.iahu.ca
+ */
+#include <tommath.h>
+
+/* AND two ints together */
+int
+mp_and (mp_int * a, mp_int * b, mp_int * c)
+{
+ int res, ix, px;
+ mp_int t, *x;
+
+ if (a->used > b->used) {
+ if ((res = mp_init_copy (&t, a)) != MP_OKAY) {
+ return res;
+ }
+ px = b->used;
+ x = b;
+ } else {
+ if ((res = mp_init_copy (&t, b)) != MP_OKAY) {
+ return res;
+ }
+ px = a->used;
+ x = a;
+ }
+
+ for (ix = 0; ix < px; ix++) {
+ t.dp[ix] &= x->dp[ix];
+ }
+
+ /* zero digits above the last from the smallest mp_int */
+ for (; ix < t.used; ix++) {
+ t.dp[ix] = 0;
+ }
+
+ mp_clamp (&t);
+ mp_exch (c, &t);
+ mp_clear (&t);
+ return MP_OKAY;
+}
View
26 bn_mp_clamp.c
@@ -0,0 +1,26 @@
+/* LibTomMath, multiple-precision integer library -- Tom St Denis
+ *
+ * LibTomMath is library that provides for multiple-precision
+ * integer arithmetic as well as number theoretic functionality.
+ *
+ * The library is 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://libtommath.iahu.ca
+ */
+#include <tommath.h>
+
+/* trim unused digits */
+void
+mp_clamp (mp_int * a)
+{
+ while (a->used > 0 && a->dp[a->used - 1] == 0)
+ --(a->used);
+ if (a->used == 0) {
+ a->sign = MP_ZPOS;
+ }
+}
View
33 bn_mp_clear.c
@@ -0,0 +1,33 @@
+/* LibTomMath, multiple-precision integer library -- Tom St Denis
+ *
+ * LibTomMath is library that provides for multiple-precision
+ * integer arithmetic as well as number theoretic functionality.
+ *
+ * The library is 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://libtommath.iahu.ca
+ */
+#include <tommath.h>
+
+/* clear one (frees) */
+void
+mp_clear (mp_int * a)
+{
+ if (a->dp != NULL) {
+
+ /* first zero the digits */
+ memset (a->dp, 0, sizeof (mp_digit) * a->used);
+
+ /* free ram */
+ free (a->dp);
+
+ /* reset members to make debugging easier */
+ a->dp = NULL;
+ a->alloc = a->used = 0;
+ }
+}
View
30 bn_mp_cmp.c
@@ -0,0 +1,30 @@
+/* LibTomMath, multiple-precision integer library -- Tom St Denis
+ *
+ * LibTomMath is library that provides for multiple-precision
+ * integer arithmetic as well as number theoretic functionality.
+ *
+ * The library is 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://libtommath.iahu.ca
+ */
+#include <tommath.h>
+
+/* compare two ints (signed)*/
+int
+mp_cmp (mp_int * a, mp_int * b)
+{
+ int res;
+ /* compare based on sign */
+ if (a->sign == MP_NEG && b->sign == MP_ZPOS) {
+ return MP_LT;
+ } else if (a->sign == MP_ZPOS && b->sign == MP_NEG) {
+ return MP_GT;
+ }
+ res = mp_cmp_mag (a, b);
+ return res;
+}
View
37 bn_mp_cmp_d.c
@@ -0,0 +1,37 @@
+/* LibTomMath, multiple-precision integer library -- Tom St Denis
+ *
+ * LibTomMath is library that provides for multiple-precision
+ * integer arithmetic as well as number theoretic functionality.
+ *
+ * The library is 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://libtommath.iahu.ca
+ */
+#include <tommath.h>
+
+/* compare a digit */
+int
+mp_cmp_d (mp_int * a, mp_digit b)
+{
+
+ if (a->sign == MP_NEG) {
+ return MP_LT;
+ }
+
+ if (a->used > 1) {
+ return MP_GT;
+ }
+
+ if (a->dp[0] > b) {
+ return MP_GT;
+ } else if (a->dp[0] < b) {
+ return MP_LT;
+ } else {
+ return MP_EQ;
+ }
+}
View
40 bn_mp_cmp_mag.c
@@ -0,0 +1,40 @@
+/* LibTomMath, multiple-precision integer library -- Tom St Denis
+ *
+ * LibTomMath is library that provides for multiple-precision
+ * integer arithmetic as well as number theoretic functionality.
+ *
+ * The library is 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://libtommath.iahu.ca
+ */
+#include <tommath.h>
+
+/* compare maginitude of two ints (unsigned) */
+int
+mp_cmp_mag (mp_int * a, mp_int * b)
+{
+ int n;
+
+
+ /* compare based on # of non-zero digits */
+ if (a->used > b->used) {
+ return MP_GT;
+ } else if (a->used < b->used) {
+ return MP_LT;
+ }
+
+ /* compare based on digits */
+ for (n = a->used - 1; n >= 0; n--) {
+ if (a->dp[n] > b->dp[n]) {
+ return MP_GT;
+ } else if (a->dp[n] < b->dp[n]) {
+ return MP_LT;
+ }
+ }
+ return MP_EQ;
+}
View
48 bn_mp_copy.c
@@ -0,0 +1,48 @@
+/* LibTomMath, multiple-precision integer library -- Tom St Denis
+ *
+ * LibTomMath is library that provides for multiple-precision
+ * integer arithmetic as well as number theoretic functionality.
+ *
+ * The library is 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://libtommath.iahu.ca
+ */
+#include <tommath.h>
+
+/* copy, b = a */
+int
+mp_copy (mp_int * a, mp_int * b)
+{
+ int res, n;
+
+
+ /* if dst == src do nothing */
+ if (a == b || a->dp == b->dp) {
+ return MP_OKAY;
+ }
+
+ /* grow dest */
+ if ((res = mp_grow (b, a->used)) != MP_OKAY) {
+ return res;
+ }
+
+ /* zero b and copy the parameters over */
+ b->used = a->used;
+ b->sign = a->sign;
+
+ /* copy all the digits */
+ for (n = 0; n < a->used; n++) {
+ b->dp[n] = a->dp[n];
+ }
+
+ /* clear high digits */
+ for (n = b->used; n < b->alloc; n++) {
+ b->dp[n] = 0;
+ }
+ return MP_OKAY;
+}
View
35 bn_mp_count_bits.c
@@ -0,0 +1,35 @@
+/* LibTomMath, multiple-precision integer library -- Tom St Denis
+ *
+ * LibTomMath is library that provides for multiple-precision
+ * integer arithmetic as well as number theoretic functionality.
+ *
+ * The library is 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://libtommath.iahu.ca
+ */
+#include <tommath.h>
+
+/* returns the number of bits in an int */
+int
+mp_count_bits (mp_int * a)
+{
+ int r;
+ mp_digit q;
+
+ if (a->used == 0) {
+ return 0;
+ }
+
+ r = (a->used - 1) * DIGIT_BIT;
+ q = a->dp[a->used - 1];
+ while (q > ((mp_digit) 0)) {
+ ++r;
+ q >>= ((mp_digit) 1);
+ }
+ return r;
+}
View
200 bn_mp_div.c
@@ -0,0 +1,200 @@
+/* LibTomMath, multiple-precision integer library -- Tom St Denis
+ *
+ * LibTomMath is library that provides for multiple-precision
+ * integer arithmetic as well as number theoretic functionality.
+ *
+ * The library is 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://libtommath.iahu.ca
+ */
+#include <tommath.h>
+
+/* integer signed division. c*b + d == a [e.g. a/b, c=quotient, d=remainder]
+ * HAC pp.598 Algorithm 14.20
+ *
+ * Note that the description in HAC is horribly incomplete. For example,
+ * it doesn't consider the case where digits are removed from 'x' in the inner
+ * loop. It also doesn't consider the case that y has fewer than three digits, etc..
+ *
+ * The overall algorithm is as described as 14.20 from HAC but fixed to treat these cases.
+*/
+int
+mp_div (mp_int * a, mp_int * b, mp_int * c, mp_int * d)
+{
+ mp_int q, x, y, t1, t2;
+ int res, n, t, i, norm, neg;
+
+
+ /* is divisor zero ? */
+ if (mp_iszero (b) == 1) {
+ return MP_VAL;
+ }
+
+ /* if a < b then q=0, r = a */
+ if (mp_cmp_mag (a, b) == MP_LT) {
+ if (d != NULL) {
+ res = mp_copy (a, d);
+ } else {
+ res = MP_OKAY;
+ }
+ if (c != NULL) {
+ mp_zero (c);
+ }
+ return res;
+ }
+
+ if ((res = mp_init_size (&q, a->used + 2)) != MP_OKAY) {
+ return res;
+ }
+ q.used = a->used + 2;
+
+ if ((res = mp_init (&t1)) != MP_OKAY) {
+ goto __Q;
+ }
+
+ if ((res = mp_init (&t2)) != MP_OKAY) {
+ goto __T1;
+ }
+
+ if ((res = mp_init_copy (&x, a)) != MP_OKAY) {
+ goto __T2;
+ }
+
+ if ((res = mp_init_copy (&y, b)) != MP_OKAY) {
+ goto __X;
+ }
+
+ /* fix the sign */
+ neg = (a->sign == b->sign) ? MP_ZPOS : MP_NEG;
+ x.sign = y.sign = MP_ZPOS;
+
+ /* normalize both x and y, ensure that y >= b/2, [b == 2^DIGIT_BIT] */
+ norm = 0;
+ while ((y.dp[y.used - 1] & (((mp_digit) 1) << (DIGIT_BIT - 1))) ==
+ ((mp_digit) 0)) {
+ ++norm;
+ if ((res = mp_mul_2d (&x, 1, &x)) != MP_OKAY) {
+ goto __Y;
+ }
+ if ((res = mp_mul_2d (&y, 1, &y)) != MP_OKAY) {
+ goto __Y;
+ }
+ }
+
+ /* note hac does 0 based, so if used==5 then its 0,1,2,3,4, e.g. use 4 */
+ n = x.used - 1;
+ t = y.used - 1;
+
+ /* step 2. while (x >= y*b^n-t) do { q[n-t] += 1; x -= y*b^{n-t} } */
+ if ((res = mp_lshd (&y, n - t)) != MP_OKAY) { /* y = y*b^{n-t} */
+ goto __Y;
+ }
+
+ while (mp_cmp (&x, &y) != MP_LT) {
+ ++(q.dp[n - t]);
+ if ((res = mp_sub (&x, &y, &x)) != MP_OKAY) {
+ goto __Y;
+ }
+ }
+
+ /* reset y by shifting it back down */
+ mp_rshd (&y, n - t);
+
+ /* step 3. for i from n down to (t + 1) */
+ for (i = n; i >= (t + 1); i--) {
+ if (i > x.alloc)
+ continue;
+
+ /* step 3.1 if xi == yt then set q{i-t-1} to b-1, otherwise set q{i-t-1} to (xi*b + x{i-1})/yt */
+ if (x.dp[i] == y.dp[t]) {
+ q.dp[i - t - 1] = ((1UL << DIGIT_BIT) - 1UL);
+ } else {
+ mp_word tmp;
+ tmp = ((mp_word) x.dp[i]) << ((mp_word) DIGIT_BIT);
+ tmp |= ((mp_word) x.dp[i - 1]);
+ tmp /= ((mp_word) y.dp[t]);
+ if (tmp > (mp_word) MP_MASK)
+ tmp = MP_MASK;
+ q.dp[i - t - 1] = (mp_digit) (tmp & (mp_word) (MP_MASK));
+ }
+
+ /* step 3.2 while (q{i-t-1} * (yt * b + y{t-1})) > xi * b^2 + xi-1 * b + xi-2 do q{i-t-1} -= 1; */
+ q.dp[i - t - 1] = (q.dp[i - t - 1] + 1) & MP_MASK;
+ do {
+ q.dp[i - t - 1] = (q.dp[i - t - 1] - 1) & MP_MASK;
+
+ /* find left hand */
+ mp_zero (&t1);
+ t1.dp[0] = (t - 1 < 0) ? 0 : y.dp[t - 1];
+ t1.dp[1] = y.dp[t];
+ t1.used = 2;
+ if ((res = mp_mul_d (&t1, q.dp[i - t - 1], &t1)) != MP_OKAY) {
+ goto __Y;
+ }
+
+ /* find right hand */
+ t2.dp[0] = (i - 2 < 0) ? 0 : x.dp[i - 2];
+ t2.dp[1] = (i - 1 < 0) ? 0 : x.dp[i - 1];
+ t2.dp[2] = x.dp[i];
+ t2.used = 3;
+ }
+ while (mp_cmp (&t1, &t2) == MP_GT);
+
+ /* step 3.3 x = x - q{i-t-1} * y * b^{i-t-1} */
+ if ((res = mp_mul_d (&y, q.dp[i - t - 1], &t1)) != MP_OKAY) {
+ goto __Y;
+ }
+
+ if ((res = mp_lshd (&t1, i - t - 1)) != MP_OKAY) {
+ goto __Y;
+ }
+
+ if ((res = mp_sub (&x, &t1, &x)) != MP_OKAY) {
+ goto __Y;
+ }
+
+ /* step 3.4 if x < 0 then { x = x + y*b^{i-t-1}; q{i-t-1} -= 1; } */
+ if (x.sign == MP_NEG) {
+ if ((res = mp_copy (&y, &t1)) != MP_OKAY) {
+ goto __Y;
+ }
+ if ((res = mp_lshd (&t1, i - t - 1)) != MP_OKAY) {
+ goto __Y;
+ }
+ if ((res = mp_add (&x, &t1, &x)) != MP_OKAY) {
+ goto __Y;
+ }
+
+ q.dp[i - t - 1] = (q.dp[i - t - 1] - 1UL) & MP_MASK;
+ }
+ }
+
+ /* now q is the quotient and x is the remainder [which we have to normalize] */
+ /* get sign before writing to c */
+ x.sign = a->sign;
+ if (c != NULL) {
+ mp_clamp (&q);
+ mp_exch (&q, c);
+ c->sign = neg;
+ }
+
+ if (d != NULL) {
+ mp_div_2d (&x, norm, &x, NULL);
+ mp_clamp (&x);
+ mp_exch (&x, d);
+ }
+
+ res = MP_OKAY;
+
+__Y:mp_clear (&y);
+__X:mp_clear (&x);
+__T2:mp_clear (&t2);
+__T1:mp_clear (&t1);
+__Q:mp_clear (&q);
+ return res;
+}
View
38 bn_mp_div_2.c
@@ -0,0 +1,38 @@
+/* LibTomMath, multiple-precision integer library -- Tom St Denis
+ *
+ * LibTomMath is library that provides for multiple-precision
+ * integer arithmetic as well as number theoretic functionality.
+ *
+ * The library is 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://libtommath.iahu.ca
+ */
+#include <tommath.h>
+
+/* b = a/2 */
+int
+mp_div_2 (mp_int * a, mp_int * b)
+{
+ mp_digit r, rr;
+ int x, res;
+
+
+ /* copy */
+ if ((res = mp_copy (a, b)) != MP_OKAY) {
+ return res;
+ }
+
+ r = 0;
+ for (x = b->used - 1; x >= 0; x--) {
+ rr = b->dp[x] & 1;
+ b->dp[x] = (b->dp[x] >> 1) | (r << (DIGIT_BIT - 1));
+ r = rr;
+ }
+ mp_clamp (b);
+ return MP_OKAY;
+}
View
78 bn_mp_div_2d.c
@@ -0,0 +1,78 @@
+/* LibTomMath, multiple-precision integer library -- Tom St Denis
+ *
+ * LibTomMath is library that provides for multiple-precision
+ * integer arithmetic as well as number theoretic functionality.
+ *
+ * The library is 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://libtommath.iahu.ca
+ */
+#include <tommath.h>
+
+/* shift right by a certain bit count (store quotient in c, remainder in d) */
+int
+mp_div_2d (mp_int * a, int b, mp_int * c, mp_int * d)
+{
+ mp_digit D, r, rr;
+ int x, res;
+ mp_int t;
+
+
+ /* if the shift count is <= 0 then we do no work */
+ if (b <= 0) {
+ res = mp_copy (a, c);
+ if (d != NULL) {
+ mp_zero (d);
+ }
+ return res;
+ }
+
+ if ((res = mp_init (&t)) != MP_OKAY) {
+ return res;
+ }
+
+ /* get the remainder */
+ if (d != NULL) {
+ if ((res = mp_mod_2d (a, b, &t)) != MP_OKAY) {
+ mp_clear (&t);
+ return res;
+ }
+ }
+
+ /* copy */
+ if ((res = mp_copy (a, c)) != MP_OKAY) {
+ mp_clear (&t);
+ return res;
+ }
+
+ /* shift by as many digits in the bit count */
+ mp_rshd (c, b / DIGIT_BIT);
+
+ /* shift any bit count < DIGIT_BIT */
+ D = (mp_digit) (b % DIGIT_BIT);
+ if (D != 0) {
+ r = 0;
+ for (x = c->used - 1; x >= 0; x--) {
+ /* get the lower bits of this word in a temp */
+ rr = c->dp[x] & ((mp_digit) ((1U << D) - 1U));
+
+ /* shift the current word and mix in the carry bits from the previous word */
+ c->dp[x] = (c->dp[x] >> D) | (r << (DIGIT_BIT - D));
+
+ /* set the carry to the carry bits of the current word found above */
+ r = rr;
+ }
+ }
+ mp_clamp (c);
+ res = MP_OKAY;
+ if (d != NULL) {
+ mp_exch (&t, d);
+ }
+ mp_clear (&t);
+ return MP_OKAY;
+}
View
44 bn_mp_div_d.c
@@ -0,0 +1,44 @@
+/* LibTomMath, multiple-precision integer library -- Tom St Denis
+ *
+ * LibTomMath is library that provides for multiple-precision
+ * integer arithmetic as well as number theoretic functionality.
+ *
+ * The library is 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://libtommath.iahu.ca
+ */
+#include <tommath.h>
+
+/* single digit division */
+int
+mp_div_d (mp_int * a, mp_digit b, mp_int * c, mp_digit * d)
+{
+ mp_int t, t2;
+ int res;
+
+
+ if ((res = mp_init (&t)) != MP_OKAY) {
+ return res;
+ }
+
+ if ((res = mp_init (&t2)) != MP_OKAY) {
+ mp_clear (&t);
+ return res;
+ }
+
+ mp_set (&t, b);
+ res = mp_div (a, &t, c, &t2);
+
+ if (d != NULL) {
+ *d = t2.dp[0];
+ }
+
+ mp_clear (&t);
+ mp_clear (&t2);
+ return res;
+}
View
25 bn_mp_exch.c
@@ -0,0 +1,25 @@
+/* LibTomMath, multiple-precision integer library -- Tom St Denis
+ *
+ * LibTomMath is library that provides for multiple-precision
+ * integer arithmetic as well as number theoretic functionality.
+ *
+ * The library is 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://libtommath.iahu.ca
+ */
+#include <tommath.h>
+
+void
+mp_exch (mp_int * a, mp_int * b)
+{
+ mp_int t;
+
+ t = *a;
+ *a = *b;
+ *b = t;
+}
View
49 bn_mp_expt_d.c
@@ -0,0 +1,49 @@
+/* LibTomMath, multiple-precision integer library -- Tom St Denis
+ *
+ * LibTomMath is library that provides for multiple-precision
+ * integer arithmetic as well as number theoretic functionality.
+ *
+ * The library is 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://libtommath.iahu.ca
+ */
+#include <tommath.h>
+
+int
+mp_expt_d (mp_int * a, mp_digit b, mp_int * c)
+{
+ int res, x;
+ mp_int g;
+
+
+ if ((res = mp_init_copy (&g, a)) != MP_OKAY) {
+ return res;
+ }
+
+ /* set initial result */
+ mp_set (c, 1);
+
+ for (x = 0; x < (int) DIGIT_BIT; x++) {
+ if ((res = mp_sqr (c, c)) != MP_OKAY) {
+ mp_clear (&g);
+ return res;
+ }
+
+ if ((b & (mp_digit) (1 << (DIGIT_BIT - 1))) != 0) {
+ if ((res = mp_mul (c, &g, c)) != MP_OKAY) {
+ mp_clear (&g);
+ return res;
+ }
+ }
+
+ b <<= 1;
+ }
+
+ mp_clear (&g);
+ return MP_OKAY;
+}
View
213 bn_mp_exptmod.c
@@ -0,0 +1,213 @@
+/* LibTomMath, multiple-precision integer library -- Tom St Denis
+ *
+ * LibTomMath is library that provides for multiple-precision
+ * integer arithmetic as well as number theoretic functionality.
+ *
+ * The library is 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://libtommath.iahu.ca
+ */
+#include <tommath.h>
+
+int
+mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y)
+{
+ mp_int M[256], res, mu;
+ mp_digit buf;
+ int err, bitbuf, bitcpy, bitcnt, mode, digidx, x, y, winsize;
+
+
+ /* if the modulus is odd use the fast method */
+ if (mp_isodd (P) == 1 && P->used > 4 && P->used < MONTGOMERY_EXPT_CUTOFF) {
+ err = mp_exptmod_fast (G, X, P, Y);
+ return err;
+ }
+
+ /* 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;
+ }
+
+ /* init G array */
+ for (x = 0; x < (1 << winsize); x++) {
+ if ((err = mp_init_size (&M[x], 1)) != MP_OKAY) {
+ for (y = 0; y < x; y++) {
+ mp_clear (&M[y]);
+ }
+ return err;
+ }
+ }
+
+ /* create mu, used for Barrett reduction */
+ if ((err = mp_init (&mu)) != MP_OKAY) {
+ goto __M;
+ }
+ if ((err = mp_reduce_setup (&mu, P)) != MP_OKAY) {
+ goto __MU;
+ }
+
+ /* create M table
+ *
+ * The M table contains powers of the input 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 __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 __MU;
+ }
+
+ for (x = 0; x < (winsize - 1); x++) {
+ if ((err =
+ mp_sqr (&M[1 << (winsize - 1)],
+ &M[1 << (winsize - 1)])) != MP_OKAY) {
+ goto __MU;
+ }
+ if ((err = mp_reduce (&M[1 << (winsize - 1)], P, &mu)) != MP_OKAY) {
+ goto __MU;
+ }
+ }
+
+ /* create upper table */
+ for (x = (1 << (winsize - 1)) + 1; x < (1 << winsize); x++) {
+ if ((err = mp_mul (&M[x - 1], &M[1], &M[x])) != MP_OKAY) {
+ goto __MU;
+ }
+ if ((err = mp_reduce (&M[x], P, &mu)) != MP_OKAY) {
+ goto __MU;
+ }
+ }
+
+ /* setup result */
+ if ((err = mp_init (&res)) != MP_OKAY) {
+ goto __MU;
+ }
+ mp_set (&res, 1);
+
+ /* set initial mode and bit cnt */
+ mode = 0;
+ bitcnt = 0;
+ buf = 0;
+ digidx = X->used - 1;
+ bitcpy = bitbuf = 0;
+
+ bitcnt = 1;
+ for (;;) {
+ /* grab next digit as required */
+ if (--bitcnt == 0) {
+ if (digidx == -1) {
+ break;
+ }
+ buf = X->dp[digidx--];
+ bitcnt = (int) DIGIT_BIT;
+ }
+
+ /* grab the next msb from the exponent */
+ y = (buf >> (DIGIT_BIT - 1)) & 1;
+ buf <<= 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 __RES;
+ }
+ if ((err = mp_reduce (&res, P, &mu)) != MP_OKAY) {
+ goto __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 multiply */
+ /* square first */
+ for (x = 0; x < winsize; x++) {
+ if ((err = mp_sqr (&res, &res)) != MP_OKAY) {
+ goto __RES;
+ }
+ if ((err = mp_reduce (&res, P, &mu)) != MP_OKAY) {
+ goto __RES;
+ }
+ }
+
+ /* then multiply */
+ if ((err = mp_mul (&res, &M[bitbuf], &res)) != MP_OKAY) {
+ goto __MU;
+ }
+ if ((err = mp_reduce (&res, P, &mu)) != MP_OKAY) {
+ goto __MU;
+ }
+
+ /* empty window and reset */
+ bitcpy = 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 __RES;
+ }
+ if ((err = mp_reduce (&res, P, &mu)) != MP_OKAY) {
+ goto __RES;
+ }
+
+ bitbuf <<= 1;
+ if ((bitbuf & (1 << winsize)) != 0) {
+ /* then multiply */
+ if ((err = mp_mul (&res, &M[1], &res)) != MP_OKAY) {
+ goto __RES;
+ }
+ if ((err = mp_reduce (&res, P, &mu)) != MP_OKAY) {
+ goto __RES;
+ }
+ }
+ }
+ }
+
+ mp_exch (&res, Y);
+ err = MP_OKAY;
+__RES:mp_clear (&res);
+__MU:mp_clear (&mu);
+__M:
+ for (x = 0; x < (1 << winsize); x++) {
+ mp_clear (&M[x]);
+ }
+ return err;
+}
View
229 bn_mp_exptmod_fast.c
@@ -0,0 +1,229 @@
+/* LibTomMath, multiple-precision integer library -- Tom St Denis
+ *
+ * LibTomMath is library that provides for multiple-precision
+ * integer arithmetic as well as number theoretic functionality.
+ *
+ * The library is 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://libtommath.iahu.ca
+ */
+#include <tommath.h>
+
+/* computes Y == G^X mod P, HAC pp.616, Algorithm 14.85
+ *
+ * Uses a left-to-right k-ary sliding window to compute the modular exponentiation.
+ * The value of k changes based on the size of the exponent.
+ *
+ * Uses Montgomery reduction
+ */
+int
+mp_exptmod_fast (mp_int * G, mp_int * X, mp_int * P, mp_int * Y)
+{
+ mp_int M[256], res;
+ mp_digit buf, mp;
+ 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;
+ }
+
+ /* init G array */
+ for (x = 0; x < (1 << winsize); x++) {
+ if ((err = mp_init_size (&M[x], 1)) != MP_OKAY) {
+ for (y = 0; y < x; y++) {
+ mp_clear (&M[y]);
+ }
+ return err;
+ }
+ }
+
+ /* now setup montgomery */
+ if ((err = mp_montgomery_setup (P, &mp)) != MP_OKAY) {
+ goto __M;
+ }
+
+ /* setup result */
+ if ((err = mp_init (&res)) != MP_OKAY) {
+ goto __RES;
+ }
+
+ /* now we need R mod m */
+ if ((err = mp_2expt (&res, P->used * DIGIT_BIT)) != MP_OKAY) {
+ goto __RES;
+ }
+
+ /* res = R mod m (can use modified double/subtract ...) */
+ if ((err = mp_mod (&res, P, &res)) != MP_OKAY) {
+ goto __RES;
+ }
+
+ /* create M table
+ *
+ * The M table contains powers of the input 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 __RES;
+ }
+
+ /* now set M[1] to G * R mod m */
+ if ((err = mp_mulmod (&M[1], &res, P, &M[1])) != MP_OKAY) {
+ goto __RES;
+ }
+
+ /* 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 __RES;
+ }
+
+ for (x = 0; x < (winsize - 1); x++) {
+ if ((err =
+ mp_sqr (&M[1 << (winsize - 1)],
+ &M[1 << (winsize - 1)])) != MP_OKAY) {
+ goto __RES;
+ }
+ if ((err =
+ mp_montgomery_reduce (&M[1 << (winsize - 1)], P, mp)) != MP_OKAY) {
+ goto __RES;
+ }
+ }
+
+ /* create upper table */
+ for (x = (1 << (winsize - 1)) + 1; x < (1 << winsize); x++) {
+ if ((err = mp_mul (&M[x - 1], &M[1], &M[x])) != MP_OKAY) {
+ goto __RES;
+ }
+ if ((err = mp_montgomery_reduce (&M[x], P, mp)) != MP_OKAY) {
+ goto __RES;
+ }
+ }
+
+ /* set initial mode and bit cnt */
+ mode = 0;
+ bitcnt = 0;
+ buf = 0;
+ digidx = X->used - 1;
+ bitcpy = bitbuf = 0;
+
+ bitcnt = 1;
+ for (;;) {
+ /* grab next digit as required */
+ if (--bitcnt == 0) {
+ if (digidx == -1) {
+ break;
+ }
+ buf = X->dp[digidx--];
+ bitcnt = (int) DIGIT_BIT;
+ }
+
+ /* grab the next msb from the exponent */
+ y = (buf >> (DIGIT_BIT - 1)) & 1;
+ buf <<= 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 __RES;
+ }
+ if ((err = mp_montgomery_reduce (&res, P, mp)) != MP_OKAY) {
+ goto __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 multiply */
+ /* square first */
+ for (x = 0; x < winsize; x++) {
+ if ((err = mp_sqr (&res, &res)) != MP_OKAY) {
+ goto __RES;
+ }
+ if ((err = mp_montgomery_reduce (&res, P, mp)) != MP_OKAY) {
+ goto __RES;
+ }
+ }
+
+ /* then multiply */
+ if ((err = mp_mul (&res, &M[bitbuf], &res)) != MP_OKAY) {
+ goto __RES;
+ }
+ if ((err = mp_montgomery_reduce (&res, P, mp)) != MP_OKAY) {
+ goto __RES;
+ }
+
+ /* empty window and reset */
+ bitcpy = 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 __RES;
+ }
+ if ((err = mp_montgomery_reduce (&res, P, mp)) != MP_OKAY) {
+ goto __RES;
+ }
+
+ bitbuf <<= 1;
+ if ((bitbuf & (1 << winsize)) != 0) {
+ /* then multiply */
+ if ((err = mp_mul (&res, &M[1], &res)) != MP_OKAY) {
+ goto __RES;
+ }
+ if ((err = mp_montgomery_reduce (&res, P, mp)) != MP_OKAY) {
+ goto __RES;
+ }
+ }
+ }
+ }
+
+ /* fixup result */
+ if ((err = mp_montgomery_reduce (&res, P, mp)) != MP_OKAY) {
+ goto __RES;
+ }
+
+ mp_exch (&res, Y);
+ err = MP_OKAY;
+__RES:mp_clear (&res);
+__M:
+ for (x = 0; x < (1 << winsize); x++) {
+ mp_clear (&M[x]);
+ }
+ return err;
+}
View
118 bn_mp_gcd.c
@@ -0,0 +1,118 @@
+/* LibTomMath, multiple-precision integer library -- Tom St Denis
+ *
+ * LibTomMath is library that provides for multiple-precision
+ * integer arithmetic as well as number theoretic functionality.
+ *
+ * The library is 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://libtommath.iahu.ca
+ */
+#include <tommath.h>
+
+/* Greatest Common Divisor using the binary method [Algorithm B, page 338, vol2 of TAOCP]
+ */
+int
+mp_gcd (mp_int * a, mp_int * b, mp_int * c)
+{
+ mp_int u, v, t;
+ int k, res, neg;
+
+
+ /* either zero than gcd is the largest */
+ if (mp_iszero (a) == 1 && mp_iszero (b) == 0) {
+ return mp_copy (b, c);
+ }
+ if (mp_iszero (a) == 0 && mp_iszero (b) == 1) {
+ return mp_copy (a, c);
+ }
+ if (mp_iszero (a) == 1 && mp_iszero (b) == 1) {
+ mp_set (c, 1);
+ return MP_OKAY;
+ }
+
+ /* if both are negative they share (-1) as a common divisor */
+ neg = (a->sign == b->sign) ? a->sign : MP_ZPOS;
+
+ if ((res = mp_init_copy (&u, a)) != MP_OKAY) {
+ return res;
+ }
+
+ if ((res = mp_init_copy (&v, b)) != MP_OKAY) {
+ goto __U;
+ }
+
+ /* must be positive for the remainder of the algorithm */
+ u.sign = v.sign = MP_ZPOS;
+
+ if ((res = mp_init (&t)) != MP_OKAY) {
+ goto __V;
+ }
+
+ /* B1. Find power of two */
+ k = 0;
+ while ((u.dp[0] & 1) == 0 && (v.dp[0] & 1) == 0) {
+ ++k;
+ if ((res = mp_div_2d (&u, 1, &u, NULL)) != MP_OKAY) {
+ goto __T;
+ }
+ if ((res = mp_div_2d (&v, 1, &v, NULL)) != MP_OKAY) {
+ goto __T;
+ }
+ }
+
+ /* B2. Initialize */
+ if ((u.dp[0] & 1) == 1) {
+ if ((res = mp_copy (&v, &t)) != MP_OKAY) {
+ goto __T;
+ }
+ t.sign = MP_NEG;
+ } else {
+ if ((res = mp_copy (&u, &t)) != MP_OKAY) {
+ goto __T;
+ }
+ }
+
+ do {
+ /* B3 (and B4). Halve t, if even */
+ while (t.used != 0 && (t.dp[0] & 1) == 0) {
+ if ((res = mp_div_2d (&t, 1, &t, NULL)) != MP_OKAY) {
+ goto __T;
+ }
+ }
+
+ /* B5. if t>0 then u=t otherwise v=-t */
+ if (t.used != 0 && t.sign != MP_NEG) {
+ if ((res = mp_copy (&t, &u)) != MP_OKAY) {
+ goto __T;
+ }
+ } else {
+ if ((res = mp_copy (&t, &v)) != MP_OKAY) {
+ goto __T;
+ }
+ v.sign = (v.sign == MP_ZPOS) ? MP_NEG : MP_ZPOS;
+ }
+
+ /* B6. t = u - v, if t != 0 loop otherwise terminate */
+ if ((res = mp_sub (&u, &v, &t)) != MP_OKAY) {
+ goto __T;
+ }
+ }
+ while (t.used != 0);
+
+ if ((res = mp_mul_2d (&u, k, &u)) != MP_OKAY) {
+ goto __T;
+ }
+
+ mp_exch (&u, c);
+ c->sign = neg;
+ res = MP_OKAY;
+__T:mp_clear (&t);
+__V:mp_clear (&u);
+__U:mp_clear (&v);
+ return res;
+}
View
40 bn_mp_grow.c
@@ -0,0 +1,40 @@
+/* LibTomMath, multiple-precision integer library -- Tom St Denis
+ *
+ * LibTomMath is library that provides for multiple-precision
+ * integer arithmetic as well as number theoretic functionality.
+ *
+ * The library is 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://libtommath.iahu.ca
+ */
+#include <tommath.h>
+
+/* grow as required */
+int
+mp_grow (mp_int * a, int size)
+{
+ int i, n;
+
+
+ /* if the alloc size is smaller alloc more ram */
+ if (a->alloc < size) {
+ size += (MP_PREC * 2) - (size & (MP_PREC - 1)); /* ensure there are always at least MP_PREC digits extra on top */
+
+ a->dp = realloc (a->dp, sizeof (mp_digit) * size);
+ if (a->dp == NULL) {
+ return MP_MEM;
+ }
+
+ n = a->alloc;
+ a->alloc = size;
+ for (i = n; i < a->alloc; i++) {
+ a->dp[i] = 0;
+ }
+ }
+ return MP_OKAY;
+}
View
35 bn_mp_init.c
@@ -0,0 +1,35 @@
+/* LibTomMath, multiple-precision integer library -- Tom St Denis
+ *
+ * LibTomMath is library that provides for multiple-precision
+ * integer arithmetic as well as number theoretic functionality.
+ *
+ * The library is 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://libtommath.iahu.ca
+ */
+#include <tommath.h>
+
+/* init a new bigint */
+int
+mp_init (mp_int * a)
+{
+
+ /* allocate ram required and clear it */
+ a->dp = calloc (sizeof (mp_digit), MP_PREC);
+ if (a->dp == NULL) {
+ return MP_MEM;
+ }
+
+ /* set the used to zero, allocated digit to the default precision
+ * and sign to positive */
+ a->used = 0;
+ a->alloc = MP_PREC;
+ a->sign = MP_ZPOS;
+
+ return MP_OKAY;
+}
View
28 bn_mp_init_copy.c
@@ -0,0 +1,28 @@
+/* LibTomMath, multiple-precision integer library -- Tom St Denis
+ *
+ * LibTomMath is library that provides for multiple-precision
+ * integer arithmetic as well as number theoretic functionality.
+ *
+ * The library is 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://libtommath.iahu.ca
+ */
+#include <tommath.h>
+
+/* creates "a" then copies b into it */
+int
+mp_init_copy (mp_int * a, mp_int * b)
+{
+ int res;
+
+ if ((res = mp_init (a)) != MP_OKAY) {
+ return res;
+ }
+ res = mp_copy (b, a);
+ return res;
+}
View
33 bn_mp_init_size.c
@@ -0,0 +1,33 @@
+/* LibTomMath, multiple-precision integer library -- Tom St Denis
+ *