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