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