]> git.cryptolib.org Git - arm-crypto-lib.git/blob - bigint/bigint.c
switching to packed structure
[arm-crypto-lib.git] / bigint / bigint.c
1 /* bigint.c */
2 /*
3     This file is part of the ARM-Crypto-Lib.
4     Copyright (C) 2008  Daniel Otte (daniel.otte@rub.de)
5
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.
10
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.
15
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/>.
18 */
19 /**
20  * \file                bigint.c
21  * \author              Daniel Otte
22  * \date                2010-02-22
23  * 
24  * \license         GPLv3 or later
25  * 
26  */
27  
28
29 #define STRING2(x) #x
30 #define STRING(x) STRING2(x)
31 #define STR_LINE STRING(__LINE__)
32
33 #include "bigint.h"
34 #include <string.h>
35
36 #define DEBUG 0
37
38 #if DEBUG
39 #include "cli.h"
40 #include "bigint_io.h"
41 #endif
42
43 #ifndef MAX
44  #define MAX(a,b) (((a)>(b))?(a):(b))
45 #endif
46
47 #ifndef MIN
48  #define MIN(a,b) (((a)<(b))?(a):(b))
49 #endif
50
51 #define SET_FBS(a, v) do{(a)->info &=~BIGINT_FBS_MASK; (a)->info |= (v);}while(0)
52 #define GET_FBS(a)   ((a)->info&BIGINT_FBS_MASK)
53 #define SET_NEG(a)   (a)->info |= BIGINT_NEG_MASK
54 #define SET_POS(a)   (a)->info &= ~BIGINT_NEG_MASK
55 #define XCHG(a,b)    do{(a)^=(b); (b)^=(a); (a)^=(b);}while(0)
56 #define XCHG_PTR(a,b)    do{ a = (void*)(((intptr_t)(a)) ^ ((intptr_t)(b))); \
57                                  b = (void*)(((intptr_t)(a)) ^ ((intptr_t)(b))); \
58                                  a = (void*)(((intptr_t)(a)) ^ ((intptr_t)(b)));}while(0)
59
60 #define GET_SIGN(a) ((a)->info&BIGINT_NEG_MASK)
61
62 /******************************************************************************/
63 void bigint_adjust(bigint_t* a){
64         while(a->length_W!=0 && a->wordv[a->length_W-1]==0){
65                 a->length_W--;
66         }
67         if(a->length_W==0){
68                 a->info=0;
69                 return;
70         }
71         bigint_word_t t;
72         uint8_t i = BIGINT_WORD_SIZE-1;
73         t = a->wordv[a->length_W-1];
74         while((t&(1L<<(BIGINT_WORD_SIZE-1)))==0 && i){
75                 t<<=1;
76                 i--;
77         }
78         SET_FBS(a, i);
79 }
80
81 /******************************************************************************/
82
83 uint16_t bigint_length_b(const bigint_t* a){
84         if(!a->length_W || a->length_W==0){
85                 return 0;
86         }
87         return (a->length_W-1) * BIGINT_WORD_SIZE + GET_FBS(a);
88 }
89
90 /******************************************************************************/
91
92 uint16_t bigint_length_B(const bigint_t* a){
93         return a->length_W * sizeof(bigint_word_t);
94 }
95
96 /******************************************************************************/
97
98 uint32_t bigint_get_first_set_bit(const bigint_t* a){
99         if(a->length_W==0){
100                 return (uint32_t)(-1);
101         }
102         return (a->length_W-1)*sizeof(bigint_word_t)*8+GET_FBS(a);
103 }
104
105
106 /******************************************************************************/
107
108 uint32_t bigint_get_last_set_bit(const bigint_t* a){
109         uint32_t r=0;
110         uint8_t b=0;
111         bigint_word_t x=1;
112         if(a->length_W==0){
113                 return (uint32_t)(-1);
114         }
115         while(a->wordv[r]==0 && r<a->length_W){
116                 ++r;
117         }
118         if(a->wordv[r] == 0){
119                 return (uint32_t)(-1);
120         }
121         while((x&a->wordv[r])==0){
122                 ++b;
123                 x <<= 1;
124         }
125         return r*BIGINT_WORD_SIZE+b;
126 }
127
128 /******************************************************************************/
129
130 void bigint_copy(bigint_t* dest, const bigint_t* src){
131         memcpy(dest->wordv, src->wordv, src->length_W*sizeof(bigint_word_t));
132         dest->length_W = src->length_W;
133         dest->info = src->info;
134 }
135
136 /******************************************************************************/
137
138 /* this should be implemented in assembly */
139 void bigint_add_u(bigint_t* dest, const bigint_t* a, const bigint_t* b){
140         uint16_t i;
141         bigint_wordplus_t t=0LL;
142         if(a->length_W < b->length_W){
143                 XCHG_PTR(a,b);
144         }
145         for(i=0; i<b->length_W; ++i){
146 //              t = (bigint_wordplus_t)(a->wordv[i]) + (bigint_wordplus_t)(b->wordv[i]) + t;
147                 t += a->wordv[i];
148                 t += b->wordv[i];
149                 dest->wordv[i] = (bigint_word_t)t;
150                 t>>=BIGINT_WORD_SIZE;
151         }
152         for(; i<a->length_W; ++i){
153                 t += a->wordv[i];
154                 dest->wordv[i] = (bigint_word_t)t;
155                 t>>=BIGINT_WORD_SIZE;
156         }
157         if(t){
158                 dest->wordv[i++] = (bigint_word_t)t;
159         }
160         dest->length_W = i;
161         bigint_adjust(dest);
162 }
163
164 /******************************************************************************/
165
166 /* this should be implemented in assembly */
167 void bigint_add_scale_u(bigint_t* dest, const bigint_t* a, uint16_t scale){
168         if(a->length_W == 0){
169                 return;
170         }
171         if(scale == 0){
172                 bigint_add_u(dest, dest, a);
173                 return;
174         }
175         bigint_t x;
176 #if BIGINT_WORD_SIZE == 8
177         memset(dest->wordv + dest->length_W, 0, MAX(dest->length_W, a->length_W + scale) - dest->length_W);
178         x.wordv = dest->wordv + scale;
179         x.length_W = dest->length_W - scale;
180         if((int16_t)x.length_W < 0){
181                 x.length_W = 0;
182                 x.info = 0;
183         } else {
184                 x.info = dest->info;
185         }
186         bigint_add_u(&x, &x, a);
187         dest->length_W = x.length_W + scale;
188         dest->info = 0;
189         bigint_adjust(dest);
190 #else
191         bigint_t s;
192         uint16_t word_shift = scale / sizeof(bigint_word_t), byte_shift = scale % sizeof(bigint_word_t);
193         bigint_word_t bv[a->length_W + 1];
194         s.wordv = bv;
195         bv[0] = bv[a->length_W] = 0;
196         memcpy((uint8_t*)bv + byte_shift, a->wordv, a->length_W * sizeof(bigint_word_t));
197         s.length_W = a->length_W + 1;
198         bigint_adjust(&s);
199         memset(dest->wordv + dest->length_W, 0, (MAX(dest->length_W, s.length_W + word_shift) - dest->length_W) * sizeof(bigint_word_t));
200         x.wordv = dest->wordv + word_shift;
201         x.length_W = dest->length_W - word_shift;
202         if((int16_t)x.length_W < 0){
203                 x.length_W = 0;
204                 x.info = 0;
205         }else{
206                 x.info = dest->info;
207         }
208         bigint_add_u(&x, &x, &s);
209         dest->length_W = x.length_W + word_shift;
210         dest->info = 0;
211         bigint_adjust(dest);
212 #endif
213
214
215 /*      uint16_t i,j=0;
216         uint16_t scale_w;
217         bigint_word_t *dst;
218         bigint_wordplus_t t=0;
219         scale_w = (scale+sizeof(bigint_word_t)-1)/sizeof(bigint_word_t);
220         if(scale>dest->length_W*sizeof(bigint_word_t)){
221                 memset(((uint8_t*)dest->wordv)+dest->length_W*sizeof(bigint_word_t), 0, scale-dest->length_W*sizeof(bigint_word_t));
222         }
223         // a->wordv = (const uint32_t*)(((uint8_t*)a->wordv)+(scale&3));
224         dst  = dest->wordv + (scale&(sizeof(bigint_word_t)-1));
225         for(i=scale/sizeof(bigint_word_t); i<a->length_W+scale_w; ++i,++j){
226                 t += a->wordv[j];
227                 if(dest->length_W>i){
228                         t += dst[i];
229                 }
230                 dst[i] = (bigint_word_t)t;
231                 t>>=BIGINT_WORD_SIZE;
232         }
233         while(t){
234                 if(dest->length_W>i){
235                         t += dst[i];
236                 }
237                 dst[i] = (bigint_word_t)t;
238                 t>>=BIGINT_WORD_SIZE;
239                 ++i;
240         }
241         if(dest->length_W < i){
242                 dest->length_W = i;
243         }
244         bigint_adjust(dest);
245         */
246 }
247
248 /******************************************************************************/
249
250 /* this should be implemented in assembly */
251 void bigint_sub_u(bigint_t* dest, const bigint_t* a, const bigint_t* b){
252         int8_t borrow=0;
253         int8_t  r;
254         bigint_wordplus_signed_t t=0LL;
255         uint16_t i, min, max;
256         min = MIN(a->length_W, b->length_W);
257         max = MAX(a->length_W, b->length_W);
258         r = bigint_cmp_u(a,b);
259         if(r==0){
260                 bigint_set_zero(dest);
261                 return;
262         }
263         if(b->length_W==0){
264                 bigint_copy(dest, a);
265                 SET_POS(dest);
266                 return;
267         }
268         if(a->length_W==0){
269                 bigint_copy(dest, b);
270                 SET_NEG(dest);
271                 return;
272         }
273         if(r<0){
274                 bigint_sub_u(dest, b, a);
275                 SET_NEG(dest);
276                 return;
277         }
278         for(i=0; i<max; ++i){
279                 t = a->wordv[i];
280                 if(i<min){
281                         t -= b->wordv[i];
282                 }
283                 t -= borrow;
284                 dest->wordv[i]=(bigint_word_t)t;
285                 if(t<0){
286                         borrow = 1;
287                 }else{
288                         borrow = 0;
289                 }
290         }
291         SET_POS(dest);
292         dest->length_W = i;
293         bigint_adjust(dest);
294 }
295
296 /******************************************************************************/
297
298 int8_t bigint_cmp_u(const bigint_t* a, const bigint_t* b){
299         if(a->length_W > b->length_W){
300                 return 1;
301         }
302         if(a->length_W < b->length_W){
303                 return -1;
304         }
305         if(a->length_W==0){
306                 return 0;
307         }
308         uint16_t i;
309         i = a->length_W-1;
310         do{
311                 if(a->wordv[i] != b->wordv[i]){
312                         if(a->wordv[i] > b->wordv[i]){
313                                 return 1;
314                         }else{
315                                 return -1;
316                         }
317                 }
318         }while(i--);
319         return 0;
320 }
321
322 /******************************************************************************/
323
324 void bigint_add_s(bigint_t* dest, const bigint_t* a, const bigint_t* b){
325         uint8_t s;
326         s  = GET_SIGN(a)?2:0;
327         s |= GET_SIGN(b)?1:0;
328         switch(s){
329                 case 0: /* both positive */
330                         bigint_add_u(dest, a,b);
331                         SET_POS(dest);
332                         break;
333                 case 1: /* a positive, b negative */
334                         bigint_sub_u(dest, a, b);
335                         break;
336                 case 2: /* a negative, b positive */
337                         bigint_sub_u(dest, b, a);
338                         break;
339                 case 3: /* both negative */
340                         bigint_add_u(dest, a, b);
341                         SET_NEG(dest);
342                         break;
343                 default: /* how can this happen?*/
344                         break;
345         }
346 }
347
348 /******************************************************************************/
349
350 void bigint_sub_s(bigint_t* dest, const bigint_t* a, const bigint_t* b){
351         uint8_t s;
352         s  = GET_SIGN(a)?2:0;
353         s |= GET_SIGN(b)?1:0;
354         switch(s){
355                 case 0: /* both positive */
356                         bigint_sub_u(dest, a,b);
357                         break;
358                 case 1: /* a positive, b negative */
359                         bigint_add_u(dest, a, b);
360                         SET_POS(dest);
361                         break;
362                 case 2: /* a negative, b positive */
363                         bigint_add_u(dest, a, b);
364                         SET_NEG(dest);
365                         break;
366                 case 3: /* both negative */
367                         bigint_sub_u(dest, b, a);
368                         break;
369                 default: /* how can this happen?*/
370                                         break;
371         }
372
373 }
374
375 /******************************************************************************/
376
377 int8_t bigint_cmp_s(const bigint_t* a, const bigint_t* b){
378         uint8_t s;
379         if(a->length_W==0 && b->length_W==0){
380                 return 0;
381         }
382         s  = GET_SIGN(a)?2:0;
383         s |= GET_SIGN(b)?1:0;
384         switch(s){
385                 case 0: /* both positive */
386                         return bigint_cmp_u(a, b);
387                         break;
388                 case 1: /* a positive, b negative */
389                         return 1;
390                         break;
391                 case 2: /* a negative, b positive */
392                         return -1;
393                         break;
394                 case 3: /* both negative */
395                         return bigint_cmp_u(b, a);
396                         break;
397                 default: /* how can this happen?*/
398                                         break;
399         }
400         return 0; /* just to satisfy the compiler */
401 }
402
403 /******************************************************************************/
404
405 void bigint_shiftleft(bigint_t* a, uint16_t shift){
406         uint16_t byteshift, word_alloc, words_to_shift;
407         int16_t i;
408         uint8_t bitshift;
409         bigint_word_t *p;
410         bigint_wordplus_t t=0;
411         if(shift==0){
412                 return;
413         }
414         byteshift = shift/8;
415         bitshift = shift&7;
416         for(i=0;i<=byteshift/sizeof(bigint_word_t); ++i){
417                 a->wordv[a->length_W+i] = 0;
418         }
419         if(byteshift){
420                 memmove(((uint8_t*)a->wordv) + byteshift, a->wordv, a->length_W * sizeof(bigint_word_t));
421                 memset(a->wordv, 0, byteshift);
422         }
423         p = a->wordv + byteshift / sizeof(bigint_word_t);
424         words_to_shift = a->length_W + (byteshift % sizeof(bigint_word_t)?1:0);
425         word_alloc = a->length_W + (byteshift + sizeof(bigint_word_t) - 1) / sizeof(bigint_word_t) + 1;
426         a->wordv[word_alloc-1]=0;
427         if(bitshift!=0){
428                 for(i=0; i < words_to_shift; ++i){
429                         t |= ((bigint_wordplus_t)p[i])<<bitshift;
430                         p[i] = (bigint_word_t)t;
431                         t >>= BIGINT_WORD_SIZE;
432                 }
433                 p[i] = (bigint_word_t)t;
434         }
435         a->length_W = word_alloc;
436         bigint_adjust(a);
437 }
438
439 /******************************************************************************/
440
441 void bigint_shiftright(bigint_t* a, uint16_t shift){
442         uint16_t byteshift;
443         uint16_t i;
444         uint8_t bitshift;
445         bigint_wordplus_t t=0;
446         byteshift = shift/8;
447         bitshift = shift&7;
448         if(byteshift >= a->length_W * sizeof(bigint_word_t)){ /* we would shift out more than we have */
449                 bigint_set_zero(a);
450                 return;
451         }
452         if(byteshift == a->length_W * sizeof(bigint_word_t) - 1 && bitshift > GET_FBS(a)){
453                 bigint_set_zero(a);
454                 return;
455         }
456         if(byteshift){
457                 memmove(a->wordv, (uint8_t*)a->wordv + byteshift, a->length_W * sizeof(bigint_word_t) - byteshift);
458                 memset((uint8_t*)a->wordv + a->length_W * sizeof(bigint_word_t) - byteshift, 0,  byteshift);
459         }
460         byteshift /= sizeof(bigint_word_t);
461     a->length_W -= (byteshift  + sizeof(bigint_word_t) - 1) / sizeof(bigint_word_t);
462         if(bitshift != 0 && a->length_W){
463          /* shift to the right */
464                 i = a->length_W - 1;
465                 do{
466                         t |= ((bigint_wordplus_t)(a->wordv[i])) << (BIGINT_WORD_SIZE - bitshift);
467                         a->wordv[i] = (bigint_word_t)(t >> BIGINT_WORD_SIZE);
468                         t <<= BIGINT_WORD_SIZE;
469                 }while(i--);
470         }
471         bigint_adjust(a);
472 }
473
474 /******************************************************************************/
475
476 void bigint_xor(bigint_t* dest, const bigint_t* a){
477         uint16_t i;
478         for(i=0; i<a->length_W; ++i){
479                 dest->wordv[i] ^= a->wordv[i];
480         }
481         bigint_adjust(dest);
482 }
483
484 /******************************************************************************/
485
486 void bigint_set_zero(bigint_t* a){
487         a->length_W=0;
488 }
489
490 /******************************************************************************/
491
492 /* using the Karatsuba-Algorithm */
493 /* x*y = (xh*yh)*b**2n + ((xh+xl)*(yh+yl) - xh*yh - xl*yl)*b**n + yh*yl */
494 void bigint_mul_u(bigint_t* dest, const bigint_t* a, const bigint_t* b){
495         if(a->length_W==0 || b->length_W==0){
496                 bigint_set_zero(dest);
497                 return;
498         }
499         if(dest==a || dest==b){
500                 bigint_t d;
501                 bigint_word_t d_b[a->length_W+b->length_W];
502                 d.wordv = d_b;
503                 bigint_mul_u(&d, a, b);
504                 bigint_copy(dest, &d);
505                 return;
506         }
507         if(a->length_W==1 || b->length_W==1){
508                 if(a->length_W!=1){
509                         XCHG_PTR(a,b);
510                 }
511                 bigint_wordplus_t t=0;
512                 uint16_t i;
513                 bigint_word_t x = a->wordv[0];
514                 for(i=0; i < b->length_W; ++i){
515                         t += ((bigint_wordplus_t)b->wordv[i])*((bigint_wordplus_t)x);
516                         dest->wordv[i] = (bigint_word_t)t;
517                         t>>=BIGINT_WORD_SIZE;
518                 }
519                 dest->wordv[i] = (bigint_word_t)t;
520                 dest->length_W = i+1;
521                 dest->info = 0;
522                 bigint_adjust(dest);
523                 return;
524         }
525         if(a->length_W * sizeof(bigint_word_t) <= 4 && b->length_W * sizeof(bigint_word_t) <= 4){
526                 uint32_t p=0, q=0;
527                 uint64_t r;
528                 memcpy(&p, a->wordv, a->length_W*sizeof(bigint_word_t));
529                 memcpy(&q, b->wordv, b->length_W*sizeof(bigint_word_t));
530                 r = (uint64_t)p * (uint64_t)q;
531                 memcpy(dest->wordv, &r, (dest->length_W = a->length_W + b->length_W)*sizeof(bigint_word_t));
532                 bigint_adjust(dest);
533                 return;
534         }
535         bigint_set_zero(dest);
536         /* split a in xh & xl; split b in yh & yl */
537         const uint16_t n = (MAX(a->length_W, b->length_W)+1)/2;
538         bigint_t xl, xh, yl, yh;
539         xl.wordv = a->wordv;
540         yl.wordv = b->wordv;
541         if(a->length_W<=n){
542                 bigint_set_zero(&xh);
543                 xl.length_W = a->length_W;
544                 xl.info = a->info;
545         }else{
546                 xl.length_W=n;
547                 xl.info = 0;
548                 bigint_adjust(&xl);
549                 xh.wordv = &(a->wordv[n]);
550                 xh.length_W = a->length_W-n;
551                 xh.info = a->info;
552         }
553         if(b->length_W<=n){
554                 bigint_set_zero(&yh);
555                 yl.length_W = b->length_W;
556                 yl.info = b->info;
557         }else{
558                 yl.length_W=n;
559                 yl.info = 0;
560                 bigint_adjust(&yl);
561                 yh.wordv = &(b->wordv[n]);
562                 yh.length_W = b->length_W-n;
563                 yh.info = b->info;
564         }
565         /* now we have split up a and b */
566         /* remember we want to do:
567          * x*y = (xh*yh)*b**2n + ((xh+xl)*(yh+yl) - xh*yh - xl*yl)*b**n + yh*yl
568          *          5          9     2   4   3    7   5   6   1         8   1
569          */
570         bigint_word_t  tmp_b[2*n+2], m_b[2*(n+1)];
571         bigint_t tmp, tmp2, m;
572         tmp.wordv = tmp_b;
573         tmp2.wordv = &(tmp_b[n+1]);
574         m.wordv = m_b;
575
576         bigint_mul_u(dest, &xl, &yl);  /* 1: dest <= xl*yl     */
577         bigint_add_u(&tmp2, &xh, &xl); /* 2: tmp2 <= xh+xl     */
578         bigint_add_u(&tmp, &yh, &yl);  /* 3: tmp  <= yh+yl     */
579         bigint_mul_u(&m, &tmp2, &tmp); /* 4: m    <= tmp2*tmp  */
580         bigint_mul_u(&tmp, &xh, &yh);  /* 5: h    <= xh*yh     */
581         bigint_sub_u(&m, &m, dest);    /* 6: m    <= m-dest    */
582     bigint_sub_u(&m, &m, &tmp);    /* 7: m    <= m-h       */
583         bigint_add_scale_u(dest, &m, n*sizeof(bigint_word_t)); /* 8: dest <= dest+m**n*/
584         bigint_add_scale_u(dest, &tmp, 2*n*sizeof(bigint_word_t)); /* 9: dest <= dest+tmp**(2*n) */
585 }
586
587 /******************************************************************************/
588
589 void bigint_mul_s(bigint_t* dest, const bigint_t* a, const bigint_t* b){
590         uint8_t s;
591         s  = GET_SIGN(a)?2:0;
592         s |= GET_SIGN(b)?1:0;
593         switch(s){
594                 case 0: /* both positive */
595                         bigint_mul_u(dest, a,b);
596                         SET_POS(dest);
597                         break;
598                 case 1: /* a positive, b negative */
599                         bigint_mul_u(dest, a,b);
600                         SET_NEG(dest);
601                         break;
602                 case 2: /* a negative, b positive */
603                         bigint_mul_u(dest, a,b);
604                         SET_NEG(dest);
605                         break;
606                 case 3: /* both negative */
607                         bigint_mul_u(dest, a,b);
608                         SET_POS(dest);
609                         break;
610                 default: /* how can this happen?*/
611                         break;
612         }
613 }
614
615 /******************************************************************************/
616
617 /* square */
618 /* (xh*b^n+xl)^2 = xh^2*b^2n + 2*xh*xl*b^n + xl^2 */
619 void bigint_square(bigint_t* dest, const bigint_t* a){
620         if(a->length_W*sizeof(bigint_word_t)<=4){
621                 uint64_t r=0;
622                 memcpy(&r, a->wordv, a->length_W*sizeof(bigint_word_t));
623                 r = r*r;
624                 memcpy(dest->wordv, &r, 2*a->length_W*sizeof(bigint_word_t));
625                 SET_POS(dest);
626                 dest->length_W=2*a->length_W;
627                 bigint_adjust(dest);
628                 return;
629         }
630         if(dest==a){
631                 bigint_t d;
632                 bigint_word_t d_b[a->length_W*2];
633                 d.wordv = d_b;
634                 bigint_square(&d, a);
635                 bigint_copy(dest, &d);
636                 return;
637         }
638         uint16_t n;
639         n=(a->length_W+1)/2;
640         bigint_t xh, xl, tmp; /* x-high, x-low, temp */
641         bigint_word_t buffer[2*n+1];
642         xl.wordv = a->wordv;
643         xl.length_W = n;
644         xl.info = 0;
645         xh.wordv = &(a->wordv[n]);
646         xh.length_W = a->length_W-n;
647         xh.info = 0;
648         bigint_adjust(&xl);
649         bigint_adjust(&xh);
650         tmp.wordv = buffer;
651 /* (xh * b**n + xl)**2 = xh**2 * b**2n + 2 * xh * xl * b**n + xl**2 */
652
653 //      cli_putstr("\r\nDBG (a): xl: "); bigint_print_hex(&xl);
654 //      cli_putstr("\r\nDBG (b): xh: "); bigint_print_hex(&xh);
655         bigint_square(dest, &xl);
656 //      cli_putstr("\r\nDBG (1): xl**2: "); bigint_print_hex(dest);
657         bigint_square(&tmp, &xh);
658 //      cli_putstr("\r\nDBG (2): xh**2: "); bigint_print_hex(&tmp);
659         bigint_add_scale_u(dest, &tmp, 2*n*sizeof(bigint_word_t));
660 //      cli_putstr("\r\nDBG (3): xl**2 + xh**2*n**2: "); bigint_print_hex(dest);
661         bigint_mul_u(&tmp, &xl, &xh);
662 //      cli_putstr("\r\nDBG (4): xl*xh: "); bigint_print_hex(&tmp);
663         bigint_shiftleft(&tmp, 1);
664 //      cli_putstr("\r\nDBG (5): xl*xh*2: "); bigint_print_hex(&tmp);
665         bigint_add_scale_u(dest, &tmp, n*sizeof(bigint_word_t));
666 //      cli_putstr("\r\nDBG (6): x**2: "); bigint_print_hex(dest);
667 //      cli_putstr("\r\n");
668 }
669
670 /******************************************************************************/
671 void bigint_sub_u_bitscale(bigint_t* a, const bigint_t* b, uint16_t bitscale){
672         bigint_t tmp, x;
673         bigint_word_t tmp_b[b->length_W + 1];
674         const uint16_t word_shift = bitscale / BIGINT_WORD_SIZE;
675
676         if(a->length_W < b->length_W + word_shift){
677 #if DEBUG
678                 cli_putstr("\r\nDBG: *bang*\r\n");
679 #endif
680                 bigint_set_zero(a);
681                 return;
682         }
683         tmp.wordv = tmp_b;
684         bigint_copy(&tmp, b);
685         bigint_shiftleft(&tmp, bitscale % BIGINT_WORD_SIZE);
686
687         x.info = a->info;
688         x.wordv = &(a->wordv[word_shift]);
689         x.length_W = a->length_W - word_shift;
690
691         bigint_sub_u(&x, &x, &tmp);
692         bigint_adjust(a);
693         return;
694 }
695
696 /******************************************************************************/
697
698 void bigint_reduce(bigint_t* a, const bigint_t* r){
699 //      bigint_adjust((bigint_t*)r);
700         uint8_t rfbs = GET_FBS(r);
701 #if DEBUG
702         cli_putstr("\r\nDBG: (a) = "); bigint_print_hex(a);
703 #endif
704         if(r->length_W==0 || a->length_W==0){
705                 return;
706         }
707         if((r->length_W*sizeof(bigint_word_t)<=4) && (a->length_W*sizeof(bigint_word_t)<=4)){
708                 uint32_t p=0, q=0;
709                 memcpy(&p, a->wordv, a->length_W*sizeof(bigint_word_t));
710                 memcpy(&q, r->wordv, r->length_W*sizeof(bigint_word_t));
711                 p %= q;
712                 memcpy(a->wordv, &p, a->length_W*sizeof(bigint_word_t));
713                 bigint_adjust(a);
714 //              cli_putstr("\r\nDBG: (0) = "); bigint_print_hex(a);
715                 return;
716         }
717         uint16_t shift;
718         while(a->length_W > r->length_W){
719                 shift = (a->length_W - r->length_W) * 8 * sizeof(bigint_word_t) + GET_FBS(a) - rfbs - 1;
720                 /*
721                 if((a->wordv[a->length_W-1] & ((1LL<<GET_FBS(a)) - 1)) > r->wordv[r->length_W-1]){
722                         // cli_putc('~');
723                         cli_putstr("\r\n ~ [a] = ");
724                         cli_hexdump_rev(&a->wordv[a->length_W-1], 4);
725                         cli_putstr("  [r] = ");
726                         cli_hexdump_rev(&r->wordv[r->length_W-1], 4);
727                         shift += 1;
728                 }
729                 */
730 //              cli_putstr("\r\nDBG: (p) shift = "); cli_hexdump_rev(&shift, 2);
731 //              cli_putstr(" a_len = "); cli_hexdump_rev(&a->length_W, 2);
732 //              cli_putstr(" r_len = "); cli_hexdump_rev(&r->length_W, 2);
733 //              uart_flush(0);
734                 bigint_sub_u_bitscale(a, r, shift);
735 //              cli_putstr("\r\nDBG: (1) = "); bigint_print_hex(a);
736         }
737         while((GET_FBS(a) > rfbs) && (a->length_W == r->length_W)){
738                 shift = GET_FBS(a)-rfbs-1;
739 //              cli_putstr("\r\nDBG: (q) shift = "); cli_hexdump_rev(&shift, 2);
740                 bigint_sub_u_bitscale(a, r, shift);
741 //              cli_putstr("\r\nDBG: (2) = "); bigint_print_hex(a);
742         }
743         while(bigint_cmp_u(a,r)>=0){
744                 bigint_sub_u(a,a,r);
745 //              cli_putstr("\r\nDBG: (3) = "); bigint_print_hex(a);
746         }
747         bigint_adjust(a);
748 //      cli_putstr("\r\nDBG: (a) = "); bigint_print_hex(a);
749 //      cli_putstr("\r\n");
750 }
751
752 /******************************************************************************/
753
754 /* calculate dest = a**exp % r */
755 /* using square&multiply */
756 void bigint_expmod_u(bigint_t* dest, const bigint_t* a, const bigint_t* exp, const bigint_t* r){
757         if(a->length_W==0 || r->length_W==0){
758                 return;
759         }
760
761         bigint_t res, base;
762         bigint_word_t t, base_b[MAX(a->length_W,r->length_W)], res_b[r->length_W*2];
763         uint16_t i;
764         uint8_t j;
765 //      uint16_t *xaddr = &i;
766 //      cli_putstr("\r\npre-alloc (");
767 //      cli_hexdump_rev(&xaddr, 4);
768 //      cli_putstr(") ...");
769         res.wordv = res_b;
770         base.wordv = base_b;
771         bigint_copy(&base, a);
772 //      cli_putstr("\r\npost-copy");
773         bigint_reduce(&base, r);
774         res.wordv[0]=1;
775         res.length_W=1;
776         res.info = 0;
777         bigint_adjust(&res);
778         if(exp->length_W == 0){
779                 bigint_copy(dest, &res);
780                 return;
781         }
782         uint8_t flag = 0;
783         t=exp->wordv[exp->length_W - 1];
784         for(i=exp->length_W; i > 0; --i){
785                 t = exp->wordv[i - 1];
786                 for(j=BIGINT_WORD_SIZE; j > 0; --j){
787                         if(!flag){
788                                 if(t & (1<<(BIGINT_WORD_SIZE-1))){
789                                         flag = 1;
790                                 }
791                         }
792                         if(flag){
793                                 bigint_square(&res, &res);
794                                 bigint_reduce(&res, r);
795                                 if(t & (1 << (BIGINT_WORD_SIZE - 1))){
796                                         bigint_mul_u(&res, &res, &base);
797                                         bigint_reduce(&res, r);
798                                 }
799                         }
800                         t <<= 1;
801                 }
802         }
803
804 //      cli_putc('+');
805         SET_POS(&res);
806         bigint_copy(dest, &res);
807 }
808
809 /******************************************************************************/
810
811 #define cli_putstr(a)
812 #define bigint_print_hex(a)
813 #define cli_hexdump_rev(a,b)
814 #define uart_flush(a)
815
816 /* gcd <-- gcd(x,y) a*x+b*y=gcd */
817 void bigint_gcdext(bigint_t* gcd, bigint_t* a, bigint_t* b, const bigint_t* x, const bigint_t* y){
818          bigint_t g, x_, y_, u, v, a_, b_, c_, d_;
819          uint16_t i=0;
820          if(x->length_W==0 || y->length_W==0){
821                  return;
822          }
823          if(x->length_W==1 && x->wordv[0]==1){
824                  gcd->length_W = 1;
825                  gcd->wordv[0] = 1;
826                  if(a){
827                          a->length_W = 1;
828                          a->wordv[0] = 1;
829                          SET_POS(a);
830                          bigint_adjust(a);
831                  }
832                  if(b){
833                          bigint_set_zero(b);
834                  }
835                  return;
836          }
837          if(y->length_W==1 && y->wordv[0]==1){
838                  gcd->length_W = 1;
839                  gcd->wordv[0] = 1;
840                  if(b){
841                          b->length_W = 1;
842                          b->wordv[0] = 1;
843                          SET_POS(b);
844                          bigint_adjust(b);
845                  }
846                  if(a){
847                          bigint_set_zero(a);
848                  }
849                  return;
850          }
851
852          while(x->wordv[i]==0 && y->wordv[i]==0){
853                  ++i;
854          }
855          bigint_word_t g_b[i+2], x_b[x->length_W-i], y_b[y->length_W-i];
856          bigint_word_t u_b[x->length_W-i], v_b[y->length_W-i];
857          bigint_word_t a_b[y->length_W+2], c_b[y->length_W+2];
858          bigint_word_t b_b[x->length_W+2], d_b[x->length_W+2];
859
860          g.wordv = g_b;
861          x_.wordv = x_b;
862          y_.wordv = y_b;
863          memset(g_b, 0, i*sizeof(bigint_word_t));
864          g_b[i]=1;
865          g.length_W = i+1;
866          g.info=0;
867          x_.info = y_.info = 0;
868          x_.length_W = x->length_W-i;
869          y_.length_W = y->length_W-i;
870          memcpy(x_.wordv, x->wordv+i, x_.length_W*sizeof(bigint_word_t));
871          memcpy(y_.wordv, y->wordv+i, y_.length_W*sizeof(bigint_word_t));
872          for(i=0; (x_.wordv[0]&(1<<i))==0 && (y_.wordv[0]&(1<<i))==0; ++i){
873          }
874
875          bigint_adjust(&x_);
876          bigint_adjust(&y_);
877
878          if(i){
879                  bigint_shiftleft(&g, i);
880                  bigint_shiftright(&x_, i);
881                  bigint_shiftright(&y_, i);
882          }
883
884          u.wordv = u_b;
885          v.wordv = v_b;
886          a_.wordv = a_b;
887          b_.wordv = b_b;
888          c_.wordv = c_b;
889          d_.wordv = d_b;
890
891          bigint_copy(&u, &x_);
892          bigint_copy(&v, &y_);
893          a_.wordv[0] = 1;
894          a_.length_W = 1;
895          a_.info = 0;
896          d_.wordv[0] = 1;
897          d_.length_W = 1;
898          d_.info = 0;
899          bigint_set_zero(&b_);
900          bigint_set_zero(&c_);
901          do{
902                  cli_putstr("\r\nDBG (gcdext) 0");
903                  while((u.wordv[0]&1)==0){
904                          cli_putstr("\r\nDBG (gcdext) 0.1");
905                          bigint_shiftright(&u, 1);
906                          if((a_.wordv[0]&1) || (b_.wordv[0]&1)){
907                                  bigint_add_s(&a_, &a_, &y_);
908                                  bigint_sub_s(&b_, &b_, &x_);
909                          }
910                          bigint_shiftright(&a_, 1);
911                          bigint_shiftright(&b_, 1);
912                  }
913                  while((v.wordv[0]&1)==0){
914                          cli_putstr("\r\nDBG (gcdext) 0.2");
915                          bigint_shiftright(&v, 1);
916                          if((c_.wordv[0]&1) || (d_.wordv[0]&1)){
917                                  bigint_add_s(&c_, &c_, &y_);
918                                  bigint_sub_s(&d_, &d_, &x_);
919                          }
920                          bigint_shiftright(&c_, 1);
921                          bigint_shiftright(&d_, 1);
922
923                  }
924                  if(bigint_cmp_u(&u, &v)>=0){
925                         bigint_sub_u(&u, &u, &v);
926                         bigint_sub_s(&a_, &a_, &c_);
927                         bigint_sub_s(&b_, &b_, &d_);
928                  }else{
929                         bigint_sub_u(&v, &v, &u);
930                         bigint_sub_s(&c_, &c_, &a_);
931                         bigint_sub_s(&d_, &d_, &b_);
932                  }
933          }while(u.length_W);
934          if(gcd){
935                  bigint_mul_s(gcd, &v, &g);
936          }
937          if(a){
938                 bigint_copy(a, &c_);
939          }
940          if(b){
941                  bigint_copy(b, &d_);
942          }
943 }
944
945 /******************************************************************************/
946
947 void bigint_inverse(bigint_t* dest, const bigint_t* a, const bigint_t* m){
948         bigint_gcdext(NULL, dest, NULL, a, m);
949         while(dest->info&BIGINT_NEG_MASK){
950                 bigint_add_s(dest, dest, m);
951         }
952 }
953
954 /******************************************************************************/
955
956 void bigint_changeendianess(bigint_t* a){
957         uint8_t t, *p, *q;
958         p = (uint8_t*)(a->wordv);
959         q = p + a->length_W * sizeof(bigint_word_t) - 1;
960         while(p<q){
961                 t = *p;
962                 *p = *q;
963                 *q = t;
964                 ++p; --q;
965         }
966 }
967
968 /******************************************************************************/
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991