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