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);
432 /* clear the carry */
434 for (ix = 0; ix < pa; ix++) {
437 tma_mp_digit *tmpx, *tmpy;
439 /* get offsets into the two bignums */
440 ty = MIN(b->used-1, ix);
443 /* setup temp aliases */
447 /* this is the number of times the loop will iterrate, essentially
448 while (tx++ < a->used && ty-- >= 0) { ... }
450 iy = MIN(a->used-tx, ty+1);
453 for (iz = 0; iz < iy; ++iz) {
454 _W += ((tma_mp_word)*tmpx++)*((tma_mp_word)*tmpy--);
459 W[ix] = ((tma_mp_digit)_W) & MP_MASK;
461 /* make next carry */
462 _W = _W >> ((tma_mp_word)DIGIT_BIT);
470 register tma_mp_digit *tmpc;
472 for (ix = 0; ix < pa+1; ix++) {
473 /* now extract the previous digit [below the carry] */
477 /* clear unused digits [that existed in the old copy of c] */
478 for (; ix < olduse; ix++) {
491 /* End: bn_fast_s_tma_mp_mul_digs.c */
493 /* Start: bn_fast_s_tma_mp_mul_high_digs.c */
495 #ifdef BN_FAST_S_MP_MUL_HIGH_DIGS_C
496 /* LibTomMath, multiple-precision integer library -- Tom St Denis
498 * LibTomMath is a library that provides multiple-precision
499 * integer arithmetic as well as number theoretic functionality.
501 * The library was designed directly after the MPI library by
502 * Michael Fromberger but has been written from scratch with
503 * additional optimizations in place.
505 * The library is free for all purposes without any express
506 * guarantee it works.
508 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
511 /* this is a modified version of fast_s_mul_digs that only produces
512 * output digits *above* digs. See the comments for fast_s_mul_digs
513 * to see how it works.
515 * This is used in the Barrett reduction since for one of the multiplications
516 * only the higher digits were needed. This essentially halves the work.
518 * Based on Algorithm 14.12 on pp.595 of HAC.
520 int fast_s_tma_mp_mul_high_digs (tma_mp_int * a, tma_mp_int * b, tma_mp_int * c, int digs)
522 int olduse, res, pa, ix, iz;
523 tma_mp_digit W[MP_WARRAY];
526 /* grow the destination as required */
527 pa = a->used + b->used;
529 if ((res = tma_mp_grow (c, pa)) != MP_OKAY) {
534 /* number of output digits to produce */
535 pa = a->used + b->used;
540 for (ix = digs; ix < pa; ix++) {
542 tma_mp_digit *tmpx, *tmpy;
544 /* get offsets into the two bignums */
545 ty = MIN(b->used-1, ix);
548 /* setup temp aliases */
552 /* this is the number of times the loop will iterrate, essentially its
553 while (tx++ < a->used && ty-- >= 0) { ... }
555 iy = MIN(a->used-tx, ty+1);
558 for (iz = 0; iz < iy; iz++) {
559 _W += ((tma_mp_word)*tmpx++)*((tma_mp_word)*tmpy--);
563 W[ix] = ((tma_mp_digit)_W) & MP_MASK;
565 /* make next carry */
566 _W = _W >> ((tma_mp_word)DIGIT_BIT);
574 register tma_mp_digit *tmpc;
577 for (ix = digs; ix < pa; ix++) {
578 /* now extract the previous digit [below the carry] */
582 /* clear unused digits [that existed in the old copy of c] */
583 for (; ix < olduse; ix++) {
596 /* End: bn_fast_s_tma_mp_mul_high_digs.c */
598 /* Start: bn_fast_s_tma_mp_sqr.c */
600 #ifdef BN_FAST_S_MP_SQR_C
601 /* LibTomMath, multiple-precision integer library -- Tom St Denis
603 * LibTomMath is a library that provides multiple-precision
604 * integer arithmetic as well as number theoretic functionality.
606 * The library was designed directly after the MPI library by
607 * Michael Fromberger but has been written from scratch with
608 * additional optimizations in place.
610 * The library is free for all purposes without any express
611 * guarantee it works.
613 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
616 /* the jist of squaring...
617 * you do like mult except the offset of the tmpx [one that
618 * starts closer to zero] can't equal the offset of tmpy.
619 * So basically you set up iy like before then you min it with
620 * (ty-tx) so that it never happens. You double all those
621 * you add in the inner loop
623 After that loop you do the squares and add them in.
626 int fast_s_tma_mp_sqr (tma_mp_int * a, tma_mp_int * b)
628 int olduse, res, pa, ix, iz;
629 tma_mp_digit W[MP_WARRAY], *tmpx;
632 /* grow the destination as required */
633 pa = a->used + a->used;
635 if ((res = tma_mp_grow (b, pa)) != MP_OKAY) {
640 /* number of output digits to produce */
642 for (ix = 0; ix < pa; ix++) {
650 /* get offsets into the two bignums */
651 ty = MIN(a->used-1, ix);
654 /* setup temp aliases */
658 /* this is the number of times the loop will iterrate, essentially
659 while (tx++ < a->used && ty-- >= 0) { ... }
661 iy = MIN(a->used-tx, ty+1);
663 /* now for squaring tx can never equal ty
664 * we halve the distance since they approach at a rate of 2x
665 * and we have to round because odd cases need to be executed
667 iy = MIN(iy, (ty-tx+1)>>1);
670 for (iz = 0; iz < iy; iz++) {
671 _W += ((tma_mp_word)*tmpx++)*((tma_mp_word)*tmpy--);
674 /* double the inner product and add carry */
677 /* even columns have the square term in them */
679 _W += ((tma_mp_word)a->dp[ix>>1])*((tma_mp_word)a->dp[ix>>1]);
683 W[ix] = (tma_mp_digit)(_W & MP_MASK);
685 /* make next carry */
686 W1 = _W >> ((tma_mp_word)DIGIT_BIT);
691 b->used = a->used+a->used;
696 for (ix = 0; ix < pa; ix++) {
697 *tmpb++ = W[ix] & MP_MASK;
700 /* clear unused digits [that existed in the old copy of c] */
701 for (; ix < olduse; ix++) {
714 /* End: bn_fast_s_tma_mp_sqr.c */
716 /* Start: bn_tma_mp_2expt.c */
719 /* LibTomMath, multiple-precision integer library -- Tom St Denis
721 * LibTomMath is a library that provides multiple-precision
722 * integer arithmetic as well as number theoretic functionality.
724 * The library was designed directly after the MPI library by
725 * Michael Fromberger but has been written from scratch with
726 * additional optimizations in place.
728 * The library is free for all purposes without any express
729 * guarantee it works.
731 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
736 * Simple algorithm which zeroes the int, grows it then just sets one bit
740 tma_mp_2expt (tma_mp_int * a, int b)
744 /* zero a as per default */
747 /* grow a to accomodate the single bit */
748 if ((res = tma_mp_grow (a, b / DIGIT_BIT + 1)) != MP_OKAY) {
752 /* set the used count of where the bit will go */
753 a->used = b / DIGIT_BIT + 1;
755 /* put the single bit in its place */
756 a->dp[b / DIGIT_BIT] = ((tma_mp_digit)1) << (b % DIGIT_BIT);
766 /* End: bn_tma_mp_2expt.c */
768 /* Start: bn_tma_mp_abs.c */
771 /* LibTomMath, multiple-precision integer library -- Tom St Denis
773 * LibTomMath is a library that provides multiple-precision
774 * integer arithmetic as well as number theoretic functionality.
776 * The library was designed directly after the MPI library by
777 * Michael Fromberger but has been written from scratch with
778 * additional optimizations in place.
780 * The library is free for all purposes without any express
781 * guarantee it works.
783 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
788 * Simple function copies the input and fixes the sign to positive
791 tma_mp_abs (tma_mp_int * a, tma_mp_int * b)
797 if ((res = tma_mp_copy (a, b)) != MP_OKAY) {
802 /* force the sign of b to positive */
813 /* End: bn_tma_mp_abs.c */
815 /* Start: bn_tma_mp_add.c */
818 /* LibTomMath, multiple-precision integer library -- Tom St Denis
820 * LibTomMath is a library that provides multiple-precision
821 * integer arithmetic as well as number theoretic functionality.
823 * The library was designed directly after the MPI library by
824 * Michael Fromberger but has been written from scratch with
825 * additional optimizations in place.
827 * The library is free for all purposes without any express
828 * guarantee it works.
830 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
833 /* high level addition (handles signs) */
834 int tma_mp_add (tma_mp_int * a, tma_mp_int * b, tma_mp_int * c)
838 /* get sign of both inputs */
842 /* handle two cases, not four */
844 /* both positive or both negative */
845 /* add their magnitudes, copy the sign */
847 res = s_tma_mp_add (a, b, c);
849 /* one positive, the other negative */
850 /* subtract the one with the greater magnitude from */
851 /* the one of the lesser magnitude. The result gets */
852 /* the sign of the one with the greater magnitude. */
853 if (tma_mp_cmp_mag (a, b) == MP_LT) {
855 res = s_tma_mp_sub (b, a, c);
858 res = s_tma_mp_sub (a, b, c);
870 /* End: bn_tma_mp_add.c */
872 /* Start: bn_tma_mp_add_d.c */
875 /* LibTomMath, multiple-precision integer library -- Tom St Denis
877 * LibTomMath is a library that provides multiple-precision
878 * integer arithmetic as well as number theoretic functionality.
880 * The library was designed directly after the MPI library by
881 * Michael Fromberger but has been written from scratch with
882 * additional optimizations in place.
884 * The library is free for all purposes without any express
885 * guarantee it works.
887 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
890 /* single digit addition */
892 tma_mp_add_d (tma_mp_int * a, tma_mp_digit b, tma_mp_int * c)
894 int res, ix, oldused;
895 tma_mp_digit *tmpa, *tmpc, mu;
897 /* grow c as required */
898 if (c->alloc < a->used + 1) {
899 if ((res = tma_mp_grow(c, a->used + 1)) != MP_OKAY) {
904 /* if a is negative and |a| >= b, call c = |a| - b */
905 if (a->sign == MP_NEG && (a->used > 1 || a->dp[0] >= b)) {
906 /* temporarily fix sign of a */
910 res = tma_mp_sub_d(a, b, c);
913 a->sign = c->sign = MP_NEG;
921 /* old number of used digits in c */
924 /* sign always positive */
930 /* destination alias */
933 /* if a is positive */
934 if (a->sign == MP_ZPOS) {
935 /* add digit, after this we're propagating
939 mu = *tmpc >> DIGIT_BIT;
942 /* now handle rest of the digits */
943 for (ix = 1; ix < a->used; ix++) {
944 *tmpc = *tmpa++ + mu;
945 mu = *tmpc >> DIGIT_BIT;
948 /* set final carry */
953 c->used = a->used + 1;
955 /* a was negative and |a| < b */
958 /* the result is a single digit */
960 *tmpc++ = b - a->dp[0];
965 /* setup count so the clearing of oldused
966 * can fall through correctly
971 /* now zero to oldused */
972 while (ix++ < oldused) {
986 /* End: bn_tma_mp_add_d.c */
988 /* Start: bn_tma_mp_addmod.c */
990 #ifdef BN_MP_ADDMOD_C
991 /* LibTomMath, multiple-precision integer library -- Tom St Denis
993 * LibTomMath is a library that provides multiple-precision
994 * integer arithmetic as well as number theoretic functionality.
996 * The library was designed directly after the MPI library by
997 * Michael Fromberger but has been written from scratch with
998 * additional optimizations in place.
1000 * The library is free for all purposes without any express
1001 * guarantee it works.
1003 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
1006 /* d = a + b (mod c) */
1008 tma_mp_addmod (tma_mp_int * a, tma_mp_int * b, tma_mp_int * c, tma_mp_int * d)
1013 if ((res = tma_mp_init (&t)) != MP_OKAY) {
1017 if ((res = tma_mp_add (a, b, &t)) != MP_OKAY) {
1021 res = tma_mp_mod (&t, c, d);
1031 /* End: bn_tma_mp_addmod.c */
1033 /* Start: bn_tma_mp_and.c */
1036 /* LibTomMath, multiple-precision integer library -- Tom St Denis
1038 * LibTomMath is a library that provides multiple-precision
1039 * integer arithmetic as well as number theoretic functionality.
1041 * The library was designed directly after the MPI library by
1042 * Michael Fromberger but has been written from scratch with
1043 * additional optimizations in place.
1045 * The library is free for all purposes without any express
1046 * guarantee it works.
1048 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
1051 /* AND two ints together */
1053 tma_mp_and (tma_mp_int * a, tma_mp_int * b, tma_mp_int * c)
1058 if (a->used > b->used) {
1059 if ((res = tma_mp_init_copy (&t, a)) != MP_OKAY) {
1065 if ((res = tma_mp_init_copy (&t, b)) != MP_OKAY) {
1072 for (ix = 0; ix < px; ix++) {
1073 t.dp[ix] &= x->dp[ix];
1076 /* zero digits above the last from the smallest tma_mp_int */
1077 for (; ix < t.used; ix++) {
1082 tma_mp_exch (c, &t);
1092 /* End: bn_tma_mp_and.c */
1094 /* Start: bn_tma_mp_clamp.c */
1096 #ifdef BN_MP_CLAMP_C
1097 /* LibTomMath, multiple-precision integer library -- Tom St Denis
1099 * LibTomMath is a library that provides multiple-precision
1100 * integer arithmetic as well as number theoretic functionality.
1102 * The library was designed directly after the MPI library by
1103 * Michael Fromberger but has been written from scratch with
1104 * additional optimizations in place.
1106 * The library is free for all purposes without any express
1107 * guarantee it works.
1109 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
1112 /* trim unused digits
1114 * This is used to ensure that leading zero digits are
1115 * trimed and the leading "used" digit will be non-zero
1116 * Typically very fast. Also fixes the sign if there
1117 * are no more leading digits
1120 tma_mp_clamp (tma_mp_int * a)
1122 /* decrease used while the most significant digit is
1125 while (a->used > 0 && a->dp[a->used - 1] == 0) {
1129 /* reset the sign flag if used == 0 */
1140 /* End: bn_tma_mp_clamp.c */
1142 /* Start: bn_tma_mp_clear.c */
1144 #ifdef BN_MP_CLEAR_C
1145 /* LibTomMath, multiple-precision integer library -- Tom St Denis
1147 * LibTomMath is a library that provides multiple-precision
1148 * integer arithmetic as well as number theoretic functionality.
1150 * The library was designed directly after the MPI library by
1151 * Michael Fromberger but has been written from scratch with
1152 * additional optimizations in place.
1154 * The library is free for all purposes without any express
1155 * guarantee it works.
1157 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
1160 /* clear one (frees) */
1162 tma_mp_clear (tma_mp_int * a)
1166 /* only do anything if a hasn't been freed previously */
1167 if (a->dp != NULL) {
1168 /* first zero the digits */
1169 for (i = 0; i < a->used; i++) {
1176 /* reset members to make debugging easier */
1178 a->alloc = a->used = 0;
1188 /* End: bn_tma_mp_clear.c */
1190 /* Start: bn_tma_mp_clear_multi.c */
1192 #ifdef BN_MP_CLEAR_MULTI_C
1193 /* LibTomMath, multiple-precision integer library -- Tom St Denis
1195 * LibTomMath is a library that provides multiple-precision
1196 * integer arithmetic as well as number theoretic functionality.
1198 * The library was designed directly after the MPI library by
1199 * Michael Fromberger but has been written from scratch with
1200 * additional optimizations in place.
1202 * The library is free for all purposes without any express
1203 * guarantee it works.
1205 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
1209 void tma_mp_clear_multi(tma_mp_int *mp, ...)
1211 tma_mp_int* next_mp = mp;
1214 while (next_mp != NULL) {
1215 tma_mp_clear(next_mp);
1216 next_mp = va_arg(args, tma_mp_int*);
1226 /* End: bn_tma_mp_clear_multi.c */
1228 /* Start: bn_tma_mp_cmp.c */
1231 /* LibTomMath, multiple-precision integer library -- Tom St Denis
1233 * LibTomMath is a library that provides multiple-precision
1234 * integer arithmetic as well as number theoretic functionality.
1236 * The library was designed directly after the MPI library by
1237 * Michael Fromberger but has been written from scratch with
1238 * additional optimizations in place.
1240 * The library is free for all purposes without any express
1241 * guarantee it works.
1243 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
1246 /* compare two ints (signed)*/
1248 tma_mp_cmp (tma_mp_int * a, tma_mp_int * b)
1250 /* compare based on sign */
1251 if (a->sign != b->sign) {
1252 if (a->sign == MP_NEG) {
1259 /* compare digits */
1260 if (a->sign == MP_NEG) {
1261 /* if negative compare opposite direction */
1262 return tma_mp_cmp_mag(b, a);
1264 return tma_mp_cmp_mag(a, b);
1273 /* End: bn_tma_mp_cmp.c */
1275 /* Start: bn_tma_mp_cmp_d.c */
1277 #ifdef BN_MP_CMP_D_C
1278 /* LibTomMath, multiple-precision integer library -- Tom St Denis
1280 * LibTomMath is a library that provides multiple-precision
1281 * integer arithmetic as well as number theoretic functionality.
1283 * The library was designed directly after the MPI library by
1284 * Michael Fromberger but has been written from scratch with
1285 * additional optimizations in place.
1287 * The library is free for all purposes without any express
1288 * guarantee it works.
1290 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
1293 /* compare a digit */
1294 int tma_mp_cmp_d(tma_mp_int * a, tma_mp_digit b)
1296 /* compare based on sign */
1297 if (a->sign == MP_NEG) {
1301 /* compare based on magnitude */
1306 /* compare the only digit of a to b */
1309 } else if (a->dp[0] < b) {
1321 /* End: bn_tma_mp_cmp_d.c */
1323 /* Start: bn_tma_mp_cmp_mag.c */
1325 #ifdef BN_MP_CMP_MAG_C
1326 /* LibTomMath, multiple-precision integer library -- Tom St Denis
1328 * LibTomMath is a library that provides multiple-precision
1329 * integer arithmetic as well as number theoretic functionality.
1331 * The library was designed directly after the MPI library by
1332 * Michael Fromberger but has been written from scratch with
1333 * additional optimizations in place.
1335 * The library is free for all purposes without any express
1336 * guarantee it works.
1338 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
1341 /* compare maginitude of two ints (unsigned) */
1342 int tma_mp_cmp_mag (tma_mp_int * a, tma_mp_int * b)
1345 tma_mp_digit *tmpa, *tmpb;
1347 /* compare based on # of non-zero digits */
1348 if (a->used > b->used) {
1352 if (a->used < b->used) {
1357 tmpa = a->dp + (a->used - 1);
1360 tmpb = b->dp + (a->used - 1);
1362 /* compare based on digits */
1363 for (n = 0; n < a->used; ++n, --tmpa, --tmpb) {
1364 if (*tmpa > *tmpb) {
1368 if (*tmpa < *tmpb) {
1380 /* End: bn_tma_mp_cmp_mag.c */
1382 /* Start: bn_tma_mp_cnt_lsb.c */
1384 #ifdef BN_MP_CNT_LSB_C
1385 /* LibTomMath, multiple-precision integer library -- Tom St Denis
1387 * LibTomMath is a library that provides multiple-precision
1388 * integer arithmetic as well as number theoretic functionality.
1390 * The library was designed directly after the MPI library by
1391 * Michael Fromberger but has been written from scratch with
1392 * additional optimizations in place.
1394 * The library is free for all purposes without any express
1395 * guarantee it works.
1397 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
1400 static const int lnz[16] = {
1401 4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0
1404 /* Counts the number of lsbs which are zero before the first zero bit */
1405 int tma_mp_cnt_lsb(tma_mp_int *a)
1411 if (tma_mp_iszero(a) == 1) {
1415 /* scan lower digits until non-zero */
1416 for (x = 0; x < a->used && a->dp[x] == 0; x++);
1420 /* now scan this digit until a 1 is found */
1437 /* End: bn_tma_mp_cnt_lsb.c */
1439 /* Start: bn_tma_mp_copy.c */
1442 /* LibTomMath, multiple-precision integer library -- Tom St Denis
1444 * LibTomMath is a library that provides multiple-precision
1445 * integer arithmetic as well as number theoretic functionality.
1447 * The library was designed directly after the MPI library by
1448 * Michael Fromberger but has been written from scratch with
1449 * additional optimizations in place.
1451 * The library is free for all purposes without any express
1452 * guarantee it works.
1454 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
1459 tma_mp_copy (tma_mp_int * a, tma_mp_int * b)
1463 /* if dst == src do nothing */
1469 if (b->alloc < a->used) {
1470 if ((res = tma_mp_grow (b, a->used)) != MP_OKAY) {
1475 /* zero b and copy the parameters over */
1477 register tma_mp_digit *tmpa, *tmpb;
1479 /* pointer aliases */
1487 /* copy all the digits */
1488 for (n = 0; n < a->used; n++) {
1492 /* clear high digits */
1493 for (; n < b->used; n++) {
1498 /* copy used count and sign */
1509 /* End: bn_tma_mp_copy.c */
1511 /* Start: bn_tma_mp_count_bits.c */
1513 #ifdef BN_MP_COUNT_BITS_C
1514 /* LibTomMath, multiple-precision integer library -- Tom St Denis
1516 * LibTomMath is a library that provides multiple-precision
1517 * integer arithmetic as well as number theoretic functionality.
1519 * The library was designed directly after the MPI library by
1520 * Michael Fromberger but has been written from scratch with
1521 * additional optimizations in place.
1523 * The library is free for all purposes without any express
1524 * guarantee it works.
1526 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
1529 /* returns the number of bits in an int */
1531 tma_mp_count_bits (tma_mp_int * a)
1541 /* get number of digits and add that */
1542 r = (a->used - 1) * DIGIT_BIT;
1544 /* take the last digit and count the bits in it */
1545 q = a->dp[a->used - 1];
1546 while (q > ((tma_mp_digit) 0)) {
1548 q >>= ((tma_mp_digit) 1);
1558 /* End: bn_tma_mp_count_bits.c */
1560 /* Start: bn_tma_mp_div.c */
1563 /* LibTomMath, multiple-precision integer library -- Tom St Denis
1565 * LibTomMath is a library that provides multiple-precision
1566 * integer arithmetic as well as number theoretic functionality.
1568 * The library was designed directly after the MPI library by
1569 * Michael Fromberger but has been written from scratch with
1570 * additional optimizations in place.
1572 * The library is free for all purposes without any express
1573 * guarantee it works.
1575 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
1578 #ifdef BN_MP_DIV_SMALL
1580 /* slower bit-bang division... also smaller */
1581 int tma_mp_div(tma_mp_int * a, tma_mp_int * b, tma_mp_int * c, tma_mp_int * d)
1583 tma_mp_int ta, tb, tq, q;
1586 /* is divisor zero ? */
1587 if (tma_mp_iszero (b) == 1) {
1591 /* if a < b then q=0, r = a */
1592 if (tma_mp_cmp_mag (a, b) == MP_LT) {
1594 res = tma_mp_copy (a, d);
1604 /* init our temps */
1605 if ((res = tma_mp_init_multi(&ta, &tb, &tq, &q, NULL) != MP_OKAY)) {
1611 n = tma_mp_count_bits(a) - tma_mp_count_bits(b);
1612 if (((res = tma_mp_abs(a, &ta)) != MP_OKAY) ||
1613 ((res = tma_mp_abs(b, &tb)) != MP_OKAY) ||
1614 ((res = tma_mp_mul_2d(&tb, n, &tb)) != MP_OKAY) ||
1615 ((res = tma_mp_mul_2d(&tq, n, &tq)) != MP_OKAY)) {
1620 if (tma_mp_cmp(&tb, &ta) != MP_GT) {
1621 if (((res = tma_mp_sub(&ta, &tb, &ta)) != MP_OKAY) ||
1622 ((res = tma_mp_add(&q, &tq, &q)) != MP_OKAY)) {
1626 if (((res = tma_mp_div_2d(&tb, 1, &tb, NULL)) != MP_OKAY) ||
1627 ((res = tma_mp_div_2d(&tq, 1, &tq, NULL)) != MP_OKAY)) {
1632 /* now q == quotient and ta == remainder */
1634 n2 = (a->sign == b->sign ? MP_ZPOS : MP_NEG);
1637 c->sign = (tma_mp_iszero(c) == MP_YES) ? MP_ZPOS : n2;
1640 tma_mp_exch(d, &ta);
1641 d->sign = (tma_mp_iszero(d) == MP_YES) ? MP_ZPOS : n;
1644 tma_mp_clear_multi(&ta, &tb, &tq, &q, NULL);
1650 /* integer signed division.
1651 * c*b + d == a [e.g. a/b, c=quotient, d=remainder]
1652 * HAC pp.598 Algorithm 14.20
1654 * Note that the description in HAC is horribly
1655 * incomplete. For example, it doesn't consider
1656 * the case where digits are removed from 'x' in
1657 * the inner loop. It also doesn't consider the
1658 * case that y has fewer than three digits, etc..
1660 * The overall algorithm is as described as
1661 * 14.20 from HAC but fixed to treat these cases.
1663 int tma_mp_div (tma_mp_int * a, tma_mp_int * b, tma_mp_int * c, tma_mp_int * d)
1665 tma_mp_int q, x, y, t1, t2;
1666 int res, n, t, i, norm, neg;
1668 /* is divisor zero ? */
1669 if (tma_mp_iszero (b) == 1) {
1673 /* if a < b then q=0, r = a */
1674 if (tma_mp_cmp_mag (a, b) == MP_LT) {
1676 res = tma_mp_copy (a, d);
1686 if ((res = tma_mp_init_size (&q, a->used + 2)) != MP_OKAY) {
1689 q.used = a->used + 2;
1691 if ((res = tma_mp_init (&t1)) != MP_OKAY) {
1695 if ((res = tma_mp_init (&t2)) != MP_OKAY) {
1699 if ((res = tma_mp_init_copy (&x, a)) != MP_OKAY) {
1703 if ((res = tma_mp_init_copy (&y, b)) != MP_OKAY) {
1708 neg = (a->sign == b->sign) ? MP_ZPOS : MP_NEG;
1709 x.sign = y.sign = MP_ZPOS;
1711 /* normalize both x and y, ensure that y >= b/2, [b == 2**DIGIT_BIT] */
1712 norm = tma_mp_count_bits(&y) % DIGIT_BIT;
1713 if (norm < (int)(DIGIT_BIT-1)) {
1714 norm = (DIGIT_BIT-1) - norm;
1715 if ((res = tma_mp_mul_2d (&x, norm, &x)) != MP_OKAY) {
1718 if ((res = tma_mp_mul_2d (&y, norm, &y)) != MP_OKAY) {
1725 /* note hac does 0 based, so if used==5 then its 0,1,2,3,4, e.g. use 4 */
1729 /* while (x >= y*b**n-t) do { q[n-t] += 1; x -= y*b**{n-t} } */
1730 if ((res = tma_mp_lshd (&y, n - t)) != MP_OKAY) { /* y = y*b**{n-t} */
1734 while (tma_mp_cmp (&x, &y) != MP_LT) {
1736 if ((res = tma_mp_sub (&x, &y, &x)) != MP_OKAY) {
1741 /* reset y by shifting it back down */
1742 tma_mp_rshd (&y, n - t);
1744 /* step 3. for i from n down to (t + 1) */
1745 for (i = n; i >= (t + 1); i--) {
1750 /* step 3.1 if xi == yt then set q{i-t-1} to b-1,
1751 * otherwise set q{i-t-1} to (xi*b + x{i-1})/yt */
1752 if (x.dp[i] == y.dp[t]) {
1753 q.dp[i - t - 1] = ((((tma_mp_digit)1) << DIGIT_BIT) - 1);
1756 tmp = ((tma_mp_word) x.dp[i]) << ((tma_mp_word) DIGIT_BIT);
1757 tmp |= ((tma_mp_word) x.dp[i - 1]);
1758 tmp /= ((tma_mp_word) y.dp[t]);
1759 if (tmp > (tma_mp_word) MP_MASK)
1761 q.dp[i - t - 1] = (tma_mp_digit) (tmp & (tma_mp_word) (MP_MASK));
1764 /* while (q{i-t-1} * (yt * b + y{t-1})) >
1765 xi * b**2 + xi-1 * b + xi-2
1769 q.dp[i - t - 1] = (q.dp[i - t - 1] + 1) & MP_MASK;
1771 q.dp[i - t - 1] = (q.dp[i - t - 1] - 1) & MP_MASK;
1773 /* find left hand */
1775 t1.dp[0] = (t - 1 < 0) ? 0 : y.dp[t - 1];
1778 if ((res = tma_mp_mul_d (&t1, q.dp[i - t - 1], &t1)) != MP_OKAY) {
1782 /* find right hand */
1783 t2.dp[0] = (i - 2 < 0) ? 0 : x.dp[i - 2];
1784 t2.dp[1] = (i - 1 < 0) ? 0 : x.dp[i - 1];
1787 } while (tma_mp_cmp_mag(&t1, &t2) == MP_GT);
1789 /* step 3.3 x = x - q{i-t-1} * y * b**{i-t-1} */
1790 if ((res = tma_mp_mul_d (&y, q.dp[i - t - 1], &t1)) != MP_OKAY) {
1794 if ((res = tma_mp_lshd (&t1, i - t - 1)) != MP_OKAY) {
1798 if ((res = tma_mp_sub (&x, &t1, &x)) != MP_OKAY) {
1802 /* if x < 0 then { x = x + y*b**{i-t-1}; q{i-t-1} -= 1; } */
1803 if (x.sign == MP_NEG) {
1804 if ((res = tma_mp_copy (&y, &t1)) != MP_OKAY) {
1807 if ((res = tma_mp_lshd (&t1, i - t - 1)) != MP_OKAY) {
1810 if ((res = tma_mp_add (&x, &t1, &x)) != MP_OKAY) {
1814 q.dp[i - t - 1] = (q.dp[i - t - 1] - 1UL) & MP_MASK;
1818 /* now q is the quotient and x is the remainder
1819 * [which we have to normalize]
1822 /* get sign before writing to c */
1823 x.sign = x.used == 0 ? MP_ZPOS : a->sign;
1827 tma_mp_exch (&q, c);
1832 tma_mp_div_2d (&x, norm, &x, NULL);
1833 tma_mp_exch (&x, d);
1838 LBL_Y:tma_mp_clear (&y);
1839 LBL_X:tma_mp_clear (&x);
1840 LBL_T2:tma_mp_clear (&t2);
1841 LBL_T1:tma_mp_clear (&t1);
1842 LBL_Q:tma_mp_clear (&q);
1854 /* End: bn_tma_mp_div.c */
1856 /* Start: bn_tma_mp_div_2.c */
1858 #ifdef BN_MP_DIV_2_C
1859 /* LibTomMath, multiple-precision integer library -- Tom St Denis
1861 * LibTomMath is a library that provides multiple-precision
1862 * integer arithmetic as well as number theoretic functionality.
1864 * The library was designed directly after the MPI library by
1865 * Michael Fromberger but has been written from scratch with
1866 * additional optimizations in place.
1868 * The library is free for all purposes without any express
1869 * guarantee it works.
1871 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
1875 int tma_mp_div_2(tma_mp_int * a, tma_mp_int * b)
1877 int x, res, oldused;
1880 if (b->alloc < a->used) {
1881 if ((res = tma_mp_grow (b, a->used)) != MP_OKAY) {
1889 register tma_mp_digit r, rr, *tmpa, *tmpb;
1892 tmpa = a->dp + b->used - 1;
1895 tmpb = b->dp + b->used - 1;
1899 for (x = b->used - 1; x >= 0; x--) {
1900 /* get the carry for the next iteration */
1903 /* shift the current digit, add in carry and store */
1904 *tmpb-- = (*tmpa-- >> 1) | (r << (DIGIT_BIT - 1));
1906 /* forward carry to next iteration */
1910 /* zero excess digits */
1911 tmpb = b->dp + b->used;
1912 for (x = b->used; x < oldused; x++) {
1926 /* End: bn_tma_mp_div_2.c */
1928 /* Start: bn_tma_mp_div_2d.c */
1930 #ifdef BN_MP_DIV_2D_C
1931 /* LibTomMath, multiple-precision integer library -- Tom St Denis
1933 * LibTomMath is a library that provides multiple-precision
1934 * integer arithmetic as well as number theoretic functionality.
1936 * The library was designed directly after the MPI library by
1937 * Michael Fromberger but has been written from scratch with
1938 * additional optimizations in place.
1940 * The library is free for all purposes without any express
1941 * guarantee it works.
1943 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
1946 /* shift right by a certain bit count (store quotient in c, optional remainder in d) */
1947 int tma_mp_div_2d (tma_mp_int * a, int b, tma_mp_int * c, tma_mp_int * d)
1949 tma_mp_digit D, r, rr;
1954 /* if the shift count is <= 0 then we do no work */
1956 res = tma_mp_copy (a, c);
1963 if ((res = tma_mp_init (&t)) != MP_OKAY) {
1967 /* get the remainder */
1969 if ((res = tma_mp_mod_2d (a, b, &t)) != MP_OKAY) {
1976 if ((res = tma_mp_copy (a, c)) != MP_OKAY) {
1981 /* shift by as many digits in the bit count */
1982 if (b >= (int)DIGIT_BIT) {
1983 tma_mp_rshd (c, b / DIGIT_BIT);
1986 /* shift any bit count < DIGIT_BIT */
1987 D = (tma_mp_digit) (b % DIGIT_BIT);
1989 register tma_mp_digit *tmpc, mask, shift;
1992 mask = (((tma_mp_digit)1) << D) - 1;
1995 shift = DIGIT_BIT - D;
1998 tmpc = c->dp + (c->used - 1);
2002 for (x = c->used - 1; x >= 0; x--) {
2003 /* get the lower bits of this word in a temp */
2006 /* shift the current word and mix in the carry bits from the previous word */
2007 *tmpc = (*tmpc >> D) | (r << shift);
2010 /* set the carry to the carry bits of the current word found above */
2016 tma_mp_exch (&t, d);
2027 /* End: bn_tma_mp_div_2d.c */
2029 /* Start: bn_tma_mp_div_3.c */
2031 #ifdef BN_MP_DIV_3_C
2032 /* LibTomMath, multiple-precision integer library -- Tom St Denis
2034 * LibTomMath is a library that provides multiple-precision
2035 * integer arithmetic as well as number theoretic functionality.
2037 * The library was designed directly after the MPI library by
2038 * Michael Fromberger but has been written from scratch with
2039 * additional optimizations in place.
2041 * The library is free for all purposes without any express
2042 * guarantee it works.
2044 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
2047 /* divide by three (based on routine from MPI and the GMP manual) */
2049 tma_mp_div_3 (tma_mp_int * a, tma_mp_int *c, tma_mp_digit * d)
2056 /* b = 2**DIGIT_BIT / 3 */
2057 b = (((tma_mp_word)1) << ((tma_mp_word)DIGIT_BIT)) / ((tma_mp_word)3);
2059 if ((res = tma_mp_init_size(&q, a->used)) != MP_OKAY) {
2066 for (ix = a->used - 1; ix >= 0; ix--) {
2067 w = (w << ((tma_mp_word)DIGIT_BIT)) | ((tma_mp_word)a->dp[ix]);
2070 /* multiply w by [1/3] */
2071 t = (w * ((tma_mp_word)b)) >> ((tma_mp_word)DIGIT_BIT);
2073 /* now subtract 3 * [w/3] from w, to get the remainder */
2076 /* fixup the remainder as required since
2077 * the optimization is not exact.
2086 q.dp[ix] = (tma_mp_digit)t;
2089 /* [optional] store the remainder */
2091 *d = (tma_mp_digit)w;
2094 /* [optional] store the quotient */
2110 /* End: bn_tma_mp_div_3.c */
2112 /* Start: bn_tma_mp_div_d.c */
2114 #ifdef BN_MP_DIV_D_C
2115 /* LibTomMath, multiple-precision integer library -- Tom St Denis
2117 * LibTomMath is a library that provides multiple-precision
2118 * integer arithmetic as well as number theoretic functionality.
2120 * The library was designed directly after the MPI library by
2121 * Michael Fromberger but has been written from scratch with
2122 * additional optimizations in place.
2124 * The library is free for all purposes without any express
2125 * guarantee it works.
2127 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
2130 static int s_is_power_of_two(tma_mp_digit b, int *p)
2134 /* fast return if no power of two */
2135 if ((b==0) || (b & (b-1))) {
2139 for (x = 0; x < DIGIT_BIT; x++) {
2140 if (b == (((tma_mp_digit)1)<<x)) {
2148 /* single digit division (based on routine from MPI) */
2149 int tma_mp_div_d (tma_mp_int * a, tma_mp_digit b, tma_mp_int * c, tma_mp_digit * d)
2156 /* cannot divide by zero */
2162 if (b == 1 || tma_mp_iszero(a) == 1) {
2167 return tma_mp_copy(a, c);
2172 /* power of two ? */
2173 if (s_is_power_of_two(b, &ix) == 1) {
2175 *d = a->dp[0] & ((((tma_mp_digit)1)<<ix) - 1);
2178 return tma_mp_div_2d(a, ix, c, NULL);
2183 #ifdef BN_MP_DIV_3_C
2186 return tma_mp_div_3(a, c, d);
2190 /* no easy answer [c'est la vie]. Just division */
2191 if ((res = tma_mp_init_size(&q, a->used)) != MP_OKAY) {
2198 for (ix = a->used - 1; ix >= 0; ix--) {
2199 w = (w << ((tma_mp_word)DIGIT_BIT)) | ((tma_mp_word)a->dp[ix]);
2202 t = (tma_mp_digit)(w / b);
2203 w -= ((tma_mp_word)t) * ((tma_mp_word)b);
2207 q.dp[ix] = (tma_mp_digit)t;
2211 *d = (tma_mp_digit)w;
2229 /* End: bn_tma_mp_div_d.c */
2231 /* Start: bn_tma_mp_dr_is_modulus.c */
2233 #ifdef BN_MP_DR_IS_MODULUS_C
2234 /* LibTomMath, multiple-precision integer library -- Tom St Denis
2236 * LibTomMath is a library that provides multiple-precision
2237 * integer arithmetic as well as number theoretic functionality.
2239 * The library was designed directly after the MPI library by
2240 * Michael Fromberger but has been written from scratch with
2241 * additional optimizations in place.
2243 * The library is free for all purposes without any express
2244 * guarantee it works.
2246 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
2249 /* determines if a number is a valid DR modulus */
2250 int tma_mp_dr_is_modulus(tma_mp_int *a)
2254 /* must be at least two digits */
2259 /* must be of the form b**k - a [a <= b] so all
2260 * but the first digit must be equal to -1 (mod b).
2262 for (ix = 1; ix < a->used; ix++) {
2263 if (a->dp[ix] != MP_MASK) {
2276 /* End: bn_tma_mp_dr_is_modulus.c */
2278 /* Start: bn_tma_mp_dr_reduce.c */
2280 #ifdef BN_MP_DR_REDUCE_C
2281 /* LibTomMath, multiple-precision integer library -- Tom St Denis
2283 * LibTomMath is a library that provides multiple-precision
2284 * integer arithmetic as well as number theoretic functionality.
2286 * The library was designed directly after the MPI library by
2287 * Michael Fromberger but has been written from scratch with
2288 * additional optimizations in place.
2290 * The library is free for all purposes without any express
2291 * guarantee it works.
2293 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
2296 /* reduce "x" in place modulo "n" using the Diminished Radix algorithm.
2298 * Based on algorithm from the paper
2300 * "Generating Efficient Primes for Discrete Log Cryptosystems"
2301 * Chae Hoon Lim, Pil Joong Lee,
2302 * POSTECH Information Research Laboratories
2304 * The modulus must be of a special format [see manual]
2306 * Has been modified to use algorithm 7.10 from the LTM book instead
2308 * Input x must be in the range 0 <= x <= (n-1)**2
2311 tma_mp_dr_reduce (tma_mp_int * x, tma_mp_int * n, tma_mp_digit k)
2315 tma_mp_digit mu, *tmpx1, *tmpx2;
2317 /* m = digits in modulus */
2320 /* ensure that "x" has at least 2m digits */
2321 if (x->alloc < m + m) {
2322 if ((err = tma_mp_grow (x, m + m)) != MP_OKAY) {
2327 /* top of loop, this is where the code resumes if
2328 * another reduction pass is required.
2331 /* aliases for digits */
2332 /* alias for lower half of x */
2335 /* alias for upper half of x, or x/B**m */
2338 /* set carry to zero */
2341 /* compute (x mod B**m) + k * [x/B**m] inline and inplace */
2342 for (i = 0; i < m; i++) {
2343 r = ((tma_mp_word)*tmpx2++) * ((tma_mp_word)k) + *tmpx1 + mu;
2344 *tmpx1++ = (tma_mp_digit)(r & MP_MASK);
2345 mu = (tma_mp_digit)(r >> ((tma_mp_word)DIGIT_BIT));
2348 /* set final carry */
2351 /* zero words above m */
2352 for (i = m + 1; i < x->used; i++) {
2356 /* clamp, sub and return */
2359 /* if x >= n then subtract and reduce again
2360 * Each successive "recursion" makes the input smaller and smaller.
2362 if (tma_mp_cmp_mag (x, n) != MP_LT) {
2363 s_tma_mp_sub(x, n, x);
2374 /* End: bn_tma_mp_dr_reduce.c */
2376 /* Start: bn_tma_mp_dr_setup.c */
2378 #ifdef BN_MP_DR_SETUP_C
2379 /* LibTomMath, multiple-precision integer library -- Tom St Denis
2381 * LibTomMath is a library that provides multiple-precision
2382 * integer arithmetic as well as number theoretic functionality.
2384 * The library was designed directly after the MPI library by
2385 * Michael Fromberger but has been written from scratch with
2386 * additional optimizations in place.
2388 * The library is free for all purposes without any express
2389 * guarantee it works.
2391 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
2394 /* determines the setup value */
2395 void tma_mp_dr_setup(tma_mp_int *a, tma_mp_digit *d)
2397 /* the casts are required if DIGIT_BIT is one less than
2398 * the number of bits in a tma_mp_digit [e.g. DIGIT_BIT==31]
2400 *d = (tma_mp_digit)((((tma_mp_word)1) << ((tma_mp_word)DIGIT_BIT)) -
2401 ((tma_mp_word)a->dp[0]));
2410 /* End: bn_tma_mp_dr_setup.c */
2412 /* Start: bn_tma_mp_exch.c */
2415 /* LibTomMath, multiple-precision integer library -- Tom St Denis
2417 * LibTomMath is a library that provides multiple-precision
2418 * integer arithmetic as well as number theoretic functionality.
2420 * The library was designed directly after the MPI library by
2421 * Michael Fromberger but has been written from scratch with
2422 * additional optimizations in place.
2424 * The library is free for all purposes without any express
2425 * guarantee it works.
2427 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
2430 /* swap the elements of two integers, for cases where you can't simply swap the
2431 * tma_mp_int pointers around
2434 tma_mp_exch (tma_mp_int * a, tma_mp_int * b)
2448 /* End: bn_tma_mp_exch.c */
2450 /* Start: bn_tma_mp_expt_d.c */
2452 #ifdef BN_MP_EXPT_D_C
2453 /* LibTomMath, multiple-precision integer library -- Tom St Denis
2455 * LibTomMath is a library that provides multiple-precision
2456 * integer arithmetic as well as number theoretic functionality.
2458 * The library was designed directly after the MPI library by
2459 * Michael Fromberger but has been written from scratch with
2460 * additional optimizations in place.
2462 * The library is free for all purposes without any express
2463 * guarantee it works.
2465 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
2468 /* calculate c = a**b using a square-multiply algorithm */
2469 int tma_mp_expt_d (tma_mp_int * a, tma_mp_digit b, tma_mp_int * c)
2474 if ((res = tma_mp_init_copy (&g, a)) != MP_OKAY) {
2478 /* set initial result */
2481 for (x = 0; x < (int) DIGIT_BIT; x++) {
2483 if ((res = tma_mp_sqr (c, c)) != MP_OKAY) {
2488 /* if the bit is set multiply */
2489 if ((b & (tma_mp_digit) (((tma_mp_digit)1) << (DIGIT_BIT - 1))) != 0) {
2490 if ((res = tma_mp_mul (c, &g, c)) != MP_OKAY) {
2496 /* shift to next bit */
2509 /* End: bn_tma_mp_expt_d.c */
2511 /* Start: bn_tma_mp_exptmod.c */
2513 #ifdef BN_MP_EXPTMOD_C
2514 /* LibTomMath, multiple-precision integer library -- Tom St Denis
2516 * LibTomMath is a library that provides multiple-precision
2517 * integer arithmetic as well as number theoretic functionality.
2519 * The library was designed directly after the MPI library by
2520 * Michael Fromberger but has been written from scratch with
2521 * additional optimizations in place.
2523 * The library is free for all purposes without any express
2524 * guarantee it works.
2526 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
2530 /* this is a shell function that calls either the normal or Montgomery
2531 * exptmod functions. Originally the call to the montgomery code was
2532 * embedded in the normal function but that wasted alot of stack space
2533 * for nothing (since 99% of the time the Montgomery code would be called)
2535 int tma_mp_exptmod (tma_mp_int * G, tma_mp_int * X, tma_mp_int * P, tma_mp_int * Y)
2539 /* modulus P must be positive */
2540 if (P->sign == MP_NEG) {
2544 /* if exponent X is negative we have to recurse */
2545 if (X->sign == MP_NEG) {
2546 #ifdef BN_MP_INVMOD_C
2547 tma_mp_int tmpG, tmpX;
2550 /* first compute 1/G mod P */
2551 if ((err = tma_mp_init(&tmpG)) != MP_OKAY) {
2554 if ((err = tma_mp_invmod(G, P, &tmpG)) != MP_OKAY) {
2555 tma_mp_clear(&tmpG);
2560 if ((err = tma_mp_init(&tmpX)) != MP_OKAY) {
2561 tma_mp_clear(&tmpG);
2564 if ((err = tma_mp_abs(X, &tmpX)) != MP_OKAY) {
2565 tma_mp_clear_multi(&tmpG, &tmpX, NULL);
2569 /* and now compute (1/G)**|X| instead of G**X [X < 0] */
2570 err = tma_mp_exptmod(&tmpG, &tmpX, P, Y);
2571 tma_mp_clear_multi(&tmpG, &tmpX, NULL);
2579 /* modified diminished radix reduction */
2580 #if defined(BN_MP_REDUCE_IS_2K_L_C) && defined(BN_MP_REDUCE_2K_L_C) && defined(BN_S_MP_EXPTMOD_C)
2581 if (tma_mp_reduce_is_2k_l(P) == MP_YES) {
2582 return s_tma_mp_exptmod(G, X, P, Y, 1);
2586 #ifdef BN_MP_DR_IS_MODULUS_C
2587 /* is it a DR modulus? */
2588 dr = tma_mp_dr_is_modulus(P);
2594 #ifdef BN_MP_REDUCE_IS_2K_C
2595 /* if not, is it a unrestricted DR modulus? */
2597 dr = tma_mp_reduce_is_2k(P) << 1;
2601 /* if the modulus is odd or dr != 0 use the montgomery method */
2602 #ifdef BN_MP_EXPTMOD_FAST_C
2603 if (tma_mp_isodd (P) == 1 || dr != 0) {
2604 return tma_mp_exptmod_fast (G, X, P, Y, dr);
2607 #ifdef BN_S_MP_EXPTMOD_C
2608 /* otherwise use the generic Barrett reduction technique */
2609 return s_tma_mp_exptmod (G, X, P, Y, 0);
2611 /* no exptmod for evens */
2614 #ifdef BN_MP_EXPTMOD_FAST_C
2625 /* End: bn_tma_mp_exptmod.c */
2627 /* Start: bn_tma_mp_exptmod_fast.c */
2629 #ifdef BN_MP_EXPTMOD_FAST_C
2630 /* LibTomMath, multiple-precision integer library -- Tom St Denis
2632 * LibTomMath is a library that provides multiple-precision
2633 * integer arithmetic as well as number theoretic functionality.
2635 * The library was designed directly after the MPI library by
2636 * Michael Fromberger but has been written from scratch with
2637 * additional optimizations in place.
2639 * The library is free for all purposes without any express
2640 * guarantee it works.
2642 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
2645 /* computes Y == G**X mod P, HAC pp.616, Algorithm 14.85
2647 * Uses a left-to-right k-ary sliding window to compute the modular exponentiation.
2648 * The value of k changes based on the size of the exponent.
2650 * Uses Montgomery or Diminished Radix reduction [whichever appropriate]
2656 #define TAB_SIZE 256
2659 int tma_mp_exptmod_fast (tma_mp_int * G, tma_mp_int * X, tma_mp_int * P, tma_mp_int * Y, int redmode)
2661 tma_mp_int M[TAB_SIZE], res;
2662 tma_mp_digit buf, mp;
2663 int err, bitbuf, bitcpy, bitcnt, mode, digidx, x, y, winsize;
2665 /* use a pointer to the reduction algorithm. This allows us to use
2666 * one of many reduction algorithms without modding the guts of
2667 * the code with if statements everywhere.
2669 int (*redux)(tma_mp_int*,tma_mp_int*,tma_mp_digit);
2671 /* find window size */
2672 x = tma_mp_count_bits (X);
2675 } else if (x <= 36) {
2677 } else if (x <= 140) {
2679 } else if (x <= 450) {
2681 } else if (x <= 1303) {
2683 } else if (x <= 3529) {
2696 /* init first cell */
2697 if ((err = tma_mp_init(&M[1])) != MP_OKAY) {
2701 /* now init the second half of the array */
2702 for (x = 1<<(winsize-1); x < (1 << winsize); x++) {
2703 if ((err = tma_mp_init(&M[x])) != MP_OKAY) {
2704 for (y = 1<<(winsize-1); y < x; y++) {
2705 tma_mp_clear (&M[y]);
2707 tma_mp_clear(&M[1]);
2712 /* determine and setup reduction code */
2714 #ifdef BN_MP_MONTGOMERY_SETUP_C
2715 /* now setup montgomery */
2716 if ((err = tma_mp_montgomery_setup (P, &mp)) != MP_OKAY) {
2724 /* automatically pick the comba one if available (saves quite a few calls/ifs) */
2725 #ifdef BN_FAST_MP_MONTGOMERY_REDUCE_C
2726 if (((P->used * 2 + 1) < MP_WARRAY) &&
2727 P->used < (1 << ((CHAR_BIT * sizeof (tma_mp_word)) - (2 * DIGIT_BIT)))) {
2728 redux = fast_tma_mp_montgomery_reduce;
2732 #ifdef BN_MP_MONTGOMERY_REDUCE_C
2733 /* use slower baseline Montgomery method */
2734 redux = tma_mp_montgomery_reduce;
2740 } else if (redmode == 1) {
2741 #if defined(BN_MP_DR_SETUP_C) && defined(BN_MP_DR_REDUCE_C)
2742 /* setup DR reduction for moduli of the form B**k - b */
2743 tma_mp_dr_setup(P, &mp);
2744 redux = tma_mp_dr_reduce;
2750 #if defined(BN_MP_REDUCE_2K_SETUP_C) && defined(BN_MP_REDUCE_2K_C)
2751 /* setup DR reduction for moduli of the form 2**k - b */
2752 if ((err = tma_mp_reduce_2k_setup(P, &mp)) != MP_OKAY) {
2755 redux = tma_mp_reduce_2k;
2763 if ((err = tma_mp_init (&res)) != MP_OKAY) {
2771 * The first half of the table is not computed though accept for M[0] and M[1]
2775 #ifdef BN_MP_MONTGOMERY_CALC_NORMALIZATION_C
2776 /* now we need R mod m */
2777 if ((err = tma_mp_montgomery_calc_normalization (&res, P)) != MP_OKAY) {
2785 /* now set M[1] to G * R mod m */
2786 if ((err = tma_mp_mulmod (G, &res, P, &M[1])) != MP_OKAY) {
2790 tma_mp_set(&res, 1);
2791 if ((err = tma_mp_mod(G, P, &M[1])) != MP_OKAY) {
2796 /* compute the value at M[1<<(winsize-1)] by squaring M[1] (winsize-1) times */
2797 if ((err = tma_mp_copy (&M[1], &M[1 << (winsize - 1)])) != MP_OKAY) {
2801 for (x = 0; x < (winsize - 1); x++) {
2802 if ((err = tma_mp_sqr (&M[1 << (winsize - 1)], &M[1 << (winsize - 1)])) != MP_OKAY) {
2805 if ((err = redux (&M[1 << (winsize - 1)], P, mp)) != MP_OKAY) {
2810 /* create upper table */
2811 for (x = (1 << (winsize - 1)) + 1; x < (1 << winsize); x++) {
2812 if ((err = tma_mp_mul (&M[x - 1], &M[1], &M[x])) != MP_OKAY) {
2815 if ((err = redux (&M[x], P, mp)) != MP_OKAY) {
2820 /* set initial mode and bit cnt */
2824 digidx = X->used - 1;
2829 /* grab next digit as required */
2830 if (--bitcnt == 0) {
2831 /* if digidx == -1 we are out of digits so break */
2835 /* read next digit and reset bitcnt */
2836 buf = X->dp[digidx--];
2837 bitcnt = (int)DIGIT_BIT;
2840 /* grab the next msb from the exponent */
2841 y = (tma_mp_digit)(buf >> (DIGIT_BIT - 1)) & 1;
2842 buf <<= (tma_mp_digit)1;
2844 /* if the bit is zero and mode == 0 then we ignore it
2845 * These represent the leading zero bits before the first 1 bit
2846 * in the exponent. Technically this opt is not required but it
2847 * does lower the # of trivial squaring/reductions used
2849 if (mode == 0 && y == 0) {
2853 /* if the bit is zero and mode == 1 then we square */
2854 if (mode == 1 && y == 0) {
2855 if ((err = tma_mp_sqr (&res, &res)) != MP_OKAY) {
2858 if ((err = redux (&res, P, mp)) != MP_OKAY) {
2864 /* else we add it to the window */
2865 bitbuf |= (y << (winsize - ++bitcpy));
2868 if (bitcpy == winsize) {
2869 /* ok window is filled so square as required and multiply */
2871 for (x = 0; x < winsize; x++) {
2872 if ((err = tma_mp_sqr (&res, &res)) != MP_OKAY) {
2875 if ((err = redux (&res, P, mp)) != MP_OKAY) {
2881 if ((err = tma_mp_mul (&res, &M[bitbuf], &res)) != MP_OKAY) {
2884 if ((err = redux (&res, P, mp)) != MP_OKAY) {
2888 /* empty window and reset */
2895 /* if bits remain then square/multiply */
2896 if (mode == 2 && bitcpy > 0) {
2897 /* square then multiply if the bit is set */
2898 for (x = 0; x < bitcpy; x++) {
2899 if ((err = tma_mp_sqr (&res, &res)) != MP_OKAY) {
2902 if ((err = redux (&res, P, mp)) != MP_OKAY) {
2906 /* get next bit of the window */
2908 if ((bitbuf & (1 << winsize)) != 0) {
2910 if ((err = tma_mp_mul (&res, &M[1], &res)) != MP_OKAY) {
2913 if ((err = redux (&res, P, mp)) != MP_OKAY) {
2921 /* fixup result if Montgomery reduction is used
2922 * recall that any value in a Montgomery system is
2923 * actually multiplied by R mod n. So we have
2924 * to reduce one more time to cancel out the factor
2927 if ((err = redux(&res, P, mp)) != MP_OKAY) {
2932 /* swap res with Y */
2933 tma_mp_exch (&res, Y);
2935 LBL_RES:tma_mp_clear (&res);
2937 tma_mp_clear(&M[1]);
2938 for (x = 1<<(winsize-1); x < (1 << winsize); x++) {
2939 tma_mp_clear (&M[x]);
2950 /* End: bn_tma_mp_exptmod_fast.c */
2952 /* Start: bn_tma_mp_exteuclid.c */
2954 #ifdef BN_MP_EXTEUCLID_C
2955 /* LibTomMath, multiple-precision integer library -- Tom St Denis
2957 * LibTomMath is a library that provides multiple-precision
2958 * integer arithmetic as well as number theoretic functionality.
2960 * The library was designed directly after the MPI library by
2961 * Michael Fromberger but has been written from scratch with
2962 * additional optimizations in place.
2964 * The library is free for all purposes without any express
2965 * guarantee it works.
2967 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
2970 /* Extended euclidean algorithm of (a, b) produces
2973 int tma_mp_exteuclid(tma_mp_int *a, tma_mp_int *b, tma_mp_int *U1, tma_mp_int *U2, tma_mp_int *U3)
2975 tma_mp_int u1,u2,u3,v1,v2,v3,t1,t2,t3,q,tmp;
2978 if ((err = tma_mp_init_multi(&u1, &u2, &u3, &v1, &v2, &v3, &t1, &t2, &t3, &q, &tmp, NULL)) != MP_OKAY) {
2982 /* initialize, (u1,u2,u3) = (1,0,a) */
2984 if ((err = tma_mp_copy(a, &u3)) != MP_OKAY) { goto _ERR; }
2986 /* initialize, (v1,v2,v3) = (0,1,b) */
2988 if ((err = tma_mp_copy(b, &v3)) != MP_OKAY) { goto _ERR; }
2990 /* loop while v3 != 0 */
2991 while (tma_mp_iszero(&v3) == MP_NO) {
2993 if ((err = tma_mp_div(&u3, &v3, &q, NULL)) != MP_OKAY) { goto _ERR; }
2995 /* (t1,t2,t3) = (u1,u2,u3) - (v1,v2,v3)q */
2996 if ((err = tma_mp_mul(&v1, &q, &tmp)) != MP_OKAY) { goto _ERR; }
2997 if ((err = tma_mp_sub(&u1, &tmp, &t1)) != MP_OKAY) { goto _ERR; }
2998 if ((err = tma_mp_mul(&v2, &q, &tmp)) != MP_OKAY) { goto _ERR; }
2999 if ((err = tma_mp_sub(&u2, &tmp, &t2)) != MP_OKAY) { goto _ERR; }
3000 if ((err = tma_mp_mul(&v3, &q, &tmp)) != MP_OKAY) { goto _ERR; }
3001 if ((err = tma_mp_sub(&u3, &tmp, &t3)) != MP_OKAY) { goto _ERR; }
3003 /* (u1,u2,u3) = (v1,v2,v3) */
3004 if ((err = tma_mp_copy(&v1, &u1)) != MP_OKAY) { goto _ERR; }
3005 if ((err = tma_mp_copy(&v2, &u2)) != MP_OKAY) { goto _ERR; }
3006 if ((err = tma_mp_copy(&v3, &u3)) != MP_OKAY) { goto _ERR; }
3008 /* (v1,v2,v3) = (t1,t2,t3) */
3009 if ((err = tma_mp_copy(&t1, &v1)) != MP_OKAY) { goto _ERR; }
3010 if ((err = tma_mp_copy(&t2, &v2)) != MP_OKAY) { goto _ERR; }
3011 if ((err = tma_mp_copy(&t3, &v3)) != MP_OKAY) { goto _ERR; }
3014 /* make sure U3 >= 0 */
3015 if (u3.sign == MP_NEG) {
3016 tma_mp_neg(&u1, &u1);
3017 tma_mp_neg(&u2, &u2);
3018 tma_mp_neg(&u3, &u3);
3021 /* copy result out */
3022 if (U1 != NULL) { tma_mp_exch(U1, &u1); }
3023 if (U2 != NULL) { tma_mp_exch(U2, &u2); }
3024 if (U3 != NULL) { tma_mp_exch(U3, &u3); }
3027 _ERR: tma_mp_clear_multi(&u1, &u2, &u3, &v1, &v2, &v3, &t1, &t2, &t3, &q, &tmp, NULL);
3036 /* End: bn_tma_mp_exteuclid.c */
3038 /* Start: bn_tma_mp_fread.c */
3040 #ifdef BN_MP_FREAD_C
3041 /* LibTomMath, multiple-precision integer library -- Tom St Denis
3043 * LibTomMath is a library that provides multiple-precision
3044 * integer arithmetic as well as number theoretic functionality.
3046 * The library was designed directly after the MPI library by
3047 * Michael Fromberger but has been written from scratch with
3048 * additional optimizations in place.
3050 * The library is free for all purposes without any express
3051 * guarantee it works.
3053 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
3056 /* read a bigint from a file stream in ASCII */
3057 int tma_mp_fread(tma_mp_int *a, int radix, FILE *stream)
3059 int err, ch, neg, y;
3064 /* if first digit is - then set negative */
3074 /* find y in the radix map */
3075 for (y = 0; y < radix; y++) {
3076 if (tma_mp_s_rmap[y] == ch) {
3084 /* shift up and add */
3085 if ((err = tma_mp_mul_d(a, radix, a)) != MP_OKAY) {
3088 if ((err = tma_mp_add_d(a, y, a)) != MP_OKAY) {
3094 if (tma_mp_cmp_d(a, 0) != MP_EQ) {
3107 /* End: bn_tma_mp_fread.c */
3109 /* Start: bn_tma_mp_fwrite.c */
3111 #ifdef BN_MP_FWRITE_C
3112 /* LibTomMath, multiple-precision integer library -- Tom St Denis
3114 * LibTomMath is a library that provides multiple-precision
3115 * integer arithmetic as well as number theoretic functionality.
3117 * The library was designed directly after the MPI library by
3118 * Michael Fromberger but has been written from scratch with
3119 * additional optimizations in place.
3121 * The library is free for all purposes without any express
3122 * guarantee it works.
3124 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
3127 int tma_mp_fwrite(tma_mp_int *a, int radix, FILE *stream)
3132 if ((err = tma_mp_radix_size(a, radix, &len)) != MP_OKAY) {
3136 buf = OPT_CAST(char) XMALLOC (len);
3141 if ((err = tma_mp_toradix(a, buf, radix)) != MP_OKAY) {
3146 for (x = 0; x < len; x++) {
3147 if (fputc(buf[x], stream) == EOF) {
3163 /* End: bn_tma_mp_fwrite.c */
3165 /* Start: bn_tma_mp_gcd.c */
3168 /* LibTomMath, multiple-precision integer library -- Tom St Denis
3170 * LibTomMath is a library that provides multiple-precision
3171 * integer arithmetic as well as number theoretic functionality.
3173 * The library was designed directly after the MPI library by
3174 * Michael Fromberger but has been written from scratch with
3175 * additional optimizations in place.
3177 * The library is free for all purposes without any express
3178 * guarantee it works.
3180 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
3183 /* Greatest Common Divisor using the binary method */
3184 int tma_mp_gcd (tma_mp_int * a, tma_mp_int * b, tma_mp_int * c)
3187 int k, u_lsb, v_lsb, res;
3189 /* either zero than gcd is the largest */
3190 if (tma_mp_iszero (a) == MP_YES) {
3191 return tma_mp_abs (b, c);
3193 if (tma_mp_iszero (b) == MP_YES) {
3194 return tma_mp_abs (a, c);
3197 /* get copies of a and b we can modify */
3198 if ((res = tma_mp_init_copy (&u, a)) != MP_OKAY) {
3202 if ((res = tma_mp_init_copy (&v, b)) != MP_OKAY) {
3206 /* must be positive for the remainder of the algorithm */
3207 u.sign = v.sign = MP_ZPOS;
3209 /* B1. Find the common power of two for u and v */
3210 u_lsb = tma_mp_cnt_lsb(&u);
3211 v_lsb = tma_mp_cnt_lsb(&v);
3212 k = MIN(u_lsb, v_lsb);
3215 /* divide the power of two out */
3216 if ((res = tma_mp_div_2d(&u, k, &u, NULL)) != MP_OKAY) {
3220 if ((res = tma_mp_div_2d(&v, k, &v, NULL)) != MP_OKAY) {
3225 /* divide any remaining factors of two out */
3227 if ((res = tma_mp_div_2d(&u, u_lsb - k, &u, NULL)) != MP_OKAY) {
3233 if ((res = tma_mp_div_2d(&v, v_lsb - k, &v, NULL)) != MP_OKAY) {
3238 while (tma_mp_iszero(&v) == 0) {
3239 /* make sure v is the largest */
3240 if (tma_mp_cmp_mag(&u, &v) == MP_GT) {
3241 /* swap u and v to make sure v is >= u */
3242 tma_mp_exch(&u, &v);
3245 /* subtract smallest from largest */
3246 if ((res = s_tma_mp_sub(&v, &u, &v)) != MP_OKAY) {
3250 /* Divide out all factors of two */
3251 if ((res = tma_mp_div_2d(&v, tma_mp_cnt_lsb(&v), &v, NULL)) != MP_OKAY) {
3256 /* multiply by 2**k which we divided out at the beginning */
3257 if ((res = tma_mp_mul_2d (&u, k, c)) != MP_OKAY) {
3262 LBL_V:tma_mp_clear (&u);
3263 LBL_U:tma_mp_clear (&v);
3272 /* End: bn_tma_mp_gcd.c */
3274 /* Start: bn_tma_mp_get_int.c */
3276 #ifdef BN_MP_GET_INT_C
3277 /* LibTomMath, multiple-precision integer library -- Tom St Denis
3279 * LibTomMath is a library that provides multiple-precision
3280 * integer arithmetic as well as number theoretic functionality.
3282 * The library was designed directly after the MPI library by
3283 * Michael Fromberger but has been written from scratch with
3284 * additional optimizations in place.
3286 * The library is free for all purposes without any express
3287 * guarantee it works.
3289 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
3292 /* get the lower 32-bits of an tma_mp_int */
3293 unsigned long tma_mp_get_int(tma_mp_int * a)
3302 /* get number of digits of the lsb we have to read */
3303 i = MIN(a->used,(int)((sizeof(unsigned long)*CHAR_BIT+DIGIT_BIT-1)/DIGIT_BIT))-1;
3305 /* get most significant digit of result */
3309 res = (res << DIGIT_BIT) | DIGIT(a,i);
3312 /* force result to 32-bits always so it is consistent on non 32-bit platforms */
3313 return res & 0xFFFFFFFFUL;
3321 /* End: bn_tma_mp_get_int.c */
3323 /* Start: bn_tma_mp_grow.c */
3326 /* LibTomMath, multiple-precision integer library -- Tom St Denis
3328 * LibTomMath is a library that provides multiple-precision
3329 * integer arithmetic as well as number theoretic functionality.
3331 * The library was designed directly after the MPI library by
3332 * Michael Fromberger but has been written from scratch with
3333 * additional optimizations in place.
3335 * The library is free for all purposes without any express
3336 * guarantee it works.
3338 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
3341 /* grow as required */
3342 int tma_mp_grow (tma_mp_int * a, int size)
3347 /* if the alloc size is smaller alloc more ram */
3348 if (a->alloc < size) {
3349 /* ensure there are always at least MP_PREC digits extra on top */
3350 size += (MP_PREC * 2) - (size % MP_PREC);
3352 /* reallocate the array a->dp
3354 * We store the return in a temporary variable
3355 * in case the operation failed we don't want
3356 * to overwrite the dp member of a.
3358 tmp = OPT_CAST(tma_mp_digit) XREALLOC (a->dp, sizeof (tma_mp_digit) * size);
3360 /* reallocation failed but "a" is still valid [can be freed] */
3364 /* reallocation succeeded so set a->dp */
3367 /* zero excess digits */
3370 for (; i < a->alloc; i++) {
3382 /* End: bn_tma_mp_grow.c */
3384 /* Start: bn_tma_mp_init.c */
3387 /* LibTomMath, multiple-precision integer library -- Tom St Denis
3389 * LibTomMath is a library that provides multiple-precision
3390 * integer arithmetic as well as number theoretic functionality.
3392 * The library was designed directly after the MPI library by
3393 * Michael Fromberger but has been written from scratch with
3394 * additional optimizations in place.
3396 * The library is free for all purposes without any express
3397 * guarantee it works.
3399 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
3402 /* init a new tma_mp_int */
3403 int tma_mp_init (tma_mp_int * a)
3407 /* allocate memory required and clear it */
3408 a->dp = OPT_CAST(tma_mp_digit) XMALLOC (sizeof (tma_mp_digit) * MP_PREC);
3409 if (a->dp == NULL) {
3413 /* set the digits to zero */
3414 for (i = 0; i < MP_PREC; i++) {
3418 /* set the used to zero, allocated digits to the default precision
3419 * and sign to positive */
3432 /* End: bn_tma_mp_init.c */
3434 /* Start: bn_tma_mp_init_copy.c */
3436 #ifdef BN_MP_INIT_COPY_C
3437 /* LibTomMath, multiple-precision integer library -- Tom St Denis
3439 * LibTomMath is a library that provides multiple-precision
3440 * integer arithmetic as well as number theoretic functionality.
3442 * The library was designed directly after the MPI library by
3443 * Michael Fromberger but has been written from scratch with
3444 * additional optimizations in place.
3446 * The library is free for all purposes without any express
3447 * guarantee it works.
3449 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
3452 /* creates "a" then copies b into it */
3453 int tma_mp_init_copy (tma_mp_int * a, tma_mp_int * b)
3457 if ((res = tma_mp_init (a)) != MP_OKAY) {
3460 return tma_mp_copy (b, a);
3468 /* End: bn_tma_mp_init_copy.c */
3470 /* Start: bn_tma_mp_init_multi.c */
3472 #ifdef BN_MP_INIT_MULTI_C
3473 /* LibTomMath, multiple-precision integer library -- Tom St Denis
3475 * LibTomMath is a library that provides multiple-precision
3476 * integer arithmetic as well as number theoretic functionality.
3478 * The library was designed directly after the MPI library by
3479 * Michael Fromberger but has been written from scratch with
3480 * additional optimizations in place.
3482 * The library is free for all purposes without any express
3483 * guarantee it works.
3485 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
3489 int tma_mp_init_multi(tma_mp_int *mp, ...)
3491 tma_mp_err res = MP_OKAY; /* Assume ok until proven otherwise */
3492 int n = 0; /* Number of ok inits */
3493 tma_mp_int* cur_arg = mp;
3496 va_start(args, mp); /* init args to next argument from caller */
3497 while (cur_arg != NULL) {
3498 if (tma_mp_init(cur_arg) != MP_OKAY) {
3499 /* Oops - error! Back-track and tma_mp_clear what we already
3500 succeeded in init-ing, then return error.
3504 /* end the current list */
3507 /* now start cleaning up */
3509 va_start(clean_args, mp);
3511 tma_mp_clear(cur_arg);
3512 cur_arg = va_arg(clean_args, tma_mp_int*);
3519 cur_arg = va_arg(args, tma_mp_int*);
3522 return res; /* Assumed ok, if error flagged above. */
3531 /* End: bn_tma_mp_init_multi.c */
3533 /* Start: bn_tma_mp_init_set.c */
3535 #ifdef BN_MP_INIT_SET_C
3536 /* LibTomMath, multiple-precision integer library -- Tom St Denis
3538 * LibTomMath is a library that provides multiple-precision
3539 * integer arithmetic as well as number theoretic functionality.
3541 * The library was designed directly after the MPI library by
3542 * Michael Fromberger but has been written from scratch with
3543 * additional optimizations in place.
3545 * The library is free for all purposes without any express
3546 * guarantee it works.
3548 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
3551 /* initialize and set a digit */
3552 int tma_mp_init_set (tma_mp_int * a, tma_mp_digit b)
3555 if ((err = tma_mp_init(a)) != MP_OKAY) {
3567 /* End: bn_tma_mp_init_set.c */
3569 /* Start: bn_tma_mp_init_set_int.c */
3571 #ifdef BN_MP_INIT_SET_INT_C
3572 /* LibTomMath, multiple-precision integer library -- Tom St Denis
3574 * LibTomMath is a library that provides multiple-precision
3575 * integer arithmetic as well as number theoretic functionality.
3577 * The library was designed directly after the MPI library by
3578 * Michael Fromberger but has been written from scratch with
3579 * additional optimizations in place.
3581 * The library is free for all purposes without any express
3582 * guarantee it works.
3584 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
3587 /* initialize and set a digit */
3588 int tma_mp_init_set_int (tma_mp_int * a, unsigned long b)
3591 if ((err = tma_mp_init(a)) != MP_OKAY) {
3594 return tma_mp_set_int(a, b);
3602 /* End: bn_tma_mp_init_set_int.c */
3604 /* Start: bn_tma_mp_init_size.c */
3606 #ifdef BN_MP_INIT_SIZE_C
3607 /* LibTomMath, multiple-precision integer library -- Tom St Denis
3609 * LibTomMath is a library that provides multiple-precision
3610 * integer arithmetic as well as number theoretic functionality.
3612 * The library was designed directly after the MPI library by
3613 * Michael Fromberger but has been written from scratch with
3614 * additional optimizations in place.
3616 * The library is free for all purposes without any express
3617 * guarantee it works.
3619 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
3622 /* init an tma_mp_init for a given size */
3623 int tma_mp_init_size (tma_mp_int * a, int size)
3627 /* pad size so there are always extra digits */
3628 size += (MP_PREC * 2) - (size % MP_PREC);
3631 a->dp = OPT_CAST(tma_mp_digit) XMALLOC (sizeof (tma_mp_digit) * size);
3632 if (a->dp == NULL) {
3636 /* set the members */
3641 /* zero the digits */
3642 for (x = 0; x < size; x++) {
3654 /* End: bn_tma_mp_init_size.c */
3656 /* Start: bn_tma_mp_invmod.c */
3658 #ifdef BN_MP_INVMOD_C
3659 /* LibTomMath, multiple-precision integer library -- Tom St Denis
3661 * LibTomMath is a library that provides multiple-precision
3662 * integer arithmetic as well as number theoretic functionality.
3664 * The library was designed directly after the MPI library by
3665 * Michael Fromberger but has been written from scratch with
3666 * additional optimizations in place.
3668 * The library is free for all purposes without any express
3669 * guarantee it works.
3671 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
3674 /* hac 14.61, pp608 */
3675 int tma_mp_invmod (tma_mp_int * a, tma_mp_int * b, tma_mp_int * c)
3677 /* b cannot be negative */
3678 if (b->sign == MP_NEG || tma_mp_iszero(b) == 1) {
3682 #ifdef BN_FAST_MP_INVMOD_C
3683 /* if the modulus is odd we can use a faster routine instead */
3684 if (tma_mp_isodd (b) == 1) {
3685 return fast_tma_mp_invmod (a, b, c);
3689 #ifdef BN_MP_INVMOD_SLOW_C
3690 return tma_mp_invmod_slow(a, b, c);
3701 /* End: bn_tma_mp_invmod.c */
3703 /* Start: bn_tma_mp_invmod_slow.c */
3705 #ifdef BN_MP_INVMOD_SLOW_C
3706 /* LibTomMath, multiple-precision integer library -- Tom St Denis
3708 * LibTomMath is a library that provides multiple-precision
3709 * integer arithmetic as well as number theoretic functionality.
3711 * The library was designed directly after the MPI library by
3712 * Michael Fromberger but has been written from scratch with
3713 * additional optimizations in place.
3715 * The library is free for all purposes without any express
3716 * guarantee it works.
3718 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
3721 /* hac 14.61, pp608 */
3722 int tma_mp_invmod_slow (tma_mp_int * a, tma_mp_int * b, tma_mp_int * c)
3724 tma_mp_int x, y, u, v, A, B, C, D;
3727 /* b cannot be negative */
3728 if (b->sign == MP_NEG || tma_mp_iszero(b) == 1) {
3733 if ((res = tma_mp_init_multi(&x, &y, &u, &v,
3734 &A, &B, &C, &D, NULL)) != MP_OKAY) {
3739 if ((res = tma_mp_mod(a, b, &x)) != MP_OKAY) {
3742 if ((res = tma_mp_copy (b, &y)) != MP_OKAY) {
3746 /* 2. [modified] if x,y are both even then return an error! */
3747 if (tma_mp_iseven (&x) == 1 && tma_mp_iseven (&y) == 1) {
3752 /* 3. u=x, v=y, A=1, B=0, C=0,D=1 */
3753 if ((res = tma_mp_copy (&x, &u)) != MP_OKAY) {
3756 if ((res = tma_mp_copy (&y, &v)) != MP_OKAY) {
3763 /* 4. while u is even do */
3764 while (tma_mp_iseven (&u) == 1) {
3766 if ((res = tma_mp_div_2 (&u, &u)) != MP_OKAY) {
3769 /* 4.2 if A or B is odd then */
3770 if (tma_mp_isodd (&A) == 1 || tma_mp_isodd (&B) == 1) {
3771 /* A = (A+y)/2, B = (B-x)/2 */
3772 if ((res = tma_mp_add (&A, &y, &A)) != MP_OKAY) {
3775 if ((res = tma_mp_sub (&B, &x, &B)) != MP_OKAY) {
3779 /* A = A/2, B = B/2 */
3780 if ((res = tma_mp_div_2 (&A, &A)) != MP_OKAY) {
3783 if ((res = tma_mp_div_2 (&B, &B)) != MP_OKAY) {
3788 /* 5. while v is even do */
3789 while (tma_mp_iseven (&v) == 1) {
3791 if ((res = tma_mp_div_2 (&v, &v)) != MP_OKAY) {
3794 /* 5.2 if C or D is odd then */
3795 if (tma_mp_isodd (&C) == 1 || tma_mp_isodd (&D) == 1) {
3796 /* C = (C+y)/2, D = (D-x)/2 */
3797 if ((res = tma_mp_add (&C, &y, &C)) != MP_OKAY) {
3800 if ((res = tma_mp_sub (&D, &x, &D)) != MP_OKAY) {
3804 /* C = C/2, D = D/2 */
3805 if ((res = tma_mp_div_2 (&C, &C)) != MP_OKAY) {
3808 if ((res = tma_mp_div_2 (&D, &D)) != MP_OKAY) {
3813 /* 6. if u >= v then */
3814 if (tma_mp_cmp (&u, &v) != MP_LT) {
3815 /* u = u - v, A = A - C, B = B - D */
3816 if ((res = tma_mp_sub (&u, &v, &u)) != MP_OKAY) {
3820 if ((res = tma_mp_sub (&A, &C, &A)) != MP_OKAY) {
3824 if ((res = tma_mp_sub (&B, &D, &B)) != MP_OKAY) {
3828 /* v - v - u, C = C - A, D = D - B */
3829 if ((res = tma_mp_sub (&v, &u, &v)) != MP_OKAY) {
3833 if ((res = tma_mp_sub (&C, &A, &C)) != MP_OKAY) {
3837 if ((res = tma_mp_sub (&D, &B, &D)) != MP_OKAY) {
3842 /* if not zero goto step 4 */
3843 if (tma_mp_iszero (&u) == 0)
3846 /* now a = C, b = D, gcd == g*v */
3848 /* if v != 1 then there is no inverse */
3849 if (tma_mp_cmp_d (&v, 1) != MP_EQ) {
3854 /* if its too low */
3855 while (tma_mp_cmp_d(&C, 0) == MP_LT) {
3856 if ((res = tma_mp_add(&C, b, &C)) != MP_OKAY) {
3862 while (tma_mp_cmp_mag(&C, b) != MP_LT) {
3863 if ((res = tma_mp_sub(&C, b, &C)) != MP_OKAY) {
3868 /* C is now the inverse */
3869 tma_mp_exch (&C, c);
3871 LBL_ERR:tma_mp_clear_multi (&x, &y, &u, &v, &A, &B, &C, &D, NULL);
3880 /* End: bn_tma_mp_invmod_slow.c */
3882 /* Start: bn_tma_mp_is_square.c */
3884 #ifdef BN_MP_IS_SQUARE_C
3885 /* LibTomMath, multiple-precision integer library -- Tom St Denis
3887 * LibTomMath is a library that provides multiple-precision
3888 * integer arithmetic as well as number theoretic functionality.
3890 * The library was designed directly after the MPI library by
3891 * Michael Fromberger but has been written from scratch with
3892 * additional optimizations in place.
3894 * The library is free for all purposes without any express
3895 * guarantee it works.
3897 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
3900 /* Check if remainders are possible squares - fast exclude non-squares */
3901 static const char rem_128[128] = {
3902 0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
3903 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
3904 1, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
3905 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
3906 0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
3907 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
3908 1, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
3909 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1
3912 static const char rem_105[105] = {
3913 0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1,
3914 0, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1,
3915 0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1,
3916 1, 0, 1, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1,
3917 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1,
3918 1, 1, 1, 1, 0, 1, 0, 1, 1, 0, 0, 1, 1, 1, 1,
3919 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1
3922 /* Store non-zero to ret if arg is square, and zero if not */
3923 int tma_mp_is_square(tma_mp_int *arg,int *ret)
3930 /* Default to Non-square :) */
3933 if (arg->sign == MP_NEG) {
3937 /* digits used? (TSD) */
3938 if (arg->used == 0) {
3942 /* First check mod 128 (suppose that DIGIT_BIT is at least 7) */
3943 if (rem_128[127 & DIGIT(arg,0)] == 1) {
3947 /* Next check mod 105 (3*5*7) */
3948 if ((res = tma_mp_mod_d(arg,105,&c)) != MP_OKAY) {
3951 if (rem_105[c] == 1) {
3956 if ((res = tma_mp_init_set_int(&t,11L*13L*17L*19L*23L*29L*31L)) != MP_OKAY) {
3959 if ((res = tma_mp_mod(arg,&t,&t)) != MP_OKAY) {
3962 r = tma_mp_get_int(&t);
3963 /* Check for other prime modules, note it's not an ERROR but we must
3964 * free "t" so the easiest way is to goto ERR. We know that res
3965 * is already equal to MP_OKAY from the tma_mp_mod call
3967 if ( (1L<<(r%11)) & 0x5C4L ) goto ERR;
3968 if ( (1L<<(r%13)) & 0x9E4L ) goto ERR;
3969 if ( (1L<<(r%17)) & 0x5CE8L ) goto ERR;
3970 if ( (1L<<(r%19)) & 0x4F50CL ) goto ERR;
3971 if ( (1L<<(r%23)) & 0x7ACCA0L ) goto ERR;
3972 if ( (1L<<(r%29)) & 0xC2EDD0CL ) goto ERR;
3973 if ( (1L<<(r%31)) & 0x6DE2B848L ) goto ERR;
3975 /* Final check - is sqr(sqrt(arg)) == arg ? */
3976 if ((res = tma_mp_sqrt(arg,&t)) != MP_OKAY) {
3979 if ((res = tma_mp_sqr(&t,&t)) != MP_OKAY) {
3983 *ret = (tma_mp_cmp_mag(&t,arg) == MP_EQ) ? MP_YES : MP_NO;
3984 ERR:tma_mp_clear(&t);
3993 /* End: bn_tma_mp_is_square.c */
3995 /* Start: bn_tma_mp_jacobi.c */
3997 #ifdef BN_MP_JACOBI_C
3998 /* LibTomMath, multiple-precision integer library -- Tom St Denis
4000 * LibTomMath is a library that provides multiple-precision
4001 * integer arithmetic as well as number theoretic functionality.
4003 * The library was designed directly after the MPI library by
4004 * Michael Fromberger but has been written from scratch with
4005 * additional optimizations in place.
4007 * The library is free for all purposes without any express
4008 * guarantee it works.
4010 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
4013 /* computes the jacobi c = (a | n) (or Legendre if n is prime)
4014 * HAC pp. 73 Algorithm 2.149
4016 int tma_mp_jacobi (tma_mp_int * a, tma_mp_int * p, int *c)
4020 tma_mp_digit residue;
4022 /* if p <= 0 return MP_VAL */
4023 if (tma_mp_cmp_d(p, 0) != MP_GT) {
4027 /* step 1. if a == 0, return 0 */
4028 if (tma_mp_iszero (a) == 1) {
4033 /* step 2. if a == 1, return 1 */
4034 if (tma_mp_cmp_d (a, 1) == MP_EQ) {
4042 /* step 3. write a = a1 * 2**k */
4043 if ((res = tma_mp_init_copy (&a1, a)) != MP_OKAY) {
4047 if ((res = tma_mp_init (&p1)) != MP_OKAY) {
4051 /* divide out larger power of two */
4052 k = tma_mp_cnt_lsb(&a1);
4053 if ((res = tma_mp_div_2d(&a1, k, &a1, NULL)) != MP_OKAY) {
4057 /* step 4. if e is even set s=1 */
4061 /* else set s=1 if p = 1/7 (mod 8) or s=-1 if p = 3/5 (mod 8) */
4062 residue = p->dp[0] & 7;
4064 if (residue == 1 || residue == 7) {
4066 } else if (residue == 3 || residue == 5) {
4071 /* step 5. if p == 3 (mod 4) *and* a1 == 3 (mod 4) then s = -s */
4072 if ( ((p->dp[0] & 3) == 3) && ((a1.dp[0] & 3) == 3)) {
4076 /* if a1 == 1 we're done */
4077 if (tma_mp_cmp_d (&a1, 1) == MP_EQ) {
4081 if ((res = tma_mp_mod (p, &a1, &p1)) != MP_OKAY) {
4084 if ((res = tma_mp_jacobi (&p1, &a1, &r)) != MP_OKAY) {
4092 LBL_P1:tma_mp_clear (&p1);
4093 LBL_A1:tma_mp_clear (&a1);
4102 /* End: bn_tma_mp_jacobi.c */
4104 /* Start: bn_tma_mp_karatsuba_mul.c */
4106 #ifdef BN_MP_KARATSUBA_MUL_C
4107 /* LibTomMath, multiple-precision integer library -- Tom St Denis
4109 * LibTomMath is a library that provides multiple-precision
4110 * integer arithmetic as well as number theoretic functionality.
4112 * The library was designed directly after the MPI library by
4113 * Michael Fromberger but has been written from scratch with
4114 * additional optimizations in place.
4116 * The library is free for all purposes without any express
4117 * guarantee it works.
4119 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
4122 /* c = |a| * |b| using Karatsuba Multiplication using
4123 * three half size multiplications
4125 * Let B represent the radix [e.g. 2**DIGIT_BIT] and
4126 * let n represent half of the number of digits in
4129 * a = a1 * B**n + a0
4130 * b = b1 * B**n + b0
4133 a1b1 * B**2n + ((a1 + a0)(b1 + b0) - (a0b0 + a1b1)) * B + a0b0
4135 * Note that a1b1 and a0b0 are used twice and only need to be
4136 * computed once. So in total three half size (half # of
4137 * digit) multiplications are performed, a0b0, a1b1 and
4140 * Note that a multiplication of half the digits requires
4141 * 1/4th the number of single precision multiplications so in
4142 * total after one call 25% of the single precision multiplications
4143 * are saved. Note also that the call to tma_mp_mul can end up back
4144 * in this function if the a0, a1, b0, or b1 are above the threshold.
4145 * This is known as divide-and-conquer and leads to the famous
4146 * O(N**lg(3)) or O(N**1.584) work which is asymptopically lower than
4147 * the standard O(N**2) that the baseline/comba methods use.
4148 * Generally though the overhead of this method doesn't pay off
4149 * until a certain size (N ~ 80) is reached.
4151 int tma_mp_karatsuba_mul (tma_mp_int * a, tma_mp_int * b, tma_mp_int * c)
4153 tma_mp_int x0, x1, y0, y1, t1, x0y0, x1y1;
4156 /* default the return code to an error */
4159 /* min # of digits */
4160 B = MIN (a->used, b->used);
4162 /* now divide in two */
4165 /* init copy all the temps */
4166 if (tma_mp_init_size (&x0, B) != MP_OKAY)
4168 if (tma_mp_init_size (&x1, a->used - B) != MP_OKAY)
4170 if (tma_mp_init_size (&y0, B) != MP_OKAY)
4172 if (tma_mp_init_size (&y1, b->used - B) != MP_OKAY)
4176 if (tma_mp_init_size (&t1, B * 2) != MP_OKAY)
4178 if (tma_mp_init_size (&x0y0, B * 2) != MP_OKAY)
4180 if (tma_mp_init_size (&x1y1, B * 2) != MP_OKAY)
4183 /* now shift the digits */
4184 x0.used = y0.used = B;
4185 x1.used = a->used - B;
4186 y1.used = b->used - B;
4190 register tma_mp_digit *tmpa, *tmpb, *tmpx, *tmpy;
4192 /* we copy the digits directly instead of using higher level functions
4193 * since we also need to shift the digits
4200 for (x = 0; x < B; x++) {
4206 for (x = B; x < a->used; x++) {
4211 for (x = B; x < b->used; x++) {
4216 /* only need to clamp the lower words since by definition the
4217 * upper words x1/y1 must have a known number of digits
4222 /* now calc the products x0y0 and x1y1 */
4223 /* after this x0 is no longer required, free temp [x0==t2]! */
4224 if (tma_mp_mul (&x0, &y0, &x0y0) != MP_OKAY)
4225 goto X1Y1; /* x0y0 = x0*y0 */
4226 if (tma_mp_mul (&x1, &y1, &x1y1) != MP_OKAY)
4227 goto X1Y1; /* x1y1 = x1*y1 */
4229 /* now calc x1+x0 and y1+y0 */
4230 if (s_tma_mp_add (&x1, &x0, &t1) != MP_OKAY)
4231 goto X1Y1; /* t1 = x1 - x0 */
4232 if (s_tma_mp_add (&y1, &y0, &x0) != MP_OKAY)
4233 goto X1Y1; /* t2 = y1 - y0 */
4234 if (tma_mp_mul (&t1, &x0, &t1) != MP_OKAY)
4235 goto X1Y1; /* t1 = (x1 + x0) * (y1 + y0) */
4238 if (tma_mp_add (&x0y0, &x1y1, &x0) != MP_OKAY)
4239 goto X1Y1; /* t2 = x0y0 + x1y1 */
4240 if (s_tma_mp_sub (&t1, &x0, &t1) != MP_OKAY)
4241 goto X1Y1; /* t1 = (x1+x0)*(y1+y0) - (x1y1 + x0y0) */
4244 if (tma_mp_lshd (&t1, B) != MP_OKAY)
4245 goto X1Y1; /* t1 = (x0y0 + x1y1 - (x1-x0)*(y1-y0))<<B */
4246 if (tma_mp_lshd (&x1y1, B * 2) != MP_OKAY)
4247 goto X1Y1; /* x1y1 = x1y1 << 2*B */
4249 if (tma_mp_add (&x0y0, &t1, &t1) != MP_OKAY)
4250 goto X1Y1; /* t1 = x0y0 + t1 */
4251 if (tma_mp_add (&t1, &x1y1, c) != MP_OKAY)
4252 goto X1Y1; /* t1 = x0y0 + t1 + x1y1 */
4254 /* Algorithm succeeded set the return code to MP_OKAY */
4257 X1Y1:tma_mp_clear (&x1y1);
4258 X0Y0:tma_mp_clear (&x0y0);
4259 T1:tma_mp_clear (&t1);
4260 Y1:tma_mp_clear (&y1);
4261 Y0:tma_mp_clear (&y0);
4262 X1:tma_mp_clear (&x1);
4263 X0:tma_mp_clear (&x0);
4273 /* End: bn_tma_mp_karatsuba_mul.c */
4275 /* Start: bn_tma_mp_karatsuba_sqr.c */
4277 #ifdef BN_MP_KARATSUBA_SQR_C
4278 /* LibTomMath, multiple-precision integer library -- Tom St Denis
4280 * LibTomMath is a library that provides multiple-precision
4281 * integer arithmetic as well as number theoretic functionality.
4283 * The library was designed directly after the MPI library by
4284 * Michael Fromberger but has been written from scratch with
4285 * additional optimizations in place.
4287 * The library is free for all purposes without any express
4288 * guarantee it works.
4290 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
4293 /* Karatsuba squaring, computes b = a*a using three
4294 * half size squarings
4296 * See comments of karatsuba_mul for details. It
4297 * is essentially the same algorithm but merely
4298 * tuned to perform recursive squarings.
4300 int tma_mp_karatsuba_sqr (tma_mp_int * a, tma_mp_int * b)
4302 tma_mp_int x0, x1, t1, t2, x0x0, x1x1;
4307 /* min # of digits */
4310 /* now divide in two */
4313 /* init copy all the temps */
4314 if (tma_mp_init_size (&x0, B) != MP_OKAY)
4316 if (tma_mp_init_size (&x1, a->used - B) != MP_OKAY)
4320 if (tma_mp_init_size (&t1, a->used * 2) != MP_OKAY)
4322 if (tma_mp_init_size (&t2, a->used * 2) != MP_OKAY)
4324 if (tma_mp_init_size (&x0x0, B * 2) != MP_OKAY)
4326 if (tma_mp_init_size (&x1x1, (a->used - B) * 2) != MP_OKAY)
4331 register tma_mp_digit *dst, *src;
4335 /* now shift the digits */
4337 for (x = 0; x < B; x++) {
4342 for (x = B; x < a->used; x++) {
4348 x1.used = a->used - B;
4352 /* now calc the products x0*x0 and x1*x1 */
4353 if (tma_mp_sqr (&x0, &x0x0) != MP_OKAY)
4354 goto X1X1; /* x0x0 = x0*x0 */
4355 if (tma_mp_sqr (&x1, &x1x1) != MP_OKAY)
4356 goto X1X1; /* x1x1 = x1*x1 */
4358 /* now calc (x1+x0)**2 */
4359 if (s_tma_mp_add (&x1, &x0, &t1) != MP_OKAY)
4360 goto X1X1; /* t1 = x1 - x0 */
4361 if (tma_mp_sqr (&t1, &t1) != MP_OKAY)
4362 goto X1X1; /* t1 = (x1 - x0) * (x1 - x0) */
4365 if (s_tma_mp_add (&x0x0, &x1x1, &t2) != MP_OKAY)
4366 goto X1X1; /* t2 = x0x0 + x1x1 */
4367 if (s_tma_mp_sub (&t1, &t2, &t1) != MP_OKAY)
4368 goto X1X1; /* t1 = (x1+x0)**2 - (x0x0 + x1x1) */
4371 if (tma_mp_lshd (&t1, B) != MP_OKAY)
4372 goto X1X1; /* t1 = (x0x0 + x1x1 - (x1-x0)*(x1-x0))<<B */
4373 if (tma_mp_lshd (&x1x1, B * 2) != MP_OKAY)
4374 goto X1X1; /* x1x1 = x1x1 << 2*B */
4376 if (tma_mp_add (&x0x0, &t1, &t1) != MP_OKAY)
4377 goto X1X1; /* t1 = x0x0 + t1 */
4378 if (tma_mp_add (&t1, &x1x1, b) != MP_OKAY)
4379 goto X1X1; /* t1 = x0x0 + t1 + x1x1 */
4383 X1X1:tma_mp_clear (&x1x1);
4384 X0X0:tma_mp_clear (&x0x0);
4385 T2:tma_mp_clear (&t2);
4386 T1:tma_mp_clear (&t1);
4387 X1:tma_mp_clear (&x1);
4388 X0:tma_mp_clear (&x0);
4398 /* End: bn_tma_mp_karatsuba_sqr.c */
4400 /* Start: bn_tma_mp_lcm.c */
4403 /* LibTomMath, multiple-precision integer library -- Tom St Denis
4405 * LibTomMath is a library that provides multiple-precision
4406 * integer arithmetic as well as number theoretic functionality.
4408 * The library was designed directly after the MPI library by
4409 * Michael Fromberger but has been written from scratch with
4410 * additional optimizations in place.
4412 * The library is free for all purposes without any express
4413 * guarantee it works.
4415 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
4418 /* computes least common multiple as |a*b|/(a, b) */
4419 int tma_mp_lcm (tma_mp_int * a, tma_mp_int * b, tma_mp_int * c)
4425 if ((res = tma_mp_init_multi (&t1, &t2, NULL)) != MP_OKAY) {
4429 /* t1 = get the GCD of the two inputs */
4430 if ((res = tma_mp_gcd (a, b, &t1)) != MP_OKAY) {
4434 /* divide the smallest by the GCD */
4435 if (tma_mp_cmp_mag(a, b) == MP_LT) {
4436 /* store quotient in t2 such that t2 * b is the LCM */
4437 if ((res = tma_mp_div(a, &t1, &t2, NULL)) != MP_OKAY) {
4440 res = tma_mp_mul(b, &t2, c);
4442 /* store quotient in t2 such that t2 * a is the LCM */
4443 if ((res = tma_mp_div(b, &t1, &t2, NULL)) != MP_OKAY) {
4446 res = tma_mp_mul(a, &t2, c);
4449 /* fix the sign to positive */
4453 tma_mp_clear_multi (&t1, &t2, NULL);
4462 /* End: bn_tma_mp_lcm.c */
4464 /* Start: bn_tma_mp_lshd.c */
4467 /* LibTomMath, multiple-precision integer library -- Tom St Denis
4469 * LibTomMath is a library that provides multiple-precision
4470 * integer arithmetic as well as number theoretic functionality.
4472 * The library was designed directly after the MPI library by
4473 * Michael Fromberger but has been written from scratch with
4474 * additional optimizations in place.
4476 * The library is free for all purposes without any express
4477 * guarantee it works.
4479 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
4482 /* shift left a certain amount of digits */
4483 int tma_mp_lshd (tma_mp_int * a, int b)
4487 /* if its less than zero return */
4492 /* grow to fit the new digits */
4493 if (a->alloc < a->used + b) {
4494 if ((res = tma_mp_grow (a, a->used + b)) != MP_OKAY) {
4500 register tma_mp_digit *top, *bottom;
4502 /* increment the used by the shift amount then copy upwards */
4506 top = a->dp + a->used - 1;
4509 bottom = a->dp + a->used - 1 - b;
4511 /* much like tma_mp_rshd this is implemented using a sliding window
4512 * except the window goes the otherway around. Copying from
4513 * the bottom to the top. see bn_tma_mp_rshd.c for more info.
4515 for (x = a->used - 1; x >= b; x--) {
4519 /* zero the lower digits */
4521 for (x = 0; x < b; x++) {
4533 /* End: bn_tma_mp_lshd.c */
4535 /* Start: bn_tma_mp_mod.c */
4538 /* LibTomMath, multiple-precision integer library -- Tom St Denis
4540 * LibTomMath is a library that provides multiple-precision
4541 * integer arithmetic as well as number theoretic functionality.
4543 * The library was designed directly after the MPI library by
4544 * Michael Fromberger but has been written from scratch with
4545 * additional optimizations in place.
4547 * The library is free for all purposes without any express
4548 * guarantee it works.
4550 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
4553 /* c = a mod b, 0 <= c < b */
4555 tma_mp_mod (tma_mp_int * a, tma_mp_int * b, tma_mp_int * c)
4560 if ((res = tma_mp_init (&t)) != MP_OKAY) {
4564 if ((res = tma_mp_div (a, b, NULL, &t)) != MP_OKAY) {
4569 if (t.sign != b->sign) {
4570 res = tma_mp_add (b, &t, c);
4573 tma_mp_exch (&t, c);
4585 /* End: bn_tma_mp_mod.c */
4587 /* Start: bn_tma_mp_mod_2d.c */
4589 #ifdef BN_MP_MOD_2D_C
4590 /* LibTomMath, multiple-precision integer library -- Tom St Denis
4592 * LibTomMath is a library that provides multiple-precision
4593 * integer arithmetic as well as number theoretic functionality.
4595 * The library was designed directly after the MPI library by
4596 * Michael Fromberger but has been written from scratch with
4597 * additional optimizations in place.
4599 * The library is free for all purposes without any express
4600 * guarantee it works.
4602 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
4605 /* calc a value mod 2**b */
4607 tma_mp_mod_2d (tma_mp_int * a, int b, tma_mp_int * c)
4611 /* if b is <= 0 then zero the int */
4617 /* if the modulus is larger than the value than return */
4618 if (b >= (int) (a->used * DIGIT_BIT)) {
4619 res = tma_mp_copy (a, c);
4624 if ((res = tma_mp_copy (a, c)) != MP_OKAY) {
4628 /* zero digits above the last digit of the modulus */
4629 for (x = (b / DIGIT_BIT) + ((b % DIGIT_BIT) == 0 ? 0 : 1); x < c->used; x++) {
4632 /* clear the digit that is not completely outside/inside the modulus */
4633 c->dp[b / DIGIT_BIT] &=
4634 (tma_mp_digit) ((((tma_mp_digit) 1) << (((tma_mp_digit) b) % DIGIT_BIT)) - ((tma_mp_digit) 1));
4644 /* End: bn_tma_mp_mod_2d.c */
4646 /* Start: bn_tma_mp_mod_d.c */
4648 #ifdef BN_MP_MOD_D_C
4649 /* LibTomMath, multiple-precision integer library -- Tom St Denis
4651 * LibTomMath is a library that provides multiple-precision
4652 * integer arithmetic as well as number theoretic functionality.
4654 * The library was designed directly after the MPI library by
4655 * Michael Fromberger but has been written from scratch with
4656 * additional optimizations in place.
4658 * The library is free for all purposes without any express
4659 * guarantee it works.
4661 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
4665 tma_mp_mod_d (tma_mp_int * a, tma_mp_digit b, tma_mp_digit * c)
4667 return tma_mp_div_d(a, b, NULL, c);
4675 /* End: bn_tma_mp_mod_d.c */
4677 /* Start: bn_tma_mp_montgomery_calc_normalization.c */
4679 #ifdef BN_MP_MONTGOMERY_CALC_NORMALIZATION_C
4680 /* LibTomMath, multiple-precision integer library -- Tom St Denis
4682 * LibTomMath is a library that provides multiple-precision
4683 * integer arithmetic as well as number theoretic functionality.
4685 * The library was designed directly after the MPI library by
4686 * Michael Fromberger but has been written from scratch with
4687 * additional optimizations in place.
4689 * The library is free for all purposes without any express
4690 * guarantee it works.
4692 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
4696 * shifts with subtractions when the result is greater than b.
4698 * The method is slightly modified to shift B unconditionally upto just under
4699 * the leading bit of b. This saves alot of multiple precision shifting.
4701 int tma_mp_montgomery_calc_normalization (tma_mp_int * a, tma_mp_int * b)
4705 /* how many bits of last digit does b use */
4706 bits = tma_mp_count_bits (b) % DIGIT_BIT;
4709 if ((res = tma_mp_2expt (a, (b->used - 1) * DIGIT_BIT + bits - 1)) != MP_OKAY) {
4718 /* now compute C = A * B mod b */
4719 for (x = bits - 1; x < (int)DIGIT_BIT; x++) {
4720 if ((res = tma_mp_mul_2 (a, a)) != MP_OKAY) {
4723 if (tma_mp_cmp_mag (a, b) != MP_LT) {
4724 if ((res = s_tma_mp_sub (a, b, a)) != MP_OKAY) {
4738 /* End: bn_tma_mp_montgomery_calc_normalization.c */
4740 /* Start: bn_tma_mp_montgomery_reduce.c */
4742 #ifdef BN_MP_MONTGOMERY_REDUCE_C
4743 /* LibTomMath, multiple-precision integer library -- Tom St Denis
4745 * LibTomMath is a library that provides multiple-precision
4746 * integer arithmetic as well as number theoretic functionality.
4748 * The library was designed directly after the MPI library by
4749 * Michael Fromberger but has been written from scratch with
4750 * additional optimizations in place.
4752 * The library is free for all purposes without any express
4753 * guarantee it works.
4755 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
4758 /* computes xR**-1 == x (mod N) via Montgomery Reduction */
4760 tma_mp_montgomery_reduce (tma_mp_int * x, tma_mp_int * n, tma_mp_digit rho)
4765 /* can the fast reduction [comba] method be used?
4767 * Note that unlike in mul you're safely allowed *less*
4768 * than the available columns [255 per default] since carries
4769 * are fixed up in the inner loop.
4771 digs = n->used * 2 + 1;
4772 if ((digs < MP_WARRAY) &&
4774 (1 << ((CHAR_BIT * sizeof (tma_mp_word)) - (2 * DIGIT_BIT)))) {
4775 return fast_tma_mp_montgomery_reduce (x, n, rho);
4778 /* grow the input as required */
4779 if (x->alloc < digs) {
4780 if ((res = tma_mp_grow (x, digs)) != MP_OKAY) {
4786 for (ix = 0; ix < n->used; ix++) {
4787 /* mu = ai * rho mod b
4789 * The value of rho must be precalculated via
4790 * montgomery_setup() such that
4791 * it equals -1/n0 mod b this allows the
4792 * following inner loop to reduce the
4793 * input one digit at a time
4795 mu = (tma_mp_digit) (((tma_mp_word)x->dp[ix]) * ((tma_mp_word)rho) & MP_MASK);
4797 /* a = a + mu * m * b**i */
4800 register tma_mp_digit *tmpn, *tmpx, u;
4801 register tma_mp_word r;
4803 /* alias for digits of the modulus */
4806 /* alias for the digits of x [the input] */
4809 /* set the carry to zero */
4812 /* Multiply and add in place */
4813 for (iy = 0; iy < n->used; iy++) {
4814 /* compute product and sum */
4815 r = ((tma_mp_word)mu) * ((tma_mp_word)*tmpn++) +
4816 ((tma_mp_word) u) + ((tma_mp_word) * tmpx);
4819 u = (tma_mp_digit)(r >> ((tma_mp_word) DIGIT_BIT));
4822 *tmpx++ = (tma_mp_digit)(r & ((tma_mp_word) MP_MASK));
4824 /* At this point the ix'th digit of x should be zero */
4827 /* propagate carries upwards as required*/
4830 u = *tmpx >> DIGIT_BIT;
4836 /* at this point the n.used'th least
4837 * significant digits of x are all zero
4838 * which means we can shift x to the
4839 * right by n.used digits and the
4840 * residue is unchanged.
4843 /* x = x/b**n.used */
4845 tma_mp_rshd (x, n->used);
4847 /* if x >= n then x = x - n */
4848 if (tma_mp_cmp_mag (x, n) != MP_LT) {
4849 return s_tma_mp_sub (x, n, x);
4860 /* End: bn_tma_mp_montgomery_reduce.c */
4862 /* Start: bn_tma_mp_montgomery_setup.c */
4864 #ifdef BN_MP_MONTGOMERY_SETUP_C
4865 /* LibTomMath, multiple-precision integer library -- Tom St Denis
4867 * LibTomMath is a library that provides multiple-precision
4868 * integer arithmetic as well as number theoretic functionality.
4870 * The library was designed directly after the MPI library by
4871 * Michael Fromberger but has been written from scratch with
4872 * additional optimizations in place.
4874 * The library is free for all purposes without any express
4875 * guarantee it works.
4877 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
4880 /* setups the montgomery reduction stuff */
4882 tma_mp_montgomery_setup (tma_mp_int * n, tma_mp_digit * rho)
4886 /* fast inversion mod 2**k
4888 * Based on the fact that
4890 * XA = 1 (mod 2**n) => (X(2-XA)) A = 1 (mod 2**2n)
4891 * => 2*X*A - X*X*A*A = 1
4892 * => 2*(1) - (1) = 1
4900 x = (((b + 2) & 4) << 1) + b; /* here x*a==1 mod 2**4 */
4901 x *= 2 - b * x; /* here x*a==1 mod 2**8 */
4902 #if !defined(MP_8BIT)
4903 x *= 2 - b * x; /* here x*a==1 mod 2**16 */
4905 #if defined(MP_64BIT) || !(defined(MP_8BIT) || defined(MP_16BIT))
4906 x *= 2 - b * x; /* here x*a==1 mod 2**32 */
4909 x *= 2 - b * x; /* here x*a==1 mod 2**64 */
4912 /* rho = -1/m mod b */
4913 *rho = (unsigned long)(((tma_mp_word)1 << ((tma_mp_word) DIGIT_BIT)) - x) & MP_MASK;
4923 /* End: bn_tma_mp_montgomery_setup.c */
4925 /* Start: bn_tma_mp_mul.c */
4928 /* LibTomMath, multiple-precision integer library -- Tom St Denis
4930 * LibTomMath is a library that provides multiple-precision
4931 * integer arithmetic as well as number theoretic functionality.
4933 * The library was designed directly after the MPI library by
4934 * Michael Fromberger but has been written from scratch with
4935 * additional optimizations in place.
4937 * The library is free for all purposes without any express
4938 * guarantee it works.
4940 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
4943 /* high level multiplication (handles sign) */
4944 int tma_mp_mul (tma_mp_int * a, tma_mp_int * b, tma_mp_int * c)
4947 neg = (a->sign == b->sign) ? MP_ZPOS : MP_NEG;
4949 /* use Toom-Cook? */
4950 #ifdef BN_MP_TOOM_MUL_C
4951 if (MIN (a->used, b->used) >= TOOM_MUL_CUTOFF) {
4952 res = tma_mp_toom_mul(a, b, c);
4955 #ifdef BN_MP_KARATSUBA_MUL_C
4956 /* use Karatsuba? */
4957 if (MIN (a->used, b->used) >= KARATSUBA_MUL_CUTOFF) {
4958 res = tma_mp_karatsuba_mul (a, b, c);
4962 /* can we use the fast multiplier?
4964 * The fast multiplier can be used if the output will
4965 * have less than MP_WARRAY digits and the number of
4966 * digits won't affect carry propagation
4968 int digs = a->used + b->used + 1;
4970 #ifdef BN_FAST_S_MP_MUL_DIGS_C
4971 if ((digs < MP_WARRAY) &&
4972 MIN(a->used, b->used) <=
4973 (1 << ((CHAR_BIT * sizeof (tma_mp_word)) - (2 * DIGIT_BIT)))) {
4974 res = fast_s_tma_mp_mul_digs (a, b, c, digs);
4977 #ifdef BN_S_MP_MUL_DIGS_C
4978 res = s_tma_mp_mul (a, b, c); /* uses s_tma_mp_mul_digs */
4984 c->sign = (c->used > 0) ? neg : MP_ZPOS;
4993 /* End: bn_tma_mp_mul.c */
4995 /* Start: bn_tma_mp_mul_2.c */
4997 #ifdef BN_MP_MUL_2_C
4998 /* LibTomMath, multiple-precision integer library -- Tom St Denis
5000 * LibTomMath is a library that provides multiple-precision
5001 * integer arithmetic as well as number theoretic functionality.
5003 * The library was designed directly after the MPI library by
5004 * Michael Fromberger but has been written from scratch with
5005 * additional optimizations in place.
5007 * The library is free for all purposes without any express
5008 * guarantee it works.
5010 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
5014 int tma_mp_mul_2(tma_mp_int * a, tma_mp_int * b)
5016 int x, res, oldused;
5018 /* grow to accomodate result */
5019 if (b->alloc < a->used + 1) {
5020 if ((res = tma_mp_grow (b, a->used + 1)) != MP_OKAY) {
5029 register tma_mp_digit r, rr, *tmpa, *tmpb;
5031 /* alias for source */
5034 /* alias for dest */
5039 for (x = 0; x < a->used; x++) {
5041 /* get what will be the *next* carry bit from the
5042 * MSB of the current digit
5044 rr = *tmpa >> ((tma_mp_digit)(DIGIT_BIT - 1));
5046 /* now shift up this digit, add in the carry [from the previous] */
5047 *tmpb++ = ((*tmpa++ << ((tma_mp_digit)1)) | r) & MP_MASK;
5049 /* copy the carry that would be from the source
5050 * digit into the next iteration
5055 /* new leading digit? */
5057 /* add a MSB which is always 1 at this point */
5062 /* now zero any excess digits on the destination
5063 * that we didn't write to
5065 tmpb = b->dp + b->used;
5066 for (x = b->used; x < oldused; x++) {
5079 /* End: bn_tma_mp_mul_2.c */
5081 /* Start: bn_tma_mp_mul_2d.c */
5083 #ifdef BN_MP_MUL_2D_C
5084 /* LibTomMath, multiple-precision integer library -- Tom St Denis
5086 * LibTomMath is a library that provides multiple-precision
5087 * integer arithmetic as well as number theoretic functionality.
5089 * The library was designed directly after the MPI library by
5090 * Michael Fromberger but has been written from scratch with
5091 * additional optimizations in place.
5093 * The library is free for all purposes without any express
5094 * guarantee it works.
5096 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
5099 /* shift left by a certain bit count */
5100 int tma_mp_mul_2d (tma_mp_int * a, int b, tma_mp_int * c)
5107 if ((res = tma_mp_copy (a, c)) != MP_OKAY) {
5112 if (c->alloc < (int)(c->used + b/DIGIT_BIT + 1)) {
5113 if ((res = tma_mp_grow (c, c->used + b / DIGIT_BIT + 1)) != MP_OKAY) {
5118 /* shift by as many digits in the bit count */
5119 if (b >= (int)DIGIT_BIT) {
5120 if ((res = tma_mp_lshd (c, b / DIGIT_BIT)) != MP_OKAY) {
5125 /* shift any bit count < DIGIT_BIT */
5126 d = (tma_mp_digit) (b % DIGIT_BIT);
5128 register tma_mp_digit *tmpc, shift, mask, r, rr;
5131 /* bitmask for carries */
5132 mask = (((tma_mp_digit)1) << d) - 1;
5134 /* shift for msbs */
5135 shift = DIGIT_BIT - d;
5142 for (x = 0; x < c->used; x++) {
5143 /* get the higher bits of the current word */
5144 rr = (*tmpc >> shift) & mask;
5146 /* shift the current word and OR in the carry */
5147 *tmpc = ((*tmpc << d) | r) & MP_MASK;
5150 /* set the carry to the carry bits of the current word */
5154 /* set final carry */
5156 c->dp[(c->used)++] = r;
5168 /* End: bn_tma_mp_mul_2d.c */
5170 /* Start: bn_tma_mp_mul_d.c */
5172 #ifdef BN_MP_MUL_D_C
5173 /* LibTomMath, multiple-precision integer library -- Tom St Denis
5175 * LibTomMath is a library that provides multiple-precision
5176 * integer arithmetic as well as number theoretic functionality.
5178 * The library was designed directly after the MPI library by
5179 * Michael Fromberger but has been written from scratch with
5180 * additional optimizations in place.
5182 * The library is free for all purposes without any express
5183 * guarantee it works.
5185 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
5188 /* multiply by a digit */
5190 tma_mp_mul_d (tma_mp_int * a, tma_mp_digit b, tma_mp_int * c)
5192 tma_mp_digit u, *tmpa, *tmpc;
5194 int ix, res, olduse;
5196 /* make sure c is big enough to hold a*b */
5197 if (c->alloc < a->used + 1) {
5198 if ((res = tma_mp_grow (c, a->used + 1)) != MP_OKAY) {
5203 /* get the original destinations used count */
5209 /* alias for a->dp [source] */
5212 /* alias for c->dp [dest] */
5218 /* compute columns */
5219 for (ix = 0; ix < a->used; ix++) {
5220 /* compute product and carry sum for this term */
5221 r = ((tma_mp_word) u) + ((tma_mp_word)*tmpa++) * ((tma_mp_word)b);
5223 /* mask off higher bits to get a single digit */
5224 *tmpc++ = (tma_mp_digit) (r & ((tma_mp_word) MP_MASK));
5226 /* send carry into next iteration */
5227 u = (tma_mp_digit) (r >> ((tma_mp_word) DIGIT_BIT));
5230 /* store final carry [if any] and increment ix offset */
5234 /* now zero digits above the top */
5235 while (ix++ < olduse) {
5239 /* set used count */
5240 c->used = a->used + 1;
5251 /* End: bn_tma_mp_mul_d.c */
5253 /* Start: bn_tma_mp_mulmod.c */
5255 #ifdef BN_MP_MULMOD_C
5256 /* LibTomMath, multiple-precision integer library -- Tom St Denis
5258 * LibTomMath is a library that provides multiple-precision
5259 * integer arithmetic as well as number theoretic functionality.
5261 * The library was designed directly after the MPI library by
5262 * Michael Fromberger but has been written from scratch with
5263 * additional optimizations in place.
5265 * The library is free for all purposes without any express
5266 * guarantee it works.
5268 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
5271 /* d = a * b (mod c) */
5272 int tma_mp_mulmod (tma_mp_int * a, tma_mp_int * b, tma_mp_int * c, tma_mp_int * d)
5277 if ((res = tma_mp_init (&t)) != MP_OKAY) {
5281 if ((res = tma_mp_mul (a, b, &t)) != MP_OKAY) {
5285 res = tma_mp_mod (&t, c, d);
5295 /* End: bn_tma_mp_mulmod.c */
5297 /* Start: bn_tma_mp_n_root.c */
5299 #ifdef BN_MP_N_ROOT_C
5300 /* LibTomMath, multiple-precision integer library -- Tom St Denis
5302 * LibTomMath is a library that provides multiple-precision
5303 * integer arithmetic as well as number theoretic functionality.
5305 * The library was designed directly after the MPI library by
5306 * Michael Fromberger but has been written from scratch with
5307 * additional optimizations in place.
5309 * The library is free for all purposes without any express
5310 * guarantee it works.
5312 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
5315 /* find the n'th root of an integer
5317 * Result found such that (c)**b <= a and (c+1)**b > a
5319 * This algorithm uses Newton's approximation
5320 * x[i+1] = x[i] - f(x[i])/f'(x[i])
5321 * which will find the root in log(N) time where
5322 * each step involves a fair bit. This is not meant to
5323 * find huge roots [square and cube, etc].
5325 int tma_mp_n_root (tma_mp_int * a, tma_mp_digit b, tma_mp_int * c)
5327 tma_mp_int t1, t2, t3;
5330 /* input must be positive if b is even */
5331 if ((b & 1) == 0 && a->sign == MP_NEG) {
5335 if ((res = tma_mp_init (&t1)) != MP_OKAY) {
5339 if ((res = tma_mp_init (&t2)) != MP_OKAY) {
5343 if ((res = tma_mp_init (&t3)) != MP_OKAY) {
5347 /* if a is negative fudge the sign but keep track */
5352 tma_mp_set (&t2, 2);
5356 if ((res = tma_mp_copy (&t2, &t1)) != MP_OKAY) {
5360 /* t2 = t1 - ((t1**b - a) / (b * t1**(b-1))) */
5362 /* t3 = t1**(b-1) */
5363 if ((res = tma_mp_expt_d (&t1, b - 1, &t3)) != MP_OKAY) {
5369 if ((res = tma_mp_mul (&t3, &t1, &t2)) != MP_OKAY) {
5373 /* t2 = t1**b - a */
5374 if ((res = tma_mp_sub (&t2, a, &t2)) != MP_OKAY) {
5379 /* t3 = t1**(b-1) * b */
5380 if ((res = tma_mp_mul_d (&t3, b, &t3)) != MP_OKAY) {
5384 /* t3 = (t1**b - a)/(b * t1**(b-1)) */
5385 if ((res = tma_mp_div (&t2, &t3, &t3, NULL)) != MP_OKAY) {
5389 if ((res = tma_mp_sub (&t1, &t3, &t2)) != MP_OKAY) {
5392 } while (tma_mp_cmp (&t1, &t2) != MP_EQ);
5394 /* result can be off by a few so check */
5396 if ((res = tma_mp_expt_d (&t1, b, &t2)) != MP_OKAY) {
5400 if (tma_mp_cmp (&t2, a) == MP_GT) {
5401 if ((res = tma_mp_sub_d (&t1, 1, &t1)) != MP_OKAY) {
5409 /* reset the sign of a first */
5412 /* set the result */
5413 tma_mp_exch (&t1, c);
5415 /* set the sign of the result */
5420 LBL_T3:tma_mp_clear (&t3);
5421 LBL_T2:tma_mp_clear (&t2);
5422 LBL_T1:tma_mp_clear (&t1);
5431 /* End: bn_tma_mp_n_root.c */
5433 /* Start: bn_tma_mp_neg.c */
5436 /* LibTomMath, multiple-precision integer library -- Tom St Denis
5438 * LibTomMath is a library that provides multiple-precision
5439 * integer arithmetic as well as number theoretic functionality.
5441 * The library was designed directly after the MPI library by
5442 * Michael Fromberger but has been written from scratch with
5443 * additional optimizations in place.
5445 * The library is free for all purposes without any express
5446 * guarantee it works.
5448 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
5452 int tma_mp_neg (tma_mp_int * a, tma_mp_int * b)
5456 if ((res = tma_mp_copy (a, b)) != MP_OKAY) {
5461 if (tma_mp_iszero(b) != MP_YES) {
5462 b->sign = (a->sign == MP_ZPOS) ? MP_NEG : MP_ZPOS;
5475 /* End: bn_tma_mp_neg.c */
5477 /* Start: bn_tma_mp_or.c */
5480 /* LibTomMath, multiple-precision integer library -- Tom St Denis
5482 * LibTomMath is a library that provides multiple-precision
5483 * integer arithmetic as well as number theoretic functionality.
5485 * The library was designed directly after the MPI library by
5486 * Michael Fromberger but has been written from scratch with
5487 * additional optimizations in place.
5489 * The library is free for all purposes without any express
5490 * guarantee it works.
5492 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
5495 /* OR two ints together */
5496 int tma_mp_or (tma_mp_int * a, tma_mp_int * b, tma_mp_int * c)
5501 if (a->used > b->used) {
5502 if ((res = tma_mp_init_copy (&t, a)) != MP_OKAY) {
5508 if ((res = tma_mp_init_copy (&t, b)) != MP_OKAY) {
5515 for (ix = 0; ix < px; ix++) {
5516 t.dp[ix] |= x->dp[ix];
5519 tma_mp_exch (c, &t);
5529 /* End: bn_tma_mp_or.c */
5531 /* Start: bn_tma_mp_prime_fermat.c */
5533 #ifdef BN_MP_PRIME_FERMAT_C
5534 /* LibTomMath, multiple-precision integer library -- Tom St Denis
5536 * LibTomMath is a library that provides multiple-precision
5537 * integer arithmetic as well as number theoretic functionality.
5539 * The library was designed directly after the MPI library by
5540 * Michael Fromberger but has been written from scratch with
5541 * additional optimizations in place.
5543 * The library is free for all purposes without any express
5544 * guarantee it works.
5546 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
5549 /* performs one Fermat test.
5551 * If "a" were prime then b**a == b (mod a) since the order of
5552 * the multiplicative sub-group would be phi(a) = a-1. That means
5553 * it would be the same as b**(a mod (a-1)) == b**1 == b (mod a).
5555 * Sets result to 1 if the congruence holds, or zero otherwise.
5557 int tma_mp_prime_fermat (tma_mp_int * a, tma_mp_int * b, int *result)
5562 /* default to composite */
5566 if (tma_mp_cmp_d(b, 1) != MP_GT) {
5571 if ((err = tma_mp_init (&t)) != MP_OKAY) {
5575 /* compute t = b**a mod a */
5576 if ((err = tma_mp_exptmod (b, a, a, &t)) != MP_OKAY) {
5580 /* is it equal to b? */
5581 if (tma_mp_cmp (&t, b) == MP_EQ) {
5586 LBL_T:tma_mp_clear (&t);
5595 /* End: bn_tma_mp_prime_fermat.c */
5597 /* Start: bn_tma_mp_prime_is_divisible.c */
5599 #ifdef BN_MP_PRIME_IS_DIVISIBLE_C
5600 /* LibTomMath, multiple-precision integer library -- Tom St Denis
5602 * LibTomMath is a library that provides multiple-precision
5603 * integer arithmetic as well as number theoretic functionality.
5605 * The library was designed directly after the MPI library by
5606 * Michael Fromberger but has been written from scratch with
5607 * additional optimizations in place.
5609 * The library is free for all purposes without any express
5610 * guarantee it works.
5612 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
5615 /* determines if an integers is divisible by one
5616 * of the first PRIME_SIZE primes or not
5618 * sets result to 0 if not, 1 if yes
5620 int tma_mp_prime_is_divisible (tma_mp_int * a, int *result)
5625 /* default to not */
5628 for (ix = 0; ix < PRIME_SIZE; ix++) {
5629 /* what is a mod LBL_prime_tab[ix] */
5630 if ((err = tma_mp_mod_d (a, ltm_prime_tab[ix], &res)) != MP_OKAY) {
5634 /* is the residue zero? */
5649 /* End: bn_tma_mp_prime_is_divisible.c */
5651 /* Start: bn_tma_mp_prime_is_prime.c */
5653 #ifdef BN_MP_PRIME_IS_PRIME_C
5654 /* LibTomMath, multiple-precision integer library -- Tom St Denis
5656 * LibTomMath is a library that provides multiple-precision
5657 * integer arithmetic as well as number theoretic functionality.
5659 * The library was designed directly after the MPI library by
5660 * Michael Fromberger but has been written from scratch with
5661 * additional optimizations in place.
5663 * The library is free for all purposes without any express
5664 * guarantee it works.
5666 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
5669 /* performs a variable number of rounds of Miller-Rabin
5671 * Probability of error after t rounds is no more than
5674 * Sets result to 1 if probably prime, 0 otherwise
5676 int tma_mp_prime_is_prime (tma_mp_int * a, int t, int *result)
5684 /* valid value of t? */
5685 if (t <= 0 || t > PRIME_SIZE) {
5689 /* is the input equal to one of the primes in the table? */
5690 for (ix = 0; ix < PRIME_SIZE; ix++) {
5691 if (tma_mp_cmp_d(a, ltm_prime_tab[ix]) == MP_EQ) {
5697 /* first perform trial division */
5698 if ((err = tma_mp_prime_is_divisible (a, &res)) != MP_OKAY) {
5702 /* return if it was trivially divisible */
5703 if (res == MP_YES) {
5707 /* now perform the miller-rabin rounds */
5708 if ((err = tma_mp_init (&b)) != MP_OKAY) {
5712 for (ix = 0; ix < t; ix++) {
5714 tma_mp_set (&b, ltm_prime_tab[ix]);
5716 if ((err = tma_mp_prime_miller_rabin (a, &b, &res)) != MP_OKAY) {
5725 /* passed the test */
5727 LBL_B:tma_mp_clear (&b);
5736 /* End: bn_tma_mp_prime_is_prime.c */
5738 /* Start: bn_tma_mp_prime_miller_rabin.c */
5740 #ifdef BN_MP_PRIME_MILLER_RABIN_C
5741 /* LibTomMath, multiple-precision integer library -- Tom St Denis
5743 * LibTomMath is a library that provides multiple-precision
5744 * integer arithmetic as well as number theoretic functionality.
5746 * The library was designed directly after the MPI library by
5747 * Michael Fromberger but has been written from scratch with
5748 * additional optimizations in place.
5750 * The library is free for all purposes without any express
5751 * guarantee it works.
5753 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
5756 /* Miller-Rabin test of "a" to the base of "b" as described in
5757 * HAC pp. 139 Algorithm 4.24
5759 * Sets result to 0 if definitely composite or 1 if probably prime.
5760 * Randomly the chance of error is no more than 1/4 and often
5763 int tma_mp_prime_miller_rabin (tma_mp_int * a, tma_mp_int * b, int *result)
5765 tma_mp_int n1, y, r;
5772 if (tma_mp_cmp_d(b, 1) != MP_GT) {
5776 /* get n1 = a - 1 */
5777 if ((err = tma_mp_init_copy (&n1, a)) != MP_OKAY) {
5780 if ((err = tma_mp_sub_d (&n1, 1, &n1)) != MP_OKAY) {
5784 /* set 2**s * r = n1 */
5785 if ((err = tma_mp_init_copy (&r, &n1)) != MP_OKAY) {
5789 /* count the number of least significant bits
5792 s = tma_mp_cnt_lsb(&r);
5794 /* now divide n - 1 by 2**s */
5795 if ((err = tma_mp_div_2d (&r, s, &r, NULL)) != MP_OKAY) {
5799 /* compute y = b**r mod a */
5800 if ((err = tma_mp_init (&y)) != MP_OKAY) {
5803 if ((err = tma_mp_exptmod (b, &r, a, &y)) != MP_OKAY) {
5807 /* if y != 1 and y != n1 do */
5808 if (tma_mp_cmp_d (&y, 1) != MP_EQ && tma_mp_cmp (&y, &n1) != MP_EQ) {
5810 /* while j <= s-1 and y != n1 */
5811 while ((j <= (s - 1)) && tma_mp_cmp (&y, &n1) != MP_EQ) {
5812 if ((err = tma_mp_sqrmod (&y, a, &y)) != MP_OKAY) {
5816 /* if y == 1 then composite */
5817 if (tma_mp_cmp_d (&y, 1) == MP_EQ) {
5824 /* if y != n1 then composite */
5825 if (tma_mp_cmp (&y, &n1) != MP_EQ) {
5830 /* probably prime now */
5832 LBL_Y:tma_mp_clear (&y);
5833 LBL_R:tma_mp_clear (&r);
5834 LBL_N1:tma_mp_clear (&n1);
5843 /* End: bn_tma_mp_prime_miller_rabin.c */
5845 /* Start: bn_tma_mp_prime_next_prime.c */
5847 #ifdef BN_MP_PRIME_NEXT_PRIME_C
5848 /* LibTomMath, multiple-precision integer library -- Tom St Denis
5850 * LibTomMath is a library that provides multiple-precision
5851 * integer arithmetic as well as number theoretic functionality.
5853 * The library was designed directly after the MPI library by
5854 * Michael Fromberger but has been written from scratch with
5855 * additional optimizations in place.
5857 * The library is free for all purposes without any express
5858 * guarantee it works.
5860 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
5863 /* finds the next prime after the number "a" using "t" trials
5866 * bbs_style = 1 means the prime must be congruent to 3 mod 4
5868 int tma_mp_prime_next_prime(tma_mp_int *a, int t, int bbs_style)
5871 tma_mp_digit res_tab[PRIME_SIZE], step, kstep;
5874 /* ensure t is valid */
5875 if (t <= 0 || t > PRIME_SIZE) {
5879 /* force positive */
5882 /* simple algo if a is less than the largest prime in the table */
5883 if (tma_mp_cmp_d(a, ltm_prime_tab[PRIME_SIZE-1]) == MP_LT) {
5884 /* find which prime it is bigger than */
5885 for (x = PRIME_SIZE - 2; x >= 0; x--) {
5886 if (tma_mp_cmp_d(a, ltm_prime_tab[x]) != MP_LT) {
5887 if (bbs_style == 1) {
5888 /* ok we found a prime smaller or
5889 * equal [so the next is larger]
5891 * however, the prime must be
5892 * congruent to 3 mod 4
5894 if ((ltm_prime_tab[x + 1] & 3) != 3) {
5895 /* scan upwards for a prime congruent to 3 mod 4 */
5896 for (y = x + 1; y < PRIME_SIZE; y++) {
5897 if ((ltm_prime_tab[y] & 3) == 3) {
5898 tma_mp_set(a, ltm_prime_tab[y]);
5904 tma_mp_set(a, ltm_prime_tab[x + 1]);
5909 /* at this point a maybe 1 */
5910 if (tma_mp_cmp_d(a, 1) == MP_EQ) {
5914 /* fall through to the sieve */
5917 /* generate a prime congruent to 3 mod 4 or 1/3 mod 4? */
5918 if (bbs_style == 1) {
5924 /* at this point we will use a combination of a sieve and Miller-Rabin */
5926 if (bbs_style == 1) {
5927 /* if a mod 4 != 3 subtract the correct value to make it so */
5928 if ((a->dp[0] & 3) != 3) {
5929 if ((err = tma_mp_sub_d(a, (a->dp[0] & 3) + 1, a)) != MP_OKAY) { return err; };
5932 if (tma_mp_iseven(a) == 1) {
5934 if ((err = tma_mp_sub_d(a, 1, a)) != MP_OKAY) {
5940 /* generate the restable */
5941 for (x = 1; x < PRIME_SIZE; x++) {
5942 if ((err = tma_mp_mod_d(a, ltm_prime_tab[x], res_tab + x)) != MP_OKAY) {
5947 /* init temp used for Miller-Rabin Testing */
5948 if ((err = tma_mp_init(&b)) != MP_OKAY) {
5953 /* skip to the next non-trivially divisible candidate */
5956 /* y == 1 if any residue was zero [e.g. cannot be prime] */
5959 /* increase step to next candidate */
5962 /* compute the new residue without using division */
5963 for (x = 1; x < PRIME_SIZE; x++) {
5964 /* add the step to each residue */
5965 res_tab[x] += kstep;
5967 /* subtract the modulus [instead of using division] */
5968 if (res_tab[x] >= ltm_prime_tab[x]) {
5969 res_tab[x] -= ltm_prime_tab[x];
5972 /* set flag if zero */
5973 if (res_tab[x] == 0) {
5977 } while (y == 1 && step < ((((tma_mp_digit)1)<<DIGIT_BIT) - kstep));
5980 if ((err = tma_mp_add_d(a, step, a)) != MP_OKAY) {
5984 /* if didn't pass sieve and step == MAX then skip test */
5985 if (y == 1 && step >= ((((tma_mp_digit)1)<<DIGIT_BIT) - kstep)) {
5989 /* is this prime? */
5990 for (x = 0; x < t; x++) {
5991 tma_mp_set(&b, ltm_prime_tab[t]);
5992 if ((err = tma_mp_prime_miller_rabin(a, &b, &res)) != MP_OKAY) {
6000 if (res == MP_YES) {
6017 /* End: bn_tma_mp_prime_next_prime.c */
6019 /* Start: bn_tma_mp_prime_rabin_miller_trials.c */
6021 #ifdef BN_MP_PRIME_RABIN_MILLER_TRIALS_C
6022 /* LibTomMath, multiple-precision integer library -- Tom St Denis
6024 * LibTomMath is a library that provides multiple-precision
6025 * integer arithmetic as well as number theoretic functionality.
6027 * The library was designed directly after the MPI library by
6028 * Michael Fromberger but has been written from scratch with
6029 * additional optimizations in place.
6031 * The library is free for all purposes without any express
6032 * guarantee it works.
6034 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
6038 static const struct {
6051 /* returns # of RM trials required for a given bit size */
6052 int tma_mp_prime_rabin_miller_trials(int size)
6056 for (x = 0; x < (int)(sizeof(sizes)/(sizeof(sizes[0]))); x++) {
6057 if (sizes[x].k == size) {
6059 } else if (sizes[x].k > size) {
6060 return (x == 0) ? sizes[0].t : sizes[x - 1].t;
6063 return sizes[x-1].t + 1;
6073 /* End: bn_tma_mp_prime_rabin_miller_trials.c */
6075 /* Start: bn_tma_mp_prime_random_ex.c */
6077 #ifdef BN_MP_PRIME_RANDOM_EX_C
6078 /* LibTomMath, multiple-precision integer library -- Tom St Denis
6080 * LibTomMath is a library that provides multiple-precision
6081 * integer arithmetic as well as number theoretic functionality.
6083 * The library was designed directly after the MPI library by
6084 * Michael Fromberger but has been written from scratch with
6085 * additional optimizations in place.
6087 * The library is free for all purposes without any express
6088 * guarantee it works.
6090 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
6093 /* makes a truly random prime of a given size (bits),
6095 * Flags are as follows:
6097 * LTM_PRIME_BBS - make prime congruent to 3 mod 4
6098 * LTM_PRIME_SAFE - make sure (p-1)/2 is prime as well (implies LTM_PRIME_BBS)
6099 * LTM_PRIME_2MSB_OFF - make the 2nd highest bit zero
6100 * LTM_PRIME_2MSB_ON - make the 2nd highest bit one
6102 * You have to supply a callback which fills in a buffer with random bytes. "dat" is a parameter you can
6103 * have passed to the callback (e.g. a state or something). This function doesn't use "dat" itself
6108 /* This is possibly the mother of all prime generation functions, muahahahahaha! */
6109 int tma_mp_prime_random_ex(tma_mp_int *a, int t, int size, int flags, ltm_prime_callback cb, void *dat)
6111 unsigned char *tmp, maskAND, maskOR_msb, maskOR_lsb;
6112 int res, err, bsize, maskOR_msb_offset;
6114 /* sanity check the input */
6115 if (size <= 1 || t <= 0) {
6119 /* LTM_PRIME_SAFE implies LTM_PRIME_BBS */
6120 if (flags & LTM_PRIME_SAFE) {
6121 flags |= LTM_PRIME_BBS;
6124 /* calc the byte size */
6125 bsize = (size>>3) + ((size&7)?1:0);
6127 /* we need a buffer of bsize bytes */
6128 tmp = OPT_CAST(unsigned char) XMALLOC(bsize);
6133 /* calc the maskAND value for the MSbyte*/
6134 maskAND = ((size&7) == 0) ? 0xFF : (0xFF >> (8 - (size & 7)));
6136 /* calc the maskOR_msb */
6138 maskOR_msb_offset = ((size & 7) == 1) ? 1 : 0;
6139 if (flags & LTM_PRIME_2MSB_ON) {
6140 maskOR_msb |= 0x80 >> ((9 - size) & 7);
6143 /* get the maskOR_lsb */
6145 if (flags & LTM_PRIME_BBS) {
6150 /* read the bytes */
6151 if (cb(tmp, bsize, dat) != bsize) {
6156 /* work over the MSbyte */
6158 tmp[0] |= 1 << ((size - 1) & 7);
6160 /* mix in the maskORs */
6161 tmp[maskOR_msb_offset] |= maskOR_msb;
6162 tmp[bsize-1] |= maskOR_lsb;
6165 if ((err = tma_mp_read_unsigned_bin(a, tmp, bsize)) != MP_OKAY) { goto error; }
6168 if ((err = tma_mp_prime_is_prime(a, t, &res)) != MP_OKAY) { goto error; }
6173 if (flags & LTM_PRIME_SAFE) {
6174 /* see if (a-1)/2 is prime */
6175 if ((err = tma_mp_sub_d(a, 1, a)) != MP_OKAY) { goto error; }
6176 if ((err = tma_mp_div_2(a, a)) != MP_OKAY) { goto error; }
6179 if ((err = tma_mp_prime_is_prime(a, t, &res)) != MP_OKAY) { goto error; }
6181 } while (res == MP_NO);
6183 if (flags & LTM_PRIME_SAFE) {
6184 /* restore a to the original value */
6185 if ((err = tma_mp_mul_2(a, a)) != MP_OKAY) { goto error; }
6186 if ((err = tma_mp_add_d(a, 1, a)) != MP_OKAY) { goto error; }
6202 /* End: bn_tma_mp_prime_random_ex.c */
6204 /* Start: bn_tma_mp_radix_size.c */
6206 #ifdef BN_MP_RADIX_SIZE_C
6207 /* LibTomMath, multiple-precision integer library -- Tom St Denis
6209 * LibTomMath is a library that provides multiple-precision
6210 * integer arithmetic as well as number theoretic functionality.
6212 * The library was designed directly after the MPI library by
6213 * Michael Fromberger but has been written from scratch with
6214 * additional optimizations in place.
6216 * The library is free for all purposes without any express
6217 * guarantee it works.
6219 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
6222 /* returns size of ASCII reprensentation */
6223 int tma_mp_radix_size (tma_mp_int * a, int radix, int *size)
6231 /* special case for binary */
6233 *size = tma_mp_count_bits (a) + (a->sign == MP_NEG ? 1 : 0) + 1;
6237 /* make sure the radix is in range */
6238 if (radix < 2 || radix > 64) {
6242 if (tma_mp_iszero(a) == MP_YES) {
6247 /* digs is the digit count */
6250 /* if it's negative add one for the sign */
6251 if (a->sign == MP_NEG) {
6255 /* init a copy of the input */
6256 if ((res = tma_mp_init_copy (&t, a)) != MP_OKAY) {
6260 /* force temp to positive */
6263 /* fetch out all of the digits */
6264 while (tma_mp_iszero (&t) == MP_NO) {
6265 if ((res = tma_mp_div_d (&t, (tma_mp_digit) radix, &t, &d)) != MP_OKAY) {
6273 /* return digs + 1, the 1 is for the NULL byte that would be required. */
6284 /* End: bn_tma_mp_radix_size.c */
6286 /* Start: bn_tma_mp_radix_smap.c */
6288 #ifdef BN_MP_RADIX_SMAP_C
6289 /* LibTomMath, multiple-precision integer library -- Tom St Denis
6291 * LibTomMath is a library that provides multiple-precision
6292 * integer arithmetic as well as number theoretic functionality.
6294 * The library was designed directly after the MPI library by
6295 * Michael Fromberger but has been written from scratch with
6296 * additional optimizations in place.
6298 * The library is free for all purposes without any express
6299 * guarantee it works.
6301 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
6304 /* chars used in radix conversions */
6305 const char *tma_mp_s_rmap = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz+/";
6312 /* End: bn_tma_mp_radix_smap.c */
6314 /* Start: bn_tma_mp_rand.c */
6317 /* LibTomMath, multiple-precision integer library -- Tom St Denis
6319 * LibTomMath is a library that provides multiple-precision
6320 * integer arithmetic as well as number theoretic functionality.
6322 * The library was designed directly after the MPI library by
6323 * Michael Fromberger but has been written from scratch with
6324 * additional optimizations in place.
6326 * The library is free for all purposes without any express
6327 * guarantee it works.
6329 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
6332 /* makes a pseudo-random int of a given size */
6334 tma_mp_rand (tma_mp_int * a, int digits)
6344 /* first place a random non-zero digit */
6346 d = ((tma_mp_digit) abs (rand ())) & MP_MASK;
6349 if ((res = tma_mp_add_d (a, d, a)) != MP_OKAY) {
6353 while (--digits > 0) {
6354 if ((res = tma_mp_lshd (a, 1)) != MP_OKAY) {
6358 if ((res = tma_mp_add_d (a, ((tma_mp_digit) abs (rand ())), a)) != MP_OKAY) {
6371 /* End: bn_tma_mp_rand.c */
6373 /* Start: bn_tma_mp_read_radix.c */
6375 #ifdef BN_MP_READ_RADIX_C
6376 /* LibTomMath, multiple-precision integer library -- Tom St Denis
6378 * LibTomMath is a library that provides multiple-precision
6379 * integer arithmetic as well as number theoretic functionality.
6381 * The library was designed directly after the MPI library by
6382 * Michael Fromberger but has been written from scratch with
6383 * additional optimizations in place.
6385 * The library is free for all purposes without any express
6386 * guarantee it works.
6388 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
6391 /* read a string [ASCII] in a given radix */
6392 int tma_mp_read_radix (tma_mp_int * a, const char *str, int radix)
6397 /* zero the digit bignum */
6400 /* make sure the radix is ok */
6401 if (radix < 2 || radix > 64) {
6405 /* if the leading digit is a
6406 * minus set the sign to negative.
6415 /* set the integer to the default of zero */
6418 /* process each digit of the string */
6420 /* if the radix < 36 the conversion is case insensitive
6421 * this allows numbers like 1AB and 1ab to represent the same value
6424 ch = (char) ((radix < 36) ? toupper (*str) : *str);
6425 for (y = 0; y < 64; y++) {
6426 if (ch == tma_mp_s_rmap[y]) {
6431 /* if the char was found in the map
6432 * and is less than the given radix add it
6433 * to the number, otherwise exit the loop.
6436 if ((res = tma_mp_mul_d (a, (tma_mp_digit) radix, a)) != MP_OKAY) {
6439 if ((res = tma_mp_add_d (a, (tma_mp_digit) y, a)) != MP_OKAY) {
6448 /* set the sign only if a != 0 */
6449 if (tma_mp_iszero(a) != 1) {
6460 /* End: bn_tma_mp_read_radix.c */
6462 /* Start: bn_tma_mp_read_signed_bin.c */
6464 #ifdef BN_MP_READ_SIGNED_BIN_C
6465 /* LibTomMath, multiple-precision integer library -- Tom St Denis
6467 * LibTomMath is a library that provides multiple-precision
6468 * integer arithmetic as well as number theoretic functionality.
6470 * The library was designed directly after the MPI library by
6471 * Michael Fromberger but has been written from scratch with
6472 * additional optimizations in place.
6474 * The library is free for all purposes without any express
6475 * guarantee it works.
6477 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
6480 /* read signed bin, big endian, first byte is 0==positive or 1==negative */
6481 int tma_mp_read_signed_bin (tma_mp_int * a, const unsigned char *b, int c)
6485 /* read magnitude */
6486 if ((res = tma_mp_read_unsigned_bin (a, b + 1, c - 1)) != MP_OKAY) {
6490 /* first byte is 0 for positive, non-zero for negative */
6505 /* End: bn_tma_mp_read_signed_bin.c */
6507 /* Start: bn_tma_mp_read_unsigned_bin.c */
6509 #ifdef BN_MP_READ_UNSIGNED_BIN_C
6510 /* LibTomMath, multiple-precision integer library -- Tom St Denis
6512 * LibTomMath is a library that provides multiple-precision
6513 * integer arithmetic as well as number theoretic functionality.
6515 * The library was designed directly after the MPI library by
6516 * Michael Fromberger but has been written from scratch with
6517 * additional optimizations in place.
6519 * The library is free for all purposes without any express
6520 * guarantee it works.
6522 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
6525 /* reads a unsigned char array, assumes the msb is stored first [big endian] */
6526 int tma_mp_read_unsigned_bin (tma_mp_int * a, const unsigned char *b, int c)
6530 /* make sure there are at least two digits */
6532 if ((res = tma_mp_grow(a, 2)) != MP_OKAY) {
6540 /* read the bytes in */
6542 if ((res = tma_mp_mul_2d (a, 8, a)) != MP_OKAY) {
6550 a->dp[0] = (*b & MP_MASK);
6551 a->dp[1] |= ((*b++ >> 7U) & 1);
6564 /* End: bn_tma_mp_read_unsigned_bin.c */
6566 /* Start: bn_tma_mp_reduce.c */
6568 #ifdef BN_MP_REDUCE_C
6569 /* LibTomMath, multiple-precision integer library -- Tom St Denis
6571 * LibTomMath is a library that provides multiple-precision
6572 * integer arithmetic as well as number theoretic functionality.
6574 * The library was designed directly after the MPI library by
6575 * Michael Fromberger but has been written from scratch with
6576 * additional optimizations in place.
6578 * The library is free for all purposes without any express
6579 * guarantee it works.
6581 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
6584 /* reduces x mod m, assumes 0 < x < m**2, mu is
6585 * precomputed via tma_mp_reduce_setup.
6586 * From HAC pp.604 Algorithm 14.42
6588 int tma_mp_reduce (tma_mp_int * x, tma_mp_int * m, tma_mp_int * mu)
6591 int res, um = m->used;
6594 if ((res = tma_mp_init_copy (&q, x)) != MP_OKAY) {
6598 /* q1 = x / b**(k-1) */
6599 tma_mp_rshd (&q, um - 1);
6601 /* according to HAC this optimization is ok */
6602 if (((unsigned long) um) > (((tma_mp_digit)1) << (DIGIT_BIT - 1))) {
6603 if ((res = tma_mp_mul (&q, mu, &q)) != MP_OKAY) {
6607 #ifdef BN_S_MP_MUL_HIGH_DIGS_C
6608 if ((res = s_tma_mp_mul_high_digs (&q, mu, &q, um)) != MP_OKAY) {
6611 #elif defined(BN_FAST_S_MP_MUL_HIGH_DIGS_C)
6612 if ((res = fast_s_tma_mp_mul_high_digs (&q, mu, &q, um)) != MP_OKAY) {
6623 /* q3 = q2 / b**(k+1) */
6624 tma_mp_rshd (&q, um + 1);
6626 /* x = x mod b**(k+1), quick (no division) */
6627 if ((res = tma_mp_mod_2d (x, DIGIT_BIT * (um + 1), x)) != MP_OKAY) {
6631 /* q = q * m mod b**(k+1), quick (no division) */
6632 if ((res = s_tma_mp_mul_digs (&q, m, &q, um + 1)) != MP_OKAY) {
6637 if ((res = tma_mp_sub (x, &q, x)) != MP_OKAY) {
6641 /* If x < 0, add b**(k+1) to it */
6642 if (tma_mp_cmp_d (x, 0) == MP_LT) {
6644 if ((res = tma_mp_lshd (&q, um + 1)) != MP_OKAY)
6646 if ((res = tma_mp_add (x, &q, x)) != MP_OKAY)
6650 /* Back off if it's too big */
6651 while (tma_mp_cmp (x, m) != MP_LT) {
6652 if ((res = s_tma_mp_sub (x, m, x)) != MP_OKAY) {
6668 /* End: bn_tma_mp_reduce.c */
6670 /* Start: bn_tma_mp_reduce_2k.c */
6672 #ifdef BN_MP_REDUCE_2K_C
6673 /* LibTomMath, multiple-precision integer library -- Tom St Denis
6675 * LibTomMath is a library that provides multiple-precision
6676 * integer arithmetic as well as number theoretic functionality.
6678 * The library was designed directly after the MPI library by
6679 * Michael Fromberger but has been written from scratch with
6680 * additional optimizations in place.
6682 * The library is free for all purposes without any express
6683 * guarantee it works.
6685 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
6688 /* reduces a modulo n where n is of the form 2**p - d */
6689 int tma_mp_reduce_2k(tma_mp_int *a, tma_mp_int *n, tma_mp_digit d)
6694 if ((res = tma_mp_init(&q)) != MP_OKAY) {
6698 p = tma_mp_count_bits(n);
6700 /* q = a/2**p, a = a mod 2**p */
6701 if ((res = tma_mp_div_2d(a, p, &q, a)) != MP_OKAY) {
6707 if ((res = tma_mp_mul_d(&q, d, &q)) != MP_OKAY) {
6713 if ((res = s_tma_mp_add(a, &q, a)) != MP_OKAY) {
6717 if (tma_mp_cmp_mag(a, n) != MP_LT) {
6718 s_tma_mp_sub(a, n, a);
6733 /* End: bn_tma_mp_reduce_2k.c */
6735 /* Start: bn_tma_mp_reduce_2k_l.c */
6737 #ifdef BN_MP_REDUCE_2K_L_C
6738 /* LibTomMath, multiple-precision integer library -- Tom St Denis
6740 * LibTomMath is a library that provides multiple-precision
6741 * integer arithmetic as well as number theoretic functionality.
6743 * The library was designed directly after the MPI library by
6744 * Michael Fromberger but has been written from scratch with
6745 * additional optimizations in place.
6747 * The library is free for all purposes without any express
6748 * guarantee it works.
6750 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
6753 /* reduces a modulo n where n is of the form 2**p - d
6754 This differs from reduce_2k since "d" can be larger
6755 than a single digit.
6757 int tma_mp_reduce_2k_l(tma_mp_int *a, tma_mp_int *n, tma_mp_int *d)
6762 if ((res = tma_mp_init(&q)) != MP_OKAY) {
6766 p = tma_mp_count_bits(n);
6768 /* q = a/2**p, a = a mod 2**p */
6769 if ((res = tma_mp_div_2d(a, p, &q, a)) != MP_OKAY) {
6774 if ((res = tma_mp_mul(&q, d, &q)) != MP_OKAY) {
6779 if ((res = s_tma_mp_add(a, &q, a)) != MP_OKAY) {
6783 if (tma_mp_cmp_mag(a, n) != MP_LT) {
6784 s_tma_mp_sub(a, n, a);
6799 /* End: bn_tma_mp_reduce_2k_l.c */
6801 /* Start: bn_tma_mp_reduce_2k_setup.c */
6803 #ifdef BN_MP_REDUCE_2K_SETUP_C
6804 /* LibTomMath, multiple-precision integer library -- Tom St Denis
6806 * LibTomMath is a library that provides multiple-precision
6807 * integer arithmetic as well as number theoretic functionality.
6809 * The library was designed directly after the MPI library by
6810 * Michael Fromberger but has been written from scratch with
6811 * additional optimizations in place.
6813 * The library is free for all purposes without any express
6814 * guarantee it works.
6816 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
6819 /* determines the setup value */
6820 int tma_mp_reduce_2k_setup(tma_mp_int *a, tma_mp_digit *d)
6825 if ((res = tma_mp_init(&tmp)) != MP_OKAY) {
6829 p = tma_mp_count_bits(a);
6830 if ((res = tma_mp_2expt(&tmp, p)) != MP_OKAY) {
6835 if ((res = s_tma_mp_sub(&tmp, a, &tmp)) != MP_OKAY) {
6850 /* End: bn_tma_mp_reduce_2k_setup.c */
6852 /* Start: bn_tma_mp_reduce_2k_setup_l.c */
6854 #ifdef BN_MP_REDUCE_2K_SETUP_L_C
6855 /* LibTomMath, multiple-precision integer library -- Tom St Denis
6857 * LibTomMath is a library that provides multiple-precision
6858 * integer arithmetic as well as number theoretic functionality.
6860 * The library was designed directly after the MPI library by
6861 * Michael Fromberger but has been written from scratch with
6862 * additional optimizations in place.
6864 * The library is free for all purposes without any express
6865 * guarantee it works.
6867 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
6870 /* determines the setup value */
6871 int tma_mp_reduce_2k_setup_l(tma_mp_int *a, tma_mp_int *d)
6876 if ((res = tma_mp_init(&tmp)) != MP_OKAY) {
6880 if ((res = tma_mp_2expt(&tmp, tma_mp_count_bits(a))) != MP_OKAY) {
6884 if ((res = s_tma_mp_sub(&tmp, a, d)) != MP_OKAY) {
6898 /* End: bn_tma_mp_reduce_2k_setup_l.c */
6900 /* Start: bn_tma_mp_reduce_is_2k.c */
6902 #ifdef BN_MP_REDUCE_IS_2K_C
6903 /* LibTomMath, multiple-precision integer library -- Tom St Denis
6905 * LibTomMath is a library that provides multiple-precision
6906 * integer arithmetic as well as number theoretic functionality.
6908 * The library was designed directly after the MPI library by
6909 * Michael Fromberger but has been written from scratch with
6910 * additional optimizations in place.
6912 * The library is free for all purposes without any express
6913 * guarantee it works.
6915 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
6918 /* determines if tma_mp_reduce_2k can be used */
6919 int tma_mp_reduce_is_2k(tma_mp_int *a)
6926 } else if (a->used == 1) {
6928 } else if (a->used > 1) {
6929 iy = tma_mp_count_bits(a);
6933 /* Test every bit from the second digit up, must be 1 */
6934 for (ix = DIGIT_BIT; ix < iy; ix++) {
6935 if ((a->dp[iw] & iz) == 0) {
6939 if (iz > (tma_mp_digit)MP_MASK) {
6954 /* End: bn_tma_mp_reduce_is_2k.c */
6956 /* Start: bn_tma_mp_reduce_is_2k_l.c */
6958 #ifdef BN_MP_REDUCE_IS_2K_L_C
6959 /* LibTomMath, multiple-precision integer library -- Tom St Denis
6961 * LibTomMath is a library that provides multiple-precision
6962 * integer arithmetic as well as number theoretic functionality.
6964 * The library was designed directly after the MPI library by
6965 * Michael Fromberger but has been written from scratch with
6966 * additional optimizations in place.
6968 * The library is free for all purposes without any express
6969 * guarantee it works.
6971 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
6974 /* determines if reduce_2k_l can be used */
6975 int tma_mp_reduce_is_2k_l(tma_mp_int *a)
6981 } else if (a->used == 1) {
6983 } else if (a->used > 1) {
6984 /* if more than half of the digits are -1 we're sold */
6985 for (iy = ix = 0; ix < a->used; ix++) {
6986 if (a->dp[ix] == MP_MASK) {
6990 return (iy >= (a->used/2)) ? MP_YES : MP_NO;
7002 /* End: bn_tma_mp_reduce_is_2k_l.c */
7004 /* Start: bn_tma_mp_reduce_setup.c */
7006 #ifdef BN_MP_REDUCE_SETUP_C
7007 /* LibTomMath, multiple-precision integer library -- Tom St Denis
7009 * LibTomMath is a library that provides multiple-precision
7010 * integer arithmetic as well as number theoretic functionality.
7012 * The library was designed directly after the MPI library by
7013 * Michael Fromberger but has been written from scratch with
7014 * additional optimizations in place.
7016 * The library is free for all purposes without any express
7017 * guarantee it works.
7019 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
7022 /* pre-calculate the value required for Barrett reduction
7023 * For a given modulus "b" it calulates the value required in "a"
7025 int tma_mp_reduce_setup (tma_mp_int * a, tma_mp_int * b)
7029 if ((res = tma_mp_2expt (a, b->used * 2 * DIGIT_BIT)) != MP_OKAY) {
7032 return tma_mp_div (a, b, a, NULL);
7040 /* End: bn_tma_mp_reduce_setup.c */
7042 /* Start: bn_tma_mp_rshd.c */
7045 /* LibTomMath, multiple-precision integer library -- Tom St Denis
7047 * LibTomMath is a library that provides multiple-precision
7048 * integer arithmetic as well as number theoretic functionality.
7050 * The library was designed directly after the MPI library by
7051 * Michael Fromberger but has been written from scratch with
7052 * additional optimizations in place.
7054 * The library is free for all purposes without any express
7055 * guarantee it works.
7057 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
7060 /* shift right a certain amount of digits */
7061 void tma_mp_rshd (tma_mp_int * a, int b)
7065 /* if b <= 0 then ignore it */
7070 /* if b > used then simply zero it and return */
7077 register tma_mp_digit *bottom, *top;
7079 /* shift the digits down */
7084 /* top [offset into digits] */
7087 /* this is implemented as a sliding window where
7088 * the window is b-digits long and digits from
7089 * the top of the window are copied to the bottom
7093 b-2 | b-1 | b0 | b1 | b2 | ... | bb | ---->
7095 \-------------------/ ---->
7097 for (x = 0; x < (a->used - b); x++) {
7101 /* zero the top digits */
7102 for (; x < a->used; x++) {
7107 /* remove excess digits */
7116 /* End: bn_tma_mp_rshd.c */
7118 /* Start: bn_tma_mp_set.c */
7121 /* LibTomMath, multiple-precision integer library -- Tom St Denis
7123 * LibTomMath is a library that provides multiple-precision
7124 * integer arithmetic as well as number theoretic functionality.
7126 * The library was designed directly after the MPI library by
7127 * Michael Fromberger but has been written from scratch with
7128 * additional optimizations in place.
7130 * The library is free for all purposes without any express
7131 * guarantee it works.
7133 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
7136 /* set to a digit */
7137 void tma_mp_set (tma_mp_int * a, tma_mp_digit b)
7140 a->dp[0] = b & MP_MASK;
7141 a->used = (a->dp[0] != 0) ? 1 : 0;
7149 /* End: bn_tma_mp_set.c */
7151 /* Start: bn_tma_mp_set_int.c */
7153 #ifdef BN_MP_SET_INT_C
7154 /* LibTomMath, multiple-precision integer library -- Tom St Denis
7156 * LibTomMath is a library that provides multiple-precision
7157 * integer arithmetic as well as number theoretic functionality.
7159 * The library was designed directly after the MPI library by
7160 * Michael Fromberger but has been written from scratch with
7161 * additional optimizations in place.
7163 * The library is free for all purposes without any express
7164 * guarantee it works.
7166 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
7169 /* set a 32-bit const */
7170 int tma_mp_set_int (tma_mp_int * a, unsigned long b)
7176 /* set four bits at a time */
7177 for (x = 0; x < 8; x++) {
7178 /* shift the number up four bits */
7179 if ((res = tma_mp_mul_2d (a, 4, a)) != MP_OKAY) {
7183 /* OR in the top four bits of the source */
7184 a->dp[0] |= (b >> 28) & 15;
7186 /* shift the source up to the next four bits */
7189 /* ensure that digits are not clamped off */
7201 /* End: bn_tma_mp_set_int.c */
7203 /* Start: bn_tma_mp_shrink.c */
7205 #ifdef BN_MP_SHRINK_C
7206 /* LibTomMath, multiple-precision integer library -- Tom St Denis
7208 * LibTomMath is a library that provides multiple-precision
7209 * integer arithmetic as well as number theoretic functionality.
7211 * The library was designed directly after the MPI library by
7212 * Michael Fromberger but has been written from scratch with
7213 * additional optimizations in place.
7215 * The library is free for all purposes without any express
7216 * guarantee it works.
7218 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
7221 /* shrink a bignum */
7222 int tma_mp_shrink (tma_mp_int * a)
7225 if (a->alloc != a->used && a->used > 0) {
7226 if ((tmp = OPT_CAST(tma_mp_digit) XREALLOC (a->dp, sizeof (tma_mp_digit) * a->used)) == NULL) {
7240 /* End: bn_tma_mp_shrink.c */
7242 /* Start: bn_tma_mp_signed_bin_size.c */
7244 #ifdef BN_MP_SIGNED_BIN_SIZE_C
7245 /* LibTomMath, multiple-precision integer library -- Tom St Denis
7247 * LibTomMath is a library that provides multiple-precision
7248 * integer arithmetic as well as number theoretic functionality.
7250 * The library was designed directly after the MPI library by
7251 * Michael Fromberger but has been written from scratch with
7252 * additional optimizations in place.
7254 * The library is free for all purposes without any express
7255 * guarantee it works.
7257 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
7260 /* get the size for an signed equivalent */
7261 int tma_mp_signed_bin_size (tma_mp_int * a)
7263 return 1 + tma_mp_unsigned_bin_size (a);
7271 /* End: bn_tma_mp_signed_bin_size.c */
7273 /* Start: bn_tma_mp_sqr.c */
7276 /* LibTomMath, multiple-precision integer library -- Tom St Denis
7278 * LibTomMath is a library that provides multiple-precision
7279 * integer arithmetic as well as number theoretic functionality.
7281 * The library was designed directly after the MPI library by
7282 * Michael Fromberger but has been written from scratch with
7283 * additional optimizations in place.
7285 * The library is free for all purposes without any express
7286 * guarantee it works.
7288 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
7291 /* computes b = a*a */
7293 tma_mp_sqr (tma_mp_int * a, tma_mp_int * b)
7297 #ifdef BN_MP_TOOM_SQR_C
7298 /* use Toom-Cook? */
7299 if (a->used >= TOOM_SQR_CUTOFF) {
7300 res = tma_mp_toom_sqr(a, b);
7304 #ifdef BN_MP_KARATSUBA_SQR_C
7305 if (a->used >= KARATSUBA_SQR_CUTOFF) {
7306 res = tma_mp_karatsuba_sqr (a, b);
7310 #ifdef BN_FAST_S_MP_SQR_C
7311 /* can we use the fast comba multiplier? */
7312 if ((a->used * 2 + 1) < MP_WARRAY &&
7314 (1 << (sizeof(tma_mp_word) * CHAR_BIT - 2*DIGIT_BIT - 1))) {
7315 res = fast_s_tma_mp_sqr (a, b);
7318 #ifdef BN_S_MP_SQR_C
7319 res = s_tma_mp_sqr (a, b);
7333 /* End: bn_tma_mp_sqr.c */
7335 /* Start: bn_tma_mp_sqrmod.c */
7337 #ifdef BN_MP_SQRMOD_C
7338 /* LibTomMath, multiple-precision integer library -- Tom St Denis
7340 * LibTomMath is a library that provides multiple-precision
7341 * integer arithmetic as well as number theoretic functionality.
7343 * The library was designed directly after the MPI library by
7344 * Michael Fromberger but has been written from scratch with
7345 * additional optimizations in place.
7347 * The library is free for all purposes without any express
7348 * guarantee it works.
7350 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
7353 /* c = a * a (mod b) */
7355 tma_mp_sqrmod (tma_mp_int * a, tma_mp_int * b, tma_mp_int * c)
7360 if ((res = tma_mp_init (&t)) != MP_OKAY) {
7364 if ((res = tma_mp_sqr (a, &t)) != MP_OKAY) {
7368 res = tma_mp_mod (&t, b, c);
7378 /* End: bn_tma_mp_sqrmod.c */
7380 /* Start: bn_tma_mp_sqrt.c */
7383 /* LibTomMath, multiple-precision integer library -- Tom St Denis
7385 * LibTomMath is a library that provides multiple-precision
7386 * integer arithmetic as well as number theoretic functionality.
7388 * The library was designed directly after the MPI library by
7389 * Michael Fromberger but has been written from scratch with
7390 * additional optimizations in place.
7392 * The library is free for all purposes without any express
7393 * guarantee it works.
7395 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
7398 /* this function is less generic than tma_mp_n_root, simpler and faster */
7399 int tma_mp_sqrt(tma_mp_int *arg, tma_mp_int *ret)
7404 /* must be positive */
7405 if (arg->sign == MP_NEG) {
7410 if (tma_mp_iszero(arg) == MP_YES) {
7415 if ((res = tma_mp_init_copy(&t1, arg)) != MP_OKAY) {
7419 if ((res = tma_mp_init(&t2)) != MP_OKAY) {
7423 /* First approx. (not very bad for large arg) */
7424 tma_mp_rshd (&t1,t1.used/2);
7427 if ((res = tma_mp_div(arg,&t1,&t2,NULL)) != MP_OKAY) {
7430 if ((res = tma_mp_add(&t1,&t2,&t1)) != MP_OKAY) {
7433 if ((res = tma_mp_div_2(&t1,&t1)) != MP_OKAY) {
7436 /* And now t1 > sqrt(arg) */
7438 if ((res = tma_mp_div(arg,&t1,&t2,NULL)) != MP_OKAY) {
7441 if ((res = tma_mp_add(&t1,&t2,&t1)) != MP_OKAY) {
7444 if ((res = tma_mp_div_2(&t1,&t1)) != MP_OKAY) {
7447 /* t1 >= sqrt(arg) >= t2 at this point */
7448 } while (tma_mp_cmp_mag(&t1,&t2) == MP_GT);
7450 tma_mp_exch(&t1,ret);
7452 E1: tma_mp_clear(&t2);
7453 E2: tma_mp_clear(&t1);
7463 /* End: bn_tma_mp_sqrt.c */
7465 /* Start: bn_tma_mp_sub.c */
7468 /* LibTomMath, multiple-precision integer library -- Tom St Denis
7470 * LibTomMath is a library that provides multiple-precision
7471 * integer arithmetic as well as number theoretic functionality.
7473 * The library was designed directly after the MPI library by
7474 * Michael Fromberger but has been written from scratch with
7475 * additional optimizations in place.
7477 * The library is free for all purposes without any express
7478 * guarantee it works.
7480 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
7483 /* high level subtraction (handles signs) */
7485 tma_mp_sub (tma_mp_int * a, tma_mp_int * b, tma_mp_int * c)
7493 /* subtract a negative from a positive, OR */
7494 /* subtract a positive from a negative. */
7495 /* In either case, ADD their magnitudes, */
7496 /* and use the sign of the first number. */
7498 res = s_tma_mp_add (a, b, c);
7500 /* subtract a positive from a positive, OR */
7501 /* subtract a negative from a negative. */
7502 /* First, take the difference between their */
7503 /* magnitudes, then... */
7504 if (tma_mp_cmp_mag (a, b) != MP_LT) {
7505 /* Copy the sign from the first */
7507 /* The first has a larger or equal magnitude */
7508 res = s_tma_mp_sub (a, b, c);
7510 /* The result has the *opposite* sign from */
7511 /* the first number. */
7512 c->sign = (sa == MP_ZPOS) ? MP_NEG : MP_ZPOS;
7513 /* The second has a larger magnitude */
7514 res = s_tma_mp_sub (b, a, c);
7526 /* End: bn_tma_mp_sub.c */
7528 /* Start: bn_tma_mp_sub_d.c */
7530 #ifdef BN_MP_SUB_D_C
7531 /* LibTomMath, multiple-precision integer library -- Tom St Denis
7533 * LibTomMath is a library that provides multiple-precision
7534 * integer arithmetic as well as number theoretic functionality.
7536 * The library was designed directly after the MPI library by
7537 * Michael Fromberger but has been written from scratch with
7538 * additional optimizations in place.
7540 * The library is free for all purposes without any express
7541 * guarantee it works.
7543 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
7546 /* single digit subtraction */
7548 tma_mp_sub_d (tma_mp_int * a, tma_mp_digit b, tma_mp_int * c)
7550 tma_mp_digit *tmpa, *tmpc, mu;
7551 int res, ix, oldused;
7553 /* grow c as required */
7554 if (c->alloc < a->used + 1) {
7555 if ((res = tma_mp_grow(c, a->used + 1)) != MP_OKAY) {
7560 /* if a is negative just do an unsigned
7561 * addition [with fudged signs]
7563 if (a->sign == MP_NEG) {
7565 res = tma_mp_add_d(a, b, c);
7566 a->sign = c->sign = MP_NEG;
7579 /* if a <= b simply fix the single digit */
7580 if ((a->used == 1 && a->dp[0] <= b) || a->used == 0) {
7582 *tmpc++ = b - *tmpa;
7588 /* negative/1digit */
7596 /* subtract first digit */
7597 *tmpc = *tmpa++ - b;
7598 mu = *tmpc >> (sizeof(tma_mp_digit) * CHAR_BIT - 1);
7601 /* handle rest of the digits */
7602 for (ix = 1; ix < a->used; ix++) {
7603 *tmpc = *tmpa++ - mu;
7604 mu = *tmpc >> (sizeof(tma_mp_digit) * CHAR_BIT - 1);
7609 /* zero excess digits */
7610 while (ix++ < oldused) {
7623 /* End: bn_tma_mp_sub_d.c */
7625 /* Start: bn_tma_mp_submod.c */
7627 #ifdef BN_MP_SUBMOD_C
7628 /* LibTomMath, multiple-precision integer library -- Tom St Denis
7630 * LibTomMath is a library that provides multiple-precision
7631 * integer arithmetic as well as number theoretic functionality.
7633 * The library was designed directly after the MPI library by
7634 * Michael Fromberger but has been written from scratch with
7635 * additional optimizations in place.
7637 * The library is free for all purposes without any express
7638 * guarantee it works.
7640 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
7643 /* d = a - b (mod c) */
7645 tma_mp_submod (tma_mp_int * a, tma_mp_int * b, tma_mp_int * c, tma_mp_int * d)
7651 if ((res = tma_mp_init (&t)) != MP_OKAY) {
7655 if ((res = tma_mp_sub (a, b, &t)) != MP_OKAY) {
7659 res = tma_mp_mod (&t, c, d);
7669 /* End: bn_tma_mp_submod.c */
7671 /* Start: bn_tma_mp_to_signed_bin.c */
7673 #ifdef BN_MP_TO_SIGNED_BIN_C
7674 /* LibTomMath, multiple-precision integer library -- Tom St Denis
7676 * LibTomMath is a library that provides multiple-precision
7677 * integer arithmetic as well as number theoretic functionality.
7679 * The library was designed directly after the MPI library by
7680 * Michael Fromberger but has been written from scratch with
7681 * additional optimizations in place.
7683 * The library is free for all purposes without any express
7684 * guarantee it works.
7686 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
7689 /* store in signed [big endian] format */
7690 int tma_mp_to_signed_bin (tma_mp_int * a, unsigned char *b)
7694 if ((res = tma_mp_to_unsigned_bin (a, b + 1)) != MP_OKAY) {
7697 b[0] = (unsigned char) ((a->sign == MP_ZPOS) ? 0 : 1);
7706 /* End: bn_tma_mp_to_signed_bin.c */
7708 /* Start: bn_tma_mp_to_signed_bin_n.c */
7710 #ifdef BN_MP_TO_SIGNED_BIN_N_C
7711 /* LibTomMath, multiple-precision integer library -- Tom St Denis
7713 * LibTomMath is a library that provides multiple-precision
7714 * integer arithmetic as well as number theoretic functionality.
7716 * The library was designed directly after the MPI library by
7717 * Michael Fromberger but has been written from scratch with
7718 * additional optimizations in place.
7720 * The library is free for all purposes without any express
7721 * guarantee it works.
7723 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
7726 /* store in signed [big endian] format */
7727 int tma_mp_to_signed_bin_n (tma_mp_int * a, unsigned char *b, unsigned long *outlen)
7729 if (*outlen < (unsigned long)tma_mp_signed_bin_size(a)) {
7732 *outlen = tma_mp_signed_bin_size(a);
7733 return tma_mp_to_signed_bin(a, b);
7741 /* End: bn_tma_mp_to_signed_bin_n.c */
7743 /* Start: bn_tma_mp_to_unsigned_bin.c */
7745 #ifdef BN_MP_TO_UNSIGNED_BIN_C
7746 /* LibTomMath, multiple-precision integer library -- Tom St Denis
7748 * LibTomMath is a library that provides multiple-precision
7749 * integer arithmetic as well as number theoretic functionality.
7751 * The library was designed directly after the MPI library by
7752 * Michael Fromberger but has been written from scratch with
7753 * additional optimizations in place.
7755 * The library is free for all purposes without any express
7756 * guarantee it works.
7758 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
7761 /* store in unsigned [big endian] format */
7762 int tma_mp_to_unsigned_bin (tma_mp_int * a, unsigned char *b)
7767 if ((res = tma_mp_init_copy (&t, a)) != MP_OKAY) {
7772 while (tma_mp_iszero (&t) == 0) {
7774 b[x++] = (unsigned char) (t.dp[0] & 255);
7776 b[x++] = (unsigned char) (t.dp[0] | ((t.dp[1] & 0x01) << 7));
7778 if ((res = tma_mp_div_2d (&t, 8, &t, NULL)) != MP_OKAY) {
7793 /* End: bn_tma_mp_to_unsigned_bin.c */
7795 /* Start: bn_tma_mp_to_unsigned_bin_n.c */
7797 #ifdef BN_MP_TO_UNSIGNED_BIN_N_C
7798 /* LibTomMath, multiple-precision integer library -- Tom St Denis
7800 * LibTomMath is a library that provides multiple-precision
7801 * integer arithmetic as well as number theoretic functionality.
7803 * The library was designed directly after the MPI library by
7804 * Michael Fromberger but has been written from scratch with
7805 * additional optimizations in place.
7807 * The library is free for all purposes without any express
7808 * guarantee it works.
7810 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
7813 /* store in unsigned [big endian] format */
7814 int tma_mp_to_unsigned_bin_n (tma_mp_int * a, unsigned char *b, unsigned long *outlen)
7816 if (*outlen < (unsigned long)tma_mp_unsigned_bin_size(a)) {
7819 *outlen = tma_mp_unsigned_bin_size(a);
7820 return tma_mp_to_unsigned_bin(a, b);
7828 /* End: bn_tma_mp_to_unsigned_bin_n.c */
7830 /* Start: bn_tma_mp_toom_mul.c */
7832 #ifdef BN_MP_TOOM_MUL_C
7833 /* LibTomMath, multiple-precision integer library -- Tom St Denis
7835 * LibTomMath is a library that provides multiple-precision
7836 * integer arithmetic as well as number theoretic functionality.
7838 * The library was designed directly after the MPI library by
7839 * Michael Fromberger but has been written from scratch with
7840 * additional optimizations in place.
7842 * The library is free for all purposes without any express
7843 * guarantee it works.
7845 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
7848 /* multiplication using the Toom-Cook 3-way algorithm
7850 * Much more complicated than Karatsuba but has a lower
7851 * asymptotic running time of O(N**1.464). This algorithm is
7852 * only particularly useful on VERY large inputs
7853 * (we're talking 1000s of digits here...).
7855 int tma_mp_toom_mul(tma_mp_int *a, tma_mp_int *b, tma_mp_int *c)
7857 tma_mp_int w0, w1, w2, w3, w4, tmp1, tmp2, a0, a1, a2, b0, b1, b2;
7861 if ((res = tma_mp_init_multi(&w0, &w1, &w2, &w3, &w4,
7862 &a0, &a1, &a2, &b0, &b1,
7863 &b2, &tmp1, &tmp2, NULL)) != MP_OKAY) {
7868 B = MIN(a->used, b->used) / 3;
7870 /* a = a2 * B**2 + a1 * B + a0 */
7871 if ((res = tma_mp_mod_2d(a, DIGIT_BIT * B, &a0)) != MP_OKAY) {
7875 if ((res = tma_mp_copy(a, &a1)) != MP_OKAY) {
7878 tma_mp_rshd(&a1, B);
7879 tma_mp_mod_2d(&a1, DIGIT_BIT * B, &a1);
7881 if ((res = tma_mp_copy(a, &a2)) != MP_OKAY) {
7884 tma_mp_rshd(&a2, B*2);
7886 /* b = b2 * B**2 + b1 * B + b0 */
7887 if ((res = tma_mp_mod_2d(b, DIGIT_BIT * B, &b0)) != MP_OKAY) {
7891 if ((res = tma_mp_copy(b, &b1)) != MP_OKAY) {
7894 tma_mp_rshd(&b1, B);
7895 tma_mp_mod_2d(&b1, DIGIT_BIT * B, &b1);
7897 if ((res = tma_mp_copy(b, &b2)) != MP_OKAY) {
7900 tma_mp_rshd(&b2, B*2);
7903 if ((res = tma_mp_mul(&a0, &b0, &w0)) != MP_OKAY) {
7908 if ((res = tma_mp_mul(&a2, &b2, &w4)) != MP_OKAY) {
7912 /* w1 = (a2 + 2(a1 + 2a0))(b2 + 2(b1 + 2b0)) */
7913 if ((res = tma_mp_mul_2(&a0, &tmp1)) != MP_OKAY) {
7916 if ((res = tma_mp_add(&tmp1, &a1, &tmp1)) != MP_OKAY) {
7919 if ((res = tma_mp_mul_2(&tmp1, &tmp1)) != MP_OKAY) {
7922 if ((res = tma_mp_add(&tmp1, &a2, &tmp1)) != MP_OKAY) {
7926 if ((res = tma_mp_mul_2(&b0, &tmp2)) != MP_OKAY) {
7929 if ((res = tma_mp_add(&tmp2, &b1, &tmp2)) != MP_OKAY) {
7932 if ((res = tma_mp_mul_2(&tmp2, &tmp2)) != MP_OKAY) {
7935 if ((res = tma_mp_add(&tmp2, &b2, &tmp2)) != MP_OKAY) {
7939 if ((res = tma_mp_mul(&tmp1, &tmp2, &w1)) != MP_OKAY) {
7943 /* w3 = (a0 + 2(a1 + 2a2))(b0 + 2(b1 + 2b2)) */
7944 if ((res = tma_mp_mul_2(&a2, &tmp1)) != MP_OKAY) {
7947 if ((res = tma_mp_add(&tmp1, &a1, &tmp1)) != MP_OKAY) {
7950 if ((res = tma_mp_mul_2(&tmp1, &tmp1)) != MP_OKAY) {
7953 if ((res = tma_mp_add(&tmp1, &a0, &tmp1)) != MP_OKAY) {
7957 if ((res = tma_mp_mul_2(&b2, &tmp2)) != MP_OKAY) {
7960 if ((res = tma_mp_add(&tmp2, &b1, &tmp2)) != MP_OKAY) {
7963 if ((res = tma_mp_mul_2(&tmp2, &tmp2)) != MP_OKAY) {
7966 if ((res = tma_mp_add(&tmp2, &b0, &tmp2)) != MP_OKAY) {
7970 if ((res = tma_mp_mul(&tmp1, &tmp2, &w3)) != MP_OKAY) {
7975 /* w2 = (a2 + a1 + a0)(b2 + b1 + b0) */
7976 if ((res = tma_mp_add(&a2, &a1, &tmp1)) != MP_OKAY) {
7979 if ((res = tma_mp_add(&tmp1, &a0, &tmp1)) != MP_OKAY) {
7982 if ((res = tma_mp_add(&b2, &b1, &tmp2)) != MP_OKAY) {
7985 if ((res = tma_mp_add(&tmp2, &b0, &tmp2)) != MP_OKAY) {
7988 if ((res = tma_mp_mul(&tmp1, &tmp2, &w2)) != MP_OKAY) {
7992 /* now solve the matrix
8000 using 12 subtractions, 4 shifts,
8001 2 small divisions and 1 small multiplication
8005 if ((res = tma_mp_sub(&w1, &w4, &w1)) != MP_OKAY) {
8009 if ((res = tma_mp_sub(&w3, &w0, &w3)) != MP_OKAY) {
8013 if ((res = tma_mp_div_2(&w1, &w1)) != MP_OKAY) {
8017 if ((res = tma_mp_div_2(&w3, &w3)) != MP_OKAY) {
8021 if ((res = tma_mp_sub(&w2, &w0, &w2)) != MP_OKAY) {
8024 if ((res = tma_mp_sub(&w2, &w4, &w2)) != MP_OKAY) {
8028 if ((res = tma_mp_sub(&w1, &w2, &w1)) != MP_OKAY) {
8032 if ((res = tma_mp_sub(&w3, &w2, &w3)) != MP_OKAY) {
8036 if ((res = tma_mp_mul_2d(&w0, 3, &tmp1)) != MP_OKAY) {
8039 if ((res = tma_mp_sub(&w1, &tmp1, &w1)) != MP_OKAY) {
8043 if ((res = tma_mp_mul_2d(&w4, 3, &tmp1)) != MP_OKAY) {
8046 if ((res = tma_mp_sub(&w3, &tmp1, &w3)) != MP_OKAY) {
8050 if ((res = tma_mp_mul_d(&w2, 3, &w2)) != MP_OKAY) {
8053 if ((res = tma_mp_sub(&w2, &w1, &w2)) != MP_OKAY) {
8056 if ((res = tma_mp_sub(&w2, &w3, &w2)) != MP_OKAY) {
8060 if ((res = tma_mp_sub(&w1, &w2, &w1)) != MP_OKAY) {
8064 if ((res = tma_mp_sub(&w3, &w2, &w3)) != MP_OKAY) {
8068 if ((res = tma_mp_div_3(&w1, &w1, NULL)) != MP_OKAY) {
8072 if ((res = tma_mp_div_3(&w3, &w3, NULL)) != MP_OKAY) {
8076 /* at this point shift W[n] by B*n */
8077 if ((res = tma_mp_lshd(&w1, 1*B)) != MP_OKAY) {
8080 if ((res = tma_mp_lshd(&w2, 2*B)) != MP_OKAY) {
8083 if ((res = tma_mp_lshd(&w3, 3*B)) != MP_OKAY) {
8086 if ((res = tma_mp_lshd(&w4, 4*B)) != MP_OKAY) {
8090 if ((res = tma_mp_add(&w0, &w1, c)) != MP_OKAY) {
8093 if ((res = tma_mp_add(&w2, &w3, &tmp1)) != MP_OKAY) {
8096 if ((res = tma_mp_add(&w4, &tmp1, &tmp1)) != MP_OKAY) {
8099 if ((res = tma_mp_add(&tmp1, c, c)) != MP_OKAY) {
8104 tma_mp_clear_multi(&w0, &w1, &w2, &w3, &w4,
8105 &a0, &a1, &a2, &b0, &b1,
8106 &b2, &tmp1, &tmp2, NULL);
8116 /* End: bn_tma_mp_toom_mul.c */
8118 /* Start: bn_tma_mp_toom_sqr.c */
8120 #ifdef BN_MP_TOOM_SQR_C
8121 /* LibTomMath, multiple-precision integer library -- Tom St Denis
8123 * LibTomMath is a library that provides multiple-precision
8124 * integer arithmetic as well as number theoretic functionality.
8126 * The library was designed directly after the MPI library by
8127 * Michael Fromberger but has been written from scratch with
8128 * additional optimizations in place.
8130 * The library is free for all purposes without any express
8131 * guarantee it works.
8133 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
8136 /* squaring using Toom-Cook 3-way algorithm */
8138 tma_mp_toom_sqr(tma_mp_int *a, tma_mp_int *b)
8140 tma_mp_int w0, w1, w2, w3, w4, tmp1, a0, a1, a2;
8144 if ((res = tma_mp_init_multi(&w0, &w1, &w2, &w3, &w4, &a0, &a1, &a2, &tmp1, NULL)) != MP_OKAY) {
8151 /* a = a2 * B**2 + a1 * B + a0 */
8152 if ((res = tma_mp_mod_2d(a, DIGIT_BIT * B, &a0)) != MP_OKAY) {
8156 if ((res = tma_mp_copy(a, &a1)) != MP_OKAY) {
8159 tma_mp_rshd(&a1, B);
8160 tma_mp_mod_2d(&a1, DIGIT_BIT * B, &a1);
8162 if ((res = tma_mp_copy(a, &a2)) != MP_OKAY) {
8165 tma_mp_rshd(&a2, B*2);
8168 if ((res = tma_mp_sqr(&a0, &w0)) != MP_OKAY) {
8173 if ((res = tma_mp_sqr(&a2, &w4)) != MP_OKAY) {
8177 /* w1 = (a2 + 2(a1 + 2a0))**2 */
8178 if ((res = tma_mp_mul_2(&a0, &tmp1)) != MP_OKAY) {
8181 if ((res = tma_mp_add(&tmp1, &a1, &tmp1)) != MP_OKAY) {
8184 if ((res = tma_mp_mul_2(&tmp1, &tmp1)) != MP_OKAY) {
8187 if ((res = tma_mp_add(&tmp1, &a2, &tmp1)) != MP_OKAY) {
8191 if ((res = tma_mp_sqr(&tmp1, &w1)) != MP_OKAY) {
8195 /* w3 = (a0 + 2(a1 + 2a2))**2 */
8196 if ((res = tma_mp_mul_2(&a2, &tmp1)) != MP_OKAY) {
8199 if ((res = tma_mp_add(&tmp1, &a1, &tmp1)) != MP_OKAY) {
8202 if ((res = tma_mp_mul_2(&tmp1, &tmp1)) != MP_OKAY) {
8205 if ((res = tma_mp_add(&tmp1, &a0, &tmp1)) != MP_OKAY) {
8209 if ((res = tma_mp_sqr(&tmp1, &w3)) != MP_OKAY) {
8214 /* w2 = (a2 + a1 + a0)**2 */
8215 if ((res = tma_mp_add(&a2, &a1, &tmp1)) != MP_OKAY) {
8218 if ((res = tma_mp_add(&tmp1, &a0, &tmp1)) != MP_OKAY) {
8221 if ((res = tma_mp_sqr(&tmp1, &w2)) != MP_OKAY) {
8225 /* now solve the matrix
8233 using 12 subtractions, 4 shifts, 2 small divisions and 1 small multiplication.
8237 if ((res = tma_mp_sub(&w1, &w4, &w1)) != MP_OKAY) {
8241 if ((res = tma_mp_sub(&w3, &w0, &w3)) != MP_OKAY) {
8245 if ((res = tma_mp_div_2(&w1, &w1)) != MP_OKAY) {
8249 if ((res = tma_mp_div_2(&w3, &w3)) != MP_OKAY) {
8253 if ((res = tma_mp_sub(&w2, &w0, &w2)) != MP_OKAY) {
8256 if ((res = tma_mp_sub(&w2, &w4, &w2)) != MP_OKAY) {
8260 if ((res = tma_mp_sub(&w1, &w2, &w1)) != MP_OKAY) {
8264 if ((res = tma_mp_sub(&w3, &w2, &w3)) != MP_OKAY) {
8268 if ((res = tma_mp_mul_2d(&w0, 3, &tmp1)) != MP_OKAY) {
8271 if ((res = tma_mp_sub(&w1, &tmp1, &w1)) != MP_OKAY) {
8275 if ((res = tma_mp_mul_2d(&w4, 3, &tmp1)) != MP_OKAY) {
8278 if ((res = tma_mp_sub(&w3, &tmp1, &w3)) != MP_OKAY) {
8282 if ((res = tma_mp_mul_d(&w2, 3, &w2)) != MP_OKAY) {
8285 if ((res = tma_mp_sub(&w2, &w1, &w2)) != MP_OKAY) {
8288 if ((res = tma_mp_sub(&w2, &w3, &w2)) != MP_OKAY) {
8292 if ((res = tma_mp_sub(&w1, &w2, &w1)) != MP_OKAY) {
8296 if ((res = tma_mp_sub(&w3, &w2, &w3)) != MP_OKAY) {
8300 if ((res = tma_mp_div_3(&w1, &w1, NULL)) != MP_OKAY) {
8304 if ((res = tma_mp_div_3(&w3, &w3, NULL)) != MP_OKAY) {
8308 /* at this point shift W[n] by B*n */
8309 if ((res = tma_mp_lshd(&w1, 1*B)) != MP_OKAY) {
8312 if ((res = tma_mp_lshd(&w2, 2*B)) != MP_OKAY) {
8315 if ((res = tma_mp_lshd(&w3, 3*B)) != MP_OKAY) {
8318 if ((res = tma_mp_lshd(&w4, 4*B)) != MP_OKAY) {
8322 if ((res = tma_mp_add(&w0, &w1, b)) != MP_OKAY) {
8325 if ((res = tma_mp_add(&w2, &w3, &tmp1)) != MP_OKAY) {
8328 if ((res = tma_mp_add(&w4, &tmp1, &tmp1)) != MP_OKAY) {
8331 if ((res = tma_mp_add(&tmp1, b, b)) != MP_OKAY) {
8336 tma_mp_clear_multi(&w0, &w1, &w2, &w3, &w4, &a0, &a1, &a2, &tmp1, NULL);
8346 /* End: bn_tma_mp_toom_sqr.c */
8348 /* Start: bn_tma_mp_toradix.c */
8350 #ifdef BN_MP_TORADIX_C
8351 /* LibTomMath, multiple-precision integer library -- Tom St Denis
8353 * LibTomMath is a library that provides multiple-precision
8354 * integer arithmetic as well as number theoretic functionality.
8356 * The library was designed directly after the MPI library by
8357 * Michael Fromberger but has been written from scratch with
8358 * additional optimizations in place.
8360 * The library is free for all purposes without any express
8361 * guarantee it works.
8363 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
8366 /* stores a bignum as a ASCII string in a given radix (2..64) */
8367 int tma_mp_toradix (tma_mp_int * a, char *str, int radix)
8374 /* check range of the radix */
8375 if (radix < 2 || radix > 64) {
8379 /* quick out if its zero */
8380 if (tma_mp_iszero(a) == 1) {
8386 if ((res = tma_mp_init_copy (&t, a)) != MP_OKAY) {
8390 /* if it is negative output a - */
8391 if (t.sign == MP_NEG) {
8398 while (tma_mp_iszero (&t) == 0) {
8399 if ((res = tma_mp_div_d (&t, (tma_mp_digit) radix, &t, &d)) != MP_OKAY) {
8403 *str++ = tma_mp_s_rmap[d];
8407 /* reverse the digits of the string. In this case _s points
8408 * to the first digit [exluding the sign] of the number]
8410 bn_reverse ((unsigned char *)_s, digs);
8412 /* append a NULL so the string is properly terminated */
8425 /* End: bn_tma_mp_toradix.c */
8427 /* Start: bn_tma_mp_toradix_n.c */
8429 #ifdef BN_MP_TORADIX_N_C
8430 /* LibTomMath, multiple-precision integer library -- Tom St Denis
8432 * LibTomMath is a library that provides multiple-precision
8433 * integer arithmetic as well as number theoretic functionality.
8435 * The library was designed directly after the MPI library by
8436 * Michael Fromberger but has been written from scratch with
8437 * additional optimizations in place.
8439 * The library is free for all purposes without any express
8440 * guarantee it works.
8442 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
8445 /* stores a bignum as a ASCII string in a given radix (2..64)
8447 * Stores upto maxlen-1 chars and always a NULL byte
8449 int tma_mp_toradix_n(tma_mp_int * a, char *str, int radix, int maxlen)
8456 /* check range of the maxlen, radix */
8457 if (maxlen < 2 || radix < 2 || radix > 64) {
8461 /* quick out if its zero */
8462 if (tma_mp_iszero(a) == MP_YES) {
8468 if ((res = tma_mp_init_copy (&t, a)) != MP_OKAY) {
8472 /* if it is negative output a - */
8473 if (t.sign == MP_NEG) {
8474 /* we have to reverse our digits later... but not the - sign!! */
8477 /* store the flag and mark the number as positive */
8481 /* subtract a char */
8486 while (tma_mp_iszero (&t) == 0) {
8491 if ((res = tma_mp_div_d (&t, (tma_mp_digit) radix, &t, &d)) != MP_OKAY) {
8495 *str++ = tma_mp_s_rmap[d];
8499 /* reverse the digits of the string. In this case _s points
8500 * to the first digit [exluding the sign] of the number
8502 bn_reverse ((unsigned char *)_s, digs);
8504 /* append a NULL so the string is properly terminated */
8517 /* End: bn_tma_mp_toradix_n.c */
8519 /* Start: bn_tma_mp_unsigned_bin_size.c */
8521 #ifdef BN_MP_UNSIGNED_BIN_SIZE_C
8522 /* LibTomMath, multiple-precision integer library -- Tom St Denis
8524 * LibTomMath is a library that provides multiple-precision
8525 * integer arithmetic as well as number theoretic functionality.
8527 * The library was designed directly after the MPI library by
8528 * Michael Fromberger but has been written from scratch with
8529 * additional optimizations in place.
8531 * The library is free for all purposes without any express
8532 * guarantee it works.
8534 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
8537 /* get the size for an unsigned equivalent */
8538 int tma_mp_unsigned_bin_size (tma_mp_int * a)
8540 int size = tma_mp_count_bits (a);
8541 return (size / 8 + ((size & 7) != 0 ? 1 : 0));
8549 /* End: bn_tma_mp_unsigned_bin_size.c */
8551 /* Start: bn_tma_mp_xor.c */
8554 /* LibTomMath, multiple-precision integer library -- Tom St Denis
8556 * LibTomMath is a library that provides multiple-precision
8557 * integer arithmetic as well as number theoretic functionality.
8559 * The library was designed directly after the MPI library by
8560 * Michael Fromberger but has been written from scratch with
8561 * additional optimizations in place.
8563 * The library is free for all purposes without any express
8564 * guarantee it works.
8566 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
8569 /* XOR two ints together */
8571 tma_mp_xor (tma_mp_int * a, tma_mp_int * b, tma_mp_int * c)
8576 if (a->used > b->used) {
8577 if ((res = tma_mp_init_copy (&t, a)) != MP_OKAY) {
8583 if ((res = tma_mp_init_copy (&t, b)) != MP_OKAY) {
8590 for (ix = 0; ix < px; ix++) {
8591 t.dp[ix] ^= x->dp[ix];
8594 tma_mp_exch (c, &t);
8604 /* End: bn_tma_mp_xor.c */
8606 /* Start: bn_tma_mp_zero.c */
8609 /* LibTomMath, multiple-precision integer library -- Tom St Denis
8611 * LibTomMath is a library that provides multiple-precision
8612 * integer arithmetic as well as number theoretic functionality.
8614 * The library was designed directly after the MPI library by
8615 * Michael Fromberger but has been written from scratch with
8616 * additional optimizations in place.
8618 * The library is free for all purposes without any express
8619 * guarantee it works.
8621 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
8625 void tma_mp_zero (tma_mp_int * a)
8634 for (n = 0; n < a->alloc; n++) {
8644 /* End: bn_tma_mp_zero.c */
8646 /* Start: bn_prime_tab.c */
8648 #ifdef BN_PRIME_TAB_C
8649 /* LibTomMath, multiple-precision integer library -- Tom St Denis
8651 * LibTomMath is a library that provides multiple-precision
8652 * integer arithmetic as well as number theoretic functionality.
8654 * The library was designed directly after the MPI library by
8655 * Michael Fromberger but has been written from scratch with
8656 * additional optimizations in place.
8658 * The library is free for all purposes without any express
8659 * guarantee it works.
8661 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
8663 const tma_mp_digit ltm_prime_tab[] = {
8664 0x0002, 0x0003, 0x0005, 0x0007, 0x000B, 0x000D, 0x0011, 0x0013,
8665 0x0017, 0x001D, 0x001F, 0x0025, 0x0029, 0x002B, 0x002F, 0x0035,
8666 0x003B, 0x003D, 0x0043, 0x0047, 0x0049, 0x004F, 0x0053, 0x0059,
8667 0x0061, 0x0065, 0x0067, 0x006B, 0x006D, 0x0071, 0x007F,
8670 0x0089, 0x008B, 0x0095, 0x0097, 0x009D, 0x00A3, 0x00A7, 0x00AD,
8671 0x00B3, 0x00B5, 0x00BF, 0x00C1, 0x00C5, 0x00C7, 0x00D3, 0x00DF,
8672 0x00E3, 0x00E5, 0x00E9, 0x00EF, 0x00F1, 0x00FB, 0x0101, 0x0107,
8673 0x010D, 0x010F, 0x0115, 0x0119, 0x011B, 0x0125, 0x0133, 0x0137,
8675 0x0139, 0x013D, 0x014B, 0x0151, 0x015B, 0x015D, 0x0161, 0x0167,
8676 0x016F, 0x0175, 0x017B, 0x017F, 0x0185, 0x018D, 0x0191, 0x0199,
8677 0x01A3, 0x01A5, 0x01AF, 0x01B1, 0x01B7, 0x01BB, 0x01C1, 0x01C9,
8678 0x01CD, 0x01CF, 0x01D3, 0x01DF, 0x01E7, 0x01EB, 0x01F3, 0x01F7,
8679 0x01FD, 0x0209, 0x020B, 0x021D, 0x0223, 0x022D, 0x0233, 0x0239,
8680 0x023B, 0x0241, 0x024B, 0x0251, 0x0257, 0x0259, 0x025F, 0x0265,
8681 0x0269, 0x026B, 0x0277, 0x0281, 0x0283, 0x0287, 0x028D, 0x0293,
8682 0x0295, 0x02A1, 0x02A5, 0x02AB, 0x02B3, 0x02BD, 0x02C5, 0x02CF,
8684 0x02D7, 0x02DD, 0x02E3, 0x02E7, 0x02EF, 0x02F5, 0x02F9, 0x0301,
8685 0x0305, 0x0313, 0x031D, 0x0329, 0x032B, 0x0335, 0x0337, 0x033B,
8686 0x033D, 0x0347, 0x0355, 0x0359, 0x035B, 0x035F, 0x036D, 0x0371,
8687 0x0373, 0x0377, 0x038B, 0x038F, 0x0397, 0x03A1, 0x03A9, 0x03AD,
8688 0x03B3, 0x03B9, 0x03C7, 0x03CB, 0x03D1, 0x03D7, 0x03DF, 0x03E5,
8689 0x03F1, 0x03F5, 0x03FB, 0x03FD, 0x0407, 0x0409, 0x040F, 0x0419,
8690 0x041B, 0x0425, 0x0427, 0x042D, 0x043F, 0x0443, 0x0445, 0x0449,
8691 0x044F, 0x0455, 0x045D, 0x0463, 0x0469, 0x047F, 0x0481, 0x048B,
8693 0x0493, 0x049D, 0x04A3, 0x04A9, 0x04B1, 0x04BD, 0x04C1, 0x04C7,
8694 0x04CD, 0x04CF, 0x04D5, 0x04E1, 0x04EB, 0x04FD, 0x04FF, 0x0503,
8695 0x0509, 0x050B, 0x0511, 0x0515, 0x0517, 0x051B, 0x0527, 0x0529,
8696 0x052F, 0x0551, 0x0557, 0x055D, 0x0565, 0x0577, 0x0581, 0x058F,
8697 0x0593, 0x0595, 0x0599, 0x059F, 0x05A7, 0x05AB, 0x05AD, 0x05B3,
8698 0x05BF, 0x05C9, 0x05CB, 0x05CF, 0x05D1, 0x05D5, 0x05DB, 0x05E7,
8699 0x05F3, 0x05FB, 0x0607, 0x060D, 0x0611, 0x0617, 0x061F, 0x0623,
8700 0x062B, 0x062F, 0x063D, 0x0641, 0x0647, 0x0649, 0x064D, 0x0653
8709 /* End: bn_prime_tab.c */
8711 /* Start: bn_reverse.c */
8714 /* LibTomMath, multiple-precision integer library -- Tom St Denis
8716 * LibTomMath is a library that provides multiple-precision
8717 * integer arithmetic as well as number theoretic functionality.
8719 * The library was designed directly after the MPI library by
8720 * Michael Fromberger but has been written from scratch with
8721 * additional optimizations in place.
8723 * The library is free for all purposes without any express
8724 * guarantee it works.
8726 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
8729 /* reverse an array, used for radix code */
8731 bn_reverse (unsigned char *s, int len)
8752 /* End: bn_reverse.c */
8754 /* Start: bn_s_tma_mp_add.c */
8756 #ifdef BN_S_MP_ADD_C
8757 /* LibTomMath, multiple-precision integer library -- Tom St Denis
8759 * LibTomMath is a library that provides multiple-precision
8760 * integer arithmetic as well as number theoretic functionality.
8762 * The library was designed directly after the MPI library by
8763 * Michael Fromberger but has been written from scratch with
8764 * additional optimizations in place.
8766 * The library is free for all purposes without any express
8767 * guarantee it works.
8769 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
8772 /* low level addition, based on HAC pp.594, Algorithm 14.7 */
8774 s_tma_mp_add (tma_mp_int * a, tma_mp_int * b, tma_mp_int * c)
8777 int olduse, res, min, max;
8779 /* find sizes, we let |a| <= |b| which means we have to sort
8780 * them. "x" will point to the input with the most digits
8782 if (a->used > b->used) {
8793 if (c->alloc < max + 1) {
8794 if ((res = tma_mp_grow (c, max + 1)) != MP_OKAY) {
8799 /* get old used digit count and set new one */
8804 register tma_mp_digit u, *tmpa, *tmpb, *tmpc;
8807 /* alias for digit pointers */
8818 /* zero the carry */
8820 for (i = 0; i < min; i++) {
8821 /* Compute the sum at one digit, T[i] = A[i] + B[i] + U */
8822 *tmpc = *tmpa++ + *tmpb++ + u;
8824 /* U = carry bit of T[i] */
8825 u = *tmpc >> ((tma_mp_digit)DIGIT_BIT);
8827 /* take away carry bit from T[i] */
8831 /* now copy higher words if any, that is in A+B
8832 * if A or B has more digits add those in
8835 for (; i < max; i++) {
8836 /* T[i] = X[i] + U */
8837 *tmpc = x->dp[i] + u;
8839 /* U = carry bit of T[i] */
8840 u = *tmpc >> ((tma_mp_digit)DIGIT_BIT);
8842 /* take away carry bit from T[i] */
8850 /* clear digits above oldused */
8851 for (i = c->used; i < olduse; i++) {
8865 /* End: bn_s_tma_mp_add.c */
8867 /* Start: bn_s_tma_mp_exptmod.c */
8869 #ifdef BN_S_MP_EXPTMOD_C
8870 /* LibTomMath, multiple-precision integer library -- Tom St Denis
8872 * LibTomMath is a library that provides multiple-precision
8873 * integer arithmetic as well as number theoretic functionality.
8875 * The library was designed directly after the MPI library by
8876 * Michael Fromberger but has been written from scratch with
8877 * additional optimizations in place.
8879 * The library is free for all purposes without any express
8880 * guarantee it works.
8882 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
8887 #define TAB_SIZE 256
8890 int s_tma_mp_exptmod (tma_mp_int * G, tma_mp_int * X, tma_mp_int * P, tma_mp_int * Y, int redmode)
8892 tma_mp_int M[TAB_SIZE], res, mu;
8894 int err, bitbuf, bitcpy, bitcnt, mode, digidx, x, y, winsize;
8895 int (*redux)(tma_mp_int*,tma_mp_int*,tma_mp_int*);
8897 /* find window size */
8898 x = tma_mp_count_bits (X);
8901 } else if (x <= 36) {
8903 } else if (x <= 140) {
8905 } else if (x <= 450) {
8907 } else if (x <= 1303) {
8909 } else if (x <= 3529) {
8922 /* init first cell */
8923 if ((err = tma_mp_init(&M[1])) != MP_OKAY) {
8927 /* now init the second half of the array */
8928 for (x = 1<<(winsize-1); x < (1 << winsize); x++) {
8929 if ((err = tma_mp_init(&M[x])) != MP_OKAY) {
8930 for (y = 1<<(winsize-1); y < x; y++) {
8931 tma_mp_clear (&M[y]);
8933 tma_mp_clear(&M[1]);
8938 /* create mu, used for Barrett reduction */
8939 if ((err = tma_mp_init (&mu)) != MP_OKAY) {
8944 if ((err = tma_mp_reduce_setup (&mu, P)) != MP_OKAY) {
8947 redux = tma_mp_reduce;
8949 if ((err = tma_mp_reduce_2k_setup_l (P, &mu)) != MP_OKAY) {
8952 redux = tma_mp_reduce_2k_l;
8957 * The M table contains powers of the base,
8958 * e.g. M[x] = G**x mod P
8960 * The first half of the table is not
8961 * computed though accept for M[0] and M[1]
8963 if ((err = tma_mp_mod (G, P, &M[1])) != MP_OKAY) {
8967 /* compute the value at M[1<<(winsize-1)] by squaring
8968 * M[1] (winsize-1) times
8970 if ((err = tma_mp_copy (&M[1], &M[1 << (winsize - 1)])) != MP_OKAY) {
8974 for (x = 0; x < (winsize - 1); x++) {
8976 if ((err = tma_mp_sqr (&M[1 << (winsize - 1)],
8977 &M[1 << (winsize - 1)])) != MP_OKAY) {
8981 /* reduce modulo P */
8982 if ((err = redux (&M[1 << (winsize - 1)], P, &mu)) != MP_OKAY) {
8987 /* create upper table, that is M[x] = M[x-1] * M[1] (mod P)
8988 * for x = (2**(winsize - 1) + 1) to (2**winsize - 1)
8990 for (x = (1 << (winsize - 1)) + 1; x < (1 << winsize); x++) {
8991 if ((err = tma_mp_mul (&M[x - 1], &M[1], &M[x])) != MP_OKAY) {
8994 if ((err = redux (&M[x], P, &mu)) != MP_OKAY) {
9000 if ((err = tma_mp_init (&res)) != MP_OKAY) {
9003 tma_mp_set (&res, 1);
9005 /* set initial mode and bit cnt */
9009 digidx = X->used - 1;
9014 /* grab next digit as required */
9015 if (--bitcnt == 0) {
9016 /* if digidx == -1 we are out of digits */
9020 /* read next digit and reset the bitcnt */
9021 buf = X->dp[digidx--];
9022 bitcnt = (int) DIGIT_BIT;
9025 /* grab the next msb from the exponent */
9026 y = (buf >> (tma_mp_digit)(DIGIT_BIT - 1)) & 1;
9027 buf <<= (tma_mp_digit)1;
9029 /* if the bit is zero and mode == 0 then we ignore it
9030 * These represent the leading zero bits before the first 1 bit
9031 * in the exponent. Technically this opt is not required but it
9032 * does lower the # of trivial squaring/reductions used
9034 if (mode == 0 && y == 0) {
9038 /* if the bit is zero and mode == 1 then we square */
9039 if (mode == 1 && y == 0) {
9040 if ((err = tma_mp_sqr (&res, &res)) != MP_OKAY) {
9043 if ((err = redux (&res, P, &mu)) != MP_OKAY) {
9049 /* else we add it to the window */
9050 bitbuf |= (y << (winsize - ++bitcpy));
9053 if (bitcpy == winsize) {
9054 /* ok window is filled so square as required and multiply */
9056 for (x = 0; x < winsize; x++) {
9057 if ((err = tma_mp_sqr (&res, &res)) != MP_OKAY) {
9060 if ((err = redux (&res, P, &mu)) != MP_OKAY) {
9066 if ((err = tma_mp_mul (&res, &M[bitbuf], &res)) != MP_OKAY) {
9069 if ((err = redux (&res, P, &mu)) != MP_OKAY) {
9073 /* empty window and reset */
9080 /* if bits remain then square/multiply */
9081 if (mode == 2 && bitcpy > 0) {
9082 /* square then multiply if the bit is set */
9083 for (x = 0; x < bitcpy; x++) {
9084 if ((err = tma_mp_sqr (&res, &res)) != MP_OKAY) {
9087 if ((err = redux (&res, P, &mu)) != MP_OKAY) {
9092 if ((bitbuf & (1 << winsize)) != 0) {
9094 if ((err = tma_mp_mul (&res, &M[1], &res)) != MP_OKAY) {
9097 if ((err = redux (&res, P, &mu)) != MP_OKAY) {
9104 tma_mp_exch (&res, Y);
9106 LBL_RES:tma_mp_clear (&res);
9107 LBL_MU:tma_mp_clear (&mu);
9109 tma_mp_clear(&M[1]);
9110 for (x = 1<<(winsize-1); x < (1 << winsize); x++) {
9111 tma_mp_clear (&M[x]);
9121 /* End: bn_s_tma_mp_exptmod.c */
9123 /* Start: bn_s_tma_mp_mul_digs.c */
9125 #ifdef BN_S_MP_MUL_DIGS_C
9126 /* LibTomMath, multiple-precision integer library -- Tom St Denis
9128 * LibTomMath is a library that provides multiple-precision
9129 * integer arithmetic as well as number theoretic functionality.
9131 * The library was designed directly after the MPI library by
9132 * Michael Fromberger but has been written from scratch with
9133 * additional optimizations in place.
9135 * The library is free for all purposes without any express
9136 * guarantee it works.
9138 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
9141 /* multiplies |a| * |b| and only computes upto digs digits of result
9142 * HAC pp. 595, Algorithm 14.12 Modified so you can control how
9143 * many digits of output are created.
9145 int s_tma_mp_mul_digs (tma_mp_int * a, tma_mp_int * b, tma_mp_int * c, int digs)
9148 int res, pa, pb, ix, iy;
9151 tma_mp_digit tmpx, *tmpt, *tmpy;
9153 /* can we use the fast multiplier? */
9154 if (((digs) < MP_WARRAY) &&
9155 MIN (a->used, b->used) <
9156 (1 << ((CHAR_BIT * sizeof (tma_mp_word)) - (2 * DIGIT_BIT)))) {
9157 return fast_s_tma_mp_mul_digs (a, b, c, digs);
9160 if ((res = tma_mp_init_size (&t, digs)) != MP_OKAY) {
9165 /* compute the digits of the product directly */
9167 for (ix = 0; ix < pa; ix++) {
9168 /* set the carry to zero */
9171 /* limit ourselves to making digs digits of output */
9172 pb = MIN (b->used, digs - ix);
9174 /* setup some aliases */
9175 /* copy of the digit from a used within the nested loop */
9178 /* an alias for the destination shifted ix places */
9181 /* an alias for the digits of b */
9184 /* compute the columns of the output and propagate the carry */
9185 for (iy = 0; iy < pb; iy++) {
9186 /* compute the column as a tma_mp_word */
9187 r = ((tma_mp_word)*tmpt) +
9188 ((tma_mp_word)tmpx) * ((tma_mp_word)*tmpy++) +
9191 /* the new column is the lower part of the result */
9192 *tmpt++ = (tma_mp_digit) (r & ((tma_mp_word) MP_MASK));
9194 /* get the carry word from the result */
9195 u = (tma_mp_digit) (r >> ((tma_mp_word) DIGIT_BIT));
9197 /* set carry if it is placed below digs */
9198 if (ix + iy < digs) {
9204 tma_mp_exch (&t, c);
9215 /* End: bn_s_tma_mp_mul_digs.c */
9217 /* Start: bn_s_tma_mp_mul_high_digs.c */
9219 #ifdef BN_S_MP_MUL_HIGH_DIGS_C
9220 /* LibTomMath, multiple-precision integer library -- Tom St Denis
9222 * LibTomMath is a library that provides multiple-precision
9223 * integer arithmetic as well as number theoretic functionality.
9225 * The library was designed directly after the MPI library by
9226 * Michael Fromberger but has been written from scratch with
9227 * additional optimizations in place.
9229 * The library is free for all purposes without any express
9230 * guarantee it works.
9232 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
9235 /* multiplies |a| * |b| and does not compute the lower digs digits
9236 * [meant to get the higher part of the product]
9239 s_tma_mp_mul_high_digs (tma_mp_int * a, tma_mp_int * b, tma_mp_int * c, int digs)
9242 int res, pa, pb, ix, iy;
9245 tma_mp_digit tmpx, *tmpt, *tmpy;
9247 /* can we use the fast multiplier? */
9248 #ifdef BN_FAST_S_MP_MUL_HIGH_DIGS_C
9249 if (((a->used + b->used + 1) < MP_WARRAY)
9250 && MIN (a->used, b->used) < (1 << ((CHAR_BIT * sizeof (tma_mp_word)) - (2 * DIGIT_BIT)))) {
9251 return fast_s_tma_mp_mul_high_digs (a, b, c, digs);
9255 if ((res = tma_mp_init_size (&t, a->used + b->used + 1)) != MP_OKAY) {
9258 t.used = a->used + b->used + 1;
9262 for (ix = 0; ix < pa; ix++) {
9263 /* clear the carry */
9266 /* left hand side of A[ix] * B[iy] */
9269 /* alias to the address of where the digits will be stored */
9270 tmpt = &(t.dp[digs]);
9272 /* alias for where to read the right hand side from */
9273 tmpy = b->dp + (digs - ix);
9275 for (iy = digs - ix; iy < pb; iy++) {
9276 /* calculate the double precision result */
9277 r = ((tma_mp_word)*tmpt) +
9278 ((tma_mp_word)tmpx) * ((tma_mp_word)*tmpy++) +
9281 /* get the lower part */
9282 *tmpt++ = (tma_mp_digit) (r & ((tma_mp_word) MP_MASK));
9284 /* carry the carry */
9285 u = (tma_mp_digit) (r >> ((tma_mp_word) DIGIT_BIT));
9290 tma_mp_exch (&t, c);
9300 /* End: bn_s_tma_mp_mul_high_digs.c */
9302 /* Start: bn_s_tma_mp_sqr.c */
9304 #ifdef BN_S_MP_SQR_C
9305 /* LibTomMath, multiple-precision integer library -- Tom St Denis
9307 * LibTomMath is a library that provides multiple-precision
9308 * integer arithmetic as well as number theoretic functionality.
9310 * The library was designed directly after the MPI library by
9311 * Michael Fromberger but has been written from scratch with
9312 * additional optimizations in place.
9314 * The library is free for all purposes without any express
9315 * guarantee it works.
9317 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
9320 /* low level squaring, b = a*a, HAC pp.596-597, Algorithm 14.16 */
9321 int s_tma_mp_sqr (tma_mp_int * a, tma_mp_int * b)
9324 int res, ix, iy, pa;
9326 tma_mp_digit u, tmpx, *tmpt;
9329 if ((res = tma_mp_init_size (&t, 2*pa + 1)) != MP_OKAY) {
9333 /* default used is maximum possible size */
9336 for (ix = 0; ix < pa; ix++) {
9337 /* first calculate the digit at 2*ix */
9338 /* calculate double precision result */
9339 r = ((tma_mp_word) t.dp[2*ix]) +
9340 ((tma_mp_word)a->dp[ix])*((tma_mp_word)a->dp[ix]);
9342 /* store lower part in result */
9343 t.dp[ix+ix] = (tma_mp_digit) (r & ((tma_mp_word) MP_MASK));
9346 u = (tma_mp_digit)(r >> ((tma_mp_word) DIGIT_BIT));
9348 /* left hand side of A[ix] * A[iy] */
9351 /* alias for where to store the results */
9352 tmpt = t.dp + (2*ix + 1);
9354 for (iy = ix + 1; iy < pa; iy++) {
9355 /* first calculate the product */
9356 r = ((tma_mp_word)tmpx) * ((tma_mp_word)a->dp[iy]);
9358 /* now calculate the double precision result, note we use
9359 * addition instead of *2 since it's easier to optimize
9361 r = ((tma_mp_word) *tmpt) + r + r + ((tma_mp_word) u);
9363 /* store lower part */
9364 *tmpt++ = (tma_mp_digit) (r & ((tma_mp_word) MP_MASK));
9367 u = (tma_mp_digit)(r >> ((tma_mp_word) DIGIT_BIT));
9369 /* propagate upwards */
9370 while (u != ((tma_mp_digit) 0)) {
9371 r = ((tma_mp_word) *tmpt) + ((tma_mp_word) u);
9372 *tmpt++ = (tma_mp_digit) (r & ((tma_mp_word) MP_MASK));
9373 u = (tma_mp_digit)(r >> ((tma_mp_word) DIGIT_BIT));
9378 tma_mp_exch (&t, b);
9388 /* End: bn_s_tma_mp_sqr.c */
9390 /* Start: bn_s_tma_mp_sub.c */
9392 #ifdef BN_S_MP_SUB_C
9393 /* LibTomMath, multiple-precision integer library -- Tom St Denis
9395 * LibTomMath is a library that provides multiple-precision
9396 * integer arithmetic as well as number theoretic functionality.
9398 * The library was designed directly after the MPI library by
9399 * Michael Fromberger but has been written from scratch with
9400 * additional optimizations in place.
9402 * The library is free for all purposes without any express
9403 * guarantee it works.
9405 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
9408 /* low level subtraction (assumes |a| > |b|), HAC pp.595 Algorithm 14.9 */
9410 s_tma_mp_sub (tma_mp_int * a, tma_mp_int * b, tma_mp_int * c)
9412 int olduse, res, min, max;
9419 if (c->alloc < max) {
9420 if ((res = tma_mp_grow (c, max)) != MP_OKAY) {
9428 register tma_mp_digit u, *tmpa, *tmpb, *tmpc;
9431 /* alias for digit pointers */
9436 /* set carry to zero */
9438 for (i = 0; i < min; i++) {
9439 /* T[i] = A[i] - B[i] - U */
9440 *tmpc = *tmpa++ - *tmpb++ - u;
9442 /* U = carry bit of T[i]
9443 * Note this saves performing an AND operation since
9444 * if a carry does occur it will propagate all the way to the
9445 * MSB. As a result a single shift is enough to get the carry
9447 u = *tmpc >> ((tma_mp_digit)(CHAR_BIT * sizeof (tma_mp_digit) - 1));
9449 /* Clear carry from T[i] */
9453 /* now copy higher words if any, e.g. if A has more digits than B */
9454 for (; i < max; i++) {
9455 /* T[i] = A[i] - U */
9456 *tmpc = *tmpa++ - u;
9458 /* U = carry bit of T[i] */
9459 u = *tmpc >> ((tma_mp_digit)(CHAR_BIT * sizeof (tma_mp_digit) - 1));
9461 /* Clear carry from T[i] */
9465 /* clear digits above used (since we may not have grown result above) */
9466 for (i = c->used; i < olduse; i++) {
9481 /* End: bn_s_tma_mp_sub.c */
9483 /* Start: bncore.c */
9486 /* LibTomMath, multiple-precision integer library -- Tom St Denis
9488 * LibTomMath is a library that provides multiple-precision
9489 * integer arithmetic as well as number theoretic functionality.
9491 * The library was designed directly after the MPI library by
9492 * Michael Fromberger but has been written from scratch with
9493 * additional optimizations in place.
9495 * The library is free for all purposes without any express
9496 * guarantee it works.
9498 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
9501 /* Known optimal configurations
9503 CPU /Compiler /MUL CUTOFF/SQR CUTOFF
9504 -------------------------------------------------------------
9505 Intel P4 Northwood /GCC v3.4.1 / 88/ 128/LTM 0.32 ;-)
9506 AMD Athlon64 /GCC v3.4.4 / 80/ 120/LTM 0.35
9510 int KARATSUBA_MUL_CUTOFF = 80, /* Min. number of digits before Karatsuba multiplication is used. */
9511 KARATSUBA_SQR_CUTOFF = 120, /* Min. number of digits before Karatsuba squaring is used. */
9513 TOOM_MUL_CUTOFF = 350, /* no optimal values of these are known yet so set em high */
9514 TOOM_SQR_CUTOFF = 400;