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 || 1
  39 #include "cli.h"
  40 #include "uart_i.h"
  41 #include "bigint_io.h"
  42 #include <stdio.h>
  43 #endif
  44
  45 #ifndef MAX
  46  #define MAX(a,b) (((a)>(b))?(a):(b))
  47 #endif
  48
  49 #ifndef MIN
  50  #define MIN(a,b) (((a) < (b)) ? (a) : (b))
  51 #endif
  52
  53 #define SET_FBS(a, v) do{(a)->info &= ~BIGINT_FBS_MASK; (a)->info |= (v);} while(0)
  54 #define GET_FBS(a)   ((a)->info&BIGINT_FBS_MASK)
  55 #define SET_NEG(a)   (a)->info |= BIGINT_NEG_MASK
  56 #define SET_POS(a)   (a)->info &= ~BIGINT_NEG_MASK
  57 #define XCHG(a,b)    do{(a) ^= (b); (b) ^= (a); (a) ^= (b);} while(0)
  58 #define XCHG_PTR(a,b)    do{ a = (void*)(((intptr_t)(a)) ^ ((intptr_t)(b))); \
  59                                  b = (void*)(((intptr_t)(a)) ^ ((intptr_t)(b))); \
  60                                  a = (void*)(((intptr_t)(a)) ^ ((intptr_t)(b)));} while(0)
  61
  62 #define GET_SIGN(a) ((a)->info & BIGINT_NEG_MASK)
  63
  64 /******************************************************************************/
  65 void bigint_adjust(bigint_t *a){
  66         while(a->length_W!=0 && a->wordv[a->length_W-1]==0){
  67                 a->length_W--;
  68         }
  69         if(a->length_W==0){
  70                 a->info=0;
  71                 return;
  72         }
  73         bigint_word_t t;
  74         uint8_t i = BIGINT_WORD_SIZE-1;
  75         t = a->wordv[a->length_W-1];
  76         while((t&(1L<<(BIGINT_WORD_SIZE-1)))==0 && i){
  77                 t<<=1;
  78                 i--;
  79         }
  80         SET_FBS(a, i);
  81 }
  82
  83 /******************************************************************************/
  84
  85 bigint_length_t bigint_length_b(const bigint_t *a){
  86         if(!a->length_W || a->length_W==0){
  87                 return 0;
  88         }
  89         return (a->length_W-1) * BIGINT_WORD_SIZE + GET_FBS(a);
  90 }
  91
  92 /******************************************************************************/
  93
  94 bigint_length_t bigint_length_B(const bigint_t *a){
  95         return a->length_W * sizeof(bigint_word_t);
  96 }
  97
  98 /******************************************************************************/
  99
 100 uint32_t bigint_get_first_set_bit(const bigint_t *a){
 101         if(a->length_W==0){
 102                 return (uint32_t)(-1);
 103         }
 104         return (a->length_W-1)*sizeof(bigint_word_t)*8+GET_FBS(a);
 105 }
 106
 107
 108 /******************************************************************************/
 109
 110 uint32_t bigint_get_last_set_bit(const bigint_t *a){
 111         uint32_t r=0;
 112         uint8_t b=0;
 113         bigint_word_t x=1;
 114         if(a->length_W==0){
 115                 return (uint32_t)(-1);
 116         }
 117         while(a->wordv[r]==0 && r<a->length_W){
 118                 ++r;
 119         }
 120         if(a->wordv[r] == 0){
 121                 return (uint32_t)(-1);
 122         }
 123         while((x&a->wordv[r])==0){
 124                 ++b;
 125                 x <<= 1;
 126         }
 127         return r*BIGINT_WORD_SIZE+b;
 128 }
 129
 130 /******************************************************************************/
 131
 132 void bigint_copy(bigint_t *dest, const bigint_t *src){
 133     if(dest->wordv != src->wordv){
 134             memcpy(dest->wordv, src->wordv, src->length_W * sizeof(bigint_word_t));
 135     }
 136     dest->length_W = src->length_W;
 137         dest->info = src->info;
 138 }
 139
 140 /******************************************************************************/
 141
 142 /* this should be implemented in assembly */
 143 void bigint_add_u(bigint_t *dest, const bigint_t *a, const bigint_t *b){
 144         bigint_length_t i;
 145         bigint_wordplus_t t = 0LL;
 146         if (a->length_W < b->length_W) {
 147                 XCHG_PTR(a, b);
 148         }
 149         for(i = 0; i < b->length_W; ++i) {
 150                 t += a->wordv[i];
 151                 t += b->wordv[i];
 152                 dest->wordv[i] = (bigint_word_t)t;
 153                 t >>= BIGINT_WORD_SIZE;
 154         }
 155         for(; i < a->length_W; ++i){
 156                 t += a->wordv[i];
 157                 dest->wordv[i] = (bigint_word_t)t;
 158                 t >>= BIGINT_WORD_SIZE;
 159         }
 160         if(t){
 161                 dest->wordv[i++] = (bigint_word_t)t;
 162         }
 163         dest->length_W = i;
 164         SET_POS(dest);
 165         bigint_adjust(dest);
 166 }
 167
 168 /******************************************************************************/
 169
 170 /* this should be implemented in assembly */
 171 void bigint_add_scale_u(bigint_t *dest, const bigint_t *a, bigint_length_t scale){
 172         if(a->length_W == 0){
 173                 return;
 174         }
 175         if(scale == 0){
 176                 bigint_add_u(dest, dest, a);
 177                 return;
 178         }
 179         bigint_t x;
 180 #if BIGINT_WORD_SIZE == 8
 181         memset(dest->wordv + dest->length_W, 0, MAX(dest->length_W, a->length_W + scale) - dest->length_W);
 182         x.wordv = dest->wordv + scale;
 183         x.length_W = dest->length_W - scale;
 184         if((int16_t)x.length_W < 0){
 185                 x.length_W = 0;
 186                 x.info = 0;
 187         } else {
 188                 x.info = dest->info;
 189         }
 190         bigint_add_u(&x, &x, a);
 191         dest->length_W = x.length_W + scale;
 192         dest->info = 0;
 193         bigint_adjust(dest);
 194 #else
 195         bigint_t s;
 196         bigint_length_t word_shift = scale / sizeof(bigint_word_t), byte_shift = scale % sizeof(bigint_word_t);
 197         bigint_word_t bv[a->length_W + 1];
 198         s.wordv = bv;
 199         bv[0] = bv[a->length_W] = 0;
 200         memcpy((uint8_t*)bv + byte_shift, a->wordv, a->length_W * sizeof(bigint_word_t));
 201         s.length_W = a->length_W + 1;
 202         bigint_adjust(&s);
 203         memset(dest->wordv + dest->length_W, 0, (MAX(dest->length_W, s.length_W + word_shift) - dest->length_W) * sizeof(bigint_word_t));
 204         x.wordv = dest->wordv + word_shift;
 205         x.length_W = dest->length_W - word_shift;
 206         if((int16_t)x.length_W < 0){
 207                 x.length_W = 0;
 208                 x.info = 0;
 209         }else{
 210                 x.info = dest->info;
 211         }
 212         bigint_add_u(&x, &x, &s);
 213         dest->length_W = x.length_W + word_shift;
 214         dest->info = 0;
 215         bigint_adjust(dest);
 216 #endif
 217 }
 218
 219 /******************************************************************************/
 220
 221 /* this should be implemented in assembly */
 222 void bigint_sub_u(bigint_t *dest, const bigint_t *a, const bigint_t *b){
 223         int8_t borrow = 0;
 224         int8_t  r;
 225         bigint_wordplus_signed_t t = 0LL;
 226         bigint_length_t i;
 227         if(b->length_W == 0){
 228                 bigint_copy(dest, a);
 229                 SET_POS(dest);
 230                 return;
 231         }
 232         if(a->length_W == 0){
 233                 bigint_copy(dest, b);
 234                 SET_NEG(dest);
 235                 return;
 236         }
 237     r = bigint_cmp_u(a,b);
 238     if(r == 0){
 239         bigint_set_zero(dest);
 240         return;
 241     }
 242         if(r < 0){
 243                 bigint_sub_u(dest, b, a);
 244                 SET_NEG(dest);
 245                 return;
 246         }
 247         for(i = 0; i < a->length_W; ++i){
 248                 t = a->wordv[i];
 249                 if(i < b->length_W){
 250                         t -= b->wordv[i];
 251                 }
 252                 t -= borrow;
 253                 dest->wordv[i] = (bigint_word_t)t;
 254                 if(t<0){
 255                         borrow = 1;
 256                 }else{
 257                         borrow = 0;
 258                 }
 259         }
 260         SET_POS(dest);
 261         dest->length_W = i;
 262         bigint_adjust(dest);
 263 }
 264
 265 /******************************************************************************/
 266
 267 int8_t bigint_cmp_u(const bigint_t *a, const bigint_t *b){
 268         if(a->length_W > b->length_W){
 269                 return 1;
 270         }
 271         if(a->length_W < b->length_W){
 272                 return -1;
 273         }
 274         if(a->length_W == 0){
 275                 return 0;
 276         }
 277         bigint_length_t i;
 278         i = a->length_W - 1;
 279         do{
 280                 if(a->wordv[i] != b->wordv[i]){
 281                         if(a->wordv[i] > b->wordv[i]){
 282                                 return 1;
 283                         }else{
 284                                 return -1;
 285                         }
 286                 }
 287         }while(i--);
 288         return 0;
 289 }
 290
 291 /******************************************************************************/
 292
 293 void bigint_add_s(bigint_t *dest, const bigint_t *a, const bigint_t *b){
 294         uint8_t s;
 295         s  = GET_SIGN(a)?2:0;
 296         s |= GET_SIGN(b)?1:0;
 297         switch(s){
 298                 case 0: /* both positive */
 299                         bigint_add_u(dest, a,b);
 300                         SET_POS(dest);
 301                         break;
 302                 case 1: /* a positive, b negative */
 303                         bigint_sub_u(dest, a, b);
 304                         break;
 305                 case 2: /* a negative, b positive */
 306                         bigint_sub_u(dest, b, a);
 307                         break;
 308                 case 3: /* both negative */
 309                         bigint_add_u(dest, a, b);
 310                         SET_NEG(dest);
 311                         break;
 312                 default: /* how can this happen?*/
 313                         break;
 314         }
 315 }
 316
 317 /******************************************************************************/
 318
 319 void bigint_sub_s(bigint_t *dest, const bigint_t *a, const bigint_t *b){
 320         uint8_t s;
 321         s  = GET_SIGN(a)?2:0;
 322         s |= GET_SIGN(b)?1:0;
 323         switch(s){
 324                 case 0: /* both positive */
 325                         bigint_sub_u(dest, a,b);
 326                         break;
 327                 case 1: /* a positive, b negative */
 328                         bigint_add_u(dest, a, b);
 329                         SET_POS(dest);
 330                         break;
 331                 case 2: /* a negative, b positive */
 332                         bigint_add_u(dest, a, b);
 333                         SET_NEG(dest);
 334                         break;
 335                 case 3: /* both negative */
 336                         bigint_sub_u(dest, b, a);
 337                         break;
 338                 default: /* how can this happen?*/
 339                                         break;
 340         }
 341
 342 }
 343
 344 /******************************************************************************/
 345
 346 int8_t bigint_cmp_s(const bigint_t *a, const bigint_t *b){
 347         uint8_t s;
 348         if(a->length_W==0 && b->length_W==0){
 349                 return 0;
 350         }
 351         s  = GET_SIGN(a)?2:0;
 352         s |= GET_SIGN(b)?1:0;
 353         switch(s){
 354                 case 0: /* both positive */
 355                         return bigint_cmp_u(a, b);
 356                         break;
 357                 case 1: /* a positive, b negative */
 358                         return 1;
 359                         break;
 360                 case 2: /* a negative, b positive */
 361                         return -1;
 362                         break;
 363                 case 3: /* both negative */
 364                         return bigint_cmp_u(b, a);
 365                         break;
 366                 default: /* how can this happen?*/
 367                                         break;
 368         }
 369         return 0; /* just to satisfy the compiler */
 370 }
 371
 372 /******************************************************************************/
 373
 374 void bigint_shiftleft(bigint_t *a, bigint_length_t shift){
 375         bigint_length_t byteshift, words_to_shift;
 376         int16_t i;
 377         uint8_t bitshift;
 378         bigint_word_t *p;
 379         bigint_wordplus_t t = 0;
 380         if(shift == 0){
 381                 return;
 382         }
 383         byteshift = shift / 8;
 384         bitshift = shift & 7;
 385
 386         if(byteshift){
 387                 memmove(((uint8_t*)a->wordv) + byteshift, a->wordv, a->length_W * sizeof(bigint_word_t));
 388                 memset(a->wordv, 0, byteshift);
 389         }
 390         if(bitshift == 0){
 391             a->length_W += (byteshift + sizeof(bigint_word_t) - 1) / sizeof(bigint_word_t);
 392             bigint_adjust(a);
 393             return;
 394         }
 395         p = a->wordv + byteshift / sizeof(bigint_word_t);
 396         words_to_shift = a->length_W + (byteshift % sizeof(bigint_word_t)?1:0);
 397     for(i=0; i < words_to_shift; ++i){
 398         t |= ((bigint_wordplus_t)p[i]) << bitshift;
 399         p[i] = (bigint_word_t)t;
 400         t >>= BIGINT_WORD_SIZE;
 401     }
 402     if(t){
 403         p[i] = (bigint_word_t)t;
 404         a->length_W += 1;
 405     }
 406     a->length_W += (byteshift + sizeof(bigint_word_t) - 1) / sizeof(bigint_word_t);
 407         bigint_adjust(a);
 408 }
 409
 410 /******************************************************************************/
 411
 412 void bigint_shiftright(bigint_t *a, bigint_length_t shift){
 413         bigint_length_t byteshift;
 414         bigint_length_t i;
 415         uint8_t bitshift;
 416         bigint_wordplus_t t = 0;
 417         byteshift = shift / 8;
 418         bitshift = shift & 7;
 419
 420         if(byteshift >= a->length_W * sizeof(bigint_word_t)){ /* we would shift out more than we have */
 421                 bigint_set_zero(a);
 422                 return;
 423         }
 424         if(byteshift == a->length_W * sizeof(bigint_word_t) - 1 && bitshift > GET_FBS(a)){
 425                 bigint_set_zero(a);
 426                 return;
 427         }
 428
 429         if(byteshift){
 430                 memmove(a->wordv, (uint8_t*)a->wordv + byteshift, a->length_W * sizeof(bigint_word_t) - byteshift);
 431         }
 432
 433     byteshift /= sizeof(bigint_word_t); /* byteshift is now wordshift */
 434     a->length_W -= byteshift;
 435         if(bitshift != 0 && a->length_W){
 436          /* shift to the right */
 437                 i = a->length_W - 1;
 438                 do{
 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                 }while(i--);
 443         }
 444         bigint_adjust(a);
 445 }
 446
 447 /******************************************************************************/
 448
 449 void bigint_xor(bigint_t *dest, const bigint_t *a){
 450         bigint_length_t i;
 451         for(i=0; i<a->length_W; ++i){
 452                 dest->wordv[i] ^= a->wordv[i];
 453         }
 454         bigint_adjust(dest);
 455 }
 456
 457 /******************************************************************************/
 458
 459 void bigint_set_zero(bigint_t *a){
 460         a->length_W=0;
 461 }
 462
 463 /******************************************************************************/
 464
 465 /* using the Karatsuba-Algorithm */
 466 /* x*y = (xh*yh)*b**2n + ((xh+xl)*(yh+yl) - xh*yh - xl*yl)*b**n + yh*yl */
 467 void bigint_mul_u(bigint_t *dest, const bigint_t *a, const bigint_t *b){
 468         if(a->length_W == 0 || b->length_W == 0){
 469                 bigint_set_zero(dest);
 470                 return;
 471         }
 472         if(dest == a || dest == b){
 473                 bigint_t d;
 474                 bigint_word_t d_b[a->length_W + b->length_W];
 475                 d.wordv = d_b;
 476                 bigint_mul_u(&d, a, b);
 477                 bigint_copy(dest, &d);
 478                 return;
 479         }
 480         if(a->length_W == 1 || b->length_W == 1){
 481                 if(a->length_W != 1){
 482                         XCHG_PTR(a,b);
 483                 }
 484                 bigint_wordplus_t t = 0;
 485                 bigint_length_t i;
 486                 bigint_word_t x = a->wordv[0];
 487                 for(i=0; i < b->length_W; ++i){
 488                         t += ((bigint_wordplus_t)b->wordv[i]) * ((bigint_wordplus_t)x);
 489                         dest->wordv[i] = (bigint_word_t)t;
 490                         t >>= BIGINT_WORD_SIZE;
 491                 }
 492                 dest->length_W = i;
 493                 if(t){
 494                     dest->wordv[i] = (bigint_word_t)t;
 495                     dest->length_W += 1;
 496                 }
 497                 dest->info = 0;
 498                 bigint_adjust(dest);
 499                 return;
 500         }
 501         if(a->length_W * sizeof(bigint_word_t) <= 4 && b->length_W * sizeof(bigint_word_t) <= 4){
 502                 uint32_t p=0, q=0;
 503                 uint64_t r;
 504                 memcpy(&p, a->wordv, a->length_W*sizeof(bigint_word_t));
 505                 memcpy(&q, b->wordv, b->length_W*sizeof(bigint_word_t));
 506                 r = (uint64_t)p * (uint64_t)q;
 507                 memcpy(dest->wordv, &r, (dest->length_W = a->length_W + b->length_W)*sizeof(bigint_word_t));
 508                 bigint_adjust(dest);
 509                 return;
 510         }
 511         /* split a in xh & xl; split b in yh & yl */
 512         const bigint_length_t n = (MAX(a->length_W, b->length_W)+1)/2;
 513         bigint_t xl, xh, yl, yh;
 514         xl.wordv = a->wordv;
 515         yl.wordv = b->wordv;
 516         if(a->length_W<=n){
 517                 bigint_set_zero(&xh);
 518                 xl.length_W = a->length_W;
 519                 xl.info = a->info;
 520         }else{
 521                 xl.length_W=n;
 522                 xl.info = 0;
 523                 bigint_adjust(&xl);
 524                 xh.wordv = &(a->wordv[n]);
 525                 xh.length_W = a->length_W-n;
 526                 xh.info = a->info;
 527         }
 528         if(b->length_W<=n){
 529                 bigint_set_zero(&yh);
 530                 yl.length_W = b->length_W;
 531                 yl.info = b->info;
 532         }else{
 533                 yl.length_W=n;
 534                 yl.info = 0;
 535                 bigint_adjust(&yl);
 536                 yh.wordv = &(b->wordv[n]);
 537                 yh.length_W = b->length_W-n;
 538                 yh.info = b->info;
 539         }
 540         /* now we have split up a and b */
 541         /* remember we want to do:
 542          * x*y = (xh * b**n + xl) * (yh * b**n + yl)
 543          *     = (xh * yh) * b**2n + xh * b**n * yl + yh * b**n * xl + xl * yl
 544          *     = (xh * yh) * b**2n + (xh * yl + yh * xl) * b**n + xl *yl
 545          *     // xh * yl + yh * xl = (xh + yh) * (xl + yl) - xh * yh - xl * yl
 546          * x*y = (xh * yh) * b**2n + ((xh+xl)*(yh+yl) - xh*yh - xl*yl)*b**n + xl*yl
 547          *          5              9     2   4   3    7   5   6   1         8   1
 548          */
 549         bigint_word_t  tmp_b[2 * n + 2], m_b[2 * (n + 1)];
 550         bigint_t tmp, tmp2, m;
 551         tmp.wordv = tmp_b;
 552         tmp2.wordv = &(tmp_b[n + 1]);
 553         m.wordv = m_b;
 554
 555         bigint_mul_u(dest, &xl, &yl);  /* 1: dest <= xl*yl     */
 556         bigint_add_u(&tmp2, &xh, &xl); /* 2: tmp2 <= xh+xl     */
 557         bigint_add_u(&tmp, &yh, &yl);  /* 3: tmp  <= yh+yl     */
 558         bigint_mul_u(&m, &tmp2, &tmp); /* 4: m    <= tmp2*tmp  */
 559         bigint_mul_u(&tmp, &xh, &yh);  /* 5: tmp  <= xh*yh     */
 560         bigint_sub_u(&m, &m, dest);    /* 6: m    <= m-dest    */
 561     bigint_sub_u(&m, &m, &tmp);    /* 7: m    <= m-tmp     */
 562         bigint_add_scale_u(dest, &m, n * sizeof(bigint_word_t));       /* 8: dest <= dest+m**n*/
 563         bigint_add_scale_u(dest, &tmp, 2 * n * sizeof(bigint_word_t)); /* 9: dest <= dest+tmp**(2*n) */
 564 }
 565
 566 /******************************************************************************/
 567
 568 void bigint_mul_s(bigint_t *dest, const bigint_t *a, const bigint_t *b){
 569         uint8_t s;
 570         s  = GET_SIGN(a)?2:0;
 571         s |= GET_SIGN(b)?1:0;
 572         switch(s){
 573                 case 0: /* both positive */
 574                         bigint_mul_u(dest, a,b);
 575                         SET_POS(dest);
 576                         break;
 577                 case 1: /* a positive, b negative */
 578                         bigint_mul_u(dest, a,b);
 579                         SET_NEG(dest);
 580                         break;
 581                 case 2: /* a negative, b positive */
 582                         bigint_mul_u(dest, a,b);
 583                         SET_NEG(dest);
 584                         break;
 585                 case 3: /* both negative */
 586                         bigint_mul_u(dest, a,b);
 587                         SET_POS(dest);
 588                         break;
 589                 default: /* how can this happen?*/
 590                         break;
 591         }
 592 }
 593
 594 /******************************************************************************/
 595 #if OLD_SQUARE
 596
 597 #if DEBUG_SQUARE
 598 unsigned square_depth = 0;
 599 #endif
 600
 601 /* square */
 602 /* (xh*b^n+xl)^2 = xh^2*b^2n + 2*xh*xl*b^n + xl^2 */
 603 void bigint_square(bigint_t *dest, const bigint_t *a){
 604     if(a->length_W * sizeof(bigint_word_t) <= 4){
 605                 uint64_t r = 0;
 606                 memcpy(&r, a->wordv, a->length_W * sizeof(bigint_word_t));
 607                 r = r * r;
 608                 memcpy(dest->wordv, &r, 2 * a->length_W * sizeof(bigint_word_t));
 609                 SET_POS(dest);
 610                 dest->length_W = 2 * a->length_W;
 611                 bigint_adjust(dest);
 612                 return;
 613         }
 614
 615         if(dest->wordv == a->wordv){
 616                 bigint_t d;
 617                 bigint_word_t d_b[a->length_W*2];
 618                 d.wordv = d_b;
 619                 bigint_square(&d, a);
 620                 bigint_copy(dest, &d);
 621                 return;
 622         }
 623         bigint_fast_square(dest, a);
 624         return;
 625 #if DEBUG_SQUARE
 626         square_depth += 1;
 627 #endif
 628
 629         bigint_length_t n;
 630         n=(a->length_W+1)/2;
 631         bigint_t xh, xl, tmp; /* x-high, x-low, temp */
 632         bigint_word_t buffer[2*n+1];
 633         xl.wordv = a->wordv;
 634         xl.length_W = n;
 635         xl.info = 0;
 636         xh.wordv = &(a->wordv[n]);
 637         xh.length_W = a->length_W-n;
 638         xh.info = a->info;
 639         bigint_adjust(&xl);
 640         tmp.wordv = buffer;
 641 /* (xh * b**n + xl)**2 = xh**2 * b**2n + 2 * xh * xl * b**n + xl**2 */
 642 #if DEBUG_SQUARE
 643         if(square_depth == 1){
 644         cli_putstr("\r\nDBG (a): xl: "); bigint_print_hex(&xl);
 645         cli_putstr("\r\nDBG (b): xh: "); bigint_print_hex(&xh);
 646         }
 647 #endif
 648         bigint_square(dest, &xl);
 649 #if DEBUG_SQUARE
 650         if(square_depth == 1){
 651             cli_putstr("\r\nDBG (1): xl**2: "); bigint_print_hex(dest);
 652         }
 653 #endif
 654     bigint_square(&tmp, &xh);
 655 #if DEBUG_SQUARE
 656     if(square_depth == 1){
 657         cli_putstr("\r\nDBG (2): xh**2: "); bigint_print_hex(&tmp);
 658     }
 659 #endif
 660         bigint_add_scale_u(dest, &tmp, 2 * n * sizeof(bigint_word_t));
 661 #if DEBUG_SQUARE
 662         if(square_depth == 1){
 663             cli_putstr("\r\nDBG (3): xl**2 + xh**2*n**2: "); bigint_print_hex(dest);
 664         }
 665 #endif
 666         bigint_mul_u(&tmp, &xl, &xh);
 667 #if DEBUG_SQUARE
 668         if(square_depth == 1){
 669             cli_putstr("\r\nDBG (4): xl*xh: "); bigint_print_hex(&tmp);
 670         }
 671 #endif
 672         bigint_shiftleft(&tmp, 1);
 673 #if DEBUG_SQUARE
 674         if(square_depth == 1){
 675             cli_putstr("\r\nDBG (5): xl*xh*2: "); bigint_print_hex(&tmp);
 676         }
 677 #endif
 678         bigint_add_scale_u(dest, &tmp, n * sizeof(bigint_word_t));
 679 #if DEBUG_SQUARE
 680         if(square_depth == 1){
 681             cli_putstr("\r\nDBG (6): x**2: "); bigint_print_hex(dest);
 682             cli_putstr("\r\n");
 683         }
 684         square_depth -= 1;
 685 #endif
 686 }
 687
 688 #endif
 689
 690 /******************************************************************************/
 691
 692 void bigint_square(bigint_t *dest, const bigint_t *a) {
 693     union {
 694         bigint_word_t u[2];
 695         bigint_wordplus_t uv;
 696     } acc;
 697     bigint_word_t q, c1, c2;
 698     bigint_length_t i, j;
 699
 700     if(a->length_W * sizeof(bigint_word_t) <= 4){
 701         uint64_t r = 0;
 702         memcpy(&r, a->wordv, a->length_W * sizeof(bigint_word_t));
 703         r = r * r;
 704         memcpy(dest->wordv, &r, 2 * a->length_W * sizeof(bigint_word_t));
 705         SET_POS(dest);
 706         dest->length_W = 2 * a->length_W;
 707         bigint_adjust(dest);
 708         return;
 709     }
 710
 711     if(dest->wordv == a->wordv){
 712         bigint_t d;
 713         bigint_word_t d_b[a->length_W*2];
 714         d.wordv = d_b;
 715         bigint_square(&d, a);
 716         bigint_copy(dest, &d);
 717         return;
 718     }
 719
 720     memset(dest->wordv, 0, a->length_W * 2 * sizeof(bigint_word_t));
 721
 722     for (i = 0; i < a->length_W; ++i) {
 723         acc.uv = a->wordv[i] * a->wordv[i] + dest->wordv[2 * i];
 724         dest->wordv[2 * i] = acc.u[0];
 725         c1 = acc.u[1];
 726         c2 = 0;
 727         for (j = i + 1; j < a->length_W; ++j) {
 728             acc.uv = a->wordv[i] * a->wordv[j];
 729             q = acc.u[1] >> (BIGINT_WORD_SIZE - 1);
 730             acc.uv <<= 1;
 731             acc.uv += dest->wordv[i + j];
 732             q += (acc.uv < dest->wordv[i + j]);
 733             acc.uv += c1;
 734             q += (acc.uv < c1);
 735             dest->wordv[i + j] = acc.u[0];
 736             c1 = acc.u[1] + c2;
 737             c2 = q + (c1 < c2);
 738         }
 739         dest->wordv[i + a->length_W] += c1;
 740         dest->wordv[i + a->length_W + 1] = c2 + (dest->wordv[i + a->length_W] < c1);
 741     }
 742     SET_POS(dest);
 743     dest->length_W = 2 * a->length_W;
 744     bigint_adjust(dest);
 745 }
 746
 747 /******************************************************************************/
 748
 749 void bigint_sub_u_bitscale(bigint_t *a, const bigint_t *b, bigint_length_t bitscale){
 750         bigint_t tmp, x;
 751         bigint_word_t tmp_b[b->length_W + 1];
 752         const bigint_length_t word_shift = bitscale / BIGINT_WORD_SIZE;
 753
 754         if(a->length_W < b->length_W + word_shift){
 755 #if DEBUG
 756                 cli_putstr("\r\nDBG: *bang*\r\n");
 757 #endif
 758                 bigint_set_zero(a);
 759                 return;
 760         }
 761         tmp.wordv = tmp_b;
 762         bigint_copy(&tmp, b);
 763         bigint_shiftleft(&tmp, bitscale % BIGINT_WORD_SIZE);
 764
 765         x.info = a->info;
 766         x.wordv = &(a->wordv[word_shift]);
 767         x.length_W = a->length_W - word_shift;
 768
 769         bigint_sub_u(&x, &x, &tmp);
 770         bigint_adjust(a);
 771         return;
 772 }
 773
 774 /******************************************************************************/
 775
 776 void bigint_reduce(bigint_t *a, const bigint_t *r){
 777         uint8_t rfbs = GET_FBS(r);
 778         if(r->length_W==0 || a->length_W==0){
 779                 return;
 780         }
 781
 782         if(bigint_length_b(a) + 3 > bigint_length_b(r)){
 783         if((r->length_W*sizeof(bigint_word_t)<=4) && (a->length_W*sizeof(bigint_word_t)<=4)){
 784             uint32_t p=0, q=0;
 785             memcpy(&p, a->wordv, a->length_W*sizeof(bigint_word_t));
 786             memcpy(&q, r->wordv, r->length_W*sizeof(bigint_word_t));
 787             p %= q;
 788             memcpy(a->wordv, &p, a->length_W*sizeof(bigint_word_t));
 789             bigint_adjust(a);
 790             return;
 791         }
 792         bigint_length_t shift;
 793         while(a->length_W > r->length_W){
 794             shift = (a->length_W - r->length_W) * 8 * sizeof(bigint_word_t) + GET_FBS(a) - rfbs - 1;
 795             bigint_sub_u_bitscale(a, r, shift);
 796         }
 797         while((GET_FBS(a) > rfbs) && (a->length_W == r->length_W)){
 798             shift = GET_FBS(a)-rfbs-1;
 799             bigint_sub_u_bitscale(a, r, shift);
 800         }
 801         }
 802         while(bigint_cmp_u(a,r)>=0){
 803                 bigint_sub_u(a,a,r);
 804         }
 805         bigint_adjust(a);
 806 }
 807
 808 /******************************************************************************/
 809
 810 /* calculate dest = a**exp % r */
 811 /* using square&multiply */
 812 void bigint_expmod_u_sam(bigint_t *dest, const bigint_t *a, const bigint_t *exp, const bigint_t *r){
 813         if(a->length_W == 0){
 814             bigint_set_zero(dest);
 815                 return;
 816         }
 817
 818         bigint_t res, base;
 819         bigint_word_t t, base_w[MAX(a->length_W,r->length_W)], res_w[r->length_W*2];
 820         bigint_length_t i;
 821         uint8_t j;
 822 //      bigint_length_t *xaddr = &i;
 823 //      cli_putstr("\r\npre-alloc (");
 824 //      cli_hexdump_rev(&xaddr, 4);
 825 //      cli_putstr(") ...");
 826         res.wordv = res_w;
 827         base.wordv = base_w;
 828         bigint_copy(&base, a);
 829 //      cli_putstr("\r\npost-copy");
 830         bigint_reduce(&base, r);
 831         res.wordv[0] = 1;
 832         res.length_W = 1;
 833         res.info = 0;
 834         bigint_adjust(&res);
 835         if(exp->length_W == 0){
 836                 bigint_copy(dest, &res);
 837                 return;
 838         }
 839         uint8_t flag = 0;
 840         t=exp->wordv[exp->length_W - 1];
 841         for(i = exp->length_W; i > 0; --i){
 842                 t = exp->wordv[i - 1];
 843                 for(j = BIGINT_WORD_SIZE; j > 0; --j){
 844                         if(!flag){
 845                                 if(t & (1 << (BIGINT_WORD_SIZE - 1))){
 846                                         flag = 1;
 847                                 }
 848                         }
 849                         if(flag){
 850                                 bigint_square(&res, &res);
 851                                 bigint_reduce(&res, r);
 852                                 if(t & (1 << (BIGINT_WORD_SIZE - 1))){
 853                                         bigint_mul_u(&res, &res, &base);
 854                                         bigint_reduce(&res, r);
 855                                 }
 856                         }
 857                         t <<= 1;
 858                 }
 859         }
 860
 861 //      cli_putc('+');
 862         SET_POS(&res);
 863         bigint_copy(dest, &res);
 864 }
 865
 866 /******************************************************************************/
 867 /* gcd <-- gcd(x,y) a*x+b*y=gcd */
 868 void bigint_gcdext(bigint_t *gcd, bigint_t *a, bigint_t *b, const bigint_t *x, const bigint_t *y){
 869          bigint_length_t i = 0;
 870          if(x->length_W == 0 || y->length_W == 0){
 871              if(gcd){
 872                  bigint_set_zero(gcd);
 873              }
 874              if(a){
 875                  bigint_set_zero(a);
 876              }
 877          if(b){
 878              bigint_set_zero(b);
 879          }
 880                  return;
 881          }
 882          if(x->length_W == 1 && x->wordv[0] == 1){
 883              if(gcd){
 884              gcd->length_W = 1;
 885              gcd->wordv[0] = 1;
 886              gcd->info = 0;
 887              }
 888                  if(a){
 889                          a->length_W = 1;
 890                          a->wordv[0] = 1;
 891                          SET_POS(a);
 892                          bigint_adjust(a);
 893                  }
 894                  if(b){
 895                          bigint_set_zero(b);
 896                  }
 897                  return;
 898          }
 899          if(y->length_W == 1 && y->wordv[0] == 1){
 900                  if(gcd){
 901              gcd->length_W = 1;
 902              gcd->wordv[0] = 1;
 903              gcd->info = 0;
 904                  }
 905                  if(b){
 906                          b->length_W = 1;
 907                          b->wordv[0] = 1;
 908                          SET_POS(b);
 909                          bigint_adjust(b);
 910                  }
 911                  if(a){
 912                          bigint_set_zero(a);
 913                  }
 914                  return;
 915          }
 916
 917          while(x->wordv[i] == 0 && y->wordv[i] == 0){
 918                  ++i;
 919          }
 920          bigint_word_t g_b[i + 2], x_b[x->length_W - i], y_b[y->length_W - i];
 921          bigint_word_t u_b[x->length_W - i], v_b[y->length_W - i];
 922          bigint_word_t a_b[y->length_W + 2], c_b[y->length_W + 2];
 923          bigint_word_t b_b[x->length_W + 2], d_b[x->length_W + 2];
 924      bigint_t g, x_, y_, u, v, a_, b_, c_, d_;
 925
 926          g.wordv = g_b;
 927          x_.wordv = x_b;
 928          y_.wordv = y_b;
 929          memset(g_b, 0, i * sizeof(bigint_word_t));
 930          g_b[i] = 1;
 931          g.length_W = i + 1;
 932          g.info = 0;
 933          x_.info = y_.info = 0;
 934          x_.length_W = x->length_W - i;
 935          y_.length_W = y->length_W - i;
 936          memcpy(x_.wordv, x->wordv + i, x_.length_W * sizeof(bigint_word_t));
 937          memcpy(y_.wordv, y->wordv + i, y_.length_W * sizeof(bigint_word_t));
 938          for(i = 0; (x_.wordv[0] & (1 << i)) == 0 && (y_.wordv[0] & (1 << i)) == 0; ++i){
 939          }
 940
 941          bigint_adjust(&x_);
 942          bigint_adjust(&y_);
 943
 944          if(i){
 945                  bigint_shiftleft(&g, i);
 946                  bigint_shiftright(&x_, i);
 947                  bigint_shiftright(&y_, i);
 948          }
 949
 950          u.wordv = u_b;
 951          v.wordv = v_b;
 952          a_.wordv = a_b;
 953          b_.wordv = b_b;
 954          c_.wordv = c_b;
 955          d_.wordv = d_b;
 956
 957          bigint_copy(&u, &x_);
 958          bigint_copy(&v, &y_);
 959          a_.wordv[0] = 1;
 960          a_.length_W = 1;
 961          a_.info = 0;
 962          d_.wordv[0] = 1;
 963          d_.length_W = 1;
 964          d_.info = 0;
 965          bigint_set_zero(&b_);
 966          bigint_set_zero(&c_);
 967          do{
 968                  while((u.wordv[0] & 1) == 0){
 969                          bigint_shiftright(&u, 1);
 970                          if((a_.wordv[0] & 1) || (b_.wordv[0] & 1)){
 971                                  bigint_add_s(&a_, &a_, &y_);
 972                                  bigint_sub_s(&b_, &b_, &x_);
 973                          }
 974                          bigint_shiftright(&a_, 1);
 975                          bigint_shiftright(&b_, 1);
 976                  }
 977                  while((v.wordv[0] & 1) == 0){
 978                          bigint_shiftright(&v, 1);
 979                          if((c_.wordv[0] & 1) || (d_.wordv[0] & 1)){
 980                                  bigint_add_s(&c_, &c_, &y_);
 981                                  bigint_sub_s(&d_, &d_, &x_);
 982                          }
 983                          bigint_shiftright(&c_, 1);
 984                          bigint_shiftright(&d_, 1);
 985                  }
 986                  if(bigint_cmp_u(&u, &v) >= 0){
 987                         bigint_sub_u(&u, &u, &v);
 988                         bigint_sub_s(&a_, &a_, &c_);
 989                         bigint_sub_s(&b_, &b_, &d_);
 990                  }else{
 991                         bigint_sub_u(&v, &v, &u);
 992                         bigint_sub_s(&c_, &c_, &a_);
 993                         bigint_sub_s(&d_, &d_, &b_);
 994                  }
 995          }while(u.length_W);
 996          if(gcd){
 997                  bigint_mul_s(gcd, &v, &g);
 998          }
 999          if(a){
1000                 bigint_copy(a, &c_);
1001          }
1002          if(b){
1003                  bigint_copy(b, &d_);
1004          }
1005 }
1006
1007 /******************************************************************************/
1008
1009 void bigint_inverse(bigint_t *dest, const bigint_t *a, const bigint_t *m){
1010         bigint_gcdext(NULL, dest, NULL, a, m);
1011         while(dest->info&BIGINT_NEG_MASK){
1012                 bigint_add_s(dest, dest, m);
1013         }
1014 }
1015
1016 /******************************************************************************/
1017
1018 void bigint_changeendianess(bigint_t *a){
1019         uint8_t t, *p, *q;
1020         p = (uint8_t*)(a->wordv);
1021         q = p + a->length_W * sizeof(bigint_word_t) - 1;
1022         while(p<q){
1023                 t = *p;
1024                 *p = *q;
1025                 *q = t;
1026                 ++p; --q;
1027         }
1028 }
1029
1030 /******************************************************************************/
1031
1032 void bigint_mul_word_u(bigint_t *a, bigint_word_t b){
1033     bigint_wordplus_t c0 = 0, c1 = 0;
1034     bigint_length_t i;
1035
1036     if(b == 0){
1037         bigint_set_zero(a);
1038         return;
1039     }
1040
1041     for(i = 0; i < a->length_W; ++i){
1042         c1 = ((bigint_wordplus_t)(a->wordv[i])) * ((bigint_wordplus_t)b);
1043         c1 += c0;
1044         a->wordv[i] = (bigint_word_t)c1;
1045         c0 = c1 >> BIGINT_WORD_SIZE;
1046     }
1047     if(c0){
1048         a->wordv[a->length_W] = (bigint_word_t)c0;
1049         a->length_W += 1;
1050     }
1051     bigint_adjust(a);
1052 }
1053
1054 /******************************************************************************/
1055 #if 1
1056
1057 void bigint_clip(bigint_t *dest, bigint_length_t length_W){
1058     if(dest->length_W > length_W){
1059         dest->length_W = length_W;
1060     }
1061     bigint_adjust(dest);
1062 }
1063 /******************************************************************************/
1064
1065 /*
1066  * m_ = m * m'[0]
1067  */
1068
1069 void bigint_mont_mul(bigint_t *dest, const bigint_t *a, const bigint_t *b, const bigint_t *m, const bigint_t *m_){
1070     const bigint_length_t s = MAX(MAX(a->length_W, b->length_W), m->length_W);
1071     bigint_t u, t;
1072     bigint_word_t u_w[s + 2], t_w[s + 2];
1073     bigint_length_t i;
1074
1075     if (a->length_W == 0 || b->length_W == 0) {
1076         bigint_set_zero(dest);
1077         return;
1078     }
1079     u.wordv = u_w;
1080     u.info = 0;
1081     u.length_W = 0;
1082     t.wordv = t_w;
1083     for (i = 0; i < a->length_W; ++i) {
1084         bigint_copy(&t, b);
1085         bigint_mul_word_u(&t, a->wordv[i]);
1086         bigint_add_u(&u, &u, &t);
1087         bigint_copy(&t, m_);
1088         if (u.length_W != 0) {
1089             bigint_mul_word_u(&t, u.wordv[0]);
1090             bigint_add_u(&u, &u, &t);
1091         }
1092         bigint_shiftright(&u, BIGINT_WORD_SIZE);
1093     }
1094     for (; i < s; ++i) {
1095         bigint_copy(&t, m_);
1096         if (u.length_W != 0) {
1097             bigint_mul_word_u(&t, u.wordv[0]);
1098             bigint_add_u(&u, &u, &t);
1099         }
1100         bigint_shiftright(&u, BIGINT_WORD_SIZE);
1101     }
1102     bigint_reduce(&u, m);
1103     bigint_copy(dest, &u);
1104 }
1105
1106 /******************************************************************************/
1107
1108 void bigint_mont_red(bigint_t *dest, const bigint_t *a, const bigint_t *m, const bigint_t *m_){
1109     bigint_t u, t;
1110     bigint_length_t i, s = MAX(a->length_W, m->length_W);
1111     bigint_word_t u_w[s + 2], t_w[s + 2];
1112
1113     t.wordv = t_w;
1114     u.wordv = u_w;
1115     if (a->length_W == 0) {
1116         bigint_set_zero(dest);
1117         return;
1118     }
1119     bigint_copy(&u, a);
1120     for (i = 0; i < m->length_W; ++i) {
1121         bigint_copy(&t, m_);
1122         if (u.length_W != 0) {
1123             bigint_mul_word_u(&t, u.wordv[0]);
1124             bigint_add_u(&u, &u, &t);
1125         }
1126         bigint_shiftright(&u, BIGINT_WORD_SIZE);
1127     }
1128     bigint_reduce(&u, m);
1129     bigint_copy(dest, &u);
1130 }
1131
1132 /******************************************************************************/
1133 /*
1134  * m_ = m * (- m0^-1 (mod 2^W))
1135  */
1136 void bigint_mont_gen_m_(bigint_t* dest, const bigint_t* m){
1137     bigint_word_t x_w[2], m_w_0[1];
1138     bigint_t x, m_0;
1139     if (m->length_W == 0) {
1140         bigint_set_zero(dest);
1141         return;
1142     }
1143     if ((m->wordv[0] & 1) == 0) {
1144         printf_P(PSTR("ERROR: m must not be even, m = "));
1145         bigint_print_hex(m);
1146         putchar('\n');
1147         uart0_flush();
1148         return;
1149     }
1150     x.wordv = x_w;
1151     x.info = 0;
1152     x.length_W = 2;
1153     x_w[0] = 0;
1154     x_w[1] = 1;
1155     m_0.wordv = m_w_0;
1156     m_0.info = 0;
1157     m_0.length_W = 1;
1158     m_0.wordv[0] = m->wordv[0];
1159     bigint_adjust(&x);
1160     bigint_adjust(&m_0);
1161     bigint_inverse(dest, &m_0, &x);
1162     bigint_sub_s(&x, &x, dest);
1163     bigint_copy(dest, m);
1164     bigint_mul_word_u(dest, x.wordv[0]);
1165 }
1166
1167 /******************************************************************************/
1168
1169 /*
1170  * dest = a * R mod m
1171  */
1172 void bigint_mont_trans(bigint_t *dest, const bigint_t *a, const bigint_t *m){
1173     bigint_t t;
1174     bigint_word_t t_w[a->length_W + m->length_W];
1175
1176     t.wordv = t_w;
1177     memset(t_w, 0, m->length_W * sizeof(bigint_word_t));
1178     memcpy(t_w + m->length_W * sizeof(bigint_word_t), a->wordv, a->length_W * sizeof(bigint_word_t));
1179     t.info = a->info;
1180     t.length_W = a->length_W + m->length_W;
1181     bigint_reduce(&t, m);
1182     bigint_copy(dest, &t);
1183 }
1184
1185 /******************************************************************************/
1186
1187 /* calculate dest = a**exp % r */
1188 /* using square&multiply */
1189 void bigint_expmod_u_mont_sam(bigint_t *dest, const bigint_t *a, const bigint_t *exp, const bigint_t *r){
1190     if(r->length_W == 0) {
1191         return;
1192     }
1193     if(a->length_W == 0) {
1194         bigint_set_zero(dest);
1195         return;
1196     }
1197     bigint_length_t s = r->length_W;
1198     bigint_t res, m_, ax;
1199     bigint_word_t t, res_w[r->length_W*2], ax_w[MAX(s, a->length_W)];
1200     bigint_word_t m_w_[s + 1];
1201     bigint_length_t i;
1202     uint8_t j;
1203
1204     res.wordv = res_w;
1205     m_.wordv = m_w_;
1206     ax.wordv = ax_w;
1207
1208     res.wordv[0] = 1;
1209     res.length_W = 1;
1210     res.info = 0;
1211     bigint_adjust(&res);
1212     if (exp->length_W == 0) {
1213         bigint_copy(dest, &res);
1214         return;
1215     }
1216     bigint_copy(&ax, a);
1217     bigint_reduce(&ax, r);
1218     bigint_mont_trans(&ax, &ax, r);
1219     bigint_mont_trans(&res, &res, r);
1220     bigint_mont_gen_m_(&m_, r);
1221     uint8_t flag = 0;
1222     t = exp->wordv[exp->length_W - 1];
1223     for (i = exp->length_W; i > 0; --i) {
1224         t = exp->wordv[i - 1];
1225         for(j = BIGINT_WORD_SIZE; j > 0; --j){
1226             if (!flag) {
1227                 if(t & (1 << (BIGINT_WORD_SIZE - 1))){
1228                     flag = 1;
1229                 }
1230             }
1231             if (flag) {
1232                 bigint_square(&res, &res);
1233                 bigint_mont_red(&res, &res, r, &m_);
1234                 if (t & (1 << (BIGINT_WORD_SIZE - 1))) {
1235                     bigint_mont_mul(&res, &res, &ax, r, &m_);
1236                 }
1237             }
1238             t <<= 1;
1239         }
1240     }
1241     SET_POS(&res);
1242     bigint_mont_red(dest, &res, r, &m_);
1243 }
1244
1245 /******************************************************************************/
1246
1247 #endif
1248
1249 void bigint_expmod_u(bigint_t *dest, const bigint_t *a, const bigint_t *exp, const bigint_t *r){
1250     if (r->wordv[0] & 1) {
1251         bigint_expmod_u_mont_sam(dest, a, exp, r);
1252     } else {
1253         bigint_expmod_u_sam(dest, a, exp, r);
1254     }
1255 }
1256
1257
1258
1259
1260
1261
1262
1263
1264
1265
1266
1267
1268
1269
1270
1271