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