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