3 This file is part of the ARM-Crypto-Lib.
4 Copyright (C) 2008 Daniel Otte (daniel.otte@rub.de)
6 This program is free software: you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation, either version 3 of the License, or
9 (at your option) any later version.
11 This program is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
16 You should have received a copy of the GNU General Public License
17 along with this program. If not, see <http://www.gnu.org/licenses/>.
24 * \license GPLv3 or later
30 #define STRING(x) STRING2(x)
31 #define STR_LINE STRING(__LINE__)
41 #include "bigint_io.h"
46 #define MAX(a,b) (((a)>(b))?(a):(b))
50 #define MIN(a,b) (((a) < (b)) ? (a) : (b))
53 #define SET_FBS(a, v) do {(a)->info &= ~BIGINT_FBS_MASK; (a)->info |= (v);} while(0)
54 #define GET_FBS(a) ((a)->info & BIGINT_FBS_MASK)
55 #define SET_NEG(a) (a)->info |= BIGINT_NEG_MASK
56 #define SET_POS(a) (a)->info &= ~BIGINT_NEG_MASK
57 #define XCHG(a,b) do{(a) ^= (b); (b) ^= (a); (a) ^= (b);} while(0)
58 #define XCHG_PTR(a,b) do{ a = (void*)(((intptr_t)(a)) ^ ((intptr_t)(b))); \
59 b = (void*)(((intptr_t)(a)) ^ ((intptr_t)(b))); \
60 a = (void*)(((intptr_t)(a)) ^ ((intptr_t)(b)));} while(0)
62 #define GET_SIGN(a) ((a)->info & BIGINT_NEG_MASK)
64 /******************************************************************************/
65 void bigint_adjust(bigint_t *a){
66 while (a->length_W != 0 && a->wordv[a->length_W - 1] == 0) {
69 if (a->length_W == 0) {
74 uint8_t i = BIGINT_WORD_SIZE - 1;
75 t = a->wordv[a->length_W - 1];
76 while ((t & (((bigint_word_t)1) << (BIGINT_WORD_SIZE - 1))) == 0 && i) {
83 /******************************************************************************/
85 bigint_length_t bigint_length_b(const bigint_t *a){
86 if(!a->length_W || a->length_W==0){
89 return (a->length_W-1) * BIGINT_WORD_SIZE + GET_FBS(a);
92 /******************************************************************************/
94 bigint_length_t bigint_length_B(const bigint_t *a){
95 return a->length_W * sizeof(bigint_word_t);
98 /******************************************************************************/
100 int32_t bigint_get_first_set_bit(const bigint_t *a){
101 if(a->length_W == 0) {
104 return (a->length_W-1) * sizeof(bigint_word_t) * CHAR_BIT + GET_FBS(a);
108 /******************************************************************************/
110 int32_t bigint_get_last_set_bit(const bigint_t *a){
114 if (a->length_W == 0) {
117 while (a->wordv[r] == 0 && r < a->length_W) {
120 if (a->wordv[r] == 0) {
121 return (uint32_t)(-1);
123 while ((x&a->wordv[r])==0) {
127 return r * BIGINT_WORD_SIZE + b;
130 /******************************************************************************/
132 void bigint_copy(bigint_t *dest, const bigint_t *src){
133 if(dest->wordv != src->wordv){
134 memcpy(dest->wordv, src->wordv, src->length_W * sizeof(bigint_word_t));
136 dest->length_W = src->length_W;
137 dest->info = src->info;
140 /******************************************************************************/
142 /* this should be implemented in assembly */
143 void bigint_add_u(bigint_t *dest, const bigint_t *a, const bigint_t *b){
145 bigint_wordplus_t t = 0LL;
146 if (a->length_W < b->length_W) {
149 for(i = 0; i < b->length_W; ++i) {
152 dest->wordv[i] = (bigint_word_t)t;
153 t >>= BIGINT_WORD_SIZE;
155 for(; i < a->length_W; ++i){
157 dest->wordv[i] = (bigint_word_t)t;
158 t >>= BIGINT_WORD_SIZE;
161 dest->wordv[i++] = (bigint_word_t)t;
168 /******************************************************************************/
170 /* this should be implemented in assembly */
171 void bigint_add_scale_u(bigint_t *dest, const bigint_t *a, bigint_length_t scale){
172 if(a->length_W == 0){
176 bigint_add_u(dest, dest, a);
180 #if BIGINT_WORD_SIZE == 8
181 memset(dest->wordv + dest->length_W, 0, MAX(dest->length_W, a->length_W + scale) - dest->length_W);
182 x.wordv = dest->wordv + scale;
183 x.length_W = dest->length_W - scale;
184 if((int16_t)x.length_W < 0){
190 bigint_add_u(&x, &x, a);
191 dest->length_W = x.length_W + scale;
196 bigint_length_t word_shift = scale / sizeof(bigint_word_t), byte_shift = scale % sizeof(bigint_word_t);
197 bigint_word_t bv[a->length_W + 1];
199 bv[0] = bv[a->length_W] = 0;
200 memcpy((uint8_t*)bv + byte_shift, a->wordv, a->length_W * sizeof(bigint_word_t));
201 s.length_W = a->length_W + 1;
203 memset(dest->wordv + dest->length_W, 0, (MAX(dest->length_W, s.length_W + word_shift) - dest->length_W) * sizeof(bigint_word_t));
204 x.wordv = dest->wordv + word_shift;
205 x.length_W = dest->length_W - word_shift;
206 if((int16_t)x.length_W < 0){
212 bigint_add_u(&x, &x, &s);
213 dest->length_W = x.length_W + word_shift;
219 /******************************************************************************/
221 /* this should be implemented in assembly */
222 void bigint_sub_u(bigint_t *dest, const bigint_t *a, const bigint_t *b){
225 bigint_wordplus_signed_t t = 0;
227 if(b->length_W == 0){
228 bigint_copy(dest, a);
232 if(a->length_W == 0){
233 bigint_copy(dest, b);
237 r = bigint_cmp_u(a,b);
239 bigint_set_zero(dest);
243 bigint_sub_u(dest, b, a);
247 for(i = 0; i < a->length_W; ++i){
253 dest->wordv[i] = (bigint_word_t)t;
254 borrow = t < 0 ? 1 : 0;
261 /******************************************************************************/
263 int8_t bigint_cmp_u(const bigint_t *a, const bigint_t *b){
264 if(a->length_W > b->length_W){
267 if(a->length_W < b->length_W){
270 if(a->length_W == 0){
276 if(a->wordv[i] != b->wordv[i]){
277 if(a->wordv[i] > b->wordv[i]){
287 /******************************************************************************/
289 void bigint_add_s(bigint_t *dest, const bigint_t *a, const bigint_t *b){
292 s |= GET_SIGN(b)?1:0;
294 case 0: /* both positive */
295 bigint_add_u(dest, a,b);
298 case 1: /* a positive, b negative */
299 bigint_sub_u(dest, a, b);
301 case 2: /* a negative, b positive */
302 bigint_sub_u(dest, b, a);
304 case 3: /* both negative */
305 bigint_add_u(dest, a, b);
308 default: /* how can this happen?*/
313 /******************************************************************************/
315 void bigint_sub_s(bigint_t *dest, const bigint_t *a, const bigint_t *b){
318 s |= GET_SIGN(b)?1:0;
320 case 0: /* both positive */
321 bigint_sub_u(dest, a,b);
323 case 1: /* a positive, b negative */
324 bigint_add_u(dest, a, b);
327 case 2: /* a negative, b positive */
328 bigint_add_u(dest, a, b);
331 case 3: /* both negative */
332 bigint_sub_u(dest, b, a);
334 default: /* how can this happen?*/
340 /******************************************************************************/
342 int8_t bigint_cmp_s(const bigint_t *a, const bigint_t *b){
344 if(a->length_W==0 && b->length_W==0){
348 s |= GET_SIGN(b)?1:0;
350 case 0: /* both positive */
351 return bigint_cmp_u(a, b);
353 case 1: /* a positive, b negative */
356 case 2: /* a negative, b positive */
359 case 3: /* both negative */
360 return bigint_cmp_u(b, a);
362 default: /* how can this happen?*/
365 return 0; /* just to satisfy the compiler */
368 /******************************************************************************/
370 void bigint_shiftleft(bigint_t *a, bigint_length_t shift){
371 bigint_length_t byteshift, words_to_shift;
375 bigint_wordplus_t t = 0;
377 if (a->length_W == 0 || shift == 0) {
380 byteshift = shift / 8;
381 bitshift = shift & 7;
382 if (byteshift % sizeof(bigint_word_t)) {
383 a->wordv[a->length_W + byteshift / sizeof(bigint_t)] = 0;
386 memmove(((uint8_t*)a->wordv) + byteshift, a->wordv, a->length_W * sizeof(bigint_word_t));
387 memset(a->wordv, 0, byteshift);
390 a->length_W += (byteshift + sizeof(bigint_word_t) - 1) / sizeof(bigint_word_t);
394 p = &a->wordv[byteshift / sizeof(bigint_word_t)];
395 words_to_shift = a->length_W + (byteshift % sizeof(bigint_word_t) ? 1 : 0);
396 for (i = 0; i < words_to_shift; ++i) {
397 t |= ((bigint_wordplus_t)p[i]) << bitshift;
398 p[i] = (bigint_word_t)t;
399 t >>= BIGINT_WORD_SIZE;
402 p[i] = (bigint_word_t)t;
405 a->length_W += (byteshift + sizeof(bigint_word_t) - 1) / sizeof(bigint_word_t);
409 /******************************************************************************/
411 void bigint_shiftright(bigint_t *a, bigint_length_t shift){
412 bigint_length_t byteshift;
415 bigint_wordplus_t t = 0;
416 byteshift = shift / 8;
417 bitshift = shift & 7;
419 if (a->length_W == 0) {
423 if(bigint_get_first_set_bit(a) < shift){ /* we would shift out more than we have */
429 memmove(a->wordv, (uint8_t*)a->wordv + byteshift, a->length_W * sizeof(bigint_word_t) - byteshift);
430 memset((uint8_t*)&a->wordv[a->length_W] - byteshift, 0, byteshift);
431 a->length_W -= byteshift / sizeof(bigint_word_t);
435 if(bitshift != 0 && a->length_W){
436 /* shift to the right */
439 t |= ((bigint_wordplus_t)(a->wordv[i])) << (BIGINT_WORD_SIZE - bitshift);
440 a->wordv[i] = (bigint_word_t)(t >> BIGINT_WORD_SIZE);
441 t <<= BIGINT_WORD_SIZE;
447 /******************************************************************************/
449 void bigint_xor(bigint_t *dest, const bigint_t *a){
451 for(i=0; i<a->length_W; ++i){
452 dest->wordv[i] ^= a->wordv[i];
457 /******************************************************************************/
459 void bigint_set_zero(bigint_t *a){
463 /******************************************************************************/
465 /* using the Karatsuba-Algorithm */
466 /* x*y = (xh*yh)*b**2n + ((xh+xl)*(yh+yl) - xh*yh - xl*yl)*b**n + yh*yl */
467 void bigint_mul_u(bigint_t *dest, const bigint_t *a, const bigint_t *b){
468 if(a->length_W == 0 || b->length_W == 0){
469 bigint_set_zero(dest);
472 if(dest == a || dest == b){
474 bigint_word_t d_b[a->length_W + b->length_W];
476 bigint_mul_u(&d, a, b);
477 bigint_copy(dest, &d);
480 if(a->length_W == 1 || b->length_W == 1){
481 if(a->length_W != 1){
484 bigint_wordplus_t t = 0;
486 bigint_word_t x = a->wordv[0];
487 for(i=0; i < b->length_W; ++i){
488 t += ((bigint_wordplus_t)b->wordv[i]) * ((bigint_wordplus_t)x);
489 dest->wordv[i] = (bigint_word_t)t;
490 t >>= BIGINT_WORD_SIZE;
494 dest->wordv[i] = (bigint_word_t)t;
501 if(a->length_W * sizeof(bigint_word_t) <= 4 && b->length_W * sizeof(bigint_word_t) <= 4){
504 memcpy(&p, a->wordv, a->length_W*sizeof(bigint_word_t));
505 memcpy(&q, b->wordv, b->length_W*sizeof(bigint_word_t));
506 r = (uint64_t)p * (uint64_t)q;
507 memcpy(dest->wordv, &r, (dest->length_W = a->length_W + b->length_W)*sizeof(bigint_word_t));
511 /* split a in xh & xl; split b in yh & yl */
512 const bigint_length_t n = (MAX(a->length_W, b->length_W)+1)/2;
513 bigint_t xl, xh, yl, yh;
517 bigint_set_zero(&xh);
518 xl.length_W = a->length_W;
524 xh.wordv = &(a->wordv[n]);
525 xh.length_W = a->length_W-n;
529 bigint_set_zero(&yh);
530 yl.length_W = b->length_W;
536 yh.wordv = &(b->wordv[n]);
537 yh.length_W = b->length_W-n;
540 /* now we have split up a and b */
541 /* remember we want to do:
542 * x*y = (xh * b**n + xl) * (yh * b**n + yl)
543 * = (xh * yh) * b**2n + xh * b**n * yl + yh * b**n * xl + xl * yl
544 * = (xh * yh) * b**2n + (xh * yl + yh * xl) * b**n + xl *yl
545 * // xh * yl + yh * xl = (xh + yh) * (xl + yl) - xh * yh - xl * yl
546 * x*y = (xh * yh) * b**2n + ((xh+xl)*(yh+yl) - xh*yh - xl*yl)*b**n + xl*yl
547 * 5 9 2 4 3 7 5 6 1 8 1
549 bigint_word_t tmp_b[2 * n + 2], m_b[2 * (n + 1)];
550 bigint_t tmp, tmp2, m;
552 tmp2.wordv = &(tmp_b[n + 1]);
555 bigint_mul_u(dest, &xl, &yl); /* 1: dest <= xl*yl */
556 bigint_add_u(&tmp2, &xh, &xl); /* 2: tmp2 <= xh+xl */
557 bigint_add_u(&tmp, &yh, &yl); /* 3: tmp <= yh+yl */
558 bigint_mul_u(&m, &tmp2, &tmp); /* 4: m <= tmp2*tmp */
559 bigint_mul_u(&tmp, &xh, &yh); /* 5: tmp <= xh*yh */
560 bigint_sub_u(&m, &m, dest); /* 6: m <= m-dest */
561 bigint_sub_u(&m, &m, &tmp); /* 7: m <= m-tmp */
562 bigint_add_scale_u(dest, &m, n * sizeof(bigint_word_t)); /* 8: dest <= dest+m**n*/
563 bigint_add_scale_u(dest, &tmp, 2 * n * sizeof(bigint_word_t)); /* 9: dest <= dest+tmp**(2*n) */
566 /******************************************************************************/
568 void bigint_mul_s(bigint_t *dest, const bigint_t *a, const bigint_t *b){
571 s |= GET_SIGN(b)?1:0;
573 case 0: /* both positive */
574 bigint_mul_u(dest, a,b);
577 case 1: /* a positive, b negative */
578 bigint_mul_u(dest, a,b);
581 case 2: /* a negative, b positive */
582 bigint_mul_u(dest, a,b);
585 case 3: /* both negative */
586 bigint_mul_u(dest, a,b);
589 default: /* how can this happen?*/
594 /******************************************************************************/
598 unsigned square_depth = 0;
602 /* (xh*b^n+xl)^2 = xh^2*b^2n + 2*xh*xl*b^n + xl^2 */
603 void bigint_square(bigint_t *dest, const bigint_t *a){
604 if(a->length_W * sizeof(bigint_word_t) <= 4){
606 memcpy(&r, a->wordv, a->length_W * sizeof(bigint_word_t));
608 memcpy(dest->wordv, &r, 2 * a->length_W * sizeof(bigint_word_t));
610 dest->length_W = 2 * a->length_W;
615 if(dest->wordv == a->wordv){
617 bigint_word_t d_b[a->length_W*2];
619 bigint_square(&d, a);
620 bigint_copy(dest, &d);
623 bigint_fast_square(dest, a);
631 bigint_t xh, xl, tmp; /* x-high, x-low, temp */
632 bigint_word_t buffer[2*n+1];
636 xh.wordv = &(a->wordv[n]);
637 xh.length_W = a->length_W-n;
641 /* (xh * b**n + xl)**2 = xh**2 * b**2n + 2 * xh * xl * b**n + xl**2 */
643 if(square_depth == 1){
644 cli_putstr("\r\nDBG (a): xl: "); bigint_print_hex(&xl);
645 cli_putstr("\r\nDBG (b): xh: "); bigint_print_hex(&xh);
648 bigint_square(dest, &xl);
650 if(square_depth == 1){
651 cli_putstr("\r\nDBG (1): xl**2: "); bigint_print_hex(dest);
654 bigint_square(&tmp, &xh);
656 if(square_depth == 1){
657 cli_putstr("\r\nDBG (2): xh**2: "); bigint_print_hex(&tmp);
660 bigint_add_scale_u(dest, &tmp, 2 * n * sizeof(bigint_word_t));
662 if(square_depth == 1){
663 cli_putstr("\r\nDBG (3): xl**2 + xh**2*n**2: "); bigint_print_hex(dest);
666 bigint_mul_u(&tmp, &xl, &xh);
668 if(square_depth == 1){
669 cli_putstr("\r\nDBG (4): xl*xh: "); bigint_print_hex(&tmp);
672 bigint_shiftleft(&tmp, 1);
674 if(square_depth == 1){
675 cli_putstr("\r\nDBG (5): xl*xh*2: "); bigint_print_hex(&tmp);
678 bigint_add_scale_u(dest, &tmp, n * sizeof(bigint_word_t));
680 if(square_depth == 1){
681 cli_putstr("\r\nDBG (6): x**2: "); bigint_print_hex(dest);
690 /******************************************************************************/
692 void bigint_square(bigint_t *dest, const bigint_t *a) {
693 union __attribute__((packed)) {
695 bigint_wordplus_t uv;
697 bigint_word_t q, c1, c2;
698 bigint_length_t i, j;
700 if(a->length_W * sizeof(bigint_word_t) <= 4){
702 memcpy(&r, a->wordv, a->length_W * sizeof(bigint_word_t));
704 memcpy(dest->wordv, &r, 2 * a->length_W * sizeof(bigint_word_t));
706 dest->length_W = 2 * a->length_W;
711 if(dest->wordv == a->wordv){
713 bigint_word_t d_b[a->length_W * 2];
715 bigint_square(&d, a);
716 bigint_copy(dest, &d);
720 memset(dest->wordv, 0, a->length_W * 2 * sizeof(bigint_word_t));
722 for (i = 0; i < a->length_W; ++i) {
723 acc.uv = (bigint_wordplus_t)a->wordv[i] * (bigint_wordplus_t)a->wordv[i] + (bigint_wordplus_t)dest->wordv[2 * i];
724 dest->wordv[2 * i] = acc.u[0];
727 for (j = i + 1; j < a->length_W; ++j) {
728 acc.uv = (bigint_wordplus_t)a->wordv[i] * (bigint_wordplus_t)a->wordv[j];
729 q = acc.u[1] >> (BIGINT_WORD_SIZE - 1);
731 acc.uv += dest->wordv[i + j];
732 q += (acc.uv < dest->wordv[i + j]);
735 dest->wordv[i + j] = acc.u[0];
736 c1 = (bigint_wordplus_t)acc.u[1] + c2;
739 dest->wordv[i + a->length_W] += c1;
740 if (i < a->length_W - 1) {
741 dest->wordv[i + a->length_W + 1] = c2 + (dest->wordv[i + a->length_W] < c1);
745 dest->length_W = 2 * a->length_W;
749 /******************************************************************************/
751 void bigint_sub_u_bitscale(bigint_t *a, const bigint_t *b, bigint_length_t bitscale){
753 bigint_word_t tmp_b[b->length_W + 1];
754 const bigint_length_t word_shift = bitscale / BIGINT_WORD_SIZE;
756 if (a->length_W < b->length_W + word_shift) {
758 cli_putstr("\r\nDBG: *bang*\r\n");
764 bigint_copy(&tmp, b);
765 bigint_shiftleft(&tmp, bitscale % BIGINT_WORD_SIZE);
768 x.wordv = &(a->wordv[word_shift]);
769 x.length_W = a->length_W - word_shift;
771 bigint_sub_u(&x, &x, &tmp);
776 /******************************************************************************/
778 void bigint_reduce(bigint_t *a, const bigint_t *r){
779 uint8_t rfbs = GET_FBS(r);
780 if (r->length_W == 0 || a->length_W == 0) {
784 if (bigint_length_b(a) + 3 > bigint_length_b(r)) {
785 if ((r->length_W * sizeof(bigint_word_t) <= 4) && (a->length_W * sizeof(bigint_word_t) <= 4)) {
786 uint32_t p = 0, q = 0;
787 memcpy(&p, a->wordv, a->length_W * sizeof(bigint_word_t));
788 memcpy(&q, r->wordv, r->length_W * sizeof(bigint_word_t));
790 memcpy(a->wordv, &p, a->length_W * sizeof(bigint_word_t));
791 a->length_W = r->length_W;
796 while (a->length_W > r->length_W) {
797 shift = (a->length_W - r->length_W) * CHAR_BIT * sizeof(bigint_word_t) + GET_FBS(a) - rfbs - 1;
798 bigint_sub_u_bitscale(a, r, shift);
800 while ((GET_FBS(a) > rfbs) && (a->length_W == r->length_W)) {
801 shift = GET_FBS(a) - rfbs - 1;
802 bigint_sub_u_bitscale(a, r, shift);
805 while (bigint_cmp_u(a, r) >= 0) {
806 bigint_sub_u(a, a, r);
811 /******************************************************************************/
813 /* calculate dest = a**exp % r */
814 /* using square&multiply */
815 void bigint_expmod_u_sam(bigint_t *dest, const bigint_t *a, const bigint_t *exp, const bigint_t *r){
816 if(a->length_W == 0){
817 bigint_set_zero(dest);
822 bigint_word_t t, base_w[MAX(a->length_W, r->length_W)], res_w[r->length_W * 2];
827 bigint_copy(&base, a);
828 bigint_reduce(&base, r);
833 if(exp->length_W == 0){
834 bigint_copy(dest, &res);
837 bigint_copy(&res, &base);
839 t = exp->wordv[exp->length_W - 1];
840 for (i = exp->length_W; i > 0; --i) {
841 t = exp->wordv[i - 1];
842 for (j = BIGINT_WORD_SIZE; j > 0; --j) {
844 if (t & (((bigint_word_t)1) << (BIGINT_WORD_SIZE - 1))) {
848 bigint_square(&res, &res);
849 bigint_reduce(&res, r);
850 if(t & (((bigint_word_t)1) << (BIGINT_WORD_SIZE - 1))){
851 bigint_mul_u(&res, &res, &base);
852 bigint_reduce(&res, r);
860 bigint_copy(dest, &res);
863 /******************************************************************************/
864 /* gcd <-- gcd(x,y) a*x+b*y=gcd */
865 void bigint_gcdext(bigint_t *gcd, bigint_t *a, bigint_t *b, const bigint_t *x, const bigint_t *y){
866 bigint_length_t i = 0;
867 if(x->length_W == 0 || y->length_W == 0){
869 bigint_set_zero(gcd);
879 if(x->length_W == 1 && x->wordv[0] == 1){
896 if(y->length_W == 1 && y->wordv[0] == 1){
914 while(x->wordv[i] == 0 && y->wordv[i] == 0){
917 bigint_word_t g_b[i + 2], x_b[x->length_W - i], y_b[y->length_W - i];
918 bigint_word_t u_b[x->length_W - i], v_b[y->length_W - i];
919 bigint_word_t a_b[y->length_W + 2], c_b[y->length_W + 2];
920 bigint_word_t b_b[x->length_W + 2], d_b[x->length_W + 2];
921 bigint_t g, x_, y_, u, v, a_, b_, c_, d_;
926 memset(g_b, 0, i * sizeof(bigint_word_t));
930 x_.info = y_.info = 0;
931 x_.length_W = x->length_W - i;
932 y_.length_W = y->length_W - i;
933 memcpy(x_.wordv, x->wordv + i, x_.length_W * sizeof(bigint_word_t));
934 memcpy(y_.wordv, y->wordv + i, y_.length_W * sizeof(bigint_word_t));
935 for(i = 0; (x_.wordv[0] & (1 << i)) == 0 && (y_.wordv[0] & (1 << i)) == 0; ++i){
942 bigint_shiftleft(&g, i);
943 bigint_shiftright(&x_, i);
944 bigint_shiftright(&y_, i);
954 bigint_copy(&u, &x_);
955 bigint_copy(&v, &y_);
962 bigint_set_zero(&b_);
963 bigint_set_zero(&c_);
965 while ((u.wordv[0] & 1) == 0) {
966 bigint_shiftright(&u, 1);
967 if((a_.wordv[0] & 1) || (b_.wordv[0] & 1)){
968 bigint_add_s(&a_, &a_, &y_);
969 bigint_sub_s(&b_, &b_, &x_);
971 bigint_shiftright(&a_, 1);
972 bigint_shiftright(&b_, 1);
974 while ((v.wordv[0] & 1) == 0) {
975 bigint_shiftright(&v, 1);
976 if((c_.wordv[0] & 1) || (d_.wordv[0] & 1)){
977 bigint_add_s(&c_, &c_, &y_);
978 bigint_sub_s(&d_, &d_, &x_);
980 bigint_shiftright(&c_, 1);
981 bigint_shiftright(&d_, 1);
983 if(bigint_cmp_u(&u, &v) >= 0){
984 bigint_sub_u(&u, &u, &v);
985 bigint_sub_s(&a_, &a_, &c_);
986 bigint_sub_s(&b_, &b_, &d_);
988 bigint_sub_u(&v, &v, &u);
989 bigint_sub_s(&c_, &c_, &a_);
990 bigint_sub_s(&d_, &d_, &b_);
994 bigint_mul_s(gcd, &v, &g);
1000 bigint_copy(b, &d_);
1004 /******************************************************************************/
1006 void bigint_inverse(bigint_t *dest, const bigint_t *a, const bigint_t *m){
1007 bigint_gcdext(NULL, dest, NULL, a, m);
1008 while(dest->info&BIGINT_NEG_MASK){
1009 bigint_add_s(dest, dest, m);
1013 /******************************************************************************/
1015 void bigint_changeendianess(bigint_t *a){
1017 p = (uint8_t*)(a->wordv);
1018 q = p + a->length_W * sizeof(bigint_word_t) - 1;
1027 /******************************************************************************/
1029 void bigint_mul_word_u(bigint_t *a, bigint_word_t b){
1030 bigint_wordplus_t c0 = 0, c1 = 0;
1038 for(i = 0; i < a->length_W; ++i){
1039 c1 = ((bigint_wordplus_t)(a->wordv[i])) * ((bigint_wordplus_t)b);
1041 a->wordv[i] = (bigint_word_t)c1;
1042 c0 = c1 >> BIGINT_WORD_SIZE;
1045 a->wordv[a->length_W] = (bigint_word_t)c0;
1051 /******************************************************************************/
1054 void bigint_clip(bigint_t *dest, bigint_length_t length_W){
1055 if(dest->length_W > length_W){
1056 dest->length_W = length_W;
1058 bigint_adjust(dest);
1060 /******************************************************************************/
1064 * dest = (a * b) % m (?)
1067 void bigint_mont_mul(bigint_t *dest, const bigint_t *a, const bigint_t *b, const bigint_t *m, const bigint_t *m_){
1068 const bigint_length_t s = MAX(MAX(a->length_W, b->length_W), m->length_W);
1070 bigint_word_t u_w[s + 2], t_w[s + 2];
1073 if (a->length_W == 0 || b->length_W == 0) {
1074 bigint_set_zero(dest);
1081 for (i = 0; i < a->length_W; ++i) {
1083 bigint_mul_word_u(&t, a->wordv[i]);
1084 bigint_add_u(&u, &u, &t);
1085 bigint_copy(&t, m_);
1086 if (u.length_W != 0) {
1087 bigint_mul_word_u(&t, u.wordv[0]);
1088 bigint_add_u(&u, &u, &t);
1090 bigint_shiftright(&u, BIGINT_WORD_SIZE);
1092 for (; i < s; ++i) {
1093 bigint_copy(&t, m_);
1094 if (u.length_W != 0) {
1095 bigint_mul_word_u(&t, u.wordv[0]);
1096 bigint_add_u(&u, &u, &t);
1098 bigint_shiftright(&u, BIGINT_WORD_SIZE);
1100 bigint_reduce(&u, m);
1101 bigint_copy(dest, &u);
1104 /******************************************************************************/
1106 void bigint_mont_red(bigint_t *dest, const bigint_t *a, const bigint_t *m, const bigint_t *m_){
1108 bigint_length_t i, s = MAX(a->length_W, m->length_W);
1109 bigint_word_t u_w[s + 2], t_w[s + 2];
1113 if (a->length_W == 0) {
1114 bigint_set_zero(dest);
1118 for (i = 0; i < m->length_W; ++i) {
1119 bigint_copy(&t, m_);
1120 if (u.length_W != 0) {
1121 bigint_mul_word_u(&t, u.wordv[0]);
1122 bigint_add_u(&u, &u, &t);
1124 bigint_shiftright(&u, BIGINT_WORD_SIZE);
1126 bigint_reduce(&u, m);
1127 bigint_copy(dest, &u);
1130 /******************************************************************************/
1132 * m_ = m * (- m0^-1 (mod 2^W))
1134 void bigint_mont_gen_m_(bigint_t* dest, const bigint_t* m){
1135 bigint_word_t x_w[2], m_w_0[1];
1139 bigint_print_hex(m);
1142 if (m->length_W == 0) {
1143 bigint_set_zero(dest);
1146 if ((m->wordv[0] & 1) == 0) {
1147 printf_P(PSTR("ERROR: m must not be even, m = "));
1148 bigint_print_hex(m);
1161 m_0.wordv[0] = m->wordv[0];
1163 bigint_adjust(&m_0);
1166 bigint_print_hex(&m_0);
1168 bigint_print_hex(&x);
1171 bigint_inverse(dest, &m_0, &x);
1173 printf("\nm0^-1 = ");
1174 bigint_print_hex(dest);
1177 bigint_sub_s(&x, &x, dest);
1179 printf("\n-m0^-1 = ");
1180 bigint_print_hex(&x);
1183 bigint_copy(dest, m);
1184 bigint_mul_word_u(dest, x.wordv[0]);
1187 /******************************************************************************/
1190 * dest = a * R mod m
1192 void bigint_mont_trans(bigint_t *dest, const bigint_t *a, const bigint_t *m){
1194 bigint_word_t t_w[a->length_W + m->length_W];
1197 memset(t_w, 0, m->length_W * sizeof(bigint_word_t));
1198 memcpy(&t_w[m->length_W], a->wordv, a->length_W * sizeof(bigint_word_t));
1200 t.length_W = a->length_W + m->length_W;
1201 bigint_reduce(&t, m);
1202 bigint_copy(dest, &t);
1205 /******************************************************************************/
1207 /* calculate dest = a**exp % r */
1208 /* using square&multiply */
1209 void bigint_expmod_u_mont_accel(bigint_t *dest, const bigint_t *a, const bigint_t *exp, const bigint_t *r, const bigint_t *m_){
1210 if(r->length_W == 0) {
1214 bigint_length_t s = r->length_W;
1216 bigint_word_t t, res_w[r->length_W * 2], ax_w[MAX(s, a->length_W)];
1226 bigint_adjust(&res);
1227 if (exp->length_W == 0) {
1228 bigint_copy(dest, &res);
1231 bigint_copy(&ax, a);
1232 bigint_reduce(&ax, r);
1233 bigint_mont_trans(&ax, &ax, r);
1234 bigint_mont_trans(&res, &res, r);
1236 t = exp->wordv[exp->length_W - 1];
1237 for (i = exp->length_W; i > 0; --i) {
1238 t = exp->wordv[i - 1];
1239 for(j = BIGINT_WORD_SIZE; j > 0; --j){
1241 if(t & (((bigint_wordplus_t)1) << (BIGINT_WORD_SIZE - 1))){
1246 bigint_square(&res, &res);
1247 bigint_mont_red(&res, &res, r, m_);
1248 if (t & (((bigint_wordplus_t)1) << (BIGINT_WORD_SIZE - 1))) {
1249 bigint_mont_mul(&res, &res, &ax, r, m_);
1256 bigint_mont_red(dest, &res, r, m_);
1259 /******************************************************************************/
1261 void bigint_expmod_u_mont_sam(bigint_t *dest, const bigint_t *a, const bigint_t *exp, const bigint_t *r){
1262 if(r->length_W == 0) {
1265 if(a->length_W == 0) {
1266 bigint_set_zero(dest);
1270 bigint_word_t m_w_[r->length_W + 1];
1272 bigint_mont_gen_m_(&m_, r);
1273 bigint_expmod_u_mont_accel(dest, a, exp, r,&m_);
1276 /******************************************************************************/
1280 void bigint_expmod_u(bigint_t *dest, const bigint_t *a, const bigint_t *exp, const bigint_t *r){
1281 if (r->wordv[0] & 1) {
1282 bigint_expmod_u_mont_sam(dest, a, exp, r);
1284 bigint_expmod_u_sam(dest, a, exp, r);