3 * Copyright (C) 2006-2012 Sawtooth Consulting Ltd.
5 * This file is part of CyaSSL.
7 * CyaSSL is free software; you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation; either version 2 of the License, or
10 * (at your option) any later version.
12 * CyaSSL is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 * GNU General Public License for more details.
17 * You should have received a copy of the GNU General Public License
18 * along with this program; if not, write to the Free Software
19 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA
24 * Based on public domain TomsFastMath 0.10 by Tom St Denis, tomstdenis@iahu.ca,
25 * http://math.libtomcrypt.com
29 * Edited by Moisés Guimarães (moises.guimaraes@phoebus.com.br)
30 * to fit CyaSSL's needs.
37 /* in case user set USE_FAST_MATH there */
38 #include <cyassl/ctaocrypt/settings.h>
42 #include <cyassl/ctaocrypt/tfm.h>
43 #include <ctaocrypt/src/asm.c> /* will define asm MACROS or C ones */
46 /* math settings check */
47 word32 CheckRunTimeSettings(void)
53 /* math settings size check */
54 word32 CheckRunTimeFastMath(void)
62 void fp_add(fp_int *a, fp_int *b, fp_int *c)
66 /* get sign of both inputs */
70 /* handle two cases, not four */
72 /* both positive or both negative */
73 /* add their magnitudes, copy the sign */
77 /* one positive, the other negative */
78 /* subtract the one with the greater magnitude from */
79 /* the one of the lesser magnitude. The result gets */
80 /* the sign of the one with the greater magnitude. */
81 if (fp_cmp_mag (a, b) == FP_LT) {
91 /* unsigned addition */
92 void s_fp_add(fp_int *a, fp_int *b, fp_int *c)
97 y = MAX(a->used, b->used);
102 for (x = 0; x < y; x++) {
103 t += ((fp_word)a->dp[x]) + ((fp_word)b->dp[x]);
104 c->dp[x] = (fp_digit)t;
107 if (t != 0 && x < FP_SIZE) {
108 c->dp[c->used++] = (fp_digit)t;
113 for (; x < oldused; x++) {
120 void fp_sub(fp_int *a, fp_int *b, fp_int *c)
128 /* subtract a negative from a positive, OR */
129 /* subtract a positive from a negative. */
130 /* In either case, ADD their magnitudes, */
131 /* and use the sign of the first number. */
135 /* subtract a positive from a positive, OR */
136 /* subtract a negative from a negative. */
137 /* First, take the difference between their */
138 /* magnitudes, then... */
139 if (fp_cmp_mag (a, b) != FP_LT) {
140 /* Copy the sign from the first */
142 /* The first has a larger or equal magnitude */
145 /* The result has the *opposite* sign from */
146 /* the first number. */
147 c->sign = (sa == FP_ZPOS) ? FP_NEG : FP_ZPOS;
148 /* The second has a larger magnitude */
154 /* unsigned subtraction ||a|| >= ||b|| ALWAYS! */
155 void s_fp_sub(fp_int *a, fp_int *b, fp_int *c)
157 int x, oldbused, oldused;
164 for (x = 0; x < oldbused; x++) {
165 t = ((fp_word)a->dp[x]) - (((fp_word)b->dp[x]) + t);
166 c->dp[x] = (fp_digit)t;
167 t = (t >> DIGIT_BIT)&1;
169 for (; x < a->used; x++) {
170 t = ((fp_word)a->dp[x]) - t;
171 c->dp[x] = (fp_digit)t;
172 t = (t >> DIGIT_BIT);
174 for (; x < oldused; x++) {
181 void fp_mul(fp_int *A, fp_int *B, fp_int *C)
185 y = MAX(A->used, B->used);
186 yy = MIN(A->used, B->used);
188 /* call generic if we're out of range */
189 if (y + yy > FP_SIZE) {
190 fp_mul_comba(A, B, C);
194 /* pick a comba (unrolled 4/8/16/32 x or rolled) based on the size
195 of the largest input. We also want to avoid doing excess mults if the
196 inputs are not close to the next power of two. That is, for example,
197 if say y=17 then we would do (32-17)^2 = 225 unneeded multiplications
202 fp_mul_comba3(A,B,C);
208 fp_mul_comba4(A,B,C);
214 fp_mul_comba6(A,B,C);
220 fp_mul_comba7(A,B,C);
226 fp_mul_comba8(A,B,C);
232 fp_mul_comba9(A,B,C);
238 fp_mul_comba12(A,B,C);
244 fp_mul_comba17(A,B,C);
251 fp_mul_comba_small(A,B,C);
255 #if defined(TFM_MUL20)
257 fp_mul_comba20(A,B,C);
261 #if defined(TFM_MUL24)
262 if (yy >= 16 && y <= 24) {
263 fp_mul_comba24(A,B,C);
267 #if defined(TFM_MUL28)
268 if (yy >= 20 && y <= 28) {
269 fp_mul_comba28(A,B,C);
273 #if defined(TFM_MUL32)
274 if (yy >= 24 && y <= 32) {
275 fp_mul_comba32(A,B,C);
279 #if defined(TFM_MUL48)
280 if (yy >= 40 && y <= 48) {
281 fp_mul_comba48(A,B,C);
285 #if defined(TFM_MUL64)
286 if (yy >= 56 && y <= 64) {
287 fp_mul_comba64(A,B,C);
294 void fp_mul_2(fp_int * a, fp_int * b)
302 register fp_digit r, rr, *tmpa, *tmpb;
304 /* alias for source */
312 for (x = 0; x < a->used; x++) {
314 /* get what will be the *next* carry bit from the
315 * MSB of the current digit
317 rr = *tmpa >> ((fp_digit)(DIGIT_BIT - 1));
319 /* now shift up this digit, add in the carry [from the previous] */
320 *tmpb++ = ((*tmpa++ << ((fp_digit)1)) | r);
322 /* copy the carry that would be from the source
323 * digit into the next iteration
328 /* new leading digit? */
329 if (r != 0 && b->used != (FP_SIZE-1)) {
330 /* add a MSB which is always 1 at this point */
335 /* now zero any excess digits on the destination
336 * that we didn't write to
338 tmpb = b->dp + b->used;
339 for (x = b->used; x < oldused; x++) {
347 void fp_mul_d(fp_int *a, fp_digit b, fp_int *c)
356 for (x = 0; x < a->used; x++) {
357 w = ((fp_word)a->dp[x]) * ((fp_word)b) + w;
358 c->dp[x] = (fp_digit)w;
361 if (w != 0 && (a->used != FP_SIZE)) {
362 c->dp[c->used++] = (fp_digit) w;
365 for (; x < oldused; x++) {
372 void fp_mul_2d(fp_int *a, int b, fp_int *c)
374 fp_digit carry, carrytmp, shift;
380 /* handle whole digits */
381 if (b >= DIGIT_BIT) {
382 fp_lshd(c, b/DIGIT_BIT);
386 /* shift the digits */
389 shift = DIGIT_BIT - b;
390 for (x = 0; x < c->used; x++) {
391 carrytmp = c->dp[x] >> shift;
392 c->dp[x] = (c->dp[x] << b) + carry;
395 /* store last carry if room */
396 if (carry && x < FP_SIZE) {
397 c->dp[c->used++] = carry;
403 /* generic PxQ multiplier */
404 void fp_mul_comba(fp_int *A, fp_int *B, fp_int *C)
406 int ix, iy, iz, tx, ty, pa;
407 fp_digit c0, c1, c2, *tmpx, *tmpy;
413 /* get size of output and trim */
414 pa = A->used + B->used;
419 if (A == C || B == C) {
427 for (ix = 0; ix < pa; ix++) {
428 /* get offsets into the two bignums */
429 ty = MIN(ix, B->used-1);
432 /* setup temp aliases */
436 /* this is the number of times the loop will iterrate, essentially its
437 while (tx++ < a->used && ty-- >= 0) { ... }
439 iy = MIN(A->used-tx, ty+1);
443 for (iz = 0; iz < iy; ++iz) {
444 /* TAO change COMBA_ADD back to MULADD */
445 MULADD(*tmpx++, *tmpy--);
449 COMBA_STORE(dst->dp[ix]);
454 dst->sign = A->sign ^ B->sign;
459 /* a/b => cb + d == a */
460 int fp_div(fp_int *a, fp_int *b, fp_int *c, fp_int *d)
462 fp_int q, x, y, t1, t2;
463 int n, t, i, norm, neg;
465 /* is divisor zero ? */
466 if (fp_iszero (b) == 1) {
470 /* if a < b then q=0, r = a */
471 if (fp_cmp_mag (a, b) == FP_LT) {
482 q.used = a->used + 2;
490 neg = (a->sign == b->sign) ? FP_ZPOS : FP_NEG;
491 x.sign = y.sign = FP_ZPOS;
493 /* normalize both x and y, ensure that y >= b/2, [b == 2**DIGIT_BIT] */
494 norm = fp_count_bits(&y) % DIGIT_BIT;
495 if (norm < (int)(DIGIT_BIT-1)) {
496 norm = (DIGIT_BIT-1) - norm;
497 fp_mul_2d (&x, norm, &x);
498 fp_mul_2d (&y, norm, &y);
503 /* note hac does 0 based, so if used==5 then its 0,1,2,3,4, e.g. use 4 */
507 /* while (x >= y*b**n-t) do { q[n-t] += 1; x -= y*b**{n-t} } */
508 fp_lshd (&y, n - t); /* y = y*b**{n-t} */
510 while (fp_cmp (&x, &y) != FP_LT) {
515 /* reset y by shifting it back down */
518 /* step 3. for i from n down to (t + 1) */
519 for (i = n; i >= (t + 1); i--) {
524 /* step 3.1 if xi == yt then set q{i-t-1} to b-1,
525 * otherwise set q{i-t-1} to (xi*b + x{i-1})/yt */
526 if (x.dp[i] == y.dp[t]) {
527 q.dp[i - t - 1] = ((((fp_word)1) << DIGIT_BIT) - 1);
530 tmp = ((fp_word) x.dp[i]) << ((fp_word) DIGIT_BIT);
531 tmp |= ((fp_word) x.dp[i - 1]);
532 tmp /= ((fp_word)y.dp[t]);
533 q.dp[i - t - 1] = (fp_digit) (tmp);
536 /* while (q{i-t-1} * (yt * b + y{t-1})) >
537 xi * b**2 + xi-1 * b + xi-2
541 q.dp[i - t - 1] = (q.dp[i - t - 1] + 1);
543 q.dp[i - t - 1] = (q.dp[i - t - 1] - 1);
547 t1.dp[0] = (t - 1 < 0) ? 0 : y.dp[t - 1];
550 fp_mul_d (&t1, q.dp[i - t - 1], &t1);
552 /* find right hand */
553 t2.dp[0] = (i - 2 < 0) ? 0 : x.dp[i - 2];
554 t2.dp[1] = (i - 1 < 0) ? 0 : x.dp[i - 1];
557 } while (fp_cmp_mag(&t1, &t2) == FP_GT);
559 /* step 3.3 x = x - q{i-t-1} * y * b**{i-t-1} */
560 fp_mul_d (&y, q.dp[i - t - 1], &t1);
561 fp_lshd (&t1, i - t - 1);
562 fp_sub (&x, &t1, &x);
564 /* if x < 0 then { x = x + y*b**{i-t-1}; q{i-t-1} -= 1; } */
565 if (x.sign == FP_NEG) {
567 fp_lshd (&t1, i - t - 1);
568 fp_add (&x, &t1, &x);
569 q.dp[i - t - 1] = q.dp[i - t - 1] - 1;
573 /* now q is the quotient and x is the remainder
574 * [which we have to normalize]
577 /* get sign before writing to c */
578 x.sign = x.used == 0 ? FP_ZPOS : a->sign;
587 fp_div_2d (&x, norm, &x, NULL);
589 /* the following is a kludge, essentially we were seeing the right remainder but
590 with excess digits that should have been zero
592 for (i = b->used; i < x.used; i++) {
603 void fp_div_2(fp_int * a, fp_int * b)
610 register fp_digit r, rr, *tmpa, *tmpb;
613 tmpa = a->dp + b->used - 1;
616 tmpb = b->dp + b->used - 1;
620 for (x = b->used - 1; x >= 0; x--) {
621 /* get the carry for the next iteration */
624 /* shift the current digit, add in carry and store */
625 *tmpb-- = (*tmpa-- >> 1) | (r << (DIGIT_BIT - 1));
627 /* forward carry to next iteration */
631 /* zero excess digits */
632 tmpb = b->dp + b->used;
633 for (x = b->used; x < oldused; x++) {
642 void fp_div_2d(fp_int *a, int b, fp_int *c, fp_int *d)
648 /* if the shift count is <= 0 then we do no work */
659 /* get the remainder */
661 fp_mod_2d (a, b, &t);
667 /* shift by as many digits in the bit count */
668 if (b >= (int)DIGIT_BIT) {
669 fp_rshd (c, b / DIGIT_BIT);
672 /* shift any bit count < DIGIT_BIT */
673 D = (fp_digit) (b % DIGIT_BIT);
675 register fp_digit *tmpc, mask, shift;
678 mask = (((fp_digit)1) << D) - 1;
681 shift = DIGIT_BIT - D;
684 tmpc = c->dp + (c->used - 1);
688 for (x = c->used - 1; x >= 0; x--) {
689 /* get the lower bits of this word in a temp */
692 /* shift the current word and mix in the carry bits from the previous word */
693 *tmpc = (*tmpc >> D) | (r << shift);
696 /* set the carry to the carry bits of the current word found above */
706 /* c = a mod b, 0 <= c < b */
707 int fp_mod(fp_int *a, fp_int *b, fp_int *c)
713 if ((err = fp_div(a, b, NULL, &t)) != FP_OKAY) {
716 if (t.sign != b->sign) {
725 void fp_mod_2d(fp_int *a, int b, fp_int *c)
729 /* zero if count less than or equal to zero */
735 /* get copy of input */
738 /* if 2**d is larger than we just return */
739 if (b >= (DIGIT_BIT * a->used)) {
743 /* zero digits above the last digit of the modulus */
744 for (x = (b / DIGIT_BIT) + ((b % DIGIT_BIT) == 0 ? 0 : 1); x < c->used; x++) {
747 /* clear the digit that is not completely outside/inside the modulus */
748 c->dp[b / DIGIT_BIT] &= ~((fp_digit)0) >> (DIGIT_BIT - b);
752 static int fp_invmod_slow (fp_int * a, fp_int * b, fp_int * c)
754 fp_int x, y, u, v, A, B, C, D;
757 /* b cannot be negative */
758 if (b->sign == FP_NEG || fp_iszero(b) == 1) {
763 fp_init(&x); fp_init(&y);
764 fp_init(&u); fp_init(&v);
765 fp_init(&A); fp_init(&B);
766 fp_init(&C); fp_init(&D);
769 if ((res = fp_mod(a, b, &x)) != FP_OKAY) {
774 /* 2. [modified] if x,y are both even then return an error! */
775 if (fp_iseven (&x) == 1 && fp_iseven (&y) == 1) {
779 /* 3. u=x, v=y, A=1, B=0, C=0,D=1 */
786 /* 4. while u is even do */
787 while (fp_iseven (&u) == 1) {
791 /* 4.2 if A or B is odd then */
792 if (fp_isodd (&A) == 1 || fp_isodd (&B) == 1) {
793 /* A = (A+y)/2, B = (B-x)/2 */
797 /* A = A/2, B = B/2 */
802 /* 5. while v is even do */
803 while (fp_iseven (&v) == 1) {
807 /* 5.2 if C or D is odd then */
808 if (fp_isodd (&C) == 1 || fp_isodd (&D) == 1) {
809 /* C = (C+y)/2, D = (D-x)/2 */
813 /* C = C/2, D = D/2 */
818 /* 6. if u >= v then */
819 if (fp_cmp (&u, &v) != FP_LT) {
820 /* u = u - v, A = A - C, B = B - D */
825 /* v - v - u, C = C - A, D = D - B */
831 /* if not zero goto step 4 */
832 if (fp_iszero (&u) == 0)
835 /* now a = C, b = D, gcd == g*v */
837 /* if v != 1 then there is no inverse */
838 if (fp_cmp_d (&v, 1) != FP_EQ) {
843 while (fp_cmp_d(&C, 0) == FP_LT) {
848 while (fp_cmp_mag(&C, b) != FP_LT) {
852 /* C is now the inverse */
857 /* c = 1/a (mod b) for odd b only */
858 int fp_invmod(fp_int *a, fp_int *b, fp_int *c)
860 fp_int x, y, u, v, B, D;
863 /* 2. [modified] b must be odd */
864 if (fp_iseven (b) == FP_YES) {
865 return fp_invmod_slow(a,b,c);
868 /* init all our temps */
869 fp_init(&x); fp_init(&y);
870 fp_init(&u); fp_init(&v);
871 fp_init(&B); fp_init(&D);
873 /* x == modulus, y == value to invert */
876 /* we need y = |a| */
879 /* 3. u=x, v=y, A=1, B=0, C=0,D=1 */
885 /* 4. while u is even do */
886 while (fp_iseven (&u) == FP_YES) {
890 /* 4.2 if B is odd then */
891 if (fp_isodd (&B) == FP_YES) {
898 /* 5. while v is even do */
899 while (fp_iseven (&v) == FP_YES) {
903 /* 5.2 if D is odd then */
904 if (fp_isodd (&D) == FP_YES) {
912 /* 6. if u >= v then */
913 if (fp_cmp (&u, &v) != FP_LT) {
914 /* u = u - v, B = B - D */
918 /* v - v - u, D = D - B */
923 /* if not zero goto step 4 */
924 if (fp_iszero (&u) == FP_NO) {
928 /* now a = C, b = D, gcd == g*v */
930 /* if v != 1 then there is no inverse */
931 if (fp_cmp_d (&v, 1) != FP_EQ) {
935 /* b is now the inverse */
937 while (D.sign == FP_NEG) {
945 /* d = a * b (mod c) */
946 int fp_mulmod(fp_int *a, fp_int *b, fp_int *c, fp_int *d)
951 return fp_mod(&tmp, c, d);
954 #ifdef TFM_TIMING_RESISTANT
956 /* timing resistant montgomery ladder based exptmod
958 Based on work by Marc Joye, Sung-Ming Yen, "The Montgomery Powering Ladder", Cryptographic Hardware and Embedded Systems, CHES 2002
960 static int _fp_exptmod(fp_int * G, fp_int * X, fp_int * P, fp_int * Y)
964 int err, bitcnt, digidx, y;
966 /* now setup montgomery */
967 if ((err = fp_montgomery_setup (P, &mp)) != FP_OKAY) {
974 /* now we need R mod m */
975 fp_montgomery_calc_normalization (&R[0], P);
977 /* now set R[0][1] to G * R mod m */
978 if (fp_cmp_mag(P, G) != FP_GT) {
979 /* G > P so we reduce it first */
984 fp_mulmod (&R[1], &R[0], P, &R[1]);
986 /* for j = t-1 downto 0 do
987 r_!k = R0*R1; r_k = r_k^2
990 /* set initial mode and bit cnt */
993 digidx = X->used - 1;
996 /* grab next digit as required */
998 /* if digidx == -1 we are out of digits so break */
1002 /* read next digit and reset bitcnt */
1003 buf = X->dp[digidx--];
1004 bitcnt = (int)DIGIT_BIT;
1007 /* grab the next msb from the exponent */
1008 y = (fp_digit)(buf >> (DIGIT_BIT - 1)) & 1;
1009 buf <<= (fp_digit)1;
1012 fp_mul(&R[0], &R[1], &R[y^1]); fp_montgomery_reduce(&R[y^1], P, mp);
1013 fp_sqr(&R[y], &R[y]); fp_montgomery_reduce(&R[y], P, mp);
1016 fp_montgomery_reduce(&R[0], P, mp);
1024 * Some restrictions... x must be positive and < b
1026 static int _fp_exptmod(fp_int * G, fp_int * X, fp_int * P, fp_int * Y)
1030 int err, bitbuf, bitcpy, bitcnt, mode, digidx, x, y, winsize;
1032 /* find window size */
1033 x = fp_count_bits (X);
1036 } else if (x <= 36) {
1038 } else if (x <= 140) {
1040 } else if (x <= 450) {
1047 XMEMSET(M, 0, sizeof(M));
1049 /* now setup montgomery */
1050 if ((err = fp_montgomery_setup (P, &mp)) != FP_OKAY) {
1059 * The M table contains powers of the input base, e.g. M[x] = G^x mod P
1061 * The first half of the table is not computed though accept for M[0] and M[1]
1064 /* now we need R mod m */
1065 fp_montgomery_calc_normalization (&res, P);
1067 /* now set M[1] to G * R mod m */
1068 if (fp_cmp_mag(P, G) != FP_GT) {
1069 /* G > P so we reduce it first */
1070 fp_mod(G, P, &M[1]);
1074 fp_mulmod (&M[1], &res, P, &M[1]);
1076 /* compute the value at M[1<<(winsize-1)] by squaring M[1] (winsize-1) times */
1077 fp_copy (&M[1], &M[1 << (winsize - 1)]);
1078 for (x = 0; x < (winsize - 1); x++) {
1079 fp_sqr (&M[1 << (winsize - 1)], &M[1 << (winsize - 1)]);
1080 fp_montgomery_reduce (&M[1 << (winsize - 1)], P, mp);
1083 /* create upper table */
1084 for (x = (1 << (winsize - 1)) + 1; x < (1 << winsize); x++) {
1085 fp_mul(&M[x - 1], &M[1], &M[x]);
1086 fp_montgomery_reduce(&M[x], P, mp);
1089 /* set initial mode and bit cnt */
1093 digidx = X->used - 1;
1098 /* grab next digit as required */
1099 if (--bitcnt == 0) {
1100 /* if digidx == -1 we are out of digits so break */
1104 /* read next digit and reset bitcnt */
1105 buf = X->dp[digidx--];
1106 bitcnt = (int)DIGIT_BIT;
1109 /* grab the next msb from the exponent */
1110 y = (fp_digit)(buf >> (DIGIT_BIT - 1)) & 1;
1111 buf <<= (fp_digit)1;
1113 /* if the bit is zero and mode == 0 then we ignore it
1114 * These represent the leading zero bits before the first 1 bit
1115 * in the exponent. Technically this opt is not required but it
1116 * does lower the # of trivial squaring/reductions used
1118 if (mode == 0 && y == 0) {
1122 /* if the bit is zero and mode == 1 then we square */
1123 if (mode == 1 && y == 0) {
1125 fp_montgomery_reduce(&res, P, mp);
1129 /* else we add it to the window */
1130 bitbuf |= (y << (winsize - ++bitcpy));
1133 if (bitcpy == winsize) {
1134 /* ok window is filled so square as required and multiply */
1136 for (x = 0; x < winsize; x++) {
1138 fp_montgomery_reduce(&res, P, mp);
1142 fp_mul(&res, &M[bitbuf], &res);
1143 fp_montgomery_reduce(&res, P, mp);
1145 /* empty window and reset */
1152 /* if bits remain then square/multiply */
1153 if (mode == 2 && bitcpy > 0) {
1154 /* square then multiply if the bit is set */
1155 for (x = 0; x < bitcpy; x++) {
1157 fp_montgomery_reduce(&res, P, mp);
1159 /* get next bit of the window */
1161 if ((bitbuf & (1 << winsize)) != 0) {
1163 fp_mul(&res, &M[1], &res);
1164 fp_montgomery_reduce(&res, P, mp);
1169 /* fixup result if Montgomery reduction is used
1170 * recall that any value in a Montgomery system is
1171 * actually multiplied by R mod n. So we have
1172 * to reduce one more time to cancel out the factor
1175 fp_montgomery_reduce(&res, P, mp);
1177 /* swap res with Y */
1184 int fp_exptmod(fp_int * G, fp_int * X, fp_int * P, fp_int * Y)
1189 /* prevent overflows */
1190 if (P->used > (FP_SIZE/2)) {
1194 /* is X negative? */
1195 if (X->sign == FP_NEG) {
1196 /* yes, copy G and invmod it */
1198 if ((err = fp_invmod(&tmp, P, &tmp)) != FP_OKAY) {
1202 err = _fp_exptmod(&tmp, X, P, Y);
1208 /* Positive exponent so just exptmod */
1209 return _fp_exptmod(G, X, P, Y);
1213 /* computes a = 2**b */
1214 void fp_2expt(fp_int *a, int b)
1218 /* zero a as per default */
1230 /* set the used count of where the bit will go */
1233 /* put the single bit in its place */
1234 a->dp[z] = ((fp_digit)1) << (b % DIGIT_BIT);
1238 void fp_sqr(fp_int *A, fp_int *B)
1242 /* call generic if we're out of range */
1243 if (y + y > FP_SIZE) {
1248 #if defined(TFM_SQR3)
1254 #if defined(TFM_SQR4)
1260 #if defined(TFM_SQR6)
1266 #if defined(TFM_SQR7)
1272 #if defined(TFM_SQR8)
1278 #if defined(TFM_SQR9)
1284 #if defined(TFM_SQR12)
1286 fp_sqr_comba12(A,B);
1290 #if defined(TFM_SQR17)
1292 fp_sqr_comba17(A,B);
1296 #if defined(TFM_SMALL_SET)
1298 fp_sqr_comba_small(A,B);
1302 #if defined(TFM_SQR20)
1304 fp_sqr_comba20(A,B);
1308 #if defined(TFM_SQR24)
1310 fp_sqr_comba24(A,B);
1314 #if defined(TFM_SQR28)
1316 fp_sqr_comba28(A,B);
1320 #if defined(TFM_SQR32)
1322 fp_sqr_comba32(A,B);
1326 #if defined(TFM_SQR48)
1328 fp_sqr_comba48(A,B);
1332 #if defined(TFM_SQR64)
1334 fp_sqr_comba64(A,B);
1341 /* generic comba squarer */
1342 void fp_sqr_comba(fp_int *A, fp_int *B)
1345 fp_digit c0, c1, c2;
1351 /* get size of output and trim */
1352 pa = A->used + A->used;
1353 if (pa >= FP_SIZE) {
1357 /* number of output digits to produce */
1369 for (ix = 0; ix < pa; ix++) {
1371 fp_digit *tmpy, *tmpx;
1373 /* get offsets into the two bignums */
1374 ty = MIN(A->used-1, ix);
1377 /* setup temp aliases */
1381 /* this is the number of times the loop will iterrate,
1382 while (tx++ < a->used && ty-- >= 0) { ... }
1384 iy = MIN(A->used-tx, ty+1);
1386 /* now for squaring tx can never equal ty
1387 * we halve the distance since they approach
1388 * at a rate of 2x and we have to round because
1389 * odd cases need to be executed
1391 iy = MIN(iy, (ty-tx+1)>>1);
1393 /* forward carries */
1397 for (iz = 0; iz < iy; iz++) {
1398 SQRADD2(*tmpx++, *tmpy--);
1401 /* even columns have the square term in them */
1403 /* TAO change COMBA_ADD back to SQRADD */
1404 SQRADD(A->dp[ix>>1], A->dp[ix>>1]);
1408 COMBA_STORE(dst->dp[ix]);
1421 int fp_cmp(fp_int *a, fp_int *b)
1423 if (a->sign == FP_NEG && b->sign == FP_ZPOS) {
1425 } else if (a->sign == FP_ZPOS && b->sign == FP_NEG) {
1428 /* compare digits */
1429 if (a->sign == FP_NEG) {
1430 /* if negative compare opposite direction */
1431 return fp_cmp_mag(b, a);
1433 return fp_cmp_mag(a, b);
1438 /* compare against a single digit */
1439 int fp_cmp_d(fp_int *a, fp_digit b)
1441 /* compare based on sign */
1442 if ((b && a->used == 0) || a->sign == FP_NEG) {
1446 /* compare based on magnitude */
1451 /* compare the only digit of a to b */
1454 } else if (a->dp[0] < b) {
1462 int fp_cmp_mag(fp_int *a, fp_int *b)
1466 if (a->used > b->used) {
1468 } else if (a->used < b->used) {
1471 for (x = a->used - 1; x >= 0; x--) {
1472 if (a->dp[x] > b->dp[x]) {
1474 } else if (a->dp[x] < b->dp[x]) {
1482 /* setups the montgomery reduction */
1483 int fp_montgomery_setup(fp_int *a, fp_digit *rho)
1487 /* fast inversion mod 2**k
1489 * Based on the fact that
1491 * XA = 1 (mod 2**n) => (X(2-XA)) A = 1 (mod 2**2n)
1492 * => 2*X*A - X*X*A*A = 1
1493 * => 2*(1) - (1) = 1
1501 x = (((b + 2) & 4) << 1) + b; /* here x*a==1 mod 2**4 */
1502 x *= 2 - b * x; /* here x*a==1 mod 2**8 */
1503 x *= 2 - b * x; /* here x*a==1 mod 2**16 */
1504 x *= 2 - b * x; /* here x*a==1 mod 2**32 */
1506 x *= 2 - b * x; /* here x*a==1 mod 2**64 */
1509 /* rho = -1/m mod b */
1510 *rho = (fp_digit) (((fp_word) 1 << ((fp_word) DIGIT_BIT)) - ((fp_word)x));
1515 /* computes a = B**n mod b without division or multiplication useful for
1516 * normalizing numbers in a Montgomery system.
1518 void fp_montgomery_calc_normalization(fp_int *a, fp_int *b)
1522 /* how many bits of last digit does b use */
1523 bits = fp_count_bits (b) % DIGIT_BIT;
1524 if (!bits) bits = DIGIT_BIT;
1526 /* compute A = B^(n-1) * 2^(bits-1) */
1528 fp_2expt (a, (b->used - 1) * DIGIT_BIT + bits - 1);
1534 /* now compute C = A * B mod b */
1535 for (x = bits - 1; x < (int)DIGIT_BIT; x++) {
1537 if (fp_cmp_mag (a, b) != FP_LT) {
1544 #ifdef TFM_SMALL_MONT_SET
1545 #include "fp_mont_small.i"
1548 /* computes x/R == x (mod N) via Montgomery Reduction */
1549 void fp_montgomery_reduce(fp_int *a, fp_int *m, fp_digit mp)
1551 fp_digit c[FP_SIZE], *_c, *tmpm, mu;
1552 int oldused, x, y, pa;
1554 /* bail if too large */
1555 if (m->used > (FP_SIZE/2)) {
1556 (void)mu; /* shut up compiler */
1560 #ifdef TFM_SMALL_MONT_SET
1561 if (m->used <= 16) {
1562 fp_montgomery_reduce_small(a, m, mp);
1568 #if defined(USE_MEMSET)
1569 /* now zero the buff */
1570 XMEMSET(c, 0, sizeof c);
1574 /* copy the input */
1576 for (x = 0; x < oldused; x++) {
1579 #if !defined(USE_MEMSET)
1580 for (; x < 2*pa+1; x++) {
1586 for (x = 0; x < pa; x++) {
1588 /* get Mu for this round */
1593 #if (defined(TFM_SSE2) || defined(TFM_X86_64))
1594 for (; y < (pa & ~7); y += 8) {
1601 for (; y < pa; y++) {
1615 for (x = 0; x < pa+1; x++) {
1619 for (; x < oldused; x++) {
1628 /* if A >= m then A = A - m */
1629 if (fp_cmp_mag (a, m) != FP_LT) {
1634 void fp_read_unsigned_bin(fp_int *a, unsigned char *b, int c)
1639 /* If we know the endianness of this architecture, and we're using
1640 32-bit fp_digits, we can optimize this */
1641 #if (defined(ENDIAN_LITTLE) || defined(ENDIAN_BIG)) && !defined(FP_64BIT)
1642 /* But not for both simultaneously */
1643 #if defined(ENDIAN_LITTLE) && defined(ENDIAN_BIG)
1644 #error Both ENDIAN_LITTLE and ENDIAN_BIG defined.
1647 unsigned char *pd = (unsigned char *)a->dp;
1649 if ((unsigned)c > (FP_SIZE * sizeof(fp_digit))) {
1650 int excess = c - (FP_SIZE * sizeof(fp_digit));
1654 a->used = (c + sizeof(fp_digit) - 1)/sizeof(fp_digit);
1655 /* read the bytes in */
1658 /* Use Duff's device to unroll the loop. */
1659 int idx = (c - 1) & ~3;
1661 case 0: do { pd[idx+0] = *b++;
1662 case 3: pd[idx+1] = *b++;
1663 case 2: pd[idx+2] = *b++;
1664 case 1: pd[idx+3] = *b++;
1666 } while ((c -= 4) > 0);
1670 for (c -= 1; c >= 0; c -= 1) {
1676 /* read the bytes in */
1677 for (; c > 0; c--) {
1678 fp_mul_2d (a, 8, a);
1686 void fp_to_unsigned_bin(fp_int *a, unsigned char *b)
1691 fp_init_copy(&t, a);
1694 while (fp_iszero (&t) == FP_NO) {
1695 b[x++] = (unsigned char) (t.dp[0] & 255);
1696 fp_div_2d (&t, 8, &t, NULL);
1701 int fp_unsigned_bin_size(fp_int *a)
1703 int size = fp_count_bits (a);
1704 return (size / 8 + ((size & 7) != 0 ? 1 : 0));
1707 void fp_set(fp_int *a, fp_digit b)
1711 a->used = a->dp[0] ? 1 : 0;
1714 int fp_count_bits (fp_int * a)
1724 /* get number of digits and add that */
1725 r = (a->used - 1) * DIGIT_BIT;
1727 /* take the last digit and count the bits in it */
1728 q = a->dp[a->used - 1];
1729 while (q > ((fp_digit) 0)) {
1731 q >>= ((fp_digit) 1);
1736 void fp_lshd(fp_int *a, int x)
1740 /* move up and truncate as required */
1741 y = MIN(a->used + x - 1, (int)(FP_SIZE-1));
1743 /* store new size */
1747 for (; y >= x; y--) {
1748 a->dp[y] = a->dp[y-x];
1751 /* zero lower digits */
1752 for (; y >= 0; y--) {
1760 void fp_rshd(fp_int *a, int x)
1764 /* too many digits just zero and return */
1771 for (y = 0; y < a->used - x; y++) {
1772 a->dp[y] = a->dp[y+x];
1776 for (; y < a->used; y++) {
1780 /* decrement count */
1785 /* reverse an array, used for radix code */
1786 void fp_reverse (unsigned char *s, int len)
1804 void fp_sub_d(fp_int *a, fp_digit b, fp_int *c)
1812 /* CyaSSL callers from normal lib */
1814 /* init a new mp_int */
1815 int mp_init (mp_int * a)
1822 /* clear one (frees) */
1823 void mp_clear (mp_int * a)
1828 /* handle up to 6 inits */
1829 int mp_init_multi(mp_int* a, mp_int* b, mp_int* c, mp_int* d, mp_int* e, mp_int* f)
1847 /* high level addition (handles signs) */
1848 int mp_add (mp_int * a, mp_int * b, mp_int * c)
1854 /* high level subtraction (handles signs) */
1855 int mp_sub (mp_int * a, mp_int * b, mp_int * c)
1861 /* high level multiplication (handles sign) */
1862 int mp_mul (mp_int * a, mp_int * b, mp_int * c)
1868 /* d = a * b (mod c) */
1869 int mp_mulmod (mp_int * a, mp_int * b, mp_int * c, mp_int * d)
1871 return fp_mulmod(a, b, c, d);
1874 /* c = a mod b, 0 <= c < b */
1875 int mp_mod (mp_int * a, mp_int * b, mp_int * c)
1877 return fp_mod (a, b, c);
1880 /* hac 14.61, pp608 */
1881 int mp_invmod (mp_int * a, mp_int * b, mp_int * c)
1883 return fp_invmod(a, b, c);
1886 /* this is a shell function that calls either the normal or Montgomery
1887 * exptmod functions. Originally the call to the montgomery code was
1888 * embedded in the normal function but that wasted alot of stack space
1889 * for nothing (since 99% of the time the Montgomery code would be called)
1891 int mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y)
1893 return fp_exptmod(G, X, P, Y);
1896 /* compare two ints (signed)*/
1897 int mp_cmp (mp_int * a, mp_int * b)
1899 return fp_cmp(a, b);
1902 /* compare a digit */
1903 int mp_cmp_d(mp_int * a, mp_digit b)
1905 return fp_cmp_d(a, b);
1908 /* get the size for an unsigned equivalent */
1909 int mp_unsigned_bin_size (mp_int * a)
1911 return fp_unsigned_bin_size(a);
1914 /* store in unsigned [big endian] format */
1915 int mp_to_unsigned_bin (mp_int * a, unsigned char *b)
1917 fp_to_unsigned_bin(a,b);
1921 /* reads a unsigned char array, assumes the msb is stored first [big endian] */
1922 int mp_read_unsigned_bin (mp_int * a, const unsigned char *b, int c)
1924 fp_read_unsigned_bin(a, (unsigned char *)b, c);
1929 int mp_sub_d(fp_int *a, fp_digit b, fp_int *c)
1936 /* fast math conversion */
1937 int mp_copy(fp_int* a, fp_int* b)
1944 /* fast math conversion */
1945 int mp_isodd(mp_int* a)
1951 /* fast math conversion */
1952 int mp_iszero(mp_int* a)
1954 return fp_iszero(a);
1958 /* fast math conversion */
1959 int mp_count_bits (mp_int* a)
1961 return fp_count_bits(a);
1965 /* fast math wrappers */
1966 int mp_set_int(fp_int *a, fp_digit b)
1973 #if defined(CYASSL_KEY_GEN) || defined (HAVE_ECC)
1975 /* c = a * a (mod b) */
1976 int fp_sqrmod(fp_int *a, fp_int *b, fp_int *c)
1981 return fp_mod(&tmp, b, c);
1984 /* fast math conversion */
1985 int mp_sqrmod(mp_int *a, mp_int *b, mp_int *c)
1987 return fp_sqrmod(a, b, c);
1990 /* fast math conversion */
1991 int mp_montgomery_calc_normalization(mp_int *a, mp_int *b)
1993 fp_montgomery_calc_normalization(a, b);
1997 #endif /* CYASSL_KEYGEN || HAVE_ECC */
2000 #ifdef CYASSL_KEY_GEN
2002 void fp_gcd(fp_int *a, fp_int *b, fp_int *c);
2003 void fp_lcm(fp_int *a, fp_int *b, fp_int *c);
2004 int fp_isprime(fp_int *a);
2005 int fp_cnt_lsb(fp_int *a);
2007 int mp_gcd(fp_int *a, fp_int *b, fp_int *c)
2014 int mp_lcm(fp_int *a, fp_int *b, fp_int *c)
2021 int mp_prime_is_prime(mp_int* a, int t, int* result)
2024 *result = fp_isprime(a);
2030 static int s_is_power_of_two(fp_digit b, int *p)
2034 /* fast return if no power of two */
2035 if ((b==0) || (b & (b-1))) {
2039 for (x = 0; x < DIGIT_BIT; x++) {
2040 if (b == (((fp_digit)1)<<x)) {
2048 /* a/b => cb + d == a */
2049 static int fp_div_d(fp_int *a, fp_digit b, fp_int *c, fp_digit *d)
2056 /* cannot divide by zero */
2062 if (b == 1 || fp_iszero(a) == 1) {
2072 /* power of two ? */
2073 if (s_is_power_of_two(b, &ix) == 1) {
2075 *d = a->dp[0] & ((((fp_digit)1)<<ix) - 1);
2078 fp_div_2d(a, ix, c, NULL);
2083 /* no easy answer [c'est la vie]. Just division */
2089 for (ix = a->used - 1; ix >= 0; ix--) {
2090 w = (w << ((fp_word)DIGIT_BIT)) | ((fp_word)a->dp[ix]);
2093 t = (fp_digit)(w / b);
2094 w -= ((fp_word)t) * ((fp_word)b);
2098 q.dp[ix] = (fp_digit)t;
2114 /* c = a mod b, 0 <= c < b */
2115 static int fp_mod_d(fp_int *a, fp_digit b, fp_digit *c)
2117 return fp_div_d(a, b, NULL, c);
2121 /* Miller-Rabin test of "a" to the base of "b" as described in
2122 * HAC pp. 139 Algorithm 4.24
2124 * Sets result to 0 if definitely composite or 1 if probably prime.
2125 * Randomly the chance of error is no more than 1/4 and often
2128 static void fp_prime_miller_rabin (fp_int * a, fp_int * b, int *result)
2137 if (fp_cmp_d(b, 1) != FP_GT) {
2141 /* get n1 = a - 1 */
2142 fp_init_copy(&n1, a);
2143 fp_sub_d(&n1, 1, &n1);
2145 /* set 2**s * r = n1 */
2146 fp_init_copy(&r, &n1);
2148 /* count the number of least significant bits
2153 /* now divide n - 1 by 2**s */
2154 fp_div_2d (&r, s, &r, NULL);
2156 /* compute y = b**r mod a */
2158 fp_exptmod(b, &r, a, &y);
2160 /* if y != 1 and y != n1 do */
2161 if (fp_cmp_d (&y, 1) != FP_EQ && fp_cmp (&y, &n1) != FP_EQ) {
2163 /* while j <= s-1 and y != n1 */
2164 while ((j <= (s - 1)) && fp_cmp (&y, &n1) != FP_EQ) {
2165 fp_sqrmod (&y, a, &y);
2167 /* if y == 1 then composite */
2168 if (fp_cmp_d (&y, 1) == FP_EQ) {
2174 /* if y != n1 then composite */
2175 if (fp_cmp (&y, &n1) != FP_EQ) {
2180 /* probably prime now */
2186 static const fp_digit primes[256] = {
2187 0x0002, 0x0003, 0x0005, 0x0007, 0x000B, 0x000D, 0x0011, 0x0013,
2188 0x0017, 0x001D, 0x001F, 0x0025, 0x0029, 0x002B, 0x002F, 0x0035,
2189 0x003B, 0x003D, 0x0043, 0x0047, 0x0049, 0x004F, 0x0053, 0x0059,
2190 0x0061, 0x0065, 0x0067, 0x006B, 0x006D, 0x0071, 0x007F, 0x0083,
2191 0x0089, 0x008B, 0x0095, 0x0097, 0x009D, 0x00A3, 0x00A7, 0x00AD,
2192 0x00B3, 0x00B5, 0x00BF, 0x00C1, 0x00C5, 0x00C7, 0x00D3, 0x00DF,
2193 0x00E3, 0x00E5, 0x00E9, 0x00EF, 0x00F1, 0x00FB, 0x0101, 0x0107,
2194 0x010D, 0x010F, 0x0115, 0x0119, 0x011B, 0x0125, 0x0133, 0x0137,
2196 0x0139, 0x013D, 0x014B, 0x0151, 0x015B, 0x015D, 0x0161, 0x0167,
2197 0x016F, 0x0175, 0x017B, 0x017F, 0x0185, 0x018D, 0x0191, 0x0199,
2198 0x01A3, 0x01A5, 0x01AF, 0x01B1, 0x01B7, 0x01BB, 0x01C1, 0x01C9,
2199 0x01CD, 0x01CF, 0x01D3, 0x01DF, 0x01E7, 0x01EB, 0x01F3, 0x01F7,
2200 0x01FD, 0x0209, 0x020B, 0x021D, 0x0223, 0x022D, 0x0233, 0x0239,
2201 0x023B, 0x0241, 0x024B, 0x0251, 0x0257, 0x0259, 0x025F, 0x0265,
2202 0x0269, 0x026B, 0x0277, 0x0281, 0x0283, 0x0287, 0x028D, 0x0293,
2203 0x0295, 0x02A1, 0x02A5, 0x02AB, 0x02B3, 0x02BD, 0x02C5, 0x02CF,
2205 0x02D7, 0x02DD, 0x02E3, 0x02E7, 0x02EF, 0x02F5, 0x02F9, 0x0301,
2206 0x0305, 0x0313, 0x031D, 0x0329, 0x032B, 0x0335, 0x0337, 0x033B,
2207 0x033D, 0x0347, 0x0355, 0x0359, 0x035B, 0x035F, 0x036D, 0x0371,
2208 0x0373, 0x0377, 0x038B, 0x038F, 0x0397, 0x03A1, 0x03A9, 0x03AD,
2209 0x03B3, 0x03B9, 0x03C7, 0x03CB, 0x03D1, 0x03D7, 0x03DF, 0x03E5,
2210 0x03F1, 0x03F5, 0x03FB, 0x03FD, 0x0407, 0x0409, 0x040F, 0x0419,
2211 0x041B, 0x0425, 0x0427, 0x042D, 0x043F, 0x0443, 0x0445, 0x0449,
2212 0x044F, 0x0455, 0x045D, 0x0463, 0x0469, 0x047F, 0x0481, 0x048B,
2214 0x0493, 0x049D, 0x04A3, 0x04A9, 0x04B1, 0x04BD, 0x04C1, 0x04C7,
2215 0x04CD, 0x04CF, 0x04D5, 0x04E1, 0x04EB, 0x04FD, 0x04FF, 0x0503,
2216 0x0509, 0x050B, 0x0511, 0x0515, 0x0517, 0x051B, 0x0527, 0x0529,
2217 0x052F, 0x0551, 0x0557, 0x055D, 0x0565, 0x0577, 0x0581, 0x058F,
2218 0x0593, 0x0595, 0x0599, 0x059F, 0x05A7, 0x05AB, 0x05AD, 0x05B3,
2219 0x05BF, 0x05C9, 0x05CB, 0x05CF, 0x05D1, 0x05D5, 0x05DB, 0x05E7,
2220 0x05F3, 0x05FB, 0x0607, 0x060D, 0x0611, 0x0617, 0x061F, 0x0623,
2221 0x062B, 0x062F, 0x063D, 0x0641, 0x0647, 0x0649, 0x064D, 0x0653
2224 int fp_isprime(fp_int *a)
2230 /* do trial division */
2231 for (r = 0; r < 256; r++) {
2232 fp_mod_d(a, primes[r], &d);
2238 /* now do 8 miller rabins */
2240 for (r = 0; r < 8; r++) {
2241 fp_set(&b, primes[r]);
2242 fp_prime_miller_rabin(a, &b, &res);
2252 void fp_lcm(fp_int *a, fp_int *b, fp_int *c)
2259 if (fp_cmp_mag(a, b) == FP_GT) {
2260 fp_div(a, &t1, &t2, NULL);
2263 fp_div(b, &t1, &t2, NULL);
2269 static const int lnz[16] = {
2270 4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0
2273 /* Counts the number of lsbs which are zero before the first zero bit */
2274 int fp_cnt_lsb(fp_int *a)
2280 if (fp_iszero(a) == 1) {
2284 /* scan lower digits until non-zero */
2285 for (x = 0; x < a->used && a->dp[x] == 0; x++);
2289 /* now scan this digit until a 1 is found */
2302 void fp_gcd(fp_int *a, fp_int *b, fp_int *c)
2306 /* either zero than gcd is the largest */
2307 if (fp_iszero (a) == 1 && fp_iszero (b) == 0) {
2311 if (fp_iszero (a) == 0 && fp_iszero (b) == 1) {
2316 /* optimized. At this point if a == 0 then
2317 * b must equal zero too
2319 if (fp_iszero (a) == 1) {
2325 if (fp_cmp_mag(a, b) != FP_LT) {
2326 fp_init_copy(&u, a);
2327 fp_init_copy(&v, b);
2329 fp_init_copy(&u, b);
2330 fp_init_copy(&v, a);
2334 while (fp_iszero(&v) == FP_NO) {
2342 #endif /* CYASSL_KEY_GEN */
2345 #if defined(HAVE_ECC) || !defined(NO_PWDBASED)
2347 void fp_add_d(fp_int *a, fp_digit b, fp_int *c)
2354 /* external compatibility */
2355 int mp_add_d(fp_int *a, fp_digit b, fp_int *c)
2361 #endif /* HAVE_ECC || !NO_PWDBASED */
2366 /* chars used in radix conversions */
2367 const char *fp_s_rmap = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz+/";
2369 static int fp_read_radix(fp_int *a, const char *str, int radix)
2374 /* make sure the radix is ok */
2375 if (radix < 2 || radix > 64) {
2379 /* if the leading digit is a
2380 * minus set the sign to negative.
2389 /* set the integer to the default of zero */
2392 /* process each digit of the string */
2394 /* if the radix < 36 the conversion is case insensitive
2395 * this allows numbers like 1AB and 1ab to represent the same value
2398 ch = (char) ((radix < 36) ? XTOUPPER(*str) : *str);
2399 for (y = 0; y < 64; y++) {
2400 if (ch == fp_s_rmap[y]) {
2405 /* if the char was found in the map
2406 * and is less than the given radix add it
2407 * to the number, otherwise exit the loop.
2410 fp_mul_d (a, (fp_digit) radix, a);
2411 fp_add_d (a, (fp_digit) y, a);
2418 /* set the sign only if a != 0 */
2419 if (fp_iszero(a) != FP_YES) {
2425 /* fast math conversion */
2426 int mp_read_radix(mp_int *a, const char *str, int radix)
2428 return fp_read_radix(a, str, radix);
2431 /* fast math conversion */
2432 int mp_set(fp_int *a, fp_digit b)
2438 /* fast math conversion */
2439 int mp_sqr(fp_int *A, fp_int *B)
2445 /* fast math conversion */
2446 int mp_montgomery_reduce(fp_int *a, fp_int *m, fp_digit mp)
2448 fp_montgomery_reduce(a, m, mp);
2453 /* fast math conversion */
2454 int mp_montgomery_setup(fp_int *a, fp_digit *rho)
2456 return fp_montgomery_setup(a, rho);
2459 int mp_div_2(fp_int * a, fp_int * b)
2466 int mp_init_copy(fp_int * a, fp_int * b)
2474 #endif /* HAVE_ECC */
2476 #endif /* USE_FAST_MATH */