Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with
or
.
Download ZIP
Browse files

added libtommath-0.17

  • Loading branch information...
commit fd181cc841b4dde061fb62ea1d1b74ac92b0141c 1 parent 14161e8
Tom St Denis authored sjaeckel committed
Showing with 14,772 additions and 6,950 deletions.
  1. BIN  bn.pdf
  2. +2 −2 bn.tex
  3. +22 −51 bn_fast_mp_invmod.c
  4. +2 −2 bn_fast_mp_montgomery_reduce.c
  5. +20 −18 bn_fast_s_mp_mul_digs.c
  6. +12 −4 bn_fast_s_mp_mul_high_digs.c
  7. +9 −9 bn_fast_s_mp_sqr.c
  8. +13 −21 bn_mp_add.c
  9. +11 −2 bn_mp_cmp.c
  10. +6 −2 bn_mp_cmp_mag.c
  11. +4 −5 bn_mp_copy.c
  12. +11 −11 bn_mp_div.c
  13. +4 −4 bn_mp_div_2.c
  14. +5 −5 bn_mp_div_2d.c
  15. +34 −0 bn_mp_dr_is_modulus.c
  16. +17 −38 bn_mp_dr_reduce.c
  17. +25 −0 bn_mp_dr_setup.c
  18. +4 −4 bn_mp_expt_d.c
  19. +74 −33 bn_mp_exptmod.c
  20. +49 −34 bn_mp_exptmod_fast.c
  21. +5 −5 bn_mp_gcd.c
  22. +4 −4 bn_mp_grow.c
  23. +0 −1  bn_mp_init.c
  24. +28 −63 bn_mp_invmod.c
  25. +1 −1  bn_mp_jacobi.c
  26. +18 −21 bn_mp_karatsuba_mul.c
  27. +11 −11 bn_mp_karatsuba_sqr.c
  28. +4 −3 bn_mp_lshd.c
  29. +4 −4 bn_mp_montgomery_calc_normalization.c
  30. +15 −8 bn_mp_montgomery_reduce.c
  31. +15 −8 bn_mp_montgomery_setup.c
  32. +5 −5 bn_mp_mul.c
  33. +15 −22 bn_mp_mul_2.c
  34. +23 −12 bn_mp_mul_2d.c
  35. +20 −6 bn_mp_mul_d.c
  36. +64 −0 bn_mp_multi.c
  37. +2 −2 bn_mp_prime_is_divisible.c
  38. +9 −1 bn_mp_prime_is_prime.c
  39. +7 −7 bn_mp_prime_next_prime.c
  40. +12 −10 bn_mp_reduce.c
  41. +4 −3 bn_mp_rshd.c
  42. +3 −6 bn_mp_set_int.c
  43. +1 −2  bn_mp_sqr.c
  44. +19 −24 bn_mp_sub.c
  45. +4 −1 bn_prime_tab.c
  46. +77 −0 bn_radix.c
  47. +1 −1  bn_reverse.c
  48. +14 −15 bn_s_mp_add.c
  49. +21 −5 bn_s_mp_mul_digs.c
  50. +3 −3 bn_s_mp_mul_high_digs.c
  51. +5 −6 bn_s_mp_sub.c
  52. +12 −4 bncore.c
  53. +261 −0 booker.pl
  54. +34 −0 changes.txt
  55. +160 −181 demo/demo.c
  56. 0  demo/test.c
  57. +19 −2 etc/makefile
  58. +45 −0 etc/mont.c
  59. +37 −0 etc/timer.asm
  60. +37 −80 etc/tune.c
  61. +14 −23 gen.pl
  62. +13 −0 logs/README
  63. +16 −0 logs/add.log
  64. BIN  logs/addsub.png
  65. +7 −0 logs/expt.log
  66. BIN  logs/expt.png
  67. +7 −0 logs/expt_dr.log
  68. +17 −0 logs/graphs.dem
  69. +24 −0 logs/index.html
  70. +32 −0 logs/invmod.log
  71. BIN  logs/invmod.png
  72. +17 −0 logs/mult.log
  73. BIN  logs/mult.png
  74. +17 −0 logs/mult_kara.log
  75. +17 −0 logs/sqr.log
  76. +17 −0 logs/sqr_kara.log
  77. +16 −0 logs/sub.log
  78. +35 −9 makefile
  79. +2 −1  makefile.msvc
  80. +77 −57 mtest/mtest.c
  81. +17 −0 pics/makefile
  82. BIN  pics/sliding_window.TIF
  83. BIN  pics/sliding_window.sxd
  84. +6,356 −6,051 pre_gen/mpi.c
  85. +74 −37 tommath.h
  86. +2,459 −0 tommath.src
  87. +4,195 −0 tommath.tex
View
BIN  bn.pdf
Binary file not shown
View
4 bn.tex
@@ -1,7 +1,7 @@
-\documentclass[]{report}
+\documentclass[]{article}
\begin{document}
-\title{LibTomMath v0.16 \\ A Free Multiple Precision Integer Library \\ http://math.libtomcrypt.org }
+\title{LibTomMath v0.17 \\ A Free Multiple Precision Integer Library \\ http://math.libtomcrypt.org }
\author{Tom St Denis \\ tomstdenis@iahu.ca}
\maketitle
\newpage
View
73 bn_fast_mp_invmod.c
@@ -27,41 +27,18 @@ fast_mp_invmod (mp_int * a, mp_int * b, mp_int * c)
int res, neg;
/* init all our temps */
- 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;
+ if ((res = mp_init_multi(&x, &y, &u, &v, &B, &D, NULL)) != MP_OKAY) {
+ return res;
}
/* 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;
+ goto __ERR;
}
- /* we need |y| */
- if ((res = mp_abs (&y, &y)) != MP_OKAY) {
- goto __D;
+ /* we need y = |a| */
+ if ((res = mp_abs (a, &y)) != MP_OKAY) {
+ goto __ERR;
}
/* 2. [modified] if x,y are both even then return an error!
@@ -70,15 +47,15 @@ fast_mp_invmod (mp_int * a, mp_int * b, mp_int * c)
*/
if (mp_iseven (&x) == 1 && mp_iseven (&y) == 1) {
res = MP_VAL;
- goto __D;
+ goto __ERR;
}
/* 3. u=x, v=y, A=1, B=0, C=0,D=1 */
if ((res = mp_copy (&x, &u)) != MP_OKAY) {
- goto __D;
+ goto __ERR;
}
if ((res = mp_copy (&y, &v)) != MP_OKAY) {
- goto __D;
+ goto __ERR;
}
mp_set (&D, 1);
@@ -87,17 +64,17 @@ fast_mp_invmod (mp_int * a, mp_int * b, mp_int * c)
while (mp_iseven (&u) == 1) {
/* 4.1 u = u/2 */
if ((res = mp_div_2 (&u, &u)) != MP_OKAY) {
- goto __D;
+ goto __ERR;
}
/* 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;
+ goto __ERR;
}
}
/* B = B/2 */
if ((res = mp_div_2 (&B, &B)) != MP_OKAY) {
- goto __D;
+ goto __ERR;
}
}
@@ -105,18 +82,18 @@ fast_mp_invmod (mp_int * a, mp_int * b, mp_int * c)
while (mp_iseven (&v) == 1) {
/* 5.1 v = v/2 */
if ((res = mp_div_2 (&v, &v)) != MP_OKAY) {
- goto __D;
+ goto __ERR;
}
/* 5.2 if C,D are even then */
if (mp_iseven (&D) == 0) {
/* D = (D-x)/2 */
if ((res = mp_sub (&D, &x, &D)) != MP_OKAY) {
- goto __D;
+ goto __ERR;
}
}
/* D = D/2 */
if ((res = mp_div_2 (&D, &D)) != MP_OKAY) {
- goto __D;
+ goto __ERR;
}
}
@@ -124,20 +101,20 @@ fast_mp_invmod (mp_int * a, mp_int * b, mp_int * c)
if (mp_cmp (&u, &v) != MP_LT) {
/* u = u - v, B = B - D */
if ((res = mp_sub (&u, &v, &u)) != MP_OKAY) {
- goto __D;
+ goto __ERR;
}
if ((res = mp_sub (&B, &D, &B)) != MP_OKAY) {
- goto __D;
+ goto __ERR;
}
} else {
/* v - v - u, D = D - B */
if ((res = mp_sub (&v, &u, &v)) != MP_OKAY) {
- goto __D;
+ goto __ERR;
}
if ((res = mp_sub (&D, &B, &D)) != MP_OKAY) {
- goto __D;
+ goto __ERR;
}
}
@@ -151,26 +128,20 @@ fast_mp_invmod (mp_int * a, mp_int * b, mp_int * c)
/* if v != 1 then there is no inverse */
if (mp_cmp_d (&v, 1) != MP_EQ) {
res = MP_VAL;
- goto __D;
+ goto __ERR;
}
/* b is now the inverse */
neg = a->sign;
while (D.sign == MP_NEG) {
if ((res = mp_add (&D, b, &D)) != MP_OKAY) {
- goto __D;
+ goto __ERR;
}
}
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:
+__ERR:mp_clear_multi (&x, &y, &u, &v, &B, &D, NULL);
return res;
}
View
4 bn_fast_mp_montgomery_reduce.c
@@ -26,7 +26,7 @@ int
fast_mp_montgomery_reduce (mp_int * a, mp_int * m, mp_digit mp)
{
int ix, res, olduse;
- mp_word W[512];
+ mp_word W[MP_WARRAY];
/* get old used count */
olduse = a->used;
@@ -92,7 +92,7 @@ fast_mp_montgomery_reduce (mp_int * a, mp_int * m, mp_digit mp)
/* inner loop */
for (iy = 0; iy < m->used; iy++) {
- *_W++ += ((mp_word) ui) * ((mp_word) * tmpx++);
+ *_W++ += ((mp_word) ui) * ((mp_word) * tmpx++);
}
}
View
38 bn_fast_s_mp_mul_digs.c
@@ -16,14 +16,16 @@
/* 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
+ * 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.
*
- * This has been modified to produce a variable number of digits of output so
- * if say only a half-product is required you don't have to compute the upper half
- * (a feature required for fast Barrett reduction).
+ * This has been modified to produce a variable number of
+ * digits of output so if say only a half-product is required
+ * you don't have to compute the upper half (a feature
+ * required for fast Barrett reduction).
*
* Based on Algorithm 14.12 on pp.595 of HAC.
*
@@ -32,7 +34,7 @@ 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];
+ mp_word W[MP_WARRAY];
/* grow the destination as required */
if (c->alloc < digs) {
@@ -47,10 +49,9 @@ fast_s_mp_mul_digs (mp_int * a, mp_int * b, mp_int * c, int 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 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
@@ -73,14 +74,14 @@ fast_s_mp_mul_digs (mp_int * a, mp_int * b, mp_int * c, int digs)
*/
_W = W + ix;
- /* the number of digits is limited by their placement. E.g.
+ /* 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++);
+ *_W++ += ((mp_word) tmpx) * ((mp_word) * tmpy++);
}
}
@@ -97,11 +98,12 @@ fast_s_mp_mul_digs (mp_int * a, mp_int * b, mp_int * c, int digs)
* 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.
+ * 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.
*/
tmpc = c->dp;
for (ix = 1; ix < digs; ix++) {
View
16 bn_fast_s_mp_mul_high_digs.c
@@ -27,7 +27,7 @@ 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];
+ mp_word W[MP_WARRAY];
/* calculate size of product and allocate more space if required */
newused = a->used + b->used + 1;
@@ -55,15 +55,23 @@ fast_s_mp_mul_high_digs (mp_int * a, mp_int * b, mp_int * c, int digs)
/* 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]);
+ _W = W + digs;
+
+ /* skip cases below zero where ix > digs */
+ if (iy < 0) {
+ iy = abs(iy);
+ tmpy += iy;
+ _W += iy;
+ iy = 0;
+ }
/* compute column products for digits above the minimum */
for (; iy < pb; iy++) {
- *_W++ += ((mp_word) tmpx) * ((mp_word) * tmpy++);
+ *_W++ += ((mp_word) tmpx) * ((mp_word) * tmpy++);
}
}
}
View
18 bn_fast_s_mp_sqr.c
@@ -20,7 +20,7 @@
* 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.
+ * 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
@@ -33,7 +33,7 @@ int
fast_s_mp_sqr (mp_int * a, mp_int * b)
{
int olduse, newused, res, ix, pa;
- mp_word W2[512], W[512];
+ mp_word W2[MP_WARRAY], W[MP_WARRAY];
/* calculate size of product and allocate as required */
pa = a->used;
@@ -44,9 +44,9 @@ fast_s_mp_sqr (mp_int * a, mp_int * b)
}
}
- /* zero temp buffer (columns)
+ /* zero temp buffer (columns)
* Note that there are two buffers. Since squaring requires
- * a outter and inner product and the inner product 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
@@ -60,7 +60,7 @@ fast_s_mp_sqr (mp_int * a, mp_int * b)
* values in W2 are only written in even locations which means
* we can collapse the array to 256 words [and fixup the memset above]
* provided we also fix up the summations below. Ideally
- * the fixup loop should be unrolled twice to handle the even/odd
+ * the fixup loop should be unrolled twice to handle the even/odd
* cases, and then a final step to handle odd cases [e.g. newused == odd]
*
* This will not only save ~8*256 = 2KB of stack but lower the number of
@@ -71,10 +71,10 @@ fast_s_mp_sqr (mp_int * a, mp_int * b)
* the multiplication by two is done afterwards in the N loop.
*/
for (ix = 0; ix < pa; ix++) {
- /* compute the outer product
+ /* compute the outer product
*
- * Note that every outer product is computed
- * for a particular column only once which means that
+ * Note that every outer product is computed
+ * for a particular column only once which means that
* there is no need todo a double precision addition
*/
W2[ix + ix] = ((mp_word) a->dp[ix]) * ((mp_word) a->dp[ix]);
@@ -95,7 +95,7 @@ fast_s_mp_sqr (mp_int * a, mp_int * b)
/* inner products */
for (iy = ix + 1; iy < pa; iy++) {
- *_W++ += ((mp_word) tmpx) * ((mp_word) * tmpy++);
+ *_W++ += ((mp_word) tmpx) * ((mp_word) * tmpy++);
}
}
}
View
34 bn_mp_add.c
@@ -24,33 +24,25 @@ mp_add (mp_int * a, mp_int * b, mp_int * c)
sa = a->sign;
sb = b->sign;
- /* handle four cases */
- if (sa == MP_ZPOS && sb == MP_ZPOS) {
- /* both positive */
+ /* handle two cases, not four */
+ if (sa == sb) {
+ /* both positive or both negative */
+ /* add their magnitudes, copy the sign */
+ c->sign = sa;
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) */
+ } else {
+ /* one positive, the other negative */
+ /* subtract the one with the greater magnitude from */
+ /* the one of the lesser magnitude. The result gets */
+ /* the sign of the one with the greater magnitude. */
if (mp_cmp_mag (a, b) == MP_LT) {
+ c->sign = sb;
res = s_mp_sub (b, a, c);
- c->sign = MP_NEG;
} else {
+ c->sign = sa;
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
13 bn_mp_cmp.c
@@ -21,8 +21,17 @@ mp_cmp (mp_int * a, mp_int * b)
/* 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) {
+ }
+
+ if (a->sign == MP_ZPOS && b->sign == MP_NEG) {
return MP_GT;
}
- return mp_cmp_mag (a, b);
+
+ /* compare digits */
+ if (a->sign == MP_NEG) {
+ /* if negative compare opposite direction */
+ return mp_cmp_mag(b, a);
+ } else {
+ return mp_cmp_mag(a, b);
+ }
}
View
8 bn_mp_cmp_mag.c
@@ -23,7 +23,9 @@ mp_cmp_mag (mp_int * a, mp_int * b)
/* compare based on # of non-zero digits */
if (a->used > b->used) {
return MP_GT;
- } else if (a->used < b->used) {
+ }
+
+ if (a->used < b->used) {
return MP_LT;
}
@@ -31,7 +33,9 @@ mp_cmp_mag (mp_int * a, mp_int * b)
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]) {
+ }
+
+ if (a->dp[n] < b->dp[n]) {
return MP_LT;
}
}
View
9 bn_mp_copy.c
@@ -31,13 +31,10 @@ mp_copy (mp_int * a, mp_int * b)
}
/* zero b and copy the parameters over */
- b->used = a->used;
- b->sign = a->sign;
-
{
register mp_digit *tmpa, *tmpb;
- /* point aliases */
+ /* pointer aliases */
tmpa = a->dp;
tmpb = b->dp;
@@ -47,9 +44,11 @@ mp_copy (mp_int * a, mp_int * b)
}
/* clear high digits */
- for (; n < b->alloc; n++) {
+ for (; n < b->used; n++) {
*tmpb++ = 0;
}
}
+ b->used = a->used;
+ b->sign = a->sign;
return MP_OKAY;
}
View
22 bn_mp_div.c
@@ -75,7 +75,7 @@ mp_div (mp_int * a, mp_int * b, mp_int * c, mp_int * d)
/* normalize both x and y, ensure that y >= b/2, [b == 2^DIGIT_BIT] */
norm = mp_count_bits(&y) % DIGIT_BIT;
- if (norm < (DIGIT_BIT-1)) {
+ if (norm < (int)(DIGIT_BIT-1)) {
norm = (DIGIT_BIT-1) - norm;
if ((res = mp_mul_2d (&x, norm, &x)) != MP_OKAY) {
goto __Y;
@@ -86,13 +86,13 @@ mp_div (mp_int * a, mp_int * b, mp_int * c, mp_int * d)
} else {
norm = 0;
}
-
+
/* 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} */
+ if ((res = mp_lshd (&y, n - t)) != MP_OKAY) { /* y = y*b^{n-t} */
goto __Y;
}
@@ -113,14 +113,14 @@ mp_div (mp_int * a, mp_int * b, mp_int * c, mp_int * d)
/* 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);
+ q.dp[i - t - 1] = ((((mp_digit)1) << DIGIT_BIT) - 1);
} 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;
+ tmp = MP_MASK;
q.dp[i - t - 1] = (mp_digit) (tmp & (mp_word) (MP_MASK));
}
@@ -135,7 +135,7 @@ mp_div (mp_int * a, mp_int * b, mp_int * c, mp_int * d)
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;
+ goto __Y;
}
/* find right hand */
@@ -143,7 +143,7 @@ mp_div (mp_int * a, mp_int * b, mp_int * c, mp_int * d)
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);
+ } while (mp_cmp_mag(&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) {
@@ -161,19 +161,19 @@ mp_div (mp_int * a, mp_int * b, mp_int * c, mp_int * d)
/* 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;
+ goto __Y;
}
if ((res = mp_lshd (&t1, i - t - 1)) != MP_OKAY) {
- goto __Y;
+ goto __Y;
}
if ((res = mp_add (&x, &t1, &x)) != MP_OKAY) {
- goto __Y;
+ 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;
View
8 bn_mp_div_2.c
@@ -34,19 +34,19 @@ mp_div_2 (mp_int * a, mp_int * b)
/* source alias */
tmpa = a->dp + b->used - 1;
-
+
/* dest alias */
tmpb = b->dp + b->used - 1;
-
+
/* carry */
r = 0;
for (x = b->used - 1; x >= 0; x--) {
/* get the carry for the next iteration */
rr = *tmpa & 1;
-
+
/* shift the current digit, add in carry and store */
*tmpb-- = (*tmpa-- >> 1) | (r << (DIGIT_BIT - 1));
-
+
/* forward carry to next iteration */
r = rr;
}
View
10 bn_mp_div_2d.c
@@ -51,7 +51,7 @@ mp_div_2d (mp_int * a, int b, mp_int * c, mp_int * d)
}
/* shift by as many digits in the bit count */
- if (b >= DIGIT_BIT) {
+ if (b >= (int)DIGIT_BIT) {
mp_rshd (c, b / DIGIT_BIT);
}
@@ -59,13 +59,13 @@ mp_div_2d (mp_int * a, int b, mp_int * c, mp_int * d)
D = (mp_digit) (b % DIGIT_BIT);
if (D != 0) {
register mp_digit *tmpc, mask;
-
+
/* mask */
- mask = (1U << D) - 1U;
-
+ mask = (((mp_digit)1) << D) - 1;
+
/* alias */
tmpc = c->dp + (c->used - 1);
-
+
/* carry */
r = 0;
for (x = c->used - 1; x >= 0; x--) {
View
34 bn_mp_dr_is_modulus.c
@@ -0,0 +1,34 @@
+/* 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://math.libtomcrypt.org
+ */
+#include <tommath.h>
+
+/* determines if a number is a valid DR modulus */
+int mp_dr_is_modulus(mp_int *a)
+{
+ int ix;
+
+ /* must be at least two digits */
+ if (a->used < 2) {
+ return 0;
+ }
+
+ for (ix = 1; ix < a->used; ix++) {
+ if (a->dp[ix] != MP_MASK) {
+ return 0;
+ }
+ }
+ return 1;
+}
+
View
55 bn_mp_dr_reduce.c
@@ -16,7 +16,7 @@
/* reduce "a" in place modulo "b" using the Diminished Radix algorithm.
*
- * Based on algorithm from the paper
+ * Based on algorithm from the paper
*
* "Generating Efficient Primes for Discrete Log Cryptosystems"
* Chae Hoon Lim, Pil Loong Lee,
@@ -40,15 +40,15 @@ mp_dr_reduce (mp_int * a, mp_int * b, mp_digit mp)
return err;
}
}
-
+
/* alias for a->dp[i] */
tmpi = a->dp + k + k - 1;
- /* for (i = 2k - 1; i >= k; i = i - 1)
+ /* for (i = 2k - 1; i >= k; i = i - 1)
*
* This is the main loop of the reduction. Note that at the end
* the words above position k are not zeroed as expected. The end
- * result is that the digits from 0 to k-1 are the residue. So
+ * result is that the digits from 0 to k-1 are the residue. So
* we have to clear those afterwards.
*/
for (i = k + k - 1; i >= k; i = i - 1) {
@@ -57,10 +57,10 @@ mp_dr_reduce (mp_int * a, mp_int * b, mp_digit mp)
/* x[i] * mp */
r = ((mp_word) *tmpi--) * ((mp_word) mp);
- /* now add r to x[i-1:i-k]
+ /* now add r to x[i-1:i-k]
*
* First add it to the first digit x[i-k] then form the carry
- * then enter the main loop
+ * then enter the main loop
*/
j = i - k;
@@ -74,14 +74,14 @@ mp_dr_reduce (mp_int * a, mp_int * b, mp_digit mp)
mu = (r >> ((mp_word) DIGIT_BIT)) + (*tmpj >> DIGIT_BIT);
/* clear carry from a->dp[j] */
- *tmpj++ &= MP_MASK;
+ *tmpj++ &= MP_MASK;
- /* now add rest of the digits
- *
+ /* now add rest of the digits
+ *
* Note this is basically a simple single digit addition to
* a larger multiple digit number. This is optimized somewhat
* because the propagation of carries is not likely to move
- * more than a few digits.
+ * more than a few digits.
*
*/
for (++j; mu != 0 && j <= (i - 1); ++j) {
@@ -99,16 +99,16 @@ mp_dr_reduce (mp_int * a, mp_int * b, mp_digit mp)
*tmpj += mp;
mu = *tmpj >> DIGIT_BIT;
*tmpj++ &= MP_MASK;
-
+
/* now handle carries */
for (++j; mu != 0 && j <= (i - 1); j++) {
- *tmpj += mu;
- mu = *tmpj >> DIGIT_BIT;
- *tmpj++ &= MP_MASK;
+ *tmpj += mu;
+ mu = *tmpj >> DIGIT_BIT;
+ *tmpj++ &= MP_MASK;
}
}
}
-
+
/* zero words above k */
tmpi = a->dp + k;
for (i = k; i < a->used; i++) {
@@ -117,34 +117,13 @@ mp_dr_reduce (mp_int * a, mp_int * b, mp_digit mp)
/* clamp, sub and return */
mp_clamp (a);
-
+
+ /* if a >= b [b == modulus] then subtract the modulus to fix up */
if (mp_cmp_mag (a, b) != MP_LT) {
return s_mp_sub (a, b, a);
}
return MP_OKAY;
}
-/* determines if a number is a valid DR modulus */
-int mp_dr_is_modulus(mp_int *a)
-{
- int ix;
-
- /* must be at least two digits */
- if (a->used < 2) {
- return 0;
- }
-
- for (ix = 1; ix < a->used; ix++) {
- if (a->dp[ix] != MP_MASK) {
- return 0;
- }
- }
- return 1;
-}
-/* determines the setup value */
-void mp_dr_setup(mp_int *a, mp_digit *d)
-{
- *d = (1 << DIGIT_BIT) - a->dp[0];
-}
View
25 bn_mp_dr_setup.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://math.libtomcrypt.org
+ */
+#include <tommath.h>
+
+/* determines the setup value */
+void mp_dr_setup(mp_int *a, mp_digit *d)
+{
+ /* the casts are required if DIGIT_BIT is one less than
+ * the number of bits in a mp_digit [e.g. DIGIT_BIT==31]
+ */
+ *d = (mp_digit)((((mp_word)1) << ((mp_word)DIGIT_BIT)) - ((mp_word)a->dp[0]));
+}
+
View
8 bn_mp_expt_d.c
@@ -35,11 +35,11 @@ mp_expt_d (mp_int * a, mp_digit b, mp_int * c)
return res;
}
- /* if the bit is set multiply */
- if ((b & (mp_digit) (1 << (DIGIT_BIT - 1))) != 0) {
+ /* if the bit is set multiply */
+ if ((b & (mp_digit) (((mp_digit)1) << (DIGIT_BIT - 1))) != 0) {
if ((res = mp_mul (c, &g, c)) != MP_OKAY) {
- mp_clear (&g);
- return res;
+ mp_clear (&g);
+ return res;
}
}
View
107 bn_mp_exptmod.c
@@ -17,7 +17,7 @@
static int f_mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y);
/* this is a shell function that calls either the normal or Montgomery
- * exptmod functions. Originally the call to the montgomery code was
+ * exptmod functions. Originally the call to the montgomery code was
* embedded in the normal function but that wasted alot of stack space
* for nothing (since 99% of the time the Montgomery code would be called)
*/
@@ -25,10 +25,46 @@ int
mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y)
{
int dr;
-
+
+ /* modulus P must be positive */
+ if (P->sign == MP_NEG) {
+ return MP_VAL;
+ }
+
+ /* if exponent X is negative we have to recurse */
+ if (X->sign == MP_NEG) {
+ mp_int tmpG, tmpX;
+ int err;
+
+ /* first compute 1/G mod P */
+ if ((err = mp_init(&tmpG)) != MP_OKAY) {
+ return err;
+ }
+ if ((err = mp_invmod(G, P, &tmpG)) != MP_OKAY) {
+ mp_clear(&tmpG);
+ return err;
+ }
+
+ /* now get |X| */
+ if ((err = mp_init(&tmpX)) != MP_OKAY) {
+ mp_clear(&tmpG);
+ return err;
+ }
+ if ((err = mp_abs(X, &tmpX)) != MP_OKAY) {
+ mp_clear_multi(&tmpG, &tmpX, NULL);
+ return err;
+ }
+
+ /* and now compute (1/G)^|X| instead of G^X [X < 0] */
+ err = mp_exptmod(&tmpG, &tmpX, P, Y);
+ mp_clear_multi(&tmpG, &tmpX, NULL);
+ return err;
+ }
+
+
dr = mp_dr_is_modulus(P);
/* if the modulus is odd use the fast method */
- if (((mp_isodd (P) == 1 && P->used < MONTGOMERY_EXPT_CUTOFF) || dr == 1) && P->used > 4) {
+ if ((mp_isodd (P) == 1 || dr == 1) && P->used > 4) {
return mp_exptmod_fast (G, X, P, Y, dr);
} else {
return f_mp_exptmod (G, X, P, Y);
@@ -60,11 +96,17 @@ f_mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y)
winsize = 8;
}
+#ifdef MP_LOW_MEM
+ if (winsize > 5) {
+ winsize = 5;
+ }
+#endif
+
/* 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]);
+ mp_clear (&M[y]);
}
return err;
}
@@ -78,7 +120,7 @@ f_mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y)
goto __MU;
}
- /* create M table
+ /* create M table
*
* The M table contains powers of the input base, e.g. M[x] = G^x mod P
*
@@ -119,30 +161,29 @@ f_mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y)
mp_set (&res, 1);
/* set initial mode and bit cnt */
- mode = 0;
- bitcnt = 0;
- buf = 0;
+ mode = 0;
+ bitcnt = 1;
+ buf = 0;
digidx = X->used - 1;
bitcpy = bitbuf = 0;
- bitcnt = 1;
for (;;) {
/* grab next digit as required */
if (--bitcnt == 0) {
if (digidx == -1) {
- break;
+ break;
}
buf = X->dp[digidx--];
bitcnt = (int) DIGIT_BIT;
}
/* grab the next msb from the exponent */
- y = (buf >> (DIGIT_BIT - 1)) & 1;
- buf <<= 1;
+ y = (buf >> (mp_digit)(DIGIT_BIT - 1)) & 1;
+ buf <<= (mp_digit)1;
- /* if the bit is zero and mode == 0 then we ignore it
+ /* 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
+ * in the exponent. Technically this opt is not required but it
* does lower the # of trivial squaring/reductions used
*/
if (mode == 0 && y == 0)
@@ -151,10 +192,10 @@ f_mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y)
/* 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;
+ goto __RES;
}
if ((err = mp_reduce (&res, P, &mu)) != MP_OKAY) {
- goto __RES;
+ goto __RES;
}
continue;
}
@@ -167,20 +208,20 @@ f_mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y)
/* 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 __RES;
- }
- if ((err = mp_reduce (&res, P, &mu)) != MP_OKAY) {
- goto __RES;
- }
+ 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;
+ goto __MU;
}
if ((err = mp_reduce (&res, P, &mu)) != MP_OKAY) {
- goto __MU;
+ goto __MU;
}
/* empty window and reset */
@@ -194,21 +235,21 @@ f_mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y)
/* square then multiply if the bit is set */
for (x = 0; x < bitcpy; x++) {
if ((err = mp_sqr (&res, &res)) != MP_OKAY) {
- goto __RES;
+ goto __RES;
}
if ((err = mp_reduce (&res, P, &mu)) != MP_OKAY) {
- goto __RES;
+ 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;
- }
+ /* 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;
+ }
}
}
}
View
83 bn_mp_exptmod_fast.c
@@ -19,7 +19,7 @@
* 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 or Diminished Radix reduction [whichever appropriate]
+ * Uses Montgomery or Diminished Radix reduction [whichever appropriate]
*/
int
mp_exptmod_fast (mp_int * G, mp_int * X, mp_int * P, mp_int * Y, int redmode)
@@ -28,7 +28,7 @@ mp_exptmod_fast (mp_int * G, mp_int * X, mp_int * P, mp_int * Y, int redmode)
mp_digit buf, mp;
int err, bitbuf, bitcpy, bitcnt, mode, digidx, x, y, winsize;
int (*redux)(mp_int*,mp_int*,mp_digit);
-
+
/* find window size */
x = mp_count_bits (X);
if (x <= 7) {
@@ -47,22 +47,37 @@ mp_exptmod_fast (mp_int * G, mp_int * X, mp_int * P, mp_int * Y, int redmode)
winsize = 8;
}
+#ifdef MP_LOW_MEM
+ if (winsize > 5) {
+ winsize = 5;
+ }
+#endif
+
+
/* init G array */
for (x = 0; x < (1 << winsize); x++) {
if ((err = mp_init (&M[x])) != MP_OKAY) {
for (y = 0; y < x; y++) {
- mp_clear (&M[y]);
+ mp_clear (&M[y]);
}
return err;
}
}
-
+
if (redmode == 0) {
/* now setup montgomery */
if ((err = mp_montgomery_setup (P, &mp)) != MP_OKAY) {
goto __M;
}
- redux = mp_montgomery_reduce;
+
+ /* automatically pick the comba one if available (saves quite a few calls/ifs) */
+ if ( ((P->used * 2 + 1) < MP_WARRAY) &&
+ P->used < (1 << ((CHAR_BIT * sizeof (mp_word)) - (2 * DIGIT_BIT)))) {
+ redux = fast_mp_montgomery_reduce;
+ } else {
+ /* use slower baselien method */
+ redux = mp_montgomery_reduce;
+ }
} else {
/* setup DR reduction */
mp_dr_setup(P, &mp);
@@ -97,7 +112,7 @@ mp_exptmod_fast (mp_int * G, mp_int * X, mp_int * P, mp_int * Y, int redmode)
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;
@@ -123,42 +138,42 @@ mp_exptmod_fast (mp_int * G, mp_int * X, mp_int * P, mp_int * Y, int redmode)
}
/* set initial mode and bit cnt */
- mode = 0;
- bitcnt = 0;
- buf = 0;
+ mode = 0;
+ bitcnt = 1;
+ buf = 0;
digidx = X->used - 1;
bitcpy = bitbuf = 0;
- bitcnt = 1;
for (;;) {
/* grab next digit as required */
if (--bitcnt == 0) {
if (digidx == -1) {
- break;
+ break;
}
buf = X->dp[digidx--];
bitcnt = (int) DIGIT_BIT;
}
/* grab the next msb from the exponent */
- y = (buf >> (DIGIT_BIT - 1)) & 1;
- buf <<= 1;
+ y = (mp_digit)(buf >> (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)
+ 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;
+ goto __RES;
}
if ((err = redux (&res, P, mp)) != MP_OKAY) {
- goto __RES;
+ goto __RES;
}
continue;
}
@@ -171,20 +186,20 @@ mp_exptmod_fast (mp_int * G, mp_int * X, mp_int * P, mp_int * Y, int redmode)
/* 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 __RES;
- }
- if ((err = redux (&res, P, mp)) != MP_OKAY) {
- goto __RES;
- }
+ if ((err = mp_sqr (&res, &res)) != MP_OKAY) {
+ goto __RES;
+ }
+ if ((err = redux (&res, P, mp)) != MP_OKAY) {
+ goto __RES;
+ }
}
/* then multiply */
if ((err = mp_mul (&res, &M[bitbuf], &res)) != MP_OKAY) {
- goto __RES;
+ goto __RES;
}
if ((err = redux (&res, P, mp)) != MP_OKAY) {
- goto __RES;
+ goto __RES;
}
/* empty window and reset */
@@ -198,21 +213,21 @@ mp_exptmod_fast (mp_int * G, mp_int * X, mp_int * P, mp_int * Y, int redmode)
/* square then multiply if the bit is set */
for (x = 0; x < bitcpy; x++) {
if ((err = mp_sqr (&res, &res)) != MP_OKAY) {
- goto __RES;
+ goto __RES;
}
if ((err = redux (&res, P, mp)) != MP_OKAY) {
- goto __RES;
+ 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 = redux (&res, P, mp)) != MP_OKAY) {
- goto __RES;
- }
+ /* then multiply */
+ if ((err = mp_mul (&res, &M[1], &res)) != MP_OKAY) {
+ goto __RES;
+ }
+ if ((err = redux (&res, P, mp)) != MP_OKAY) {
+ goto __RES;
+ }
}
}
}
@@ -222,7 +237,7 @@ mp_exptmod_fast (mp_int * G, mp_int * X, mp_int * P, mp_int * Y, int redmode)
if ((err = mp_montgomery_reduce (&res, P, mp)) != MP_OKAY) {
goto __RES;
}
- }
+ }
mp_exch (&res, Y);
err = MP_OKAY;
View
10 bn_mp_gcd.c
@@ -82,18 +82,18 @@ mp_gcd (mp_int * a, mp_int * b, mp_int * c)
/* B3 (and B4). Halve t, if even */
while (t.used != 0 && mp_iseven(&t) == 1) {
if ((res = mp_div_2 (&t, &t)) != MP_OKAY) {
- goto __T;
+ 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;
+ goto __T;
}
} else {
if ((res = mp_copy (&t, &v)) != MP_OKAY) {
- goto __T;
+ goto __T;
}
v.sign = (v.sign == MP_ZPOS) ? MP_NEG : MP_ZPOS;
}
@@ -102,9 +102,9 @@ mp_gcd (mp_int * a, mp_int * b, mp_int * c)
if ((res = mp_sub (&u, &v, &t)) != MP_OKAY) {
goto __T;
}
- }
- while (t.used != 0);
+ } while (mp_iszero(&t) == 0);
+ /* multiply by 2^k which we divided out at the beginning */
if ((res = mp_mul_2d (&u, k, &u)) != MP_OKAY) {
goto __T;
}
View
8 bn_mp_grow.c
@@ -18,12 +18,12 @@
int
mp_grow (mp_int * a, int size)
{
- int i, n;
+ int i;
/* if the alloc size is smaller alloc more ram */
if (a->alloc < size) {
/* ensure there are always at least MP_PREC digits extra on top */
- size += (MP_PREC * 2) - (size & (MP_PREC - 1));
+ size += (MP_PREC * 2) - (size & (MP_PREC - 1));
a->dp = OPT_CAST realloc (a->dp, sizeof (mp_digit) * size);
if (a->dp == NULL) {
@@ -31,9 +31,9 @@ mp_grow (mp_int * a, int size)
}
/* zero excess digits */
- n = a->alloc;
+ i = a->alloc;
a->alloc = size;
- for (i = n; i < a->alloc; i++) {
+ for (; i < a->alloc; i++) {
a->dp[i] = 0;
}
}
View
1  bn_mp_init.c
@@ -18,7 +18,6 @@
int
mp_init (mp_int * a)
{
-
/* allocate ram required and clear it */
a->dp = OPT_CAST calloc (sizeof (mp_digit), MP_PREC);
if (a->dp == NULL) {
View
91 bn_mp_invmod.c
@@ -29,63 +29,36 @@ mp_invmod (mp_int * a, mp_int * b, mp_int * c)
if (mp_iseven (b) == 0) {
return fast_mp_invmod (a, b, c);
}
-
- 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 (&A)) != MP_OKAY) {
- goto __V;
- }
-
- if ((res = mp_init (&B)) != MP_OKAY) {
- goto __A;
- }
-
- if ((res = mp_init (&C)) != MP_OKAY) {
- goto __B;
- }
-
- if ((res = mp_init (&D)) != MP_OKAY) {
- goto __C;
+
+ /* init temps */
+ if ((res = mp_init_multi(&x, &y, &u, &v, &A, &B, &C, &D, NULL)) != MP_OKAY) {
+ return res;
}
/* x = a, y = b */
if ((res = mp_copy (a, &x)) != MP_OKAY) {
- goto __D;
+ goto __ERR;
}
if ((res = mp_copy (b, &y)) != MP_OKAY) {
- goto __D;
+ goto __ERR;
}
if ((res = mp_abs (&x, &x)) != MP_OKAY) {
- goto __D;
+ goto __ERR;
}
/* 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;
+ goto __ERR;
}
/* 3. u=x, v=y, A=1, B=0, C=0,D=1 */
if ((res = mp_copy (&x, &u)) != MP_OKAY) {
- goto __D;
+ goto __ERR;
}
if ((res = mp_copy (&y, &v)) != MP_OKAY) {
- goto __D;
+ goto __ERR;
}
mp_set (&A, 1);
mp_set (&D, 1);
@@ -96,24 +69,24 @@ mp_invmod (mp_int * a, mp_int * b, mp_int * c)
while (mp_iseven (&u) == 1) {
/* 4.1 u = u/2 */
if ((res = mp_div_2 (&u, &u)) != MP_OKAY) {
- goto __D;
+ goto __ERR;
}
/* 4.2 if A or B is odd then */
if (mp_iseven (&A) == 0 || mp_iseven (&B) == 0) {
/* A = (A+y)/2, B = (B-x)/2 */
if ((res = mp_add (&A, &y, &A)) != MP_OKAY) {
- goto __D;
+ goto __ERR;
}
if ((res = mp_sub (&B, &x, &B)) != MP_OKAY) {
- goto __D;
+ goto __ERR;
}
}
/* A = A/2, B = B/2 */
if ((res = mp_div_2 (&A, &A)) != MP_OKAY) {
- goto __D;
+ goto __ERR;
}
if ((res = mp_div_2 (&B, &B)) != MP_OKAY) {
- goto __D;
+ goto __ERR;
}
}
@@ -122,24 +95,24 @@ mp_invmod (mp_int * a, mp_int * b, mp_int * c)
while (mp_iseven (&v) == 1) {
/* 5.1 v = v/2 */
if ((res = mp_div_2 (&v, &v)) != MP_OKAY) {
- goto __D;
+ goto __ERR;
}
/* 5.2 if C,D are even then */
if (mp_iseven (&C) == 0 || mp_iseven (&D) == 0) {
/* C = (C+y)/2, D = (D-x)/2 */
if ((res = mp_add (&C, &y, &C)) != MP_OKAY) {
- goto __D;
+ goto __ERR;
}
if ((res = mp_sub (&D, &x, &D)) != MP_OKAY) {
- goto __D;
+ goto __ERR;
}
}
/* C = C/2, D = D/2 */
if ((res = mp_div_2 (&C, &C)) != MP_OKAY) {
- goto __D;
+ goto __ERR;
}
if ((res = mp_div_2 (&D, &D)) != MP_OKAY) {
- goto __D;
+ goto __ERR;
}
}
@@ -147,28 +120,28 @@ mp_invmod (mp_int * a, mp_int * b, mp_int * c)
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;
+ goto __ERR;
}
if ((res = mp_sub (&A, &C, &A)) != MP_OKAY) {
- goto __D;
+ goto __ERR;
}
if ((res = mp_sub (&B, &D, &B)) != MP_OKAY) {
- goto __D;
+ goto __ERR;
}
} else {
/* v - v - u, C = C - A, D = D - B */
if ((res = mp_sub (&v, &u, &v)) != MP_OKAY) {
- goto __D;
+ goto __ERR;
}
if ((res = mp_sub (&C, &A, &C)) != MP_OKAY) {
- goto __D;
+ goto __ERR;
}
if ((res = mp_sub (&D, &B, &D)) != MP_OKAY) {
- goto __D;
+ goto __ERR;
}
}
@@ -181,21 +154,13 @@ mp_invmod (mp_int * a, mp_int * b, mp_int * c)
/* if v != 1 then there is no inverse */
if (mp_cmp_d (&v, 1) != MP_EQ) {
res = MP_VAL;
- goto __D;
+ goto __ERR;
}
/* a is now the inverse */
mp_exch (&C, c);
res = MP_OKAY;
-__D:mp_clear (&D);
-__C:mp_clear (&C);
-__B:mp_clear (&B);
-__A:mp_clear (&A);
-__V:mp_clear (&v);
-__U:mp_clear (&u);
-__Y:mp_clear (&y);
-__X:mp_clear (&x);
-__ERR:
+__ERR:mp_clear_multi (&x, &y, &u, &v, &A, &B, &C, &D, NULL);
return res;
}
View
2  bn_mp_jacobi.c
@@ -14,7 +14,7 @@
*/
#include <tommath.h>
-/* computes the jacobi c = (a | n) (or Legendre if b is prime)
+/* computes the jacobi c = (a | n) (or Legendre if n is prime)
* HAC pp. 73 Algorithm 2.149
*/
int
View
39 bn_mp_karatsuba_mul.c
@@ -36,7 +36,7 @@
int
mp_karatsuba_mul (mp_int * a, mp_int * b, mp_int * c)
{
- mp_int x0, x1, y0, y1, t1, t2, x0y0, x1y1;
+ mp_int x0, x1, y0, y1, t1, x0y0, x1y1;
int B, err;
err = MP_MEM;
@@ -60,10 +60,8 @@ mp_karatsuba_mul (mp_int * a, mp_int * b, mp_int * c)
/* init temps */
if (mp_init_size (&t1, B * 2) != MP_OKAY)
goto Y1;
- if (mp_init_size (&t2, B * 2) != MP_OKAY)
- goto T1;
if (mp_init_size (&x0y0, B * 2) != MP_OKAY)
- goto T2;
+ goto T1;
if (mp_init_size (&x1y1, B * 2) != MP_OKAY)
goto X0Y0;
@@ -110,41 +108,40 @@ mp_karatsuba_mul (mp_int * a, mp_int * b, mp_int * c)
mp_clamp (&y0);
/* now calc the products x0y0 and x1y1 */
- if (mp_mul (&x0, &y0, &x0y0) != MP_OKAY)
- goto X1Y1; /* x0y0 = x0*y0 */
+ if (mp_mul (&x0, &y0, &x0y0) != MP_OKAY) /* after this x0 is no longer required, free temp [x0==t2]! */
+ goto X1Y1; /* x0y0 = x0*y0 */
if (mp_mul (&x1, &y1, &x1y1) != MP_OKAY)
- goto X1Y1; /* x1y1 = x1*y1 */
+ goto X1Y1; /* x1y1 = x1*y1 */
/* now calc x1-x0 and y1-y0 */
if (mp_sub (&x1, &x0, &t1) != MP_OKAY)
- goto X1Y1; /* t1 = x1 - x0 */
- if (mp_sub (&y1, &y0, &t2) != MP_OKAY)
- goto X1Y1; /* t2 = y1 - y0 */
- if (mp_mul (&t1, &t2, &t1) != MP_OKAY)
- goto X1Y1; /* t1 = (x1 - x0) * (y1 - y0) */
+ goto X1Y1; /* t1 = x1 - x0 */
+ if (mp_sub (&y1, &y0, &x0) != MP_OKAY)
+ goto X1Y1; /* t2 = y1 - y0 */
+ if (mp_mul (&t1, &x0, &t1) != MP_OKAY)
+ goto X1Y1; /* t1 = (x1 - x0) * (y1 - y0) */
/* add x0y0 */
- if (mp_add (&x0y0, &x1y1, &t2) != MP_OKAY)
- goto X1Y1; /* t2 = x0y0 + x1y1 */
- if (mp_sub (&t2, &t1, &t1) != MP_OKAY)
- goto X1Y1; /* t1 = x0y0 + x1y1 - (x1-x0)*(y1-y0) */
+ if (mp_add (&x0y0, &x1y1, &x0) != MP_OKAY)
+ goto X1Y1; /* t2 = x0y0 + x1y1 */
+ if (mp_sub (&x0, &t1, &t1) != MP_OKAY)
+ goto X1Y1; /* t1 = x0y0 + x1y1 - (x1-x0)*(y1-y0) */
/* shift by B */
if (mp_lshd (&t1, B) != MP_OKAY)
- goto X1Y1; /* t1 = (x0y0 + x1y1 - (x1-x0)*(y1-y0))<<B */
+ goto X1Y1; /* t1 = (x0y0 + x1y1 - (x1-x0)*(y1-y0))<<B */
if (mp_lshd (&x1y1, B * 2) != MP_OKAY)
- goto X1Y1; /* x1y1 = x1y1 << 2*B */
+ goto X1Y1; /* x1y1 = x1y1 << 2*B */
if (mp_add (&x0y0, &t1, &t1) != MP_OKAY)
- goto X1Y1; /* t1 = x0y0 + t1 */
+ goto X1Y1; /* t1 = x0y0 + t1 */
if (mp_add (&t1, &x1y1, c) != MP_OKAY)
- goto X1Y1; /* t1 = x0y0 + t1 + x1y1 */
+ goto X1Y1; /* t1 = x0y0 + t1 + x1y1 */
err = MP_OKAY;
X1Y1:mp_clear (&x1y1);
X0Y0:mp_clear (&x0y0);
-T2:mp_clear (&t2);
T1:mp_clear (&t1);
Y1:mp_clear (&y1);
Y0:mp_clear (&y0);
View
22 bn_mp_karatsuba_sqr.c
@@ -74,32 +74,32 @@ mp_karatsuba_sqr (mp_int * a, mp_int * b)
/* now calc the products x0*x0 and x1*x1 */
if (mp_sqr (&x0, &x0x0) != MP_OKAY)
- goto X1X1; /* x0x0 = x0*x0 */
+ goto X1X1; /* x0x0 = x0*x0 */
if (mp_sqr (&x1, &x1x1) != MP_OKAY)
- goto X1X1; /* x1x1 = x1*x1 */
+ goto X1X1; /* x1x1 = x1*x1 */
- /* now calc x1-x0 and y1-y0 */
+ /* now calc (x1-x0)^2 */
if (mp_sub (&x1, &x0, &t1) != MP_OKAY)
- goto X1X1; /* t1 = x1 - x0 */
+ goto X1X1; /* t1 = x1 - x0 */
if (mp_sqr (&t1, &t1) != MP_OKAY)
- goto X1X1; /* t1 = (x1 - x0) * (y1 - y0) */
+ goto X1X1; /* t1 = (x1 - x0) * (x1 - x0) */
/* add x0y0 */
if (s_mp_add (&x0x0, &x1x1, &t2) != MP_OKAY)
- goto X1X1; /* t2 = x0y0 + x1y1 */
+ goto X1X1; /* t2 = x0y0 + x1y1 */
if (mp_sub (&t2, &t1, &t1) != MP_OKAY)
- goto X1X1; /* t1 = x0y0 + x1y1 - (x1-x0)*(y1-y0) */
+ goto X1X1; /* t1 = x0y0 + x1y1 - (x1-x0)*(y1-y0) */
/* shift by B */
if (mp_lshd (&t1, B) != MP_OKAY)
- goto X1X1; /* t1 = (x0y0 + x1y1 - (x1-x0)*(y1-y0))<<B */
+ goto X1X1; /* t1 = (x0y0 + x1y1 - (x1-x0)*(y1-y0))<<B */
if (mp_lshd (&x1x1, B * 2) != MP_OKAY)
- goto X1X1; /* x1y1 = x1y1 << 2*B */
+ goto X1X1; /* x1y1 = x1y1 << 2*B */
if (mp_add (&x0x0, &t1, &t1) != MP_OKAY)
- goto X1X1; /* t1 = x0y0 + t1 */
+ goto X1X1; /* t1 = x0y0 + t1 */
if (mp_add (&t1, &x1x1, b) != MP_OKAY)
- goto X1X1; /* t1 = x0y0 + t1 + x1y1 */
+ goto X1X1; /* t1 = x0y0 + t1 + x1y1 */
err = MP_OKAY;
View
7 bn_mp_lshd.c
@@ -20,15 +20,16 @@ mp_lshd (mp_int * a, int b)
{
int x, res;
-
/* if its less than zero return */
if (b <= 0) {
return MP_OKAY;
}
/* grow to fit the new digits */
- if ((res = mp_grow (a, a->used + b)) != MP_OKAY) {
- return res;
+ if (a->alloc < a->used + b) {
+ if ((res = mp_grow (a, a->used + b)) != MP_OKAY) {
+ return res;
+ }
}
{
View
8 bn_mp_montgomery_calc_normalization.c
@@ -15,10 +15,10 @@
#include <tommath.h>
/* calculates a = B^n mod b for Montgomery reduction
- * Where B is the base [e.g. 2^DIGIT_BIT].
+ * Where B is the base [e.g. 2^DIGIT_BIT].
* B^n mod b is computed by first computing
* A = B^(n-1) which doesn't require a reduction but a simple OR.
- * then C = A * B = B^n is computed by performing upto DIGIT_BIT
+ * then C = A * B = B^n is computed by performing upto DIGIT_BIT
* shifts with subtractions when the result is greater than b.
*
* The method is slightly modified to shift B unconditionally upto just under
@@ -38,13 +38,13 @@ mp_montgomery_calc_normalization (mp_int * a, mp_int * b)
}
/* now compute C = A * B mod b */
- for (x = bits - 1; x < DIGIT_BIT; x++) {
+ for (x = bits - 1; x < (int)DIGIT_BIT; x++) {
if ((res = mp_mul_2 (a, a)) != MP_OKAY) {
return res;
}
if (mp_cmp_mag (a, b) != MP_LT) {
if ((res = s_mp_sub (a, b, a)) != MP_OKAY) {
- return res;
+ return res;
}
}
}
View
23 bn_mp_montgomery_reduce.c
@@ -21,12 +21,19 @@ mp_montgomery_reduce (mp_int * a, mp_int * m, mp_digit mp)
int ix, res, digs;
mp_digit ui;
+ /* can the fast reduction [comba] method be used?
+ *
+ * Note that unlike in mp_mul you're safely allowed *less*
+ * than the available columns [255 per default] since carries
+ * are fixed up in the inner loop.
+ */
digs = m->used * 2 + 1;
- if ((digs < 512)
- && digs < (1 << ((CHAR_BIT * sizeof (mp_word)) - (2 * DIGIT_BIT)))) {
+ if ((digs < MP_WARRAY)
+ && m->used < (1 << ((CHAR_BIT * sizeof (mp_word)) - (2 * DIGIT_BIT)))) {
return fast_mp_montgomery_reduce (a, m, mp);
}
+ /* grow the input as required */
if (a->alloc < m->used * 2 + 1) {
if ((res = mp_grow (a, m->used * 2 + 1)) != MP_OKAY) {
return res;
@@ -50,15 +57,15 @@ mp_montgomery_reduce (mp_int * a, mp_int * m, mp_digit mp)
mu = 0;
for (iy = 0; iy < m->used; iy++) {
- r = ((mp_word) ui) * ((mp_word) * tmpx++) + ((mp_word) mu) + ((mp_word) * tmpy);
- mu = (r >> ((mp_word) DIGIT_BIT));
- *tmpy++ = (r & ((mp_word) MP_MASK));
+ r = ((mp_word) ui) * ((mp_word) * tmpx++) + ((mp_word) mu) + ((mp_word) * tmpy);
+ mu = (r >> ((mp_word) DIGIT_BIT));
+ *tmpy++ = (r & ((mp_word) MP_MASK));
}
/* propagate carries */
while (mu) {
- *tmpy += mu;
- mu = (*tmpy >> DIGIT_BIT) & 1;
- *tmpy++ &= MP_MASK;
+ *tmpy += mu;
+ mu = (*tmpy >> DIGIT_BIT) & 1;
+ *tmpy++ &= MP_MASK;
}
}
}
View
23 bn_mp_montgomery_setup.c
@@ -18,11 +18,11 @@
int
mp_montgomery_setup (mp_int * a, mp_digit * mp)
{
- unsigned long x, b;
+ mp_digit x, b;
-/* fast inversion mod 2^32
+/* fast inversion mod 2^k
*
- * Based on the fact that
+ * Based on the fact that
*
* XA = 1 (mod 2^n) => (X(2-XA)) A = 1 (mod 2^2n)
* => 2*X*A - X*X*A*A = 1
@@ -34,13 +34,20 @@ mp_montgomery_setup (mp_int * a, mp_digit * mp)
return MP_VAL;
}
- x = (((b + 2) & 4) << 1) + b; /* here x*a==1 mod 2^4 */
- x *= 2 - b * x; /* here x*a==1 mod 2^8 */
- x *= 2 - b * x; /* here x*a==1 mod 2^16; each step doubles the nb of bits */
- x *= 2 - b * x; /* here x*a==1 mod 2^32 */
+ x = (((b + 2) & 4) << 1) + b; /* here x*a==1 mod 2^4 */
+ x *= 2 - b * x; /* here x*a==1 mod 2^8 */
+#if !defined(MP_8BIT)
+ x *= 2 - b * x; /* here x*a==1 mod 2^16; each step doubles the nb of bits */
+#endif
+#if defined(MP_64BIT) || !(defined(MP_8BIT) || defined(MP_16BIT))
+ x *= 2 - b * x; /* here x*a==1 mod 2^32 */
+#endif
+#ifdef MP_64BIT
+ x *= 2 - b * x; /* here x*a==1 mod 2^64 */
+#endif
/* t = -1/m mod b */
- *mp = ((mp_digit) 1 << ((mp_digit) DIGIT_BIT)) - (x & MP_MASK);
+ *mp = (((mp_digit) 1 << ((mp_digit) DIGIT_BIT)) - x) & MP_MASK;
return MP_OKAY;
}
View
10 bn_mp_mul.c
@@ -24,15 +24,15 @@ mp_mul (mp_int * a, mp_int * b, mp_int * c)
res = mp_karatsuba_mul (a, b, c);
} else {
- /* can we use the fast multiplier?
+ /* can we use the fast multiplier?
*
- * The fast multiplier can be used if the output will have less than
- * 512 digits and the number of digits won't affect carry propagation
+ * The fast multiplier can be used if the output will have less than
+ * MP_WARRAY digits and the number of digits won't affect carry propagation
*/
int digs = a->used + b->used + 1;
- if ((digs < 512)
- && digs < (1 << ((CHAR_BIT * sizeof (mp_word)) - (2 * DIGIT_BIT)))) {
+ if ((digs < MP_WARRAY)
+ && MIN(a->used, b->used) <= (1 << ((CHAR_BIT * sizeof (mp_word)) - (2 * DIGIT_BIT)))) {
res = fast_s_mp_mul_digs (a, b, c, digs);
} else {
res = s_mp_mul (a, b, c);
View
37 bn_mp_mul_2.c
@@ -20,10 +20,9 @@ mp_mul_2 (mp_int * a, mp_int * b)
{
int x, res, oldused;
- /* Optimization: should copy and shift at the same time */
-
- if (b->alloc < a->used) {
- if ((res = mp_grow (b, a->used)) != MP_OKAY) {
+ /* grow to accomodate result */
+ if (b->alloc < a->used + 1) {
+ if ((res = mp_grow (b, a->used + 1)) != MP_OKAY) {
return res;
}
}
@@ -31,7 +30,6 @@ mp_mul_2 (mp_int * a, mp_int * b)
oldused = b->used;
b->used = a->used;
- /* shift any bit count < DIGIT_BIT */
{
register mp_digit r, rr, *tmpa, *tmpb;
@@ -43,37 +41,32 @@ mp_mul_2 (mp_int * a, mp_int * b)
/* carry */
r = 0;
- for (x = 0; x < b->used; x++) {
+ for (x = 0; x < a->used; x++) {
- /* get what will be the *next* carry bit from the MSB of the current digit */
- rr = *tmpa >> (DIGIT_BIT - 1);
+ /* get what will be the *next* carry bit from the
+ * MSB of the current digit
+ */
+ rr = *tmpa >> ((mp_digit)(DIGIT_BIT - 1));
/* now shift up this digit, add in the carry [from the previous] */
- *tmpb++ = ((*tmpa++ << 1) | r) & MP_MASK;
+ *tmpb++ = ((*tmpa++ << ((mp_digit)1)) | r) & MP_MASK;
- /* copy the carry that would be from the source digit into the next iteration */
+ /* copy the carry that would be from the source
+ * digit into the next iteration
+ */
r = rr;
}
/* new leading digit? */
if (r != 0) {
- /* do we have to grow to accomodate the new digit? */
- if (b->alloc == b->used) {
- if ((res = mp_grow (b, b->used + 1)) != MP_OKAY) {
- return res;
- }
-
- /* after the grow *tmpb is no longer valid so we have to reset it!
- * (this bug took me about 17 minutes to find...!)
- */
- tmpb = b->dp + b->used;
- }
/* add a MSB which is always 1 at this point */
*tmpb = 1;
++b->used;
}
- /* now zero any excess digits on the destination that we didn't write to */
+ /* now zero any excess digits on the destination
+ * that we didn't write to
+ */
tmpb = b->dp + b->used;
for (x = b->used; x < oldused; x++) {
*tmpb++ = 0;
View
35 bn_mp_mul_2d.c
@@ -14,24 +14,34 @@
*/
#include <tommath.h>
+/* NOTE: This routine requires updating. For instance the c->used = c->alloc bit
+ is wrong. We should just shift c->used digits then set the carry as c->dp[c->used] = carry
+
+ To be fixed for LTM 0.18
+ */
+
/* shift left by a certain bit count */
int
mp_mul_2d (mp_int * a, int b, mp_int * c)
{
- mp_digit d, r, rr;
- int x, res;
+ mp_digit d;
+ int res;
/* copy */
- if ((res = mp_copy (a, c)) != MP_OKAY) {
- return res;
+ if (a != c) {
+ if ((res = mp_copy (a, c)) != MP_OKAY) {
+ return res;
+ }
}
- if ((res = mp_grow (c, c->used + b / DIGIT_BIT + 1)) != MP_OKAY) {
- return res;
+ if (c->alloc < (int)(c->used + b/DIGIT_BIT + 2)) {
+ if ((res = mp_grow (c, c->used + b / DIGIT_BIT + 2)) != MP_OKAY) {
+ return res;
+ }
}
/* shift by as many digits in the bit count */
- if (b >= DIGIT_BIT) {
+ if (b >= (int)DIGIT_BIT) {
if ((res = mp_lshd (c, b / DIGIT_BIT)) != MP_OKAY) {
return res;
}
@@ -41,14 +51,15 @@ mp_mul_2d (mp_int * a, int b, mp_int * c)
/* shift any bit count < DIGIT_BIT */
d = (mp_digit) (b % DIGIT_BIT);
if (d != 0) {
- register mp_digit *tmpc, mask;
-
+ register mp_digit *tmpc, mask, r, rr;
+ register int x;
+
/* bitmask for carries */
- mask = (1U << d) - 1U;
-
+ mask = (((mp_digit)1) << d) - 1;
+
/* alias */
tmpc = c->dp;
-
+
/* carry */
r = 0;
for (x = 0; x < c->used; x++) {
View
26 bn_mp_mul_d.c
@@ -20,6 +20,7 @@ mp_mul_d (mp_int * a, mp_digit b, mp_int * c)
{
int res, pa, olduse;
+ /* make sure c is big enough to hold a*b */
pa = a->used;
if (c->alloc < pa + 1) {
if ((res = mp_grow (c, pa + 1)) != MP_OKAY) {
@@ -27,7 +28,10 @@ mp_mul_d (mp_int * a, mp_digit b, mp_int * c)
}
}
+ /* get the original destinations used count */
olduse = c->used;
+
+ /* set the new temporary used count */
c->used = pa + 1;
{
@@ -35,21 +39,31 @@ mp_mul_d (mp_int * a, mp_digit b, mp_int * c)
register mp_word r;
register int ix;
- tmpc = c->dp + c->used;
- for (ix = c->used; ix < olduse; ix++) {
- *tmpc++ = 0;
- }
-
+ /* alias for a->dp [source] */
tmpa = a->dp;
+
+ /* alias for c->dp [dest] */
tmpc = c->dp;
+ /* zero carry */
u = 0;
for (ix = 0; ix < pa; ix++) {
+ /* compute product and carry sum for this term */
r = ((mp_word) u) + ((mp_word) * tmpa++) * ((mp_word) b);
+
+ /* mask off higher bits to get a single digit */
*tmpc++ = (mp_digit) (r & ((mp_word) MP_MASK));
+
+ /* send carry into next iteration */
u = (mp_digit) (r >> ((mp_word) DIGIT_BIT));
}
- *tmpc = u;
+ /* store final carry [if any] */
+ *tmpc++ = u;
+
+ /* now zero digits above the top */
+ for (; pa < olduse; pa++) {
+ *tmpc++ = 0;
+ }
}
mp_clamp (c);
View
64 bn_mp_multi.c
@@ -0,0 +1,64 @@
+/* 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://math.libtomcrypt.org
+ */
+#include <tommath.h>
+#include <stdarg.h>
+
+int mp_init_multi(mp_int *mp, ...)
+{
+ mp_err res = MP_OKAY; /* Assume ok until proven otherwise */
+ int n = 0; /* Number of ok inits */
+ mp_int* cur_arg = mp;
+ va_list args;
+
+ va_start(args, mp); /* init args to next argument from caller */
+ while (cur_arg != NULL) {
+ if (mp_init(cur_arg) != MP_OKAY) {
+ /* Oops - error! Back-track and mp_clear what we already
+ succeeded in init-ing, then return error.
+ */
+ va_list clean_args;
+
+ /* end the current list */
+ va_end(args);
+
+ /* now start cleaning up */