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@gmail.com, http://libtom.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 *tma_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";
52 /* Start: bn_fast_tma_mp_invmod.c */
54 #ifdef BN_FAST_MP_INVMOD_C
55 /* LibTomMath, multiple-precision integer library -- Tom St Denis
57 * LibTomMath is a library that provides multiple-precision
58 * integer arithmetic as well as number theoretic functionality.
60 * The library was designed directly after the MPI library by
61 * Michael Fromberger but has been written from scratch with
62 * additional optimizations in place.
64 * The library is free for all purposes without any express
67 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
70 /* computes the modular inverse via binary extended euclidean algorithm,
71 * that is c = 1/a mod b
73 * Based on slow invmod except this is optimized for the case where b is
74 * odd as per HAC Note 14.64 on pp. 610
76 int fast_tma_mp_invmod (tma_mp_int * a, tma_mp_int * b, tma_mp_int * c)
78 tma_mp_int x, y, u, v, B, D;
81 /* 2. [modified] b must be odd */
82 if (tma_mp_iseven (b) == 1) {
86 /* init all our temps */
87 if ((res = tma_mp_init_multi(&x, &y, &u, &v, &B, &D, NULL)) != MP_OKAY) {
91 /* x == modulus, y == value to invert */
92 if ((res = tma_mp_copy (b, &x)) != MP_OKAY) {
97 if ((res = tma_mp_mod (a, b, &y)) != MP_OKAY) {
101 /* 3. u=x, v=y, A=1, B=0, C=0,D=1 */
102 if ((res = tma_mp_copy (&x, &u)) != MP_OKAY) {
105 if ((res = tma_mp_copy (&y, &v)) != MP_OKAY) {
111 /* 4. while u is even do */
112 while (tma_mp_iseven (&u) == 1) {
114 if ((res = tma_mp_div_2 (&u, &u)) != MP_OKAY) {
117 /* 4.2 if B is odd then */
118 if (tma_mp_isodd (&B) == 1) {
119 if ((res = tma_mp_sub (&B, &x, &B)) != MP_OKAY) {
124 if ((res = tma_mp_div_2 (&B, &B)) != MP_OKAY) {
129 /* 5. while v is even do */
130 while (tma_mp_iseven (&v) == 1) {
132 if ((res = tma_mp_div_2 (&v, &v)) != MP_OKAY) {
135 /* 5.2 if D is odd then */
136 if (tma_mp_isodd (&D) == 1) {
138 if ((res = tma_mp_sub (&D, &x, &D)) != MP_OKAY) {
143 if ((res = tma_mp_div_2 (&D, &D)) != MP_OKAY) {
148 /* 6. if u >= v then */
149 if (tma_mp_cmp (&u, &v) != MP_LT) {
150 /* u = u - v, B = B - D */
151 if ((res = tma_mp_sub (&u, &v, &u)) != MP_OKAY) {
155 if ((res = tma_mp_sub (&B, &D, &B)) != MP_OKAY) {
159 /* v - v - u, D = D - B */
160 if ((res = tma_mp_sub (&v, &u, &v)) != MP_OKAY) {
164 if ((res = tma_mp_sub (&D, &B, &D)) != MP_OKAY) {
169 /* if not zero goto step 4 */
170 if (tma_mp_iszero (&u) == 0) {
174 /* now a = C, b = D, gcd == g*v */
176 /* if v != 1 then there is no inverse */
177 if (tma_mp_cmp_d (&v, 1) != MP_EQ) {
182 /* b is now the inverse */
184 while (D.sign == MP_NEG) {
185 if ((res = tma_mp_add (&D, b, &D)) != MP_OKAY) {
193 LBL_ERR:tma_mp_clear_multi (&x, &y, &u, &v, &B, &D, NULL);
202 /* End: bn_fast_tma_mp_invmod.c */
204 /* Start: bn_fast_tma_mp_montgomery_reduce.c */
206 #ifdef BN_FAST_MP_MONTGOMERY_REDUCE_C
207 /* LibTomMath, multiple-precision integer library -- Tom St Denis
209 * LibTomMath is a library that provides multiple-precision
210 * integer arithmetic as well as number theoretic functionality.
212 * The library was designed directly after the MPI library by
213 * Michael Fromberger but has been written from scratch with
214 * additional optimizations in place.
216 * The library is free for all purposes without any express
217 * guarantee it works.
219 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
222 /* computes xR**-1 == x (mod N) via Montgomery Reduction
224 * This is an optimized implementation of montgomery_reduce
225 * which uses the comba method to quickly calculate the columns of the
228 * Based on Algorithm 14.32 on pp.601 of HAC.
230 int fast_tma_mp_montgomery_reduce (tma_mp_int * x, tma_mp_int * n, tma_mp_digit rho)
233 tma_mp_word W[MP_WARRAY];
235 /* get old used count */
238 /* grow a as required */
239 if (x->alloc < n->used + 1) {
240 if ((res = tma_mp_grow (x, n->used + 1)) != MP_OKAY) {
245 /* first we have to get the digits of the input into
246 * an array of double precision words W[...]
249 register tma_mp_word *_W;
250 register tma_mp_digit *tmpx;
252 /* alias for the W[] array */
255 /* alias for the digits of x*/
258 /* copy the digits of a into W[0..a->used-1] */
259 for (ix = 0; ix < x->used; ix++) {
263 /* zero the high words of W[a->used..m->used*2] */
264 for (; ix < n->used * 2 + 1; ix++) {
269 /* now we proceed to zero successive digits
270 * from the least significant upwards
272 for (ix = 0; ix < n->used; ix++) {
273 /* mu = ai * m' mod b
275 * We avoid a double precision multiplication (which isn't required)
276 * by casting the value down to a tma_mp_digit. Note this requires
277 * that W[ix-1] have the carry cleared (see after the inner loop)
279 register tma_mp_digit mu;
280 mu = (tma_mp_digit) (((W[ix] & MP_MASK) * rho) & MP_MASK);
282 /* a = a + mu * m * b**i
284 * This is computed in place and on the fly. The multiplication
285 * by b**i is handled by offseting which columns the results
288 * Note the comba method normally doesn't handle carries in the
289 * inner loop In this case we fix the carry from the previous
290 * column since the Montgomery reduction requires digits of the
291 * result (so far) [see above] to work. This is
292 * handled by fixing up one carry after the inner loop. The
293 * carry fixups are done in order so after these loops the
294 * first m->used words of W[] have the carries fixed
298 register tma_mp_digit *tmpn;
299 register tma_mp_word *_W;
301 /* alias for the digits of the modulus */
304 /* Alias for the columns set by an offset of ix */
308 for (iy = 0; iy < n->used; iy++) {
309 *_W++ += ((tma_mp_word)mu) * ((tma_mp_word)*tmpn++);
313 /* now fix carry for next digit, W[ix+1] */
314 W[ix + 1] += W[ix] >> ((tma_mp_word) DIGIT_BIT);
317 /* now we have to propagate the carries and
318 * shift the words downward [all those least
319 * significant digits we zeroed].
322 register tma_mp_digit *tmpx;
323 register tma_mp_word *_W, *_W1;
325 /* nox fix rest of carries */
327 /* alias for current word */
330 /* alias for next word, where the carry goes */
333 for (; ix <= n->used * 2 + 1; ix++) {
334 *_W++ += *_W1++ >> ((tma_mp_word) DIGIT_BIT);
337 /* copy out, A = A/b**n
339 * The result is A/b**n but instead of converting from an
340 * array of tma_mp_word to tma_mp_digit than calling tma_mp_rshd
341 * we just copy them in the right order
344 /* alias for destination word */
347 /* alias for shifted double precision result */
350 for (ix = 0; ix < n->used + 1; ix++) {
351 *tmpx++ = (tma_mp_digit)(*_W++ & ((tma_mp_word) MP_MASK));
354 /* zero oldused digits, if the input a was larger than
355 * m->used+1 we'll have to clear the digits
357 for (; ix < olduse; ix++) {
362 /* set the max used and clamp */
363 x->used = n->used + 1;
366 /* if A >= m then A = A - m */
367 if (tma_mp_cmp_mag (x, n) != MP_LT) {
368 return s_tma_mp_sub (x, n, x);
378 /* End: bn_fast_tma_mp_montgomery_reduce.c */
380 /* Start: bn_fast_s_tma_mp_mul_digs.c */
382 #ifdef BN_FAST_S_MP_MUL_DIGS_C
383 /* LibTomMath, multiple-precision integer library -- Tom St Denis
385 * LibTomMath is a library that provides multiple-precision
386 * integer arithmetic as well as number theoretic functionality.
388 * The library was designed directly after the MPI library by
389 * Michael Fromberger but has been written from scratch with
390 * additional optimizations in place.
392 * The library is free for all purposes without any express
393 * guarantee it works.
395 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
398 /* Fast (comba) multiplier
400 * This is the fast column-array [comba] multiplier. It is
401 * designed to compute the columns of the product first
402 * then handle the carries afterwards. This has the effect
403 * of making the nested loops that compute the columns very
404 * simple and schedulable on super-scalar processors.
406 * This has been modified to produce a variable number of
407 * digits of output so if say only a half-product is required
408 * you don't have to compute the upper half (a feature
409 * required for fast Barrett reduction).
411 * Based on Algorithm 14.12 on pp.595 of HAC.
414 int fast_s_tma_mp_mul_digs (tma_mp_int * a, tma_mp_int * b, tma_mp_int * c, int digs)
416 int olduse, res, pa, ix, iz;
417 tma_mp_digit W[MP_WARRAY];
418 register tma_mp_word _W;
420 /* grow the destination as required */
421 if (c->alloc < digs) {
422 if ((res = tma_mp_grow (c, digs)) != MP_OKAY) {
427 /* number of output digits to produce */
428 pa = MIN(digs, a->used + b->used);
430 /* clear the carry */
432 for (ix = 0; ix < pa; ix++) {
435 tma_mp_digit *tmpx, *tmpy;
437 /* get offsets into the two bignums */
438 ty = MIN(b->used-1, ix);
441 /* setup temp aliases */
445 /* this is the number of times the loop will iterrate, essentially
446 while (tx++ < a->used && ty-- >= 0) { ... }
448 iy = MIN(a->used-tx, ty+1);
451 for (iz = 0; iz < iy; ++iz) {
452 _W += ((tma_mp_word)*tmpx++)*((tma_mp_word)*tmpy--);
457 W[ix] = ((tma_mp_digit)_W) & MP_MASK;
459 /* make next carry */
460 _W = _W >> ((tma_mp_word)DIGIT_BIT);
468 register tma_mp_digit *tmpc;
470 for (ix = 0; ix < pa+1; ix++) {
471 /* now extract the previous digit [below the carry] */
475 /* clear unused digits [that existed in the old copy of c] */
476 for (; ix < olduse; ix++) {
489 /* End: bn_fast_s_tma_mp_mul_digs.c */
491 /* Start: bn_fast_s_tma_mp_mul_high_digs.c */
493 #ifdef BN_FAST_S_MP_MUL_HIGH_DIGS_C
494 /* LibTomMath, multiple-precision integer library -- Tom St Denis
496 * LibTomMath is a library that provides multiple-precision
497 * integer arithmetic as well as number theoretic functionality.
499 * The library was designed directly after the MPI library by
500 * Michael Fromberger but has been written from scratch with
501 * additional optimizations in place.
503 * The library is free for all purposes without any express
504 * guarantee it works.
506 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
509 /* this is a modified version of fast_s_mul_digs that only produces
510 * output digits *above* digs. See the comments for fast_s_mul_digs
511 * to see how it works.
513 * This is used in the Barrett reduction since for one of the multiplications
514 * only the higher digits were needed. This essentially halves the work.
516 * Based on Algorithm 14.12 on pp.595 of HAC.
518 int fast_s_tma_mp_mul_high_digs (tma_mp_int * a, tma_mp_int * b, tma_mp_int * c, int digs)
520 int olduse, res, pa, ix, iz;
521 tma_mp_digit W[MP_WARRAY];
524 /* grow the destination as required */
525 pa = a->used + b->used;
527 if ((res = tma_mp_grow (c, pa)) != MP_OKAY) {
532 /* number of output digits to produce */
533 pa = a->used + b->used;
535 for (ix = digs; ix < pa; ix++) {
537 tma_mp_digit *tmpx, *tmpy;
539 /* get offsets into the two bignums */
540 ty = MIN(b->used-1, ix);
543 /* setup temp aliases */
547 /* this is the number of times the loop will iterrate, essentially its
548 while (tx++ < a->used && ty-- >= 0) { ... }
550 iy = MIN(a->used-tx, ty+1);
553 for (iz = 0; iz < iy; iz++) {
554 _W += ((tma_mp_word)*tmpx++)*((tma_mp_word)*tmpy--);
558 W[ix] = ((tma_mp_digit)_W) & MP_MASK;
560 /* make next carry */
561 _W = _W >> ((tma_mp_word)DIGIT_BIT);
569 register tma_mp_digit *tmpc;
572 for (ix = digs; ix < pa; ix++) {
573 /* now extract the previous digit [below the carry] */
577 /* clear unused digits [that existed in the old copy of c] */
578 for (; ix < olduse; ix++) {
591 /* End: bn_fast_s_tma_mp_mul_high_digs.c */
593 /* Start: bn_fast_s_tma_mp_sqr.c */
595 #ifdef BN_FAST_S_MP_SQR_C
596 /* LibTomMath, multiple-precision integer library -- Tom St Denis
598 * LibTomMath is a library that provides multiple-precision
599 * integer arithmetic as well as number theoretic functionality.
601 * The library was designed directly after the MPI library by
602 * Michael Fromberger but has been written from scratch with
603 * additional optimizations in place.
605 * The library is free for all purposes without any express
606 * guarantee it works.
608 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
611 /* the jist of squaring...
612 * you do like mult except the offset of the tmpx [one that
613 * starts closer to zero] can't equal the offset of tmpy.
614 * So basically you set up iy like before then you min it with
615 * (ty-tx) so that it never happens. You double all those
616 * you add in the inner loop
618 After that loop you do the squares and add them in.
621 int fast_s_tma_mp_sqr (tma_mp_int * a, tma_mp_int * b)
623 int olduse, res, pa, ix, iz;
624 tma_mp_digit W[MP_WARRAY], *tmpx;
627 /* grow the destination as required */
628 pa = a->used + a->used;
630 if ((res = tma_mp_grow (b, pa)) != MP_OKAY) {
635 /* number of output digits to produce */
637 for (ix = 0; ix < pa; ix++) {
645 /* get offsets into the two bignums */
646 ty = MIN(a->used-1, ix);
649 /* setup temp aliases */
653 /* this is the number of times the loop will iterrate, essentially
654 while (tx++ < a->used && ty-- >= 0) { ... }
656 iy = MIN(a->used-tx, ty+1);
658 /* now for squaring tx can never equal ty
659 * we halve the distance since they approach at a rate of 2x
660 * and we have to round because odd cases need to be executed
662 iy = MIN(iy, (ty-tx+1)>>1);
665 for (iz = 0; iz < iy; iz++) {
666 _W += ((tma_mp_word)*tmpx++)*((tma_mp_word)*tmpy--);
669 /* double the inner product and add carry */
672 /* even columns have the square term in them */
674 _W += ((tma_mp_word)a->dp[ix>>1])*((tma_mp_word)a->dp[ix>>1]);
678 W[ix] = (tma_mp_digit)(_W & MP_MASK);
680 /* make next carry */
681 W1 = _W >> ((tma_mp_word)DIGIT_BIT);
686 b->used = a->used+a->used;
691 for (ix = 0; ix < pa; ix++) {
692 *tmpb++ = W[ix] & MP_MASK;
695 /* clear unused digits [that existed in the old copy of c] */
696 for (; ix < olduse; ix++) {
709 /* End: bn_fast_s_tma_mp_sqr.c */
711 /* Start: bn_tma_mp_2expt.c */
714 /* LibTomMath, multiple-precision integer library -- Tom St Denis
716 * LibTomMath is a library that provides multiple-precision
717 * integer arithmetic as well as number theoretic functionality.
719 * The library was designed directly after the MPI library by
720 * Michael Fromberger but has been written from scratch with
721 * additional optimizations in place.
723 * The library is free for all purposes without any express
724 * guarantee it works.
726 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
731 * Simple algorithm which zeroes the int, grows it then just sets one bit
735 tma_mp_2expt (tma_mp_int * a, int b)
739 /* zero a as per default */
742 /* grow a to accomodate the single bit */
743 if ((res = tma_mp_grow (a, b / DIGIT_BIT + 1)) != MP_OKAY) {
747 /* set the used count of where the bit will go */
748 a->used = b / DIGIT_BIT + 1;
750 /* put the single bit in its place */
751 a->dp[b / DIGIT_BIT] = ((tma_mp_digit)1) << (b % DIGIT_BIT);
761 /* End: bn_tma_mp_2expt.c */
763 /* Start: bn_tma_mp_abs.c */
766 /* LibTomMath, multiple-precision integer library -- Tom St Denis
768 * LibTomMath is a library that provides multiple-precision
769 * integer arithmetic as well as number theoretic functionality.
771 * The library was designed directly after the MPI library by
772 * Michael Fromberger but has been written from scratch with
773 * additional optimizations in place.
775 * The library is free for all purposes without any express
776 * guarantee it works.
778 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
783 * Simple function copies the input and fixes the sign to positive
786 tma_mp_abs (tma_mp_int * a, tma_mp_int * b)
792 if ((res = tma_mp_copy (a, b)) != MP_OKAY) {
797 /* force the sign of b to positive */
808 /* End: bn_tma_mp_abs.c */
810 /* Start: bn_tma_mp_add.c */
813 /* LibTomMath, multiple-precision integer library -- Tom St Denis
815 * LibTomMath is a library that provides multiple-precision
816 * integer arithmetic as well as number theoretic functionality.
818 * The library was designed directly after the MPI library by
819 * Michael Fromberger but has been written from scratch with
820 * additional optimizations in place.
822 * The library is free for all purposes without any express
823 * guarantee it works.
825 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
828 /* high level addition (handles signs) */
829 int tma_mp_add (tma_mp_int * a, tma_mp_int * b, tma_mp_int * c)
833 /* get sign of both inputs */
837 /* handle two cases, not four */
839 /* both positive or both negative */
840 /* add their magnitudes, copy the sign */
842 res = s_tma_mp_add (a, b, c);
844 /* one positive, the other negative */
845 /* subtract the one with the greater magnitude from */
846 /* the one of the lesser magnitude. The result gets */
847 /* the sign of the one with the greater magnitude. */
848 if (tma_mp_cmp_mag (a, b) == MP_LT) {
850 res = s_tma_mp_sub (b, a, c);
853 res = s_tma_mp_sub (a, b, c);
865 /* End: bn_tma_mp_add.c */
867 /* Start: bn_tma_mp_add_d.c */
870 /* LibTomMath, multiple-precision integer library -- Tom St Denis
872 * LibTomMath is a library that provides multiple-precision
873 * integer arithmetic as well as number theoretic functionality.
875 * The library was designed directly after the MPI library by
876 * Michael Fromberger but has been written from scratch with
877 * additional optimizations in place.
879 * The library is free for all purposes without any express
880 * guarantee it works.
882 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
885 /* single digit addition */
887 tma_mp_add_d (tma_mp_int * a, tma_mp_digit b, tma_mp_int * c)
889 int res, ix, oldused;
890 tma_mp_digit *tmpa, *tmpc, mu;
892 /* grow c as required */
893 if (c->alloc < a->used + 1) {
894 if ((res = tma_mp_grow(c, a->used + 1)) != MP_OKAY) {
899 /* if a is negative and |a| >= b, call c = |a| - b */
900 if (a->sign == MP_NEG && (a->used > 1 || a->dp[0] >= b)) {
901 /* temporarily fix sign of a */
905 res = tma_mp_sub_d(a, b, c);
908 a->sign = c->sign = MP_NEG;
916 /* old number of used digits in c */
919 /* sign always positive */
925 /* destination alias */
928 /* if a is positive */
929 if (a->sign == MP_ZPOS) {
930 /* add digit, after this we're propagating
934 mu = *tmpc >> DIGIT_BIT;
937 /* now handle rest of the digits */
938 for (ix = 1; ix < a->used; ix++) {
939 *tmpc = *tmpa++ + mu;
940 mu = *tmpc >> DIGIT_BIT;
943 /* set final carry */
948 c->used = a->used + 1;
950 /* a was negative and |a| < b */
953 /* the result is a single digit */
955 *tmpc++ = b - a->dp[0];
960 /* setup count so the clearing of oldused
961 * can fall through correctly
966 /* now zero to oldused */
967 while (ix++ < oldused) {
981 /* End: bn_tma_mp_add_d.c */
983 /* Start: bn_tma_mp_addmod.c */
985 #ifdef BN_MP_ADDMOD_C
986 /* LibTomMath, multiple-precision integer library -- Tom St Denis
988 * LibTomMath is a library that provides multiple-precision
989 * integer arithmetic as well as number theoretic functionality.
991 * The library was designed directly after the MPI library by
992 * Michael Fromberger but has been written from scratch with
993 * additional optimizations in place.
995 * The library is free for all purposes without any express
996 * guarantee it works.
998 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
1001 /* d = a + b (mod c) */
1003 tma_mp_addmod (tma_mp_int * a, tma_mp_int * b, tma_mp_int * c, tma_mp_int * d)
1008 if ((res = tma_mp_init (&t)) != MP_OKAY) {
1012 if ((res = tma_mp_add (a, b, &t)) != MP_OKAY) {
1016 res = tma_mp_mod (&t, c, d);
1026 /* End: bn_tma_mp_addmod.c */
1028 /* Start: bn_tma_mp_and.c */
1031 /* LibTomMath, multiple-precision integer library -- Tom St Denis
1033 * LibTomMath is a library that provides multiple-precision
1034 * integer arithmetic as well as number theoretic functionality.
1036 * The library was designed directly after the MPI library by
1037 * Michael Fromberger but has been written from scratch with
1038 * additional optimizations in place.
1040 * The library is free for all purposes without any express
1041 * guarantee it works.
1043 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
1046 /* AND two ints together */
1048 tma_mp_and (tma_mp_int * a, tma_mp_int * b, tma_mp_int * c)
1053 if (a->used > b->used) {
1054 if ((res = tma_mp_init_copy (&t, a)) != MP_OKAY) {
1060 if ((res = tma_mp_init_copy (&t, b)) != MP_OKAY) {
1067 for (ix = 0; ix < px; ix++) {
1068 t.dp[ix] &= x->dp[ix];
1071 /* zero digits above the last from the smallest tma_mp_int */
1072 for (; ix < t.used; ix++) {
1077 tma_mp_exch (c, &t);
1087 /* End: bn_tma_mp_and.c */
1089 /* Start: bn_tma_mp_clamp.c */
1091 #ifdef BN_MP_CLAMP_C
1092 /* LibTomMath, multiple-precision integer library -- Tom St Denis
1094 * LibTomMath is a library that provides multiple-precision
1095 * integer arithmetic as well as number theoretic functionality.
1097 * The library was designed directly after the MPI library by
1098 * Michael Fromberger but has been written from scratch with
1099 * additional optimizations in place.
1101 * The library is free for all purposes without any express
1102 * guarantee it works.
1104 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
1107 /* trim unused digits
1109 * This is used to ensure that leading zero digits are
1110 * trimed and the leading "used" digit will be non-zero
1111 * Typically very fast. Also fixes the sign if there
1112 * are no more leading digits
1115 tma_mp_clamp (tma_mp_int * a)
1117 /* decrease used while the most significant digit is
1120 while (a->used > 0 && a->dp[a->used - 1] == 0) {
1124 /* reset the sign flag if used == 0 */
1135 /* End: bn_tma_mp_clamp.c */
1137 /* Start: bn_tma_mp_clear.c */
1139 #ifdef BN_MP_CLEAR_C
1140 /* LibTomMath, multiple-precision integer library -- Tom St Denis
1142 * LibTomMath is a library that provides multiple-precision
1143 * integer arithmetic as well as number theoretic functionality.
1145 * The library was designed directly after the MPI library by
1146 * Michael Fromberger but has been written from scratch with
1147 * additional optimizations in place.
1149 * The library is free for all purposes without any express
1150 * guarantee it works.
1152 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
1155 /* clear one (frees) */
1157 tma_mp_clear (tma_mp_int * a)
1161 /* only do anything if a hasn't been freed previously */
1162 if (a->dp != NULL) {
1163 /* first zero the digits */
1164 for (i = 0; i < a->used; i++) {
1171 /* reset members to make debugging easier */
1173 a->alloc = a->used = 0;
1183 /* End: bn_tma_mp_clear.c */
1185 /* Start: bn_tma_mp_clear_multi.c */
1187 #ifdef BN_MP_CLEAR_MULTI_C
1188 /* LibTomMath, multiple-precision integer library -- Tom St Denis
1190 * LibTomMath is a library that provides multiple-precision
1191 * integer arithmetic as well as number theoretic functionality.
1193 * The library was designed directly after the MPI library by
1194 * Michael Fromberger but has been written from scratch with
1195 * additional optimizations in place.
1197 * The library is free for all purposes without any express
1198 * guarantee it works.
1200 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
1204 void tma_mp_clear_multi(tma_mp_int *mp, ...)
1206 tma_mp_int* next_mp = mp;
1209 while (next_mp != NULL) {
1210 tma_mp_clear(next_mp);
1211 next_mp = va_arg(args, tma_mp_int*);
1221 /* End: bn_tma_mp_clear_multi.c */
1223 /* Start: bn_tma_mp_cmp.c */
1226 /* LibTomMath, multiple-precision integer library -- Tom St Denis
1228 * LibTomMath is a library that provides multiple-precision
1229 * integer arithmetic as well as number theoretic functionality.
1231 * The library was designed directly after the MPI library by
1232 * Michael Fromberger but has been written from scratch with
1233 * additional optimizations in place.
1235 * The library is free for all purposes without any express
1236 * guarantee it works.
1238 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
1241 /* compare two ints (signed)*/
1243 tma_mp_cmp (tma_mp_int * a, tma_mp_int * b)
1245 /* compare based on sign */
1246 if (a->sign != b->sign) {
1247 if (a->sign == MP_NEG) {
1254 /* compare digits */
1255 if (a->sign == MP_NEG) {
1256 /* if negative compare opposite direction */
1257 return tma_mp_cmp_mag(b, a);
1259 return tma_mp_cmp_mag(a, b);
1268 /* End: bn_tma_mp_cmp.c */
1270 /* Start: bn_tma_mp_cmp_d.c */
1272 #ifdef BN_MP_CMP_D_C
1273 /* LibTomMath, multiple-precision integer library -- Tom St Denis
1275 * LibTomMath is a library that provides multiple-precision
1276 * integer arithmetic as well as number theoretic functionality.
1278 * The library was designed directly after the MPI library by
1279 * Michael Fromberger but has been written from scratch with
1280 * additional optimizations in place.
1282 * The library is free for all purposes without any express
1283 * guarantee it works.
1285 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
1288 /* compare a digit */
1289 int tma_mp_cmp_d(tma_mp_int * a, tma_mp_digit b)
1291 /* compare based on sign */
1292 if (a->sign == MP_NEG) {
1296 /* compare based on magnitude */
1301 /* compare the only digit of a to b */
1304 } else if (a->dp[0] < b) {
1316 /* End: bn_tma_mp_cmp_d.c */
1318 /* Start: bn_tma_mp_cmp_mag.c */
1320 #ifdef BN_MP_CMP_MAG_C
1321 /* LibTomMath, multiple-precision integer library -- Tom St Denis
1323 * LibTomMath is a library that provides multiple-precision
1324 * integer arithmetic as well as number theoretic functionality.
1326 * The library was designed directly after the MPI library by
1327 * Michael Fromberger but has been written from scratch with
1328 * additional optimizations in place.
1330 * The library is free for all purposes without any express
1331 * guarantee it works.
1333 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
1336 /* compare maginitude of two ints (unsigned) */
1337 int tma_mp_cmp_mag (tma_mp_int * a, tma_mp_int * b)
1340 tma_mp_digit *tmpa, *tmpb;
1342 /* compare based on # of non-zero digits */
1343 if (a->used > b->used) {
1347 if (a->used < b->used) {
1352 tmpa = a->dp + (a->used - 1);
1355 tmpb = b->dp + (a->used - 1);
1357 /* compare based on digits */
1358 for (n = 0; n < a->used; ++n, --tmpa, --tmpb) {
1359 if (*tmpa > *tmpb) {
1363 if (*tmpa < *tmpb) {
1375 /* End: bn_tma_mp_cmp_mag.c */
1377 /* Start: bn_tma_mp_cnt_lsb.c */
1379 #ifdef BN_MP_CNT_LSB_C
1380 /* LibTomMath, multiple-precision integer library -- Tom St Denis
1382 * LibTomMath is a library that provides multiple-precision
1383 * integer arithmetic as well as number theoretic functionality.
1385 * The library was designed directly after the MPI library by
1386 * Michael Fromberger but has been written from scratch with
1387 * additional optimizations in place.
1389 * The library is free for all purposes without any express
1390 * guarantee it works.
1392 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
1395 static const int lnz[16] = {
1396 4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0
1399 /* Counts the number of lsbs which are zero before the first zero bit */
1400 int tma_mp_cnt_lsb(tma_mp_int *a)
1406 if (tma_mp_iszero(a) == 1) {
1410 /* scan lower digits until non-zero */
1411 for (x = 0; x < a->used && a->dp[x] == 0; x++);
1415 /* now scan this digit until a 1 is found */
1432 /* End: bn_tma_mp_cnt_lsb.c */
1434 /* Start: bn_tma_mp_copy.c */
1437 /* LibTomMath, multiple-precision integer library -- Tom St Denis
1439 * LibTomMath is a library that provides multiple-precision
1440 * integer arithmetic as well as number theoretic functionality.
1442 * The library was designed directly after the MPI library by
1443 * Michael Fromberger but has been written from scratch with
1444 * additional optimizations in place.
1446 * The library is free for all purposes without any express
1447 * guarantee it works.
1449 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
1454 tma_mp_copy (tma_mp_int * a, tma_mp_int * b)
1458 /* if dst == src do nothing */
1464 if (b->alloc < a->used) {
1465 if ((res = tma_mp_grow (b, a->used)) != MP_OKAY) {
1470 /* zero b and copy the parameters over */
1472 register tma_mp_digit *tmpa, *tmpb;
1474 /* pointer aliases */
1482 /* copy all the digits */
1483 for (n = 0; n < a->used; n++) {
1487 /* clear high digits */
1488 for (; n < b->used; n++) {
1493 /* copy used count and sign */
1504 /* End: bn_tma_mp_copy.c */
1506 /* Start: bn_tma_mp_count_bits.c */
1508 #ifdef BN_MP_COUNT_BITS_C
1509 /* LibTomMath, multiple-precision integer library -- Tom St Denis
1511 * LibTomMath is a library that provides multiple-precision
1512 * integer arithmetic as well as number theoretic functionality.
1514 * The library was designed directly after the MPI library by
1515 * Michael Fromberger but has been written from scratch with
1516 * additional optimizations in place.
1518 * The library is free for all purposes without any express
1519 * guarantee it works.
1521 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
1524 /* returns the number of bits in an int */
1526 tma_mp_count_bits (tma_mp_int * a)
1536 /* get number of digits and add that */
1537 r = (a->used - 1) * DIGIT_BIT;
1539 /* take the last digit and count the bits in it */
1540 q = a->dp[a->used - 1];
1541 while (q > ((tma_mp_digit) 0)) {
1543 q >>= ((tma_mp_digit) 1);
1553 /* End: bn_tma_mp_count_bits.c */
1555 /* Start: bn_tma_mp_div.c */
1558 /* LibTomMath, multiple-precision integer library -- Tom St Denis
1560 * LibTomMath is a library that provides multiple-precision
1561 * integer arithmetic as well as number theoretic functionality.
1563 * The library was designed directly after the MPI library by
1564 * Michael Fromberger but has been written from scratch with
1565 * additional optimizations in place.
1567 * The library is free for all purposes without any express
1568 * guarantee it works.
1570 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
1573 #ifdef BN_MP_DIV_SMALL
1575 /* slower bit-bang division... also smaller */
1576 int tma_mp_div(tma_mp_int * a, tma_mp_int * b, tma_mp_int * c, tma_mp_int * d)
1578 tma_mp_int ta, tb, tq, q;
1581 /* is divisor zero ? */
1582 if (tma_mp_iszero (b) == 1) {
1586 /* if a < b then q=0, r = a */
1587 if (tma_mp_cmp_mag (a, b) == MP_LT) {
1589 res = tma_mp_copy (a, d);
1599 /* init our temps */
1600 if ((res = tma_mp_init_multi(&ta, &tb, &tq, &q, NULL) != MP_OKAY)) {
1606 n = tma_mp_count_bits(a) - tma_mp_count_bits(b);
1607 if (((res = tma_mp_abs(a, &ta)) != MP_OKAY) ||
1608 ((res = tma_mp_abs(b, &tb)) != MP_OKAY) ||
1609 ((res = tma_mp_mul_2d(&tb, n, &tb)) != MP_OKAY) ||
1610 ((res = tma_mp_mul_2d(&tq, n, &tq)) != MP_OKAY)) {
1615 if (tma_mp_cmp(&tb, &ta) != MP_GT) {
1616 if (((res = tma_mp_sub(&ta, &tb, &ta)) != MP_OKAY) ||
1617 ((res = tma_mp_add(&q, &tq, &q)) != MP_OKAY)) {
1621 if (((res = tma_mp_div_2d(&tb, 1, &tb, NULL)) != MP_OKAY) ||
1622 ((res = tma_mp_div_2d(&tq, 1, &tq, NULL)) != MP_OKAY)) {
1627 /* now q == quotient and ta == remainder */
1629 n2 = (a->sign == b->sign ? MP_ZPOS : MP_NEG);
1632 c->sign = (tma_mp_iszero(c) == MP_YES) ? MP_ZPOS : n2;
1635 tma_mp_exch(d, &ta);
1636 d->sign = (tma_mp_iszero(d) == MP_YES) ? MP_ZPOS : n;
1639 tma_mp_clear_multi(&ta, &tb, &tq, &q, NULL);
1645 /* integer signed division.
1646 * c*b + d == a [e.g. a/b, c=quotient, d=remainder]
1647 * HAC pp.598 Algorithm 14.20
1649 * Note that the description in HAC is horribly
1650 * incomplete. For example, it doesn't consider
1651 * the case where digits are removed from 'x' in
1652 * the inner loop. It also doesn't consider the
1653 * case that y has fewer than three digits, etc..
1655 * The overall algorithm is as described as
1656 * 14.20 from HAC but fixed to treat these cases.
1658 int tma_mp_div (tma_mp_int * a, tma_mp_int * b, tma_mp_int * c, tma_mp_int * d)
1660 tma_mp_int q, x, y, t1, t2;
1661 int res, n, t, i, norm, neg;
1663 /* is divisor zero ? */
1664 if (tma_mp_iszero (b) == 1) {
1668 /* if a < b then q=0, r = a */
1669 if (tma_mp_cmp_mag (a, b) == MP_LT) {
1671 res = tma_mp_copy (a, d);
1681 if ((res = tma_mp_init_size (&q, a->used + 2)) != MP_OKAY) {
1684 q.used = a->used + 2;
1686 if ((res = tma_mp_init (&t1)) != MP_OKAY) {
1690 if ((res = tma_mp_init (&t2)) != MP_OKAY) {
1694 if ((res = tma_mp_init_copy (&x, a)) != MP_OKAY) {
1698 if ((res = tma_mp_init_copy (&y, b)) != MP_OKAY) {
1703 neg = (a->sign == b->sign) ? MP_ZPOS : MP_NEG;
1704 x.sign = y.sign = MP_ZPOS;
1706 /* normalize both x and y, ensure that y >= b/2, [b == 2**DIGIT_BIT] */
1707 norm = tma_mp_count_bits(&y) % DIGIT_BIT;
1708 if (norm < (int)(DIGIT_BIT-1)) {
1709 norm = (DIGIT_BIT-1) - norm;
1710 if ((res = tma_mp_mul_2d (&x, norm, &x)) != MP_OKAY) {
1713 if ((res = tma_mp_mul_2d (&y, norm, &y)) != MP_OKAY) {
1720 /* note hac does 0 based, so if used==5 then its 0,1,2,3,4, e.g. use 4 */
1724 /* while (x >= y*b**n-t) do { q[n-t] += 1; x -= y*b**{n-t} } */
1725 if ((res = tma_mp_lshd (&y, n - t)) != MP_OKAY) { /* y = y*b**{n-t} */
1729 while (tma_mp_cmp (&x, &y) != MP_LT) {
1731 if ((res = tma_mp_sub (&x, &y, &x)) != MP_OKAY) {
1736 /* reset y by shifting it back down */
1737 tma_mp_rshd (&y, n - t);
1739 /* step 3. for i from n down to (t + 1) */
1740 for (i = n; i >= (t + 1); i--) {
1745 /* step 3.1 if xi == yt then set q{i-t-1} to b-1,
1746 * otherwise set q{i-t-1} to (xi*b + x{i-1})/yt */
1747 if (x.dp[i] == y.dp[t]) {
1748 q.dp[i - t - 1] = ((((tma_mp_digit)1) << DIGIT_BIT) - 1);
1751 tmp = ((tma_mp_word) x.dp[i]) << ((tma_mp_word) DIGIT_BIT);
1752 tmp |= ((tma_mp_word) x.dp[i - 1]);
1753 tmp /= ((tma_mp_word) y.dp[t]);
1754 if (tmp > (tma_mp_word) MP_MASK)
1756 q.dp[i - t - 1] = (tma_mp_digit) (tmp & (tma_mp_word) (MP_MASK));
1759 /* while (q{i-t-1} * (yt * b + y{t-1})) >
1760 xi * b**2 + xi-1 * b + xi-2
1764 q.dp[i - t - 1] = (q.dp[i - t - 1] + 1) & MP_MASK;
1766 q.dp[i - t - 1] = (q.dp[i - t - 1] - 1) & MP_MASK;
1768 /* find left hand */
1770 t1.dp[0] = (t - 1 < 0) ? 0 : y.dp[t - 1];
1773 if ((res = tma_mp_mul_d (&t1, q.dp[i - t - 1], &t1)) != MP_OKAY) {
1777 /* find right hand */
1778 t2.dp[0] = (i - 2 < 0) ? 0 : x.dp[i - 2];
1779 t2.dp[1] = (i - 1 < 0) ? 0 : x.dp[i - 1];
1782 } while (tma_mp_cmp_mag(&t1, &t2) == MP_GT);
1784 /* step 3.3 x = x - q{i-t-1} * y * b**{i-t-1} */
1785 if ((res = tma_mp_mul_d (&y, q.dp[i - t - 1], &t1)) != MP_OKAY) {
1789 if ((res = tma_mp_lshd (&t1, i - t - 1)) != MP_OKAY) {
1793 if ((res = tma_mp_sub (&x, &t1, &x)) != MP_OKAY) {
1797 /* if x < 0 then { x = x + y*b**{i-t-1}; q{i-t-1} -= 1; } */
1798 if (x.sign == MP_NEG) {
1799 if ((res = tma_mp_copy (&y, &t1)) != MP_OKAY) {
1802 if ((res = tma_mp_lshd (&t1, i - t - 1)) != MP_OKAY) {
1805 if ((res = tma_mp_add (&x, &t1, &x)) != MP_OKAY) {
1809 q.dp[i - t - 1] = (q.dp[i - t - 1] - 1UL) & MP_MASK;
1813 /* now q is the quotient and x is the remainder
1814 * [which we have to normalize]
1817 /* get sign before writing to c */
1818 x.sign = x.used == 0 ? MP_ZPOS : a->sign;
1822 tma_mp_exch (&q, c);
1827 tma_mp_div_2d (&x, norm, &x, NULL);
1828 tma_mp_exch (&x, d);
1833 LBL_Y:tma_mp_clear (&y);
1834 LBL_X:tma_mp_clear (&x);
1835 LBL_T2:tma_mp_clear (&t2);
1836 LBL_T1:tma_mp_clear (&t1);
1837 LBL_Q:tma_mp_clear (&q);
1849 /* End: bn_tma_mp_div.c */
1851 /* Start: bn_tma_mp_div_2.c */
1853 #ifdef BN_MP_DIV_2_C
1854 /* LibTomMath, multiple-precision integer library -- Tom St Denis
1856 * LibTomMath is a library that provides multiple-precision
1857 * integer arithmetic as well as number theoretic functionality.
1859 * The library was designed directly after the MPI library by
1860 * Michael Fromberger but has been written from scratch with
1861 * additional optimizations in place.
1863 * The library is free for all purposes without any express
1864 * guarantee it works.
1866 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
1870 int tma_mp_div_2(tma_mp_int * a, tma_mp_int * b)
1872 int x, res, oldused;
1875 if (b->alloc < a->used) {
1876 if ((res = tma_mp_grow (b, a->used)) != MP_OKAY) {
1884 register tma_mp_digit r, rr, *tmpa, *tmpb;
1887 tmpa = a->dp + b->used - 1;
1890 tmpb = b->dp + b->used - 1;
1894 for (x = b->used - 1; x >= 0; x--) {
1895 /* get the carry for the next iteration */
1898 /* shift the current digit, add in carry and store */
1899 *tmpb-- = (*tmpa-- >> 1) | (r << (DIGIT_BIT - 1));
1901 /* forward carry to next iteration */
1905 /* zero excess digits */
1906 tmpb = b->dp + b->used;
1907 for (x = b->used; x < oldused; x++) {
1921 /* End: bn_tma_mp_div_2.c */
1923 /* Start: bn_tma_mp_div_2d.c */
1925 #ifdef BN_MP_DIV_2D_C
1926 /* LibTomMath, multiple-precision integer library -- Tom St Denis
1928 * LibTomMath is a library that provides multiple-precision
1929 * integer arithmetic as well as number theoretic functionality.
1931 * The library was designed directly after the MPI library by
1932 * Michael Fromberger but has been written from scratch with
1933 * additional optimizations in place.
1935 * The library is free for all purposes without any express
1936 * guarantee it works.
1938 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
1941 /* shift right by a certain bit count (store quotient in c, optional remainder in d) */
1942 int tma_mp_div_2d (tma_mp_int * a, int b, tma_mp_int * c, tma_mp_int * d)
1944 tma_mp_digit D, r, rr;
1949 /* if the shift count is <= 0 then we do no work */
1951 res = tma_mp_copy (a, c);
1958 if ((res = tma_mp_init (&t)) != MP_OKAY) {
1962 /* get the remainder */
1964 if ((res = tma_mp_mod_2d (a, b, &t)) != MP_OKAY) {
1971 if ((res = tma_mp_copy (a, c)) != MP_OKAY) {
1976 /* shift by as many digits in the bit count */
1977 if (b >= (int)DIGIT_BIT) {
1978 tma_mp_rshd (c, b / DIGIT_BIT);
1981 /* shift any bit count < DIGIT_BIT */
1982 D = (tma_mp_digit) (b % DIGIT_BIT);
1984 register tma_mp_digit *tmpc, mask, shift;
1987 mask = (((tma_mp_digit)1) << D) - 1;
1990 shift = DIGIT_BIT - D;
1993 tmpc = c->dp + (c->used - 1);
1997 for (x = c->used - 1; x >= 0; x--) {
1998 /* get the lower bits of this word in a temp */
2001 /* shift the current word and mix in the carry bits from the previous word */
2002 *tmpc = (*tmpc >> D) | (r << shift);
2005 /* set the carry to the carry bits of the current word found above */
2011 tma_mp_exch (&t, d);
2022 /* End: bn_tma_mp_div_2d.c */
2024 /* Start: bn_tma_mp_div_3.c */
2026 #ifdef BN_MP_DIV_3_C
2027 /* LibTomMath, multiple-precision integer library -- Tom St Denis
2029 * LibTomMath is a library that provides multiple-precision
2030 * integer arithmetic as well as number theoretic functionality.
2032 * The library was designed directly after the MPI library by
2033 * Michael Fromberger but has been written from scratch with
2034 * additional optimizations in place.
2036 * The library is free for all purposes without any express
2037 * guarantee it works.
2039 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
2042 /* divide by three (based on routine from MPI and the GMP manual) */
2044 tma_mp_div_3 (tma_mp_int * a, tma_mp_int *c, tma_mp_digit * d)
2051 /* b = 2**DIGIT_BIT / 3 */
2052 b = (((tma_mp_word)1) << ((tma_mp_word)DIGIT_BIT)) / ((tma_mp_word)3);
2054 if ((res = tma_mp_init_size(&q, a->used)) != MP_OKAY) {
2061 for (ix = a->used - 1; ix >= 0; ix--) {
2062 w = (w << ((tma_mp_word)DIGIT_BIT)) | ((tma_mp_word)a->dp[ix]);
2065 /* multiply w by [1/3] */
2066 t = (w * ((tma_mp_word)b)) >> ((tma_mp_word)DIGIT_BIT);
2068 /* now subtract 3 * [w/3] from w, to get the remainder */
2071 /* fixup the remainder as required since
2072 * the optimization is not exact.
2081 q.dp[ix] = (tma_mp_digit)t;
2084 /* [optional] store the remainder */
2086 *d = (tma_mp_digit)w;
2089 /* [optional] store the quotient */
2105 /* End: bn_tma_mp_div_3.c */
2107 /* Start: bn_tma_mp_div_d.c */
2109 #ifdef BN_MP_DIV_D_C
2110 /* LibTomMath, multiple-precision integer library -- Tom St Denis
2112 * LibTomMath is a library that provides multiple-precision
2113 * integer arithmetic as well as number theoretic functionality.
2115 * The library was designed directly after the MPI library by
2116 * Michael Fromberger but has been written from scratch with
2117 * additional optimizations in place.
2119 * The library is free for all purposes without any express
2120 * guarantee it works.
2122 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
2125 static int s_is_power_of_two(tma_mp_digit b, int *p)
2129 /* fast return if no power of two */
2130 if ((b==0) || (b & (b-1))) {
2134 for (x = 0; x < DIGIT_BIT; x++) {
2135 if (b == (((tma_mp_digit)1)<<x)) {
2143 /* single digit division (based on routine from MPI) */
2144 int tma_mp_div_d (tma_mp_int * a, tma_mp_digit b, tma_mp_int * c, tma_mp_digit * d)
2151 /* cannot divide by zero */
2157 if (b == 1 || tma_mp_iszero(a) == 1) {
2162 return tma_mp_copy(a, c);
2167 /* power of two ? */
2168 if (s_is_power_of_two(b, &ix) == 1) {
2170 *d = a->dp[0] & ((((tma_mp_digit)1)<<ix) - 1);
2173 return tma_mp_div_2d(a, ix, c, NULL);
2178 #ifdef BN_MP_DIV_3_C
2181 return tma_mp_div_3(a, c, d);
2185 /* no easy answer [c'est la vie]. Just division */
2186 if ((res = tma_mp_init_size(&q, a->used)) != MP_OKAY) {
2193 for (ix = a->used - 1; ix >= 0; ix--) {
2194 w = (w << ((tma_mp_word)DIGIT_BIT)) | ((tma_mp_word)a->dp[ix]);
2197 t = (tma_mp_digit)(w / b);
2198 w -= ((tma_mp_word)t) * ((tma_mp_word)b);
2202 q.dp[ix] = (tma_mp_digit)t;
2206 *d = (tma_mp_digit)w;
2224 /* End: bn_tma_mp_div_d.c */
2226 /* Start: bn_tma_mp_dr_is_modulus.c */
2228 #ifdef BN_MP_DR_IS_MODULUS_C
2229 /* LibTomMath, multiple-precision integer library -- Tom St Denis
2231 * LibTomMath is a library that provides multiple-precision
2232 * integer arithmetic as well as number theoretic functionality.
2234 * The library was designed directly after the MPI library by
2235 * Michael Fromberger but has been written from scratch with
2236 * additional optimizations in place.
2238 * The library is free for all purposes without any express
2239 * guarantee it works.
2241 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
2244 /* determines if a number is a valid DR modulus */
2245 int tma_mp_dr_is_modulus(tma_mp_int *a)
2249 /* must be at least two digits */
2254 /* must be of the form b**k - a [a <= b] so all
2255 * but the first digit must be equal to -1 (mod b).
2257 for (ix = 1; ix < a->used; ix++) {
2258 if (a->dp[ix] != MP_MASK) {
2271 /* End: bn_tma_mp_dr_is_modulus.c */
2273 /* Start: bn_tma_mp_dr_reduce.c */
2275 #ifdef BN_MP_DR_REDUCE_C
2276 /* LibTomMath, multiple-precision integer library -- Tom St Denis
2278 * LibTomMath is a library that provides multiple-precision
2279 * integer arithmetic as well as number theoretic functionality.
2281 * The library was designed directly after the MPI library by
2282 * Michael Fromberger but has been written from scratch with
2283 * additional optimizations in place.
2285 * The library is free for all purposes without any express
2286 * guarantee it works.
2288 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
2291 /* reduce "x" in place modulo "n" using the Diminished Radix algorithm.
2293 * Based on algorithm from the paper
2295 * "Generating Efficient Primes for Discrete Log Cryptosystems"
2296 * Chae Hoon Lim, Pil Joong Lee,
2297 * POSTECH Information Research Laboratories
2299 * The modulus must be of a special format [see manual]
2301 * Has been modified to use algorithm 7.10 from the LTM book instead
2303 * Input x must be in the range 0 <= x <= (n-1)**2
2306 tma_mp_dr_reduce (tma_mp_int * x, tma_mp_int * n, tma_mp_digit k)
2310 tma_mp_digit mu, *tmpx1, *tmpx2;
2312 /* m = digits in modulus */
2315 /* ensure that "x" has at least 2m digits */
2316 if (x->alloc < m + m) {
2317 if ((err = tma_mp_grow (x, m + m)) != MP_OKAY) {
2322 /* top of loop, this is where the code resumes if
2323 * another reduction pass is required.
2326 /* aliases for digits */
2327 /* alias for lower half of x */
2330 /* alias for upper half of x, or x/B**m */
2333 /* set carry to zero */
2336 /* compute (x mod B**m) + k * [x/B**m] inline and inplace */
2337 for (i = 0; i < m; i++) {
2338 r = ((tma_mp_word)*tmpx2++) * ((tma_mp_word)k) + *tmpx1 + mu;
2339 *tmpx1++ = (tma_mp_digit)(r & MP_MASK);
2340 mu = (tma_mp_digit)(r >> ((tma_mp_word)DIGIT_BIT));
2343 /* set final carry */
2346 /* zero words above m */
2347 for (i = m + 1; i < x->used; i++) {
2351 /* clamp, sub and return */
2354 /* if x >= n then subtract and reduce again
2355 * Each successive "recursion" makes the input smaller and smaller.
2357 if (tma_mp_cmp_mag (x, n) != MP_LT) {
2358 s_tma_mp_sub(x, n, x);
2369 /* End: bn_tma_mp_dr_reduce.c */
2371 /* Start: bn_tma_mp_dr_setup.c */
2373 #ifdef BN_MP_DR_SETUP_C
2374 /* LibTomMath, multiple-precision integer library -- Tom St Denis
2376 * LibTomMath is a library that provides multiple-precision
2377 * integer arithmetic as well as number theoretic functionality.
2379 * The library was designed directly after the MPI library by
2380 * Michael Fromberger but has been written from scratch with
2381 * additional optimizations in place.
2383 * The library is free for all purposes without any express
2384 * guarantee it works.
2386 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
2389 /* determines the setup value */
2390 void tma_mp_dr_setup(tma_mp_int *a, tma_mp_digit *d)
2392 /* the casts are required if DIGIT_BIT is one less than
2393 * the number of bits in a tma_mp_digit [e.g. DIGIT_BIT==31]
2395 *d = (tma_mp_digit)((((tma_mp_word)1) << ((tma_mp_word)DIGIT_BIT)) -
2396 ((tma_mp_word)a->dp[0]));
2405 /* End: bn_tma_mp_dr_setup.c */
2407 /* Start: bn_tma_mp_exch.c */
2410 /* LibTomMath, multiple-precision integer library -- Tom St Denis
2412 * LibTomMath is a library that provides multiple-precision
2413 * integer arithmetic as well as number theoretic functionality.
2415 * The library was designed directly after the MPI library by
2416 * Michael Fromberger but has been written from scratch with
2417 * additional optimizations in place.
2419 * The library is free for all purposes without any express
2420 * guarantee it works.
2422 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
2425 /* swap the elements of two integers, for cases where you can't simply swap the
2426 * tma_mp_int pointers around
2429 tma_mp_exch (tma_mp_int * a, tma_mp_int * b)
2443 /* End: bn_tma_mp_exch.c */
2445 /* Start: bn_tma_mp_expt_d.c */
2447 #ifdef BN_MP_EXPT_D_C
2448 /* LibTomMath, multiple-precision integer library -- Tom St Denis
2450 * LibTomMath is a library that provides multiple-precision
2451 * integer arithmetic as well as number theoretic functionality.
2453 * The library was designed directly after the MPI library by
2454 * Michael Fromberger but has been written from scratch with
2455 * additional optimizations in place.
2457 * The library is free for all purposes without any express
2458 * guarantee it works.
2460 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
2463 /* calculate c = a**b using a square-multiply algorithm */
2464 int tma_mp_expt_d (tma_mp_int * a, tma_mp_digit b, tma_mp_int * c)
2469 if ((res = tma_mp_init_copy (&g, a)) != MP_OKAY) {
2473 /* set initial result */
2476 for (x = 0; x < (int) DIGIT_BIT; x++) {
2478 if ((res = tma_mp_sqr (c, c)) != MP_OKAY) {
2483 /* if the bit is set multiply */
2484 if ((b & (tma_mp_digit) (((tma_mp_digit)1) << (DIGIT_BIT - 1))) != 0) {
2485 if ((res = tma_mp_mul (c, &g, c)) != MP_OKAY) {
2491 /* shift to next bit */
2504 /* End: bn_tma_mp_expt_d.c */
2506 /* Start: bn_tma_mp_exptmod.c */
2508 #ifdef BN_MP_EXPTMOD_C
2509 /* LibTomMath, multiple-precision integer library -- Tom St Denis
2511 * LibTomMath is a library that provides multiple-precision
2512 * integer arithmetic as well as number theoretic functionality.
2514 * The library was designed directly after the MPI library by
2515 * Michael Fromberger but has been written from scratch with
2516 * additional optimizations in place.
2518 * The library is free for all purposes without any express
2519 * guarantee it works.
2521 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
2525 /* this is a shell function that calls either the normal or Montgomery
2526 * exptmod functions. Originally the call to the montgomery code was
2527 * embedded in the normal function but that wasted alot of stack space
2528 * for nothing (since 99% of the time the Montgomery code would be called)
2530 int tma_mp_exptmod (tma_mp_int * G, tma_mp_int * X, tma_mp_int * P, tma_mp_int * Y)
2534 /* modulus P must be positive */
2535 if (P->sign == MP_NEG) {
2539 /* if exponent X is negative we have to recurse */
2540 if (X->sign == MP_NEG) {
2541 #ifdef BN_MP_INVMOD_C
2542 tma_mp_int tmpG, tmpX;
2545 /* first compute 1/G mod P */
2546 if ((err = tma_mp_init(&tmpG)) != MP_OKAY) {
2549 if ((err = tma_mp_invmod(G, P, &tmpG)) != MP_OKAY) {
2550 tma_mp_clear(&tmpG);
2555 if ((err = tma_mp_init(&tmpX)) != MP_OKAY) {
2556 tma_mp_clear(&tmpG);
2559 if ((err = tma_mp_abs(X, &tmpX)) != MP_OKAY) {
2560 tma_mp_clear_multi(&tmpG, &tmpX, NULL);
2564 /* and now compute (1/G)**|X| instead of G**X [X < 0] */
2565 err = tma_mp_exptmod(&tmpG, &tmpX, P, Y);
2566 tma_mp_clear_multi(&tmpG, &tmpX, NULL);
2574 /* modified diminished radix reduction */
2575 #if defined(BN_MP_REDUCE_IS_2K_L_C) && defined(BN_MP_REDUCE_2K_L_C) && defined(BN_S_MP_EXPTMOD_C)
2576 if (tma_mp_reduce_is_2k_l(P) == MP_YES) {
2577 return s_tma_mp_exptmod(G, X, P, Y, 1);
2581 #ifdef BN_MP_DR_IS_MODULUS_C
2582 /* is it a DR modulus? */
2583 dr = tma_mp_dr_is_modulus(P);
2589 #ifdef BN_MP_REDUCE_IS_2K_C
2590 /* if not, is it a unrestricted DR modulus? */
2592 dr = tma_mp_reduce_is_2k(P) << 1;
2596 /* if the modulus is odd or dr != 0 use the montgomery method */
2597 #ifdef BN_MP_EXPTMOD_FAST_C
2598 if (tma_mp_isodd (P) == 1 || dr != 0) {
2599 return tma_mp_exptmod_fast (G, X, P, Y, dr);
2602 #ifdef BN_S_MP_EXPTMOD_C
2603 /* otherwise use the generic Barrett reduction technique */
2604 return s_tma_mp_exptmod (G, X, P, Y, 0);
2606 /* no exptmod for evens */
2609 #ifdef BN_MP_EXPTMOD_FAST_C
2620 /* End: bn_tma_mp_exptmod.c */
2622 /* Start: bn_tma_mp_exptmod_fast.c */
2624 #ifdef BN_MP_EXPTMOD_FAST_C
2625 /* LibTomMath, multiple-precision integer library -- Tom St Denis
2627 * LibTomMath is a library that provides multiple-precision
2628 * integer arithmetic as well as number theoretic functionality.
2630 * The library was designed directly after the MPI library by
2631 * Michael Fromberger but has been written from scratch with
2632 * additional optimizations in place.
2634 * The library is free for all purposes without any express
2635 * guarantee it works.
2637 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
2640 /* computes Y == G**X mod P, HAC pp.616, Algorithm 14.85
2642 * Uses a left-to-right k-ary sliding window to compute the modular exponentiation.
2643 * The value of k changes based on the size of the exponent.
2645 * Uses Montgomery or Diminished Radix reduction [whichever appropriate]
2651 #define TAB_SIZE 256
2654 int tma_mp_exptmod_fast (tma_mp_int * G, tma_mp_int * X, tma_mp_int * P, tma_mp_int * Y, int redmode)
2656 tma_mp_int M[TAB_SIZE], res;
2657 tma_mp_digit buf, mp;
2658 int err, bitbuf, bitcpy, bitcnt, mode, digidx, x, y, winsize;
2660 /* use a pointer to the reduction algorithm. This allows us to use
2661 * one of many reduction algorithms without modding the guts of
2662 * the code with if statements everywhere.
2664 int (*redux)(tma_mp_int*,tma_mp_int*,tma_mp_digit);
2666 /* find window size */
2667 x = tma_mp_count_bits (X);
2670 } else if (x <= 36) {
2672 } else if (x <= 140) {
2674 } else if (x <= 450) {
2676 } else if (x <= 1303) {
2678 } else if (x <= 3529) {
2691 /* init first cell */
2692 if ((err = tma_mp_init(&M[1])) != MP_OKAY) {
2696 /* now init the second half of the array */
2697 for (x = 1<<(winsize-1); x < (1 << winsize); x++) {
2698 if ((err = tma_mp_init(&M[x])) != MP_OKAY) {
2699 for (y = 1<<(winsize-1); y < x; y++) {
2700 tma_mp_clear (&M[y]);
2702 tma_mp_clear(&M[1]);
2707 /* determine and setup reduction code */
2709 #ifdef BN_MP_MONTGOMERY_SETUP_C
2710 /* now setup montgomery */
2711 if ((err = tma_mp_montgomery_setup (P, &mp)) != MP_OKAY) {
2719 /* automatically pick the comba one if available (saves quite a few calls/ifs) */
2720 #ifdef BN_FAST_MP_MONTGOMERY_REDUCE_C
2721 if (((P->used * 2 + 1) < MP_WARRAY) &&
2722 P->used < (1 << ((CHAR_BIT * sizeof (tma_mp_word)) - (2 * DIGIT_BIT)))) {
2723 redux = fast_tma_mp_montgomery_reduce;
2727 #ifdef BN_MP_MONTGOMERY_REDUCE_C
2728 /* use slower baseline Montgomery method */
2729 redux = tma_mp_montgomery_reduce;
2735 } else if (redmode == 1) {
2736 #if defined(BN_MP_DR_SETUP_C) && defined(BN_MP_DR_REDUCE_C)
2737 /* setup DR reduction for moduli of the form B**k - b */
2738 tma_mp_dr_setup(P, &mp);
2739 redux = tma_mp_dr_reduce;
2745 #if defined(BN_MP_REDUCE_2K_SETUP_C) && defined(BN_MP_REDUCE_2K_C)
2746 /* setup DR reduction for moduli of the form 2**k - b */
2747 if ((err = tma_mp_reduce_2k_setup(P, &mp)) != MP_OKAY) {
2750 redux = tma_mp_reduce_2k;
2758 if ((err = tma_mp_init (&res)) != MP_OKAY) {
2766 * The first half of the table is not computed though accept for M[0] and M[1]
2770 #ifdef BN_MP_MONTGOMERY_CALC_NORMALIZATION_C
2771 /* now we need R mod m */
2772 if ((err = tma_mp_montgomery_calc_normalization (&res, P)) != MP_OKAY) {
2780 /* now set M[1] to G * R mod m */
2781 if ((err = tma_mp_mulmod (G, &res, P, &M[1])) != MP_OKAY) {
2785 tma_mp_set(&res, 1);
2786 if ((err = tma_mp_mod(G, P, &M[1])) != MP_OKAY) {
2791 /* compute the value at M[1<<(winsize-1)] by squaring M[1] (winsize-1) times */
2792 if ((err = tma_mp_copy (&M[1], &M[1 << (winsize - 1)])) != MP_OKAY) {
2796 for (x = 0; x < (winsize - 1); x++) {
2797 if ((err = tma_mp_sqr (&M[1 << (winsize - 1)], &M[1 << (winsize - 1)])) != MP_OKAY) {
2800 if ((err = redux (&M[1 << (winsize - 1)], P, mp)) != MP_OKAY) {
2805 /* create upper table */
2806 for (x = (1 << (winsize - 1)) + 1; x < (1 << winsize); x++) {
2807 if ((err = tma_mp_mul (&M[x - 1], &M[1], &M[x])) != MP_OKAY) {
2810 if ((err = redux (&M[x], P, mp)) != MP_OKAY) {
2815 /* set initial mode and bit cnt */
2819 digidx = X->used - 1;
2824 /* grab next digit as required */
2825 if (--bitcnt == 0) {
2826 /* if digidx == -1 we are out of digits so break */
2830 /* read next digit and reset bitcnt */
2831 buf = X->dp[digidx--];
2832 bitcnt = (int)DIGIT_BIT;
2835 /* grab the next msb from the exponent */
2836 y = (tma_mp_digit)(buf >> (DIGIT_BIT - 1)) & 1;
2837 buf <<= (tma_mp_digit)1;
2839 /* if the bit is zero and mode == 0 then we ignore it
2840 * These represent the leading zero bits before the first 1 bit
2841 * in the exponent. Technically this opt is not required but it
2842 * does lower the # of trivial squaring/reductions used
2844 if (mode == 0 && y == 0) {
2848 /* if the bit is zero and mode == 1 then we square */
2849 if (mode == 1 && y == 0) {
2850 if ((err = tma_mp_sqr (&res, &res)) != MP_OKAY) {
2853 if ((err = redux (&res, P, mp)) != MP_OKAY) {
2859 /* else we add it to the window */
2860 bitbuf |= (y << (winsize - ++bitcpy));
2863 if (bitcpy == winsize) {
2864 /* ok window is filled so square as required and multiply */
2866 for (x = 0; x < winsize; x++) {
2867 if ((err = tma_mp_sqr (&res, &res)) != MP_OKAY) {
2870 if ((err = redux (&res, P, mp)) != MP_OKAY) {
2876 if ((err = tma_mp_mul (&res, &M[bitbuf], &res)) != MP_OKAY) {
2879 if ((err = redux (&res, P, mp)) != MP_OKAY) {
2883 /* empty window and reset */
2890 /* if bits remain then square/multiply */
2891 if (mode == 2 && bitcpy > 0) {
2892 /* square then multiply if the bit is set */
2893 for (x = 0; x < bitcpy; x++) {
2894 if ((err = tma_mp_sqr (&res, &res)) != MP_OKAY) {
2897 if ((err = redux (&res, P, mp)) != MP_OKAY) {
2901 /* get next bit of the window */
2903 if ((bitbuf & (1 << winsize)) != 0) {
2905 if ((err = tma_mp_mul (&res, &M[1], &res)) != MP_OKAY) {
2908 if ((err = redux (&res, P, mp)) != MP_OKAY) {
2916 /* fixup result if Montgomery reduction is used
2917 * recall that any value in a Montgomery system is
2918 * actually multiplied by R mod n. So we have
2919 * to reduce one more time to cancel out the factor
2922 if ((err = redux(&res, P, mp)) != MP_OKAY) {
2927 /* swap res with Y */
2928 tma_mp_exch (&res, Y);
2930 LBL_RES:tma_mp_clear (&res);
2932 tma_mp_clear(&M[1]);
2933 for (x = 1<<(winsize-1); x < (1 << winsize); x++) {
2934 tma_mp_clear (&M[x]);
2945 /* End: bn_tma_mp_exptmod_fast.c */
2947 /* Start: bn_tma_mp_exteuclid.c */
2949 #ifdef BN_MP_EXTEUCLID_C
2950 /* LibTomMath, multiple-precision integer library -- Tom St Denis
2952 * LibTomMath is a library that provides multiple-precision
2953 * integer arithmetic as well as number theoretic functionality.
2955 * The library was designed directly after the MPI library by
2956 * Michael Fromberger but has been written from scratch with
2957 * additional optimizations in place.
2959 * The library is free for all purposes without any express
2960 * guarantee it works.
2962 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
2965 /* Extended euclidean algorithm of (a, b) produces
2968 int tma_mp_exteuclid(tma_mp_int *a, tma_mp_int *b, tma_mp_int *U1, tma_mp_int *U2, tma_mp_int *U3)
2970 tma_mp_int u1,u2,u3,v1,v2,v3,t1,t2,t3,q,tmp;
2973 if ((err = tma_mp_init_multi(&u1, &u2, &u3, &v1, &v2, &v3, &t1, &t2, &t3, &q, &tmp, NULL)) != MP_OKAY) {
2977 /* initialize, (u1,u2,u3) = (1,0,a) */
2979 if ((err = tma_mp_copy(a, &u3)) != MP_OKAY) { goto _ERR; }
2981 /* initialize, (v1,v2,v3) = (0,1,b) */
2983 if ((err = tma_mp_copy(b, &v3)) != MP_OKAY) { goto _ERR; }
2985 /* loop while v3 != 0 */
2986 while (tma_mp_iszero(&v3) == MP_NO) {
2988 if ((err = tma_mp_div(&u3, &v3, &q, NULL)) != MP_OKAY) { goto _ERR; }
2990 /* (t1,t2,t3) = (u1,u2,u3) - (v1,v2,v3)q */
2991 if ((err = tma_mp_mul(&v1, &q, &tmp)) != MP_OKAY) { goto _ERR; }
2992 if ((err = tma_mp_sub(&u1, &tmp, &t1)) != MP_OKAY) { goto _ERR; }
2993 if ((err = tma_mp_mul(&v2, &q, &tmp)) != MP_OKAY) { goto _ERR; }
2994 if ((err = tma_mp_sub(&u2, &tmp, &t2)) != MP_OKAY) { goto _ERR; }
2995 if ((err = tma_mp_mul(&v3, &q, &tmp)) != MP_OKAY) { goto _ERR; }
2996 if ((err = tma_mp_sub(&u3, &tmp, &t3)) != MP_OKAY) { goto _ERR; }
2998 /* (u1,u2,u3) = (v1,v2,v3) */
2999 if ((err = tma_mp_copy(&v1, &u1)) != MP_OKAY) { goto _ERR; }
3000 if ((err = tma_mp_copy(&v2, &u2)) != MP_OKAY) { goto _ERR; }
3001 if ((err = tma_mp_copy(&v3, &u3)) != MP_OKAY) { goto _ERR; }
3003 /* (v1,v2,v3) = (t1,t2,t3) */
3004 if ((err = tma_mp_copy(&t1, &v1)) != MP_OKAY) { goto _ERR; }
3005 if ((err = tma_mp_copy(&t2, &v2)) != MP_OKAY) { goto _ERR; }
3006 if ((err = tma_mp_copy(&t3, &v3)) != MP_OKAY) { goto _ERR; }
3009 /* make sure U3 >= 0 */
3010 if (u3.sign == MP_NEG) {
3011 tma_mp_neg(&u1, &u1);
3012 tma_mp_neg(&u2, &u2);
3013 tma_mp_neg(&u3, &u3);
3016 /* copy result out */
3017 if (U1 != NULL) { tma_mp_exch(U1, &u1); }
3018 if (U2 != NULL) { tma_mp_exch(U2, &u2); }
3019 if (U3 != NULL) { tma_mp_exch(U3, &u3); }
3022 _ERR: tma_mp_clear_multi(&u1, &u2, &u3, &v1, &v2, &v3, &t1, &t2, &t3, &q, &tmp, NULL);
3031 /* End: bn_tma_mp_exteuclid.c */
3033 /* Start: bn_tma_mp_fread.c */
3035 #ifdef BN_MP_FREAD_C
3036 /* LibTomMath, multiple-precision integer library -- Tom St Denis
3038 * LibTomMath is a library that provides multiple-precision
3039 * integer arithmetic as well as number theoretic functionality.
3041 * The library was designed directly after the MPI library by
3042 * Michael Fromberger but has been written from scratch with
3043 * additional optimizations in place.
3045 * The library is free for all purposes without any express
3046 * guarantee it works.
3048 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
3051 /* read a bigint from a file stream in ASCII */
3052 int tma_mp_fread(tma_mp_int *a, int radix, FILE *stream)
3054 int err, ch, neg, y;
3059 /* if first digit is - then set negative */
3069 /* find y in the radix map */
3070 for (y = 0; y < radix; y++) {
3071 if (tma_mp_s_rmap[y] == ch) {
3079 /* shift up and add */
3080 if ((err = tma_mp_mul_d(a, radix, a)) != MP_OKAY) {
3083 if ((err = tma_mp_add_d(a, y, a)) != MP_OKAY) {
3089 if (tma_mp_cmp_d(a, 0) != MP_EQ) {
3102 /* End: bn_tma_mp_fread.c */
3104 /* Start: bn_tma_mp_fwrite.c */
3106 #ifdef BN_MP_FWRITE_C
3107 /* LibTomMath, multiple-precision integer library -- Tom St Denis
3109 * LibTomMath is a library that provides multiple-precision
3110 * integer arithmetic as well as number theoretic functionality.
3112 * The library was designed directly after the MPI library by
3113 * Michael Fromberger but has been written from scratch with
3114 * additional optimizations in place.
3116 * The library is free for all purposes without any express
3117 * guarantee it works.
3119 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
3122 int tma_mp_fwrite(tma_mp_int *a, int radix, FILE *stream)
3127 if ((err = tma_mp_radix_size(a, radix, &len)) != MP_OKAY) {
3131 buf = OPT_CAST(char) XMALLOC (len);
3136 if ((err = tma_mp_toradix(a, buf, radix)) != MP_OKAY) {
3141 for (x = 0; x < len; x++) {
3142 if (fputc(buf[x], stream) == EOF) {
3158 /* End: bn_tma_mp_fwrite.c */
3160 /* Start: bn_tma_mp_gcd.c */
3163 /* LibTomMath, multiple-precision integer library -- Tom St Denis
3165 * LibTomMath is a library that provides multiple-precision
3166 * integer arithmetic as well as number theoretic functionality.
3168 * The library was designed directly after the MPI library by
3169 * Michael Fromberger but has been written from scratch with
3170 * additional optimizations in place.
3172 * The library is free for all purposes without any express
3173 * guarantee it works.
3175 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
3178 /* Greatest Common Divisor using the binary method */
3179 int tma_mp_gcd (tma_mp_int * a, tma_mp_int * b, tma_mp_int * c)
3182 int k, u_lsb, v_lsb, res;
3184 /* either zero than gcd is the largest */
3185 if (tma_mp_iszero (a) == MP_YES) {
3186 return tma_mp_abs (b, c);
3188 if (tma_mp_iszero (b) == MP_YES) {
3189 return tma_mp_abs (a, c);
3192 /* get copies of a and b we can modify */
3193 if ((res = tma_mp_init_copy (&u, a)) != MP_OKAY) {
3197 if ((res = tma_mp_init_copy (&v, b)) != MP_OKAY) {
3201 /* must be positive for the remainder of the algorithm */
3202 u.sign = v.sign = MP_ZPOS;
3204 /* B1. Find the common power of two for u and v */
3205 u_lsb = tma_mp_cnt_lsb(&u);
3206 v_lsb = tma_mp_cnt_lsb(&v);
3207 k = MIN(u_lsb, v_lsb);
3210 /* divide the power of two out */
3211 if ((res = tma_mp_div_2d(&u, k, &u, NULL)) != MP_OKAY) {
3215 if ((res = tma_mp_div_2d(&v, k, &v, NULL)) != MP_OKAY) {
3220 /* divide any remaining factors of two out */
3222 if ((res = tma_mp_div_2d(&u, u_lsb - k, &u, NULL)) != MP_OKAY) {
3228 if ((res = tma_mp_div_2d(&v, v_lsb - k, &v, NULL)) != MP_OKAY) {
3233 while (tma_mp_iszero(&v) == 0) {
3234 /* make sure v is the largest */
3235 if (tma_mp_cmp_mag(&u, &v) == MP_GT) {
3236 /* swap u and v to make sure v is >= u */
3237 tma_mp_exch(&u, &v);
3240 /* subtract smallest from largest */
3241 if ((res = s_tma_mp_sub(&v, &u, &v)) != MP_OKAY) {
3245 /* Divide out all factors of two */
3246 if ((res = tma_mp_div_2d(&v, tma_mp_cnt_lsb(&v), &v, NULL)) != MP_OKAY) {
3251 /* multiply by 2**k which we divided out at the beginning */
3252 if ((res = tma_mp_mul_2d (&u, k, c)) != MP_OKAY) {
3257 LBL_V:tma_mp_clear (&u);
3258 LBL_U:tma_mp_clear (&v);
3267 /* End: bn_tma_mp_gcd.c */
3269 /* Start: bn_tma_mp_get_int.c */
3271 #ifdef BN_MP_GET_INT_C
3272 /* LibTomMath, multiple-precision integer library -- Tom St Denis
3274 * LibTomMath is a library that provides multiple-precision
3275 * integer arithmetic as well as number theoretic functionality.
3277 * The library was designed directly after the MPI library by
3278 * Michael Fromberger but has been written from scratch with
3279 * additional optimizations in place.
3281 * The library is free for all purposes without any express
3282 * guarantee it works.
3284 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
3287 /* get the lower 32-bits of an tma_mp_int */
3288 unsigned long tma_mp_get_int(tma_mp_int * a)
3297 /* get number of digits of the lsb we have to read */
3298 i = MIN(a->used,(int)((sizeof(unsigned long)*CHAR_BIT+DIGIT_BIT-1)/DIGIT_BIT))-1;
3300 /* get most significant digit of result */
3304 res = (res << DIGIT_BIT) | DIGIT(a,i);
3307 /* force result to 32-bits always so it is consistent on non 32-bit platforms */
3308 return res & 0xFFFFFFFFUL;
3316 /* End: bn_tma_mp_get_int.c */
3318 /* Start: bn_tma_mp_grow.c */
3321 /* LibTomMath, multiple-precision integer library -- Tom St Denis
3323 * LibTomMath is a library that provides multiple-precision
3324 * integer arithmetic as well as number theoretic functionality.
3326 * The library was designed directly after the MPI library by
3327 * Michael Fromberger but has been written from scratch with
3328 * additional optimizations in place.
3330 * The library is free for all purposes without any express
3331 * guarantee it works.
3333 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
3336 /* grow as required */
3337 int tma_mp_grow (tma_mp_int * a, int size)
3342 /* if the alloc size is smaller alloc more ram */
3343 if (a->alloc < size) {
3344 /* ensure there are always at least MP_PREC digits extra on top */
3345 size += (MP_PREC * 2) - (size % MP_PREC);
3347 /* reallocate the array a->dp
3349 * We store the return in a temporary variable
3350 * in case the operation failed we don't want
3351 * to overwrite the dp member of a.
3353 tmp = OPT_CAST(tma_mp_digit) XREALLOC (a->dp, sizeof (tma_mp_digit) * size);
3355 /* reallocation failed but "a" is still valid [can be freed] */
3359 /* reallocation succeeded so set a->dp */
3362 /* zero excess digits */
3365 for (; i < a->alloc; i++) {
3377 /* End: bn_tma_mp_grow.c */
3379 /* Start: bn_tma_mp_init.c */
3382 /* LibTomMath, multiple-precision integer library -- Tom St Denis
3384 * LibTomMath is a library that provides multiple-precision
3385 * integer arithmetic as well as number theoretic functionality.
3387 * The library was designed directly after the MPI library by
3388 * Michael Fromberger but has been written from scratch with
3389 * additional optimizations in place.
3391 * The library is free for all purposes without any express
3392 * guarantee it works.
3394 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
3397 /* init a new tma_mp_int */
3398 int tma_mp_init (tma_mp_int * a)
3402 /* allocate memory required and clear it */
3403 a->dp = OPT_CAST(tma_mp_digit) XMALLOC (sizeof (tma_mp_digit) * MP_PREC);
3404 if (a->dp == NULL) {
3408 /* set the digits to zero */
3409 for (i = 0; i < MP_PREC; i++) {
3413 /* set the used to zero, allocated digits to the default precision
3414 * and sign to positive */
3427 /* End: bn_tma_mp_init.c */
3429 /* Start: bn_tma_mp_init_copy.c */
3431 #ifdef BN_MP_INIT_COPY_C
3432 /* LibTomMath, multiple-precision integer library -- Tom St Denis
3434 * LibTomMath is a library that provides multiple-precision
3435 * integer arithmetic as well as number theoretic functionality.
3437 * The library was designed directly after the MPI library by
3438 * Michael Fromberger but has been written from scratch with
3439 * additional optimizations in place.
3441 * The library is free for all purposes without any express
3442 * guarantee it works.
3444 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
3447 /* creates "a" then copies b into it */
3448 int tma_mp_init_copy (tma_mp_int * a, tma_mp_int * b)
3452 if ((res = tma_mp_init (a)) != MP_OKAY) {
3455 return tma_mp_copy (b, a);
3463 /* End: bn_tma_mp_init_copy.c */
3465 /* Start: bn_tma_mp_init_multi.c */
3467 #ifdef BN_MP_INIT_MULTI_C
3468 /* LibTomMath, multiple-precision integer library -- Tom St Denis
3470 * LibTomMath is a library that provides multiple-precision
3471 * integer arithmetic as well as number theoretic functionality.
3473 * The library was designed directly after the MPI library by
3474 * Michael Fromberger but has been written from scratch with
3475 * additional optimizations in place.
3477 * The library is free for all purposes without any express
3478 * guarantee it works.
3480 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
3484 int tma_mp_init_multi(tma_mp_int *mp, ...)
3486 tma_mp_err res = MP_OKAY; /* Assume ok until proven otherwise */
3487 int n = 0; /* Number of ok inits */
3488 tma_mp_int* cur_arg = mp;
3491 va_start(args, mp); /* init args to next argument from caller */
3492 while (cur_arg != NULL) {
3493 if (tma_mp_init(cur_arg) != MP_OKAY) {
3494 /* Oops - error! Back-track and tma_mp_clear what we already
3495 succeeded in init-ing, then return error.
3499 /* end the current list */
3502 /* now start cleaning up */
3504 va_start(clean_args, mp);
3506 tma_mp_clear(cur_arg);
3507 cur_arg = va_arg(clean_args, tma_mp_int*);
3514 cur_arg = va_arg(args, tma_mp_int*);
3517 return res; /* Assumed ok, if error flagged above. */
3526 /* End: bn_tma_mp_init_multi.c */
3528 /* Start: bn_tma_mp_init_set.c */
3530 #ifdef BN_MP_INIT_SET_C
3531 /* LibTomMath, multiple-precision integer library -- Tom St Denis
3533 * LibTomMath is a library that provides multiple-precision
3534 * integer arithmetic as well as number theoretic functionality.
3536 * The library was designed directly after the MPI library by
3537 * Michael Fromberger but has been written from scratch with
3538 * additional optimizations in place.
3540 * The library is free for all purposes without any express
3541 * guarantee it works.
3543 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
3546 /* initialize and set a digit */
3547 int tma_mp_init_set (tma_mp_int * a, tma_mp_digit b)
3550 if ((err = tma_mp_init(a)) != MP_OKAY) {
3562 /* End: bn_tma_mp_init_set.c */
3564 /* Start: bn_tma_mp_init_set_int.c */
3566 #ifdef BN_MP_INIT_SET_INT_C
3567 /* LibTomMath, multiple-precision integer library -- Tom St Denis
3569 * LibTomMath is a library that provides multiple-precision
3570 * integer arithmetic as well as number theoretic functionality.
3572 * The library was designed directly after the MPI library by
3573 * Michael Fromberger but has been written from scratch with
3574 * additional optimizations in place.
3576 * The library is free for all purposes without any express
3577 * guarantee it works.
3579 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
3582 /* initialize and set a digit */
3583 int tma_mp_init_set_int (tma_mp_int * a, unsigned long b)
3586 if ((err = tma_mp_init(a)) != MP_OKAY) {
3589 return tma_mp_set_int(a, b);
3597 /* End: bn_tma_mp_init_set_int.c */
3599 /* Start: bn_tma_mp_init_size.c */
3601 #ifdef BN_MP_INIT_SIZE_C
3602 /* LibTomMath, multiple-precision integer library -- Tom St Denis
3604 * LibTomMath is a library that provides multiple-precision
3605 * integer arithmetic as well as number theoretic functionality.
3607 * The library was designed directly after the MPI library by
3608 * Michael Fromberger but has been written from scratch with
3609 * additional optimizations in place.
3611 * The library is free for all purposes without any express
3612 * guarantee it works.
3614 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
3617 /* init an tma_mp_init for a given size */
3618 int tma_mp_init_size (tma_mp_int * a, int size)
3622 /* pad size so there are always extra digits */
3623 size += (MP_PREC * 2) - (size % MP_PREC);
3626 a->dp = OPT_CAST(tma_mp_digit) XMALLOC (sizeof (tma_mp_digit) * size);
3627 if (a->dp == NULL) {
3631 /* set the members */
3636 /* zero the digits */
3637 for (x = 0; x < size; x++) {
3649 /* End: bn_tma_mp_init_size.c */
3651 /* Start: bn_tma_mp_invmod.c */
3653 #ifdef BN_MP_INVMOD_C
3654 /* LibTomMath, multiple-precision integer library -- Tom St Denis
3656 * LibTomMath is a library that provides multiple-precision
3657 * integer arithmetic as well as number theoretic functionality.
3659 * The library was designed directly after the MPI library by
3660 * Michael Fromberger but has been written from scratch with
3661 * additional optimizations in place.
3663 * The library is free for all purposes without any express
3664 * guarantee it works.
3666 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
3669 /* hac 14.61, pp608 */
3670 int tma_mp_invmod (tma_mp_int * a, tma_mp_int * b, tma_mp_int * c)
3672 /* b cannot be negative */
3673 if (b->sign == MP_NEG || tma_mp_iszero(b) == 1) {
3677 #ifdef BN_FAST_MP_INVMOD_C
3678 /* if the modulus is odd we can use a faster routine instead */
3679 if (tma_mp_isodd (b) == 1) {
3680 return fast_tma_mp_invmod (a, b, c);
3684 #ifdef BN_MP_INVMOD_SLOW_C
3685 return tma_mp_invmod_slow(a, b, c);
3696 /* End: bn_tma_mp_invmod.c */
3698 /* Start: bn_tma_mp_invmod_slow.c */
3700 #ifdef BN_MP_INVMOD_SLOW_C
3701 /* LibTomMath, multiple-precision integer library -- Tom St Denis
3703 * LibTomMath is a library that provides multiple-precision
3704 * integer arithmetic as well as number theoretic functionality.
3706 * The library was designed directly after the MPI library by
3707 * Michael Fromberger but has been written from scratch with
3708 * additional optimizations in place.
3710 * The library is free for all purposes without any express
3711 * guarantee it works.
3713 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
3716 /* hac 14.61, pp608 */
3717 int tma_mp_invmod_slow (tma_mp_int * a, tma_mp_int * b, tma_mp_int * c)
3719 tma_mp_int x, y, u, v, A, B, C, D;
3722 /* b cannot be negative */
3723 if (b->sign == MP_NEG || tma_mp_iszero(b) == 1) {
3728 if ((res = tma_mp_init_multi(&x, &y, &u, &v,
3729 &A, &B, &C, &D, NULL)) != MP_OKAY) {
3734 if ((res = tma_mp_mod(a, b, &x)) != MP_OKAY) {
3737 if ((res = tma_mp_copy (b, &y)) != MP_OKAY) {
3741 /* 2. [modified] if x,y are both even then return an error! */
3742 if (tma_mp_iseven (&x) == 1 && tma_mp_iseven (&y) == 1) {
3747 /* 3. u=x, v=y, A=1, B=0, C=0,D=1 */
3748 if ((res = tma_mp_copy (&x, &u)) != MP_OKAY) {
3751 if ((res = tma_mp_copy (&y, &v)) != MP_OKAY) {
3758 /* 4. while u is even do */
3759 while (tma_mp_iseven (&u) == 1) {
3761 if ((res = tma_mp_div_2 (&u, &u)) != MP_OKAY) {
3764 /* 4.2 if A or B is odd then */
3765 if (tma_mp_isodd (&A) == 1 || tma_mp_isodd (&B) == 1) {
3766 /* A = (A+y)/2, B = (B-x)/2 */
3767 if ((res = tma_mp_add (&A, &y, &A)) != MP_OKAY) {
3770 if ((res = tma_mp_sub (&B, &x, &B)) != MP_OKAY) {
3774 /* A = A/2, B = B/2 */
3775 if ((res = tma_mp_div_2 (&A, &A)) != MP_OKAY) {
3778 if ((res = tma_mp_div_2 (&B, &B)) != MP_OKAY) {
3783 /* 5. while v is even do */
3784 while (tma_mp_iseven (&v) == 1) {
3786 if ((res = tma_mp_div_2 (&v, &v)) != MP_OKAY) {
3789 /* 5.2 if C or D is odd then */
3790 if (tma_mp_isodd (&C) == 1 || tma_mp_isodd (&D) == 1) {
3791 /* C = (C+y)/2, D = (D-x)/2 */
3792 if ((res = tma_mp_add (&C, &y, &C)) != MP_OKAY) {
3795 if ((res = tma_mp_sub (&D, &x, &D)) != MP_OKAY) {
3799 /* C = C/2, D = D/2 */
3800 if ((res = tma_mp_div_2 (&C, &C)) != MP_OKAY) {
3803 if ((res = tma_mp_div_2 (&D, &D)) != MP_OKAY) {
3808 /* 6. if u >= v then */
3809 if (tma_mp_cmp (&u, &v) != MP_LT) {
3810 /* u = u - v, A = A - C, B = B - D */
3811 if ((res = tma_mp_sub (&u, &v, &u)) != MP_OKAY) {
3815 if ((res = tma_mp_sub (&A, &C, &A)) != MP_OKAY) {
3819 if ((res = tma_mp_sub (&B, &D, &B)) != MP_OKAY) {
3823 /* v - v - u, C = C - A, D = D - B */
3824 if ((res = tma_mp_sub (&v, &u, &v)) != MP_OKAY) {
3828 if ((res = tma_mp_sub (&C, &A, &C)) != MP_OKAY) {
3832 if ((res = tma_mp_sub (&D, &B, &D)) != MP_OKAY) {
3837 /* if not zero goto step 4 */
3838 if (tma_mp_iszero (&u) == 0)
3841 /* now a = C, b = D, gcd == g*v */
3843 /* if v != 1 then there is no inverse */
3844 if (tma_mp_cmp_d (&v, 1) != MP_EQ) {
3849 /* if its too low */
3850 while (tma_mp_cmp_d(&C, 0) == MP_LT) {
3851 if ((res = tma_mp_add(&C, b, &C)) != MP_OKAY) {
3857 while (tma_mp_cmp_mag(&C, b) != MP_LT) {
3858 if ((res = tma_mp_sub(&C, b, &C)) != MP_OKAY) {
3863 /* C is now the inverse */
3864 tma_mp_exch (&C, c);
3866 LBL_ERR:tma_mp_clear_multi (&x, &y, &u, &v, &A, &B, &C, &D, NULL);
3875 /* End: bn_tma_mp_invmod_slow.c */
3877 /* Start: bn_tma_mp_is_square.c */
3879 #ifdef BN_MP_IS_SQUARE_C
3880 /* LibTomMath, multiple-precision integer library -- Tom St Denis
3882 * LibTomMath is a library that provides multiple-precision
3883 * integer arithmetic as well as number theoretic functionality.
3885 * The library was designed directly after the MPI library by
3886 * Michael Fromberger but has been written from scratch with
3887 * additional optimizations in place.
3889 * The library is free for all purposes without any express
3890 * guarantee it works.
3892 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
3895 /* Check if remainders are possible squares - fast exclude non-squares */
3896 static const char rem_128[128] = {
3897 0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
3898 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
3899 1, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
3900 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
3901 0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
3902 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
3903 1, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
3904 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1
3907 static const char rem_105[105] = {
3908 0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1,
3909 0, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1,
3910 0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1,
3911 1, 0, 1, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1,
3912 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1,
3913 1, 1, 1, 1, 0, 1, 0, 1, 1, 0, 0, 1, 1, 1, 1,
3914 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1
3917 /* Store non-zero to ret if arg is square, and zero if not */
3918 int tma_mp_is_square(tma_mp_int *arg,int *ret)
3925 /* Default to Non-square :) */
3928 if (arg->sign == MP_NEG) {
3932 /* digits used? (TSD) */
3933 if (arg->used == 0) {
3937 /* First check mod 128 (suppose that DIGIT_BIT is at least 7) */
3938 if (rem_128[127 & DIGIT(arg,0)] == 1) {
3942 /* Next check mod 105 (3*5*7) */
3943 if ((res = tma_mp_mod_d(arg,105,&c)) != MP_OKAY) {
3946 if (rem_105[c] == 1) {
3951 if ((res = tma_mp_init_set_int(&t,11L*13L*17L*19L*23L*29L*31L)) != MP_OKAY) {
3954 if ((res = tma_mp_mod(arg,&t,&t)) != MP_OKAY) {
3957 r = tma_mp_get_int(&t);
3958 /* Check for other prime modules, note it's not an ERROR but we must
3959 * free "t" so the easiest way is to goto ERR. We know that res
3960 * is already equal to MP_OKAY from the tma_mp_mod call
3962 if ( (1L<<(r%11)) & 0x5C4L ) goto ERR;
3963 if ( (1L<<(r%13)) & 0x9E4L ) goto ERR;
3964 if ( (1L<<(r%17)) & 0x5CE8L ) goto ERR;
3965 if ( (1L<<(r%19)) & 0x4F50CL ) goto ERR;
3966 if ( (1L<<(r%23)) & 0x7ACCA0L ) goto ERR;
3967 if ( (1L<<(r%29)) & 0xC2EDD0CL ) goto ERR;
3968 if ( (1L<<(r%31)) & 0x6DE2B848L ) goto ERR;
3970 /* Final check - is sqr(sqrt(arg)) == arg ? */
3971 if ((res = tma_mp_sqrt(arg,&t)) != MP_OKAY) {
3974 if ((res = tma_mp_sqr(&t,&t)) != MP_OKAY) {
3978 *ret = (tma_mp_cmp_mag(&t,arg) == MP_EQ) ? MP_YES : MP_NO;
3979 ERR:tma_mp_clear(&t);
3988 /* End: bn_tma_mp_is_square.c */
3990 /* Start: bn_tma_mp_jacobi.c */
3992 #ifdef BN_MP_JACOBI_C
3993 /* LibTomMath, multiple-precision integer library -- Tom St Denis
3995 * LibTomMath is a library that provides multiple-precision
3996 * integer arithmetic as well as number theoretic functionality.
3998 * The library was designed directly after the MPI library by
3999 * Michael Fromberger but has been written from scratch with
4000 * additional optimizations in place.
4002 * The library is free for all purposes without any express
4003 * guarantee it works.
4005 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
4008 /* computes the jacobi c = (a | n) (or Legendre if n is prime)
4009 * HAC pp. 73 Algorithm 2.149
4011 int tma_mp_jacobi (tma_mp_int * a, tma_mp_int * p, int *c)
4015 tma_mp_digit residue;
4017 /* if p <= 0 return MP_VAL */
4018 if (tma_mp_cmp_d(p, 0) != MP_GT) {
4022 /* step 1. if a == 0, return 0 */
4023 if (tma_mp_iszero (a) == 1) {
4028 /* step 2. if a == 1, return 1 */
4029 if (tma_mp_cmp_d (a, 1) == MP_EQ) {
4037 /* step 3. write a = a1 * 2**k */
4038 if ((res = tma_mp_init_copy (&a1, a)) != MP_OKAY) {
4042 if ((res = tma_mp_init (&p1)) != MP_OKAY) {
4046 /* divide out larger power of two */
4047 k = tma_mp_cnt_lsb(&a1);
4048 if ((res = tma_mp_div_2d(&a1, k, &a1, NULL)) != MP_OKAY) {
4052 /* step 4. if e is even set s=1 */
4056 /* else set s=1 if p = 1/7 (mod 8) or s=-1 if p = 3/5 (mod 8) */
4057 residue = p->dp[0] & 7;
4059 if (residue == 1 || residue == 7) {
4061 } else if (residue == 3 || residue == 5) {
4066 /* step 5. if p == 3 (mod 4) *and* a1 == 3 (mod 4) then s = -s */
4067 if ( ((p->dp[0] & 3) == 3) && ((a1.dp[0] & 3) == 3)) {
4071 /* if a1 == 1 we're done */
4072 if (tma_mp_cmp_d (&a1, 1) == MP_EQ) {
4076 if ((res = tma_mp_mod (p, &a1, &p1)) != MP_OKAY) {
4079 if ((res = tma_mp_jacobi (&p1, &a1, &r)) != MP_OKAY) {
4087 LBL_P1:tma_mp_clear (&p1);
4088 LBL_A1:tma_mp_clear (&a1);
4097 /* End: bn_tma_mp_jacobi.c */
4099 /* Start: bn_tma_mp_karatsuba_mul.c */
4101 #ifdef BN_MP_KARATSUBA_MUL_C
4102 /* LibTomMath, multiple-precision integer library -- Tom St Denis
4104 * LibTomMath is a library that provides multiple-precision
4105 * integer arithmetic as well as number theoretic functionality.
4107 * The library was designed directly after the MPI library by
4108 * Michael Fromberger but has been written from scratch with
4109 * additional optimizations in place.
4111 * The library is free for all purposes without any express
4112 * guarantee it works.
4114 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
4117 /* c = |a| * |b| using Karatsuba Multiplication using
4118 * three half size multiplications
4120 * Let B represent the radix [e.g. 2**DIGIT_BIT] and
4121 * let n represent half of the number of digits in
4124 * a = a1 * B**n + a0
4125 * b = b1 * B**n + b0
4128 a1b1 * B**2n + ((a1 + a0)(b1 + b0) - (a0b0 + a1b1)) * B + a0b0
4130 * Note that a1b1 and a0b0 are used twice and only need to be
4131 * computed once. So in total three half size (half # of
4132 * digit) multiplications are performed, a0b0, a1b1 and
4135 * Note that a multiplication of half the digits requires
4136 * 1/4th the number of single precision multiplications so in
4137 * total after one call 25% of the single precision multiplications
4138 * are saved. Note also that the call to tma_mp_mul can end up back
4139 * in this function if the a0, a1, b0, or b1 are above the threshold.
4140 * This is known as divide-and-conquer and leads to the famous
4141 * O(N**lg(3)) or O(N**1.584) work which is asymptopically lower than
4142 * the standard O(N**2) that the baseline/comba methods use.
4143 * Generally though the overhead of this method doesn't pay off
4144 * until a certain size (N ~ 80) is reached.
4146 int tma_mp_karatsuba_mul (tma_mp_int * a, tma_mp_int * b, tma_mp_int * c)
4148 tma_mp_int x0, x1, y0, y1, t1, x0y0, x1y1;
4151 /* default the return code to an error */
4154 /* min # of digits */
4155 B = MIN (a->used, b->used);
4157 /* now divide in two */
4160 /* init copy all the temps */
4161 if (tma_mp_init_size (&x0, B) != MP_OKAY)
4163 if (tma_mp_init_size (&x1, a->used - B) != MP_OKAY)
4165 if (tma_mp_init_size (&y0, B) != MP_OKAY)
4167 if (tma_mp_init_size (&y1, b->used - B) != MP_OKAY)
4171 if (tma_mp_init_size (&t1, B * 2) != MP_OKAY)
4173 if (tma_mp_init_size (&x0y0, B * 2) != MP_OKAY)
4175 if (tma_mp_init_size (&x1y1, B * 2) != MP_OKAY)
4178 /* now shift the digits */
4179 x0.used = y0.used = B;
4180 x1.used = a->used - B;
4181 y1.used = b->used - B;
4185 register tma_mp_digit *tmpa, *tmpb, *tmpx, *tmpy;
4187 /* we copy the digits directly instead of using higher level functions
4188 * since we also need to shift the digits
4195 for (x = 0; x < B; x++) {
4201 for (x = B; x < a->used; x++) {
4206 for (x = B; x < b->used; x++) {
4211 /* only need to clamp the lower words since by definition the
4212 * upper words x1/y1 must have a known number of digits
4217 /* now calc the products x0y0 and x1y1 */
4218 /* after this x0 is no longer required, free temp [x0==t2]! */
4219 if (tma_mp_mul (&x0, &y0, &x0y0) != MP_OKAY)
4220 goto X1Y1; /* x0y0 = x0*y0 */
4221 if (tma_mp_mul (&x1, &y1, &x1y1) != MP_OKAY)
4222 goto X1Y1; /* x1y1 = x1*y1 */
4224 /* now calc x1+x0 and y1+y0 */
4225 if (s_tma_mp_add (&x1, &x0, &t1) != MP_OKAY)
4226 goto X1Y1; /* t1 = x1 - x0 */
4227 if (s_tma_mp_add (&y1, &y0, &x0) != MP_OKAY)
4228 goto X1Y1; /* t2 = y1 - y0 */
4229 if (tma_mp_mul (&t1, &x0, &t1) != MP_OKAY)
4230 goto X1Y1; /* t1 = (x1 + x0) * (y1 + y0) */
4233 if (tma_mp_add (&x0y0, &x1y1, &x0) != MP_OKAY)
4234 goto X1Y1; /* t2 = x0y0 + x1y1 */
4235 if (s_tma_mp_sub (&t1, &x0, &t1) != MP_OKAY)
4236 goto X1Y1; /* t1 = (x1+x0)*(y1+y0) - (x1y1 + x0y0) */
4239 if (tma_mp_lshd (&t1, B) != MP_OKAY)
4240 goto X1Y1; /* t1 = (x0y0 + x1y1 - (x1-x0)*(y1-y0))<<B */
4241 if (tma_mp_lshd (&x1y1, B * 2) != MP_OKAY)
4242 goto X1Y1; /* x1y1 = x1y1 << 2*B */
4244 if (tma_mp_add (&x0y0, &t1, &t1) != MP_OKAY)
4245 goto X1Y1; /* t1 = x0y0 + t1 */
4246 if (tma_mp_add (&t1, &x1y1, c) != MP_OKAY)
4247 goto X1Y1; /* t1 = x0y0 + t1 + x1y1 */
4249 /* Algorithm succeeded set the return code to MP_OKAY */
4252 X1Y1:tma_mp_clear (&x1y1);
4253 X0Y0:tma_mp_clear (&x0y0);
4254 T1:tma_mp_clear (&t1);
4255 Y1:tma_mp_clear (&y1);
4256 Y0:tma_mp_clear (&y0);
4257 X1:tma_mp_clear (&x1);
4258 X0:tma_mp_clear (&x0);
4268 /* End: bn_tma_mp_karatsuba_mul.c */
4270 /* Start: bn_tma_mp_karatsuba_sqr.c */
4272 #ifdef BN_MP_KARATSUBA_SQR_C
4273 /* LibTomMath, multiple-precision integer library -- Tom St Denis
4275 * LibTomMath is a library that provides multiple-precision
4276 * integer arithmetic as well as number theoretic functionality.
4278 * The library was designed directly after the MPI library by
4279 * Michael Fromberger but has been written from scratch with
4280 * additional optimizations in place.
4282 * The library is free for all purposes without any express
4283 * guarantee it works.
4285 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
4288 /* Karatsuba squaring, computes b = a*a using three
4289 * half size squarings
4291 * See comments of karatsuba_mul for details. It
4292 * is essentially the same algorithm but merely
4293 * tuned to perform recursive squarings.
4295 int tma_mp_karatsuba_sqr (tma_mp_int * a, tma_mp_int * b)
4297 tma_mp_int x0, x1, t1, t2, x0x0, x1x1;
4302 /* min # of digits */
4305 /* now divide in two */
4308 /* init copy all the temps */
4309 if (tma_mp_init_size (&x0, B) != MP_OKAY)
4311 if (tma_mp_init_size (&x1, a->used - B) != MP_OKAY)
4315 if (tma_mp_init_size (&t1, a->used * 2) != MP_OKAY)
4317 if (tma_mp_init_size (&t2, a->used * 2) != MP_OKAY)
4319 if (tma_mp_init_size (&x0x0, B * 2) != MP_OKAY)
4321 if (tma_mp_init_size (&x1x1, (a->used - B) * 2) != MP_OKAY)
4326 register tma_mp_digit *dst, *src;
4330 /* now shift the digits */
4332 for (x = 0; x < B; x++) {
4337 for (x = B; x < a->used; x++) {
4343 x1.used = a->used - B;
4347 /* now calc the products x0*x0 and x1*x1 */
4348 if (tma_mp_sqr (&x0, &x0x0) != MP_OKAY)
4349 goto X1X1; /* x0x0 = x0*x0 */
4350 if (tma_mp_sqr (&x1, &x1x1) != MP_OKAY)
4351 goto X1X1; /* x1x1 = x1*x1 */
4353 /* now calc (x1+x0)**2 */
4354 if (s_tma_mp_add (&x1, &x0, &t1) != MP_OKAY)
4355 goto X1X1; /* t1 = x1 - x0 */
4356 if (tma_mp_sqr (&t1, &t1) != MP_OKAY)
4357 goto X1X1; /* t1 = (x1 - x0) * (x1 - x0) */
4360 if (s_tma_mp_add (&x0x0, &x1x1, &t2) != MP_OKAY)
4361 goto X1X1; /* t2 = x0x0 + x1x1 */
4362 if (s_tma_mp_sub (&t1, &t2, &t1) != MP_OKAY)
4363 goto X1X1; /* t1 = (x1+x0)**2 - (x0x0 + x1x1) */
4366 if (tma_mp_lshd (&t1, B) != MP_OKAY)
4367 goto X1X1; /* t1 = (x0x0 + x1x1 - (x1-x0)*(x1-x0))<<B */
4368 if (tma_mp_lshd (&x1x1, B * 2) != MP_OKAY)
4369 goto X1X1; /* x1x1 = x1x1 << 2*B */
4371 if (tma_mp_add (&x0x0, &t1, &t1) != MP_OKAY)
4372 goto X1X1; /* t1 = x0x0 + t1 */
4373 if (tma_mp_add (&t1, &x1x1, b) != MP_OKAY)
4374 goto X1X1; /* t1 = x0x0 + t1 + x1x1 */
4378 X1X1:tma_mp_clear (&x1x1);
4379 X0X0:tma_mp_clear (&x0x0);
4380 T2:tma_mp_clear (&t2);
4381 T1:tma_mp_clear (&t1);
4382 X1:tma_mp_clear (&x1);
4383 X0:tma_mp_clear (&x0);
4393 /* End: bn_tma_mp_karatsuba_sqr.c */
4395 /* Start: bn_tma_mp_lcm.c */
4398 /* LibTomMath, multiple-precision integer library -- Tom St Denis
4400 * LibTomMath is a library that provides multiple-precision
4401 * integer arithmetic as well as number theoretic functionality.
4403 * The library was designed directly after the MPI library by
4404 * Michael Fromberger but has been written from scratch with
4405 * additional optimizations in place.
4407 * The library is free for all purposes without any express
4408 * guarantee it works.
4410 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
4413 /* computes least common multiple as |a*b|/(a, b) */
4414 int tma_mp_lcm (tma_mp_int * a, tma_mp_int * b, tma_mp_int * c)
4420 if ((res = tma_mp_init_multi (&t1, &t2, NULL)) != MP_OKAY) {
4424 /* t1 = get the GCD of the two inputs */
4425 if ((res = tma_mp_gcd (a, b, &t1)) != MP_OKAY) {
4429 /* divide the smallest by the GCD */
4430 if (tma_mp_cmp_mag(a, b) == MP_LT) {
4431 /* store quotient in t2 such that t2 * b is the LCM */
4432 if ((res = tma_mp_div(a, &t1, &t2, NULL)) != MP_OKAY) {
4435 res = tma_mp_mul(b, &t2, c);
4437 /* store quotient in t2 such that t2 * a is the LCM */
4438 if ((res = tma_mp_div(b, &t1, &t2, NULL)) != MP_OKAY) {
4441 res = tma_mp_mul(a, &t2, c);
4444 /* fix the sign to positive */
4448 tma_mp_clear_multi (&t1, &t2, NULL);
4457 /* End: bn_tma_mp_lcm.c */
4459 /* Start: bn_tma_mp_lshd.c */
4462 /* LibTomMath, multiple-precision integer library -- Tom St Denis
4464 * LibTomMath is a library that provides multiple-precision
4465 * integer arithmetic as well as number theoretic functionality.
4467 * The library was designed directly after the MPI library by
4468 * Michael Fromberger but has been written from scratch with
4469 * additional optimizations in place.
4471 * The library is free for all purposes without any express
4472 * guarantee it works.
4474 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
4477 /* shift left a certain amount of digits */
4478 int tma_mp_lshd (tma_mp_int * a, int b)
4482 /* if its less than zero return */
4487 /* grow to fit the new digits */
4488 if (a->alloc < a->used + b) {
4489 if ((res = tma_mp_grow (a, a->used + b)) != MP_OKAY) {
4495 register tma_mp_digit *top, *bottom;
4497 /* increment the used by the shift amount then copy upwards */
4501 top = a->dp + a->used - 1;
4504 bottom = a->dp + a->used - 1 - b;
4506 /* much like tma_mp_rshd this is implemented using a sliding window
4507 * except the window goes the otherway around. Copying from
4508 * the bottom to the top. see bn_tma_mp_rshd.c for more info.
4510 for (x = a->used - 1; x >= b; x--) {
4514 /* zero the lower digits */
4516 for (x = 0; x < b; x++) {
4528 /* End: bn_tma_mp_lshd.c */
4530 /* Start: bn_tma_mp_mod.c */
4533 /* LibTomMath, multiple-precision integer library -- Tom St Denis
4535 * LibTomMath is a library that provides multiple-precision
4536 * integer arithmetic as well as number theoretic functionality.
4538 * The library was designed directly after the MPI library by
4539 * Michael Fromberger but has been written from scratch with
4540 * additional optimizations in place.
4542 * The library is free for all purposes without any express
4543 * guarantee it works.
4545 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
4548 /* c = a mod b, 0 <= c < b */
4550 tma_mp_mod (tma_mp_int * a, tma_mp_int * b, tma_mp_int * c)
4555 if ((res = tma_mp_init (&t)) != MP_OKAY) {
4559 if ((res = tma_mp_div (a, b, NULL, &t)) != MP_OKAY) {
4564 if (t.sign != b->sign) {
4565 res = tma_mp_add (b, &t, c);
4568 tma_mp_exch (&t, c);
4580 /* End: bn_tma_mp_mod.c */
4582 /* Start: bn_tma_mp_mod_2d.c */
4584 #ifdef BN_MP_MOD_2D_C
4585 /* LibTomMath, multiple-precision integer library -- Tom St Denis
4587 * LibTomMath is a library that provides multiple-precision
4588 * integer arithmetic as well as number theoretic functionality.
4590 * The library was designed directly after the MPI library by
4591 * Michael Fromberger but has been written from scratch with
4592 * additional optimizations in place.
4594 * The library is free for all purposes without any express
4595 * guarantee it works.
4597 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
4600 /* calc a value mod 2**b */
4602 tma_mp_mod_2d (tma_mp_int * a, int b, tma_mp_int * c)
4606 /* if b is <= 0 then zero the int */
4612 /* if the modulus is larger than the value than return */
4613 if (b >= (int) (a->used * DIGIT_BIT)) {
4614 res = tma_mp_copy (a, c);
4619 if ((res = tma_mp_copy (a, c)) != MP_OKAY) {
4623 /* zero digits above the last digit of the modulus */
4624 for (x = (b / DIGIT_BIT) + ((b % DIGIT_BIT) == 0 ? 0 : 1); x < c->used; x++) {
4627 /* clear the digit that is not completely outside/inside the modulus */
4628 c->dp[b / DIGIT_BIT] &=
4629 (tma_mp_digit) ((((tma_mp_digit) 1) << (((tma_mp_digit) b) % DIGIT_BIT)) - ((tma_mp_digit) 1));
4639 /* End: bn_tma_mp_mod_2d.c */
4641 /* Start: bn_tma_mp_mod_d.c */
4643 #ifdef BN_MP_MOD_D_C
4644 /* LibTomMath, multiple-precision integer library -- Tom St Denis
4646 * LibTomMath is a library that provides multiple-precision
4647 * integer arithmetic as well as number theoretic functionality.
4649 * The library was designed directly after the MPI library by
4650 * Michael Fromberger but has been written from scratch with
4651 * additional optimizations in place.
4653 * The library is free for all purposes without any express
4654 * guarantee it works.
4656 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
4660 tma_mp_mod_d (tma_mp_int * a, tma_mp_digit b, tma_mp_digit * c)
4662 return tma_mp_div_d(a, b, NULL, c);
4670 /* End: bn_tma_mp_mod_d.c */
4672 /* Start: bn_tma_mp_montgomery_calc_normalization.c */
4674 #ifdef BN_MP_MONTGOMERY_CALC_NORMALIZATION_C
4675 /* LibTomMath, multiple-precision integer library -- Tom St Denis
4677 * LibTomMath is a library that provides multiple-precision
4678 * integer arithmetic as well as number theoretic functionality.
4680 * The library was designed directly after the MPI library by
4681 * Michael Fromberger but has been written from scratch with
4682 * additional optimizations in place.
4684 * The library is free for all purposes without any express
4685 * guarantee it works.
4687 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
4691 * shifts with subtractions when the result is greater than b.
4693 * The method is slightly modified to shift B unconditionally upto just under
4694 * the leading bit of b. This saves alot of multiple precision shifting.
4696 int tma_mp_montgomery_calc_normalization (tma_mp_int * a, tma_mp_int * b)
4700 /* how many bits of last digit does b use */
4701 bits = tma_mp_count_bits (b) % DIGIT_BIT;
4704 if ((res = tma_mp_2expt (a, (b->used - 1) * DIGIT_BIT + bits - 1)) != MP_OKAY) {
4713 /* now compute C = A * B mod b */
4714 for (x = bits - 1; x < (int)DIGIT_BIT; x++) {
4715 if ((res = tma_mp_mul_2 (a, a)) != MP_OKAY) {
4718 if (tma_mp_cmp_mag (a, b) != MP_LT) {
4719 if ((res = s_tma_mp_sub (a, b, a)) != MP_OKAY) {
4733 /* End: bn_tma_mp_montgomery_calc_normalization.c */
4735 /* Start: bn_tma_mp_montgomery_reduce.c */
4737 #ifdef BN_MP_MONTGOMERY_REDUCE_C
4738 /* LibTomMath, multiple-precision integer library -- Tom St Denis
4740 * LibTomMath is a library that provides multiple-precision
4741 * integer arithmetic as well as number theoretic functionality.
4743 * The library was designed directly after the MPI library by
4744 * Michael Fromberger but has been written from scratch with
4745 * additional optimizations in place.
4747 * The library is free for all purposes without any express
4748 * guarantee it works.
4750 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
4753 /* computes xR**-1 == x (mod N) via Montgomery Reduction */
4755 tma_mp_montgomery_reduce (tma_mp_int * x, tma_mp_int * n, tma_mp_digit rho)
4760 /* can the fast reduction [comba] method be used?
4762 * Note that unlike in mul you're safely allowed *less*
4763 * than the available columns [255 per default] since carries
4764 * are fixed up in the inner loop.
4766 digs = n->used * 2 + 1;
4767 if ((digs < MP_WARRAY) &&
4769 (1 << ((CHAR_BIT * sizeof (tma_mp_word)) - (2 * DIGIT_BIT)))) {
4770 return fast_tma_mp_montgomery_reduce (x, n, rho);
4773 /* grow the input as required */
4774 if (x->alloc < digs) {
4775 if ((res = tma_mp_grow (x, digs)) != MP_OKAY) {
4781 for (ix = 0; ix < n->used; ix++) {
4782 /* mu = ai * rho mod b
4784 * The value of rho must be precalculated via
4785 * montgomery_setup() such that
4786 * it equals -1/n0 mod b this allows the
4787 * following inner loop to reduce the
4788 * input one digit at a time
4790 mu = (tma_mp_digit) (((tma_mp_word)x->dp[ix]) * ((tma_mp_word)rho) & MP_MASK);
4792 /* a = a + mu * m * b**i */
4795 register tma_mp_digit *tmpn, *tmpx, u;
4796 register tma_mp_word r;
4798 /* alias for digits of the modulus */
4801 /* alias for the digits of x [the input] */
4804 /* set the carry to zero */
4807 /* Multiply and add in place */
4808 for (iy = 0; iy < n->used; iy++) {
4809 /* compute product and sum */
4810 r = ((tma_mp_word)mu) * ((tma_mp_word)*tmpn++) +
4811 ((tma_mp_word) u) + ((tma_mp_word) * tmpx);
4814 u = (tma_mp_digit)(r >> ((tma_mp_word) DIGIT_BIT));
4817 *tmpx++ = (tma_mp_digit)(r & ((tma_mp_word) MP_MASK));
4819 /* At this point the ix'th digit of x should be zero */
4822 /* propagate carries upwards as required*/
4825 u = *tmpx >> DIGIT_BIT;
4831 /* at this point the n.used'th least
4832 * significant digits of x are all zero
4833 * which means we can shift x to the
4834 * right by n.used digits and the
4835 * residue is unchanged.
4838 /* x = x/b**n.used */
4840 tma_mp_rshd (x, n->used);
4842 /* if x >= n then x = x - n */
4843 if (tma_mp_cmp_mag (x, n) != MP_LT) {
4844 return s_tma_mp_sub (x, n, x);
4855 /* End: bn_tma_mp_montgomery_reduce.c */
4857 /* Start: bn_tma_mp_montgomery_setup.c */
4859 #ifdef BN_MP_MONTGOMERY_SETUP_C
4860 /* LibTomMath, multiple-precision integer library -- Tom St Denis
4862 * LibTomMath is a library that provides multiple-precision
4863 * integer arithmetic as well as number theoretic functionality.
4865 * The library was designed directly after the MPI library by
4866 * Michael Fromberger but has been written from scratch with
4867 * additional optimizations in place.
4869 * The library is free for all purposes without any express
4870 * guarantee it works.
4872 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
4875 /* setups the montgomery reduction stuff */
4877 tma_mp_montgomery_setup (tma_mp_int * n, tma_mp_digit * rho)
4881 /* fast inversion mod 2**k
4883 * Based on the fact that
4885 * XA = 1 (mod 2**n) => (X(2-XA)) A = 1 (mod 2**2n)
4886 * => 2*X*A - X*X*A*A = 1
4887 * => 2*(1) - (1) = 1
4895 x = (((b + 2) & 4) << 1) + b; /* here x*a==1 mod 2**4 */
4896 x *= 2 - b * x; /* here x*a==1 mod 2**8 */
4897 #if !defined(MP_8BIT)
4898 x *= 2 - b * x; /* here x*a==1 mod 2**16 */
4900 #if defined(MP_64BIT) || !(defined(MP_8BIT) || defined(MP_16BIT))
4901 x *= 2 - b * x; /* here x*a==1 mod 2**32 */
4904 x *= 2 - b * x; /* here x*a==1 mod 2**64 */
4907 /* rho = -1/m mod b */
4908 *rho = (unsigned long)(((tma_mp_word)1 << ((tma_mp_word) DIGIT_BIT)) - x) & MP_MASK;
4918 /* End: bn_tma_mp_montgomery_setup.c */
4920 /* Start: bn_tma_mp_mul.c */
4923 /* LibTomMath, multiple-precision integer library -- Tom St Denis
4925 * LibTomMath is a library that provides multiple-precision
4926 * integer arithmetic as well as number theoretic functionality.
4928 * The library was designed directly after the MPI library by
4929 * Michael Fromberger but has been written from scratch with
4930 * additional optimizations in place.
4932 * The library is free for all purposes without any express
4933 * guarantee it works.
4935 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
4938 /* high level multiplication (handles sign) */
4939 int tma_mp_mul (tma_mp_int * a, tma_mp_int * b, tma_mp_int * c)
4942 neg = (a->sign == b->sign) ? MP_ZPOS : MP_NEG;
4944 /* use Toom-Cook? */
4945 #ifdef BN_MP_TOOM_MUL_C
4946 if (MIN (a->used, b->used) >= TOOM_MUL_CUTOFF) {
4947 res = tma_mp_toom_mul(a, b, c);
4950 #ifdef BN_MP_KARATSUBA_MUL_C
4951 /* use Karatsuba? */
4952 if (MIN (a->used, b->used) >= KARATSUBA_MUL_CUTOFF) {
4953 res = tma_mp_karatsuba_mul (a, b, c);
4957 /* can we use the fast multiplier?
4959 * The fast multiplier can be used if the output will
4960 * have less than MP_WARRAY digits and the number of
4961 * digits won't affect carry propagation
4963 int digs = a->used + b->used + 1;
4965 #ifdef BN_FAST_S_MP_MUL_DIGS_C
4966 if ((digs < MP_WARRAY) &&
4967 MIN(a->used, b->used) <=
4968 (1 << ((CHAR_BIT * sizeof (tma_mp_word)) - (2 * DIGIT_BIT)))) {
4969 res = fast_s_tma_mp_mul_digs (a, b, c, digs);
4972 #ifdef BN_S_MP_MUL_DIGS_C
4973 res = s_tma_mp_mul (a, b, c); /* uses s_tma_mp_mul_digs */
4979 c->sign = (c->used > 0) ? neg : MP_ZPOS;
4988 /* End: bn_tma_mp_mul.c */
4990 /* Start: bn_tma_mp_mul_2.c */
4992 #ifdef BN_MP_MUL_2_C
4993 /* LibTomMath, multiple-precision integer library -- Tom St Denis
4995 * LibTomMath is a library that provides multiple-precision
4996 * integer arithmetic as well as number theoretic functionality.
4998 * The library was designed directly after the MPI library by
4999 * Michael Fromberger but has been written from scratch with
5000 * additional optimizations in place.
5002 * The library is free for all purposes without any express
5003 * guarantee it works.
5005 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
5009 int tma_mp_mul_2(tma_mp_int * a, tma_mp_int * b)
5011 int x, res, oldused;
5013 /* grow to accomodate result */
5014 if (b->alloc < a->used + 1) {
5015 if ((res = tma_mp_grow (b, a->used + 1)) != MP_OKAY) {
5024 register tma_mp_digit r, rr, *tmpa, *tmpb;
5026 /* alias for source */
5029 /* alias for dest */
5034 for (x = 0; x < a->used; x++) {
5036 /* get what will be the *next* carry bit from the
5037 * MSB of the current digit
5039 rr = *tmpa >> ((tma_mp_digit)(DIGIT_BIT - 1));
5041 /* now shift up this digit, add in the carry [from the previous] */
5042 *tmpb++ = ((*tmpa++ << ((tma_mp_digit)1)) | r) & MP_MASK;
5044 /* copy the carry that would be from the source
5045 * digit into the next iteration
5050 /* new leading digit? */
5052 /* add a MSB which is always 1 at this point */
5057 /* now zero any excess digits on the destination
5058 * that we didn't write to
5060 tmpb = b->dp + b->used;
5061 for (x = b->used; x < oldused; x++) {
5074 /* End: bn_tma_mp_mul_2.c */
5076 /* Start: bn_tma_mp_mul_2d.c */
5078 #ifdef BN_MP_MUL_2D_C
5079 /* LibTomMath, multiple-precision integer library -- Tom St Denis
5081 * LibTomMath is a library that provides multiple-precision
5082 * integer arithmetic as well as number theoretic functionality.
5084 * The library was designed directly after the MPI library by
5085 * Michael Fromberger but has been written from scratch with
5086 * additional optimizations in place.
5088 * The library is free for all purposes without any express
5089 * guarantee it works.
5091 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
5094 /* shift left by a certain bit count */
5095 int tma_mp_mul_2d (tma_mp_int * a, int b, tma_mp_int * c)
5102 if ((res = tma_mp_copy (a, c)) != MP_OKAY) {
5107 if (c->alloc < (int)(c->used + b/DIGIT_BIT + 1)) {
5108 if ((res = tma_mp_grow (c, c->used + b / DIGIT_BIT + 1)) != MP_OKAY) {
5113 /* shift by as many digits in the bit count */
5114 if (b >= (int)DIGIT_BIT) {
5115 if ((res = tma_mp_lshd (c, b / DIGIT_BIT)) != MP_OKAY) {
5120 /* shift any bit count < DIGIT_BIT */
5121 d = (tma_mp_digit) (b % DIGIT_BIT);
5123 register tma_mp_digit *tmpc, shift, mask, r, rr;
5126 /* bitmask for carries */
5127 mask = (((tma_mp_digit)1) << d) - 1;
5129 /* shift for msbs */
5130 shift = DIGIT_BIT - d;
5137 for (x = 0; x < c->used; x++) {
5138 /* get the higher bits of the current word */
5139 rr = (*tmpc >> shift) & mask;
5141 /* shift the current word and OR in the carry */
5142 *tmpc = ((*tmpc << d) | r) & MP_MASK;
5145 /* set the carry to the carry bits of the current word */
5149 /* set final carry */
5151 c->dp[(c->used)++] = r;
5163 /* End: bn_tma_mp_mul_2d.c */
5165 /* Start: bn_tma_mp_mul_d.c */
5167 #ifdef BN_MP_MUL_D_C
5168 /* LibTomMath, multiple-precision integer library -- Tom St Denis
5170 * LibTomMath is a library that provides multiple-precision
5171 * integer arithmetic as well as number theoretic functionality.
5173 * The library was designed directly after the MPI library by
5174 * Michael Fromberger but has been written from scratch with
5175 * additional optimizations in place.
5177 * The library is free for all purposes without any express
5178 * guarantee it works.
5180 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
5183 /* multiply by a digit */
5185 tma_mp_mul_d (tma_mp_int * a, tma_mp_digit b, tma_mp_int * c)
5187 tma_mp_digit u, *tmpa, *tmpc;
5189 int ix, res, olduse;
5191 /* make sure c is big enough to hold a*b */
5192 if (c->alloc < a->used + 1) {
5193 if ((res = tma_mp_grow (c, a->used + 1)) != MP_OKAY) {
5198 /* get the original destinations used count */
5204 /* alias for a->dp [source] */
5207 /* alias for c->dp [dest] */
5213 /* compute columns */
5214 for (ix = 0; ix < a->used; ix++) {
5215 /* compute product and carry sum for this term */
5216 r = ((tma_mp_word) u) + ((tma_mp_word)*tmpa++) * ((tma_mp_word)b);
5218 /* mask off higher bits to get a single digit */
5219 *tmpc++ = (tma_mp_digit) (r & ((tma_mp_word) MP_MASK));
5221 /* send carry into next iteration */
5222 u = (tma_mp_digit) (r >> ((tma_mp_word) DIGIT_BIT));
5225 /* store final carry [if any] and increment ix offset */
5229 /* now zero digits above the top */
5230 while (ix++ < olduse) {
5234 /* set used count */
5235 c->used = a->used + 1;
5246 /* End: bn_tma_mp_mul_d.c */
5248 /* Start: bn_tma_mp_mulmod.c */
5250 #ifdef BN_MP_MULMOD_C
5251 /* LibTomMath, multiple-precision integer library -- Tom St Denis
5253 * LibTomMath is a library that provides multiple-precision
5254 * integer arithmetic as well as number theoretic functionality.
5256 * The library was designed directly after the MPI library by
5257 * Michael Fromberger but has been written from scratch with
5258 * additional optimizations in place.
5260 * The library is free for all purposes without any express
5261 * guarantee it works.
5263 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
5266 /* d = a * b (mod c) */
5267 int tma_mp_mulmod (tma_mp_int * a, tma_mp_int * b, tma_mp_int * c, tma_mp_int * d)
5272 if ((res = tma_mp_init (&t)) != MP_OKAY) {
5276 if ((res = tma_mp_mul (a, b, &t)) != MP_OKAY) {
5280 res = tma_mp_mod (&t, c, d);
5290 /* End: bn_tma_mp_mulmod.c */
5292 /* Start: bn_tma_mp_n_root.c */
5294 #ifdef BN_MP_N_ROOT_C
5295 /* LibTomMath, multiple-precision integer library -- Tom St Denis
5297 * LibTomMath is a library that provides multiple-precision
5298 * integer arithmetic as well as number theoretic functionality.
5300 * The library was designed directly after the MPI library by
5301 * Michael Fromberger but has been written from scratch with
5302 * additional optimizations in place.
5304 * The library is free for all purposes without any express
5305 * guarantee it works.
5307 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
5310 /* find the n'th root of an integer
5312 * Result found such that (c)**b <= a and (c+1)**b > a
5314 * This algorithm uses Newton's approximation
5315 * x[i+1] = x[i] - f(x[i])/f'(x[i])
5316 * which will find the root in log(N) time where
5317 * each step involves a fair bit. This is not meant to
5318 * find huge roots [square and cube, etc].
5320 int tma_mp_n_root (tma_mp_int * a, tma_mp_digit b, tma_mp_int * c)
5322 tma_mp_int t1, t2, t3;
5325 /* input must be positive if b is even */
5326 if ((b & 1) == 0 && a->sign == MP_NEG) {
5330 if ((res = tma_mp_init (&t1)) != MP_OKAY) {
5334 if ((res = tma_mp_init (&t2)) != MP_OKAY) {
5338 if ((res = tma_mp_init (&t3)) != MP_OKAY) {
5342 /* if a is negative fudge the sign but keep track */
5347 tma_mp_set (&t2, 2);
5351 if ((res = tma_mp_copy (&t2, &t1)) != MP_OKAY) {
5355 /* t2 = t1 - ((t1**b - a) / (b * t1**(b-1))) */
5357 /* t3 = t1**(b-1) */
5358 if ((res = tma_mp_expt_d (&t1, b - 1, &t3)) != MP_OKAY) {
5364 if ((res = tma_mp_mul (&t3, &t1, &t2)) != MP_OKAY) {
5368 /* t2 = t1**b - a */
5369 if ((res = tma_mp_sub (&t2, a, &t2)) != MP_OKAY) {
5374 /* t3 = t1**(b-1) * b */
5375 if ((res = tma_mp_mul_d (&t3, b, &t3)) != MP_OKAY) {
5379 /* t3 = (t1**b - a)/(b * t1**(b-1)) */
5380 if ((res = tma_mp_div (&t2, &t3, &t3, NULL)) != MP_OKAY) {
5384 if ((res = tma_mp_sub (&t1, &t3, &t2)) != MP_OKAY) {
5387 } while (tma_mp_cmp (&t1, &t2) != MP_EQ);
5389 /* result can be off by a few so check */
5391 if ((res = tma_mp_expt_d (&t1, b, &t2)) != MP_OKAY) {
5395 if (tma_mp_cmp (&t2, a) == MP_GT) {
5396 if ((res = tma_mp_sub_d (&t1, 1, &t1)) != MP_OKAY) {
5404 /* reset the sign of a first */
5407 /* set the result */
5408 tma_mp_exch (&t1, c);
5410 /* set the sign of the result */
5415 LBL_T3:tma_mp_clear (&t3);
5416 LBL_T2:tma_mp_clear (&t2);
5417 LBL_T1:tma_mp_clear (&t1);
5426 /* End: bn_tma_mp_n_root.c */
5428 /* Start: bn_tma_mp_neg.c */
5431 /* LibTomMath, multiple-precision integer library -- Tom St Denis
5433 * LibTomMath is a library that provides multiple-precision
5434 * integer arithmetic as well as number theoretic functionality.
5436 * The library was designed directly after the MPI library by
5437 * Michael Fromberger but has been written from scratch with
5438 * additional optimizations in place.
5440 * The library is free for all purposes without any express
5441 * guarantee it works.
5443 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
5447 int tma_mp_neg (tma_mp_int * a, tma_mp_int * b)
5451 if ((res = tma_mp_copy (a, b)) != MP_OKAY) {
5456 if (tma_mp_iszero(b) != MP_YES) {
5457 b->sign = (a->sign == MP_ZPOS) ? MP_NEG : MP_ZPOS;
5470 /* End: bn_tma_mp_neg.c */
5472 /* Start: bn_tma_mp_or.c */
5475 /* LibTomMath, multiple-precision integer library -- Tom St Denis
5477 * LibTomMath is a library that provides multiple-precision
5478 * integer arithmetic as well as number theoretic functionality.
5480 * The library was designed directly after the MPI library by
5481 * Michael Fromberger but has been written from scratch with
5482 * additional optimizations in place.
5484 * The library is free for all purposes without any express
5485 * guarantee it works.
5487 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
5490 /* OR two ints together */
5491 int tma_mp_or (tma_mp_int * a, tma_mp_int * b, tma_mp_int * c)
5496 if (a->used > b->used) {
5497 if ((res = tma_mp_init_copy (&t, a)) != MP_OKAY) {
5503 if ((res = tma_mp_init_copy (&t, b)) != MP_OKAY) {
5510 for (ix = 0; ix < px; ix++) {
5511 t.dp[ix] |= x->dp[ix];
5514 tma_mp_exch (c, &t);
5524 /* End: bn_tma_mp_or.c */
5526 /* Start: bn_tma_mp_prime_fermat.c */
5528 #ifdef BN_MP_PRIME_FERMAT_C
5529 /* LibTomMath, multiple-precision integer library -- Tom St Denis
5531 * LibTomMath is a library that provides multiple-precision
5532 * integer arithmetic as well as number theoretic functionality.
5534 * The library was designed directly after the MPI library by
5535 * Michael Fromberger but has been written from scratch with
5536 * additional optimizations in place.
5538 * The library is free for all purposes without any express
5539 * guarantee it works.
5541 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
5544 /* performs one Fermat test.
5546 * If "a" were prime then b**a == b (mod a) since the order of
5547 * the multiplicative sub-group would be phi(a) = a-1. That means
5548 * it would be the same as b**(a mod (a-1)) == b**1 == b (mod a).
5550 * Sets result to 1 if the congruence holds, or zero otherwise.
5552 int tma_mp_prime_fermat (tma_mp_int * a, tma_mp_int * b, int *result)
5557 /* default to composite */
5561 if (tma_mp_cmp_d(b, 1) != MP_GT) {
5566 if ((err = tma_mp_init (&t)) != MP_OKAY) {
5570 /* compute t = b**a mod a */
5571 if ((err = tma_mp_exptmod (b, a, a, &t)) != MP_OKAY) {
5575 /* is it equal to b? */
5576 if (tma_mp_cmp (&t, b) == MP_EQ) {
5581 LBL_T:tma_mp_clear (&t);
5590 /* End: bn_tma_mp_prime_fermat.c */
5592 /* Start: bn_tma_mp_prime_is_divisible.c */
5594 #ifdef BN_MP_PRIME_IS_DIVISIBLE_C
5595 /* LibTomMath, multiple-precision integer library -- Tom St Denis
5597 * LibTomMath is a library that provides multiple-precision
5598 * integer arithmetic as well as number theoretic functionality.
5600 * The library was designed directly after the MPI library by
5601 * Michael Fromberger but has been written from scratch with
5602 * additional optimizations in place.
5604 * The library is free for all purposes without any express
5605 * guarantee it works.
5607 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
5610 /* determines if an integers is divisible by one
5611 * of the first PRIME_SIZE primes or not
5613 * sets result to 0 if not, 1 if yes
5615 int tma_mp_prime_is_divisible (tma_mp_int * a, int *result)
5620 /* default to not */
5623 for (ix = 0; ix < PRIME_SIZE; ix++) {
5624 /* what is a mod LBL_prime_tab[ix] */
5625 if ((err = tma_mp_mod_d (a, ltm_prime_tab[ix], &res)) != MP_OKAY) {
5629 /* is the residue zero? */
5644 /* End: bn_tma_mp_prime_is_divisible.c */
5646 /* Start: bn_tma_mp_prime_is_prime.c */
5648 #ifdef BN_MP_PRIME_IS_PRIME_C
5649 /* LibTomMath, multiple-precision integer library -- Tom St Denis
5651 * LibTomMath is a library that provides multiple-precision
5652 * integer arithmetic as well as number theoretic functionality.
5654 * The library was designed directly after the MPI library by
5655 * Michael Fromberger but has been written from scratch with
5656 * additional optimizations in place.
5658 * The library is free for all purposes without any express
5659 * guarantee it works.
5661 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
5664 /* performs a variable number of rounds of Miller-Rabin
5666 * Probability of error after t rounds is no more than
5669 * Sets result to 1 if probably prime, 0 otherwise
5671 int tma_mp_prime_is_prime (tma_mp_int * a, int t, int *result)
5679 /* valid value of t? */
5680 if (t <= 0 || t > PRIME_SIZE) {
5684 /* is the input equal to one of the primes in the table? */
5685 for (ix = 0; ix < PRIME_SIZE; ix++) {
5686 if (tma_mp_cmp_d(a, ltm_prime_tab[ix]) == MP_EQ) {
5692 /* first perform trial division */
5693 if ((err = tma_mp_prime_is_divisible (a, &res)) != MP_OKAY) {
5697 /* return if it was trivially divisible */
5698 if (res == MP_YES) {
5702 /* now perform the miller-rabin rounds */
5703 if ((err = tma_mp_init (&b)) != MP_OKAY) {
5707 for (ix = 0; ix < t; ix++) {
5709 tma_mp_set (&b, ltm_prime_tab[ix]);
5711 if ((err = tma_mp_prime_miller_rabin (a, &b, &res)) != MP_OKAY) {
5720 /* passed the test */
5722 LBL_B:tma_mp_clear (&b);
5731 /* End: bn_tma_mp_prime_is_prime.c */
5733 /* Start: bn_tma_mp_prime_miller_rabin.c */
5735 #ifdef BN_MP_PRIME_MILLER_RABIN_C
5736 /* LibTomMath, multiple-precision integer library -- Tom St Denis
5738 * LibTomMath is a library that provides multiple-precision
5739 * integer arithmetic as well as number theoretic functionality.
5741 * The library was designed directly after the MPI library by
5742 * Michael Fromberger but has been written from scratch with
5743 * additional optimizations in place.
5745 * The library is free for all purposes without any express
5746 * guarantee it works.
5748 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
5751 /* Miller-Rabin test of "a" to the base of "b" as described in
5752 * HAC pp. 139 Algorithm 4.24
5754 * Sets result to 0 if definitely composite or 1 if probably prime.
5755 * Randomly the chance of error is no more than 1/4 and often
5758 int tma_mp_prime_miller_rabin (tma_mp_int * a, tma_mp_int * b, int *result)
5760 tma_mp_int n1, y, r;
5767 if (tma_mp_cmp_d(b, 1) != MP_GT) {
5771 /* get n1 = a - 1 */
5772 if ((err = tma_mp_init_copy (&n1, a)) != MP_OKAY) {
5775 if ((err = tma_mp_sub_d (&n1, 1, &n1)) != MP_OKAY) {
5779 /* set 2**s * r = n1 */
5780 if ((err = tma_mp_init_copy (&r, &n1)) != MP_OKAY) {
5784 /* count the number of least significant bits
5787 s = tma_mp_cnt_lsb(&r);
5789 /* now divide n - 1 by 2**s */
5790 if ((err = tma_mp_div_2d (&r, s, &r, NULL)) != MP_OKAY) {
5794 /* compute y = b**r mod a */
5795 if ((err = tma_mp_init (&y)) != MP_OKAY) {
5798 if ((err = tma_mp_exptmod (b, &r, a, &y)) != MP_OKAY) {
5802 /* if y != 1 and y != n1 do */
5803 if (tma_mp_cmp_d (&y, 1) != MP_EQ && tma_mp_cmp (&y, &n1) != MP_EQ) {
5805 /* while j <= s-1 and y != n1 */
5806 while ((j <= (s - 1)) && tma_mp_cmp (&y, &n1) != MP_EQ) {
5807 if ((err = tma_mp_sqrmod (&y, a, &y)) != MP_OKAY) {
5811 /* if y == 1 then composite */
5812 if (tma_mp_cmp_d (&y, 1) == MP_EQ) {
5819 /* if y != n1 then composite */
5820 if (tma_mp_cmp (&y, &n1) != MP_EQ) {
5825 /* probably prime now */
5827 LBL_Y:tma_mp_clear (&y);
5828 LBL_R:tma_mp_clear (&r);
5829 LBL_N1:tma_mp_clear (&n1);
5838 /* End: bn_tma_mp_prime_miller_rabin.c */
5840 /* Start: bn_tma_mp_prime_next_prime.c */
5842 #ifdef BN_MP_PRIME_NEXT_PRIME_C
5843 /* LibTomMath, multiple-precision integer library -- Tom St Denis
5845 * LibTomMath is a library that provides multiple-precision
5846 * integer arithmetic as well as number theoretic functionality.
5848 * The library was designed directly after the MPI library by
5849 * Michael Fromberger but has been written from scratch with
5850 * additional optimizations in place.
5852 * The library is free for all purposes without any express
5853 * guarantee it works.
5855 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
5858 /* finds the next prime after the number "a" using "t" trials
5861 * bbs_style = 1 means the prime must be congruent to 3 mod 4
5863 int tma_mp_prime_next_prime(tma_mp_int *a, int t, int bbs_style)
5866 tma_mp_digit res_tab[PRIME_SIZE], step, kstep;
5869 /* ensure t is valid */
5870 if (t <= 0 || t > PRIME_SIZE) {
5874 /* force positive */
5877 /* simple algo if a is less than the largest prime in the table */
5878 if (tma_mp_cmp_d(a, ltm_prime_tab[PRIME_SIZE-1]) == MP_LT) {
5879 /* find which prime it is bigger than */
5880 for (x = PRIME_SIZE - 2; x >= 0; x--) {
5881 if (tma_mp_cmp_d(a, ltm_prime_tab[x]) != MP_LT) {
5882 if (bbs_style == 1) {
5883 /* ok we found a prime smaller or
5884 * equal [so the next is larger]
5886 * however, the prime must be
5887 * congruent to 3 mod 4
5889 if ((ltm_prime_tab[x + 1] & 3) != 3) {
5890 /* scan upwards for a prime congruent to 3 mod 4 */
5891 for (y = x + 1; y < PRIME_SIZE; y++) {
5892 if ((ltm_prime_tab[y] & 3) == 3) {
5893 tma_mp_set(a, ltm_prime_tab[y]);
5899 tma_mp_set(a, ltm_prime_tab[x + 1]);
5904 /* at this point a maybe 1 */
5905 if (tma_mp_cmp_d(a, 1) == MP_EQ) {
5909 /* fall through to the sieve */
5912 /* generate a prime congruent to 3 mod 4 or 1/3 mod 4? */
5913 if (bbs_style == 1) {
5919 /* at this point we will use a combination of a sieve and Miller-Rabin */
5921 if (bbs_style == 1) {
5922 /* if a mod 4 != 3 subtract the correct value to make it so */
5923 if ((a->dp[0] & 3) != 3) {
5924 if ((err = tma_mp_sub_d(a, (a->dp[0] & 3) + 1, a)) != MP_OKAY) { return err; };
5927 if (tma_mp_iseven(a) == 1) {
5929 if ((err = tma_mp_sub_d(a, 1, a)) != MP_OKAY) {
5935 /* generate the restable */
5936 for (x = 1; x < PRIME_SIZE; x++) {
5937 if ((err = tma_mp_mod_d(a, ltm_prime_tab[x], res_tab + x)) != MP_OKAY) {
5942 /* init temp used for Miller-Rabin Testing */
5943 if ((err = tma_mp_init(&b)) != MP_OKAY) {
5948 /* skip to the next non-trivially divisible candidate */
5951 /* y == 1 if any residue was zero [e.g. cannot be prime] */
5954 /* increase step to next candidate */
5957 /* compute the new residue without using division */
5958 for (x = 1; x < PRIME_SIZE; x++) {
5959 /* add the step to each residue */
5960 res_tab[x] += kstep;
5962 /* subtract the modulus [instead of using division] */
5963 if (res_tab[x] >= ltm_prime_tab[x]) {
5964 res_tab[x] -= ltm_prime_tab[x];
5967 /* set flag if zero */
5968 if (res_tab[x] == 0) {
5972 } while (y == 1 && step < ((((tma_mp_digit)1)<<DIGIT_BIT) - kstep));
5975 if ((err = tma_mp_add_d(a, step, a)) != MP_OKAY) {
5979 /* if didn't pass sieve and step == MAX then skip test */
5980 if (y == 1 && step >= ((((tma_mp_digit)1)<<DIGIT_BIT) - kstep)) {
5984 /* is this prime? */
5985 for (x = 0; x < t; x++) {
5986 tma_mp_set(&b, ltm_prime_tab[t]);
5987 if ((err = tma_mp_prime_miller_rabin(a, &b, &res)) != MP_OKAY) {
5995 if (res == MP_YES) {
6012 /* End: bn_tma_mp_prime_next_prime.c */
6014 /* Start: bn_tma_mp_prime_rabin_miller_trials.c */
6016 #ifdef BN_MP_PRIME_RABIN_MILLER_TRIALS_C
6017 /* LibTomMath, multiple-precision integer library -- Tom St Denis
6019 * LibTomMath is a library that provides multiple-precision
6020 * integer arithmetic as well as number theoretic functionality.
6022 * The library was designed directly after the MPI library by
6023 * Michael Fromberger but has been written from scratch with
6024 * additional optimizations in place.
6026 * The library is free for all purposes without any express
6027 * guarantee it works.
6029 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
6033 static const struct {
6046 /* returns # of RM trials required for a given bit size */
6047 int tma_mp_prime_rabin_miller_trials(int size)
6051 for (x = 0; x < (int)(sizeof(sizes)/(sizeof(sizes[0]))); x++) {
6052 if (sizes[x].k == size) {
6054 } else if (sizes[x].k > size) {
6055 return (x == 0) ? sizes[0].t : sizes[x - 1].t;
6058 return sizes[x-1].t + 1;
6068 /* End: bn_tma_mp_prime_rabin_miller_trials.c */
6070 /* Start: bn_tma_mp_prime_random_ex.c */
6072 #ifdef BN_MP_PRIME_RANDOM_EX_C
6073 /* LibTomMath, multiple-precision integer library -- Tom St Denis
6075 * LibTomMath is a library that provides multiple-precision
6076 * integer arithmetic as well as number theoretic functionality.
6078 * The library was designed directly after the MPI library by
6079 * Michael Fromberger but has been written from scratch with
6080 * additional optimizations in place.
6082 * The library is free for all purposes without any express
6083 * guarantee it works.
6085 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
6088 /* makes a truly random prime of a given size (bits),
6090 * Flags are as follows:
6092 * LTM_PRIME_BBS - make prime congruent to 3 mod 4
6093 * LTM_PRIME_SAFE - make sure (p-1)/2 is prime as well (implies LTM_PRIME_BBS)
6094 * LTM_PRIME_2MSB_OFF - make the 2nd highest bit zero
6095 * LTM_PRIME_2MSB_ON - make the 2nd highest bit one
6097 * You have to supply a callback which fills in a buffer with random bytes. "dat" is a parameter you can
6098 * have passed to the callback (e.g. a state or something). This function doesn't use "dat" itself
6103 /* This is possibly the mother of all prime generation functions, muahahahahaha! */
6104 int tma_mp_prime_random_ex(tma_mp_int *a, int t, int size, int flags, ltm_prime_callback cb, void *dat)
6106 unsigned char *tmp, maskAND, maskOR_msb, maskOR_lsb;
6107 int res, err, bsize, maskOR_msb_offset;
6109 /* sanity check the input */
6110 if (size <= 1 || t <= 0) {
6114 /* LTM_PRIME_SAFE implies LTM_PRIME_BBS */
6115 if (flags & LTM_PRIME_SAFE) {
6116 flags |= LTM_PRIME_BBS;
6119 /* calc the byte size */
6120 bsize = (size>>3) + ((size&7)?1:0);
6122 /* we need a buffer of bsize bytes */
6123 tmp = OPT_CAST(unsigned char) XMALLOC(bsize);
6128 /* calc the maskAND value for the MSbyte*/
6129 maskAND = ((size&7) == 0) ? 0xFF : (0xFF >> (8 - (size & 7)));
6131 /* calc the maskOR_msb */
6133 maskOR_msb_offset = ((size & 7) == 1) ? 1 : 0;
6134 if (flags & LTM_PRIME_2MSB_ON) {
6135 maskOR_msb |= 0x80 >> ((9 - size) & 7);
6138 /* get the maskOR_lsb */
6140 if (flags & LTM_PRIME_BBS) {
6145 /* read the bytes */
6146 if (cb(tmp, bsize, dat) != bsize) {
6151 /* work over the MSbyte */
6153 tmp[0] |= 1 << ((size - 1) & 7);
6155 /* mix in the maskORs */
6156 tmp[maskOR_msb_offset] |= maskOR_msb;
6157 tmp[bsize-1] |= maskOR_lsb;
6160 if ((err = tma_mp_read_unsigned_bin(a, tmp, bsize)) != MP_OKAY) { goto error; }
6163 if ((err = tma_mp_prime_is_prime(a, t, &res)) != MP_OKAY) { goto error; }
6168 if (flags & LTM_PRIME_SAFE) {
6169 /* see if (a-1)/2 is prime */
6170 if ((err = tma_mp_sub_d(a, 1, a)) != MP_OKAY) { goto error; }
6171 if ((err = tma_mp_div_2(a, a)) != MP_OKAY) { goto error; }
6174 if ((err = tma_mp_prime_is_prime(a, t, &res)) != MP_OKAY) { goto error; }
6176 } while (res == MP_NO);
6178 if (flags & LTM_PRIME_SAFE) {
6179 /* restore a to the original value */
6180 if ((err = tma_mp_mul_2(a, a)) != MP_OKAY) { goto error; }
6181 if ((err = tma_mp_add_d(a, 1, a)) != MP_OKAY) { goto error; }
6197 /* End: bn_tma_mp_prime_random_ex.c */
6199 /* Start: bn_tma_mp_radix_size.c */
6201 #ifdef BN_MP_RADIX_SIZE_C
6202 /* LibTomMath, multiple-precision integer library -- Tom St Denis
6204 * LibTomMath is a library that provides multiple-precision
6205 * integer arithmetic as well as number theoretic functionality.
6207 * The library was designed directly after the MPI library by
6208 * Michael Fromberger but has been written from scratch with
6209 * additional optimizations in place.
6211 * The library is free for all purposes without any express
6212 * guarantee it works.
6214 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
6217 /* returns size of ASCII reprensentation */
6218 int tma_mp_radix_size (tma_mp_int * a, int radix, int *size)
6226 /* special case for binary */
6228 *size = tma_mp_count_bits (a) + (a->sign == MP_NEG ? 1 : 0) + 1;
6232 /* make sure the radix is in range */
6233 if (radix < 2 || radix > 64) {
6237 if (tma_mp_iszero(a) == MP_YES) {
6242 /* digs is the digit count */
6245 /* if it's negative add one for the sign */
6246 if (a->sign == MP_NEG) {
6250 /* init a copy of the input */
6251 if ((res = tma_mp_init_copy (&t, a)) != MP_OKAY) {
6255 /* force temp to positive */
6258 /* fetch out all of the digits */
6259 while (tma_mp_iszero (&t) == MP_NO) {
6260 if ((res = tma_mp_div_d (&t, (tma_mp_digit) radix, &t, &d)) != MP_OKAY) {
6268 /* return digs + 1, the 1 is for the NULL byte that would be required. */
6279 /* End: bn_tma_mp_radix_size.c */
6281 /* Start: bn_tma_mp_radix_smap.c */
6283 #ifdef BN_MP_RADIX_SMAP_C
6284 /* LibTomMath, multiple-precision integer library -- Tom St Denis
6286 * LibTomMath is a library that provides multiple-precision
6287 * integer arithmetic as well as number theoretic functionality.
6289 * The library was designed directly after the MPI library by
6290 * Michael Fromberger but has been written from scratch with
6291 * additional optimizations in place.
6293 * The library is free for all purposes without any express
6294 * guarantee it works.
6296 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
6299 /* chars used in radix conversions */
6300 const char *tma_mp_s_rmap = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz+/";
6307 /* End: bn_tma_mp_radix_smap.c */
6309 /* Start: bn_tma_mp_rand.c */
6312 /* LibTomMath, multiple-precision integer library -- Tom St Denis
6314 * LibTomMath is a library that provides multiple-precision
6315 * integer arithmetic as well as number theoretic functionality.
6317 * The library was designed directly after the MPI library by
6318 * Michael Fromberger but has been written from scratch with
6319 * additional optimizations in place.
6321 * The library is free for all purposes without any express
6322 * guarantee it works.
6324 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
6327 /* makes a pseudo-random int of a given size */
6329 tma_mp_rand (tma_mp_int * a, int digits)
6339 /* first place a random non-zero digit */
6341 d = ((tma_mp_digit) abs (rand ())) & MP_MASK;
6344 if ((res = tma_mp_add_d (a, d, a)) != MP_OKAY) {
6348 while (--digits > 0) {
6349 if ((res = tma_mp_lshd (a, 1)) != MP_OKAY) {
6353 if ((res = tma_mp_add_d (a, ((tma_mp_digit) abs (rand ())), a)) != MP_OKAY) {
6366 /* End: bn_tma_mp_rand.c */
6368 /* Start: bn_tma_mp_read_radix.c */
6370 #ifdef BN_MP_READ_RADIX_C
6371 /* LibTomMath, multiple-precision integer library -- Tom St Denis
6373 * LibTomMath is a library that provides multiple-precision
6374 * integer arithmetic as well as number theoretic functionality.
6376 * The library was designed directly after the MPI library by
6377 * Michael Fromberger but has been written from scratch with
6378 * additional optimizations in place.
6380 * The library is free for all purposes without any express
6381 * guarantee it works.
6383 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
6386 /* read a string [ASCII] in a given radix */
6387 int tma_mp_read_radix (tma_mp_int * a, const char *str, int radix)
6392 /* zero the digit bignum */
6395 /* make sure the radix is ok */
6396 if (radix < 2 || radix > 64) {
6400 /* if the leading digit is a
6401 * minus set the sign to negative.
6410 /* set the integer to the default of zero */
6413 /* process each digit of the string */
6415 /* if the radix < 36 the conversion is case insensitive
6416 * this allows numbers like 1AB and 1ab to represent the same value
6419 ch = (char) ((radix < 36) ? toupper (*str) : *str);
6420 for (y = 0; y < 64; y++) {
6421 if (ch == tma_mp_s_rmap[y]) {
6426 /* if the char was found in the map
6427 * and is less than the given radix add it
6428 * to the number, otherwise exit the loop.
6431 if ((res = tma_mp_mul_d (a, (tma_mp_digit) radix, a)) != MP_OKAY) {
6434 if ((res = tma_mp_add_d (a, (tma_mp_digit) y, a)) != MP_OKAY) {
6443 /* set the sign only if a != 0 */
6444 if (tma_mp_iszero(a) != 1) {
6455 /* End: bn_tma_mp_read_radix.c */
6457 /* Start: bn_tma_mp_read_signed_bin.c */
6459 #ifdef BN_MP_READ_SIGNED_BIN_C
6460 /* LibTomMath, multiple-precision integer library -- Tom St Denis
6462 * LibTomMath is a library that provides multiple-precision
6463 * integer arithmetic as well as number theoretic functionality.
6465 * The library was designed directly after the MPI library by
6466 * Michael Fromberger but has been written from scratch with
6467 * additional optimizations in place.
6469 * The library is free for all purposes without any express
6470 * guarantee it works.
6472 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
6475 /* read signed bin, big endian, first byte is 0==positive or 1==negative */
6476 int tma_mp_read_signed_bin (tma_mp_int * a, const unsigned char *b, int c)
6480 /* read magnitude */
6481 if ((res = tma_mp_read_unsigned_bin (a, b + 1, c - 1)) != MP_OKAY) {
6485 /* first byte is 0 for positive, non-zero for negative */
6500 /* End: bn_tma_mp_read_signed_bin.c */
6502 /* Start: bn_tma_mp_read_unsigned_bin.c */
6504 #ifdef BN_MP_READ_UNSIGNED_BIN_C
6505 /* LibTomMath, multiple-precision integer library -- Tom St Denis
6507 * LibTomMath is a library that provides multiple-precision
6508 * integer arithmetic as well as number theoretic functionality.
6510 * The library was designed directly after the MPI library by
6511 * Michael Fromberger but has been written from scratch with
6512 * additional optimizations in place.
6514 * The library is free for all purposes without any express
6515 * guarantee it works.
6517 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
6520 /* reads a unsigned char array, assumes the msb is stored first [big endian] */
6521 int tma_mp_read_unsigned_bin (tma_mp_int * a, const unsigned char *b, int c)
6525 /* make sure there are at least two digits */
6527 if ((res = tma_mp_grow(a, 2)) != MP_OKAY) {
6535 /* read the bytes in */
6537 if ((res = tma_mp_mul_2d (a, 8, a)) != MP_OKAY) {
6545 a->dp[0] = (*b & MP_MASK);
6546 a->dp[1] |= ((*b++ >> 7U) & 1);
6559 /* End: bn_tma_mp_read_unsigned_bin.c */
6561 /* Start: bn_tma_mp_reduce.c */
6563 #ifdef BN_MP_REDUCE_C
6564 /* LibTomMath, multiple-precision integer library -- Tom St Denis
6566 * LibTomMath is a library that provides multiple-precision
6567 * integer arithmetic as well as number theoretic functionality.
6569 * The library was designed directly after the MPI library by
6570 * Michael Fromberger but has been written from scratch with
6571 * additional optimizations in place.
6573 * The library is free for all purposes without any express
6574 * guarantee it works.
6576 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
6579 /* reduces x mod m, assumes 0 < x < m**2, mu is
6580 * precomputed via tma_mp_reduce_setup.
6581 * From HAC pp.604 Algorithm 14.42
6583 int tma_mp_reduce (tma_mp_int * x, tma_mp_int * m, tma_mp_int * mu)
6586 int res, um = m->used;
6589 if ((res = tma_mp_init_copy (&q, x)) != MP_OKAY) {
6593 /* q1 = x / b**(k-1) */
6594 tma_mp_rshd (&q, um - 1);
6596 /* according to HAC this optimization is ok */
6597 if (((unsigned long) um) > (((tma_mp_digit)1) << (DIGIT_BIT - 1))) {
6598 if ((res = tma_mp_mul (&q, mu, &q)) != MP_OKAY) {
6602 #ifdef BN_S_MP_MUL_HIGH_DIGS_C
6603 if ((res = s_tma_mp_mul_high_digs (&q, mu, &q, um)) != MP_OKAY) {
6606 #elif defined(BN_FAST_S_MP_MUL_HIGH_DIGS_C)
6607 if ((res = fast_s_tma_mp_mul_high_digs (&q, mu, &q, um)) != MP_OKAY) {
6618 /* q3 = q2 / b**(k+1) */
6619 tma_mp_rshd (&q, um + 1);
6621 /* x = x mod b**(k+1), quick (no division) */
6622 if ((res = tma_mp_mod_2d (x, DIGIT_BIT * (um + 1), x)) != MP_OKAY) {
6626 /* q = q * m mod b**(k+1), quick (no division) */
6627 if ((res = s_tma_mp_mul_digs (&q, m, &q, um + 1)) != MP_OKAY) {
6632 if ((res = tma_mp_sub (x, &q, x)) != MP_OKAY) {
6636 /* If x < 0, add b**(k+1) to it */
6637 if (tma_mp_cmp_d (x, 0) == MP_LT) {
6639 if ((res = tma_mp_lshd (&q, um + 1)) != MP_OKAY)
6641 if ((res = tma_mp_add (x, &q, x)) != MP_OKAY)
6645 /* Back off if it's too big */
6646 while (tma_mp_cmp (x, m) != MP_LT) {
6647 if ((res = s_tma_mp_sub (x, m, x)) != MP_OKAY) {
6663 /* End: bn_tma_mp_reduce.c */
6665 /* Start: bn_tma_mp_reduce_2k.c */
6667 #ifdef BN_MP_REDUCE_2K_C
6668 /* LibTomMath, multiple-precision integer library -- Tom St Denis
6670 * LibTomMath is a library that provides multiple-precision
6671 * integer arithmetic as well as number theoretic functionality.
6673 * The library was designed directly after the MPI library by
6674 * Michael Fromberger but has been written from scratch with
6675 * additional optimizations in place.
6677 * The library is free for all purposes without any express
6678 * guarantee it works.
6680 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
6683 /* reduces a modulo n where n is of the form 2**p - d */
6684 int tma_mp_reduce_2k(tma_mp_int *a, tma_mp_int *n, tma_mp_digit d)
6689 if ((res = tma_mp_init(&q)) != MP_OKAY) {
6693 p = tma_mp_count_bits(n);
6695 /* q = a/2**p, a = a mod 2**p */
6696 if ((res = tma_mp_div_2d(a, p, &q, a)) != MP_OKAY) {
6702 if ((res = tma_mp_mul_d(&q, d, &q)) != MP_OKAY) {
6708 if ((res = s_tma_mp_add(a, &q, a)) != MP_OKAY) {
6712 if (tma_mp_cmp_mag(a, n) != MP_LT) {
6713 s_tma_mp_sub(a, n, a);
6728 /* End: bn_tma_mp_reduce_2k.c */
6730 /* Start: bn_tma_mp_reduce_2k_l.c */
6732 #ifdef BN_MP_REDUCE_2K_L_C
6733 /* LibTomMath, multiple-precision integer library -- Tom St Denis
6735 * LibTomMath is a library that provides multiple-precision
6736 * integer arithmetic as well as number theoretic functionality.
6738 * The library was designed directly after the MPI library by
6739 * Michael Fromberger but has been written from scratch with
6740 * additional optimizations in place.
6742 * The library is free for all purposes without any express
6743 * guarantee it works.
6745 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
6748 /* reduces a modulo n where n is of the form 2**p - d
6749 This differs from reduce_2k since "d" can be larger
6750 than a single digit.
6752 int tma_mp_reduce_2k_l(tma_mp_int *a, tma_mp_int *n, tma_mp_int *d)
6757 if ((res = tma_mp_init(&q)) != MP_OKAY) {
6761 p = tma_mp_count_bits(n);
6763 /* q = a/2**p, a = a mod 2**p */
6764 if ((res = tma_mp_div_2d(a, p, &q, a)) != MP_OKAY) {
6769 if ((res = tma_mp_mul(&q, d, &q)) != MP_OKAY) {
6774 if ((res = s_tma_mp_add(a, &q, a)) != MP_OKAY) {
6778 if (tma_mp_cmp_mag(a, n) != MP_LT) {
6779 s_tma_mp_sub(a, n, a);
6794 /* End: bn_tma_mp_reduce_2k_l.c */
6796 /* Start: bn_tma_mp_reduce_2k_setup.c */
6798 #ifdef BN_MP_REDUCE_2K_SETUP_C
6799 /* LibTomMath, multiple-precision integer library -- Tom St Denis
6801 * LibTomMath is a library that provides multiple-precision
6802 * integer arithmetic as well as number theoretic functionality.
6804 * The library was designed directly after the MPI library by
6805 * Michael Fromberger but has been written from scratch with
6806 * additional optimizations in place.
6808 * The library is free for all purposes without any express
6809 * guarantee it works.
6811 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
6814 /* determines the setup value */
6815 int tma_mp_reduce_2k_setup(tma_mp_int *a, tma_mp_digit *d)
6820 if ((res = tma_mp_init(&tmp)) != MP_OKAY) {
6824 p = tma_mp_count_bits(a);
6825 if ((res = tma_mp_2expt(&tmp, p)) != MP_OKAY) {
6830 if ((res = s_tma_mp_sub(&tmp, a, &tmp)) != MP_OKAY) {
6845 /* End: bn_tma_mp_reduce_2k_setup.c */
6847 /* Start: bn_tma_mp_reduce_2k_setup_l.c */
6849 #ifdef BN_MP_REDUCE_2K_SETUP_L_C
6850 /* LibTomMath, multiple-precision integer library -- Tom St Denis
6852 * LibTomMath is a library that provides multiple-precision
6853 * integer arithmetic as well as number theoretic functionality.
6855 * The library was designed directly after the MPI library by
6856 * Michael Fromberger but has been written from scratch with
6857 * additional optimizations in place.
6859 * The library is free for all purposes without any express
6860 * guarantee it works.
6862 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
6865 /* determines the setup value */
6866 int tma_mp_reduce_2k_setup_l(tma_mp_int *a, tma_mp_int *d)
6871 if ((res = tma_mp_init(&tmp)) != MP_OKAY) {
6875 if ((res = tma_mp_2expt(&tmp, tma_mp_count_bits(a))) != MP_OKAY) {
6879 if ((res = s_tma_mp_sub(&tmp, a, d)) != MP_OKAY) {
6893 /* End: bn_tma_mp_reduce_2k_setup_l.c */
6895 /* Start: bn_tma_mp_reduce_is_2k.c */
6897 #ifdef BN_MP_REDUCE_IS_2K_C
6898 /* LibTomMath, multiple-precision integer library -- Tom St Denis
6900 * LibTomMath is a library that provides multiple-precision
6901 * integer arithmetic as well as number theoretic functionality.
6903 * The library was designed directly after the MPI library by
6904 * Michael Fromberger but has been written from scratch with
6905 * additional optimizations in place.
6907 * The library is free for all purposes without any express
6908 * guarantee it works.
6910 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
6913 /* determines if tma_mp_reduce_2k can be used */
6914 int tma_mp_reduce_is_2k(tma_mp_int *a)
6921 } else if (a->used == 1) {
6923 } else if (a->used > 1) {
6924 iy = tma_mp_count_bits(a);
6928 /* Test every bit from the second digit up, must be 1 */
6929 for (ix = DIGIT_BIT; ix < iy; ix++) {
6930 if ((a->dp[iw] & iz) == 0) {
6934 if (iz > (tma_mp_digit)MP_MASK) {
6949 /* End: bn_tma_mp_reduce_is_2k.c */
6951 /* Start: bn_tma_mp_reduce_is_2k_l.c */
6953 #ifdef BN_MP_REDUCE_IS_2K_L_C
6954 /* LibTomMath, multiple-precision integer library -- Tom St Denis
6956 * LibTomMath is a library that provides multiple-precision
6957 * integer arithmetic as well as number theoretic functionality.
6959 * The library was designed directly after the MPI library by
6960 * Michael Fromberger but has been written from scratch with
6961 * additional optimizations in place.
6963 * The library is free for all purposes without any express
6964 * guarantee it works.
6966 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
6969 /* determines if reduce_2k_l can be used */
6970 int tma_mp_reduce_is_2k_l(tma_mp_int *a)
6976 } else if (a->used == 1) {
6978 } else if (a->used > 1) {
6979 /* if more than half of the digits are -1 we're sold */
6980 for (iy = ix = 0; ix < a->used; ix++) {
6981 if (a->dp[ix] == MP_MASK) {
6985 return (iy >= (a->used/2)) ? MP_YES : MP_NO;
6997 /* End: bn_tma_mp_reduce_is_2k_l.c */
6999 /* Start: bn_tma_mp_reduce_setup.c */
7001 #ifdef BN_MP_REDUCE_SETUP_C
7002 /* LibTomMath, multiple-precision integer library -- Tom St Denis
7004 * LibTomMath is a library that provides multiple-precision
7005 * integer arithmetic as well as number theoretic functionality.
7007 * The library was designed directly after the MPI library by
7008 * Michael Fromberger but has been written from scratch with
7009 * additional optimizations in place.
7011 * The library is free for all purposes without any express
7012 * guarantee it works.
7014 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
7017 /* pre-calculate the value required for Barrett reduction
7018 * For a given modulus "b" it calulates the value required in "a"
7020 int tma_mp_reduce_setup (tma_mp_int * a, tma_mp_int * b)
7024 if ((res = tma_mp_2expt (a, b->used * 2 * DIGIT_BIT)) != MP_OKAY) {
7027 return tma_mp_div (a, b, a, NULL);
7035 /* End: bn_tma_mp_reduce_setup.c */
7037 /* Start: bn_tma_mp_rshd.c */
7040 /* LibTomMath, multiple-precision integer library -- Tom St Denis
7042 * LibTomMath is a library that provides multiple-precision
7043 * integer arithmetic as well as number theoretic functionality.
7045 * The library was designed directly after the MPI library by
7046 * Michael Fromberger but has been written from scratch with
7047 * additional optimizations in place.
7049 * The library is free for all purposes without any express
7050 * guarantee it works.
7052 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
7055 /* shift right a certain amount of digits */
7056 void tma_mp_rshd (tma_mp_int * a, int b)
7060 /* if b <= 0 then ignore it */
7065 /* if b > used then simply zero it and return */
7072 register tma_mp_digit *bottom, *top;
7074 /* shift the digits down */
7079 /* top [offset into digits] */
7082 /* this is implemented as a sliding window where
7083 * the window is b-digits long and digits from
7084 * the top of the window are copied to the bottom
7088 b-2 | b-1 | b0 | b1 | b2 | ... | bb | ---->
7090 \-------------------/ ---->
7092 for (x = 0; x < (a->used - b); x++) {
7096 /* zero the top digits */
7097 for (; x < a->used; x++) {
7102 /* remove excess digits */
7111 /* End: bn_tma_mp_rshd.c */
7113 /* Start: bn_tma_mp_set.c */
7116 /* LibTomMath, multiple-precision integer library -- Tom St Denis
7118 * LibTomMath is a library that provides multiple-precision
7119 * integer arithmetic as well as number theoretic functionality.
7121 * The library was designed directly after the MPI library by
7122 * Michael Fromberger but has been written from scratch with
7123 * additional optimizations in place.
7125 * The library is free for all purposes without any express
7126 * guarantee it works.
7128 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
7131 /* set to a digit */
7132 void tma_mp_set (tma_mp_int * a, tma_mp_digit b)
7135 a->dp[0] = b & MP_MASK;
7136 a->used = (a->dp[0] != 0) ? 1 : 0;
7144 /* End: bn_tma_mp_set.c */
7146 /* Start: bn_tma_mp_set_int.c */
7148 #ifdef BN_MP_SET_INT_C
7149 /* LibTomMath, multiple-precision integer library -- Tom St Denis
7151 * LibTomMath is a library that provides multiple-precision
7152 * integer arithmetic as well as number theoretic functionality.
7154 * The library was designed directly after the MPI library by
7155 * Michael Fromberger but has been written from scratch with
7156 * additional optimizations in place.
7158 * The library is free for all purposes without any express
7159 * guarantee it works.
7161 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
7164 /* set a 32-bit const */
7165 int tma_mp_set_int (tma_mp_int * a, unsigned long b)
7171 /* set four bits at a time */
7172 for (x = 0; x < 8; x++) {
7173 /* shift the number up four bits */
7174 if ((res = tma_mp_mul_2d (a, 4, a)) != MP_OKAY) {
7178 /* OR in the top four bits of the source */
7179 a->dp[0] |= (b >> 28) & 15;
7181 /* shift the source up to the next four bits */
7184 /* ensure that digits are not clamped off */
7196 /* End: bn_tma_mp_set_int.c */
7198 /* Start: bn_tma_mp_shrink.c */
7200 #ifdef BN_MP_SHRINK_C
7201 /* LibTomMath, multiple-precision integer library -- Tom St Denis
7203 * LibTomMath is a library that provides multiple-precision
7204 * integer arithmetic as well as number theoretic functionality.
7206 * The library was designed directly after the MPI library by
7207 * Michael Fromberger but has been written from scratch with
7208 * additional optimizations in place.
7210 * The library is free for all purposes without any express
7211 * guarantee it works.
7213 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
7216 /* shrink a bignum */
7217 int tma_mp_shrink (tma_mp_int * a)
7220 if (a->alloc != a->used && a->used > 0) {
7221 if ((tmp = OPT_CAST(tma_mp_digit) XREALLOC (a->dp, sizeof (tma_mp_digit) * a->used)) == NULL) {
7235 /* End: bn_tma_mp_shrink.c */
7237 /* Start: bn_tma_mp_signed_bin_size.c */
7239 #ifdef BN_MP_SIGNED_BIN_SIZE_C
7240 /* LibTomMath, multiple-precision integer library -- Tom St Denis
7242 * LibTomMath is a library that provides multiple-precision
7243 * integer arithmetic as well as number theoretic functionality.
7245 * The library was designed directly after the MPI library by
7246 * Michael Fromberger but has been written from scratch with
7247 * additional optimizations in place.
7249 * The library is free for all purposes without any express
7250 * guarantee it works.
7252 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
7255 /* get the size for an signed equivalent */
7256 int tma_mp_signed_bin_size (tma_mp_int * a)
7258 return 1 + tma_mp_unsigned_bin_size (a);
7266 /* End: bn_tma_mp_signed_bin_size.c */
7268 /* Start: bn_tma_mp_sqr.c */
7271 /* LibTomMath, multiple-precision integer library -- Tom St Denis
7273 * LibTomMath is a library that provides multiple-precision
7274 * integer arithmetic as well as number theoretic functionality.
7276 * The library was designed directly after the MPI library by
7277 * Michael Fromberger but has been written from scratch with
7278 * additional optimizations in place.
7280 * The library is free for all purposes without any express
7281 * guarantee it works.
7283 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
7286 /* computes b = a*a */
7288 tma_mp_sqr (tma_mp_int * a, tma_mp_int * b)
7292 #ifdef BN_MP_TOOM_SQR_C
7293 /* use Toom-Cook? */
7294 if (a->used >= TOOM_SQR_CUTOFF) {
7295 res = tma_mp_toom_sqr(a, b);
7299 #ifdef BN_MP_KARATSUBA_SQR_C
7300 if (a->used >= KARATSUBA_SQR_CUTOFF) {
7301 res = tma_mp_karatsuba_sqr (a, b);
7305 #ifdef BN_FAST_S_MP_SQR_C
7306 /* can we use the fast comba multiplier? */
7307 if ((a->used * 2 + 1) < MP_WARRAY &&
7309 (1 << (sizeof(tma_mp_word) * CHAR_BIT - 2*DIGIT_BIT - 1))) {
7310 res = fast_s_tma_mp_sqr (a, b);
7313 #ifdef BN_S_MP_SQR_C
7314 res = s_tma_mp_sqr (a, b);
7328 /* End: bn_tma_mp_sqr.c */
7330 /* Start: bn_tma_mp_sqrmod.c */
7332 #ifdef BN_MP_SQRMOD_C
7333 /* LibTomMath, multiple-precision integer library -- Tom St Denis
7335 * LibTomMath is a library that provides multiple-precision
7336 * integer arithmetic as well as number theoretic functionality.
7338 * The library was designed directly after the MPI library by
7339 * Michael Fromberger but has been written from scratch with
7340 * additional optimizations in place.
7342 * The library is free for all purposes without any express
7343 * guarantee it works.
7345 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
7348 /* c = a * a (mod b) */
7350 tma_mp_sqrmod (tma_mp_int * a, tma_mp_int * b, tma_mp_int * c)
7355 if ((res = tma_mp_init (&t)) != MP_OKAY) {
7359 if ((res = tma_mp_sqr (a, &t)) != MP_OKAY) {
7363 res = tma_mp_mod (&t, b, c);
7373 /* End: bn_tma_mp_sqrmod.c */
7375 /* Start: bn_tma_mp_sqrt.c */
7378 /* LibTomMath, multiple-precision integer library -- Tom St Denis
7380 * LibTomMath is a library that provides multiple-precision
7381 * integer arithmetic as well as number theoretic functionality.
7383 * The library was designed directly after the MPI library by
7384 * Michael Fromberger but has been written from scratch with
7385 * additional optimizations in place.
7387 * The library is free for all purposes without any express
7388 * guarantee it works.
7390 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
7393 /* this function is less generic than tma_mp_n_root, simpler and faster */
7394 int tma_mp_sqrt(tma_mp_int *arg, tma_mp_int *ret)
7399 /* must be positive */
7400 if (arg->sign == MP_NEG) {
7405 if (tma_mp_iszero(arg) == MP_YES) {
7410 if ((res = tma_mp_init_copy(&t1, arg)) != MP_OKAY) {
7414 if ((res = tma_mp_init(&t2)) != MP_OKAY) {
7418 /* First approx. (not very bad for large arg) */
7419 tma_mp_rshd (&t1,t1.used/2);
7422 if ((res = tma_mp_div(arg,&t1,&t2,NULL)) != MP_OKAY) {
7425 if ((res = tma_mp_add(&t1,&t2,&t1)) != MP_OKAY) {
7428 if ((res = tma_mp_div_2(&t1,&t1)) != MP_OKAY) {
7431 /* And now t1 > sqrt(arg) */
7433 if ((res = tma_mp_div(arg,&t1,&t2,NULL)) != MP_OKAY) {
7436 if ((res = tma_mp_add(&t1,&t2,&t1)) != MP_OKAY) {
7439 if ((res = tma_mp_div_2(&t1,&t1)) != MP_OKAY) {
7442 /* t1 >= sqrt(arg) >= t2 at this point */
7443 } while (tma_mp_cmp_mag(&t1,&t2) == MP_GT);
7445 tma_mp_exch(&t1,ret);
7447 E1: tma_mp_clear(&t2);
7448 E2: tma_mp_clear(&t1);
7458 /* End: bn_tma_mp_sqrt.c */
7460 /* Start: bn_tma_mp_sub.c */
7463 /* LibTomMath, multiple-precision integer library -- Tom St Denis
7465 * LibTomMath is a library that provides multiple-precision
7466 * integer arithmetic as well as number theoretic functionality.
7468 * The library was designed directly after the MPI library by
7469 * Michael Fromberger but has been written from scratch with
7470 * additional optimizations in place.
7472 * The library is free for all purposes without any express
7473 * guarantee it works.
7475 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
7478 /* high level subtraction (handles signs) */
7480 tma_mp_sub (tma_mp_int * a, tma_mp_int * b, tma_mp_int * c)
7488 /* subtract a negative from a positive, OR */
7489 /* subtract a positive from a negative. */
7490 /* In either case, ADD their magnitudes, */
7491 /* and use the sign of the first number. */
7493 res = s_tma_mp_add (a, b, c);
7495 /* subtract a positive from a positive, OR */
7496 /* subtract a negative from a negative. */
7497 /* First, take the difference between their */
7498 /* magnitudes, then... */
7499 if (tma_mp_cmp_mag (a, b) != MP_LT) {
7500 /* Copy the sign from the first */
7502 /* The first has a larger or equal magnitude */
7503 res = s_tma_mp_sub (a, b, c);
7505 /* The result has the *opposite* sign from */
7506 /* the first number. */
7507 c->sign = (sa == MP_ZPOS) ? MP_NEG : MP_ZPOS;
7508 /* The second has a larger magnitude */
7509 res = s_tma_mp_sub (b, a, c);
7521 /* End: bn_tma_mp_sub.c */
7523 /* Start: bn_tma_mp_sub_d.c */
7525 #ifdef BN_MP_SUB_D_C
7526 /* LibTomMath, multiple-precision integer library -- Tom St Denis
7528 * LibTomMath is a library that provides multiple-precision
7529 * integer arithmetic as well as number theoretic functionality.
7531 * The library was designed directly after the MPI library by
7532 * Michael Fromberger but has been written from scratch with
7533 * additional optimizations in place.
7535 * The library is free for all purposes without any express
7536 * guarantee it works.
7538 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
7541 /* single digit subtraction */
7543 tma_mp_sub_d (tma_mp_int * a, tma_mp_digit b, tma_mp_int * c)
7545 tma_mp_digit *tmpa, *tmpc, mu;
7546 int res, ix, oldused;
7548 /* grow c as required */
7549 if (c->alloc < a->used + 1) {
7550 if ((res = tma_mp_grow(c, a->used + 1)) != MP_OKAY) {
7555 /* if a is negative just do an unsigned
7556 * addition [with fudged signs]
7558 if (a->sign == MP_NEG) {
7560 res = tma_mp_add_d(a, b, c);
7561 a->sign = c->sign = MP_NEG;
7574 /* if a <= b simply fix the single digit */
7575 if ((a->used == 1 && a->dp[0] <= b) || a->used == 0) {
7577 *tmpc++ = b - *tmpa;
7583 /* negative/1digit */
7591 /* subtract first digit */
7592 *tmpc = *tmpa++ - b;
7593 mu = *tmpc >> (sizeof(tma_mp_digit) * CHAR_BIT - 1);
7596 /* handle rest of the digits */
7597 for (ix = 1; ix < a->used; ix++) {
7598 *tmpc = *tmpa++ - mu;
7599 mu = *tmpc >> (sizeof(tma_mp_digit) * CHAR_BIT - 1);
7604 /* zero excess digits */
7605 while (ix++ < oldused) {
7618 /* End: bn_tma_mp_sub_d.c */
7620 /* Start: bn_tma_mp_submod.c */
7622 #ifdef BN_MP_SUBMOD_C
7623 /* LibTomMath, multiple-precision integer library -- Tom St Denis
7625 * LibTomMath is a library that provides multiple-precision
7626 * integer arithmetic as well as number theoretic functionality.
7628 * The library was designed directly after the MPI library by
7629 * Michael Fromberger but has been written from scratch with
7630 * additional optimizations in place.
7632 * The library is free for all purposes without any express
7633 * guarantee it works.
7635 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
7638 /* d = a - b (mod c) */
7640 tma_mp_submod (tma_mp_int * a, tma_mp_int * b, tma_mp_int * c, tma_mp_int * d)
7646 if ((res = tma_mp_init (&t)) != MP_OKAY) {
7650 if ((res = tma_mp_sub (a, b, &t)) != MP_OKAY) {
7654 res = tma_mp_mod (&t, c, d);
7664 /* End: bn_tma_mp_submod.c */
7666 /* Start: bn_tma_mp_to_signed_bin.c */
7668 #ifdef BN_MP_TO_SIGNED_BIN_C
7669 /* LibTomMath, multiple-precision integer library -- Tom St Denis
7671 * LibTomMath is a library that provides multiple-precision
7672 * integer arithmetic as well as number theoretic functionality.
7674 * The library was designed directly after the MPI library by
7675 * Michael Fromberger but has been written from scratch with
7676 * additional optimizations in place.
7678 * The library is free for all purposes without any express
7679 * guarantee it works.
7681 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
7684 /* store in signed [big endian] format */
7685 int tma_mp_to_signed_bin (tma_mp_int * a, unsigned char *b)
7689 if ((res = tma_mp_to_unsigned_bin (a, b + 1)) != MP_OKAY) {
7692 b[0] = (unsigned char) ((a->sign == MP_ZPOS) ? 0 : 1);
7701 /* End: bn_tma_mp_to_signed_bin.c */
7703 /* Start: bn_tma_mp_to_signed_bin_n.c */
7705 #ifdef BN_MP_TO_SIGNED_BIN_N_C
7706 /* LibTomMath, multiple-precision integer library -- Tom St Denis
7708 * LibTomMath is a library that provides multiple-precision
7709 * integer arithmetic as well as number theoretic functionality.
7711 * The library was designed directly after the MPI library by
7712 * Michael Fromberger but has been written from scratch with
7713 * additional optimizations in place.
7715 * The library is free for all purposes without any express
7716 * guarantee it works.
7718 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
7721 /* store in signed [big endian] format */
7722 int tma_mp_to_signed_bin_n (tma_mp_int * a, unsigned char *b, unsigned long *outlen)
7724 if (*outlen < (unsigned long)tma_mp_signed_bin_size(a)) {
7727 *outlen = tma_mp_signed_bin_size(a);
7728 return tma_mp_to_signed_bin(a, b);
7736 /* End: bn_tma_mp_to_signed_bin_n.c */
7738 /* Start: bn_tma_mp_to_unsigned_bin.c */
7740 #ifdef BN_MP_TO_UNSIGNED_BIN_C
7741 /* LibTomMath, multiple-precision integer library -- Tom St Denis
7743 * LibTomMath is a library that provides multiple-precision
7744 * integer arithmetic as well as number theoretic functionality.
7746 * The library was designed directly after the MPI library by
7747 * Michael Fromberger but has been written from scratch with
7748 * additional optimizations in place.
7750 * The library is free for all purposes without any express
7751 * guarantee it works.
7753 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
7756 /* store in unsigned [big endian] format */
7757 int tma_mp_to_unsigned_bin (tma_mp_int * a, unsigned char *b)
7762 if ((res = tma_mp_init_copy (&t, a)) != MP_OKAY) {
7767 while (tma_mp_iszero (&t) == 0) {
7769 b[x++] = (unsigned char) (t.dp[0] & 255);
7771 b[x++] = (unsigned char) (t.dp[0] | ((t.dp[1] & 0x01) << 7));
7773 if ((res = tma_mp_div_2d (&t, 8, &t, NULL)) != MP_OKAY) {
7788 /* End: bn_tma_mp_to_unsigned_bin.c */
7790 /* Start: bn_tma_mp_to_unsigned_bin_n.c */
7792 #ifdef BN_MP_TO_UNSIGNED_BIN_N_C
7793 /* LibTomMath, multiple-precision integer library -- Tom St Denis
7795 * LibTomMath is a library that provides multiple-precision
7796 * integer arithmetic as well as number theoretic functionality.
7798 * The library was designed directly after the MPI library by
7799 * Michael Fromberger but has been written from scratch with
7800 * additional optimizations in place.
7802 * The library is free for all purposes without any express
7803 * guarantee it works.
7805 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
7808 /* store in unsigned [big endian] format */
7809 int tma_mp_to_unsigned_bin_n (tma_mp_int * a, unsigned char *b, unsigned long *outlen)
7811 if (*outlen < (unsigned long)tma_mp_unsigned_bin_size(a)) {
7814 *outlen = tma_mp_unsigned_bin_size(a);
7815 return tma_mp_to_unsigned_bin(a, b);
7823 /* End: bn_tma_mp_to_unsigned_bin_n.c */
7825 /* Start: bn_tma_mp_toom_mul.c */
7827 #ifdef BN_MP_TOOM_MUL_C
7828 /* LibTomMath, multiple-precision integer library -- Tom St Denis
7830 * LibTomMath is a library that provides multiple-precision
7831 * integer arithmetic as well as number theoretic functionality.
7833 * The library was designed directly after the MPI library by
7834 * Michael Fromberger but has been written from scratch with
7835 * additional optimizations in place.
7837 * The library is free for all purposes without any express
7838 * guarantee it works.
7840 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
7843 /* multiplication using the Toom-Cook 3-way algorithm
7845 * Much more complicated than Karatsuba but has a lower
7846 * asymptotic running time of O(N**1.464). This algorithm is
7847 * only particularly useful on VERY large inputs
7848 * (we're talking 1000s of digits here...).
7850 int tma_mp_toom_mul(tma_mp_int *a, tma_mp_int *b, tma_mp_int *c)
7852 tma_mp_int w0, w1, w2, w3, w4, tmp1, tmp2, a0, a1, a2, b0, b1, b2;
7856 if ((res = tma_mp_init_multi(&w0, &w1, &w2, &w3, &w4,
7857 &a0, &a1, &a2, &b0, &b1,
7858 &b2, &tmp1, &tmp2, NULL)) != MP_OKAY) {
7863 B = MIN(a->used, b->used) / 3;
7865 /* a = a2 * B**2 + a1 * B + a0 */
7866 if ((res = tma_mp_mod_2d(a, DIGIT_BIT * B, &a0)) != MP_OKAY) {
7870 if ((res = tma_mp_copy(a, &a1)) != MP_OKAY) {
7873 tma_mp_rshd(&a1, B);
7874 tma_mp_mod_2d(&a1, DIGIT_BIT * B, &a1);
7876 if ((res = tma_mp_copy(a, &a2)) != MP_OKAY) {
7879 tma_mp_rshd(&a2, B*2);
7881 /* b = b2 * B**2 + b1 * B + b0 */
7882 if ((res = tma_mp_mod_2d(b, DIGIT_BIT * B, &b0)) != MP_OKAY) {
7886 if ((res = tma_mp_copy(b, &b1)) != MP_OKAY) {
7889 tma_mp_rshd(&b1, B);
7890 tma_mp_mod_2d(&b1, DIGIT_BIT * B, &b1);
7892 if ((res = tma_mp_copy(b, &b2)) != MP_OKAY) {
7895 tma_mp_rshd(&b2, B*2);
7898 if ((res = tma_mp_mul(&a0, &b0, &w0)) != MP_OKAY) {
7903 if ((res = tma_mp_mul(&a2, &b2, &w4)) != MP_OKAY) {
7907 /* w1 = (a2 + 2(a1 + 2a0))(b2 + 2(b1 + 2b0)) */
7908 if ((res = tma_mp_mul_2(&a0, &tmp1)) != MP_OKAY) {
7911 if ((res = tma_mp_add(&tmp1, &a1, &tmp1)) != MP_OKAY) {
7914 if ((res = tma_mp_mul_2(&tmp1, &tmp1)) != MP_OKAY) {
7917 if ((res = tma_mp_add(&tmp1, &a2, &tmp1)) != MP_OKAY) {
7921 if ((res = tma_mp_mul_2(&b0, &tmp2)) != MP_OKAY) {
7924 if ((res = tma_mp_add(&tmp2, &b1, &tmp2)) != MP_OKAY) {
7927 if ((res = tma_mp_mul_2(&tmp2, &tmp2)) != MP_OKAY) {
7930 if ((res = tma_mp_add(&tmp2, &b2, &tmp2)) != MP_OKAY) {
7934 if ((res = tma_mp_mul(&tmp1, &tmp2, &w1)) != MP_OKAY) {
7938 /* w3 = (a0 + 2(a1 + 2a2))(b0 + 2(b1 + 2b2)) */
7939 if ((res = tma_mp_mul_2(&a2, &tmp1)) != MP_OKAY) {
7942 if ((res = tma_mp_add(&tmp1, &a1, &tmp1)) != MP_OKAY) {
7945 if ((res = tma_mp_mul_2(&tmp1, &tmp1)) != MP_OKAY) {
7948 if ((res = tma_mp_add(&tmp1, &a0, &tmp1)) != MP_OKAY) {
7952 if ((res = tma_mp_mul_2(&b2, &tmp2)) != MP_OKAY) {
7955 if ((res = tma_mp_add(&tmp2, &b1, &tmp2)) != MP_OKAY) {
7958 if ((res = tma_mp_mul_2(&tmp2, &tmp2)) != MP_OKAY) {
7961 if ((res = tma_mp_add(&tmp2, &b0, &tmp2)) != MP_OKAY) {
7965 if ((res = tma_mp_mul(&tmp1, &tmp2, &w3)) != MP_OKAY) {
7970 /* w2 = (a2 + a1 + a0)(b2 + b1 + b0) */
7971 if ((res = tma_mp_add(&a2, &a1, &tmp1)) != MP_OKAY) {
7974 if ((res = tma_mp_add(&tmp1, &a0, &tmp1)) != MP_OKAY) {
7977 if ((res = tma_mp_add(&b2, &b1, &tmp2)) != MP_OKAY) {
7980 if ((res = tma_mp_add(&tmp2, &b0, &tmp2)) != MP_OKAY) {
7983 if ((res = tma_mp_mul(&tmp1, &tmp2, &w2)) != MP_OKAY) {
7987 /* now solve the matrix
7995 using 12 subtractions, 4 shifts,
7996 2 small divisions and 1 small multiplication
8000 if ((res = tma_mp_sub(&w1, &w4, &w1)) != MP_OKAY) {
8004 if ((res = tma_mp_sub(&w3, &w0, &w3)) != MP_OKAY) {
8008 if ((res = tma_mp_div_2(&w1, &w1)) != MP_OKAY) {
8012 if ((res = tma_mp_div_2(&w3, &w3)) != MP_OKAY) {
8016 if ((res = tma_mp_sub(&w2, &w0, &w2)) != MP_OKAY) {
8019 if ((res = tma_mp_sub(&w2, &w4, &w2)) != MP_OKAY) {
8023 if ((res = tma_mp_sub(&w1, &w2, &w1)) != MP_OKAY) {
8027 if ((res = tma_mp_sub(&w3, &w2, &w3)) != MP_OKAY) {
8031 if ((res = tma_mp_mul_2d(&w0, 3, &tmp1)) != MP_OKAY) {
8034 if ((res = tma_mp_sub(&w1, &tmp1, &w1)) != MP_OKAY) {
8038 if ((res = tma_mp_mul_2d(&w4, 3, &tmp1)) != MP_OKAY) {
8041 if ((res = tma_mp_sub(&w3, &tmp1, &w3)) != MP_OKAY) {
8045 if ((res = tma_mp_mul_d(&w2, 3, &w2)) != MP_OKAY) {
8048 if ((res = tma_mp_sub(&w2, &w1, &w2)) != MP_OKAY) {
8051 if ((res = tma_mp_sub(&w2, &w3, &w2)) != MP_OKAY) {
8055 if ((res = tma_mp_sub(&w1, &w2, &w1)) != MP_OKAY) {
8059 if ((res = tma_mp_sub(&w3, &w2, &w3)) != MP_OKAY) {
8063 if ((res = tma_mp_div_3(&w1, &w1, NULL)) != MP_OKAY) {
8067 if ((res = tma_mp_div_3(&w3, &w3, NULL)) != MP_OKAY) {
8071 /* at this point shift W[n] by B*n */
8072 if ((res = tma_mp_lshd(&w1, 1*B)) != MP_OKAY) {
8075 if ((res = tma_mp_lshd(&w2, 2*B)) != MP_OKAY) {
8078 if ((res = tma_mp_lshd(&w3, 3*B)) != MP_OKAY) {
8081 if ((res = tma_mp_lshd(&w4, 4*B)) != MP_OKAY) {
8085 if ((res = tma_mp_add(&w0, &w1, c)) != MP_OKAY) {
8088 if ((res = tma_mp_add(&w2, &w3, &tmp1)) != MP_OKAY) {
8091 if ((res = tma_mp_add(&w4, &tmp1, &tmp1)) != MP_OKAY) {
8094 if ((res = tma_mp_add(&tmp1, c, c)) != MP_OKAY) {
8099 tma_mp_clear_multi(&w0, &w1, &w2, &w3, &w4,
8100 &a0, &a1, &a2, &b0, &b1,
8101 &b2, &tmp1, &tmp2, NULL);
8111 /* End: bn_tma_mp_toom_mul.c */
8113 /* Start: bn_tma_mp_toom_sqr.c */
8115 #ifdef BN_MP_TOOM_SQR_C
8116 /* LibTomMath, multiple-precision integer library -- Tom St Denis
8118 * LibTomMath is a library that provides multiple-precision
8119 * integer arithmetic as well as number theoretic functionality.
8121 * The library was designed directly after the MPI library by
8122 * Michael Fromberger but has been written from scratch with
8123 * additional optimizations in place.
8125 * The library is free for all purposes without any express
8126 * guarantee it works.
8128 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
8131 /* squaring using Toom-Cook 3-way algorithm */
8133 tma_mp_toom_sqr(tma_mp_int *a, tma_mp_int *b)
8135 tma_mp_int w0, w1, w2, w3, w4, tmp1, a0, a1, a2;
8139 if ((res = tma_mp_init_multi(&w0, &w1, &w2, &w3, &w4, &a0, &a1, &a2, &tmp1, NULL)) != MP_OKAY) {
8146 /* a = a2 * B**2 + a1 * B + a0 */
8147 if ((res = tma_mp_mod_2d(a, DIGIT_BIT * B, &a0)) != MP_OKAY) {
8151 if ((res = tma_mp_copy(a, &a1)) != MP_OKAY) {
8154 tma_mp_rshd(&a1, B);
8155 tma_mp_mod_2d(&a1, DIGIT_BIT * B, &a1);
8157 if ((res = tma_mp_copy(a, &a2)) != MP_OKAY) {
8160 tma_mp_rshd(&a2, B*2);
8163 if ((res = tma_mp_sqr(&a0, &w0)) != MP_OKAY) {
8168 if ((res = tma_mp_sqr(&a2, &w4)) != MP_OKAY) {
8172 /* w1 = (a2 + 2(a1 + 2a0))**2 */
8173 if ((res = tma_mp_mul_2(&a0, &tmp1)) != MP_OKAY) {
8176 if ((res = tma_mp_add(&tmp1, &a1, &tmp1)) != MP_OKAY) {
8179 if ((res = tma_mp_mul_2(&tmp1, &tmp1)) != MP_OKAY) {
8182 if ((res = tma_mp_add(&tmp1, &a2, &tmp1)) != MP_OKAY) {
8186 if ((res = tma_mp_sqr(&tmp1, &w1)) != MP_OKAY) {
8190 /* w3 = (a0 + 2(a1 + 2a2))**2 */
8191 if ((res = tma_mp_mul_2(&a2, &tmp1)) != MP_OKAY) {
8194 if ((res = tma_mp_add(&tmp1, &a1, &tmp1)) != MP_OKAY) {
8197 if ((res = tma_mp_mul_2(&tmp1, &tmp1)) != MP_OKAY) {
8200 if ((res = tma_mp_add(&tmp1, &a0, &tmp1)) != MP_OKAY) {
8204 if ((res = tma_mp_sqr(&tmp1, &w3)) != MP_OKAY) {
8209 /* w2 = (a2 + a1 + a0)**2 */
8210 if ((res = tma_mp_add(&a2, &a1, &tmp1)) != MP_OKAY) {
8213 if ((res = tma_mp_add(&tmp1, &a0, &tmp1)) != MP_OKAY) {
8216 if ((res = tma_mp_sqr(&tmp1, &w2)) != MP_OKAY) {
8220 /* now solve the matrix
8228 using 12 subtractions, 4 shifts, 2 small divisions and 1 small multiplication.
8232 if ((res = tma_mp_sub(&w1, &w4, &w1)) != MP_OKAY) {
8236 if ((res = tma_mp_sub(&w3, &w0, &w3)) != MP_OKAY) {
8240 if ((res = tma_mp_div_2(&w1, &w1)) != MP_OKAY) {
8244 if ((res = tma_mp_div_2(&w3, &w3)) != MP_OKAY) {
8248 if ((res = tma_mp_sub(&w2, &w0, &w2)) != MP_OKAY) {
8251 if ((res = tma_mp_sub(&w2, &w4, &w2)) != MP_OKAY) {
8255 if ((res = tma_mp_sub(&w1, &w2, &w1)) != MP_OKAY) {
8259 if ((res = tma_mp_sub(&w3, &w2, &w3)) != MP_OKAY) {
8263 if ((res = tma_mp_mul_2d(&w0, 3, &tmp1)) != MP_OKAY) {
8266 if ((res = tma_mp_sub(&w1, &tmp1, &w1)) != MP_OKAY) {
8270 if ((res = tma_mp_mul_2d(&w4, 3, &tmp1)) != MP_OKAY) {
8273 if ((res = tma_mp_sub(&w3, &tmp1, &w3)) != MP_OKAY) {
8277 if ((res = tma_mp_mul_d(&w2, 3, &w2)) != MP_OKAY) {
8280 if ((res = tma_mp_sub(&w2, &w1, &w2)) != MP_OKAY) {
8283 if ((res = tma_mp_sub(&w2, &w3, &w2)) != MP_OKAY) {
8287 if ((res = tma_mp_sub(&w1, &w2, &w1)) != MP_OKAY) {
8291 if ((res = tma_mp_sub(&w3, &w2, &w3)) != MP_OKAY) {
8295 if ((res = tma_mp_div_3(&w1, &w1, NULL)) != MP_OKAY) {
8299 if ((res = tma_mp_div_3(&w3, &w3, NULL)) != MP_OKAY) {
8303 /* at this point shift W[n] by B*n */
8304 if ((res = tma_mp_lshd(&w1, 1*B)) != MP_OKAY) {
8307 if ((res = tma_mp_lshd(&w2, 2*B)) != MP_OKAY) {
8310 if ((res = tma_mp_lshd(&w3, 3*B)) != MP_OKAY) {
8313 if ((res = tma_mp_lshd(&w4, 4*B)) != MP_OKAY) {
8317 if ((res = tma_mp_add(&w0, &w1, b)) != MP_OKAY) {
8320 if ((res = tma_mp_add(&w2, &w3, &tmp1)) != MP_OKAY) {
8323 if ((res = tma_mp_add(&w4, &tmp1, &tmp1)) != MP_OKAY) {
8326 if ((res = tma_mp_add(&tmp1, b, b)) != MP_OKAY) {
8331 tma_mp_clear_multi(&w0, &w1, &w2, &w3, &w4, &a0, &a1, &a2, &tmp1, NULL);
8341 /* End: bn_tma_mp_toom_sqr.c */
8343 /* Start: bn_tma_mp_toradix.c */
8345 #ifdef BN_MP_TORADIX_C
8346 /* LibTomMath, multiple-precision integer library -- Tom St Denis
8348 * LibTomMath is a library that provides multiple-precision
8349 * integer arithmetic as well as number theoretic functionality.
8351 * The library was designed directly after the MPI library by
8352 * Michael Fromberger but has been written from scratch with
8353 * additional optimizations in place.
8355 * The library is free for all purposes without any express
8356 * guarantee it works.
8358 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
8361 /* stores a bignum as a ASCII string in a given radix (2..64) */
8362 int tma_mp_toradix (tma_mp_int * a, char *str, int radix)
8369 /* check range of the radix */
8370 if (radix < 2 || radix > 64) {
8374 /* quick out if its zero */
8375 if (tma_mp_iszero(a) == 1) {
8381 if ((res = tma_mp_init_copy (&t, a)) != MP_OKAY) {
8385 /* if it is negative output a - */
8386 if (t.sign == MP_NEG) {
8393 while (tma_mp_iszero (&t) == 0) {
8394 if ((res = tma_mp_div_d (&t, (tma_mp_digit) radix, &t, &d)) != MP_OKAY) {
8398 *str++ = tma_mp_s_rmap[d];
8402 /* reverse the digits of the string. In this case _s points
8403 * to the first digit [exluding the sign] of the number]
8405 bn_reverse ((unsigned char *)_s, digs);
8407 /* append a NULL so the string is properly terminated */
8420 /* End: bn_tma_mp_toradix.c */
8422 /* Start: bn_tma_mp_toradix_n.c */
8424 #ifdef BN_MP_TORADIX_N_C
8425 /* LibTomMath, multiple-precision integer library -- Tom St Denis
8427 * LibTomMath is a library that provides multiple-precision
8428 * integer arithmetic as well as number theoretic functionality.
8430 * The library was designed directly after the MPI library by
8431 * Michael Fromberger but has been written from scratch with
8432 * additional optimizations in place.
8434 * The library is free for all purposes without any express
8435 * guarantee it works.
8437 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
8440 /* stores a bignum as a ASCII string in a given radix (2..64)
8442 * Stores upto maxlen-1 chars and always a NULL byte
8444 int tma_mp_toradix_n(tma_mp_int * a, char *str, int radix, int maxlen)
8451 /* check range of the maxlen, radix */
8452 if (maxlen < 2 || radix < 2 || radix > 64) {
8456 /* quick out if its zero */
8457 if (tma_mp_iszero(a) == MP_YES) {
8463 if ((res = tma_mp_init_copy (&t, a)) != MP_OKAY) {
8467 /* if it is negative output a - */
8468 if (t.sign == MP_NEG) {
8469 /* we have to reverse our digits later... but not the - sign!! */
8472 /* store the flag and mark the number as positive */
8476 /* subtract a char */
8481 while (tma_mp_iszero (&t) == 0) {
8486 if ((res = tma_mp_div_d (&t, (tma_mp_digit) radix, &t, &d)) != MP_OKAY) {
8490 *str++ = tma_mp_s_rmap[d];
8494 /* reverse the digits of the string. In this case _s points
8495 * to the first digit [exluding the sign] of the number
8497 bn_reverse ((unsigned char *)_s, digs);
8499 /* append a NULL so the string is properly terminated */
8512 /* End: bn_tma_mp_toradix_n.c */
8514 /* Start: bn_tma_mp_unsigned_bin_size.c */
8516 #ifdef BN_MP_UNSIGNED_BIN_SIZE_C
8517 /* LibTomMath, multiple-precision integer library -- Tom St Denis
8519 * LibTomMath is a library that provides multiple-precision
8520 * integer arithmetic as well as number theoretic functionality.
8522 * The library was designed directly after the MPI library by
8523 * Michael Fromberger but has been written from scratch with
8524 * additional optimizations in place.
8526 * The library is free for all purposes without any express
8527 * guarantee it works.
8529 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
8532 /* get the size for an unsigned equivalent */
8533 int tma_mp_unsigned_bin_size (tma_mp_int * a)
8535 int size = tma_mp_count_bits (a);
8536 return (size / 8 + ((size & 7) != 0 ? 1 : 0));
8544 /* End: bn_tma_mp_unsigned_bin_size.c */
8546 /* Start: bn_tma_mp_xor.c */
8549 /* LibTomMath, multiple-precision integer library -- Tom St Denis
8551 * LibTomMath is a library that provides multiple-precision
8552 * integer arithmetic as well as number theoretic functionality.
8554 * The library was designed directly after the MPI library by
8555 * Michael Fromberger but has been written from scratch with
8556 * additional optimizations in place.
8558 * The library is free for all purposes without any express
8559 * guarantee it works.
8561 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
8564 /* XOR two ints together */
8566 tma_mp_xor (tma_mp_int * a, tma_mp_int * b, tma_mp_int * c)
8571 if (a->used > b->used) {
8572 if ((res = tma_mp_init_copy (&t, a)) != MP_OKAY) {
8578 if ((res = tma_mp_init_copy (&t, b)) != MP_OKAY) {
8585 for (ix = 0; ix < px; ix++) {
8586 t.dp[ix] ^= x->dp[ix];
8589 tma_mp_exch (c, &t);
8599 /* End: bn_tma_mp_xor.c */
8601 /* Start: bn_tma_mp_zero.c */
8604 /* LibTomMath, multiple-precision integer library -- Tom St Denis
8606 * LibTomMath is a library that provides multiple-precision
8607 * integer arithmetic as well as number theoretic functionality.
8609 * The library was designed directly after the MPI library by
8610 * Michael Fromberger but has been written from scratch with
8611 * additional optimizations in place.
8613 * The library is free for all purposes without any express
8614 * guarantee it works.
8616 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
8620 void tma_mp_zero (tma_mp_int * a)
8629 for (n = 0; n < a->alloc; n++) {
8639 /* End: bn_tma_mp_zero.c */
8641 /* Start: bn_prime_tab.c */
8643 #ifdef BN_PRIME_TAB_C
8644 /* LibTomMath, multiple-precision integer library -- Tom St Denis
8646 * LibTomMath is a library that provides multiple-precision
8647 * integer arithmetic as well as number theoretic functionality.
8649 * The library was designed directly after the MPI library by
8650 * Michael Fromberger but has been written from scratch with
8651 * additional optimizations in place.
8653 * The library is free for all purposes without any express
8654 * guarantee it works.
8656 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
8658 const tma_mp_digit ltm_prime_tab[] = {
8659 0x0002, 0x0003, 0x0005, 0x0007, 0x000B, 0x000D, 0x0011, 0x0013,
8660 0x0017, 0x001D, 0x001F, 0x0025, 0x0029, 0x002B, 0x002F, 0x0035,
8661 0x003B, 0x003D, 0x0043, 0x0047, 0x0049, 0x004F, 0x0053, 0x0059,
8662 0x0061, 0x0065, 0x0067, 0x006B, 0x006D, 0x0071, 0x007F,
8665 0x0089, 0x008B, 0x0095, 0x0097, 0x009D, 0x00A3, 0x00A7, 0x00AD,
8666 0x00B3, 0x00B5, 0x00BF, 0x00C1, 0x00C5, 0x00C7, 0x00D3, 0x00DF,
8667 0x00E3, 0x00E5, 0x00E9, 0x00EF, 0x00F1, 0x00FB, 0x0101, 0x0107,
8668 0x010D, 0x010F, 0x0115, 0x0119, 0x011B, 0x0125, 0x0133, 0x0137,
8670 0x0139, 0x013D, 0x014B, 0x0151, 0x015B, 0x015D, 0x0161, 0x0167,
8671 0x016F, 0x0175, 0x017B, 0x017F, 0x0185, 0x018D, 0x0191, 0x0199,
8672 0x01A3, 0x01A5, 0x01AF, 0x01B1, 0x01B7, 0x01BB, 0x01C1, 0x01C9,
8673 0x01CD, 0x01CF, 0x01D3, 0x01DF, 0x01E7, 0x01EB, 0x01F3, 0x01F7,
8674 0x01FD, 0x0209, 0x020B, 0x021D, 0x0223, 0x022D, 0x0233, 0x0239,
8675 0x023B, 0x0241, 0x024B, 0x0251, 0x0257, 0x0259, 0x025F, 0x0265,
8676 0x0269, 0x026B, 0x0277, 0x0281, 0x0283, 0x0287, 0x028D, 0x0293,
8677 0x0295, 0x02A1, 0x02A5, 0x02AB, 0x02B3, 0x02BD, 0x02C5, 0x02CF,
8679 0x02D7, 0x02DD, 0x02E3, 0x02E7, 0x02EF, 0x02F5, 0x02F9, 0x0301,
8680 0x0305, 0x0313, 0x031D, 0x0329, 0x032B, 0x0335, 0x0337, 0x033B,
8681 0x033D, 0x0347, 0x0355, 0x0359, 0x035B, 0x035F, 0x036D, 0x0371,
8682 0x0373, 0x0377, 0x038B, 0x038F, 0x0397, 0x03A1, 0x03A9, 0x03AD,
8683 0x03B3, 0x03B9, 0x03C7, 0x03CB, 0x03D1, 0x03D7, 0x03DF, 0x03E5,
8684 0x03F1, 0x03F5, 0x03FB, 0x03FD, 0x0407, 0x0409, 0x040F, 0x0419,
8685 0x041B, 0x0425, 0x0427, 0x042D, 0x043F, 0x0443, 0x0445, 0x0449,
8686 0x044F, 0x0455, 0x045D, 0x0463, 0x0469, 0x047F, 0x0481, 0x048B,
8688 0x0493, 0x049D, 0x04A3, 0x04A9, 0x04B1, 0x04BD, 0x04C1, 0x04C7,
8689 0x04CD, 0x04CF, 0x04D5, 0x04E1, 0x04EB, 0x04FD, 0x04FF, 0x0503,
8690 0x0509, 0x050B, 0x0511, 0x0515, 0x0517, 0x051B, 0x0527, 0x0529,
8691 0x052F, 0x0551, 0x0557, 0x055D, 0x0565, 0x0577, 0x0581, 0x058F,
8692 0x0593, 0x0595, 0x0599, 0x059F, 0x05A7, 0x05AB, 0x05AD, 0x05B3,
8693 0x05BF, 0x05C9, 0x05CB, 0x05CF, 0x05D1, 0x05D5, 0x05DB, 0x05E7,
8694 0x05F3, 0x05FB, 0x0607, 0x060D, 0x0611, 0x0617, 0x061F, 0x0623,
8695 0x062B, 0x062F, 0x063D, 0x0641, 0x0647, 0x0649, 0x064D, 0x0653
8704 /* End: bn_prime_tab.c */
8706 /* Start: bn_reverse.c */
8709 /* LibTomMath, multiple-precision integer library -- Tom St Denis
8711 * LibTomMath is a library that provides multiple-precision
8712 * integer arithmetic as well as number theoretic functionality.
8714 * The library was designed directly after the MPI library by
8715 * Michael Fromberger but has been written from scratch with
8716 * additional optimizations in place.
8718 * The library is free for all purposes without any express
8719 * guarantee it works.
8721 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
8724 /* reverse an array, used for radix code */
8726 bn_reverse (unsigned char *s, int len)
8747 /* End: bn_reverse.c */
8749 /* Start: bn_s_tma_mp_add.c */
8751 #ifdef BN_S_MP_ADD_C
8752 /* LibTomMath, multiple-precision integer library -- Tom St Denis
8754 * LibTomMath is a library that provides multiple-precision
8755 * integer arithmetic as well as number theoretic functionality.
8757 * The library was designed directly after the MPI library by
8758 * Michael Fromberger but has been written from scratch with
8759 * additional optimizations in place.
8761 * The library is free for all purposes without any express
8762 * guarantee it works.
8764 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
8767 /* low level addition, based on HAC pp.594, Algorithm 14.7 */
8769 s_tma_mp_add (tma_mp_int * a, tma_mp_int * b, tma_mp_int * c)
8772 int olduse, res, min, max;
8774 /* find sizes, we let |a| <= |b| which means we have to sort
8775 * them. "x" will point to the input with the most digits
8777 if (a->used > b->used) {
8788 if (c->alloc < max + 1) {
8789 if ((res = tma_mp_grow (c, max + 1)) != MP_OKAY) {
8794 /* get old used digit count and set new one */
8799 register tma_mp_digit u, *tmpa, *tmpb, *tmpc;
8802 /* alias for digit pointers */
8813 /* zero the carry */
8815 for (i = 0; i < min; i++) {
8816 /* Compute the sum at one digit, T[i] = A[i] + B[i] + U */
8817 *tmpc = *tmpa++ + *tmpb++ + u;
8819 /* U = carry bit of T[i] */
8820 u = *tmpc >> ((tma_mp_digit)DIGIT_BIT);
8822 /* take away carry bit from T[i] */
8826 /* now copy higher words if any, that is in A+B
8827 * if A or B has more digits add those in
8830 for (; i < max; i++) {
8831 /* T[i] = X[i] + U */
8832 *tmpc = x->dp[i] + u;
8834 /* U = carry bit of T[i] */
8835 u = *tmpc >> ((tma_mp_digit)DIGIT_BIT);
8837 /* take away carry bit from T[i] */
8845 /* clear digits above oldused */
8846 for (i = c->used; i < olduse; i++) {
8860 /* End: bn_s_tma_mp_add.c */
8862 /* Start: bn_s_tma_mp_exptmod.c */
8864 #ifdef BN_S_MP_EXPTMOD_C
8865 /* LibTomMath, multiple-precision integer library -- Tom St Denis
8867 * LibTomMath is a library that provides multiple-precision
8868 * integer arithmetic as well as number theoretic functionality.
8870 * The library was designed directly after the MPI library by
8871 * Michael Fromberger but has been written from scratch with
8872 * additional optimizations in place.
8874 * The library is free for all purposes without any express
8875 * guarantee it works.
8877 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
8882 #define TAB_SIZE 256
8885 int s_tma_mp_exptmod (tma_mp_int * G, tma_mp_int * X, tma_mp_int * P, tma_mp_int * Y, int redmode)
8887 tma_mp_int M[TAB_SIZE], res, mu;
8889 int err, bitbuf, bitcpy, bitcnt, mode, digidx, x, y, winsize;
8890 int (*redux)(tma_mp_int*,tma_mp_int*,tma_mp_int*);
8892 /* find window size */
8893 x = tma_mp_count_bits (X);
8896 } else if (x <= 36) {
8898 } else if (x <= 140) {
8900 } else if (x <= 450) {
8902 } else if (x <= 1303) {
8904 } else if (x <= 3529) {
8917 /* init first cell */
8918 if ((err = tma_mp_init(&M[1])) != MP_OKAY) {
8922 /* now init the second half of the array */
8923 for (x = 1<<(winsize-1); x < (1 << winsize); x++) {
8924 if ((err = tma_mp_init(&M[x])) != MP_OKAY) {
8925 for (y = 1<<(winsize-1); y < x; y++) {
8926 tma_mp_clear (&M[y]);
8928 tma_mp_clear(&M[1]);
8933 /* create mu, used for Barrett reduction */
8934 if ((err = tma_mp_init (&mu)) != MP_OKAY) {
8939 if ((err = tma_mp_reduce_setup (&mu, P)) != MP_OKAY) {
8942 redux = tma_mp_reduce;
8944 if ((err = tma_mp_reduce_2k_setup_l (P, &mu)) != MP_OKAY) {
8947 redux = tma_mp_reduce_2k_l;
8952 * The M table contains powers of the base,
8953 * e.g. M[x] = G**x mod P
8955 * The first half of the table is not
8956 * computed though accept for M[0] and M[1]
8958 if ((err = tma_mp_mod (G, P, &M[1])) != MP_OKAY) {
8962 /* compute the value at M[1<<(winsize-1)] by squaring
8963 * M[1] (winsize-1) times
8965 if ((err = tma_mp_copy (&M[1], &M[1 << (winsize - 1)])) != MP_OKAY) {
8969 for (x = 0; x < (winsize - 1); x++) {
8971 if ((err = tma_mp_sqr (&M[1 << (winsize - 1)],
8972 &M[1 << (winsize - 1)])) != MP_OKAY) {
8976 /* reduce modulo P */
8977 if ((err = redux (&M[1 << (winsize - 1)], P, &mu)) != MP_OKAY) {
8982 /* create upper table, that is M[x] = M[x-1] * M[1] (mod P)
8983 * for x = (2**(winsize - 1) + 1) to (2**winsize - 1)
8985 for (x = (1 << (winsize - 1)) + 1; x < (1 << winsize); x++) {
8986 if ((err = tma_mp_mul (&M[x - 1], &M[1], &M[x])) != MP_OKAY) {
8989 if ((err = redux (&M[x], P, &mu)) != MP_OKAY) {
8995 if ((err = tma_mp_init (&res)) != MP_OKAY) {
8998 tma_mp_set (&res, 1);
9000 /* set initial mode and bit cnt */
9004 digidx = X->used - 1;
9009 /* grab next digit as required */
9010 if (--bitcnt == 0) {
9011 /* if digidx == -1 we are out of digits */
9015 /* read next digit and reset the bitcnt */
9016 buf = X->dp[digidx--];
9017 bitcnt = (int) DIGIT_BIT;
9020 /* grab the next msb from the exponent */
9021 y = (buf >> (tma_mp_digit)(DIGIT_BIT - 1)) & 1;
9022 buf <<= (tma_mp_digit)1;
9024 /* if the bit is zero and mode == 0 then we ignore it
9025 * These represent the leading zero bits before the first 1 bit
9026 * in the exponent. Technically this opt is not required but it
9027 * does lower the # of trivial squaring/reductions used
9029 if (mode == 0 && y == 0) {
9033 /* if the bit is zero and mode == 1 then we square */
9034 if (mode == 1 && y == 0) {
9035 if ((err = tma_mp_sqr (&res, &res)) != MP_OKAY) {
9038 if ((err = redux (&res, P, &mu)) != MP_OKAY) {
9044 /* else we add it to the window */
9045 bitbuf |= (y << (winsize - ++bitcpy));
9048 if (bitcpy == winsize) {
9049 /* ok window is filled so square as required and multiply */
9051 for (x = 0; x < winsize; x++) {
9052 if ((err = tma_mp_sqr (&res, &res)) != MP_OKAY) {
9055 if ((err = redux (&res, P, &mu)) != MP_OKAY) {
9061 if ((err = tma_mp_mul (&res, &M[bitbuf], &res)) != MP_OKAY) {
9064 if ((err = redux (&res, P, &mu)) != MP_OKAY) {
9068 /* empty window and reset */
9075 /* if bits remain then square/multiply */
9076 if (mode == 2 && bitcpy > 0) {
9077 /* square then multiply if the bit is set */
9078 for (x = 0; x < bitcpy; x++) {
9079 if ((err = tma_mp_sqr (&res, &res)) != MP_OKAY) {
9082 if ((err = redux (&res, P, &mu)) != MP_OKAY) {
9087 if ((bitbuf & (1 << winsize)) != 0) {
9089 if ((err = tma_mp_mul (&res, &M[1], &res)) != MP_OKAY) {
9092 if ((err = redux (&res, P, &mu)) != MP_OKAY) {
9099 tma_mp_exch (&res, Y);
9101 LBL_RES:tma_mp_clear (&res);
9102 LBL_MU:tma_mp_clear (&mu);
9104 tma_mp_clear(&M[1]);
9105 for (x = 1<<(winsize-1); x < (1 << winsize); x++) {
9106 tma_mp_clear (&M[x]);
9116 /* End: bn_s_tma_mp_exptmod.c */
9118 /* Start: bn_s_tma_mp_mul_digs.c */
9120 #ifdef BN_S_MP_MUL_DIGS_C
9121 /* LibTomMath, multiple-precision integer library -- Tom St Denis
9123 * LibTomMath is a library that provides multiple-precision
9124 * integer arithmetic as well as number theoretic functionality.
9126 * The library was designed directly after the MPI library by
9127 * Michael Fromberger but has been written from scratch with
9128 * additional optimizations in place.
9130 * The library is free for all purposes without any express
9131 * guarantee it works.
9133 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
9136 /* multiplies |a| * |b| and only computes upto digs digits of result
9137 * HAC pp. 595, Algorithm 14.12 Modified so you can control how
9138 * many digits of output are created.
9140 int s_tma_mp_mul_digs (tma_mp_int * a, tma_mp_int * b, tma_mp_int * c, int digs)
9143 int res, pa, pb, ix, iy;
9146 tma_mp_digit tmpx, *tmpt, *tmpy;
9148 /* can we use the fast multiplier? */
9149 if (((digs) < MP_WARRAY) &&
9150 MIN (a->used, b->used) <
9151 (1 << ((CHAR_BIT * sizeof (tma_mp_word)) - (2 * DIGIT_BIT)))) {
9152 return fast_s_tma_mp_mul_digs (a, b, c, digs);
9155 if ((res = tma_mp_init_size (&t, digs)) != MP_OKAY) {
9160 /* compute the digits of the product directly */
9162 for (ix = 0; ix < pa; ix++) {
9163 /* set the carry to zero */
9166 /* limit ourselves to making digs digits of output */
9167 pb = MIN (b->used, digs - ix);
9169 /* setup some aliases */
9170 /* copy of the digit from a used within the nested loop */
9173 /* an alias for the destination shifted ix places */
9176 /* an alias for the digits of b */
9179 /* compute the columns of the output and propagate the carry */
9180 for (iy = 0; iy < pb; iy++) {
9181 /* compute the column as a tma_mp_word */
9182 r = ((tma_mp_word)*tmpt) +
9183 ((tma_mp_word)tmpx) * ((tma_mp_word)*tmpy++) +
9186 /* the new column is the lower part of the result */
9187 *tmpt++ = (tma_mp_digit) (r & ((tma_mp_word) MP_MASK));
9189 /* get the carry word from the result */
9190 u = (tma_mp_digit) (r >> ((tma_mp_word) DIGIT_BIT));
9192 /* set carry if it is placed below digs */
9193 if (ix + iy < digs) {
9199 tma_mp_exch (&t, c);
9210 /* End: bn_s_tma_mp_mul_digs.c */
9212 /* Start: bn_s_tma_mp_mul_high_digs.c */
9214 #ifdef BN_S_MP_MUL_HIGH_DIGS_C
9215 /* LibTomMath, multiple-precision integer library -- Tom St Denis
9217 * LibTomMath is a library that provides multiple-precision
9218 * integer arithmetic as well as number theoretic functionality.
9220 * The library was designed directly after the MPI library by
9221 * Michael Fromberger but has been written from scratch with
9222 * additional optimizations in place.
9224 * The library is free for all purposes without any express
9225 * guarantee it works.
9227 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
9230 /* multiplies |a| * |b| and does not compute the lower digs digits
9231 * [meant to get the higher part of the product]
9234 s_tma_mp_mul_high_digs (tma_mp_int * a, tma_mp_int * b, tma_mp_int * c, int digs)
9237 int res, pa, pb, ix, iy;
9240 tma_mp_digit tmpx, *tmpt, *tmpy;
9242 /* can we use the fast multiplier? */
9243 #ifdef BN_FAST_S_MP_MUL_HIGH_DIGS_C
9244 if (((a->used + b->used + 1) < MP_WARRAY)
9245 && MIN (a->used, b->used) < (1 << ((CHAR_BIT * sizeof (tma_mp_word)) - (2 * DIGIT_BIT)))) {
9246 return fast_s_tma_mp_mul_high_digs (a, b, c, digs);
9250 if ((res = tma_mp_init_size (&t, a->used + b->used + 1)) != MP_OKAY) {
9253 t.used = a->used + b->used + 1;
9257 for (ix = 0; ix < pa; ix++) {
9258 /* clear the carry */
9261 /* left hand side of A[ix] * B[iy] */
9264 /* alias to the address of where the digits will be stored */
9265 tmpt = &(t.dp[digs]);
9267 /* alias for where to read the right hand side from */
9268 tmpy = b->dp + (digs - ix);
9270 for (iy = digs - ix; iy < pb; iy++) {
9271 /* calculate the double precision result */
9272 r = ((tma_mp_word)*tmpt) +
9273 ((tma_mp_word)tmpx) * ((tma_mp_word)*tmpy++) +
9276 /* get the lower part */
9277 *tmpt++ = (tma_mp_digit) (r & ((tma_mp_word) MP_MASK));
9279 /* carry the carry */
9280 u = (tma_mp_digit) (r >> ((tma_mp_word) DIGIT_BIT));
9285 tma_mp_exch (&t, c);
9295 /* End: bn_s_tma_mp_mul_high_digs.c */
9297 /* Start: bn_s_tma_mp_sqr.c */
9299 #ifdef BN_S_MP_SQR_C
9300 /* LibTomMath, multiple-precision integer library -- Tom St Denis
9302 * LibTomMath is a library that provides multiple-precision
9303 * integer arithmetic as well as number theoretic functionality.
9305 * The library was designed directly after the MPI library by
9306 * Michael Fromberger but has been written from scratch with
9307 * additional optimizations in place.
9309 * The library is free for all purposes without any express
9310 * guarantee it works.
9312 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
9315 /* low level squaring, b = a*a, HAC pp.596-597, Algorithm 14.16 */
9316 int s_tma_mp_sqr (tma_mp_int * a, tma_mp_int * b)
9319 int res, ix, iy, pa;
9321 tma_mp_digit u, tmpx, *tmpt;
9324 if ((res = tma_mp_init_size (&t, 2*pa + 1)) != MP_OKAY) {
9328 /* default used is maximum possible size */
9331 for (ix = 0; ix < pa; ix++) {
9332 /* first calculate the digit at 2*ix */
9333 /* calculate double precision result */
9334 r = ((tma_mp_word) t.dp[2*ix]) +
9335 ((tma_mp_word)a->dp[ix])*((tma_mp_word)a->dp[ix]);
9337 /* store lower part in result */
9338 t.dp[ix+ix] = (tma_mp_digit) (r & ((tma_mp_word) MP_MASK));
9341 u = (tma_mp_digit)(r >> ((tma_mp_word) DIGIT_BIT));
9343 /* left hand side of A[ix] * A[iy] */
9346 /* alias for where to store the results */
9347 tmpt = t.dp + (2*ix + 1);
9349 for (iy = ix + 1; iy < pa; iy++) {
9350 /* first calculate the product */
9351 r = ((tma_mp_word)tmpx) * ((tma_mp_word)a->dp[iy]);
9353 /* now calculate the double precision result, note we use
9354 * addition instead of *2 since it's easier to optimize
9356 r = ((tma_mp_word) *tmpt) + r + r + ((tma_mp_word) u);
9358 /* store lower part */
9359 *tmpt++ = (tma_mp_digit) (r & ((tma_mp_word) MP_MASK));
9362 u = (tma_mp_digit)(r >> ((tma_mp_word) DIGIT_BIT));
9364 /* propagate upwards */
9365 while (u != ((tma_mp_digit) 0)) {
9366 r = ((tma_mp_word) *tmpt) + ((tma_mp_word) u);
9367 *tmpt++ = (tma_mp_digit) (r & ((tma_mp_word) MP_MASK));
9368 u = (tma_mp_digit)(r >> ((tma_mp_word) DIGIT_BIT));
9373 tma_mp_exch (&t, b);
9383 /* End: bn_s_tma_mp_sqr.c */
9385 /* Start: bn_s_tma_mp_sub.c */
9387 #ifdef BN_S_MP_SUB_C
9388 /* LibTomMath, multiple-precision integer library -- Tom St Denis
9390 * LibTomMath is a library that provides multiple-precision
9391 * integer arithmetic as well as number theoretic functionality.
9393 * The library was designed directly after the MPI library by
9394 * Michael Fromberger but has been written from scratch with
9395 * additional optimizations in place.
9397 * The library is free for all purposes without any express
9398 * guarantee it works.
9400 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
9403 /* low level subtraction (assumes |a| > |b|), HAC pp.595 Algorithm 14.9 */
9405 s_tma_mp_sub (tma_mp_int * a, tma_mp_int * b, tma_mp_int * c)
9407 int olduse, res, min, max;
9414 if (c->alloc < max) {
9415 if ((res = tma_mp_grow (c, max)) != MP_OKAY) {
9423 register tma_mp_digit u, *tmpa, *tmpb, *tmpc;
9426 /* alias for digit pointers */
9431 /* set carry to zero */
9433 for (i = 0; i < min; i++) {
9434 /* T[i] = A[i] - B[i] - U */
9435 *tmpc = *tmpa++ - *tmpb++ - u;
9437 /* U = carry bit of T[i]
9438 * Note this saves performing an AND operation since
9439 * if a carry does occur it will propagate all the way to the
9440 * MSB. As a result a single shift is enough to get the carry
9442 u = *tmpc >> ((tma_mp_digit)(CHAR_BIT * sizeof (tma_mp_digit) - 1));
9444 /* Clear carry from T[i] */
9448 /* now copy higher words if any, e.g. if A has more digits than B */
9449 for (; i < max; i++) {
9450 /* T[i] = A[i] - U */
9451 *tmpc = *tmpa++ - u;
9453 /* U = carry bit of T[i] */
9454 u = *tmpc >> ((tma_mp_digit)(CHAR_BIT * sizeof (tma_mp_digit) - 1));
9456 /* Clear carry from T[i] */
9460 /* clear digits above used (since we may not have grown result above) */
9461 for (i = c->used; i < olduse; i++) {
9476 /* End: bn_s_tma_mp_sub.c */
9478 /* Start: bncore.c */
9481 /* LibTomMath, multiple-precision integer library -- Tom St Denis
9483 * LibTomMath is a library that provides multiple-precision
9484 * integer arithmetic as well as number theoretic functionality.
9486 * The library was designed directly after the MPI library by
9487 * Michael Fromberger but has been written from scratch with
9488 * additional optimizations in place.
9490 * The library is free for all purposes without any express
9491 * guarantee it works.
9493 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
9496 /* Known optimal configurations
9498 CPU /Compiler /MUL CUTOFF/SQR CUTOFF
9499 -------------------------------------------------------------
9500 Intel P4 Northwood /GCC v3.4.1 / 88/ 128/LTM 0.32 ;-)
9501 AMD Athlon64 /GCC v3.4.4 / 80/ 120/LTM 0.35
9505 int KARATSUBA_MUL_CUTOFF = 80, /* Min. number of digits before Karatsuba multiplication is used. */
9506 KARATSUBA_SQR_CUTOFF = 120, /* Min. number of digits before Karatsuba squaring is used. */
9508 TOOM_MUL_CUTOFF = 350, /* no optimal values of these are known yet so set em high */
9509 TOOM_SQR_CUTOFF = 400;