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*)(((bigint_ptr_int_t)(a)) ^ ((bigint_ptr_int_t)(b))); \
  57                                  b = (void*)(((bigint_ptr_int_t)(a)) ^ ((bigint_ptr_int_t)(b))); \
  58                                  a = (void*)(((bigint_ptr_int_t)(a)) ^ ((bigint_ptr_int_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_B!=0 && a->wordv[a->length_B-1]==0){
  65                 a->length_B--;
  66         }
  67         if(a->length_B==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_B-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(bigint_t* a){
  84         if(!a->length_B || a->length_B==0){
  85                 return 0;
  86         }
  87         return (a->length_B-1) * BIGINT_WORD_SIZE + GET_FBS(a);
  88 }
  89
  90 /******************************************************************************/
  91
  92 uint16_t bigint_length_B(bigint_t* a){
  93         return a->length_B * sizeof(bigint_word_t);
  94 }
  95
  96 /******************************************************************************/
  97
  98 uint32_t bigint_get_first_set_bit(bigint_t* a){
  99         if(a->length_B==0){
 100                 return (uint32_t)(-1);
 101         }
 102         return (a->length_B-1)*sizeof(bigint_word_t)*8+GET_FBS(a);
 103 }
 104
 105
 106 /******************************************************************************/
 107
 108 uint32_t bigint_get_last_set_bit(bigint_t* a){
 109         uint32_t r=0;
 110         uint8_t b=0;
 111         bigint_word_t x=1;
 112         if(a->length_B==0){
 113                 return (uint32_t)(-1);
 114         }
 115         while(a->wordv[r]==0 && r<a->length_B){
 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_B*sizeof(bigint_word_t));
 132         dest->length_B = src->length_B;
 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_B < b->length_B){
 143                 XCHG_PTR(a,b);
 144         }
 145         for(i=0; i<b->length_B; ++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_B; ++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_B = 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_B == 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_B, 0, MAX(dest->length_B, a->length_B + scale) - dest->length_B);
 178         x.wordv = dest->wordv + scale;
 179         x.length_B = dest->length_B - scale;
 180         if((int16_t)x.length_B < 0){
 181                 x.length_B = 0;
 182                 x.info = 0;
 183         } else {
 184                 x.info = dest->info;
 185         }
 186         bigint_add_u(&x, &x, a);
 187         dest->length_B = x.length_B + 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_B + 1];
 194         s.wordv = bv;
 195         bv[0] = bv[a->length_B] = 0;
 196         memcpy((uint8_t*)bv + byte_shift, a->wordv, a->length_B * sizeof(bigint_word_t));
 197         s.length_B = a->length_B + 1;
 198         bigint_adjust(&s);
 199         memset(dest->wordv + dest->length_B, 0, (MAX(dest->length_B, s.length_B + word_shift) - dest->length_B) * sizeof(bigint_word_t));
 200         x.wordv = dest->wordv + word_shift;
 201         x.length_B = dest->length_B - word_shift;
 202         if((int16_t)x.length_B < 0){
 203                 x.length_B = 0;
 204                 x.info = 0;
 205         }else{
 206                 x.info = dest->info;
 207         }
 208         bigint_add_u(&x, &x, &s);
 209         dest->length_B = x.length_B + 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_B*sizeof(bigint_word_t)){
 221                 memset(((uint8_t*)dest->wordv)+dest->length_B*sizeof(bigint_word_t), 0, scale-dest->length_B*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_B+scale_w; ++i,++j){
 226                 t += a->wordv[j];
 227                 if(dest->length_B>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_B>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_B < i){
 242                 dest->length_B = 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_B, b->length_B);
 257         max = MAX(a->length_B, b->length_B);
 258         r = bigint_cmp_u(a,b);
 259         if(r==0){
 260                 bigint_set_zero(dest);
 261                 return;
 262         }
 263         if(b->length_B==0){
 264                 bigint_copy(dest, a);
 265                 SET_POS(dest);
 266                 return;
 267         }
 268         if(a->length_B==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_B = 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_B > b->length_B){
 300                 return 1;
 301         }
 302         if(a->length_B < b->length_B){
 303                 return -1;
 304         }
 305         if(a->length_B==0){
 306                 return 0;
 307         }
 308         uint16_t i;
 309         i = a->length_B-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_B==0 && b->length_B==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;
 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_B+i] = 0;
 418         }
 419         if(byteshift){
 420                 memmove(((uint8_t*)a->wordv)+byteshift, a->wordv, a->length_B*sizeof(bigint_word_t));
 421                 memset(a->wordv, 0, byteshift);
 422         }
 423         p = (bigint_word_t*)(((uint8_t*)a->wordv)+byteshift);
 424         word_alloc = a->length_B+(byteshift+sizeof(bigint_word_t)-1)/sizeof(bigint_word_t)+1;
 425         a->wordv[word_alloc-1]=0;
 426         if(bitshift!=0){
 427                 for(i=0; i<a->length_B; ++i){
 428                         t |= ((bigint_wordplus_t)p[i])<<bitshift;
 429                         p[i] = (bigint_word_t)t;
 430                         t >>= BIGINT_WORD_SIZE;
 431                 }
 432                 p[i] = (bigint_word_t)t;
 433         }
 434         a->length_B = word_alloc;
 435         bigint_adjust(a);
 436 }
 437
 438 /******************************************************************************/
 439
 440 void bigint_shiftright(bigint_t* a, uint16_t shift){
 441         uint16_t byteshift;
 442         uint16_t i;
 443         uint8_t bitshift;
 444         bigint_wordplus_t t=0;
 445         byteshift = shift/8;
 446         bitshift = shift&7;
 447         if(byteshift >= a->length_B*sizeof(bigint_word_t)){ /* we would shift out more than we have */
 448                 bigint_set_zero(a);
 449                 return;
 450         }
 451         if(byteshift == a->length_B*sizeof(bigint_word_t)-1 && bitshift>GET_FBS(a)){
 452                 bigint_set_zero(a);
 453                 return;
 454         }
 455         if(byteshift){
 456                 memmove(a->wordv, (uint8_t*)a->wordv+byteshift, a->length_B-byteshift);
 457                 memset((uint8_t*)a->wordv+a->length_B-byteshift, 0,  byteshift);
 458         }
 459         byteshift /= sizeof(bigint_word_t);
 460         if(bitshift!=0){
 461          /* shift to the right */
 462                 for(i=a->length_B-byteshift-1; i>0; --i){
 463                         t |= ((bigint_wordplus_t)(a->wordv[i]))<<(BIGINT_WORD_SIZE-bitshift);
 464                         a->wordv[i] = (bigint_word_t)(t>>BIGINT_WORD_SIZE);
 465                         t <<= BIGINT_WORD_SIZE;
 466                 }
 467                 t |= ((bigint_wordplus_t)(a->wordv[0]))<<(BIGINT_WORD_SIZE-bitshift);
 468                 a->wordv[0] = (bigint_word_t)(t>>BIGINT_WORD_SIZE);
 469         }
 470     a->length_B -= ((shift/8)+sizeof(bigint_word_t)-1)/sizeof(bigint_word_t);
 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_B; ++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_B=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_B==0 || b->length_B==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_B+b->length_B];
 502                 d.wordv = d_b;
 503                 bigint_mul_u(&d, a, b);
 504                 bigint_copy(dest, &d);
 505                 return;
 506         }
 507         if(a->length_B==1 || b->length_B==1){
 508                 if(a->length_B!=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_B; ++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_B = i+1;
 521                 dest->info = 0;
 522                 bigint_adjust(dest);
 523                 return;
 524         }
 525         if(a->length_B * sizeof(bigint_word_t) <= 4 && b->length_B * sizeof(bigint_word_t) <= 4){
 526                 uint32_t p=0, q=0;
 527                 uint64_t r;
 528                 memcpy(&p, a->wordv, a->length_B*sizeof(bigint_word_t));
 529                 memcpy(&q, b->wordv, b->length_B*sizeof(bigint_word_t));
 530                 r = (uint64_t)p * (uint64_t)q;
 531                 memcpy(dest->wordv, &r, (dest->length_B = a->length_B + b->length_B)*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_B, b->length_B)+1)/2;
 538         bigint_t xl, xh, yl, yh;
 539         xl.wordv = a->wordv;
 540         yl.wordv = b->wordv;
 541         if(a->length_B<=n){
 542                 bigint_set_zero(&xh);
 543                 xl.length_B = a->length_B;
 544                 xl.info = a->info;
 545         }else{
 546                 xl.length_B=n;
 547                 xl.info = 0;
 548                 bigint_adjust(&xl);
 549                 xh.wordv = &(a->wordv[n]);
 550                 xh.length_B = a->length_B-n;
 551                 xh.info = a->info;
 552         }
 553         if(b->length_B<=n){
 554                 bigint_set_zero(&yh);
 555                 yl.length_B = b->length_B;
 556                 yl.info = b->info;
 557         }else{
 558                 yl.length_B=n;
 559                 yl.info = 0;
 560                 bigint_adjust(&yl);
 561                 yh.wordv = &(b->wordv[n]);
 562                 yh.length_B = b->length_B-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_B*sizeof(bigint_word_t)<=4){
 621                 uint64_t r=0;
 622                 memcpy(&r, a->wordv, a->length_B*sizeof(bigint_word_t));
 623                 r = r*r;
 624                 memcpy(dest->wordv, &r, 2*a->length_B*sizeof(bigint_word_t));
 625                 SET_POS(dest);
 626                 dest->length_B=2*a->length_B;
 627                 bigint_adjust(dest);
 628                 return;
 629         }
 630         if(dest==a){
 631                 bigint_t d;
 632                 bigint_word_t d_b[a->length_B*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_B+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_B = n;
 644         xl.info = 0;
 645         xh.wordv = &(a->wordv[n]);
 646         xh.length_B = a->length_B-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_B + 1];
 674         const uint16_t word_shift = bitscale / BIGINT_WORD_SIZE;
 675
 676         if(a->length_B < b->length_B + 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_B = a->length_B - 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_B==0 || a->length_B==0){
 705                 return;
 706         }
 707         if((r->length_B*sizeof(bigint_word_t)<=4) && (a->length_B*sizeof(bigint_word_t)<=4)){
 708                 uint32_t p=0, q=0;
 709                 memcpy(&p, a->wordv, a->length_B*sizeof(bigint_word_t));
 710                 memcpy(&q, r->wordv, r->length_B*sizeof(bigint_word_t));
 711                 p %= q;
 712                 memcpy(a->wordv, &p, a->length_B*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_B > r->length_B){
 719                 shift = (a->length_B - r->length_B) * 8 * sizeof(bigint_word_t) + GET_FBS(a) - rfbs - 1;
 720                 /*
 721                 if((a->wordv[a->length_B-1] & ((1LL<<GET_FBS(a)) - 1)) > r->wordv[r->length_B-1]){
 722                         // cli_putc('~');
 723                         cli_putstr("\r\n ~ [a] = ");
 724                         cli_hexdump_rev(&a->wordv[a->length_B-1], 4);
 725                         cli_putstr("  [r] = ");
 726                         cli_hexdump_rev(&r->wordv[r->length_B-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_B, 2);
 732 //              cli_putstr(" r_len = "); cli_hexdump_rev(&r->length_B, 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_B == r->length_B)){
 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_B==0 || r->length_B==0){
 758                 return;
 759         }
 760
 761         bigint_t res, base;
 762         bigint_word_t t, base_b[MAX(a->length_B,r->length_B)], res_b[r->length_B*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_B=1;
 776         res.info = 0;
 777         bigint_adjust(&res);
 778         if(exp->length_B == 0){
 779                 bigint_copy(dest, &res);
 780                 return;
 781         }
 782         uint8_t flag = 0;
 783         t=exp->wordv[exp->length_B - 1];
 784         for(i=exp->length_B; 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_B==0 || y->length_B==0){
 821                  return;
 822          }
 823          if(x->length_B==1 && x->wordv[0]==1){
 824                  gcd->length_B = 1;
 825                  gcd->wordv[0] = 1;
 826                  if(a){
 827                          a->length_B = 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_B==1 && y->wordv[0]==1){
 838                  gcd->length_B = 1;
 839                  gcd->wordv[0] = 1;
 840                  if(b){
 841                          b->length_B = 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_B-i], y_b[y->length_B-i];
 856          bigint_word_t u_b[x->length_B-i], v_b[y->length_B-i];
 857          bigint_word_t a_b[y->length_B+2], c_b[y->length_B+2];
 858          bigint_word_t b_b[x->length_B+2], d_b[x->length_B+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_B = i+1;
 866          g.info=0;
 867          x_.info = y_.info = 0;
 868          x_.length_B = x->length_B-i;
 869          y_.length_B = y->length_B-i;
 870          memcpy(x_.wordv, x->wordv+i, x_.length_B*sizeof(bigint_word_t));
 871          memcpy(y_.wordv, y->wordv+i, y_.length_B*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_B = 1;
 895          a_.info = 0;
 896          d_.wordv[0] = 1;
 897          d_.length_B = 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_B);
 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 = ((uint8_t*)p)+a->length_B*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