1 /* Start: bn_error.c */
4 /* LibTomMath, multiple-precision integer library -- Tom St Denis
6 * LibTomMath is a library that provides multiple-precision
7 * integer arithmetic as well as number theoretic functionality.
9 * The library was designed directly after the MPI library by
10 * Michael Fromberger but has been written from scratch with
11 * additional optimizations in place.
13 * The library is free for all purposes without any express
16 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
23 { MP_OKAY, "Successful" },
24 { MP_MEM, "Out of heap" },
25 { MP_VAL, "Value out of range" }
28 /* return a char * string for a given code */
29 char *mp_error_to_string(int code)
33 /* scan the lookup table for the given message */
34 for (x = 0; x < (int)(sizeof(msgs) / sizeof(msgs[0])); x++) {
35 if (msgs[x].code == code) {
40 /* generic reply for invalid code */
41 return "Invalid error code";
48 /* Start: bn_fast_mp_invmod.c */
50 #ifdef BN_FAST_MP_INVMOD_C
51 /* LibTomMath, multiple-precision integer library -- Tom St Denis
53 * LibTomMath is a library that provides multiple-precision
54 * integer arithmetic as well as number theoretic functionality.
56 * The library was designed directly after the MPI library by
57 * Michael Fromberger but has been written from scratch with
58 * additional optimizations in place.
60 * The library is free for all purposes without any express
63 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
66 /* computes the modular inverse via binary extended euclidean algorithm,
67 * that is c = 1/a mod b
69 * Based on slow invmod except this is optimized for the case where b is
70 * odd as per HAC Note 14.64 on pp. 610
72 int fast_mp_invmod (mp_int * a, mp_int * b, mp_int * c)
74 mp_int x, y, u, v, B, D;
77 /* 2. [modified] b must be odd */
78 if (mp_iseven (b) == 1) {
82 /* init all our temps */
83 if ((res = mp_init_multi(&x, &y, &u, &v, &B, &D, NULL)) != MP_OKAY) {
87 /* x == modulus, y == value to invert */
88 if ((res = mp_copy (b, &x)) != MP_OKAY) {
93 if ((res = mp_mod (a, b, &y)) != MP_OKAY) {
97 /* 3. u=x, v=y, A=1, B=0, C=0,D=1 */
98 if ((res = mp_copy (&x, &u)) != MP_OKAY) {
101 if ((res = mp_copy (&y, &v)) != MP_OKAY) {
107 /* 4. while u is even do */
108 while (mp_iseven (&u) == 1) {
110 if ((res = mp_div_2 (&u, &u)) != MP_OKAY) {
113 /* 4.2 if B is odd then */
114 if (mp_isodd (&B) == 1) {
115 if ((res = mp_sub (&B, &x, &B)) != MP_OKAY) {
120 if ((res = mp_div_2 (&B, &B)) != MP_OKAY) {
125 /* 5. while v is even do */
126 while (mp_iseven (&v) == 1) {
128 if ((res = mp_div_2 (&v, &v)) != MP_OKAY) {
131 /* 5.2 if D is odd then */
132 if (mp_isodd (&D) == 1) {
134 if ((res = mp_sub (&D, &x, &D)) != MP_OKAY) {
139 if ((res = mp_div_2 (&D, &D)) != MP_OKAY) {
144 /* 6. if u >= v then */
145 if (mp_cmp (&u, &v) != MP_LT) {
146 /* u = u - v, B = B - D */
147 if ((res = mp_sub (&u, &v, &u)) != MP_OKAY) {
151 if ((res = mp_sub (&B, &D, &B)) != MP_OKAY) {
155 /* v - v - u, D = D - B */
156 if ((res = mp_sub (&v, &u, &v)) != MP_OKAY) {
160 if ((res = mp_sub (&D, &B, &D)) != MP_OKAY) {
165 /* if not zero goto step 4 */
166 if (mp_iszero (&u) == 0) {
170 /* now a = C, b = D, gcd == g*v */
172 /* if v != 1 then there is no inverse */
173 if (mp_cmp_d (&v, 1) != MP_EQ) {
178 /* b is now the inverse */
180 while (D.sign == MP_NEG) {
181 if ((res = mp_add (&D, b, &D)) != MP_OKAY) {
189 LBL_ERR:mp_clear_multi (&x, &y, &u, &v, &B, &D, NULL);
194 /* End: bn_fast_mp_invmod.c */
196 /* Start: bn_fast_mp_montgomery_reduce.c */
198 #ifdef BN_FAST_MP_MONTGOMERY_REDUCE_C
199 /* LibTomMath, multiple-precision integer library -- Tom St Denis
201 * LibTomMath is a library that provides multiple-precision
202 * integer arithmetic as well as number theoretic functionality.
204 * The library was designed directly after the MPI library by
205 * Michael Fromberger but has been written from scratch with
206 * additional optimizations in place.
208 * The library is free for all purposes without any express
209 * guarantee it works.
211 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
214 /* computes xR**-1 == x (mod N) via Montgomery Reduction
216 * This is an optimized implementation of montgomery_reduce
217 * which uses the comba method to quickly calculate the columns of the
220 * Based on Algorithm 14.32 on pp.601 of HAC.
222 int fast_mp_montgomery_reduce (mp_int * x, mp_int * n, mp_digit rho)
225 mp_word W[MP_WARRAY];
227 /* get old used count */
230 /* grow a as required */
231 if (x->alloc < n->used + 1) {
232 if ((res = mp_grow (x, n->used + 1)) != MP_OKAY) {
237 /* first we have to get the digits of the input into
238 * an array of double precision words W[...]
241 register mp_word *_W;
242 register mp_digit *tmpx;
244 /* alias for the W[] array */
247 /* alias for the digits of x*/
250 /* copy the digits of a into W[0..a->used-1] */
251 for (ix = 0; ix < x->used; ix++) {
255 /* zero the high words of W[a->used..m->used*2] */
256 for (; ix < n->used * 2 + 1; ix++) {
261 /* now we proceed to zero successive digits
262 * from the least significant upwards
264 for (ix = 0; ix < n->used; ix++) {
265 /* mu = ai * m' mod b
267 * We avoid a double precision multiplication (which isn't required)
268 * by casting the value down to a mp_digit. Note this requires
269 * that W[ix-1] have the carry cleared (see after the inner loop)
271 register mp_digit mu;
272 mu = (mp_digit) (((W[ix] & MP_MASK) * rho) & MP_MASK);
274 /* a = a + mu * m * b**i
276 * This is computed in place and on the fly. The multiplication
277 * by b**i is handled by offseting which columns the results
280 * Note the comba method normally doesn't handle carries in the
281 * inner loop In this case we fix the carry from the previous
282 * column since the Montgomery reduction requires digits of the
283 * result (so far) [see above] to work. This is
284 * handled by fixing up one carry after the inner loop. The
285 * carry fixups are done in order so after these loops the
286 * first m->used words of W[] have the carries fixed
290 register mp_digit *tmpn;
291 register mp_word *_W;
293 /* alias for the digits of the modulus */
296 /* Alias for the columns set by an offset of ix */
300 for (iy = 0; iy < n->used; iy++) {
301 *_W++ += ((mp_word)mu) * ((mp_word)*tmpn++);
305 /* now fix carry for next digit, W[ix+1] */
306 W[ix + 1] += W[ix] >> ((mp_word) DIGIT_BIT);
309 /* now we have to propagate the carries and
310 * shift the words downward [all those least
311 * significant digits we zeroed].
314 register mp_digit *tmpx;
315 register mp_word *_W, *_W1;
317 /* nox fix rest of carries */
319 /* alias for current word */
322 /* alias for next word, where the carry goes */
325 for (; ix <= n->used * 2 + 1; ix++) {
326 *_W++ += *_W1++ >> ((mp_word) DIGIT_BIT);
329 /* copy out, A = A/b**n
331 * The result is A/b**n but instead of converting from an
332 * array of mp_word to mp_digit than calling mp_rshd
333 * we just copy them in the right order
336 /* alias for destination word */
339 /* alias for shifted double precision result */
342 for (ix = 0; ix < n->used + 1; ix++) {
343 *tmpx++ = (mp_digit)(*_W++ & ((mp_word) MP_MASK));
346 /* zero oldused digits, if the input a was larger than
347 * m->used+1 we'll have to clear the digits
349 for (; ix < olduse; ix++) {
354 /* set the max used and clamp */
355 x->used = n->used + 1;
358 /* if A >= m then A = A - m */
359 if (mp_cmp_mag (x, n) != MP_LT) {
360 return s_mp_sub (x, n, x);
366 /* End: bn_fast_mp_montgomery_reduce.c */
368 /* Start: bn_fast_s_mp_mul_digs.c */
370 #ifdef BN_FAST_S_MP_MUL_DIGS_C
371 /* LibTomMath, multiple-precision integer library -- Tom St Denis
373 * LibTomMath is a library that provides multiple-precision
374 * integer arithmetic as well as number theoretic functionality.
376 * The library was designed directly after the MPI library by
377 * Michael Fromberger but has been written from scratch with
378 * additional optimizations in place.
380 * The library is free for all purposes without any express
381 * guarantee it works.
383 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
386 /* Fast (comba) multiplier
388 * This is the fast column-array [comba] multiplier. It is
389 * designed to compute the columns of the product first
390 * then handle the carries afterwards. This has the effect
391 * of making the nested loops that compute the columns very
392 * simple and schedulable on super-scalar processors.
394 * This has been modified to produce a variable number of
395 * digits of output so if say only a half-product is required
396 * you don't have to compute the upper half (a feature
397 * required for fast Barrett reduction).
399 * Based on Algorithm 14.12 on pp.595 of HAC.
402 int fast_s_mp_mul_digs (mp_int * a, mp_int * b, mp_int * c, int digs)
404 int olduse, res, pa, ix, iz;
405 mp_digit W[MP_WARRAY];
408 /* grow the destination as required */
409 if (c->alloc < digs) {
410 if ((res = mp_grow (c, digs)) != MP_OKAY) {
415 /* number of output digits to produce */
416 pa = MIN(digs, a->used + b->used);
418 /* clear the carry */
420 for (ix = 0; ix < pa; ix++) {
423 mp_digit *tmpx, *tmpy;
425 /* get offsets into the two bignums */
426 ty = MIN(b->used-1, ix);
429 /* setup temp aliases */
433 /* this is the number of times the loop will iterrate, essentially
434 while (tx++ < a->used && ty-- >= 0) { ... }
436 iy = MIN(a->used-tx, ty+1);
439 for (iz = 0; iz < iy; ++iz) {
440 _W += ((mp_word)*tmpx++)*((mp_word)*tmpy--);
444 W[ix] = ((mp_digit)_W) & MP_MASK;
446 /* make next carry */
447 _W = _W >> ((mp_word)DIGIT_BIT);
450 /* store final carry */
451 W[ix] = (mp_digit)(_W & MP_MASK);
458 register mp_digit *tmpc;
460 for (ix = 0; ix < pa+1; ix++) {
461 /* now extract the previous digit [below the carry] */
465 /* clear unused digits [that existed in the old copy of c] */
466 for (; ix < olduse; ix++) {
475 /* End: bn_fast_s_mp_mul_digs.c */
477 /* Start: bn_fast_s_mp_mul_high_digs.c */
479 #ifdef BN_FAST_S_MP_MUL_HIGH_DIGS_C
480 /* LibTomMath, multiple-precision integer library -- Tom St Denis
482 * LibTomMath is a library that provides multiple-precision
483 * integer arithmetic as well as number theoretic functionality.
485 * The library was designed directly after the MPI library by
486 * Michael Fromberger but has been written from scratch with
487 * additional optimizations in place.
489 * The library is free for all purposes without any express
490 * guarantee it works.
492 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
495 /* this is a modified version of fast_s_mul_digs that only produces
496 * output digits *above* digs. See the comments for fast_s_mul_digs
497 * to see how it works.
499 * This is used in the Barrett reduction since for one of the multiplications
500 * only the higher digits were needed. This essentially halves the work.
502 * Based on Algorithm 14.12 on pp.595 of HAC.
504 int fast_s_mp_mul_high_digs (mp_int * a, mp_int * b, mp_int * c, int digs)
506 int olduse, res, pa, ix, iz;
507 mp_digit W[MP_WARRAY];
510 /* grow the destination as required */
511 pa = a->used + b->used;
513 if ((res = mp_grow (c, pa)) != MP_OKAY) {
518 /* number of output digits to produce */
519 pa = a->used + b->used;
521 for (ix = digs; ix < pa; ix++) {
523 mp_digit *tmpx, *tmpy;
525 /* get offsets into the two bignums */
526 ty = MIN(b->used-1, ix);
529 /* setup temp aliases */
533 /* this is the number of times the loop will iterrate, essentially its
534 while (tx++ < a->used && ty-- >= 0) { ... }
536 iy = MIN(a->used-tx, ty+1);
539 for (iz = 0; iz < iy; iz++) {
540 _W += ((mp_word)*tmpx++)*((mp_word)*tmpy--);
544 W[ix] = ((mp_digit)_W) & MP_MASK;
546 /* make next carry */
547 _W = _W >> ((mp_word)DIGIT_BIT);
550 /* store final carry */
551 W[ix] = (mp_digit)(_W & MP_MASK);
558 register mp_digit *tmpc;
561 for (ix = digs; ix <= pa; ix++) {
562 /* now extract the previous digit [below the carry] */
566 /* clear unused digits [that existed in the old copy of c] */
567 for (; ix < olduse; ix++) {
576 /* End: bn_fast_s_mp_mul_high_digs.c */
578 /* Start: bn_fast_s_mp_sqr.c */
580 #ifdef BN_FAST_S_MP_SQR_C
581 /* LibTomMath, multiple-precision integer library -- Tom St Denis
583 * LibTomMath is a library that provides multiple-precision
584 * integer arithmetic as well as number theoretic functionality.
586 * The library was designed directly after the MPI library by
587 * Michael Fromberger but has been written from scratch with
588 * additional optimizations in place.
590 * The library is free for all purposes without any express
591 * guarantee it works.
593 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
596 /* the jist of squaring...
597 * you do like mult except the offset of the tmpx [one that
598 * starts closer to zero] can't equal the offset of tmpy.
599 * So basically you set up iy like before then you min it with
600 * (ty-tx) so that it never happens. You double all those
601 * you add in the inner loop
603 After that loop you do the squares and add them in.
606 int fast_s_mp_sqr (mp_int * a, mp_int * b)
608 int olduse, res, pa, ix, iz;
609 mp_digit W[MP_WARRAY], *tmpx;
612 /* grow the destination as required */
613 pa = a->used + a->used;
615 if ((res = mp_grow (b, pa)) != MP_OKAY) {
620 /* number of output digits to produce */
622 for (ix = 0; ix < pa; ix++) {
630 /* get offsets into the two bignums */
631 ty = MIN(a->used-1, ix);
634 /* setup temp aliases */
638 /* this is the number of times the loop will iterrate, essentially
639 while (tx++ < a->used && ty-- >= 0) { ... }
641 iy = MIN(a->used-tx, ty+1);
643 /* now for squaring tx can never equal ty
644 * we halve the distance since they approach at a rate of 2x
645 * and we have to round because odd cases need to be executed
647 iy = MIN(iy, (ty-tx+1)>>1);
650 for (iz = 0; iz < iy; iz++) {
651 _W += ((mp_word)*tmpx++)*((mp_word)*tmpy--);
654 /* double the inner product and add carry */
657 /* even columns have the square term in them */
659 _W += ((mp_word)a->dp[ix>>1])*((mp_word)a->dp[ix>>1]);
663 W[ix] = (mp_digit)(_W & MP_MASK);
665 /* make next carry */
666 W1 = _W >> ((mp_word)DIGIT_BIT);
671 b->used = a->used+a->used;
676 for (ix = 0; ix < pa; ix++) {
677 *tmpb++ = W[ix] & MP_MASK;
680 /* clear unused digits [that existed in the old copy of c] */
681 for (; ix < olduse; ix++) {
690 /* End: bn_fast_s_mp_sqr.c */
692 /* Start: bn_mp_2expt.c */
695 /* LibTomMath, multiple-precision integer library -- Tom St Denis
697 * LibTomMath is a library that provides multiple-precision
698 * integer arithmetic as well as number theoretic functionality.
700 * The library was designed directly after the MPI library by
701 * Michael Fromberger but has been written from scratch with
702 * additional optimizations in place.
704 * The library is free for all purposes without any express
705 * guarantee it works.
707 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
712 * Simple algorithm which zeroes the int, grows it then just sets one bit
716 mp_2expt (mp_int * a, int b)
720 /* zero a as per default */
723 /* grow a to accomodate the single bit */
724 if ((res = mp_grow (a, b / DIGIT_BIT + 1)) != MP_OKAY) {
728 /* set the used count of where the bit will go */
729 a->used = b / DIGIT_BIT + 1;
731 /* put the single bit in its place */
732 a->dp[b / DIGIT_BIT] = ((mp_digit)1) << (b % DIGIT_BIT);
738 /* End: bn_mp_2expt.c */
740 /* Start: bn_mp_abs.c */
743 /* LibTomMath, multiple-precision integer library -- Tom St Denis
745 * LibTomMath is a library that provides multiple-precision
746 * integer arithmetic as well as number theoretic functionality.
748 * The library was designed directly after the MPI library by
749 * Michael Fromberger but has been written from scratch with
750 * additional optimizations in place.
752 * The library is free for all purposes without any express
753 * guarantee it works.
755 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
760 * Simple function copies the input and fixes the sign to positive
763 mp_abs (mp_int * a, mp_int * b)
769 if ((res = mp_copy (a, b)) != MP_OKAY) {
774 /* force the sign of b to positive */
781 /* End: bn_mp_abs.c */
783 /* Start: bn_mp_add.c */
786 /* LibTomMath, multiple-precision integer library -- Tom St Denis
788 * LibTomMath is a library that provides multiple-precision
789 * integer arithmetic as well as number theoretic functionality.
791 * The library was designed directly after the MPI library by
792 * Michael Fromberger but has been written from scratch with
793 * additional optimizations in place.
795 * The library is free for all purposes without any express
796 * guarantee it works.
798 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
801 /* high level addition (handles signs) */
802 int mp_add (mp_int * a, mp_int * b, mp_int * c)
806 /* get sign of both inputs */
810 /* handle two cases, not four */
812 /* both positive or both negative */
813 /* add their magnitudes, copy the sign */
815 res = s_mp_add (a, b, c);
817 /* one positive, the other negative */
818 /* subtract the one with the greater magnitude from */
819 /* the one of the lesser magnitude. The result gets */
820 /* the sign of the one with the greater magnitude. */
821 if (mp_cmp_mag (a, b) == MP_LT) {
823 res = s_mp_sub (b, a, c);
826 res = s_mp_sub (a, b, c);
834 /* End: bn_mp_add.c */
836 /* Start: bn_mp_add_d.c */
839 /* LibTomMath, multiple-precision integer library -- Tom St Denis
841 * LibTomMath is a library that provides multiple-precision
842 * integer arithmetic as well as number theoretic functionality.
844 * The library was designed directly after the MPI library by
845 * Michael Fromberger but has been written from scratch with
846 * additional optimizations in place.
848 * The library is free for all purposes without any express
849 * guarantee it works.
851 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
854 /* single digit addition */
856 mp_add_d (mp_int * a, mp_digit b, mp_int * c)
858 int res, ix, oldused;
859 mp_digit *tmpa, *tmpc, mu;
861 /* grow c as required */
862 if (c->alloc < a->used + 1) {
863 if ((res = mp_grow(c, a->used + 1)) != MP_OKAY) {
868 /* if a is negative and |a| >= b, call c = |a| - b */
869 if (a->sign == MP_NEG && (a->used > 1 || a->dp[0] >= b)) {
870 /* temporarily fix sign of a */
874 res = mp_sub_d(a, b, c);
877 a->sign = c->sign = MP_NEG;
882 /* old number of used digits in c */
885 /* sign always positive */
891 /* destination alias */
894 /* if a is positive */
895 if (a->sign == MP_ZPOS) {
896 /* add digit, after this we're propagating
900 mu = *tmpc >> DIGIT_BIT;
903 /* now handle rest of the digits */
904 for (ix = 1; ix < a->used; ix++) {
905 *tmpc = *tmpa++ + mu;
906 mu = *tmpc >> DIGIT_BIT;
909 /* set final carry */
914 c->used = a->used + 1;
916 /* a was negative and |a| < b */
919 /* the result is a single digit */
921 *tmpc++ = b - a->dp[0];
926 /* setup count so the clearing of oldused
927 * can fall through correctly
932 /* now zero to oldused */
933 while (ix++ < oldused) {
943 /* End: bn_mp_add_d.c */
945 /* Start: bn_mp_addmod.c */
947 #ifdef BN_MP_ADDMOD_C
948 /* LibTomMath, multiple-precision integer library -- Tom St Denis
950 * LibTomMath is a library that provides multiple-precision
951 * integer arithmetic as well as number theoretic functionality.
953 * The library was designed directly after the MPI library by
954 * Michael Fromberger but has been written from scratch with
955 * additional optimizations in place.
957 * The library is free for all purposes without any express
958 * guarantee it works.
960 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
963 /* d = a + b (mod c) */
965 mp_addmod (mp_int * a, mp_int * b, mp_int * c, mp_int * d)
970 if ((res = mp_init (&t)) != MP_OKAY) {
974 if ((res = mp_add (a, b, &t)) != MP_OKAY) {
978 res = mp_mod (&t, c, d);
984 /* End: bn_mp_addmod.c */
986 /* Start: bn_mp_and.c */
989 /* LibTomMath, multiple-precision integer library -- Tom St Denis
991 * LibTomMath is a library that provides multiple-precision
992 * integer arithmetic as well as number theoretic functionality.
994 * The library was designed directly after the MPI library by
995 * Michael Fromberger but has been written from scratch with
996 * additional optimizations in place.
998 * The library is free for all purposes without any express
999 * guarantee it works.
1001 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
1004 /* AND two ints together */
1006 mp_and (mp_int * a, mp_int * b, mp_int * c)
1011 if (a->used > b->used) {
1012 if ((res = mp_init_copy (&t, a)) != MP_OKAY) {
1018 if ((res = mp_init_copy (&t, b)) != MP_OKAY) {
1025 for (ix = 0; ix < px; ix++) {
1026 t.dp[ix] &= x->dp[ix];
1029 /* zero digits above the last from the smallest mp_int */
1030 for (; ix < t.used; ix++) {
1041 /* End: bn_mp_and.c */
1043 /* Start: bn_mp_clamp.c */
1045 #ifdef BN_MP_CLAMP_C
1046 /* LibTomMath, multiple-precision integer library -- Tom St Denis
1048 * LibTomMath is a library that provides multiple-precision
1049 * integer arithmetic as well as number theoretic functionality.
1051 * The library was designed directly after the MPI library by
1052 * Michael Fromberger but has been written from scratch with
1053 * additional optimizations in place.
1055 * The library is free for all purposes without any express
1056 * guarantee it works.
1058 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
1061 /* trim unused digits
1063 * This is used to ensure that leading zero digits are
1064 * trimed and the leading "used" digit will be non-zero
1065 * Typically very fast. Also fixes the sign if there
1066 * are no more leading digits
1069 mp_clamp (mp_int * a)
1071 /* decrease used while the most significant digit is
1074 while (a->used > 0 && a->dp[a->used - 1] == 0) {
1078 /* reset the sign flag if used == 0 */
1085 /* End: bn_mp_clamp.c */
1087 /* Start: bn_mp_clear.c */
1089 #ifdef BN_MP_CLEAR_C
1090 /* LibTomMath, multiple-precision integer library -- Tom St Denis
1092 * LibTomMath is a library that provides multiple-precision
1093 * integer arithmetic as well as number theoretic functionality.
1095 * The library was designed directly after the MPI library by
1096 * Michael Fromberger but has been written from scratch with
1097 * additional optimizations in place.
1099 * The library is free for all purposes without any express
1100 * guarantee it works.
1102 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
1105 /* clear one (frees) */
1107 mp_clear (mp_int * a)
1111 /* only do anything if a hasn't been freed previously */
1112 if (a->dp != NULL) {
1113 /* first zero the digits */
1114 for (i = 0; i < a->used; i++) {
1121 /* reset members to make debugging easier */
1123 a->alloc = a->used = 0;
1129 /* End: bn_mp_clear.c */
1131 /* Start: bn_mp_clear_multi.c */
1133 #ifdef BN_MP_CLEAR_MULTI_C
1134 /* LibTomMath, multiple-precision integer library -- Tom St Denis
1136 * LibTomMath is a library that provides multiple-precision
1137 * integer arithmetic as well as number theoretic functionality.
1139 * The library was designed directly after the MPI library by
1140 * Michael Fromberger but has been written from scratch with
1141 * additional optimizations in place.
1143 * The library is free for all purposes without any express
1144 * guarantee it works.
1146 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
1150 void mp_clear_multi(mp_int *mp, ...)
1152 mp_int* next_mp = mp;
1155 while (next_mp != NULL) {
1157 next_mp = va_arg(args, mp_int*);
1163 /* End: bn_mp_clear_multi.c */
1165 /* Start: bn_mp_cmp.c */
1168 /* LibTomMath, multiple-precision integer library -- Tom St Denis
1170 * LibTomMath is a library that provides multiple-precision
1171 * integer arithmetic as well as number theoretic functionality.
1173 * The library was designed directly after the MPI library by
1174 * Michael Fromberger but has been written from scratch with
1175 * additional optimizations in place.
1177 * The library is free for all purposes without any express
1178 * guarantee it works.
1180 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
1183 /* compare two ints (signed)*/
1185 mp_cmp (mp_int * a, mp_int * b)
1187 /* compare based on sign */
1188 if (a->sign != b->sign) {
1189 if (a->sign == MP_NEG) {
1196 /* compare digits */
1197 if (a->sign == MP_NEG) {
1198 /* if negative compare opposite direction */
1199 return mp_cmp_mag(b, a);
1201 return mp_cmp_mag(a, b);
1206 /* End: bn_mp_cmp.c */
1208 /* Start: bn_mp_cmp_d.c */
1210 #ifdef BN_MP_CMP_D_C
1211 /* LibTomMath, multiple-precision integer library -- Tom St Denis
1213 * LibTomMath is a library that provides multiple-precision
1214 * integer arithmetic as well as number theoretic functionality.
1216 * The library was designed directly after the MPI library by
1217 * Michael Fromberger but has been written from scratch with
1218 * additional optimizations in place.
1220 * The library is free for all purposes without any express
1221 * guarantee it works.
1223 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
1226 /* compare a digit */
1227 int mp_cmp_d(mp_int * a, mp_digit b)
1229 /* compare based on sign */
1230 if (a->sign == MP_NEG) {
1234 /* compare based on magnitude */
1239 /* compare the only digit of a to b */
1242 } else if (a->dp[0] < b) {
1250 /* End: bn_mp_cmp_d.c */
1252 /* Start: bn_mp_cmp_mag.c */
1254 #ifdef BN_MP_CMP_MAG_C
1255 /* LibTomMath, multiple-precision integer library -- Tom St Denis
1257 * LibTomMath is a library that provides multiple-precision
1258 * integer arithmetic as well as number theoretic functionality.
1260 * The library was designed directly after the MPI library by
1261 * Michael Fromberger but has been written from scratch with
1262 * additional optimizations in place.
1264 * The library is free for all purposes without any express
1265 * guarantee it works.
1267 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
1270 /* compare maginitude of two ints (unsigned) */
1271 int mp_cmp_mag (mp_int * a, mp_int * b)
1274 mp_digit *tmpa, *tmpb;
1276 /* compare based on # of non-zero digits */
1277 if (a->used > b->used) {
1281 if (a->used < b->used) {
1286 tmpa = a->dp + (a->used - 1);
1289 tmpb = b->dp + (a->used - 1);
1291 /* compare based on digits */
1292 for (n = 0; n < a->used; ++n, --tmpa, --tmpb) {
1293 if (*tmpa > *tmpb) {
1297 if (*tmpa < *tmpb) {
1305 /* End: bn_mp_cmp_mag.c */
1307 /* Start: bn_mp_cnt_lsb.c */
1309 #ifdef BN_MP_CNT_LSB_C
1310 /* LibTomMath, multiple-precision integer library -- Tom St Denis
1312 * LibTomMath is a library that provides multiple-precision
1313 * integer arithmetic as well as number theoretic functionality.
1315 * The library was designed directly after the MPI library by
1316 * Michael Fromberger but has been written from scratch with
1317 * additional optimizations in place.
1319 * The library is free for all purposes without any express
1320 * guarantee it works.
1322 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
1325 static const int lnz[16] = {
1326 4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0
1329 /* Counts the number of lsbs which are zero before the first zero bit */
1330 int mp_cnt_lsb(mp_int *a)
1336 if (mp_iszero(a) == 1) {
1340 /* scan lower digits until non-zero */
1341 for (x = 0; x < a->used && a->dp[x] == 0; x++);
1345 /* now scan this digit until a 1 is found */
1358 /* End: bn_mp_cnt_lsb.c */
1360 /* Start: bn_mp_copy.c */
1363 /* LibTomMath, multiple-precision integer library -- Tom St Denis
1365 * LibTomMath is a library that provides multiple-precision
1366 * integer arithmetic as well as number theoretic functionality.
1368 * The library was designed directly after the MPI library by
1369 * Michael Fromberger but has been written from scratch with
1370 * additional optimizations in place.
1372 * The library is free for all purposes without any express
1373 * guarantee it works.
1375 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
1380 mp_copy (mp_int * a, mp_int * b)
1384 /* if dst == src do nothing */
1390 if (b->alloc < a->used) {
1391 if ((res = mp_grow (b, a->used)) != MP_OKAY) {
1396 /* zero b and copy the parameters over */
1398 register mp_digit *tmpa, *tmpb;
1400 /* pointer aliases */
1408 /* copy all the digits */
1409 for (n = 0; n < a->used; n++) {
1413 /* clear high digits */
1414 for (; n < b->used; n++) {
1419 /* copy used count and sign */
1426 /* End: bn_mp_copy.c */
1428 /* Start: bn_mp_count_bits.c */
1430 #ifdef BN_MP_COUNT_BITS_C
1431 /* LibTomMath, multiple-precision integer library -- Tom St Denis
1433 * LibTomMath is a library that provides multiple-precision
1434 * integer arithmetic as well as number theoretic functionality.
1436 * The library was designed directly after the MPI library by
1437 * Michael Fromberger but has been written from scratch with
1438 * additional optimizations in place.
1440 * The library is free for all purposes without any express
1441 * guarantee it works.
1443 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
1446 /* returns the number of bits in an int */
1448 mp_count_bits (mp_int * a)
1458 /* get number of digits and add that */
1459 r = (a->used - 1) * DIGIT_BIT;
1461 /* take the last digit and count the bits in it */
1462 q = a->dp[a->used - 1];
1463 while (q > ((mp_digit) 0)) {
1465 q >>= ((mp_digit) 1);
1471 /* End: bn_mp_count_bits.c */
1473 /* Start: bn_mp_div.c */
1476 /* LibTomMath, multiple-precision integer library -- Tom St Denis
1478 * LibTomMath is a library that provides multiple-precision
1479 * integer arithmetic as well as number theoretic functionality.
1481 * The library was designed directly after the MPI library by
1482 * Michael Fromberger but has been written from scratch with
1483 * additional optimizations in place.
1485 * The library is free for all purposes without any express
1486 * guarantee it works.
1488 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
1491 #ifdef BN_MP_DIV_SMALL
1493 /* slower bit-bang division... also smaller */
1494 int mp_div(mp_int * a, mp_int * b, mp_int * c, mp_int * d)
1496 mp_int ta, tb, tq, q;
1499 /* is divisor zero ? */
1500 if (mp_iszero (b) == 1) {
1504 /* if a < b then q=0, r = a */
1505 if (mp_cmp_mag (a, b) == MP_LT) {
1507 res = mp_copy (a, d);
1517 /* init our temps */
1518 if ((res = mp_init_multi(&ta, &tb, &tq, &q, NULL) != MP_OKAY)) {
1524 n = mp_count_bits(a) - mp_count_bits(b);
1525 if (((res = mp_abs(a, &ta)) != MP_OKAY) ||
1526 ((res = mp_abs(b, &tb)) != MP_OKAY) ||
1527 ((res = mp_mul_2d(&tb, n, &tb)) != MP_OKAY) ||
1528 ((res = mp_mul_2d(&tq, n, &tq)) != MP_OKAY)) {
1533 if (mp_cmp(&tb, &ta) != MP_GT) {
1534 if (((res = mp_sub(&ta, &tb, &ta)) != MP_OKAY) ||
1535 ((res = mp_add(&q, &tq, &q)) != MP_OKAY)) {
1539 if (((res = mp_div_2d(&tb, 1, &tb, NULL)) != MP_OKAY) ||
1540 ((res = mp_div_2d(&tq, 1, &tq, NULL)) != MP_OKAY)) {
1545 /* now q == quotient and ta == remainder */
1547 n2 = (a->sign == b->sign ? MP_ZPOS : MP_NEG);
1550 c->sign = (mp_iszero(c) == MP_YES) ? MP_ZPOS : n2;
1554 d->sign = (mp_iszero(d) == MP_YES) ? MP_ZPOS : n;
1557 mp_clear_multi(&ta, &tb, &tq, &q, NULL);
1563 /* integer signed division.
1564 * c*b + d == a [e.g. a/b, c=quotient, d=remainder]
1565 * HAC pp.598 Algorithm 14.20
1567 * Note that the description in HAC is horribly
1568 * incomplete. For example, it doesn't consider
1569 * the case where digits are removed from 'x' in
1570 * the inner loop. It also doesn't consider the
1571 * case that y has fewer than three digits, etc..
1573 * The overall algorithm is as described as
1574 * 14.20 from HAC but fixed to treat these cases.
1576 int mp_div (mp_int * a, mp_int * b, mp_int * c, mp_int * d)
1578 mp_int q, x, y, t1, t2;
1579 int res, n, t, i, norm, neg;
1581 /* is divisor zero ? */
1582 if (mp_iszero (b) == 1) {
1586 /* if a < b then q=0, r = a */
1587 if (mp_cmp_mag (a, b) == MP_LT) {
1589 res = mp_copy (a, d);
1599 if ((res = mp_init_size (&q, a->used + 2)) != MP_OKAY) {
1602 q.used = a->used + 2;
1604 if ((res = mp_init (&t1)) != MP_OKAY) {
1608 if ((res = mp_init (&t2)) != MP_OKAY) {
1612 if ((res = mp_init_copy (&x, a)) != MP_OKAY) {
1616 if ((res = mp_init_copy (&y, b)) != MP_OKAY) {
1621 neg = (a->sign == b->sign) ? MP_ZPOS : MP_NEG;
1622 x.sign = y.sign = MP_ZPOS;
1624 /* normalize both x and y, ensure that y >= b/2, [b == 2**DIGIT_BIT] */
1625 norm = mp_count_bits(&y) % DIGIT_BIT;
1626 if (norm < (int)(DIGIT_BIT-1)) {
1627 norm = (DIGIT_BIT-1) - norm;
1628 if ((res = mp_mul_2d (&x, norm, &x)) != MP_OKAY) {
1631 if ((res = mp_mul_2d (&y, norm, &y)) != MP_OKAY) {
1638 /* note hac does 0 based, so if used==5 then its 0,1,2,3,4, e.g. use 4 */
1642 /* while (x >= y*b**n-t) do { q[n-t] += 1; x -= y*b**{n-t} } */
1643 if ((res = mp_lshd (&y, n - t)) != MP_OKAY) { /* y = y*b**{n-t} */
1647 while (mp_cmp (&x, &y) != MP_LT) {
1649 if ((res = mp_sub (&x, &y, &x)) != MP_OKAY) {
1654 /* reset y by shifting it back down */
1655 mp_rshd (&y, n - t);
1657 /* step 3. for i from n down to (t + 1) */
1658 for (i = n; i >= (t + 1); i--) {
1663 /* step 3.1 if xi == yt then set q{i-t-1} to b-1,
1664 * otherwise set q{i-t-1} to (xi*b + x{i-1})/yt */
1665 if (x.dp[i] == y.dp[t]) {
1666 q.dp[i - t - 1] = ((((mp_digit)1) << DIGIT_BIT) - 1);
1669 tmp = ((mp_word) x.dp[i]) << ((mp_word) DIGIT_BIT);
1670 tmp |= ((mp_word) x.dp[i - 1]);
1671 tmp /= ((mp_word) y.dp[t]);
1672 if (tmp > (mp_word) MP_MASK)
1674 q.dp[i - t - 1] = (mp_digit) (tmp & (mp_word) (MP_MASK));
1677 /* while (q{i-t-1} * (yt * b + y{t-1})) >
1678 xi * b**2 + xi-1 * b + xi-2
1682 q.dp[i - t - 1] = (q.dp[i - t - 1] + 1) & MP_MASK;
1684 q.dp[i - t - 1] = (q.dp[i - t - 1] - 1) & MP_MASK;
1686 /* find left hand */
1688 t1.dp[0] = (t - 1 < 0) ? 0 : y.dp[t - 1];
1691 if ((res = mp_mul_d (&t1, q.dp[i - t - 1], &t1)) != MP_OKAY) {
1695 /* find right hand */
1696 t2.dp[0] = (i - 2 < 0) ? 0 : x.dp[i - 2];
1697 t2.dp[1] = (i - 1 < 0) ? 0 : x.dp[i - 1];
1700 } while (mp_cmp_mag(&t1, &t2) == MP_GT);
1702 /* step 3.3 x = x - q{i-t-1} * y * b**{i-t-1} */
1703 if ((res = mp_mul_d (&y, q.dp[i - t - 1], &t1)) != MP_OKAY) {
1707 if ((res = mp_lshd (&t1, i - t - 1)) != MP_OKAY) {
1711 if ((res = mp_sub (&x, &t1, &x)) != MP_OKAY) {
1715 /* if x < 0 then { x = x + y*b**{i-t-1}; q{i-t-1} -= 1; } */
1716 if (x.sign == MP_NEG) {
1717 if ((res = mp_copy (&y, &t1)) != MP_OKAY) {
1720 if ((res = mp_lshd (&t1, i - t - 1)) != MP_OKAY) {
1723 if ((res = mp_add (&x, &t1, &x)) != MP_OKAY) {
1727 q.dp[i - t - 1] = (q.dp[i - t - 1] - 1UL) & MP_MASK;
1731 /* now q is the quotient and x is the remainder
1732 * [which we have to normalize]
1735 /* get sign before writing to c */
1736 x.sign = x.used == 0 ? MP_ZPOS : a->sign;
1745 mp_div_2d (&x, norm, &x, NULL);
1751 LBL_Y:mp_clear (&y);
1752 LBL_X:mp_clear (&x);
1753 LBL_T2:mp_clear (&t2);
1754 LBL_T1:mp_clear (&t1);
1755 LBL_Q:mp_clear (&q);
1763 /* End: bn_mp_div.c */
1765 /* Start: bn_mp_div_2.c */
1767 #ifdef BN_MP_DIV_2_C
1768 /* LibTomMath, multiple-precision integer library -- Tom St Denis
1770 * LibTomMath is a library that provides multiple-precision
1771 * integer arithmetic as well as number theoretic functionality.
1773 * The library was designed directly after the MPI library by
1774 * Michael Fromberger but has been written from scratch with
1775 * additional optimizations in place.
1777 * The library is free for all purposes without any express
1778 * guarantee it works.
1780 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
1784 int mp_div_2(mp_int * a, mp_int * b)
1786 int x, res, oldused;
1789 if (b->alloc < a->used) {
1790 if ((res = mp_grow (b, a->used)) != MP_OKAY) {
1798 register mp_digit r, rr, *tmpa, *tmpb;
1801 tmpa = a->dp + b->used - 1;
1804 tmpb = b->dp + b->used - 1;
1808 for (x = b->used - 1; x >= 0; x--) {
1809 /* get the carry for the next iteration */
1812 /* shift the current digit, add in carry and store */
1813 *tmpb-- = (*tmpa-- >> 1) | (r << (DIGIT_BIT - 1));
1815 /* forward carry to next iteration */
1819 /* zero excess digits */
1820 tmpb = b->dp + b->used;
1821 for (x = b->used; x < oldused; x++) {
1831 /* End: bn_mp_div_2.c */
1833 /* Start: bn_mp_div_2d.c */
1835 #ifdef BN_MP_DIV_2D_C
1836 /* LibTomMath, multiple-precision integer library -- Tom St Denis
1838 * LibTomMath is a library that provides multiple-precision
1839 * integer arithmetic as well as number theoretic functionality.
1841 * The library was designed directly after the MPI library by
1842 * Michael Fromberger but has been written from scratch with
1843 * additional optimizations in place.
1845 * The library is free for all purposes without any express
1846 * guarantee it works.
1848 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
1851 /* shift right by a certain bit count (store quotient in c, optional remainder in d) */
1852 int mp_div_2d (mp_int * a, int b, mp_int * c, mp_int * d)
1859 /* if the shift count is <= 0 then we do no work */
1861 res = mp_copy (a, c);
1868 if ((res = mp_init (&t)) != MP_OKAY) {
1872 /* get the remainder */
1874 if ((res = mp_mod_2d (a, b, &t)) != MP_OKAY) {
1881 if ((res = mp_copy (a, c)) != MP_OKAY) {
1886 /* shift by as many digits in the bit count */
1887 if (b >= (int)DIGIT_BIT) {
1888 mp_rshd (c, b / DIGIT_BIT);
1891 /* shift any bit count < DIGIT_BIT */
1892 D = (mp_digit) (b % DIGIT_BIT);
1894 register mp_digit *tmpc, mask, shift;
1897 mask = (((mp_digit)1) << D) - 1;
1900 shift = DIGIT_BIT - D;
1903 tmpc = c->dp + (c->used - 1);
1907 for (x = c->used - 1; x >= 0; x--) {
1908 /* get the lower bits of this word in a temp */
1911 /* shift the current word and mix in the carry bits from the previous word */
1912 *tmpc = (*tmpc >> D) | (r << shift);
1915 /* set the carry to the carry bits of the current word found above */
1928 /* End: bn_mp_div_2d.c */
1930 /* Start: bn_mp_div_3.c */
1932 #ifdef BN_MP_DIV_3_C
1933 /* LibTomMath, multiple-precision integer library -- Tom St Denis
1935 * LibTomMath is a library that provides multiple-precision
1936 * integer arithmetic as well as number theoretic functionality.
1938 * The library was designed directly after the MPI library by
1939 * Michael Fromberger but has been written from scratch with
1940 * additional optimizations in place.
1942 * The library is free for all purposes without any express
1943 * guarantee it works.
1945 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
1948 /* divide by three (based on routine from MPI and the GMP manual) */
1950 mp_div_3 (mp_int * a, mp_int *c, mp_digit * d)
1957 /* b = 2**DIGIT_BIT / 3 */
1958 b = (((mp_word)1) << ((mp_word)DIGIT_BIT)) / ((mp_word)3);
1960 if ((res = mp_init_size(&q, a->used)) != MP_OKAY) {
1967 for (ix = a->used - 1; ix >= 0; ix--) {
1968 w = (w << ((mp_word)DIGIT_BIT)) | ((mp_word)a->dp[ix]);
1971 /* multiply w by [1/3] */
1972 t = (w * ((mp_word)b)) >> ((mp_word)DIGIT_BIT);
1974 /* now subtract 3 * [w/3] from w, to get the remainder */
1977 /* fixup the remainder as required since
1978 * the optimization is not exact.
1987 q.dp[ix] = (mp_digit)t;
1990 /* [optional] store the remainder */
1995 /* [optional] store the quotient */
2007 /* End: bn_mp_div_3.c */
2009 /* Start: bn_mp_div_d.c */
2011 #ifdef BN_MP_DIV_D_C
2012 /* LibTomMath, multiple-precision integer library -- Tom St Denis
2014 * LibTomMath is a library that provides multiple-precision
2015 * integer arithmetic as well as number theoretic functionality.
2017 * The library was designed directly after the MPI library by
2018 * Michael Fromberger but has been written from scratch with
2019 * additional optimizations in place.
2021 * The library is free for all purposes without any express
2022 * guarantee it works.
2024 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
2027 static int s_is_power_of_two(mp_digit b, int *p)
2031 for (x = 1; x < DIGIT_BIT; x++) {
2032 if (b == (((mp_digit)1)<<x)) {
2040 /* single digit division (based on routine from MPI) */
2041 int mp_div_d (mp_int * a, mp_digit b, mp_int * c, mp_digit * d)
2048 /* cannot divide by zero */
2054 if (b == 1 || mp_iszero(a) == 1) {
2059 return mp_copy(a, c);
2064 /* power of two ? */
2065 if (s_is_power_of_two(b, &ix) == 1) {
2067 *d = a->dp[0] & ((((mp_digit)1)<<ix) - 1);
2070 return mp_div_2d(a, ix, c, NULL);
2075 #ifdef BN_MP_DIV_3_C
2078 return mp_div_3(a, c, d);
2082 /* no easy answer [c'est la vie]. Just division */
2083 if ((res = mp_init_size(&q, a->used)) != MP_OKAY) {
2090 for (ix = a->used - 1; ix >= 0; ix--) {
2091 w = (w << ((mp_word)DIGIT_BIT)) | ((mp_word)a->dp[ix]);
2094 t = (mp_digit)(w / b);
2095 w -= ((mp_word)t) * ((mp_word)b);
2099 q.dp[ix] = (mp_digit)t;
2117 /* End: bn_mp_div_d.c */
2119 /* Start: bn_mp_dr_is_modulus.c */
2121 #ifdef BN_MP_DR_IS_MODULUS_C
2122 /* LibTomMath, multiple-precision integer library -- Tom St Denis
2124 * LibTomMath is a library that provides multiple-precision
2125 * integer arithmetic as well as number theoretic functionality.
2127 * The library was designed directly after the MPI library by
2128 * Michael Fromberger but has been written from scratch with
2129 * additional optimizations in place.
2131 * The library is free for all purposes without any express
2132 * guarantee it works.
2134 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
2137 /* determines if a number is a valid DR modulus */
2138 int mp_dr_is_modulus(mp_int *a)
2142 /* must be at least two digits */
2147 /* must be of the form b**k - a [a <= b] so all
2148 * but the first digit must be equal to -1 (mod b).
2150 for (ix = 1; ix < a->used; ix++) {
2151 if (a->dp[ix] != MP_MASK) {
2160 /* End: bn_mp_dr_is_modulus.c */
2162 /* Start: bn_mp_dr_reduce.c */
2164 #ifdef BN_MP_DR_REDUCE_C
2165 /* LibTomMath, multiple-precision integer library -- Tom St Denis
2167 * LibTomMath is a library that provides multiple-precision
2168 * integer arithmetic as well as number theoretic functionality.
2170 * The library was designed directly after the MPI library by
2171 * Michael Fromberger but has been written from scratch with
2172 * additional optimizations in place.
2174 * The library is free for all purposes without any express
2175 * guarantee it works.
2177 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
2180 /* reduce "x" in place modulo "n" using the Diminished Radix algorithm.
2182 * Based on algorithm from the paper
2184 * "Generating Efficient Primes for Discrete Log Cryptosystems"
2185 * Chae Hoon Lim, Pil Joong Lee,
2186 * POSTECH Information Research Laboratories
2188 * The modulus must be of a special format [see manual]
2190 * Has been modified to use algorithm 7.10 from the LTM book instead
2192 * Input x must be in the range 0 <= x <= (n-1)**2
2195 mp_dr_reduce (mp_int * x, mp_int * n, mp_digit k)
2199 mp_digit mu, *tmpx1, *tmpx2;
2201 /* m = digits in modulus */
2204 /* ensure that "x" has at least 2m digits */
2205 if (x->alloc < m + m) {
2206 if ((err = mp_grow (x, m + m)) != MP_OKAY) {
2211 /* top of loop, this is where the code resumes if
2212 * another reduction pass is required.
2215 /* aliases for digits */
2216 /* alias for lower half of x */
2219 /* alias for upper half of x, or x/B**m */
2222 /* set carry to zero */
2225 /* compute (x mod B**m) + k * [x/B**m] inline and inplace */
2226 for (i = 0; i < m; i++) {
2227 r = ((mp_word)*tmpx2++) * ((mp_word)k) + *tmpx1 + mu;
2228 *tmpx1++ = (mp_digit)(r & MP_MASK);
2229 mu = (mp_digit)(r >> ((mp_word)DIGIT_BIT));
2232 /* set final carry */
2235 /* zero words above m */
2236 for (i = m + 1; i < x->used; i++) {
2240 /* clamp, sub and return */
2243 /* if x >= n then subtract and reduce again
2244 * Each successive "recursion" makes the input smaller and smaller.
2246 if (mp_cmp_mag (x, n) != MP_LT) {
2254 /* End: bn_mp_dr_reduce.c */
2256 /* Start: bn_mp_dr_setup.c */
2258 #ifdef BN_MP_DR_SETUP_C
2259 /* LibTomMath, multiple-precision integer library -- Tom St Denis
2261 * LibTomMath is a library that provides multiple-precision
2262 * integer arithmetic as well as number theoretic functionality.
2264 * The library was designed directly after the MPI library by
2265 * Michael Fromberger but has been written from scratch with
2266 * additional optimizations in place.
2268 * The library is free for all purposes without any express
2269 * guarantee it works.
2271 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
2274 /* determines the setup value */
2275 void mp_dr_setup(mp_int *a, mp_digit *d)
2277 /* the casts are required if DIGIT_BIT is one less than
2278 * the number of bits in a mp_digit [e.g. DIGIT_BIT==31]
2280 *d = (mp_digit)((((mp_word)1) << ((mp_word)DIGIT_BIT)) -
2281 ((mp_word)a->dp[0]));
2286 /* End: bn_mp_dr_setup.c */
2288 /* Start: bn_mp_exch.c */
2291 /* LibTomMath, multiple-precision integer library -- Tom St Denis
2293 * LibTomMath is a library that provides multiple-precision
2294 * integer arithmetic as well as number theoretic functionality.
2296 * The library was designed directly after the MPI library by
2297 * Michael Fromberger but has been written from scratch with
2298 * additional optimizations in place.
2300 * The library is free for all purposes without any express
2301 * guarantee it works.
2303 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
2306 /* swap the elements of two integers, for cases where you can't simply swap the
2307 * mp_int pointers around
2310 mp_exch (mp_int * a, mp_int * b)
2320 /* End: bn_mp_exch.c */
2322 /* Start: bn_mp_expt_d.c */
2324 #ifdef BN_MP_EXPT_D_C
2325 /* LibTomMath, multiple-precision integer library -- Tom St Denis
2327 * LibTomMath is a library that provides multiple-precision
2328 * integer arithmetic as well as number theoretic functionality.
2330 * The library was designed directly after the MPI library by
2331 * Michael Fromberger but has been written from scratch with
2332 * additional optimizations in place.
2334 * The library is free for all purposes without any express
2335 * guarantee it works.
2337 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
2340 /* calculate c = a**b using a square-multiply algorithm */
2341 int mp_expt_d (mp_int * a, mp_digit b, mp_int * c)
2346 if ((res = mp_init_copy (&g, a)) != MP_OKAY) {
2350 /* set initial result */
2353 for (x = 0; x < (int) DIGIT_BIT; x++) {
2355 if ((res = mp_sqr (c, c)) != MP_OKAY) {
2360 /* if the bit is set multiply */
2361 if ((b & (mp_digit) (((mp_digit)1) << (DIGIT_BIT - 1))) != 0) {
2362 if ((res = mp_mul (c, &g, c)) != MP_OKAY) {
2368 /* shift to next bit */
2377 /* End: bn_mp_expt_d.c */
2379 /* Start: bn_mp_exptmod.c */
2381 #ifdef BN_MP_EXPTMOD_C
2382 /* LibTomMath, multiple-precision integer library -- Tom St Denis
2384 * LibTomMath is a library that provides multiple-precision
2385 * integer arithmetic as well as number theoretic functionality.
2387 * The library was designed directly after the MPI library by
2388 * Michael Fromberger but has been written from scratch with
2389 * additional optimizations in place.
2391 * The library is free for all purposes without any express
2392 * guarantee it works.
2394 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
2398 /* this is a shell function that calls either the normal or Montgomery
2399 * exptmod functions. Originally the call to the montgomery code was
2400 * embedded in the normal function but that wasted alot of stack space
2401 * for nothing (since 99% of the time the Montgomery code would be called)
2403 int mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y)
2407 /* modulus P must be positive */
2408 if (P->sign == MP_NEG) {
2412 /* if exponent X is negative we have to recurse */
2413 if (X->sign == MP_NEG) {
2414 #ifdef BN_MP_INVMOD_C
2418 /* first compute 1/G mod P */
2419 if ((err = mp_init(&tmpG)) != MP_OKAY) {
2422 if ((err = mp_invmod(G, P, &tmpG)) != MP_OKAY) {
2428 if ((err = mp_init(&tmpX)) != MP_OKAY) {
2432 if ((err = mp_abs(X, &tmpX)) != MP_OKAY) {
2433 mp_clear_multi(&tmpG, &tmpX, NULL);
2437 /* and now compute (1/G)**|X| instead of G**X [X < 0] */
2438 err = mp_exptmod(&tmpG, &tmpX, P, Y);
2439 mp_clear_multi(&tmpG, &tmpX, NULL);
2447 /* modified diminished radix reduction */
2448 #if defined(BN_MP_REDUCE_IS_2K_L_C) && defined(BN_MP_REDUCE_2K_L_C)
2449 if (mp_reduce_is_2k_l(P) == MP_YES) {
2450 return s_mp_exptmod(G, X, P, Y, 1);
2454 #ifdef BN_MP_DR_IS_MODULUS_C
2455 /* is it a DR modulus? */
2456 dr = mp_dr_is_modulus(P);
2462 #ifdef BN_MP_REDUCE_IS_2K_C
2463 /* if not, is it a unrestricted DR modulus? */
2465 dr = mp_reduce_is_2k(P) << 1;
2469 /* if the modulus is odd or dr != 0 use the montgomery method */
2470 #ifdef BN_MP_EXPTMOD_FAST_C
2471 if (mp_isodd (P) == 1 || dr != 0) {
2472 return mp_exptmod_fast (G, X, P, Y, dr);
2475 #ifdef BN_S_MP_EXPTMOD_C
2476 /* otherwise use the generic Barrett reduction technique */
2477 return s_mp_exptmod (G, X, P, Y, 0);
2479 /* no exptmod for evens */
2482 #ifdef BN_MP_EXPTMOD_FAST_C
2489 /* End: bn_mp_exptmod.c */
2491 /* Start: bn_mp_exptmod_fast.c */
2493 #ifdef BN_MP_EXPTMOD_FAST_C
2494 /* LibTomMath, multiple-precision integer library -- Tom St Denis
2496 * LibTomMath is a library that provides multiple-precision
2497 * integer arithmetic as well as number theoretic functionality.
2499 * The library was designed directly after the MPI library by
2500 * Michael Fromberger but has been written from scratch with
2501 * additional optimizations in place.
2503 * The library is free for all purposes without any express
2504 * guarantee it works.
2506 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
2509 /* computes Y == G**X mod P, HAC pp.616, Algorithm 14.85
2511 * Uses a left-to-right k-ary sliding window to compute the modular exponentiation.
2512 * The value of k changes based on the size of the exponent.
2514 * Uses Montgomery or Diminished Radix reduction [whichever appropriate]
2520 #define TAB_SIZE 256
2523 int mp_exptmod_fast (mp_int * G, mp_int * X, mp_int * P, mp_int * Y, int redmode)
2525 mp_int M[TAB_SIZE], res;
2527 int err, bitbuf, bitcpy, bitcnt, mode, digidx, x, y, winsize;
2529 /* use a pointer to the reduction algorithm. This allows us to use
2530 * one of many reduction algorithms without modding the guts of
2531 * the code with if statements everywhere.
2533 int (*redux)(mp_int*,mp_int*,mp_digit);
2535 /* find window size */
2536 x = mp_count_bits (X);
2539 } else if (x <= 36) {
2541 } else if (x <= 140) {
2543 } else if (x <= 450) {
2545 } else if (x <= 1303) {
2547 } else if (x <= 3529) {
2560 /* init first cell */
2561 if ((err = mp_init(&M[1])) != MP_OKAY) {
2565 /* now init the second half of the array */
2566 for (x = 1<<(winsize-1); x < (1 << winsize); x++) {
2567 if ((err = mp_init(&M[x])) != MP_OKAY) {
2568 for (y = 1<<(winsize-1); y < x; y++) {
2576 /* determine and setup reduction code */
2578 #ifdef BN_MP_MONTGOMERY_SETUP_C
2579 /* now setup montgomery */
2580 if ((err = mp_montgomery_setup (P, &mp)) != MP_OKAY) {
2588 /* automatically pick the comba one if available (saves quite a few calls/ifs) */
2589 #ifdef BN_FAST_MP_MONTGOMERY_REDUCE_C
2590 if (((P->used * 2 + 1) < MP_WARRAY) &&
2591 P->used < (1 << ((CHAR_BIT * sizeof (mp_word)) - (2 * DIGIT_BIT)))) {
2592 redux = fast_mp_montgomery_reduce;
2596 #ifdef BN_MP_MONTGOMERY_REDUCE_C
2597 /* use slower baseline Montgomery method */
2598 redux = mp_montgomery_reduce;
2604 } else if (redmode == 1) {
2605 #if defined(BN_MP_DR_SETUP_C) && defined(BN_MP_DR_REDUCE_C)
2606 /* setup DR reduction for moduli of the form B**k - b */
2607 mp_dr_setup(P, &mp);
2608 redux = mp_dr_reduce;
2614 #if defined(BN_MP_REDUCE_2K_SETUP_C) && defined(BN_MP_REDUCE_2K_C)
2615 /* setup DR reduction for moduli of the form 2**k - b */
2616 if ((err = mp_reduce_2k_setup(P, &mp)) != MP_OKAY) {
2619 redux = mp_reduce_2k;
2627 if ((err = mp_init (&res)) != MP_OKAY) {
2635 * The first half of the table is not computed though accept for M[0] and M[1]
2639 #ifdef BN_MP_MONTGOMERY_CALC_NORMALIZATION_C
2640 /* now we need R mod m */
2641 if ((err = mp_montgomery_calc_normalization (&res, P)) != MP_OKAY) {
2649 /* now set M[1] to G * R mod m */
2650 if ((err = mp_mulmod (G, &res, P, &M[1])) != MP_OKAY) {
2655 if ((err = mp_mod(G, P, &M[1])) != MP_OKAY) {
2660 /* compute the value at M[1<<(winsize-1)] by squaring M[1] (winsize-1) times */
2661 if ((err = mp_copy (&M[1], &M[1 << (winsize - 1)])) != MP_OKAY) {
2665 for (x = 0; x < (winsize - 1); x++) {
2666 if ((err = mp_sqr (&M[1 << (winsize - 1)], &M[1 << (winsize - 1)])) != MP_OKAY) {
2669 if ((err = redux (&M[1 << (winsize - 1)], P, mp)) != MP_OKAY) {
2674 /* create upper table */
2675 for (x = (1 << (winsize - 1)) + 1; x < (1 << winsize); x++) {
2676 if ((err = mp_mul (&M[x - 1], &M[1], &M[x])) != MP_OKAY) {
2679 if ((err = redux (&M[x], P, mp)) != MP_OKAY) {
2684 /* set initial mode and bit cnt */
2688 digidx = X->used - 1;
2693 /* grab next digit as required */
2694 if (--bitcnt == 0) {
2695 /* if digidx == -1 we are out of digits so break */
2699 /* read next digit and reset bitcnt */
2700 buf = X->dp[digidx--];
2701 bitcnt = (int)DIGIT_BIT;
2704 /* grab the next msb from the exponent */
2705 y = (mp_digit)(buf >> (DIGIT_BIT - 1)) & 1;
2706 buf <<= (mp_digit)1;
2708 /* if the bit is zero and mode == 0 then we ignore it
2709 * These represent the leading zero bits before the first 1 bit
2710 * in the exponent. Technically this opt is not required but it
2711 * does lower the # of trivial squaring/reductions used
2713 if (mode == 0 && y == 0) {
2717 /* if the bit is zero and mode == 1 then we square */
2718 if (mode == 1 && y == 0) {
2719 if ((err = mp_sqr (&res, &res)) != MP_OKAY) {
2722 if ((err = redux (&res, P, mp)) != MP_OKAY) {
2728 /* else we add it to the window */
2729 bitbuf |= (y << (winsize - ++bitcpy));
2732 if (bitcpy == winsize) {
2733 /* ok window is filled so square as required and multiply */
2735 for (x = 0; x < winsize; x++) {
2736 if ((err = mp_sqr (&res, &res)) != MP_OKAY) {
2739 if ((err = redux (&res, P, mp)) != MP_OKAY) {
2745 if ((err = mp_mul (&res, &M[bitbuf], &res)) != MP_OKAY) {
2748 if ((err = redux (&res, P, mp)) != MP_OKAY) {
2752 /* empty window and reset */
2759 /* if bits remain then square/multiply */
2760 if (mode == 2 && bitcpy > 0) {
2761 /* square then multiply if the bit is set */
2762 for (x = 0; x < bitcpy; x++) {
2763 if ((err = mp_sqr (&res, &res)) != MP_OKAY) {
2766 if ((err = redux (&res, P, mp)) != MP_OKAY) {
2770 /* get next bit of the window */
2772 if ((bitbuf & (1 << winsize)) != 0) {
2774 if ((err = mp_mul (&res, &M[1], &res)) != MP_OKAY) {
2777 if ((err = redux (&res, P, mp)) != MP_OKAY) {
2785 /* fixup result if Montgomery reduction is used
2786 * recall that any value in a Montgomery system is
2787 * actually multiplied by R mod n. So we have
2788 * to reduce one more time to cancel out the factor
2791 if ((err = redux(&res, P, mp)) != MP_OKAY) {
2796 /* swap res with Y */
2799 LBL_RES:mp_clear (&res);
2802 for (x = 1<<(winsize-1); x < (1 << winsize); x++) {
2810 /* End: bn_mp_exptmod_fast.c */
2812 /* Start: bn_mp_exteuclid.c */
2814 #ifdef BN_MP_EXTEUCLID_C
2815 /* LibTomMath, multiple-precision integer library -- Tom St Denis
2817 * LibTomMath is a library that provides multiple-precision
2818 * integer arithmetic as well as number theoretic functionality.
2820 * The library was designed directly after the MPI library by
2821 * Michael Fromberger but has been written from scratch with
2822 * additional optimizations in place.
2824 * The library is free for all purposes without any express
2825 * guarantee it works.
2827 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
2830 /* Extended euclidean algorithm of (a, b) produces
2833 int mp_exteuclid(mp_int *a, mp_int *b, mp_int *U1, mp_int *U2, mp_int *U3)
2835 mp_int u1,u2,u3,v1,v2,v3,t1,t2,t3,q,tmp;
2838 if ((err = mp_init_multi(&u1, &u2, &u3, &v1, &v2, &v3, &t1, &t2, &t3, &q, &tmp, NULL)) != MP_OKAY) {
2842 /* initialize, (u1,u2,u3) = (1,0,a) */
2844 if ((err = mp_copy(a, &u3)) != MP_OKAY) { goto _ERR; }
2846 /* initialize, (v1,v2,v3) = (0,1,b) */
2848 if ((err = mp_copy(b, &v3)) != MP_OKAY) { goto _ERR; }
2850 /* loop while v3 != 0 */
2851 while (mp_iszero(&v3) == MP_NO) {
2853 if ((err = mp_div(&u3, &v3, &q, NULL)) != MP_OKAY) { goto _ERR; }
2855 /* (t1,t2,t3) = (u1,u2,u3) - (v1,v2,v3)q */
2856 if ((err = mp_mul(&v1, &q, &tmp)) != MP_OKAY) { goto _ERR; }
2857 if ((err = mp_sub(&u1, &tmp, &t1)) != MP_OKAY) { goto _ERR; }
2858 if ((err = mp_mul(&v2, &q, &tmp)) != MP_OKAY) { goto _ERR; }
2859 if ((err = mp_sub(&u2, &tmp, &t2)) != MP_OKAY) { goto _ERR; }
2860 if ((err = mp_mul(&v3, &q, &tmp)) != MP_OKAY) { goto _ERR; }
2861 if ((err = mp_sub(&u3, &tmp, &t3)) != MP_OKAY) { goto _ERR; }
2863 /* (u1,u2,u3) = (v1,v2,v3) */
2864 if ((err = mp_copy(&v1, &u1)) != MP_OKAY) { goto _ERR; }
2865 if ((err = mp_copy(&v2, &u2)) != MP_OKAY) { goto _ERR; }
2866 if ((err = mp_copy(&v3, &u3)) != MP_OKAY) { goto _ERR; }
2868 /* (v1,v2,v3) = (t1,t2,t3) */
2869 if ((err = mp_copy(&t1, &v1)) != MP_OKAY) { goto _ERR; }
2870 if ((err = mp_copy(&t2, &v2)) != MP_OKAY) { goto _ERR; }
2871 if ((err = mp_copy(&t3, &v3)) != MP_OKAY) { goto _ERR; }
2874 /* make sure U3 >= 0 */
2875 if (u3.sign == MP_NEG) {
2881 /* copy result out */
2882 if (U1 != NULL) { mp_exch(U1, &u1); }
2883 if (U2 != NULL) { mp_exch(U2, &u2); }
2884 if (U3 != NULL) { mp_exch(U3, &u3); }
2887 _ERR: mp_clear_multi(&u1, &u2, &u3, &v1, &v2, &v3, &t1, &t2, &t3, &q, &tmp, NULL);
2892 /* End: bn_mp_exteuclid.c */
2894 /* Start: bn_mp_fread.c */
2896 #ifdef BN_MP_FREAD_C
2897 /* LibTomMath, multiple-precision integer library -- Tom St Denis
2899 * LibTomMath is a library that provides multiple-precision
2900 * integer arithmetic as well as number theoretic functionality.
2902 * The library was designed directly after the MPI library by
2903 * Michael Fromberger but has been written from scratch with
2904 * additional optimizations in place.
2906 * The library is free for all purposes without any express
2907 * guarantee it works.
2909 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
2912 /* read a bigint from a file stream in ASCII */
2913 int mp_fread(mp_int *a, int radix, FILE *stream)
2915 int err, ch, neg, y;
2920 /* if first digit is - then set negative */
2930 /* find y in the radix map */
2931 for (y = 0; y < radix; y++) {
2932 if (mp_s_rmap[y] == ch) {
2940 /* shift up and add */
2941 if ((err = mp_mul_d(a, radix, a)) != MP_OKAY) {
2944 if ((err = mp_add_d(a, y, a)) != MP_OKAY) {
2950 if (mp_cmp_d(a, 0) != MP_EQ) {
2959 /* End: bn_mp_fread.c */
2961 /* Start: bn_mp_fwrite.c */
2963 #ifdef BN_MP_FWRITE_C
2964 /* LibTomMath, multiple-precision integer library -- Tom St Denis
2966 * LibTomMath is a library that provides multiple-precision
2967 * integer arithmetic as well as number theoretic functionality.
2969 * The library was designed directly after the MPI library by
2970 * Michael Fromberger but has been written from scratch with
2971 * additional optimizations in place.
2973 * The library is free for all purposes without any express
2974 * guarantee it works.
2976 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
2979 int mp_fwrite(mp_int *a, int radix, FILE *stream)
2984 if ((err = mp_radix_size(a, radix, &len)) != MP_OKAY) {
2988 buf = OPT_CAST(char) XMALLOC (len);
2993 if ((err = mp_toradix(a, buf, radix)) != MP_OKAY) {
2998 for (x = 0; x < len; x++) {
2999 if (fputc(buf[x], stream) == EOF) {
3011 /* End: bn_mp_fwrite.c */
3013 /* Start: bn_mp_gcd.c */
3016 /* LibTomMath, multiple-precision integer library -- Tom St Denis
3018 * LibTomMath is a library that provides multiple-precision
3019 * integer arithmetic as well as number theoretic functionality.
3021 * The library was designed directly after the MPI library by
3022 * Michael Fromberger but has been written from scratch with
3023 * additional optimizations in place.
3025 * The library is free for all purposes without any express
3026 * guarantee it works.
3028 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
3031 /* Greatest Common Divisor using the binary method */
3032 int mp_gcd (mp_int * a, mp_int * b, mp_int * c)
3035 int k, u_lsb, v_lsb, res;
3037 /* either zero than gcd is the largest */
3038 if (mp_iszero (a) == 1 && mp_iszero (b) == 0) {
3039 return mp_abs (b, c);
3041 if (mp_iszero (a) == 0 && mp_iszero (b) == 1) {
3042 return mp_abs (a, c);
3045 /* optimized. At this point if a == 0 then
3046 * b must equal zero too
3048 if (mp_iszero (a) == 1) {
3053 /* get copies of a and b we can modify */
3054 if ((res = mp_init_copy (&u, a)) != MP_OKAY) {
3058 if ((res = mp_init_copy (&v, b)) != MP_OKAY) {
3062 /* must be positive for the remainder of the algorithm */
3063 u.sign = v.sign = MP_ZPOS;
3065 /* B1. Find the common power of two for u and v */
3066 u_lsb = mp_cnt_lsb(&u);
3067 v_lsb = mp_cnt_lsb(&v);
3068 k = MIN(u_lsb, v_lsb);
3071 /* divide the power of two out */
3072 if ((res = mp_div_2d(&u, k, &u, NULL)) != MP_OKAY) {
3076 if ((res = mp_div_2d(&v, k, &v, NULL)) != MP_OKAY) {
3081 /* divide any remaining factors of two out */
3083 if ((res = mp_div_2d(&u, u_lsb - k, &u, NULL)) != MP_OKAY) {
3089 if ((res = mp_div_2d(&v, v_lsb - k, &v, NULL)) != MP_OKAY) {
3094 while (mp_iszero(&v) == 0) {
3095 /* make sure v is the largest */
3096 if (mp_cmp_mag(&u, &v) == MP_GT) {
3097 /* swap u and v to make sure v is >= u */
3101 /* subtract smallest from largest */
3102 if ((res = s_mp_sub(&v, &u, &v)) != MP_OKAY) {
3106 /* Divide out all factors of two */
3107 if ((res = mp_div_2d(&v, mp_cnt_lsb(&v), &v, NULL)) != MP_OKAY) {
3112 /* multiply by 2**k which we divided out at the beginning */
3113 if ((res = mp_mul_2d (&u, k, c)) != MP_OKAY) {
3118 LBL_V:mp_clear (&u);
3119 LBL_U:mp_clear (&v);
3124 /* End: bn_mp_gcd.c */
3126 /* Start: bn_mp_get_int.c */
3128 #ifdef BN_MP_GET_INT_C
3129 /* LibTomMath, multiple-precision integer library -- Tom St Denis
3131 * LibTomMath is a library that provides multiple-precision
3132 * integer arithmetic as well as number theoretic functionality.
3134 * The library was designed directly after the MPI library by
3135 * Michael Fromberger but has been written from scratch with
3136 * additional optimizations in place.
3138 * The library is free for all purposes without any express
3139 * guarantee it works.
3141 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
3144 /* get the lower 32-bits of an mp_int */
3145 unsigned long mp_get_int(mp_int * a)
3154 /* get number of digits of the lsb we have to read */
3155 i = MIN(a->used,(int)((sizeof(unsigned long)*CHAR_BIT+DIGIT_BIT-1)/DIGIT_BIT))-1;
3157 /* get most significant digit of result */
3161 res = (res << DIGIT_BIT) | DIGIT(a,i);
3164 /* force result to 32-bits always so it is consistent on non 32-bit platforms */
3165 return res & 0xFFFFFFFFUL;
3169 /* End: bn_mp_get_int.c */
3171 /* Start: bn_mp_grow.c */
3174 /* LibTomMath, multiple-precision integer library -- Tom St Denis
3176 * LibTomMath is a library that provides multiple-precision
3177 * integer arithmetic as well as number theoretic functionality.
3179 * The library was designed directly after the MPI library by
3180 * Michael Fromberger but has been written from scratch with
3181 * additional optimizations in place.
3183 * The library is free for all purposes without any express
3184 * guarantee it works.
3186 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
3189 /* grow as required */
3190 int mp_grow (mp_int * a, int size)
3195 /* if the alloc size is smaller alloc more ram */
3196 if (a->alloc < size) {
3197 /* ensure there are always at least MP_PREC digits extra on top */
3198 size += (MP_PREC * 2) - (size % MP_PREC);
3200 /* reallocate the array a->dp
3202 * We store the return in a temporary variable
3203 * in case the operation failed we don't want
3204 * to overwrite the dp member of a.
3206 tmp = OPT_CAST(mp_digit) XREALLOC (a->dp, sizeof (mp_digit) * size);
3208 /* reallocation failed but "a" is still valid [can be freed] */
3212 /* reallocation succeeded so set a->dp */
3215 /* zero excess digits */
3218 for (; i < a->alloc; i++) {
3226 /* End: bn_mp_grow.c */
3228 /* Start: bn_mp_init.c */
3231 /* LibTomMath, multiple-precision integer library -- Tom St Denis
3233 * LibTomMath is a library that provides multiple-precision
3234 * integer arithmetic as well as number theoretic functionality.
3236 * The library was designed directly after the MPI library by
3237 * Michael Fromberger but has been written from scratch with
3238 * additional optimizations in place.
3240 * The library is free for all purposes without any express
3241 * guarantee it works.
3243 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
3246 /* init a new mp_int */
3247 int mp_init (mp_int * a)
3251 /* allocate memory required and clear it */
3252 a->dp = OPT_CAST(mp_digit) XMALLOC (sizeof (mp_digit) * MP_PREC);
3253 if (a->dp == NULL) {
3257 /* set the digits to zero */
3258 for (i = 0; i < MP_PREC; i++) {
3262 /* set the used to zero, allocated digits to the default precision
3263 * and sign to positive */
3272 /* End: bn_mp_init.c */
3274 /* Start: bn_mp_init_copy.c */
3276 #ifdef BN_MP_INIT_COPY_C
3277 /* LibTomMath, multiple-precision integer library -- Tom St Denis
3279 * LibTomMath is a library that provides multiple-precision
3280 * integer arithmetic as well as number theoretic functionality.
3282 * The library was designed directly after the MPI library by
3283 * Michael Fromberger but has been written from scratch with
3284 * additional optimizations in place.
3286 * The library is free for all purposes without any express
3287 * guarantee it works.
3289 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
3292 /* creates "a" then copies b into it */
3293 int mp_init_copy (mp_int * a, mp_int * b)
3297 if ((res = mp_init (a)) != MP_OKAY) {
3300 return mp_copy (b, a);
3304 /* End: bn_mp_init_copy.c */
3306 /* Start: bn_mp_init_multi.c */
3308 #ifdef BN_MP_INIT_MULTI_C
3309 /* LibTomMath, multiple-precision integer library -- Tom St Denis
3311 * LibTomMath is a library that provides multiple-precision
3312 * integer arithmetic as well as number theoretic functionality.
3314 * The library was designed directly after the MPI library by
3315 * Michael Fromberger but has been written from scratch with
3316 * additional optimizations in place.
3318 * The library is free for all purposes without any express
3319 * guarantee it works.
3321 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
3325 int mp_init_multi(mp_int *mp, ...)
3327 mp_err res = MP_OKAY; /* Assume ok until proven otherwise */
3328 int n = 0; /* Number of ok inits */
3329 mp_int* cur_arg = mp;
3332 va_start(args, mp); /* init args to next argument from caller */
3333 while (cur_arg != NULL) {
3334 if (mp_init(cur_arg) != MP_OKAY) {
3335 /* Oops - error! Back-track and mp_clear what we already
3336 succeeded in init-ing, then return error.
3340 /* end the current list */
3343 /* now start cleaning up */
3345 va_start(clean_args, mp);
3348 cur_arg = va_arg(clean_args, mp_int*);
3355 cur_arg = va_arg(args, mp_int*);
3358 return res; /* Assumed ok, if error flagged above. */
3363 /* End: bn_mp_init_multi.c */
3365 /* Start: bn_mp_init_set.c */
3367 #ifdef BN_MP_INIT_SET_C
3368 /* LibTomMath, multiple-precision integer library -- Tom St Denis
3370 * LibTomMath is a library that provides multiple-precision
3371 * integer arithmetic as well as number theoretic functionality.
3373 * The library was designed directly after the MPI library by
3374 * Michael Fromberger but has been written from scratch with
3375 * additional optimizations in place.
3377 * The library is free for all purposes without any express
3378 * guarantee it works.
3380 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
3383 /* initialize and set a digit */
3384 int mp_init_set (mp_int * a, mp_digit b)
3387 if ((err = mp_init(a)) != MP_OKAY) {
3395 /* End: bn_mp_init_set.c */
3397 /* Start: bn_mp_init_set_int.c */
3399 #ifdef BN_MP_INIT_SET_INT_C
3400 /* LibTomMath, multiple-precision integer library -- Tom St Denis
3402 * LibTomMath is a library that provides multiple-precision
3403 * integer arithmetic as well as number theoretic functionality.
3405 * The library was designed directly after the MPI library by
3406 * Michael Fromberger but has been written from scratch with
3407 * additional optimizations in place.
3409 * The library is free for all purposes without any express
3410 * guarantee it works.
3412 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
3415 /* initialize and set a digit */
3416 int mp_init_set_int (mp_int * a, unsigned long b)
3419 if ((err = mp_init(a)) != MP_OKAY) {
3422 return mp_set_int(a, b);
3426 /* End: bn_mp_init_set_int.c */
3428 /* Start: bn_mp_init_size.c */
3430 #ifdef BN_MP_INIT_SIZE_C
3431 /* LibTomMath, multiple-precision integer library -- Tom St Denis
3433 * LibTomMath is a library that provides multiple-precision
3434 * integer arithmetic as well as number theoretic functionality.
3436 * The library was designed directly after the MPI library by
3437 * Michael Fromberger but has been written from scratch with
3438 * additional optimizations in place.
3440 * The library is free for all purposes without any express
3441 * guarantee it works.
3443 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
3446 /* init an mp_init for a given size */
3447 int mp_init_size (mp_int * a, int size)
3451 /* pad size so there are always extra digits */
3452 size += (MP_PREC * 2) - (size % MP_PREC);
3455 a->dp = OPT_CAST(mp_digit) XMALLOC (sizeof (mp_digit) * size);
3456 if (a->dp == NULL) {
3460 /* set the members */
3465 /* zero the digits */
3466 for (x = 0; x < size; x++) {
3474 /* End: bn_mp_init_size.c */
3476 /* Start: bn_mp_invmod.c */
3478 #ifdef BN_MP_INVMOD_C
3479 /* LibTomMath, multiple-precision integer library -- Tom St Denis
3481 * LibTomMath is a library that provides multiple-precision
3482 * integer arithmetic as well as number theoretic functionality.
3484 * The library was designed directly after the MPI library by
3485 * Michael Fromberger but has been written from scratch with
3486 * additional optimizations in place.
3488 * The library is free for all purposes without any express
3489 * guarantee it works.
3491 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
3494 /* hac 14.61, pp608 */
3495 int mp_invmod (mp_int * a, mp_int * b, mp_int * c)
3497 /* b cannot be negative */
3498 if (b->sign == MP_NEG || mp_iszero(b) == 1) {
3502 #ifdef BN_FAST_MP_INVMOD_C
3503 /* if the modulus is odd we can use a faster routine instead */
3504 if (mp_isodd (b) == 1) {
3505 return fast_mp_invmod (a, b, c);
3509 #ifdef BN_MP_INVMOD_SLOW_C
3510 return mp_invmod_slow(a, b, c);
3517 /* End: bn_mp_invmod.c */
3519 /* Start: bn_mp_invmod_slow.c */
3521 #ifdef BN_MP_INVMOD_SLOW_C
3522 /* LibTomMath, multiple-precision integer library -- Tom St Denis
3524 * LibTomMath is a library that provides multiple-precision
3525 * integer arithmetic as well as number theoretic functionality.
3527 * The library was designed directly after the MPI library by
3528 * Michael Fromberger but has been written from scratch with
3529 * additional optimizations in place.
3531 * The library is free for all purposes without any express
3532 * guarantee it works.
3534 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
3537 /* hac 14.61, pp608 */
3538 int mp_invmod_slow (mp_int * a, mp_int * b, mp_int * c)
3540 mp_int x, y, u, v, A, B, C, D;
3543 /* b cannot be negative */
3544 if (b->sign == MP_NEG || mp_iszero(b) == 1) {
3549 if ((res = mp_init_multi(&x, &y, &u, &v,
3550 &A, &B, &C, &D, NULL)) != MP_OKAY) {
3555 if ((res = mp_mod(a, b, &x)) != MP_OKAY) {
3558 if ((res = mp_copy (b, &y)) != MP_OKAY) {
3562 /* 2. [modified] if x,y are both even then return an error! */
3563 if (mp_iseven (&x) == 1 && mp_iseven (&y) == 1) {
3568 /* 3. u=x, v=y, A=1, B=0, C=0,D=1 */
3569 if ((res = mp_copy (&x, &u)) != MP_OKAY) {
3572 if ((res = mp_copy (&y, &v)) != MP_OKAY) {
3579 /* 4. while u is even do */
3580 while (mp_iseven (&u) == 1) {
3582 if ((res = mp_div_2 (&u, &u)) != MP_OKAY) {
3585 /* 4.2 if A or B is odd then */
3586 if (mp_isodd (&A) == 1 || mp_isodd (&B) == 1) {
3587 /* A = (A+y)/2, B = (B-x)/2 */
3588 if ((res = mp_add (&A, &y, &A)) != MP_OKAY) {
3591 if ((res = mp_sub (&B, &x, &B)) != MP_OKAY) {
3595 /* A = A/2, B = B/2 */
3596 if ((res = mp_div_2 (&A, &A)) != MP_OKAY) {
3599 if ((res = mp_div_2 (&B, &B)) != MP_OKAY) {
3604 /* 5. while v is even do */
3605 while (mp_iseven (&v) == 1) {
3607 if ((res = mp_div_2 (&v, &v)) != MP_OKAY) {
3610 /* 5.2 if C or D is odd then */
3611 if (mp_isodd (&C) == 1 || mp_isodd (&D) == 1) {
3612 /* C = (C+y)/2, D = (D-x)/2 */
3613 if ((res = mp_add (&C, &y, &C)) != MP_OKAY) {
3616 if ((res = mp_sub (&D, &x, &D)) != MP_OKAY) {
3620 /* C = C/2, D = D/2 */
3621 if ((res = mp_div_2 (&C, &C)) != MP_OKAY) {
3624 if ((res = mp_div_2 (&D, &D)) != MP_OKAY) {
3629 /* 6. if u >= v then */
3630 if (mp_cmp (&u, &v) != MP_LT) {
3631 /* u = u - v, A = A - C, B = B - D */
3632 if ((res = mp_sub (&u, &v, &u)) != MP_OKAY) {
3636 if ((res = mp_sub (&A, &C, &A)) != MP_OKAY) {
3640 if ((res = mp_sub (&B, &D, &B)) != MP_OKAY) {
3644 /* v - v - u, C = C - A, D = D - B */
3645 if ((res = mp_sub (&v, &u, &v)) != MP_OKAY) {
3649 if ((res = mp_sub (&C, &A, &C)) != MP_OKAY) {
3653 if ((res = mp_sub (&D, &B, &D)) != MP_OKAY) {
3658 /* if not zero goto step 4 */
3659 if (mp_iszero (&u) == 0)
3662 /* now a = C, b = D, gcd == g*v */
3664 /* if v != 1 then there is no inverse */
3665 if (mp_cmp_d (&v, 1) != MP_EQ) {
3670 /* if its too low */
3671 while (mp_cmp_d(&C, 0) == MP_LT) {
3672 if ((res = mp_add(&C, b, &C)) != MP_OKAY) {
3678 while (mp_cmp_mag(&C, b) != MP_LT) {
3679 if ((res = mp_sub(&C, b, &C)) != MP_OKAY) {
3684 /* C is now the inverse */
3687 LBL_ERR:mp_clear_multi (&x, &y, &u, &v, &A, &B, &C, &D, NULL);
3692 /* End: bn_mp_invmod_slow.c */
3694 /* Start: bn_mp_is_square.c */
3696 #ifdef BN_MP_IS_SQUARE_C
3697 /* LibTomMath, multiple-precision integer library -- Tom St Denis
3699 * LibTomMath is a library that provides multiple-precision
3700 * integer arithmetic as well as number theoretic functionality.
3702 * The library was designed directly after the MPI library by
3703 * Michael Fromberger but has been written from scratch with
3704 * additional optimizations in place.
3706 * The library is free for all purposes without any express
3707 * guarantee it works.
3709 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
3712 /* Check if remainders are possible squares - fast exclude non-squares */
3713 static const char rem_128[128] = {
3714 0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
3715 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
3716 1, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
3717 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
3718 0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
3719 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
3720 1, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
3721 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1
3724 static const char rem_105[105] = {
3725 0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1,
3726 0, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1,
3727 0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1,
3728 1, 0, 1, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1,
3729 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1,
3730 1, 1, 1, 1, 0, 1, 0, 1, 1, 0, 0, 1, 1, 1, 1,
3731 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1
3734 /* Store non-zero to ret if arg is square, and zero if not */
3735 int mp_is_square(mp_int *arg,int *ret)
3742 /* Default to Non-square :) */
3745 if (arg->sign == MP_NEG) {
3749 /* digits used? (TSD) */
3750 if (arg->used == 0) {
3754 /* First check mod 128 (suppose that DIGIT_BIT is at least 7) */
3755 if (rem_128[127 & DIGIT(arg,0)] == 1) {
3759 /* Next check mod 105 (3*5*7) */
3760 if ((res = mp_mod_d(arg,105,&c)) != MP_OKAY) {
3763 if (rem_105[c] == 1) {
3768 if ((res = mp_init_set_int(&t,11L*13L*17L*19L*23L*29L*31L)) != MP_OKAY) {
3771 if ((res = mp_mod(arg,&t,&t)) != MP_OKAY) {
3775 /* Check for other prime modules, note it's not an ERROR but we must
3776 * free "t" so the easiest way is to goto ERR. We know that res
3777 * is already equal to MP_OKAY from the mp_mod call
3779 if ( (1L<<(r%11)) & 0x5C4L ) goto ERR;
3780 if ( (1L<<(r%13)) & 0x9E4L ) goto ERR;
3781 if ( (1L<<(r%17)) & 0x5CE8L ) goto ERR;
3782 if ( (1L<<(r%19)) & 0x4F50CL ) goto ERR;
3783 if ( (1L<<(r%23)) & 0x7ACCA0L ) goto ERR;
3784 if ( (1L<<(r%29)) & 0xC2EDD0CL ) goto ERR;
3785 if ( (1L<<(r%31)) & 0x6DE2B848L ) goto ERR;
3787 /* Final check - is sqr(sqrt(arg)) == arg ? */
3788 if ((res = mp_sqrt(arg,&t)) != MP_OKAY) {
3791 if ((res = mp_sqr(&t,&t)) != MP_OKAY) {
3795 *ret = (mp_cmp_mag(&t,arg) == MP_EQ) ? MP_YES : MP_NO;
3801 /* End: bn_mp_is_square.c */
3803 /* Start: bn_mp_jacobi.c */
3805 #ifdef BN_MP_JACOBI_C
3806 /* LibTomMath, multiple-precision integer library -- Tom St Denis
3808 * LibTomMath is a library that provides multiple-precision
3809 * integer arithmetic as well as number theoretic functionality.
3811 * The library was designed directly after the MPI library by
3812 * Michael Fromberger but has been written from scratch with
3813 * additional optimizations in place.
3815 * The library is free for all purposes without any express
3816 * guarantee it works.
3818 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
3821 /* computes the jacobi c = (a | n) (or Legendre if n is prime)
3822 * HAC pp. 73 Algorithm 2.149
3824 int mp_jacobi (mp_int * a, mp_int * p, int *c)
3830 /* if p <= 0 return MP_VAL */
3831 if (mp_cmp_d(p, 0) != MP_GT) {
3835 /* step 1. if a == 0, return 0 */
3836 if (mp_iszero (a) == 1) {
3841 /* step 2. if a == 1, return 1 */
3842 if (mp_cmp_d (a, 1) == MP_EQ) {
3850 /* step 3. write a = a1 * 2**k */
3851 if ((res = mp_init_copy (&a1, a)) != MP_OKAY) {
3855 if ((res = mp_init (&p1)) != MP_OKAY) {
3859 /* divide out larger power of two */
3860 k = mp_cnt_lsb(&a1);
3861 if ((res = mp_div_2d(&a1, k, &a1, NULL)) != MP_OKAY) {
3865 /* step 4. if e is even set s=1 */
3869 /* else set s=1 if p = 1/7 (mod 8) or s=-1 if p = 3/5 (mod 8) */
3870 residue = p->dp[0] & 7;
3872 if (residue == 1 || residue == 7) {
3874 } else if (residue == 3 || residue == 5) {
3879 /* step 5. if p == 3 (mod 4) *and* a1 == 3 (mod 4) then s = -s */
3880 if ( ((p->dp[0] & 3) == 3) && ((a1.dp[0] & 3) == 3)) {
3884 /* if a1 == 1 we're done */
3885 if (mp_cmp_d (&a1, 1) == MP_EQ) {
3889 if ((res = mp_mod (p, &a1, &p1)) != MP_OKAY) {
3892 if ((res = mp_jacobi (&p1, &a1, &r)) != MP_OKAY) {
3900 LBL_P1:mp_clear (&p1);
3901 LBL_A1:mp_clear (&a1);
3906 /* End: bn_mp_jacobi.c */
3908 /* Start: bn_mp_karatsuba_mul.c */
3910 #ifdef BN_MP_KARATSUBA_MUL_C
3911 /* LibTomMath, multiple-precision integer library -- Tom St Denis
3913 * LibTomMath is a library that provides multiple-precision
3914 * integer arithmetic as well as number theoretic functionality.
3916 * The library was designed directly after the MPI library by
3917 * Michael Fromberger but has been written from scratch with
3918 * additional optimizations in place.
3920 * The library is free for all purposes without any express
3921 * guarantee it works.
3923 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
3926 /* c = |a| * |b| using Karatsuba Multiplication using
3927 * three half size multiplications
3929 * Let B represent the radix [e.g. 2**DIGIT_BIT] and
3930 * let n represent half of the number of digits in
3933 * a = a1 * B**n + a0
3934 * b = b1 * B**n + b0
3937 a1b1 * B**2n + ((a1 - a0)(b1 - b0) + a0b0 + a1b1) * B + a0b0
3939 * Note that a1b1 and a0b0 are used twice and only need to be
3940 * computed once. So in total three half size (half # of
3941 * digit) multiplications are performed, a0b0, a1b1 and
3944 * Note that a multiplication of half the digits requires
3945 * 1/4th the number of single precision multiplications so in
3946 * total after one call 25% of the single precision multiplications
3947 * are saved. Note also that the call to mp_mul can end up back
3948 * in this function if the a0, a1, b0, or b1 are above the threshold.
3949 * This is known as divide-and-conquer and leads to the famous
3950 * O(N**lg(3)) or O(N**1.584) work which is asymptopically lower than
3951 * the standard O(N**2) that the baseline/comba methods use.
3952 * Generally though the overhead of this method doesn't pay off
3953 * until a certain size (N ~ 80) is reached.
3955 int mp_karatsuba_mul (mp_int * a, mp_int * b, mp_int * c)
3957 mp_int x0, x1, y0, y1, t1, x0y0, x1y1;
3960 /* default the return code to an error */
3963 /* min # of digits */
3964 B = MIN (a->used, b->used);
3966 /* now divide in two */
3969 /* init copy all the temps */
3970 if (mp_init_size (&x0, B) != MP_OKAY)
3972 if (mp_init_size (&x1, a->used - B) != MP_OKAY)
3974 if (mp_init_size (&y0, B) != MP_OKAY)
3976 if (mp_init_size (&y1, b->used - B) != MP_OKAY)
3980 if (mp_init_size (&t1, B * 2) != MP_OKAY)
3982 if (mp_init_size (&x0y0, B * 2) != MP_OKAY)
3984 if (mp_init_size (&x1y1, B * 2) != MP_OKAY)
3987 /* now shift the digits */
3988 x0.used = y0.used = B;
3989 x1.used = a->used - B;
3990 y1.used = b->used - B;
3994 register mp_digit *tmpa, *tmpb, *tmpx, *tmpy;
3996 /* we copy the digits directly instead of using higher level functions
3997 * since we also need to shift the digits
4004 for (x = 0; x < B; x++) {
4010 for (x = B; x < a->used; x++) {
4015 for (x = B; x < b->used; x++) {
4020 /* only need to clamp the lower words since by definition the
4021 * upper words x1/y1 must have a known number of digits
4026 /* now calc the products x0y0 and x1y1 */
4027 /* after this x0 is no longer required, free temp [x0==t2]! */
4028 if (mp_mul (&x0, &y0, &x0y0) != MP_OKAY)
4029 goto X1Y1; /* x0y0 = x0*y0 */
4030 if (mp_mul (&x1, &y1, &x1y1) != MP_OKAY)
4031 goto X1Y1; /* x1y1 = x1*y1 */
4033 /* now calc x1-x0 and y1-y0 */
4034 if (mp_sub (&x1, &x0, &t1) != MP_OKAY)
4035 goto X1Y1; /* t1 = x1 - x0 */
4036 if (mp_sub (&y1, &y0, &x0) != MP_OKAY)
4037 goto X1Y1; /* t2 = y1 - y0 */
4038 if (mp_mul (&t1, &x0, &t1) != MP_OKAY)
4039 goto X1Y1; /* t1 = (x1 - x0) * (y1 - y0) */
4042 if (mp_add (&x0y0, &x1y1, &x0) != MP_OKAY)
4043 goto X1Y1; /* t2 = x0y0 + x1y1 */
4044 if (mp_sub (&x0, &t1, &t1) != MP_OKAY)
4045 goto X1Y1; /* t1 = x0y0 + x1y1 - (x1-x0)*(y1-y0) */
4048 if (mp_lshd (&t1, B) != MP_OKAY)
4049 goto X1Y1; /* t1 = (x0y0 + x1y1 - (x1-x0)*(y1-y0))<<B */
4050 if (mp_lshd (&x1y1, B * 2) != MP_OKAY)
4051 goto X1Y1; /* x1y1 = x1y1 << 2*B */
4053 if (mp_add (&x0y0, &t1, &t1) != MP_OKAY)
4054 goto X1Y1; /* t1 = x0y0 + t1 */
4055 if (mp_add (&t1, &x1y1, c) != MP_OKAY)
4056 goto X1Y1; /* t1 = x0y0 + t1 + x1y1 */
4058 /* Algorithm succeeded set the return code to MP_OKAY */
4061 X1Y1:mp_clear (&x1y1);
4062 X0Y0:mp_clear (&x0y0);
4073 /* End: bn_mp_karatsuba_mul.c */
4075 /* Start: bn_mp_karatsuba_sqr.c */
4077 #ifdef BN_MP_KARATSUBA_SQR_C
4078 /* LibTomMath, multiple-precision integer library -- Tom St Denis
4080 * LibTomMath is a library that provides multiple-precision
4081 * integer arithmetic as well as number theoretic functionality.
4083 * The library was designed directly after the MPI library by
4084 * Michael Fromberger but has been written from scratch with
4085 * additional optimizations in place.
4087 * The library is free for all purposes without any express
4088 * guarantee it works.
4090 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
4093 /* Karatsuba squaring, computes b = a*a using three
4094 * half size squarings
4096 * See comments of karatsuba_mul for details. It
4097 * is essentially the same algorithm but merely
4098 * tuned to perform recursive squarings.
4100 int mp_karatsuba_sqr (mp_int * a, mp_int * b)
4102 mp_int x0, x1, t1, t2, x0x0, x1x1;
4107 /* min # of digits */
4110 /* now divide in two */
4113 /* init copy all the temps */
4114 if (mp_init_size (&x0, B) != MP_OKAY)
4116 if (mp_init_size (&x1, a->used - B) != MP_OKAY)
4120 if (mp_init_size (&t1, a->used * 2) != MP_OKAY)
4122 if (mp_init_size (&t2, a->used * 2) != MP_OKAY)
4124 if (mp_init_size (&x0x0, B * 2) != MP_OKAY)
4126 if (mp_init_size (&x1x1, (a->used - B) * 2) != MP_OKAY)
4131 register mp_digit *dst, *src;
4135 /* now shift the digits */
4137 for (x = 0; x < B; x++) {
4142 for (x = B; x < a->used; x++) {
4148 x1.used = a->used - B;
4152 /* now calc the products x0*x0 and x1*x1 */
4153 if (mp_sqr (&x0, &x0x0) != MP_OKAY)
4154 goto X1X1; /* x0x0 = x0*x0 */
4155 if (mp_sqr (&x1, &x1x1) != MP_OKAY)
4156 goto X1X1; /* x1x1 = x1*x1 */
4158 /* now calc (x1-x0)**2 */
4159 if (mp_sub (&x1, &x0, &t1) != MP_OKAY)
4160 goto X1X1; /* t1 = x1 - x0 */
4161 if (mp_sqr (&t1, &t1) != MP_OKAY)
4162 goto X1X1; /* t1 = (x1 - x0) * (x1 - x0) */
4165 if (s_mp_add (&x0x0, &x1x1, &t2) != MP_OKAY)
4166 goto X1X1; /* t2 = x0x0 + x1x1 */
4167 if (mp_sub (&t2, &t1, &t1) != MP_OKAY)
4168 goto X1X1; /* t1 = x0x0 + x1x1 - (x1-x0)*(x1-x0) */
4171 if (mp_lshd (&t1, B) != MP_OKAY)
4172 goto X1X1; /* t1 = (x0x0 + x1x1 - (x1-x0)*(x1-x0))<<B */
4173 if (mp_lshd (&x1x1, B * 2) != MP_OKAY)
4174 goto X1X1; /* x1x1 = x1x1 << 2*B */
4176 if (mp_add (&x0x0, &t1, &t1) != MP_OKAY)
4177 goto X1X1; /* t1 = x0x0 + t1 */
4178 if (mp_add (&t1, &x1x1, b) != MP_OKAY)
4179 goto X1X1; /* t1 = x0x0 + t1 + x1x1 */
4183 X1X1:mp_clear (&x1x1);
4184 X0X0:mp_clear (&x0x0);
4194 /* End: bn_mp_karatsuba_sqr.c */
4196 /* Start: bn_mp_lcm.c */
4199 /* LibTomMath, multiple-precision integer library -- Tom St Denis
4201 * LibTomMath is a library that provides multiple-precision
4202 * integer arithmetic as well as number theoretic functionality.
4204 * The library was designed directly after the MPI library by
4205 * Michael Fromberger but has been written from scratch with
4206 * additional optimizations in place.
4208 * The library is free for all purposes without any express
4209 * guarantee it works.
4211 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
4214 /* computes least common multiple as |a*b|/(a, b) */
4215 int mp_lcm (mp_int * a, mp_int * b, mp_int * c)
4221 if ((res = mp_init_multi (&t1, &t2, NULL)) != MP_OKAY) {
4225 /* t1 = get the GCD of the two inputs */
4226 if ((res = mp_gcd (a, b, &t1)) != MP_OKAY) {
4230 /* divide the smallest by the GCD */
4231 if (mp_cmp_mag(a, b) == MP_LT) {
4232 /* store quotient in t2 such that t2 * b is the LCM */
4233 if ((res = mp_div(a, &t1, &t2, NULL)) != MP_OKAY) {
4236 res = mp_mul(b, &t2, c);
4238 /* store quotient in t2 such that t2 * a is the LCM */
4239 if ((res = mp_div(b, &t1, &t2, NULL)) != MP_OKAY) {
4242 res = mp_mul(a, &t2, c);
4245 /* fix the sign to positive */
4249 mp_clear_multi (&t1, &t2, NULL);
4254 /* End: bn_mp_lcm.c */
4256 /* Start: bn_mp_lshd.c */
4259 /* LibTomMath, multiple-precision integer library -- Tom St Denis
4261 * LibTomMath is a library that provides multiple-precision
4262 * integer arithmetic as well as number theoretic functionality.
4264 * The library was designed directly after the MPI library by
4265 * Michael Fromberger but has been written from scratch with
4266 * additional optimizations in place.
4268 * The library is free for all purposes without any express
4269 * guarantee it works.
4271 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
4274 /* shift left a certain amount of digits */
4275 int mp_lshd (mp_int * a, int b)
4279 /* if its less than zero return */
4284 /* grow to fit the new digits */
4285 if (a->alloc < a->used + b) {
4286 if ((res = mp_grow (a, a->used + b)) != MP_OKAY) {
4292 register mp_digit *top, *bottom;
4294 /* increment the used by the shift amount then copy upwards */
4298 top = a->dp + a->used - 1;
4301 bottom = a->dp + a->used - 1 - b;
4303 /* much like mp_rshd this is implemented using a sliding window
4304 * except the window goes the otherway around. Copying from
4305 * the bottom to the top. see bn_mp_rshd.c for more info.
4307 for (x = a->used - 1; x >= b; x--) {
4311 /* zero the lower digits */
4313 for (x = 0; x < b; x++) {
4321 /* End: bn_mp_lshd.c */
4323 /* Start: bn_mp_mod.c */
4326 /* LibTomMath, multiple-precision integer library -- Tom St Denis
4328 * LibTomMath is a library that provides multiple-precision
4329 * integer arithmetic as well as number theoretic functionality.
4331 * The library was designed directly after the MPI library by
4332 * Michael Fromberger but has been written from scratch with
4333 * additional optimizations in place.
4335 * The library is free for all purposes without any express
4336 * guarantee it works.
4338 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
4341 /* c = a mod b, 0 <= c < b */
4343 mp_mod (mp_int * a, mp_int * b, mp_int * c)
4348 if ((res = mp_init (&t)) != MP_OKAY) {
4352 if ((res = mp_div (a, b, NULL, &t)) != MP_OKAY) {
4357 if (t.sign != b->sign) {
4358 res = mp_add (b, &t, c);
4369 /* End: bn_mp_mod.c */
4371 /* Start: bn_mp_mod_2d.c */
4373 #ifdef BN_MP_MOD_2D_C
4374 /* LibTomMath, multiple-precision integer library -- Tom St Denis
4376 * LibTomMath is a library that provides multiple-precision
4377 * integer arithmetic as well as number theoretic functionality.
4379 * The library was designed directly after the MPI library by
4380 * Michael Fromberger but has been written from scratch with
4381 * additional optimizations in place.
4383 * The library is free for all purposes without any express
4384 * guarantee it works.
4386 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
4389 /* calc a value mod 2**b */
4391 mp_mod_2d (mp_int * a, int b, mp_int * c)
4395 /* if b is <= 0 then zero the int */
4401 /* if the modulus is larger than the value than return */
4402 if (b >= (int) (a->used * DIGIT_BIT)) {
4403 res = mp_copy (a, c);
4408 if ((res = mp_copy (a, c)) != MP_OKAY) {
4412 /* zero digits above the last digit of the modulus */
4413 for (x = (b / DIGIT_BIT) + ((b % DIGIT_BIT) == 0 ? 0 : 1); x < c->used; x++) {
4416 /* clear the digit that is not completely outside/inside the modulus */
4417 c->dp[b / DIGIT_BIT] &=
4418 (mp_digit) ((((mp_digit) 1) << (((mp_digit) b) % DIGIT_BIT)) - ((mp_digit) 1));
4424 /* End: bn_mp_mod_2d.c */
4426 /* Start: bn_mp_mod_d.c */
4428 #ifdef BN_MP_MOD_D_C
4429 /* LibTomMath, multiple-precision integer library -- Tom St Denis
4431 * LibTomMath is a library that provides multiple-precision
4432 * integer arithmetic as well as number theoretic functionality.
4434 * The library was designed directly after the MPI library by
4435 * Michael Fromberger but has been written from scratch with
4436 * additional optimizations in place.
4438 * The library is free for all purposes without any express
4439 * guarantee it works.
4441 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
4445 mp_mod_d (mp_int * a, mp_digit b, mp_digit * c)
4447 return mp_div_d(a, b, NULL, c);
4451 /* End: bn_mp_mod_d.c */
4453 /* Start: bn_mp_montgomery_calc_normalization.c */
4455 #ifdef BN_MP_MONTGOMERY_CALC_NORMALIZATION_C
4456 /* LibTomMath, multiple-precision integer library -- Tom St Denis
4458 * LibTomMath is a library that provides multiple-precision
4459 * integer arithmetic as well as number theoretic functionality.
4461 * The library was designed directly after the MPI library by
4462 * Michael Fromberger but has been written from scratch with
4463 * additional optimizations in place.
4465 * The library is free for all purposes without any express
4466 * guarantee it works.
4468 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
4472 * shifts with subtractions when the result is greater than b.
4474 * The method is slightly modified to shift B unconditionally upto just under
4475 * the leading bit of b. This saves alot of multiple precision shifting.
4477 int mp_montgomery_calc_normalization (mp_int * a, mp_int * b)
4481 /* how many bits of last digit does b use */
4482 bits = mp_count_bits (b) % DIGIT_BIT;
4485 if ((res = mp_2expt (a, (b->used - 1) * DIGIT_BIT + bits - 1)) != MP_OKAY) {
4494 /* now compute C = A * B mod b */
4495 for (x = bits - 1; x < (int)DIGIT_BIT; x++) {
4496 if ((res = mp_mul_2 (a, a)) != MP_OKAY) {
4499 if (mp_cmp_mag (a, b) != MP_LT) {
4500 if ((res = s_mp_sub (a, b, a)) != MP_OKAY) {
4510 /* End: bn_mp_montgomery_calc_normalization.c */
4512 /* Start: bn_mp_montgomery_reduce.c */
4514 #ifdef BN_MP_MONTGOMERY_REDUCE_C
4515 /* LibTomMath, multiple-precision integer library -- Tom St Denis
4517 * LibTomMath is a library that provides multiple-precision
4518 * integer arithmetic as well as number theoretic functionality.
4520 * The library was designed directly after the MPI library by
4521 * Michael Fromberger but has been written from scratch with
4522 * additional optimizations in place.
4524 * The library is free for all purposes without any express
4525 * guarantee it works.
4527 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
4530 /* computes xR**-1 == x (mod N) via Montgomery Reduction */
4532 mp_montgomery_reduce (mp_int * x, mp_int * n, mp_digit rho)
4537 /* can the fast reduction [comba] method be used?
4539 * Note that unlike in mul you're safely allowed *less*
4540 * than the available columns [255 per default] since carries
4541 * are fixed up in the inner loop.
4543 digs = n->used * 2 + 1;
4544 if ((digs < MP_WARRAY) &&
4546 (1 << ((CHAR_BIT * sizeof (mp_word)) - (2 * DIGIT_BIT)))) {
4547 return fast_mp_montgomery_reduce (x, n, rho);
4550 /* grow the input as required */
4551 if (x->alloc < digs) {
4552 if ((res = mp_grow (x, digs)) != MP_OKAY) {
4558 for (ix = 0; ix < n->used; ix++) {
4559 /* mu = ai * rho mod b
4561 * The value of rho must be precalculated via
4562 * montgomery_setup() such that
4563 * it equals -1/n0 mod b this allows the
4564 * following inner loop to reduce the
4565 * input one digit at a time
4567 mu = (mp_digit) (((mp_word)x->dp[ix]) * ((mp_word)rho) & MP_MASK);
4569 /* a = a + mu * m * b**i */
4572 register mp_digit *tmpn, *tmpx, u;
4575 /* alias for digits of the modulus */
4578 /* alias for the digits of x [the input] */
4581 /* set the carry to zero */
4584 /* Multiply and add in place */
4585 for (iy = 0; iy < n->used; iy++) {
4586 /* compute product and sum */
4587 r = ((mp_word)mu) * ((mp_word)*tmpn++) +
4588 ((mp_word) u) + ((mp_word) * tmpx);
4591 u = (mp_digit)(r >> ((mp_word) DIGIT_BIT));
4594 *tmpx++ = (mp_digit)(r & ((mp_word) MP_MASK));
4596 /* At this point the ix'th digit of x should be zero */
4599 /* propagate carries upwards as required*/
4602 u = *tmpx >> DIGIT_BIT;
4608 /* at this point the n.used'th least
4609 * significant digits of x are all zero
4610 * which means we can shift x to the
4611 * right by n.used digits and the
4612 * residue is unchanged.
4615 /* x = x/b**n.used */
4617 mp_rshd (x, n->used);
4619 /* if x >= n then x = x - n */
4620 if (mp_cmp_mag (x, n) != MP_LT) {
4621 return s_mp_sub (x, n, x);
4628 /* End: bn_mp_montgomery_reduce.c */
4630 /* Start: bn_mp_montgomery_setup.c */
4632 #ifdef BN_MP_MONTGOMERY_SETUP_C
4633 /* LibTomMath, multiple-precision integer library -- Tom St Denis
4635 * LibTomMath is a library that provides multiple-precision
4636 * integer arithmetic as well as number theoretic functionality.
4638 * The library was designed directly after the MPI library by
4639 * Michael Fromberger but has been written from scratch with
4640 * additional optimizations in place.
4642 * The library is free for all purposes without any express
4643 * guarantee it works.
4645 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
4648 /* setups the montgomery reduction stuff */
4650 mp_montgomery_setup (mp_int * n, mp_digit * rho)
4654 /* fast inversion mod 2**k
4656 * Based on the fact that
4658 * XA = 1 (mod 2**n) => (X(2-XA)) A = 1 (mod 2**2n)
4659 * => 2*X*A - X*X*A*A = 1
4660 * => 2*(1) - (1) = 1
4668 x = (((b + 2) & 4) << 1) + b; /* here x*a==1 mod 2**4 */
4669 x *= 2 - b * x; /* here x*a==1 mod 2**8 */
4670 #if !defined(MP_8BIT)
4671 x *= 2 - b * x; /* here x*a==1 mod 2**16 */
4673 #if defined(MP_64BIT) || !(defined(MP_8BIT) || defined(MP_16BIT))
4674 x *= 2 - b * x; /* here x*a==1 mod 2**32 */
4677 x *= 2 - b * x; /* here x*a==1 mod 2**64 */
4680 /* rho = -1/m mod b */
4681 *rho = (((mp_word)1 << ((mp_word) DIGIT_BIT)) - x) & MP_MASK;
4687 /* End: bn_mp_montgomery_setup.c */
4689 /* Start: bn_mp_mul.c */
4692 /* LibTomMath, multiple-precision integer library -- Tom St Denis
4694 * LibTomMath is a library that provides multiple-precision
4695 * integer arithmetic as well as number theoretic functionality.
4697 * The library was designed directly after the MPI library by
4698 * Michael Fromberger but has been written from scratch with
4699 * additional optimizations in place.
4701 * The library is free for all purposes without any express
4702 * guarantee it works.
4704 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
4707 /* high level multiplication (handles sign) */
4708 int mp_mul (mp_int * a, mp_int * b, mp_int * c)
4711 neg = (a->sign == b->sign) ? MP_ZPOS : MP_NEG;
4713 /* use Toom-Cook? */
4714 #ifdef BN_MP_TOOM_MUL_C
4715 if (MIN (a->used, b->used) >= TOOM_MUL_CUTOFF) {
4716 res = mp_toom_mul(a, b, c);
4719 #ifdef BN_MP_KARATSUBA_MUL_C
4720 /* use Karatsuba? */
4721 if (MIN (a->used, b->used) >= KARATSUBA_MUL_CUTOFF) {
4722 res = mp_karatsuba_mul (a, b, c);
4726 /* can we use the fast multiplier?
4728 * The fast multiplier can be used if the output will
4729 * have less than MP_WARRAY digits and the number of
4730 * digits won't affect carry propagation
4732 int digs = a->used + b->used + 1;
4734 #ifdef BN_FAST_S_MP_MUL_DIGS_C
4735 if ((digs < MP_WARRAY) &&
4736 MIN(a->used, b->used) <=
4737 (1 << ((CHAR_BIT * sizeof (mp_word)) - (2 * DIGIT_BIT)))) {
4738 res = fast_s_mp_mul_digs (a, b, c, digs);
4741 #ifdef BN_S_MP_MUL_DIGS_C
4742 res = s_mp_mul (a, b, c); /* uses s_mp_mul_digs */
4748 c->sign = (c->used > 0) ? neg : MP_ZPOS;
4753 /* End: bn_mp_mul.c */
4755 /* Start: bn_mp_mul_2.c */
4757 #ifdef BN_MP_MUL_2_C
4758 /* LibTomMath, multiple-precision integer library -- Tom St Denis
4760 * LibTomMath is a library that provides multiple-precision
4761 * integer arithmetic as well as number theoretic functionality.
4763 * The library was designed directly after the MPI library by
4764 * Michael Fromberger but has been written from scratch with
4765 * additional optimizations in place.
4767 * The library is free for all purposes without any express
4768 * guarantee it works.
4770 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
4774 int mp_mul_2(mp_int * a, mp_int * b)
4776 int x, res, oldused;
4778 /* grow to accomodate result */
4779 if (b->alloc < a->used + 1) {
4780 if ((res = mp_grow (b, a->used + 1)) != MP_OKAY) {
4789 register mp_digit r, rr, *tmpa, *tmpb;
4791 /* alias for source */
4794 /* alias for dest */
4799 for (x = 0; x < a->used; x++) {
4801 /* get what will be the *next* carry bit from the
4802 * MSB of the current digit
4804 rr = *tmpa >> ((mp_digit)(DIGIT_BIT - 1));
4806 /* now shift up this digit, add in the carry [from the previous] */
4807 *tmpb++ = ((*tmpa++ << ((mp_digit)1)) | r) & MP_MASK;
4809 /* copy the carry that would be from the source
4810 * digit into the next iteration
4815 /* new leading digit? */
4817 /* add a MSB which is always 1 at this point */
4822 /* now zero any excess digits on the destination
4823 * that we didn't write to
4825 tmpb = b->dp + b->used;
4826 for (x = b->used; x < oldused; x++) {
4835 /* End: bn_mp_mul_2.c */
4837 /* Start: bn_mp_mul_2d.c */
4839 #ifdef BN_MP_MUL_2D_C
4840 /* LibTomMath, multiple-precision integer library -- Tom St Denis
4842 * LibTomMath is a library that provides multiple-precision
4843 * integer arithmetic as well as number theoretic functionality.
4845 * The library was designed directly after the MPI library by
4846 * Michael Fromberger but has been written from scratch with
4847 * additional optimizations in place.
4849 * The library is free for all purposes without any express
4850 * guarantee it works.
4852 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
4855 /* shift left by a certain bit count */
4856 int mp_mul_2d (mp_int * a, int b, mp_int * c)
4863 if ((res = mp_copy (a, c)) != MP_OKAY) {
4868 if (c->alloc < (int)(c->used + b/DIGIT_BIT + 1)) {
4869 if ((res = mp_grow (c, c->used + b / DIGIT_BIT + 1)) != MP_OKAY) {
4874 /* shift by as many digits in the bit count */
4875 if (b >= (int)DIGIT_BIT) {
4876 if ((res = mp_lshd (c, b / DIGIT_BIT)) != MP_OKAY) {
4881 /* shift any bit count < DIGIT_BIT */
4882 d = (mp_digit) (b % DIGIT_BIT);
4884 register mp_digit *tmpc, shift, mask, r, rr;
4887 /* bitmask for carries */
4888 mask = (((mp_digit)1) << d) - 1;
4890 /* shift for msbs */
4891 shift = DIGIT_BIT - d;
4898 for (x = 0; x < c->used; x++) {
4899 /* get the higher bits of the current word */
4900 rr = (*tmpc >> shift) & mask;
4902 /* shift the current word and OR in the carry */
4903 *tmpc = ((*tmpc << d) | r) & MP_MASK;
4906 /* set the carry to the carry bits of the current word */
4910 /* set final carry */
4912 c->dp[(c->used)++] = r;
4920 /* End: bn_mp_mul_2d.c */
4922 /* Start: bn_mp_mul_d.c */
4924 #ifdef BN_MP_MUL_D_C
4925 /* LibTomMath, multiple-precision integer library -- Tom St Denis
4927 * LibTomMath is a library that provides multiple-precision
4928 * integer arithmetic as well as number theoretic functionality.
4930 * The library was designed directly after the MPI library by
4931 * Michael Fromberger but has been written from scratch with
4932 * additional optimizations in place.
4934 * The library is free for all purposes without any express
4935 * guarantee it works.
4937 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
4940 /* multiply by a digit */
4942 mp_mul_d (mp_int * a, mp_digit b, mp_int * c)
4944 mp_digit u, *tmpa, *tmpc;
4946 int ix, res, olduse;
4948 /* make sure c is big enough to hold a*b */
4949 if (c->alloc < a->used + 1) {
4950 if ((res = mp_grow (c, a->used + 1)) != MP_OKAY) {
4955 /* get the original destinations used count */
4961 /* alias for a->dp [source] */
4964 /* alias for c->dp [dest] */
4970 /* compute columns */
4971 for (ix = 0; ix < a->used; ix++) {
4972 /* compute product and carry sum for this term */
4973 r = ((mp_word) u) + ((mp_word)*tmpa++) * ((mp_word)b);
4975 /* mask off higher bits to get a single digit */
4976 *tmpc++ = (mp_digit) (r & ((mp_word) MP_MASK));
4978 /* send carry into next iteration */
4979 u = (mp_digit) (r >> ((mp_word) DIGIT_BIT));
4982 /* store final carry [if any] and increment ix offset */
4986 /* now zero digits above the top */
4987 while (ix++ < olduse) {
4991 /* set used count */
4992 c->used = a->used + 1;
4999 /* End: bn_mp_mul_d.c */
5001 /* Start: bn_mp_mulmod.c */
5003 #ifdef BN_MP_MULMOD_C
5004 /* LibTomMath, multiple-precision integer library -- Tom St Denis
5006 * LibTomMath is a library that provides multiple-precision
5007 * integer arithmetic as well as number theoretic functionality.
5009 * The library was designed directly after the MPI library by
5010 * Michael Fromberger but has been written from scratch with
5011 * additional optimizations in place.
5013 * The library is free for all purposes without any express
5014 * guarantee it works.
5016 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
5019 /* d = a * b (mod c) */
5021 mp_mulmod (mp_int * a, mp_int * b, mp_int * c, mp_int * d)
5026 if ((res = mp_init (&t)) != MP_OKAY) {
5030 if ((res = mp_mul (a, b, &t)) != MP_OKAY) {
5034 res = mp_mod (&t, c, d);
5040 /* End: bn_mp_mulmod.c */
5042 /* Start: bn_mp_n_root.c */
5044 #ifdef BN_MP_N_ROOT_C
5045 /* LibTomMath, multiple-precision integer library -- Tom St Denis
5047 * LibTomMath is a library that provides multiple-precision
5048 * integer arithmetic as well as number theoretic functionality.
5050 * The library was designed directly after the MPI library by
5051 * Michael Fromberger but has been written from scratch with
5052 * additional optimizations in place.
5054 * The library is free for all purposes without any express
5055 * guarantee it works.
5057 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
5060 /* find the n'th root of an integer
5062 * Result found such that (c)**b <= a and (c+1)**b > a
5064 * This algorithm uses Newton's approximation
5065 * x[i+1] = x[i] - f(x[i])/f'(x[i])
5066 * which will find the root in log(N) time where
5067 * each step involves a fair bit. This is not meant to
5068 * find huge roots [square and cube, etc].
5070 int mp_n_root (mp_int * a, mp_digit b, mp_int * c)
5075 /* input must be positive if b is even */
5076 if ((b & 1) == 0 && a->sign == MP_NEG) {
5080 if ((res = mp_init (&t1)) != MP_OKAY) {
5084 if ((res = mp_init (&t2)) != MP_OKAY) {
5088 if ((res = mp_init (&t3)) != MP_OKAY) {
5092 /* if a is negative fudge the sign but keep track */
5101 if ((res = mp_copy (&t2, &t1)) != MP_OKAY) {
5105 /* t2 = t1 - ((t1**b - a) / (b * t1**(b-1))) */
5107 /* t3 = t1**(b-1) */
5108 if ((res = mp_expt_d (&t1, b - 1, &t3)) != MP_OKAY) {
5114 if ((res = mp_mul (&t3, &t1, &t2)) != MP_OKAY) {
5118 /* t2 = t1**b - a */
5119 if ((res = mp_sub (&t2, a, &t2)) != MP_OKAY) {
5124 /* t3 = t1**(b-1) * b */
5125 if ((res = mp_mul_d (&t3, b, &t3)) != MP_OKAY) {
5129 /* t3 = (t1**b - a)/(b * t1**(b-1)) */
5130 if ((res = mp_div (&t2, &t3, &t3, NULL)) != MP_OKAY) {
5134 if ((res = mp_sub (&t1, &t3, &t2)) != MP_OKAY) {
5137 } while (mp_cmp (&t1, &t2) != MP_EQ);
5139 /* result can be off by a few so check */
5141 if ((res = mp_expt_d (&t1, b, &t2)) != MP_OKAY) {
5145 if (mp_cmp (&t2, a) == MP_GT) {
5146 if ((res = mp_sub_d (&t1, 1, &t1)) != MP_OKAY) {
5154 /* reset the sign of a first */
5157 /* set the result */
5160 /* set the sign of the result */
5165 LBL_T3:mp_clear (&t3);
5166 LBL_T2:mp_clear (&t2);
5167 LBL_T1:mp_clear (&t1);
5172 /* End: bn_mp_n_root.c */
5174 /* Start: bn_mp_neg.c */
5177 /* LibTomMath, multiple-precision integer library -- Tom St Denis
5179 * LibTomMath is a library that provides multiple-precision
5180 * integer arithmetic as well as number theoretic functionality.
5182 * The library was designed directly after the MPI library by
5183 * Michael Fromberger but has been written from scratch with
5184 * additional optimizations in place.
5186 * The library is free for all purposes without any express
5187 * guarantee it works.
5189 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
5193 int mp_neg (mp_int * a, mp_int * b)
5197 if ((res = mp_copy (a, b)) != MP_OKAY) {
5202 if (mp_iszero(b) != MP_YES) {
5203 b->sign = (a->sign == MP_ZPOS) ? MP_NEG : MP_ZPOS;
5212 /* End: bn_mp_neg.c */
5214 /* Start: bn_mp_or.c */
5217 /* LibTomMath, multiple-precision integer library -- Tom St Denis
5219 * LibTomMath is a library that provides multiple-precision
5220 * integer arithmetic as well as number theoretic functionality.
5222 * The library was designed directly after the MPI library by
5223 * Michael Fromberger but has been written from scratch with
5224 * additional optimizations in place.
5226 * The library is free for all purposes without any express
5227 * guarantee it works.
5229 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
5232 /* OR two ints together */
5233 int mp_or (mp_int * a, mp_int * b, mp_int * c)
5238 if (a->used > b->used) {
5239 if ((res = mp_init_copy (&t, a)) != MP_OKAY) {
5245 if ((res = mp_init_copy (&t, b)) != MP_OKAY) {
5252 for (ix = 0; ix < px; ix++) {
5253 t.dp[ix] |= x->dp[ix];
5262 /* End: bn_mp_or.c */
5264 /* Start: bn_mp_prime_fermat.c */
5266 #ifdef BN_MP_PRIME_FERMAT_C
5267 /* LibTomMath, multiple-precision integer library -- Tom St Denis
5269 * LibTomMath is a library that provides multiple-precision
5270 * integer arithmetic as well as number theoretic functionality.
5272 * The library was designed directly after the MPI library by
5273 * Michael Fromberger but has been written from scratch with
5274 * additional optimizations in place.
5276 * The library is free for all purposes without any express
5277 * guarantee it works.
5279 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
5282 /* performs one Fermat test.
5284 * If "a" were prime then b**a == b (mod a) since the order of
5285 * the multiplicative sub-group would be phi(a) = a-1. That means
5286 * it would be the same as b**(a mod (a-1)) == b**1 == b (mod a).
5288 * Sets result to 1 if the congruence holds, or zero otherwise.
5290 int mp_prime_fermat (mp_int * a, mp_int * b, int *result)
5295 /* default to composite */
5299 if (mp_cmp_d(b, 1) != MP_GT) {
5304 if ((err = mp_init (&t)) != MP_OKAY) {
5308 /* compute t = b**a mod a */
5309 if ((err = mp_exptmod (b, a, a, &t)) != MP_OKAY) {
5313 /* is it equal to b? */
5314 if (mp_cmp (&t, b) == MP_EQ) {
5319 LBL_T:mp_clear (&t);
5324 /* End: bn_mp_prime_fermat.c */
5326 /* Start: bn_mp_prime_is_divisible.c */
5328 #ifdef BN_MP_PRIME_IS_DIVISIBLE_C
5329 /* LibTomMath, multiple-precision integer library -- Tom St Denis
5331 * LibTomMath is a library that provides multiple-precision
5332 * integer arithmetic as well as number theoretic functionality.
5334 * The library was designed directly after the MPI library by
5335 * Michael Fromberger but has been written from scratch with
5336 * additional optimizations in place.
5338 * The library is free for all purposes without any express
5339 * guarantee it works.
5341 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
5344 /* determines if an integers is divisible by one
5345 * of the first PRIME_SIZE primes or not
5347 * sets result to 0 if not, 1 if yes
5349 int mp_prime_is_divisible (mp_int * a, int *result)
5354 /* default to not */
5357 for (ix = 0; ix < PRIME_SIZE; ix++) {
5358 /* what is a mod LBL_prime_tab[ix] */
5359 if ((err = mp_mod_d (a, ltm_prime_tab[ix], &res)) != MP_OKAY) {
5363 /* is the residue zero? */
5374 /* End: bn_mp_prime_is_divisible.c */
5376 /* Start: bn_mp_prime_is_prime.c */
5378 #ifdef BN_MP_PRIME_IS_PRIME_C
5379 /* LibTomMath, multiple-precision integer library -- Tom St Denis
5381 * LibTomMath is a library that provides multiple-precision
5382 * integer arithmetic as well as number theoretic functionality.
5384 * The library was designed directly after the MPI library by
5385 * Michael Fromberger but has been written from scratch with
5386 * additional optimizations in place.
5388 * The library is free for all purposes without any express
5389 * guarantee it works.
5391 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
5394 /* performs a variable number of rounds of Miller-Rabin
5396 * Probability of error after t rounds is no more than
5399 * Sets result to 1 if probably prime, 0 otherwise
5401 int mp_prime_is_prime (mp_int * a, int t, int *result)
5409 /* valid value of t? */
5410 if (t <= 0 || t > PRIME_SIZE) {
5414 /* is the input equal to one of the primes in the table? */
5415 for (ix = 0; ix < PRIME_SIZE; ix++) {
5416 if (mp_cmp_d(a, ltm_prime_tab[ix]) == MP_EQ) {
5422 /* first perform trial division */
5423 if ((err = mp_prime_is_divisible (a, &res)) != MP_OKAY) {
5427 /* return if it was trivially divisible */
5428 if (res == MP_YES) {
5432 /* now perform the miller-rabin rounds */
5433 if ((err = mp_init (&b)) != MP_OKAY) {
5437 for (ix = 0; ix < t; ix++) {
5439 mp_set (&b, ltm_prime_tab[ix]);
5441 if ((err = mp_prime_miller_rabin (a, &b, &res)) != MP_OKAY) {
5450 /* passed the test */
5452 LBL_B:mp_clear (&b);
5457 /* End: bn_mp_prime_is_prime.c */
5459 /* Start: bn_mp_prime_miller_rabin.c */
5461 #ifdef BN_MP_PRIME_MILLER_RABIN_C
5462 /* LibTomMath, multiple-precision integer library -- Tom St Denis
5464 * LibTomMath is a library that provides multiple-precision
5465 * integer arithmetic as well as number theoretic functionality.
5467 * The library was designed directly after the MPI library by
5468 * Michael Fromberger but has been written from scratch with
5469 * additional optimizations in place.
5471 * The library is free for all purposes without any express
5472 * guarantee it works.
5474 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
5477 /* Miller-Rabin test of "a" to the base of "b" as described in
5478 * HAC pp. 139 Algorithm 4.24
5480 * Sets result to 0 if definitely composite or 1 if probably prime.
5481 * Randomly the chance of error is no more than 1/4 and often
5484 int mp_prime_miller_rabin (mp_int * a, mp_int * b, int *result)
5493 if (mp_cmp_d(b, 1) != MP_GT) {
5497 /* get n1 = a - 1 */
5498 if ((err = mp_init_copy (&n1, a)) != MP_OKAY) {
5501 if ((err = mp_sub_d (&n1, 1, &n1)) != MP_OKAY) {
5505 /* set 2**s * r = n1 */
5506 if ((err = mp_init_copy (&r, &n1)) != MP_OKAY) {
5510 /* count the number of least significant bits
5515 /* now divide n - 1 by 2**s */
5516 if ((err = mp_div_2d (&r, s, &r, NULL)) != MP_OKAY) {
5520 /* compute y = b**r mod a */
5521 if ((err = mp_init (&y)) != MP_OKAY) {
5524 if ((err = mp_exptmod (b, &r, a, &y)) != MP_OKAY) {
5528 /* if y != 1 and y != n1 do */
5529 if (mp_cmp_d (&y, 1) != MP_EQ && mp_cmp (&y, &n1) != MP_EQ) {
5531 /* while j <= s-1 and y != n1 */
5532 while ((j <= (s - 1)) && mp_cmp (&y, &n1) != MP_EQ) {
5533 if ((err = mp_sqrmod (&y, a, &y)) != MP_OKAY) {
5537 /* if y == 1 then composite */
5538 if (mp_cmp_d (&y, 1) == MP_EQ) {
5545 /* if y != n1 then composite */
5546 if (mp_cmp (&y, &n1) != MP_EQ) {
5551 /* probably prime now */
5553 LBL_Y:mp_clear (&y);
5554 LBL_R:mp_clear (&r);
5555 LBL_N1:mp_clear (&n1);
5560 /* End: bn_mp_prime_miller_rabin.c */
5562 /* Start: bn_mp_prime_next_prime.c */
5564 #ifdef BN_MP_PRIME_NEXT_PRIME_C
5565 /* LibTomMath, multiple-precision integer library -- Tom St Denis
5567 * LibTomMath is a library that provides multiple-precision
5568 * integer arithmetic as well as number theoretic functionality.
5570 * The library was designed directly after the MPI library by
5571 * Michael Fromberger but has been written from scratch with
5572 * additional optimizations in place.
5574 * The library is free for all purposes without any express
5575 * guarantee it works.
5577 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
5580 /* finds the next prime after the number "a" using "t" trials
5583 * bbs_style = 1 means the prime must be congruent to 3 mod 4
5585 int mp_prime_next_prime(mp_int *a, int t, int bbs_style)
5588 mp_digit res_tab[PRIME_SIZE], step, kstep;
5591 /* ensure t is valid */
5592 if (t <= 0 || t > PRIME_SIZE) {
5596 /* force positive */
5599 /* simple algo if a is less than the largest prime in the table */
5600 if (mp_cmp_d(a, ltm_prime_tab[PRIME_SIZE-1]) == MP_LT) {
5601 /* find which prime it is bigger than */
5602 for (x = PRIME_SIZE - 2; x >= 0; x--) {
5603 if (mp_cmp_d(a, ltm_prime_tab[x]) != MP_LT) {
5604 if (bbs_style == 1) {
5605 /* ok we found a prime smaller or
5606 * equal [so the next is larger]
5608 * however, the prime must be
5609 * congruent to 3 mod 4
5611 if ((ltm_prime_tab[x + 1] & 3) != 3) {
5612 /* scan upwards for a prime congruent to 3 mod 4 */
5613 for (y = x + 1; y < PRIME_SIZE; y++) {
5614 if ((ltm_prime_tab[y] & 3) == 3) {
5615 mp_set(a, ltm_prime_tab[y]);
5621 mp_set(a, ltm_prime_tab[x + 1]);
5626 /* at this point a maybe 1 */
5627 if (mp_cmp_d(a, 1) == MP_EQ) {
5631 /* fall through to the sieve */
5634 /* generate a prime congruent to 3 mod 4 or 1/3 mod 4? */
5635 if (bbs_style == 1) {
5641 /* at this point we will use a combination of a sieve and Miller-Rabin */
5643 if (bbs_style == 1) {
5644 /* if a mod 4 != 3 subtract the correct value to make it so */
5645 if ((a->dp[0] & 3) != 3) {
5646 if ((err = mp_sub_d(a, (a->dp[0] & 3) + 1, a)) != MP_OKAY) { return err; };
5649 if (mp_iseven(a) == 1) {
5651 if ((err = mp_sub_d(a, 1, a)) != MP_OKAY) {
5657 /* generate the restable */
5658 for (x = 1; x < PRIME_SIZE; x++) {
5659 if ((err = mp_mod_d(a, ltm_prime_tab[x], res_tab + x)) != MP_OKAY) {
5664 /* init temp used for Miller-Rabin Testing */
5665 if ((err = mp_init(&b)) != MP_OKAY) {
5670 /* skip to the next non-trivially divisible candidate */
5673 /* y == 1 if any residue was zero [e.g. cannot be prime] */
5676 /* increase step to next candidate */
5679 /* compute the new residue without using division */
5680 for (x = 1; x < PRIME_SIZE; x++) {
5681 /* add the step to each residue */
5682 res_tab[x] += kstep;
5684 /* subtract the modulus [instead of using division] */
5685 if (res_tab[x] >= ltm_prime_tab[x]) {
5686 res_tab[x] -= ltm_prime_tab[x];
5689 /* set flag if zero */
5690 if (res_tab[x] == 0) {
5694 } while (y == 1 && step < ((((mp_digit)1)<<DIGIT_BIT) - kstep));
5697 if ((err = mp_add_d(a, step, a)) != MP_OKAY) {
5701 /* if didn't pass sieve and step == MAX then skip test */
5702 if (y == 1 && step >= ((((mp_digit)1)<<DIGIT_BIT) - kstep)) {
5706 /* is this prime? */
5707 for (x = 0; x < t; x++) {
5708 mp_set(&b, ltm_prime_tab[t]);
5709 if ((err = mp_prime_miller_rabin(a, &b, &res)) != MP_OKAY) {
5717 if (res == MP_YES) {
5730 /* End: bn_mp_prime_next_prime.c */
5732 /* Start: bn_mp_prime_rabin_miller_trials.c */
5734 #ifdef BN_MP_PRIME_RABIN_MILLER_TRIALS_C
5735 /* LibTomMath, multiple-precision integer library -- Tom St Denis
5737 * LibTomMath is a library that provides multiple-precision
5738 * integer arithmetic as well as number theoretic functionality.
5740 * The library was designed directly after the MPI library by
5741 * Michael Fromberger but has been written from scratch with
5742 * additional optimizations in place.
5744 * The library is free for all purposes without any express
5745 * guarantee it works.
5747 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
5751 static const struct {
5764 /* returns # of RM trials required for a given bit size */
5765 int mp_prime_rabin_miller_trials(int size)
5769 for (x = 0; x < (int)(sizeof(sizes)/(sizeof(sizes[0]))); x++) {
5770 if (sizes[x].k == size) {
5772 } else if (sizes[x].k > size) {
5773 return (x == 0) ? sizes[0].t : sizes[x - 1].t;
5776 return sizes[x-1].t + 1;
5782 /* End: bn_mp_prime_rabin_miller_trials.c */
5784 /* Start: bn_mp_prime_random_ex.c */
5786 #ifdef BN_MP_PRIME_RANDOM_EX_C
5787 /* LibTomMath, multiple-precision integer library -- Tom St Denis
5789 * LibTomMath is a library that provides multiple-precision
5790 * integer arithmetic as well as number theoretic functionality.
5792 * The library was designed directly after the MPI library by
5793 * Michael Fromberger but has been written from scratch with
5794 * additional optimizations in place.
5796 * The library is free for all purposes without any express
5797 * guarantee it works.
5799 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
5802 /* makes a truly random prime of a given size (bits),
5804 * Flags are as follows:
5806 * LTM_PRIME_BBS - make prime congruent to 3 mod 4
5807 * LTM_PRIME_SAFE - make sure (p-1)/2 is prime as well (implies LTM_PRIME_BBS)
5808 * LTM_PRIME_2MSB_OFF - make the 2nd highest bit zero
5809 * LTM_PRIME_2MSB_ON - make the 2nd highest bit one
5811 * You have to supply a callback which fills in a buffer with random bytes. "dat" is a parameter you can
5812 * have passed to the callback (e.g. a state or something). This function doesn't use "dat" itself
5817 /* This is possibly the mother of all prime generation functions, muahahahahaha! */
5818 int mp_prime_random_ex(mp_int *a, int t, int size, int flags, ltm_prime_callback cb, void *dat)
5820 unsigned char *tmp, maskAND, maskOR_msb, maskOR_lsb;
5821 int res, err, bsize, maskOR_msb_offset;
5823 /* sanity check the input */
5824 if (size <= 1 || t <= 0) {
5828 /* LTM_PRIME_SAFE implies LTM_PRIME_BBS */
5829 if (flags & LTM_PRIME_SAFE) {
5830 flags |= LTM_PRIME_BBS;
5833 /* calc the byte size */
5834 bsize = (size>>3) + ((size&7)?1:0);
5836 /* we need a buffer of bsize bytes */
5837 tmp = OPT_CAST(unsigned char) XMALLOC(bsize);
5842 /* calc the maskAND value for the MSbyte*/
5843 maskAND = ((size&7) == 0) ? 0xFF : (0xFF >> (8 - (size & 7)));
5845 /* calc the maskOR_msb */
5847 maskOR_msb_offset = ((size & 7) == 1) ? 1 : 0;
5848 if (flags & LTM_PRIME_2MSB_ON) {
5849 maskOR_msb |= 1 << ((size - 2) & 7);
5850 } else if (flags & LTM_PRIME_2MSB_OFF) {
5851 maskAND &= ~(1 << ((size - 2) & 7));
5854 /* get the maskOR_lsb */
5856 if (flags & LTM_PRIME_BBS) {
5861 /* read the bytes */
5862 if (cb(tmp, bsize, dat) != bsize) {
5867 /* work over the MSbyte */
5869 tmp[0] |= 1 << ((size - 1) & 7);
5871 /* mix in the maskORs */
5872 tmp[maskOR_msb_offset] |= maskOR_msb;
5873 tmp[bsize-1] |= maskOR_lsb;
5876 if ((err = mp_read_unsigned_bin(a, tmp, bsize)) != MP_OKAY) { goto error; }
5879 if ((err = mp_prime_is_prime(a, t, &res)) != MP_OKAY) { goto error; }
5884 if (flags & LTM_PRIME_SAFE) {
5885 /* see if (a-1)/2 is prime */
5886 if ((err = mp_sub_d(a, 1, a)) != MP_OKAY) { goto error; }
5887 if ((err = mp_div_2(a, a)) != MP_OKAY) { goto error; }
5890 if ((err = mp_prime_is_prime(a, t, &res)) != MP_OKAY) { goto error; }
5892 } while (res == MP_NO);
5894 if (flags & LTM_PRIME_SAFE) {
5895 /* restore a to the original value */
5896 if ((err = mp_mul_2(a, a)) != MP_OKAY) { goto error; }
5897 if ((err = mp_add_d(a, 1, a)) != MP_OKAY) { goto error; }
5909 /* End: bn_mp_prime_random_ex.c */
5911 /* Start: bn_mp_radix_size.c */
5913 #ifdef BN_MP_RADIX_SIZE_C
5914 /* LibTomMath, multiple-precision integer library -- Tom St Denis
5916 * LibTomMath is a library that provides multiple-precision
5917 * integer arithmetic as well as number theoretic functionality.
5919 * The library was designed directly after the MPI library by
5920 * Michael Fromberger but has been written from scratch with
5921 * additional optimizations in place.
5923 * The library is free for all purposes without any express
5924 * guarantee it works.
5926 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
5929 /* returns size of ASCII reprensentation */
5930 int mp_radix_size (mp_int * a, int radix, int *size)
5938 /* special case for binary */
5940 *size = mp_count_bits (a) + (a->sign == MP_NEG ? 1 : 0) + 1;
5944 /* make sure the radix is in range */
5945 if (radix < 2 || radix > 64) {
5949 if (mp_iszero(a) == MP_YES) {
5954 /* digs is the digit count */
5957 /* if it's negative add one for the sign */
5958 if (a->sign == MP_NEG) {
5962 /* init a copy of the input */
5963 if ((res = mp_init_copy (&t, a)) != MP_OKAY) {
5967 /* force temp to positive */
5970 /* fetch out all of the digits */
5971 while (mp_iszero (&t) == MP_NO) {
5972 if ((res = mp_div_d (&t, (mp_digit) radix, &t, &d)) != MP_OKAY) {
5980 /* return digs + 1, the 1 is for the NULL byte that would be required. */
5987 /* End: bn_mp_radix_size.c */
5989 /* Start: bn_mp_radix_smap.c */
5991 #ifdef BN_MP_RADIX_SMAP_C
5992 /* LibTomMath, multiple-precision integer library -- Tom St Denis
5994 * LibTomMath is a library that provides multiple-precision
5995 * integer arithmetic as well as number theoretic functionality.
5997 * The library was designed directly after the MPI library by
5998 * Michael Fromberger but has been written from scratch with
5999 * additional optimizations in place.
6001 * The library is free for all purposes without any express
6002 * guarantee it works.
6004 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
6007 /* chars used in radix conversions */
6008 const char *mp_s_rmap = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz+/";
6011 /* End: bn_mp_radix_smap.c */
6013 /* Start: bn_mp_rand.c */
6016 /* LibTomMath, multiple-precision integer library -- Tom St Denis
6018 * LibTomMath is a library that provides multiple-precision
6019 * integer arithmetic as well as number theoretic functionality.
6021 * The library was designed directly after the MPI library by
6022 * Michael Fromberger but has been written from scratch with
6023 * additional optimizations in place.
6025 * The library is free for all purposes without any express
6026 * guarantee it works.
6028 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
6031 /* makes a pseudo-random int of a given size */
6033 mp_rand (mp_int * a, int digits)
6043 /* first place a random non-zero digit */
6045 d = ((mp_digit) abs (rand ())) & MP_MASK;
6048 if ((res = mp_add_d (a, d, a)) != MP_OKAY) {
6052 while (--digits > 0) {
6053 if ((res = mp_lshd (a, 1)) != MP_OKAY) {
6057 if ((res = mp_add_d (a, ((mp_digit) abs (rand ())), a)) != MP_OKAY) {
6066 /* End: bn_mp_rand.c */
6068 /* Start: bn_mp_read_radix.c */
6070 #ifdef BN_MP_READ_RADIX_C
6071 /* LibTomMath, multiple-precision integer library -- Tom St Denis
6073 * LibTomMath is a library that provides multiple-precision
6074 * integer arithmetic as well as number theoretic functionality.
6076 * The library was designed directly after the MPI library by
6077 * Michael Fromberger but has been written from scratch with
6078 * additional optimizations in place.
6080 * The library is free for all purposes without any express
6081 * guarantee it works.
6083 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
6086 /* read a string [ASCII] in a given radix */
6087 int mp_read_radix (mp_int * a, const char *str, int radix)
6092 /* make sure the radix is ok */
6093 if (radix < 2 || radix > 64) {
6097 /* if the leading digit is a
6098 * minus set the sign to negative.
6107 /* set the integer to the default of zero */
6110 /* process each digit of the string */
6112 /* if the radix < 36 the conversion is case insensitive
6113 * this allows numbers like 1AB and 1ab to represent the same value
6116 ch = (char) ((radix < 36) ? toupper (*str) : *str);
6117 for (y = 0; y < 64; y++) {
6118 if (ch == mp_s_rmap[y]) {
6123 /* if the char was found in the map
6124 * and is less than the given radix add it
6125 * to the number, otherwise exit the loop.
6128 if ((res = mp_mul_d (a, (mp_digit) radix, a)) != MP_OKAY) {
6131 if ((res = mp_add_d (a, (mp_digit) y, a)) != MP_OKAY) {
6140 /* set the sign only if a != 0 */
6141 if (mp_iszero(a) != 1) {
6148 /* End: bn_mp_read_radix.c */
6150 /* Start: bn_mp_read_signed_bin.c */
6152 #ifdef BN_MP_READ_SIGNED_BIN_C
6153 /* LibTomMath, multiple-precision integer library -- Tom St Denis
6155 * LibTomMath is a library that provides multiple-precision
6156 * integer arithmetic as well as number theoretic functionality.
6158 * The library was designed directly after the MPI library by
6159 * Michael Fromberger but has been written from scratch with
6160 * additional optimizations in place.
6162 * The library is free for all purposes without any express
6163 * guarantee it works.
6165 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
6168 /* read signed bin, big endian, first byte is 0==positive or 1==negative */
6170 mp_read_signed_bin (mp_int * a, unsigned char *b, int c)
6174 /* read magnitude */
6175 if ((res = mp_read_unsigned_bin (a, b + 1, c - 1)) != MP_OKAY) {
6179 /* first byte is 0 for positive, non-zero for negative */
6190 /* End: bn_mp_read_signed_bin.c */
6192 /* Start: bn_mp_read_unsigned_bin.c */
6194 #ifdef BN_MP_READ_UNSIGNED_BIN_C
6195 /* LibTomMath, multiple-precision integer library -- Tom St Denis
6197 * LibTomMath is a library that provides multiple-precision
6198 * integer arithmetic as well as number theoretic functionality.
6200 * The library was designed directly after the MPI library by
6201 * Michael Fromberger but has been written from scratch with
6202 * additional optimizations in place.
6204 * The library is free for all purposes without any express
6205 * guarantee it works.
6207 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
6210 /* reads a unsigned char array, assumes the msb is stored first [big endian] */
6212 mp_read_unsigned_bin (mp_int * a, unsigned char *b, int c)
6216 /* make sure there are at least two digits */
6218 if ((res = mp_grow(a, 2)) != MP_OKAY) {
6226 /* read the bytes in */
6228 if ((res = mp_mul_2d (a, 8, a)) != MP_OKAY) {
6236 a->dp[0] = (*b & MP_MASK);
6237 a->dp[1] |= ((*b++ >> 7U) & 1);
6246 /* End: bn_mp_read_unsigned_bin.c */
6248 /* Start: bn_mp_reduce.c */
6250 #ifdef BN_MP_REDUCE_C
6251 /* LibTomMath, multiple-precision integer library -- Tom St Denis
6253 * LibTomMath is a library that provides multiple-precision
6254 * integer arithmetic as well as number theoretic functionality.
6256 * The library was designed directly after the MPI library by
6257 * Michael Fromberger but has been written from scratch with
6258 * additional optimizations in place.
6260 * The library is free for all purposes without any express
6261 * guarantee it works.
6263 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
6266 /* reduces x mod m, assumes 0 < x < m**2, mu is
6267 * precomputed via mp_reduce_setup.
6268 * From HAC pp.604 Algorithm 14.42
6270 int mp_reduce (mp_int * x, mp_int * m, mp_int * mu)
6273 int res, um = m->used;
6276 if ((res = mp_init_copy (&q, x)) != MP_OKAY) {
6280 /* q1 = x / b**(k-1) */
6281 mp_rshd (&q, um - 1);
6283 /* according to HAC this optimization is ok */
6284 if (((unsigned long) um) > (((mp_digit)1) << (DIGIT_BIT - 1))) {
6285 if ((res = mp_mul (&q, mu, &q)) != MP_OKAY) {
6289 #ifdef BN_S_MP_MUL_HIGH_DIGS_C
6290 if ((res = s_mp_mul_high_digs (&q, mu, &q, um)) != MP_OKAY) {
6293 #elif defined(BN_FAST_S_MP_MUL_HIGH_DIGS_C)
6294 if ((res = fast_s_mp_mul_high_digs (&q, mu, &q, um)) != MP_OKAY) {
6305 /* q3 = q2 / b**(k+1) */
6306 mp_rshd (&q, um + 1);
6308 /* x = x mod b**(k+1), quick (no division) */
6309 if ((res = mp_mod_2d (x, DIGIT_BIT * (um + 1), x)) != MP_OKAY) {
6313 /* q = q * m mod b**(k+1), quick (no division) */
6314 if ((res = s_mp_mul_digs (&q, m, &q, um + 1)) != MP_OKAY) {
6319 if ((res = mp_sub (x, &q, x)) != MP_OKAY) {
6323 /* If x < 0, add b**(k+1) to it */
6324 if (mp_cmp_d (x, 0) == MP_LT) {
6326 if ((res = mp_lshd (&q, um + 1)) != MP_OKAY)
6328 if ((res = mp_add (x, &q, x)) != MP_OKAY)
6332 /* Back off if it's too big */
6333 while (mp_cmp (x, m) != MP_LT) {
6334 if ((res = s_mp_sub (x, m, x)) != MP_OKAY) {
6346 /* End: bn_mp_reduce.c */
6348 /* Start: bn_mp_reduce_2k.c */
6350 #ifdef BN_MP_REDUCE_2K_C
6351 /* LibTomMath, multiple-precision integer library -- Tom St Denis
6353 * LibTomMath is a library that provides multiple-precision
6354 * integer arithmetic as well as number theoretic functionality.
6356 * The library was designed directly after the MPI library by
6357 * Michael Fromberger but has been written from scratch with
6358 * additional optimizations in place.
6360 * The library is free for all purposes without any express
6361 * guarantee it works.
6363 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
6366 /* reduces a modulo n where n is of the form 2**p - d */
6367 int mp_reduce_2k(mp_int *a, mp_int *n, mp_digit d)
6372 if ((res = mp_init(&q)) != MP_OKAY) {
6376 p = mp_count_bits(n);
6378 /* q = a/2**p, a = a mod 2**p */
6379 if ((res = mp_div_2d(a, p, &q, a)) != MP_OKAY) {
6385 if ((res = mp_mul_d(&q, d, &q)) != MP_OKAY) {
6391 if ((res = s_mp_add(a, &q, a)) != MP_OKAY) {
6395 if (mp_cmp_mag(a, n) != MP_LT) {
6407 /* End: bn_mp_reduce_2k.c */
6409 /* Start: bn_mp_reduce_2k_l.c */
6411 #ifdef BN_MP_REDUCE_2K_L_C
6412 /* LibTomMath, multiple-precision integer library -- Tom St Denis
6414 * LibTomMath is a library that provides multiple-precision
6415 * integer arithmetic as well as number theoretic functionality.
6417 * The library was designed directly after the MPI library by
6418 * Michael Fromberger but has been written from scratch with
6419 * additional optimizations in place.
6421 * The library is free for all purposes without any express
6422 * guarantee it works.
6424 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
6427 /* reduces a modulo n where n is of the form 2**p - d
6428 This differs from reduce_2k since "d" can be larger
6429 than a single digit.
6431 int mp_reduce_2k_l(mp_int *a, mp_int *n, mp_int *d)
6436 if ((res = mp_init(&q)) != MP_OKAY) {
6440 p = mp_count_bits(n);
6442 /* q = a/2**p, a = a mod 2**p */
6443 if ((res = mp_div_2d(a, p, &q, a)) != MP_OKAY) {
6448 if ((res = mp_mul(&q, d, &q)) != MP_OKAY) {
6453 if ((res = s_mp_add(a, &q, a)) != MP_OKAY) {
6457 if (mp_cmp_mag(a, n) != MP_LT) {
6469 /* End: bn_mp_reduce_2k_l.c */
6471 /* Start: bn_mp_reduce_2k_setup.c */
6473 #ifdef BN_MP_REDUCE_2K_SETUP_C
6474 /* LibTomMath, multiple-precision integer library -- Tom St Denis
6476 * LibTomMath is a library that provides multiple-precision
6477 * integer arithmetic as well as number theoretic functionality.
6479 * The library was designed directly after the MPI library by
6480 * Michael Fromberger but has been written from scratch with
6481 * additional optimizations in place.
6483 * The library is free for all purposes without any express
6484 * guarantee it works.
6486 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
6489 /* determines the setup value */
6490 int mp_reduce_2k_setup(mp_int *a, mp_digit *d)
6495 if ((res = mp_init(&tmp)) != MP_OKAY) {
6499 p = mp_count_bits(a);
6500 if ((res = mp_2expt(&tmp, p)) != MP_OKAY) {
6505 if ((res = s_mp_sub(&tmp, a, &tmp)) != MP_OKAY) {
6516 /* End: bn_mp_reduce_2k_setup.c */
6518 /* Start: bn_mp_reduce_2k_setup_l.c */
6520 #ifdef BN_MP_REDUCE_2K_SETUP_L_C
6521 /* LibTomMath, multiple-precision integer library -- Tom St Denis
6523 * LibTomMath is a library that provides multiple-precision
6524 * integer arithmetic as well as number theoretic functionality.
6526 * The library was designed directly after the MPI library by
6527 * Michael Fromberger but has been written from scratch with
6528 * additional optimizations in place.
6530 * The library is free for all purposes without any express
6531 * guarantee it works.
6533 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
6536 /* determines the setup value */
6537 int mp_reduce_2k_setup_l(mp_int *a, mp_int *d)
6542 if ((res = mp_init(&tmp)) != MP_OKAY) {
6546 if ((res = mp_2expt(&tmp, mp_count_bits(a))) != MP_OKAY) {
6550 if ((res = s_mp_sub(&tmp, a, d)) != MP_OKAY) {
6560 /* End: bn_mp_reduce_2k_setup_l.c */
6562 /* Start: bn_mp_reduce_is_2k.c */
6564 #ifdef BN_MP_REDUCE_IS_2K_C
6565 /* LibTomMath, multiple-precision integer library -- Tom St Denis
6567 * LibTomMath is a library that provides multiple-precision
6568 * integer arithmetic as well as number theoretic functionality.
6570 * The library was designed directly after the MPI library by
6571 * Michael Fromberger but has been written from scratch with
6572 * additional optimizations in place.
6574 * The library is free for all purposes without any express
6575 * guarantee it works.
6577 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
6580 /* determines if mp_reduce_2k can be used */
6581 int mp_reduce_is_2k(mp_int *a)
6588 } else if (a->used == 1) {
6590 } else if (a->used > 1) {
6591 iy = mp_count_bits(a);
6595 /* Test every bit from the second digit up, must be 1 */
6596 for (ix = DIGIT_BIT; ix < iy; ix++) {
6597 if ((a->dp[iw] & iz) == 0) {
6601 if (iz > (mp_digit)MP_MASK) {
6612 /* End: bn_mp_reduce_is_2k.c */
6614 /* Start: bn_mp_reduce_is_2k_l.c */
6616 #ifdef BN_MP_REDUCE_IS_2K_L_C
6617 /* LibTomMath, multiple-precision integer library -- Tom St Denis
6619 * LibTomMath is a library that provides multiple-precision
6620 * integer arithmetic as well as number theoretic functionality.
6622 * The library was designed directly after the MPI library by
6623 * Michael Fromberger but has been written from scratch with
6624 * additional optimizations in place.
6626 * The library is free for all purposes without any express
6627 * guarantee it works.
6629 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
6632 /* determines if reduce_2k_l can be used */
6633 int mp_reduce_is_2k_l(mp_int *a)
6639 } else if (a->used == 1) {
6641 } else if (a->used > 1) {
6642 /* if more than half of the digits are -1 we're sold */
6643 for (iy = ix = 0; ix < a->used; ix++) {
6644 if (a->dp[ix] == MP_MASK) {
6648 return (iy >= (a->used/2)) ? MP_YES : MP_NO;
6656 /* End: bn_mp_reduce_is_2k_l.c */
6658 /* Start: bn_mp_reduce_setup.c */
6660 #ifdef BN_MP_REDUCE_SETUP_C
6661 /* LibTomMath, multiple-precision integer library -- Tom St Denis
6663 * LibTomMath is a library that provides multiple-precision
6664 * integer arithmetic as well as number theoretic functionality.
6666 * The library was designed directly after the MPI library by
6667 * Michael Fromberger but has been written from scratch with
6668 * additional optimizations in place.
6670 * The library is free for all purposes without any express
6671 * guarantee it works.
6673 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
6676 /* pre-calculate the value required for Barrett reduction
6677 * For a given modulus "b" it calulates the value required in "a"
6679 int mp_reduce_setup (mp_int * a, mp_int * b)
6683 if ((res = mp_2expt (a, b->used * 2 * DIGIT_BIT)) != MP_OKAY) {
6686 return mp_div (a, b, a, NULL);
6690 /* End: bn_mp_reduce_setup.c */
6692 /* Start: bn_mp_rshd.c */
6695 /* LibTomMath, multiple-precision integer library -- Tom St Denis
6697 * LibTomMath is a library that provides multiple-precision
6698 * integer arithmetic as well as number theoretic functionality.
6700 * The library was designed directly after the MPI library by
6701 * Michael Fromberger but has been written from scratch with
6702 * additional optimizations in place.
6704 * The library is free for all purposes without any express
6705 * guarantee it works.
6707 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
6710 /* shift right a certain amount of digits */
6711 void mp_rshd (mp_int * a, int b)
6715 /* if b <= 0 then ignore it */
6720 /* if b > used then simply zero it and return */
6727 register mp_digit *bottom, *top;
6729 /* shift the digits down */
6734 /* top [offset into digits] */
6737 /* this is implemented as a sliding window where
6738 * the window is b-digits long and digits from
6739 * the top of the window are copied to the bottom
6743 b-2 | b-1 | b0 | b1 | b2 | ... | bb | ---->
6745 \-------------------/ ---->
6747 for (x = 0; x < (a->used - b); x++) {
6751 /* zero the top digits */
6752 for (; x < a->used; x++) {
6757 /* remove excess digits */
6762 /* End: bn_mp_rshd.c */
6764 /* Start: bn_mp_set.c */
6767 /* LibTomMath, multiple-precision integer library -- Tom St Denis
6769 * LibTomMath is a library that provides multiple-precision
6770 * integer arithmetic as well as number theoretic functionality.
6772 * The library was designed directly after the MPI library by
6773 * Michael Fromberger but has been written from scratch with
6774 * additional optimizations in place.
6776 * The library is free for all purposes without any express
6777 * guarantee it works.
6779 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
6782 /* set to a digit */
6783 void mp_set (mp_int * a, mp_digit b)
6786 a->dp[0] = b & MP_MASK;
6787 a->used = (a->dp[0] != 0) ? 1 : 0;
6791 /* End: bn_mp_set.c */
6793 /* Start: bn_mp_set_int.c */
6795 #ifdef BN_MP_SET_INT_C
6796 /* LibTomMath, multiple-precision integer library -- Tom St Denis
6798 * LibTomMath is a library that provides multiple-precision
6799 * integer arithmetic as well as number theoretic functionality.
6801 * The library was designed directly after the MPI library by
6802 * Michael Fromberger but has been written from scratch with
6803 * additional optimizations in place.
6805 * The library is free for all purposes without any express
6806 * guarantee it works.
6808 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
6811 /* set a 32-bit const */
6812 int mp_set_int (mp_int * a, unsigned long b)
6818 /* set four bits at a time */
6819 for (x = 0; x < 8; x++) {
6820 /* shift the number up four bits */
6821 if ((res = mp_mul_2d (a, 4, a)) != MP_OKAY) {
6825 /* OR in the top four bits of the source */
6826 a->dp[0] |= (b >> 28) & 15;
6828 /* shift the source up to the next four bits */
6831 /* ensure that digits are not clamped off */
6839 /* End: bn_mp_set_int.c */
6841 /* Start: bn_mp_shrink.c */
6843 #ifdef BN_MP_SHRINK_C
6844 /* LibTomMath, multiple-precision integer library -- Tom St Denis
6846 * LibTomMath is a library that provides multiple-precision
6847 * integer arithmetic as well as number theoretic functionality.
6849 * The library was designed directly after the MPI library by
6850 * Michael Fromberger but has been written from scratch with
6851 * additional optimizations in place.
6853 * The library is free for all purposes without any express
6854 * guarantee it works.
6856 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
6859 /* shrink a bignum */
6860 int mp_shrink (mp_int * a)
6863 if (a->alloc != a->used && a->used > 0) {
6864 if ((tmp = OPT_CAST(mp_digit) XREALLOC (a->dp, sizeof (mp_digit) * a->used)) == NULL) {
6874 /* End: bn_mp_shrink.c */
6876 /* Start: bn_mp_signed_bin_size.c */
6878 #ifdef BN_MP_SIGNED_BIN_SIZE_C
6879 /* LibTomMath, multiple-precision integer library -- Tom St Denis
6881 * LibTomMath is a library that provides multiple-precision
6882 * integer arithmetic as well as number theoretic functionality.
6884 * The library was designed directly after the MPI library by
6885 * Michael Fromberger but has been written from scratch with
6886 * additional optimizations in place.
6888 * The library is free for all purposes without any express
6889 * guarantee it works.
6891 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
6894 /* get the size for an signed equivalent */
6895 int mp_signed_bin_size (mp_int * a)
6897 return 1 + mp_unsigned_bin_size (a);
6901 /* End: bn_mp_signed_bin_size.c */
6903 /* Start: bn_mp_sqr.c */
6906 /* LibTomMath, multiple-precision integer library -- Tom St Denis
6908 * LibTomMath is a library that provides multiple-precision
6909 * integer arithmetic as well as number theoretic functionality.
6911 * The library was designed directly after the MPI library by
6912 * Michael Fromberger but has been written from scratch with
6913 * additional optimizations in place.
6915 * The library is free for all purposes without any express
6916 * guarantee it works.
6918 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
6921 /* computes b = a*a */
6923 mp_sqr (mp_int * a, mp_int * b)
6927 #ifdef BN_MP_TOOM_SQR_C
6928 /* use Toom-Cook? */
6929 if (a->used >= TOOM_SQR_CUTOFF) {
6930 res = mp_toom_sqr(a, b);
6934 #ifdef BN_MP_KARATSUBA_SQR_C
6935 if (a->used >= KARATSUBA_SQR_CUTOFF) {
6936 res = mp_karatsuba_sqr (a, b);
6940 #ifdef BN_FAST_S_MP_SQR_C
6941 /* can we use the fast comba multiplier? */
6942 if ((a->used * 2 + 1) < MP_WARRAY &&
6944 (1 << (sizeof(mp_word) * CHAR_BIT - 2*DIGIT_BIT - 1))) {
6945 res = fast_s_mp_sqr (a, b);
6948 #ifdef BN_S_MP_SQR_C
6949 res = s_mp_sqr (a, b);
6959 /* End: bn_mp_sqr.c */
6961 /* Start: bn_mp_sqrmod.c */
6963 #ifdef BN_MP_SQRMOD_C
6964 /* LibTomMath, multiple-precision integer library -- Tom St Denis
6966 * LibTomMath is a library that provides multiple-precision
6967 * integer arithmetic as well as number theoretic functionality.
6969 * The library was designed directly after the MPI library by
6970 * Michael Fromberger but has been written from scratch with
6971 * additional optimizations in place.
6973 * The library is free for all purposes without any express
6974 * guarantee it works.
6976 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
6979 /* c = a * a (mod b) */
6981 mp_sqrmod (mp_int * a, mp_int * b, mp_int * c)
6986 if ((res = mp_init (&t)) != MP_OKAY) {
6990 if ((res = mp_sqr (a, &t)) != MP_OKAY) {
6994 res = mp_mod (&t, b, c);
7000 /* End: bn_mp_sqrmod.c */
7002 /* Start: bn_mp_sqrt.c */
7005 /* LibTomMath, multiple-precision integer library -- Tom St Denis
7007 * LibTomMath is a library that provides multiple-precision
7008 * integer arithmetic as well as number theoretic functionality.
7010 * The library was designed directly after the MPI library by
7011 * Michael Fromberger but has been written from scratch with
7012 * additional optimizations in place.
7014 * The library is free for all purposes without any express
7015 * guarantee it works.
7017 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
7020 /* this function is less generic than mp_n_root, simpler and faster */
7021 int mp_sqrt(mp_int *arg, mp_int *ret)
7026 /* must be positive */
7027 if (arg->sign == MP_NEG) {
7032 if (mp_iszero(arg) == MP_YES) {
7037 if ((res = mp_init_copy(&t1, arg)) != MP_OKAY) {
7041 if ((res = mp_init(&t2)) != MP_OKAY) {
7045 /* First approx. (not very bad for large arg) */
7046 mp_rshd (&t1,t1.used/2);
7049 if ((res = mp_div(arg,&t1,&t2,NULL)) != MP_OKAY) {
7052 if ((res = mp_add(&t1,&t2,&t1)) != MP_OKAY) {
7055 if ((res = mp_div_2(&t1,&t1)) != MP_OKAY) {
7058 /* And now t1 > sqrt(arg) */
7060 if ((res = mp_div(arg,&t1,&t2,NULL)) != MP_OKAY) {
7063 if ((res = mp_add(&t1,&t2,&t1)) != MP_OKAY) {
7066 if ((res = mp_div_2(&t1,&t1)) != MP_OKAY) {
7069 /* t1 >= sqrt(arg) >= t2 at this point */
7070 } while (mp_cmp_mag(&t1,&t2) == MP_GT);
7081 /* End: bn_mp_sqrt.c */
7083 /* Start: bn_mp_sub.c */
7086 /* LibTomMath, multiple-precision integer library -- Tom St Denis
7088 * LibTomMath is a library that provides multiple-precision
7089 * integer arithmetic as well as number theoretic functionality.
7091 * The library was designed directly after the MPI library by
7092 * Michael Fromberger but has been written from scratch with
7093 * additional optimizations in place.
7095 * The library is free for all purposes without any express
7096 * guarantee it works.
7098 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
7101 /* high level subtraction (handles signs) */
7103 mp_sub (mp_int * a, mp_int * b, mp_int * c)
7111 /* subtract a negative from a positive, OR */
7112 /* subtract a positive from a negative. */
7113 /* In either case, ADD their magnitudes, */
7114 /* and use the sign of the first number. */
7116 res = s_mp_add (a, b, c);
7118 /* subtract a positive from a positive, OR */
7119 /* subtract a negative from a negative. */
7120 /* First, take the difference between their */
7121 /* magnitudes, then... */
7122 if (mp_cmp_mag (a, b) != MP_LT) {
7123 /* Copy the sign from the first */
7125 /* The first has a larger or equal magnitude */
7126 res = s_mp_sub (a, b, c);
7128 /* The result has the *opposite* sign from */
7129 /* the first number. */
7130 c->sign = (sa == MP_ZPOS) ? MP_NEG : MP_ZPOS;
7131 /* The second has a larger magnitude */
7132 res = s_mp_sub (b, a, c);
7140 /* End: bn_mp_sub.c */
7142 /* Start: bn_mp_sub_d.c */
7144 #ifdef BN_MP_SUB_D_C
7145 /* LibTomMath, multiple-precision integer library -- Tom St Denis
7147 * LibTomMath is a library that provides multiple-precision
7148 * integer arithmetic as well as number theoretic functionality.
7150 * The library was designed directly after the MPI library by
7151 * Michael Fromberger but has been written from scratch with
7152 * additional optimizations in place.
7154 * The library is free for all purposes without any express
7155 * guarantee it works.
7157 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
7160 /* single digit subtraction */
7162 mp_sub_d (mp_int * a, mp_digit b, mp_int * c)
7164 mp_digit *tmpa, *tmpc, mu;
7165 int res, ix, oldused;
7167 /* grow c as required */
7168 if (c->alloc < a->used + 1) {
7169 if ((res = mp_grow(c, a->used + 1)) != MP_OKAY) {
7174 /* if a is negative just do an unsigned
7175 * addition [with fudged signs]
7177 if (a->sign == MP_NEG) {
7179 res = mp_add_d(a, b, c);
7180 a->sign = c->sign = MP_NEG;
7189 /* if a <= b simply fix the single digit */
7190 if ((a->used == 1 && a->dp[0] <= b) || a->used == 0) {
7192 *tmpc++ = b - *tmpa;
7198 /* negative/1digit */
7206 /* subtract first digit */
7207 *tmpc = *tmpa++ - b;
7208 mu = *tmpc >> (sizeof(mp_digit) * CHAR_BIT - 1);
7211 /* handle rest of the digits */
7212 for (ix = 1; ix < a->used; ix++) {
7213 *tmpc = *tmpa++ - mu;
7214 mu = *tmpc >> (sizeof(mp_digit) * CHAR_BIT - 1);
7219 /* zero excess digits */
7220 while (ix++ < oldused) {
7229 /* End: bn_mp_sub_d.c */
7231 /* Start: bn_mp_submod.c */
7233 #ifdef BN_MP_SUBMOD_C
7234 /* LibTomMath, multiple-precision integer library -- Tom St Denis
7236 * LibTomMath is a library that provides multiple-precision
7237 * integer arithmetic as well as number theoretic functionality.
7239 * The library was designed directly after the MPI library by
7240 * Michael Fromberger but has been written from scratch with
7241 * additional optimizations in place.
7243 * The library is free for all purposes without any express
7244 * guarantee it works.
7246 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
7249 /* d = a - b (mod c) */
7251 mp_submod (mp_int * a, mp_int * b, mp_int * c, mp_int * d)
7257 if ((res = mp_init (&t)) != MP_OKAY) {
7261 if ((res = mp_sub (a, b, &t)) != MP_OKAY) {
7265 res = mp_mod (&t, c, d);
7271 /* End: bn_mp_submod.c */
7273 /* Start: bn_mp_to_signed_bin.c */
7275 #ifdef BN_MP_TO_SIGNED_BIN_C
7276 /* LibTomMath, multiple-precision integer library -- Tom St Denis
7278 * LibTomMath is a library that provides multiple-precision
7279 * integer arithmetic as well as number theoretic functionality.
7281 * The library was designed directly after the MPI library by
7282 * Michael Fromberger but has been written from scratch with
7283 * additional optimizations in place.
7285 * The library is free for all purposes without any express
7286 * guarantee it works.
7288 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
7291 /* store in signed [big endian] format */
7292 int mp_to_signed_bin (mp_int * a, unsigned char *b)
7296 if ((res = mp_to_unsigned_bin (a, b + 1)) != MP_OKAY) {
7299 b[0] = (unsigned char) ((a->sign == MP_ZPOS) ? 0 : 1);
7304 /* End: bn_mp_to_signed_bin.c */
7306 /* Start: bn_mp_to_signed_bin_n.c */
7308 #ifdef BN_MP_TO_SIGNED_BIN_N_C
7309 /* LibTomMath, multiple-precision integer library -- Tom St Denis
7311 * LibTomMath is a library that provides multiple-precision
7312 * integer arithmetic as well as number theoretic functionality.
7314 * The library was designed directly after the MPI library by
7315 * Michael Fromberger but has been written from scratch with
7316 * additional optimizations in place.
7318 * The library is free for all purposes without any express
7319 * guarantee it works.
7321 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
7324 /* store in signed [big endian] format */
7325 int mp_to_signed_bin_n (mp_int * a, unsigned char *b, unsigned long *outlen)
7327 if (*outlen < (unsigned long)mp_signed_bin_size(a)) {
7330 *outlen = mp_signed_bin_size(a);
7331 return mp_to_signed_bin(a, b);
7335 /* End: bn_mp_to_signed_bin_n.c */
7337 /* Start: bn_mp_to_unsigned_bin.c */
7339 #ifdef BN_MP_TO_UNSIGNED_BIN_C
7340 /* LibTomMath, multiple-precision integer library -- Tom St Denis
7342 * LibTomMath is a library that provides multiple-precision
7343 * integer arithmetic as well as number theoretic functionality.
7345 * The library was designed directly after the MPI library by
7346 * Michael Fromberger but has been written from scratch with
7347 * additional optimizations in place.
7349 * The library is free for all purposes without any express
7350 * guarantee it works.
7352 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
7355 /* store in unsigned [big endian] format */
7356 int mp_to_unsigned_bin (mp_int * a, unsigned char *b)
7361 if ((res = mp_init_copy (&t, a)) != MP_OKAY) {
7366 while (mp_iszero (&t) == 0) {
7368 b[x++] = (unsigned char) (t.dp[0] & 255);
7370 b[x++] = (unsigned char) (t.dp[0] | ((t.dp[1] & 0x01) << 7));
7372 if ((res = mp_div_2d (&t, 8, &t, NULL)) != MP_OKAY) {
7383 /* End: bn_mp_to_unsigned_bin.c */
7385 /* Start: bn_mp_to_unsigned_bin_n.c */
7387 #ifdef BN_MP_TO_UNSIGNED_BIN_N_C
7388 /* LibTomMath, multiple-precision integer library -- Tom St Denis
7390 * LibTomMath is a library that provides multiple-precision
7391 * integer arithmetic as well as number theoretic functionality.
7393 * The library was designed directly after the MPI library by
7394 * Michael Fromberger but has been written from scratch with
7395 * additional optimizations in place.
7397 * The library is free for all purposes without any express
7398 * guarantee it works.
7400 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
7403 /* store in unsigned [big endian] format */
7404 int mp_to_unsigned_bin_n (mp_int * a, unsigned char *b, unsigned long *outlen)
7406 if (*outlen < (unsigned long)mp_unsigned_bin_size(a)) {
7409 *outlen = mp_unsigned_bin_size(a);
7410 return mp_to_unsigned_bin(a, b);
7414 /* End: bn_mp_to_unsigned_bin_n.c */
7416 /* Start: bn_mp_toom_mul.c */
7418 #ifdef BN_MP_TOOM_MUL_C
7419 /* LibTomMath, multiple-precision integer library -- Tom St Denis
7421 * LibTomMath is a library that provides multiple-precision
7422 * integer arithmetic as well as number theoretic functionality.
7424 * The library was designed directly after the MPI library by
7425 * Michael Fromberger but has been written from scratch with
7426 * additional optimizations in place.
7428 * The library is free for all purposes without any express
7429 * guarantee it works.
7431 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
7434 /* multiplication using the Toom-Cook 3-way algorithm
7436 * Much more complicated than Karatsuba but has a lower
7437 * asymptotic running time of O(N**1.464). This algorithm is
7438 * only particularly useful on VERY large inputs
7439 * (we're talking 1000s of digits here...).
7441 int mp_toom_mul(mp_int *a, mp_int *b, mp_int *c)
7443 mp_int w0, w1, w2, w3, w4, tmp1, tmp2, a0, a1, a2, b0, b1, b2;
7447 if ((res = mp_init_multi(&w0, &w1, &w2, &w3, &w4,
7448 &a0, &a1, &a2, &b0, &b1,
7449 &b2, &tmp1, &tmp2, NULL)) != MP_OKAY) {
7454 B = MIN(a->used, b->used) / 3;
7456 /* a = a2 * B**2 + a1 * B + a0 */
7457 if ((res = mp_mod_2d(a, DIGIT_BIT * B, &a0)) != MP_OKAY) {
7461 if ((res = mp_copy(a, &a1)) != MP_OKAY) {
7465 mp_mod_2d(&a1, DIGIT_BIT * B, &a1);
7467 if ((res = mp_copy(a, &a2)) != MP_OKAY) {
7472 /* b = b2 * B**2 + b1 * B + b0 */
7473 if ((res = mp_mod_2d(b, DIGIT_BIT * B, &b0)) != MP_OKAY) {
7477 if ((res = mp_copy(b, &b1)) != MP_OKAY) {
7481 mp_mod_2d(&b1, DIGIT_BIT * B, &b1);
7483 if ((res = mp_copy(b, &b2)) != MP_OKAY) {
7489 if ((res = mp_mul(&a0, &b0, &w0)) != MP_OKAY) {
7494 if ((res = mp_mul(&a2, &b2, &w4)) != MP_OKAY) {
7498 /* w1 = (a2 + 2(a1 + 2a0))(b2 + 2(b1 + 2b0)) */
7499 if ((res = mp_mul_2(&a0, &tmp1)) != MP_OKAY) {
7502 if ((res = mp_add(&tmp1, &a1, &tmp1)) != MP_OKAY) {
7505 if ((res = mp_mul_2(&tmp1, &tmp1)) != MP_OKAY) {
7508 if ((res = mp_add(&tmp1, &a2, &tmp1)) != MP_OKAY) {
7512 if ((res = mp_mul_2(&b0, &tmp2)) != MP_OKAY) {
7515 if ((res = mp_add(&tmp2, &b1, &tmp2)) != MP_OKAY) {
7518 if ((res = mp_mul_2(&tmp2, &tmp2)) != MP_OKAY) {
7521 if ((res = mp_add(&tmp2, &b2, &tmp2)) != MP_OKAY) {
7525 if ((res = mp_mul(&tmp1, &tmp2, &w1)) != MP_OKAY) {
7529 /* w3 = (a0 + 2(a1 + 2a2))(b0 + 2(b1 + 2b2)) */
7530 if ((res = mp_mul_2(&a2, &tmp1)) != MP_OKAY) {
7533 if ((res = mp_add(&tmp1, &a1, &tmp1)) != MP_OKAY) {
7536 if ((res = mp_mul_2(&tmp1, &tmp1)) != MP_OKAY) {
7539 if ((res = mp_add(&tmp1, &a0, &tmp1)) != MP_OKAY) {
7543 if ((res = mp_mul_2(&b2, &tmp2)) != MP_OKAY) {
7546 if ((res = mp_add(&tmp2, &b1, &tmp2)) != MP_OKAY) {
7549 if ((res = mp_mul_2(&tmp2, &tmp2)) != MP_OKAY) {
7552 if ((res = mp_add(&tmp2, &b0, &tmp2)) != MP_OKAY) {
7556 if ((res = mp_mul(&tmp1, &tmp2, &w3)) != MP_OKAY) {
7561 /* w2 = (a2 + a1 + a0)(b2 + b1 + b0) */
7562 if ((res = mp_add(&a2, &a1, &tmp1)) != MP_OKAY) {
7565 if ((res = mp_add(&tmp1, &a0, &tmp1)) != MP_OKAY) {
7568 if ((res = mp_add(&b2, &b1, &tmp2)) != MP_OKAY) {
7571 if ((res = mp_add(&tmp2, &b0, &tmp2)) != MP_OKAY) {
7574 if ((res = mp_mul(&tmp1, &tmp2, &w2)) != MP_OKAY) {
7578 /* now solve the matrix
7586 using 12 subtractions, 4 shifts,
7587 2 small divisions and 1 small multiplication
7591 if ((res = mp_sub(&w1, &w4, &w1)) != MP_OKAY) {
7595 if ((res = mp_sub(&w3, &w0, &w3)) != MP_OKAY) {
7599 if ((res = mp_div_2(&w1, &w1)) != MP_OKAY) {
7603 if ((res = mp_div_2(&w3, &w3)) != MP_OKAY) {
7607 if ((res = mp_sub(&w2, &w0, &w2)) != MP_OKAY) {
7610 if ((res = mp_sub(&w2, &w4, &w2)) != MP_OKAY) {
7614 if ((res = mp_sub(&w1, &w2, &w1)) != MP_OKAY) {
7618 if ((res = mp_sub(&w3, &w2, &w3)) != MP_OKAY) {
7622 if ((res = mp_mul_2d(&w0, 3, &tmp1)) != MP_OKAY) {
7625 if ((res = mp_sub(&w1, &tmp1, &w1)) != MP_OKAY) {
7629 if ((res = mp_mul_2d(&w4, 3, &tmp1)) != MP_OKAY) {
7632 if ((res = mp_sub(&w3, &tmp1, &w3)) != MP_OKAY) {
7636 if ((res = mp_mul_d(&w2, 3, &w2)) != MP_OKAY) {
7639 if ((res = mp_sub(&w2, &w1, &w2)) != MP_OKAY) {
7642 if ((res = mp_sub(&w2, &w3, &w2)) != MP_OKAY) {
7646 if ((res = mp_sub(&w1, &w2, &w1)) != MP_OKAY) {
7650 if ((res = mp_sub(&w3, &w2, &w3)) != MP_OKAY) {
7654 if ((res = mp_div_3(&w1, &w1, NULL)) != MP_OKAY) {
7658 if ((res = mp_div_3(&w3, &w3, NULL)) != MP_OKAY) {
7662 /* at this point shift W[n] by B*n */
7663 if ((res = mp_lshd(&w1, 1*B)) != MP_OKAY) {
7666 if ((res = mp_lshd(&w2, 2*B)) != MP_OKAY) {
7669 if ((res = mp_lshd(&w3, 3*B)) != MP_OKAY) {
7672 if ((res = mp_lshd(&w4, 4*B)) != MP_OKAY) {
7676 if ((res = mp_add(&w0, &w1, c)) != MP_OKAY) {
7679 if ((res = mp_add(&w2, &w3, &tmp1)) != MP_OKAY) {
7682 if ((res = mp_add(&w4, &tmp1, &tmp1)) != MP_OKAY) {
7685 if ((res = mp_add(&tmp1, c, c)) != MP_OKAY) {
7690 mp_clear_multi(&w0, &w1, &w2, &w3, &w4,
7691 &a0, &a1, &a2, &b0, &b1,
7692 &b2, &tmp1, &tmp2, NULL);
7698 /* End: bn_mp_toom_mul.c */
7700 /* Start: bn_mp_toom_sqr.c */
7702 #ifdef BN_MP_TOOM_SQR_C
7703 /* LibTomMath, multiple-precision integer library -- Tom St Denis
7705 * LibTomMath is a library that provides multiple-precision
7706 * integer arithmetic as well as number theoretic functionality.
7708 * The library was designed directly after the MPI library by
7709 * Michael Fromberger but has been written from scratch with
7710 * additional optimizations in place.
7712 * The library is free for all purposes without any express
7713 * guarantee it works.
7715 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
7718 /* squaring using Toom-Cook 3-way algorithm */
7720 mp_toom_sqr(mp_int *a, mp_int *b)
7722 mp_int w0, w1, w2, w3, w4, tmp1, a0, a1, a2;
7726 if ((res = mp_init_multi(&w0, &w1, &w2, &w3, &w4, &a0, &a1, &a2, &tmp1, NULL)) != MP_OKAY) {
7733 /* a = a2 * B**2 + a1 * B + a0 */
7734 if ((res = mp_mod_2d(a, DIGIT_BIT * B, &a0)) != MP_OKAY) {
7738 if ((res = mp_copy(a, &a1)) != MP_OKAY) {
7742 mp_mod_2d(&a1, DIGIT_BIT * B, &a1);
7744 if ((res = mp_copy(a, &a2)) != MP_OKAY) {
7750 if ((res = mp_sqr(&a0, &w0)) != MP_OKAY) {
7755 if ((res = mp_sqr(&a2, &w4)) != MP_OKAY) {
7759 /* w1 = (a2 + 2(a1 + 2a0))**2 */
7760 if ((res = mp_mul_2(&a0, &tmp1)) != MP_OKAY) {
7763 if ((res = mp_add(&tmp1, &a1, &tmp1)) != MP_OKAY) {
7766 if ((res = mp_mul_2(&tmp1, &tmp1)) != MP_OKAY) {
7769 if ((res = mp_add(&tmp1, &a2, &tmp1)) != MP_OKAY) {
7773 if ((res = mp_sqr(&tmp1, &w1)) != MP_OKAY) {
7777 /* w3 = (a0 + 2(a1 + 2a2))**2 */
7778 if ((res = mp_mul_2(&a2, &tmp1)) != MP_OKAY) {
7781 if ((res = mp_add(&tmp1, &a1, &tmp1)) != MP_OKAY) {
7784 if ((res = mp_mul_2(&tmp1, &tmp1)) != MP_OKAY) {
7787 if ((res = mp_add(&tmp1, &a0, &tmp1)) != MP_OKAY) {
7791 if ((res = mp_sqr(&tmp1, &w3)) != MP_OKAY) {
7796 /* w2 = (a2 + a1 + a0)**2 */
7797 if ((res = mp_add(&a2, &a1, &tmp1)) != MP_OKAY) {
7800 if ((res = mp_add(&tmp1, &a0, &tmp1)) != MP_OKAY) {
7803 if ((res = mp_sqr(&tmp1, &w2)) != MP_OKAY) {
7807 /* now solve the matrix
7815 using 12 subtractions, 4 shifts, 2 small divisions and 1 small multiplication.
7819 if ((res = mp_sub(&w1, &w4, &w1)) != MP_OKAY) {
7823 if ((res = mp_sub(&w3, &w0, &w3)) != MP_OKAY) {
7827 if ((res = mp_div_2(&w1, &w1)) != MP_OKAY) {
7831 if ((res = mp_div_2(&w3, &w3)) != MP_OKAY) {
7835 if ((res = mp_sub(&w2, &w0, &w2)) != MP_OKAY) {
7838 if ((res = mp_sub(&w2, &w4, &w2)) != MP_OKAY) {
7842 if ((res = mp_sub(&w1, &w2, &w1)) != MP_OKAY) {
7846 if ((res = mp_sub(&w3, &w2, &w3)) != MP_OKAY) {
7850 if ((res = mp_mul_2d(&w0, 3, &tmp1)) != MP_OKAY) {
7853 if ((res = mp_sub(&w1, &tmp1, &w1)) != MP_OKAY) {
7857 if ((res = mp_mul_2d(&w4, 3, &tmp1)) != MP_OKAY) {
7860 if ((res = mp_sub(&w3, &tmp1, &w3)) != MP_OKAY) {
7864 if ((res = mp_mul_d(&w2, 3, &w2)) != MP_OKAY) {
7867 if ((res = mp_sub(&w2, &w1, &w2)) != MP_OKAY) {
7870 if ((res = mp_sub(&w2, &w3, &w2)) != MP_OKAY) {
7874 if ((res = mp_sub(&w1, &w2, &w1)) != MP_OKAY) {
7878 if ((res = mp_sub(&w3, &w2, &w3)) != MP_OKAY) {
7882 if ((res = mp_div_3(&w1, &w1, NULL)) != MP_OKAY) {
7886 if ((res = mp_div_3(&w3, &w3, NULL)) != MP_OKAY) {
7890 /* at this point shift W[n] by B*n */
7891 if ((res = mp_lshd(&w1, 1*B)) != MP_OKAY) {
7894 if ((res = mp_lshd(&w2, 2*B)) != MP_OKAY) {
7897 if ((res = mp_lshd(&w3, 3*B)) != MP_OKAY) {
7900 if ((res = mp_lshd(&w4, 4*B)) != MP_OKAY) {
7904 if ((res = mp_add(&w0, &w1, b)) != MP_OKAY) {
7907 if ((res = mp_add(&w2, &w3, &tmp1)) != MP_OKAY) {
7910 if ((res = mp_add(&w4, &tmp1, &tmp1)) != MP_OKAY) {
7913 if ((res = mp_add(&tmp1, b, b)) != MP_OKAY) {
7918 mp_clear_multi(&w0, &w1, &w2, &w3, &w4, &a0, &a1, &a2, &tmp1, NULL);
7924 /* End: bn_mp_toom_sqr.c */
7926 /* Start: bn_mp_toradix.c */
7928 #ifdef BN_MP_TORADIX_C
7929 /* LibTomMath, multiple-precision integer library -- Tom St Denis
7931 * LibTomMath is a library that provides multiple-precision
7932 * integer arithmetic as well as number theoretic functionality.
7934 * The library was designed directly after the MPI library by
7935 * Michael Fromberger but has been written from scratch with
7936 * additional optimizations in place.
7938 * The library is free for all purposes without any express
7939 * guarantee it works.
7941 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
7944 /* stores a bignum as a ASCII string in a given radix (2..64) */
7945 int mp_toradix (mp_int * a, char *str, int radix)
7952 /* check range of the radix */
7953 if (radix < 2 || radix > 64) {
7957 /* quick out if its zero */
7958 if (mp_iszero(a) == 1) {
7964 if ((res = mp_init_copy (&t, a)) != MP_OKAY) {
7968 /* if it is negative output a - */
7969 if (t.sign == MP_NEG) {
7976 while (mp_iszero (&t) == 0) {
7977 if ((res = mp_div_d (&t, (mp_digit) radix, &t, &d)) != MP_OKAY) {
7981 *str++ = mp_s_rmap[d];
7985 /* reverse the digits of the string. In this case _s points
7986 * to the first digit [exluding the sign] of the number]
7988 bn_reverse ((unsigned char *)_s, digs);
7990 /* append a NULL so the string is properly terminated */
7999 /* End: bn_mp_toradix.c */
8001 /* Start: bn_mp_toradix_n.c */
8003 #ifdef BN_MP_TORADIX_N_C
8004 /* LibTomMath, multiple-precision integer library -- Tom St Denis
8006 * LibTomMath is a library that provides multiple-precision
8007 * integer arithmetic as well as number theoretic functionality.
8009 * The library was designed directly after the MPI library by
8010 * Michael Fromberger but has been written from scratch with
8011 * additional optimizations in place.
8013 * The library is free for all purposes without any express
8014 * guarantee it works.
8016 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
8019 /* stores a bignum as a ASCII string in a given radix (2..64)
8021 * Stores upto maxlen-1 chars and always a NULL byte
8023 int mp_toradix_n(mp_int * a, char *str, int radix, int maxlen)
8030 /* check range of the maxlen, radix */
8031 if (maxlen < 3 || radix < 2 || radix > 64) {
8035 /* quick out if its zero */
8036 if (mp_iszero(a) == 1) {
8042 if ((res = mp_init_copy (&t, a)) != MP_OKAY) {
8046 /* if it is negative output a - */
8047 if (t.sign == MP_NEG) {
8048 /* we have to reverse our digits later... but not the - sign!! */
8051 /* store the flag and mark the number as positive */
8055 /* subtract a char */
8060 while (mp_iszero (&t) == 0) {
8061 if ((res = mp_div_d (&t, (mp_digit) radix, &t, &d)) != MP_OKAY) {
8065 *str++ = mp_s_rmap[d];
8068 if (--maxlen == 1) {
8074 /* reverse the digits of the string. In this case _s points
8075 * to the first digit [exluding the sign] of the number]
8077 bn_reverse ((unsigned char *)_s, digs);
8079 /* append a NULL so the string is properly terminated */
8088 /* End: bn_mp_toradix_n.c */
8090 /* Start: bn_mp_unsigned_bin_size.c */
8092 #ifdef BN_MP_UNSIGNED_BIN_SIZE_C
8093 /* LibTomMath, multiple-precision integer library -- Tom St Denis
8095 * LibTomMath is a library that provides multiple-precision
8096 * integer arithmetic as well as number theoretic functionality.
8098 * The library was designed directly after the MPI library by
8099 * Michael Fromberger but has been written from scratch with
8100 * additional optimizations in place.
8102 * The library is free for all purposes without any express
8103 * guarantee it works.
8105 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
8108 /* get the size for an unsigned equivalent */
8109 int mp_unsigned_bin_size (mp_int * a)
8111 int size = mp_count_bits (a);
8112 return (size / 8 + ((size & 7) != 0 ? 1 : 0));
8116 /* End: bn_mp_unsigned_bin_size.c */
8118 /* Start: bn_mp_xor.c */
8121 /* LibTomMath, multiple-precision integer library -- Tom St Denis
8123 * LibTomMath is a library that provides multiple-precision
8124 * integer arithmetic as well as number theoretic functionality.
8126 * The library was designed directly after the MPI library by
8127 * Michael Fromberger but has been written from scratch with
8128 * additional optimizations in place.
8130 * The library is free for all purposes without any express
8131 * guarantee it works.
8133 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
8136 /* XOR two ints together */
8138 mp_xor (mp_int * a, mp_int * b, mp_int * c)
8143 if (a->used > b->used) {
8144 if ((res = mp_init_copy (&t, a)) != MP_OKAY) {
8150 if ((res = mp_init_copy (&t, b)) != MP_OKAY) {
8157 for (ix = 0; ix < px; ix++) {
8158 t.dp[ix] ^= x->dp[ix];
8167 /* End: bn_mp_xor.c */
8169 /* Start: bn_mp_zero.c */
8172 /* LibTomMath, multiple-precision integer library -- Tom St Denis
8174 * LibTomMath is a library that provides multiple-precision
8175 * integer arithmetic as well as number theoretic functionality.
8177 * The library was designed directly after the MPI library by
8178 * Michael Fromberger but has been written from scratch with
8179 * additional optimizations in place.
8181 * The library is free for all purposes without any express
8182 * guarantee it works.
8184 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
8188 void mp_zero (mp_int * a)
8197 for (n = 0; n < a->alloc; n++) {
8203 /* End: bn_mp_zero.c */
8205 /* Start: bn_prime_tab.c */
8207 #ifdef BN_PRIME_TAB_C
8208 /* LibTomMath, multiple-precision integer library -- Tom St Denis
8210 * LibTomMath is a library that provides multiple-precision
8211 * integer arithmetic as well as number theoretic functionality.
8213 * The library was designed directly after the MPI library by
8214 * Michael Fromberger but has been written from scratch with
8215 * additional optimizations in place.
8217 * The library is free for all purposes without any express
8218 * guarantee it works.
8220 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
8222 const mp_digit ltm_prime_tab[] = {
8223 0x0002, 0x0003, 0x0005, 0x0007, 0x000B, 0x000D, 0x0011, 0x0013,
8224 0x0017, 0x001D, 0x001F, 0x0025, 0x0029, 0x002B, 0x002F, 0x0035,
8225 0x003B, 0x003D, 0x0043, 0x0047, 0x0049, 0x004F, 0x0053, 0x0059,
8226 0x0061, 0x0065, 0x0067, 0x006B, 0x006D, 0x0071, 0x007F,
8229 0x0089, 0x008B, 0x0095, 0x0097, 0x009D, 0x00A3, 0x00A7, 0x00AD,
8230 0x00B3, 0x00B5, 0x00BF, 0x00C1, 0x00C5, 0x00C7, 0x00D3, 0x00DF,
8231 0x00E3, 0x00E5, 0x00E9, 0x00EF, 0x00F1, 0x00FB, 0x0101, 0x0107,
8232 0x010D, 0x010F, 0x0115, 0x0119, 0x011B, 0x0125, 0x0133, 0x0137,
8234 0x0139, 0x013D, 0x014B, 0x0151, 0x015B, 0x015D, 0x0161, 0x0167,
8235 0x016F, 0x0175, 0x017B, 0x017F, 0x0185, 0x018D, 0x0191, 0x0199,
8236 0x01A3, 0x01A5, 0x01AF, 0x01B1, 0x01B7, 0x01BB, 0x01C1, 0x01C9,
8237 0x01CD, 0x01CF, 0x01D3, 0x01DF, 0x01E7, 0x01EB, 0x01F3, 0x01F7,
8238 0x01FD, 0x0209, 0x020B, 0x021D, 0x0223, 0x022D, 0x0233, 0x0239,
8239 0x023B, 0x0241, 0x024B, 0x0251, 0x0257, 0x0259, 0x025F, 0x0265,
8240 0x0269, 0x026B, 0x0277, 0x0281, 0x0283, 0x0287, 0x028D, 0x0293,
8241 0x0295, 0x02A1, 0x02A5, 0x02AB, 0x02B3, 0x02BD, 0x02C5, 0x02CF,
8243 0x02D7, 0x02DD, 0x02E3, 0x02E7, 0x02EF, 0x02F5, 0x02F9, 0x0301,
8244 0x0305, 0x0313, 0x031D, 0x0329, 0x032B, 0x0335, 0x0337, 0x033B,
8245 0x033D, 0x0347, 0x0355, 0x0359, 0x035B, 0x035F, 0x036D, 0x0371,
8246 0x0373, 0x0377, 0x038B, 0x038F, 0x0397, 0x03A1, 0x03A9, 0x03AD,
8247 0x03B3, 0x03B9, 0x03C7, 0x03CB, 0x03D1, 0x03D7, 0x03DF, 0x03E5,
8248 0x03F1, 0x03F5, 0x03FB, 0x03FD, 0x0407, 0x0409, 0x040F, 0x0419,
8249 0x041B, 0x0425, 0x0427, 0x042D, 0x043F, 0x0443, 0x0445, 0x0449,
8250 0x044F, 0x0455, 0x045D, 0x0463, 0x0469, 0x047F, 0x0481, 0x048B,
8252 0x0493, 0x049D, 0x04A3, 0x04A9, 0x04B1, 0x04BD, 0x04C1, 0x04C7,
8253 0x04CD, 0x04CF, 0x04D5, 0x04E1, 0x04EB, 0x04FD, 0x04FF, 0x0503,
8254 0x0509, 0x050B, 0x0511, 0x0515, 0x0517, 0x051B, 0x0527, 0x0529,
8255 0x052F, 0x0551, 0x0557, 0x055D, 0x0565, 0x0577, 0x0581, 0x058F,
8256 0x0593, 0x0595, 0x0599, 0x059F, 0x05A7, 0x05AB, 0x05AD, 0x05B3,
8257 0x05BF, 0x05C9, 0x05CB, 0x05CF, 0x05D1, 0x05D5, 0x05DB, 0x05E7,
8258 0x05F3, 0x05FB, 0x0607, 0x060D, 0x0611, 0x0617, 0x061F, 0x0623,
8259 0x062B, 0x062F, 0x063D, 0x0641, 0x0647, 0x0649, 0x064D, 0x0653
8264 /* End: bn_prime_tab.c */
8266 /* Start: bn_reverse.c */
8269 /* LibTomMath, multiple-precision integer library -- Tom St Denis
8271 * LibTomMath is a library that provides multiple-precision
8272 * integer arithmetic as well as number theoretic functionality.
8274 * The library was designed directly after the MPI library by
8275 * Michael Fromberger but has been written from scratch with
8276 * additional optimizations in place.
8278 * The library is free for all purposes without any express
8279 * guarantee it works.
8281 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
8284 /* reverse an array, used for radix code */
8286 bn_reverse (unsigned char *s, int len)
8303 /* End: bn_reverse.c */
8305 /* Start: bn_s_mp_add.c */
8307 #ifdef BN_S_MP_ADD_C
8308 /* LibTomMath, multiple-precision integer library -- Tom St Denis
8310 * LibTomMath is a library that provides multiple-precision
8311 * integer arithmetic as well as number theoretic functionality.
8313 * The library was designed directly after the MPI library by
8314 * Michael Fromberger but has been written from scratch with
8315 * additional optimizations in place.
8317 * The library is free for all purposes without any express
8318 * guarantee it works.
8320 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
8323 /* low level addition, based on HAC pp.594, Algorithm 14.7 */
8325 s_mp_add (mp_int * a, mp_int * b, mp_int * c)
8328 int olduse, res, min, max;
8330 /* find sizes, we let |a| <= |b| which means we have to sort
8331 * them. "x" will point to the input with the most digits
8333 if (a->used > b->used) {
8344 if (c->alloc < max + 1) {
8345 if ((res = mp_grow (c, max + 1)) != MP_OKAY) {
8350 /* get old used digit count and set new one */
8355 register mp_digit u, *tmpa, *tmpb, *tmpc;
8358 /* alias for digit pointers */
8369 /* zero the carry */
8371 for (i = 0; i < min; i++) {
8372 /* Compute the sum at one digit, T[i] = A[i] + B[i] + U */
8373 *tmpc = *tmpa++ + *tmpb++ + u;
8375 /* U = carry bit of T[i] */
8376 u = *tmpc >> ((mp_digit)DIGIT_BIT);
8378 /* take away carry bit from T[i] */
8382 /* now copy higher words if any, that is in A+B
8383 * if A or B has more digits add those in
8386 for (; i < max; i++) {
8387 /* T[i] = X[i] + U */
8388 *tmpc = x->dp[i] + u;
8390 /* U = carry bit of T[i] */
8391 u = *tmpc >> ((mp_digit)DIGIT_BIT);
8393 /* take away carry bit from T[i] */
8401 /* clear digits above oldused */
8402 for (i = c->used; i < olduse; i++) {
8412 /* End: bn_s_mp_add.c */
8414 /* Start: bn_s_mp_exptmod.c */
8416 #ifdef BN_S_MP_EXPTMOD_C
8417 /* LibTomMath, multiple-precision integer library -- Tom St Denis
8419 * LibTomMath is a library that provides multiple-precision
8420 * integer arithmetic as well as number theoretic functionality.
8422 * The library was designed directly after the MPI library by
8423 * Michael Fromberger but has been written from scratch with
8424 * additional optimizations in place.
8426 * The library is free for all purposes without any express
8427 * guarantee it works.
8429 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
8435 #define TAB_SIZE 256
8438 int s_mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y, int redmode)
8440 mp_int M[TAB_SIZE], res, mu;
8442 int err, bitbuf, bitcpy, bitcnt, mode, digidx, x, y, winsize;
8443 int (*redux)(mp_int*,mp_int*,mp_int*);
8445 /* find window size */
8446 x = mp_count_bits (X);
8449 } else if (x <= 36) {
8451 } else if (x <= 140) {
8453 } else if (x <= 450) {
8455 } else if (x <= 1303) {
8457 } else if (x <= 3529) {
8470 /* init first cell */
8471 if ((err = mp_init(&M[1])) != MP_OKAY) {
8475 /* now init the second half of the array */
8476 for (x = 1<<(winsize-1); x < (1 << winsize); x++) {
8477 if ((err = mp_init(&M[x])) != MP_OKAY) {
8478 for (y = 1<<(winsize-1); y < x; y++) {
8486 /* create mu, used for Barrett reduction */
8487 if ((err = mp_init (&mu)) != MP_OKAY) {
8492 if ((err = mp_reduce_setup (&mu, P)) != MP_OKAY) {
8497 if ((err = mp_reduce_2k_setup_l (P, &mu)) != MP_OKAY) {
8500 redux = mp_reduce_2k_l;
8505 * The M table contains powers of the base,
8506 * e.g. M[x] = G**x mod P
8508 * The first half of the table is not
8509 * computed though accept for M[0] and M[1]
8511 if ((err = mp_mod (G, P, &M[1])) != MP_OKAY) {
8515 /* compute the value at M[1<<(winsize-1)] by squaring
8516 * M[1] (winsize-1) times
8518 if ((err = mp_copy (&M[1], &M[1 << (winsize - 1)])) != MP_OKAY) {
8522 for (x = 0; x < (winsize - 1); x++) {
8524 if ((err = mp_sqr (&M[1 << (winsize - 1)],
8525 &M[1 << (winsize - 1)])) != MP_OKAY) {
8529 /* reduce modulo P */
8530 if ((err = redux (&M[1 << (winsize - 1)], P, &mu)) != MP_OKAY) {
8535 /* create upper table, that is M[x] = M[x-1] * M[1] (mod P)
8536 * for x = (2**(winsize - 1) + 1) to (2**winsize - 1)
8538 for (x = (1 << (winsize - 1)) + 1; x < (1 << winsize); x++) {
8539 if ((err = mp_mul (&M[x - 1], &M[1], &M[x])) != MP_OKAY) {
8542 if ((err = redux (&M[x], P, &mu)) != MP_OKAY) {
8548 if ((err = mp_init (&res)) != MP_OKAY) {
8553 /* set initial mode and bit cnt */
8557 digidx = X->used - 1;
8562 /* grab next digit as required */
8563 if (--bitcnt == 0) {
8564 /* if digidx == -1 we are out of digits */
8568 /* read next digit and reset the bitcnt */
8569 buf = X->dp[digidx--];
8570 bitcnt = (int) DIGIT_BIT;
8573 /* grab the next msb from the exponent */
8574 y = (buf >> (mp_digit)(DIGIT_BIT - 1)) & 1;
8575 buf <<= (mp_digit)1;
8577 /* if the bit is zero and mode == 0 then we ignore it
8578 * These represent the leading zero bits before the first 1 bit
8579 * in the exponent. Technically this opt is not required but it
8580 * does lower the # of trivial squaring/reductions used
8582 if (mode == 0 && y == 0) {
8586 /* if the bit is zero and mode == 1 then we square */
8587 if (mode == 1 && y == 0) {
8588 if ((err = mp_sqr (&res, &res)) != MP_OKAY) {
8591 if ((err = redux (&res, P, &mu)) != MP_OKAY) {
8597 /* else we add it to the window */
8598 bitbuf |= (y << (winsize - ++bitcpy));
8601 if (bitcpy == winsize) {
8602 /* ok window is filled so square as required and multiply */
8604 for (x = 0; x < winsize; x++) {
8605 if ((err = mp_sqr (&res, &res)) != MP_OKAY) {
8608 if ((err = redux (&res, P, &mu)) != MP_OKAY) {
8614 if ((err = mp_mul (&res, &M[bitbuf], &res)) != MP_OKAY) {
8617 if ((err = redux (&res, P, &mu)) != MP_OKAY) {
8621 /* empty window and reset */
8628 /* if bits remain then square/multiply */
8629 if (mode == 2 && bitcpy > 0) {
8630 /* square then multiply if the bit is set */
8631 for (x = 0; x < bitcpy; x++) {
8632 if ((err = mp_sqr (&res, &res)) != MP_OKAY) {
8635 if ((err = redux (&res, P, &mu)) != MP_OKAY) {
8640 if ((bitbuf & (1 << winsize)) != 0) {
8642 if ((err = mp_mul (&res, &M[1], &res)) != MP_OKAY) {
8645 if ((err = redux (&res, P, &mu)) != MP_OKAY) {
8654 LBL_RES:mp_clear (&res);
8655 LBL_MU:mp_clear (&mu);
8658 for (x = 1<<(winsize-1); x < (1 << winsize); x++) {
8665 /* End: bn_s_mp_exptmod.c */
8667 /* Start: bn_s_mp_mul_digs.c */
8669 #ifdef BN_S_MP_MUL_DIGS_C
8670 /* LibTomMath, multiple-precision integer library -- Tom St Denis
8672 * LibTomMath is a library that provides multiple-precision
8673 * integer arithmetic as well as number theoretic functionality.
8675 * The library was designed directly after the MPI library by
8676 * Michael Fromberger but has been written from scratch with
8677 * additional optimizations in place.
8679 * The library is free for all purposes without any express
8680 * guarantee it works.
8682 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
8685 /* multiplies |a| * |b| and only computes upto digs digits of result
8686 * HAC pp. 595, Algorithm 14.12 Modified so you can control how
8687 * many digits of output are created.
8689 int s_mp_mul_digs (mp_int * a, mp_int * b, mp_int * c, int digs)
8692 int res, pa, pb, ix, iy;
8695 mp_digit tmpx, *tmpt, *tmpy;
8697 /* can we use the fast multiplier? */
8698 if (((digs) < MP_WARRAY) &&
8699 MIN (a->used, b->used) <
8700 (1 << ((CHAR_BIT * sizeof (mp_word)) - (2 * DIGIT_BIT)))) {
8701 return fast_s_mp_mul_digs (a, b, c, digs);
8704 if ((res = mp_init_size (&t, digs)) != MP_OKAY) {
8709 /* compute the digits of the product directly */
8711 for (ix = 0; ix < pa; ix++) {
8712 /* set the carry to zero */
8715 /* limit ourselves to making digs digits of output */
8716 pb = MIN (b->used, digs - ix);
8718 /* setup some aliases */
8719 /* copy of the digit from a used within the nested loop */
8722 /* an alias for the destination shifted ix places */
8725 /* an alias for the digits of b */
8728 /* compute the columns of the output and propagate the carry */
8729 for (iy = 0; iy < pb; iy++) {
8730 /* compute the column as a mp_word */
8731 r = ((mp_word)*tmpt) +
8732 ((mp_word)tmpx) * ((mp_word)*tmpy++) +
8735 /* the new column is the lower part of the result */
8736 *tmpt++ = (mp_digit) (r & ((mp_word) MP_MASK));
8738 /* get the carry word from the result */
8739 u = (mp_digit) (r >> ((mp_word) DIGIT_BIT));
8741 /* set carry if it is placed below digs */
8742 if (ix + iy < digs) {
8755 /* End: bn_s_mp_mul_digs.c */
8757 /* Start: bn_s_mp_mul_high_digs.c */
8759 #ifdef BN_S_MP_MUL_HIGH_DIGS_C
8760 /* LibTomMath, multiple-precision integer library -- Tom St Denis
8762 * LibTomMath is a library that provides multiple-precision
8763 * integer arithmetic as well as number theoretic functionality.
8765 * The library was designed directly after the MPI library by
8766 * Michael Fromberger but has been written from scratch with
8767 * additional optimizations in place.
8769 * The library is free for all purposes without any express
8770 * guarantee it works.
8772 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
8775 /* multiplies |a| * |b| and does not compute the lower digs digits
8776 * [meant to get the higher part of the product]
8779 s_mp_mul_high_digs (mp_int * a, mp_int * b, mp_int * c, int digs)
8782 int res, pa, pb, ix, iy;
8785 mp_digit tmpx, *tmpt, *tmpy;
8787 /* can we use the fast multiplier? */
8788 #ifdef BN_FAST_S_MP_MUL_HIGH_DIGS_C
8789 if (((a->used + b->used + 1) < MP_WARRAY)
8790 && MIN (a->used, b->used) < (1 << ((CHAR_BIT * sizeof (mp_word)) - (2 * DIGIT_BIT)))) {
8791 return fast_s_mp_mul_high_digs (a, b, c, digs);
8795 if ((res = mp_init_size (&t, a->used + b->used + 1)) != MP_OKAY) {
8798 t.used = a->used + b->used + 1;
8802 for (ix = 0; ix < pa; ix++) {
8803 /* clear the carry */
8806 /* left hand side of A[ix] * B[iy] */
8809 /* alias to the address of where the digits will be stored */
8810 tmpt = &(t.dp[digs]);
8812 /* alias for where to read the right hand side from */
8813 tmpy = b->dp + (digs - ix);
8815 for (iy = digs - ix; iy < pb; iy++) {
8816 /* calculate the double precision result */
8817 r = ((mp_word)*tmpt) +
8818 ((mp_word)tmpx) * ((mp_word)*tmpy++) +
8821 /* get the lower part */
8822 *tmpt++ = (mp_digit) (r & ((mp_word) MP_MASK));
8824 /* carry the carry */
8825 u = (mp_digit) (r >> ((mp_word) DIGIT_BIT));
8836 /* End: bn_s_mp_mul_high_digs.c */
8838 /* Start: bn_s_mp_sqr.c */
8840 #ifdef BN_S_MP_SQR_C
8841 /* LibTomMath, multiple-precision integer library -- Tom St Denis
8843 * LibTomMath is a library that provides multiple-precision
8844 * integer arithmetic as well as number theoretic functionality.
8846 * The library was designed directly after the MPI library by
8847 * Michael Fromberger but has been written from scratch with
8848 * additional optimizations in place.
8850 * The library is free for all purposes without any express
8851 * guarantee it works.
8853 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
8856 /* low level squaring, b = a*a, HAC pp.596-597, Algorithm 14.16 */
8857 int s_mp_sqr (mp_int * a, mp_int * b)
8860 int res, ix, iy, pa;
8862 mp_digit u, tmpx, *tmpt;
8865 if ((res = mp_init_size (&t, 2*pa + 1)) != MP_OKAY) {
8869 /* default used is maximum possible size */
8872 for (ix = 0; ix < pa; ix++) {
8873 /* first calculate the digit at 2*ix */
8874 /* calculate double precision result */
8875 r = ((mp_word) t.dp[2*ix]) +
8876 ((mp_word)a->dp[ix])*((mp_word)a->dp[ix]);
8878 /* store lower part in result */
8879 t.dp[ix+ix] = (mp_digit) (r & ((mp_word) MP_MASK));
8882 u = (mp_digit)(r >> ((mp_word) DIGIT_BIT));
8884 /* left hand side of A[ix] * A[iy] */
8887 /* alias for where to store the results */
8888 tmpt = t.dp + (2*ix + 1);
8890 for (iy = ix + 1; iy < pa; iy++) {
8891 /* first calculate the product */
8892 r = ((mp_word)tmpx) * ((mp_word)a->dp[iy]);
8894 /* now calculate the double precision result, note we use
8895 * addition instead of *2 since it's easier to optimize
8897 r = ((mp_word) *tmpt) + r + r + ((mp_word) u);
8899 /* store lower part */
8900 *tmpt++ = (mp_digit) (r & ((mp_word) MP_MASK));
8903 u = (mp_digit)(r >> ((mp_word) DIGIT_BIT));
8905 /* propagate upwards */
8906 while (u != ((mp_digit) 0)) {
8907 r = ((mp_word) *tmpt) + ((mp_word) u);
8908 *tmpt++ = (mp_digit) (r & ((mp_word) MP_MASK));
8909 u = (mp_digit)(r >> ((mp_word) DIGIT_BIT));
8920 /* End: bn_s_mp_sqr.c */
8922 /* Start: bn_s_mp_sub.c */
8924 #ifdef BN_S_MP_SUB_C
8925 /* LibTomMath, multiple-precision integer library -- Tom St Denis
8927 * LibTomMath is a library that provides multiple-precision
8928 * integer arithmetic as well as number theoretic functionality.
8930 * The library was designed directly after the MPI library by
8931 * Michael Fromberger but has been written from scratch with
8932 * additional optimizations in place.
8934 * The library is free for all purposes without any express
8935 * guarantee it works.
8937 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
8940 /* low level subtraction (assumes |a| > |b|), HAC pp.595 Algorithm 14.9 */
8942 s_mp_sub (mp_int * a, mp_int * b, mp_int * c)
8944 int olduse, res, min, max;
8951 if (c->alloc < max) {
8952 if ((res = mp_grow (c, max)) != MP_OKAY) {
8960 register mp_digit u, *tmpa, *tmpb, *tmpc;
8963 /* alias for digit pointers */
8968 /* set carry to zero */
8970 for (i = 0; i < min; i++) {
8971 /* T[i] = A[i] - B[i] - U */
8972 *tmpc = *tmpa++ - *tmpb++ - u;
8974 /* U = carry bit of T[i]
8975 * Note this saves performing an AND operation since
8976 * if a carry does occur it will propagate all the way to the
8977 * MSB. As a result a single shift is enough to get the carry
8979 u = *tmpc >> ((mp_digit)(CHAR_BIT * sizeof (mp_digit) - 1));
8981 /* Clear carry from T[i] */
8985 /* now copy higher words if any, e.g. if A has more digits than B */
8986 for (; i < max; i++) {
8987 /* T[i] = A[i] - U */
8988 *tmpc = *tmpa++ - u;
8990 /* U = carry bit of T[i] */
8991 u = *tmpc >> ((mp_digit)(CHAR_BIT * sizeof (mp_digit) - 1));
8993 /* Clear carry from T[i] */
8997 /* clear digits above used (since we may not have grown result above) */
8998 for (i = c->used; i < olduse; i++) {
9009 /* End: bn_s_mp_sub.c */
9011 /* Start: bncore.c */
9014 /* LibTomMath, multiple-precision integer library -- Tom St Denis
9016 * LibTomMath is a library that provides multiple-precision
9017 * integer arithmetic as well as number theoretic functionality.
9019 * The library was designed directly after the MPI library by
9020 * Michael Fromberger but has been written from scratch with
9021 * additional optimizations in place.
9023 * The library is free for all purposes without any express
9024 * guarantee it works.
9026 * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
9029 /* Known optimal configurations
9031 CPU /Compiler /MUL CUTOFF/SQR CUTOFF
9032 -------------------------------------------------------------
9033 Intel P4 Northwood /GCC v3.4.1 / 88/ 128/LTM 0.32 ;-)
9034 AMD Athlon64 /GCC v3.4.4 / 74/ 124/LTM 0.34
9038 int KARATSUBA_MUL_CUTOFF = 74, /* Min. number of digits before Karatsuba multiplication is used. */
9039 KARATSUBA_SQR_CUTOFF = 124, /* Min. number of digits before Karatsuba squaring is used. */
9041 TOOM_MUL_CUTOFF = 350, /* no optimal values of these are known yet so set em high */
9042 TOOM_SQR_CUTOFF = 400;