bigint/bigint.c

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