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 & (((bigint_word_t)1) << (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 int32_t bigint_get_first_set_bit(const bigint_t *a){
 101         if(a->length_W == 0) {
 102                 return -1;
 103         }
 104         return (a->length_W-1) * sizeof(bigint_word_t) * CHAR_BIT + GET_FBS(a);
 105 }
 106
 107
 108 /******************************************************************************/
 109
 110 int32_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 -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 = 0;
 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                 borrow = t < 0 ? 1 : 0;
 255         }
 256         SET_POS(dest);
 257         dest->length_W = i;
 258         bigint_adjust(dest);
 259 }
 260
 261 /******************************************************************************/
 262
 263 int8_t bigint_cmp_u(const bigint_t *a, const bigint_t *b){
 264         if(a->length_W > b->length_W){
 265                 return 1;
 266         }
 267         if(a->length_W < b->length_W){
 268                 return -1;
 269         }
 270         if(a->length_W == 0){
 271                 return 0;
 272         }
 273         bigint_length_t i;
 274         i = a->length_W - 1;
 275         do{
 276                 if(a->wordv[i] != b->wordv[i]){
 277                         if(a->wordv[i] > b->wordv[i]){
 278                                 return 1;
 279                         }else{
 280                                 return -1;
 281                         }
 282                 }
 283         }while(i--);
 284         return 0;
 285 }
 286
 287 /******************************************************************************/
 288
 289 void bigint_add_s(bigint_t *dest, const bigint_t *a, const bigint_t *b){
 290         uint8_t s;
 291         s  = GET_SIGN(a)?2:0;
 292         s |= GET_SIGN(b)?1:0;
 293         switch(s){
 294                 case 0: /* both positive */
 295                         bigint_add_u(dest, a,b);
 296                         SET_POS(dest);
 297                         break;
 298                 case 1: /* a positive, b negative */
 299                         bigint_sub_u(dest, a, b);
 300                         break;
 301                 case 2: /* a negative, b positive */
 302                         bigint_sub_u(dest, b, a);
 303                         break;
 304                 case 3: /* both negative */
 305                         bigint_add_u(dest, a, b);
 306                         SET_NEG(dest);
 307                         break;
 308                 default: /* how can this happen?*/
 309                         break;
 310         }
 311 }
 312
 313 /******************************************************************************/
 314
 315 void bigint_sub_s(bigint_t *dest, const bigint_t *a, const bigint_t *b){
 316         uint8_t s;
 317         s  = GET_SIGN(a)?2:0;
 318         s |= GET_SIGN(b)?1:0;
 319         switch(s){
 320                 case 0: /* both positive */
 321                         bigint_sub_u(dest, a,b);
 322                         break;
 323                 case 1: /* a positive, b negative */
 324                         bigint_add_u(dest, a, b);
 325                         SET_POS(dest);
 326                         break;
 327                 case 2: /* a negative, b positive */
 328                         bigint_add_u(dest, a, b);
 329                         SET_NEG(dest);
 330                         break;
 331                 case 3: /* both negative */
 332                         bigint_sub_u(dest, b, a);
 333                         break;
 334                 default: /* how can this happen?*/
 335                                         break;
 336         }
 337
 338 }
 339
 340 /******************************************************************************/
 341
 342 int8_t bigint_cmp_s(const bigint_t *a, const bigint_t *b){
 343         uint8_t s;
 344         if(a->length_W==0 && b->length_W==0){
 345                 return 0;
 346         }
 347         s  = GET_SIGN(a)?2:0;
 348         s |= GET_SIGN(b)?1:0;
 349         switch(s){
 350                 case 0: /* both positive */
 351                         return bigint_cmp_u(a, b);
 352                         break;
 353                 case 1: /* a positive, b negative */
 354                         return 1;
 355                         break;
 356                 case 2: /* a negative, b positive */
 357                         return -1;
 358                         break;
 359                 case 3: /* both negative */
 360                         return bigint_cmp_u(b, a);
 361                         break;
 362                 default: /* how can this happen?*/
 363                                         break;
 364         }
 365         return 0; /* just to satisfy the compiler */
 366 }
 367
 368 /******************************************************************************/
 369
 370 void bigint_shiftleft(bigint_t *a, bigint_length_t shift){
 371         bigint_length_t byteshift, words_to_shift;
 372         int16_t i;
 373         uint8_t bitshift;
 374         bigint_word_t *p;
 375         bigint_wordplus_t t = 0;
 376
 377         if (a->length_W == 0 || shift == 0) {
 378                 return;
 379         }
 380         byteshift = shift / 8;
 381         bitshift = shift & 7;
 382     if (byteshift % sizeof(bigint_word_t)) {
 383         a->wordv[a->length_W + byteshift / sizeof(bigint_t)] = 0;
 384     }
 385         if (byteshift) {
 386                 memmove(((uint8_t*)a->wordv) + byteshift, a->wordv, a->length_W * sizeof(bigint_word_t));
 387                 memset(a->wordv, 0, byteshift);
 388         }
 389         if (bitshift == 0) {
 390             a->length_W += (byteshift + sizeof(bigint_word_t) - 1) / sizeof(bigint_word_t);
 391             bigint_adjust(a);
 392             return;
 393         }
 394         p = &a->wordv[byteshift / sizeof(bigint_word_t)];
 395         words_to_shift = a->length_W + (byteshift % sizeof(bigint_word_t) ? 1 : 0);
 396     for (i = 0; i < words_to_shift; ++i) {
 397         t |= ((bigint_wordplus_t)p[i]) << bitshift;
 398         p[i] = (bigint_word_t)t;
 399         t >>= BIGINT_WORD_SIZE;
 400     }
 401     if (t) {
 402         p[i] = (bigint_word_t)t;
 403         a->length_W += 1;
 404     }
 405     a->length_W += (byteshift + sizeof(bigint_word_t) - 1) / sizeof(bigint_word_t);
 406     bigint_adjust(a);
 407 }
 408
 409 /******************************************************************************/
 410
 411 void bigint_shiftright(bigint_t *a, bigint_length_t shift){
 412         bigint_length_t byteshift;
 413         bigint_length_t i;
 414         uint8_t bitshift;
 415         bigint_wordplus_t t = 0;
 416         byteshift = shift / 8;
 417         bitshift = shift & 7;
 418
 419         if (a->length_W == 0) {
 420             return;
 421         }
 422
 423         if(bigint_get_first_set_bit(a) < shift){ /* we would shift out more than we have */
 424                 bigint_set_zero(a);
 425                 return;
 426         }
 427
 428         if(byteshift){
 429                 memmove(a->wordv, (uint8_t*)a->wordv + byteshift, a->length_W * sizeof(bigint_word_t) - byteshift);
 430                 memset((uint8_t*)&a->wordv[a->length_W] - byteshift, 0, byteshift);
 431             a->length_W -= byteshift / sizeof(bigint_word_t);
 432         }
 433
 434
 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 __attribute__((packed)) {
 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 = (bigint_wordplus_t)a->wordv[i] * (bigint_wordplus_t)a->wordv[i] + (bigint_wordplus_t)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 = (bigint_wordplus_t)a->wordv[i] * (bigint_wordplus_t)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 = (bigint_wordplus_t)acc.u[1] + c2;
 737             c2 = q + (c1 < c2);
 738         }
 739         dest->wordv[i + a->length_W] += c1;
 740         if (i < a->length_W - 1) {
 741             dest->wordv[i + a->length_W + 1] = c2 + (dest->wordv[i + a->length_W] < c1);
 742         }
 743     }
 744     dest->info = 0;
 745     dest->length_W = 2 * a->length_W;
 746     bigint_adjust(dest);
 747 }
 748
 749 /******************************************************************************/
 750
 751 void bigint_sub_u_bitscale(bigint_t *a, const bigint_t *b, bigint_length_t bitscale){
 752         bigint_t tmp, x;
 753         bigint_word_t tmp_b[b->length_W + 1];
 754         const bigint_length_t word_shift = bitscale / BIGINT_WORD_SIZE;
 755
 756         if (a->length_W < b->length_W + word_shift) {
 757 #if DEBUG
 758                 cli_putstr("\r\nDBG: *bang*\r\n");
 759 #endif
 760                 bigint_set_zero(a);
 761                 return;
 762         }
 763         tmp.wordv = tmp_b;
 764         bigint_copy(&tmp, b);
 765         bigint_shiftleft(&tmp, bitscale % BIGINT_WORD_SIZE);
 766
 767         x.info = a->info;
 768         x.wordv = &(a->wordv[word_shift]);
 769         x.length_W = a->length_W - word_shift;
 770
 771         bigint_sub_u(&x, &x, &tmp);
 772         bigint_adjust(a);
 773         return;
 774 }
 775
 776 /******************************************************************************/
 777
 778 void bigint_reduce(bigint_t *a, const bigint_t *r){
 779         uint8_t rfbs = GET_FBS(r);
 780         if (r->length_W == 0 || a->length_W == 0) {
 781                 return;
 782         }
 783
 784         if (bigint_length_b(a) + 3 > bigint_length_b(r)) {
 785         if ((r->length_W * sizeof(bigint_word_t) <= 4) && (a->length_W * sizeof(bigint_word_t) <= 4)) {
 786             uint32_t p = 0, q = 0;
 787             memcpy(&p, a->wordv, a->length_W * sizeof(bigint_word_t));
 788             memcpy(&q, r->wordv, r->length_W * sizeof(bigint_word_t));
 789             p %= q;
 790             memcpy(a->wordv, &p, a->length_W * sizeof(bigint_word_t));
 791             a->length_W = r->length_W;
 792             bigint_adjust(a);
 793             return;
 794         }
 795         unsigned shift;
 796         while (a->length_W > r->length_W) {
 797             shift = (a->length_W - r->length_W) * CHAR_BIT * sizeof(bigint_word_t) + GET_FBS(a) - rfbs - 1;
 798             bigint_sub_u_bitscale(a, r, shift);
 799         }
 800         while ((GET_FBS(a) > rfbs) && (a->length_W == r->length_W)) {
 801             shift = GET_FBS(a) - rfbs - 1;
 802             bigint_sub_u_bitscale(a, r, shift);
 803         }
 804         }
 805         while (bigint_cmp_u(a, r) >= 0) {
 806                 bigint_sub_u(a, a, r);
 807         }
 808         bigint_adjust(a);
 809 }
 810
 811 /******************************************************************************/
 812
 813 /* calculate dest = a**exp % r */
 814 /* using square&multiply */
 815 void bigint_expmod_u_sam(bigint_t *dest, const bigint_t *a, const bigint_t *exp, const bigint_t *r){
 816         if(a->length_W == 0){
 817             bigint_set_zero(dest);
 818                 return;
 819         }
 820
 821         bigint_t res, base;
 822         bigint_word_t t, base_w[MAX(a->length_W, r->length_W)], res_w[r->length_W * 2];
 823         bigint_length_t i;
 824         uint8_t j;
 825         res.wordv = res_w;
 826         base.wordv = base_w;
 827         bigint_copy(&base, a);
 828         bigint_reduce(&base, r);
 829         res.wordv[0] = 1;
 830         res.length_W = 1;
 831         res.info = 0;
 832         bigint_adjust(&res);
 833         if(exp->length_W == 0){
 834                 bigint_copy(dest, &res);
 835                 return;
 836         }
 837         bigint_copy(&res, &base);
 838         uint8_t flag = 0;
 839         t = exp->wordv[exp->length_W - 1];
 840         for (i = exp->length_W; i > 0; --i) {
 841                 t = exp->wordv[i - 1];
 842                 for (j = BIGINT_WORD_SIZE; j > 0; --j) {
 843                         if (!flag) {
 844                                 if (t & (((bigint_word_t)1) << (BIGINT_WORD_SIZE - 1))) {
 845                                         flag = 1;
 846                                 }
 847                         } else {
 848                                 bigint_square(&res, &res);
 849                                 bigint_reduce(&res, r);
 850                                 if(t & (((bigint_word_t)1) << (BIGINT_WORD_SIZE - 1))){
 851                                         bigint_mul_u(&res, &res, &base);
 852                                         bigint_reduce(&res, r);
 853                                 }
 854                         }
 855                         t <<= 1;
 856                 }
 857         }
 858
 859         SET_POS(&res);
 860         bigint_copy(dest, &res);
 861 }
 862
 863 /******************************************************************************/
 864 /* gcd <-- gcd(x,y) a*x+b*y=gcd */
 865 void bigint_gcdext(bigint_t *gcd, bigint_t *a, bigint_t *b, const bigint_t *x, const bigint_t *y){
 866          bigint_length_t i = 0;
 867          if(x->length_W == 0 || y->length_W == 0){
 868              if(gcd){
 869                  bigint_set_zero(gcd);
 870              }
 871              if(a){
 872                  bigint_set_zero(a);
 873              }
 874          if(b){
 875              bigint_set_zero(b);
 876          }
 877                  return;
 878          }
 879          if(x->length_W == 1 && x->wordv[0] == 1){
 880              if(gcd){
 881              gcd->length_W = 1;
 882              gcd->wordv[0] = 1;
 883              gcd->info = 0;
 884              }
 885                  if(a){
 886                          a->length_W = 1;
 887                          a->wordv[0] = 1;
 888                          SET_POS(a);
 889                          bigint_adjust(a);
 890                  }
 891                  if(b){
 892                          bigint_set_zero(b);
 893                  }
 894                  return;
 895          }
 896          if(y->length_W == 1 && y->wordv[0] == 1){
 897                  if(gcd){
 898              gcd->length_W = 1;
 899              gcd->wordv[0] = 1;
 900              gcd->info = 0;
 901                  }
 902                  if(b){
 903                          b->length_W = 1;
 904                          b->wordv[0] = 1;
 905                          SET_POS(b);
 906                          bigint_adjust(b);
 907                  }
 908                  if(a){
 909                          bigint_set_zero(a);
 910                  }
 911                  return;
 912          }
 913
 914          while(x->wordv[i] == 0 && y->wordv[i] == 0){
 915                  ++i;
 916          }
 917          bigint_word_t g_b[i + 2], x_b[x->length_W - i], y_b[y->length_W - i];
 918          bigint_word_t u_b[x->length_W - i], v_b[y->length_W - i];
 919          bigint_word_t a_b[y->length_W + 2], c_b[y->length_W + 2];
 920          bigint_word_t b_b[x->length_W + 2], d_b[x->length_W + 2];
 921      bigint_t g, x_, y_, u, v, a_, b_, c_, d_;
 922
 923          g.wordv = g_b;
 924          x_.wordv = x_b;
 925          y_.wordv = y_b;
 926          memset(g_b, 0, i * sizeof(bigint_word_t));
 927          g_b[i] = 1;
 928          g.length_W = i + 1;
 929          g.info = 0;
 930          x_.info = y_.info = 0;
 931          x_.length_W = x->length_W - i;
 932          y_.length_W = y->length_W - i;
 933          memcpy(x_.wordv, x->wordv + i, x_.length_W * sizeof(bigint_word_t));
 934          memcpy(y_.wordv, y->wordv + i, y_.length_W * sizeof(bigint_word_t));
 935          for(i = 0; (x_.wordv[0] & (1 << i)) == 0 && (y_.wordv[0] & (1 << i)) == 0; ++i){
 936          }
 937
 938          bigint_adjust(&x_);
 939          bigint_adjust(&y_);
 940
 941          if(i){
 942                  bigint_shiftleft(&g, i);
 943                  bigint_shiftright(&x_, i);
 944                  bigint_shiftright(&y_, i);
 945          }
 946
 947          u.wordv = u_b;
 948          v.wordv = v_b;
 949          a_.wordv = a_b;
 950          b_.wordv = b_b;
 951          c_.wordv = c_b;
 952          d_.wordv = d_b;
 953
 954          bigint_copy(&u, &x_);
 955          bigint_copy(&v, &y_);
 956          a_.wordv[0] = 1;
 957          a_.length_W = 1;
 958          a_.info = 0;
 959          d_.wordv[0] = 1;
 960          d_.length_W = 1;
 961          d_.info = 0;
 962          bigint_set_zero(&b_);
 963          bigint_set_zero(&c_);
 964          do {
 965                  while ((u.wordv[0] & 1) == 0) {
 966                          bigint_shiftright(&u, 1);
 967                          if((a_.wordv[0] & 1) || (b_.wordv[0] & 1)){
 968                                  bigint_add_s(&a_, &a_, &y_);
 969                                  bigint_sub_s(&b_, &b_, &x_);
 970                          }
 971                          bigint_shiftright(&a_, 1);
 972                          bigint_shiftright(&b_, 1);
 973                  }
 974                  while ((v.wordv[0] & 1) == 0) {
 975                      bigint_shiftright(&v, 1);
 976                          if((c_.wordv[0] & 1) || (d_.wordv[0] & 1)){
 977                                  bigint_add_s(&c_, &c_, &y_);
 978                                  bigint_sub_s(&d_, &d_, &x_);
 979                          }
 980                          bigint_shiftright(&c_, 1);
 981                          bigint_shiftright(&d_, 1);
 982                  }
 983                  if(bigint_cmp_u(&u, &v) >= 0){
 984                         bigint_sub_u(&u, &u, &v);
 985                         bigint_sub_s(&a_, &a_, &c_);
 986                         bigint_sub_s(&b_, &b_, &d_);
 987                  }else{
 988                         bigint_sub_u(&v, &v, &u);
 989                         bigint_sub_s(&c_, &c_, &a_);
 990                         bigint_sub_s(&d_, &d_, &b_);
 991                  }
 992          } while(u.length_W);
 993          if(gcd){
 994                  bigint_mul_s(gcd, &v, &g);
 995          }
 996          if(a){
 997                 bigint_copy(a, &c_);
 998          }
 999          if(b){
1000                  bigint_copy(b, &d_);
1001          }
1002 }
1003
1004 /******************************************************************************/
1005
1006 void bigint_inverse(bigint_t *dest, const bigint_t *a, const bigint_t *m){
1007         bigint_gcdext(NULL, dest, NULL, a, m);
1008         while(dest->info&BIGINT_NEG_MASK){
1009                 bigint_add_s(dest, dest, m);
1010         }
1011 }
1012
1013 /******************************************************************************/
1014
1015 void bigint_changeendianess(bigint_t *a){
1016         uint8_t t, *p, *q;
1017         p = (uint8_t*)(a->wordv);
1018         q = p + a->length_W * sizeof(bigint_word_t) - 1;
1019         while(p<q){
1020                 t = *p;
1021                 *p = *q;
1022                 *q = t;
1023                 ++p; --q;
1024         }
1025 }
1026
1027 /******************************************************************************/
1028
1029 void bigint_mul_word_u(bigint_t *a, bigint_word_t b){
1030     bigint_wordplus_t c0 = 0, c1 = 0;
1031     bigint_length_t i;
1032
1033     if(b == 0){
1034         bigint_set_zero(a);
1035         return;
1036     }
1037
1038     for(i = 0; i < a->length_W; ++i){
1039         c1 = ((bigint_wordplus_t)(a->wordv[i])) * ((bigint_wordplus_t)b);
1040         c1 += c0;
1041         a->wordv[i] = (bigint_word_t)c1;
1042         c0 = c1 >> BIGINT_WORD_SIZE;
1043     }
1044     if(c0){
1045         a->wordv[a->length_W] = (bigint_word_t)c0;
1046         a->length_W += 1;
1047     }
1048     bigint_adjust(a);
1049 }
1050
1051 /******************************************************************************/
1052 #if 1
1053
1054 void bigint_clip(bigint_t *dest, bigint_length_t length_W){
1055     if(dest->length_W > length_W){
1056         dest->length_W = length_W;
1057     }
1058     bigint_adjust(dest);
1059 }
1060 /******************************************************************************/
1061
1062 /*
1063  * m_ = m * m'[0]
1064  * dest = (a * b) % m (?)
1065  */
1066
1067 void bigint_mont_mul(bigint_t *dest, const bigint_t *a, const bigint_t *b, const bigint_t *m, const bigint_t *m_){
1068     const bigint_length_t s = MAX(MAX(a->length_W, b->length_W), m->length_W);
1069     bigint_t u, t;
1070     bigint_word_t u_w[s + 2], t_w[s + 2];
1071     bigint_length_t i;
1072
1073     if (a->length_W == 0 || b->length_W == 0) {
1074         bigint_set_zero(dest);
1075         return;
1076     }
1077     u.wordv = u_w;
1078     u.info = 0;
1079     u.length_W = 0;
1080     t.wordv = t_w;
1081     for (i = 0; i < a->length_W; ++i) {
1082         bigint_copy(&t, b);
1083         bigint_mul_word_u(&t, a->wordv[i]);
1084         bigint_add_u(&u, &u, &t);
1085         bigint_copy(&t, m_);
1086         if (u.length_W != 0) {
1087             bigint_mul_word_u(&t, u.wordv[0]);
1088             bigint_add_u(&u, &u, &t);
1089         }
1090         bigint_shiftright(&u, BIGINT_WORD_SIZE);
1091     }
1092     for (; i < s; ++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(&u, BIGINT_WORD_SIZE);
1099     }
1100     bigint_reduce(&u, m);
1101     bigint_copy(dest, &u);
1102 }
1103
1104 /******************************************************************************/
1105
1106 void bigint_mont_red(bigint_t *dest, const bigint_t *a, const bigint_t *m, const bigint_t *m_){
1107     bigint_t u, t;
1108     bigint_length_t i, s = MAX(a->length_W, m->length_W);
1109     bigint_word_t u_w[s + 2], t_w[s + 2];
1110
1111     t.wordv = t_w;
1112     u.wordv = u_w;
1113     if (a->length_W == 0) {
1114         bigint_set_zero(dest);
1115         return;
1116     }
1117     bigint_copy(&u, a);
1118     for (i = 0; i < m->length_W; ++i) {
1119         bigint_copy(&t, m_);
1120         if (u.length_W != 0) {
1121             bigint_mul_word_u(&t, u.wordv[0]);
1122             bigint_add_u(&u, &u, &t);
1123         }
1124         bigint_shiftright(&u, BIGINT_WORD_SIZE);
1125     }
1126     bigint_reduce(&u, m);
1127     bigint_copy(dest, &u);
1128 }
1129
1130 /******************************************************************************/
1131 /*
1132  * m_ = m * (- m0^-1 (mod 2^W))
1133  */
1134 void bigint_mont_gen_m_(bigint_t* dest, const bigint_t* m){
1135     bigint_word_t x_w[2], m_w_0[1];
1136     bigint_t x, m_0;
1137 #if 0
1138     printf("\nm = ");
1139     bigint_print_hex(m);
1140     uart0_flush();
1141 #endif
1142     if (m->length_W == 0) {
1143         bigint_set_zero(dest);
1144         return;
1145     }
1146     if ((m->wordv[0] & 1) == 0) {
1147         printf_P(PSTR("ERROR: m must not be even, m = "));
1148         bigint_print_hex(m);
1149         putchar('\n');
1150         uart0_flush();
1151         return;
1152     }
1153     x.wordv = x_w;
1154     x.info = 0;
1155     x.length_W = 2;
1156     x_w[0] = 0;
1157     x_w[1] = 1;
1158     m_0.wordv = m_w_0;
1159     m_0.info = 0;
1160     m_0.length_W = 1;
1161     m_0.wordv[0] = m->wordv[0];
1162     bigint_adjust(&x);
1163     bigint_adjust(&m_0);
1164 #if 0
1165     printf("\nm0 = ");
1166     bigint_print_hex(&m_0);
1167     printf("\nx = ");
1168     bigint_print_hex(&x);
1169     uart0_flush();
1170 #endif
1171     bigint_inverse(dest, &m_0, &x);
1172 #if 0
1173     printf("\nm0^-1 = ");
1174     bigint_print_hex(dest);
1175     uart0_flush();
1176 #endif
1177     bigint_sub_s(&x, &x, dest);
1178 #if 0
1179     printf("\n-m0^-1 = ");
1180     bigint_print_hex(&x);
1181     uart0_flush();
1182 #endif
1183     bigint_copy(dest, m);
1184     bigint_mul_word_u(dest, x.wordv[0]);
1185 }
1186
1187 /******************************************************************************/
1188
1189 /*
1190  * dest = a * R mod m
1191  */
1192 void bigint_mont_trans(bigint_t *dest, const bigint_t *a, const bigint_t *m){
1193     bigint_t t;
1194     bigint_word_t t_w[a->length_W + m->length_W];
1195
1196     t.wordv = t_w;
1197     memset(t_w, 0, m->length_W * sizeof(bigint_word_t));
1198     memcpy(&t_w[m->length_W], a->wordv, a->length_W * sizeof(bigint_word_t));
1199     t.info = a->info;
1200     t.length_W = a->length_W + m->length_W;
1201     bigint_reduce(&t, m);
1202     bigint_copy(dest, &t);
1203 }
1204
1205 /******************************************************************************/
1206
1207 /* calculate dest = a**exp % r */
1208 /* using square&multiply */
1209 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_){
1210     if(r->length_W == 0) {
1211         return;
1212     }
1213
1214     bigint_length_t s = r->length_W;
1215     bigint_t res, ax;
1216     bigint_word_t t, res_w[r->length_W * 2], ax_w[MAX(s, a->length_W)];
1217     bigint_length_t i;
1218     uint8_t j;
1219
1220     res.wordv = res_w;
1221     ax.wordv = ax_w;
1222
1223     res.wordv[0] = 1;
1224     res.length_W = 1;
1225     res.info = 0;
1226     bigint_adjust(&res);
1227     if (exp->length_W == 0) {
1228         bigint_copy(dest, &res);
1229         return;
1230     }
1231     bigint_copy(&ax, a);
1232     bigint_reduce(&ax, r);
1233     bigint_mont_trans(&ax, &ax, r);
1234     bigint_mont_trans(&res, &res, r);
1235     uint8_t flag = 0;
1236     t = exp->wordv[exp->length_W - 1];
1237     for (i = exp->length_W; i > 0; --i) {
1238         t = exp->wordv[i - 1];
1239         for(j = BIGINT_WORD_SIZE; j > 0; --j){
1240             if (!flag) {
1241                 if(t & (((bigint_wordplus_t)1) << (BIGINT_WORD_SIZE - 1))){
1242                     flag = 1;
1243                 }
1244             }
1245             if (flag) {
1246                 bigint_square(&res, &res);
1247                 bigint_mont_red(&res, &res, r, m_);
1248                 if (t & (((bigint_wordplus_t)1) << (BIGINT_WORD_SIZE - 1))) {
1249                     bigint_mont_mul(&res, &res, &ax, r, m_);
1250                 }
1251             }
1252             t <<= 1;
1253         }
1254     }
1255     SET_POS(&res);
1256     bigint_mont_red(dest, &res, r, m_);
1257 }
1258
1259 /******************************************************************************/
1260
1261 void bigint_expmod_u_mont_sam(bigint_t *dest, const bigint_t *a, const bigint_t *exp, const bigint_t *r){
1262     if(r->length_W == 0) {
1263         return;
1264     }
1265     if(a->length_W == 0) {
1266         bigint_set_zero(dest);
1267         return;
1268     }
1269     bigint_t m_;
1270     bigint_word_t m_w_[r->length_W + 1];
1271     m_.wordv = m_w_;
1272     bigint_mont_gen_m_(&m_, r);
1273     bigint_expmod_u_mont_accel(dest, a, exp, r,&m_);
1274 }
1275
1276 /******************************************************************************/
1277
1278 #endif
1279
1280 void bigint_expmod_u(bigint_t *dest, const bigint_t *a, const bigint_t *exp, const bigint_t *r){
1281     if (r->wordv[0] & 1) {
1282         bigint_expmod_u_mont_sam(dest, a, exp, r);
1283     } else {
1284         bigint_expmod_u_sam(dest, a, exp, r);
1285     }
1286 }
1287
1288
1289
1290
1291
1292
1293
1294
1295
1296
1297
1298
1299
1300
1301
1302