bigint/bigint.c

   1 /* bigint.c */
   2 /*
   3     This file is part of the ARM-Crypto-Lib.
   4     Copyright (C) 2008  Daniel Otte (daniel.otte@rub.de)
   5
   6     This program is free software: you can redistribute it and/or modify
   7     it under the terms of the GNU General Public License as published by
   8     the Free Software Foundation, either version 3 of the License, or
   9     (at your option) any later version.
  10
  11     This program is distributed in the hope that it will be useful,
  12     but WITHOUT ANY WARRANTY; without even the implied warranty of
  13     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14     GNU General Public License for more details.
  15
  16     You should have received a copy of the GNU General Public License
  17     along with this program.  If not, see <http://www.gnu.org/licenses/>.
  18 */
  19 /**
  20  * \file                bigint.c
  21  * \author              Daniel Otte
  22  * \date                2010-02-22
  23  *
  24  * \license         GPLv3 or later
  25  *
  26  */
  27
  28
  29 #define STRING2(x) #x
  30 #define STRING(x) STRING2(x)
  31 #define STR_LINE STRING(__LINE__)
  32
  33 #include "bigint.h"
  34 #include <string.h>
  35
  36 #define DEBUG 0
  37
  38 #if DEBUG || 1
  39 #include "cli.h"
  40 #include "bigint_io.h"
  41 #include <stdio.h>
  42 #endif
  43
  44 #ifndef MAX
  45  #define MAX(a,b) (((a)>(b))?(a):(b))
  46 #endif
  47
  48 #ifndef MIN
  49  #define MIN(a,b) (((a)<(b))?(a):(b))
  50 #endif
  51
  52 #define SET_FBS(a, v) do{(a)->info &=~BIGINT_FBS_MASK; (a)->info |= (v);}while(0)
  53 #define GET_FBS(a)   ((a)->info&BIGINT_FBS_MASK)
  54 #define SET_NEG(a)   (a)->info |= BIGINT_NEG_MASK
  55 #define SET_POS(a)   (a)->info &= ~BIGINT_NEG_MASK
  56 #define XCHG(a,b)    do{(a)^=(b); (b)^=(a); (a)^=(b);}while(0)
  57 #define XCHG_PTR(a,b)    do{ a = (void*)(((intptr_t)(a)) ^ ((intptr_t)(b))); \
  58                                  b = (void*)(((intptr_t)(a)) ^ ((intptr_t)(b))); \
  59                                  a = (void*)(((intptr_t)(a)) ^ ((intptr_t)(b)));}while(0)
  60
  61 #define GET_SIGN(a) ((a)->info&BIGINT_NEG_MASK)
  62
  63 /******************************************************************************/
  64 void bigint_adjust(bigint_t *a){
  65         while(a->length_W!=0 && a->wordv[a->length_W-1]==0){
  66                 a->length_W--;
  67         }
  68         if(a->length_W==0){
  69                 a->info=0;
  70                 return;
  71         }
  72         bigint_word_t t;
  73         uint8_t i = BIGINT_WORD_SIZE-1;
  74         t = a->wordv[a->length_W-1];
  75         while((t&(1L<<(BIGINT_WORD_SIZE-1)))==0 && i){
  76                 t<<=1;
  77                 i--;
  78         }
  79         SET_FBS(a, i);
  80 }
  81
  82 /******************************************************************************/
  83
  84 uint16_t bigint_length_b(const bigint_t *a){
  85         if(!a->length_W || a->length_W==0){
  86                 return 0;
  87         }
  88         return (a->length_W-1) * BIGINT_WORD_SIZE + GET_FBS(a);
  89 }
  90
  91 /******************************************************************************/
  92
  93 uint16_t bigint_length_B(const bigint_t *a){
  94         return a->length_W * sizeof(bigint_word_t);
  95 }
  96
  97 /******************************************************************************/
  98
  99 uint32_t bigint_get_first_set_bit(const bigint_t *a){
 100         if(a->length_W==0){
 101                 return (uint32_t)(-1);
 102         }
 103         return (a->length_W-1)*sizeof(bigint_word_t)*8+GET_FBS(a);
 104 }
 105
 106
 107 /******************************************************************************/
 108
 109 uint32_t bigint_get_last_set_bit(const bigint_t *a){
 110         uint32_t r=0;
 111         uint8_t b=0;
 112         bigint_word_t x=1;
 113         if(a->length_W==0){
 114                 return (uint32_t)(-1);
 115         }
 116         while(a->wordv[r]==0 && r<a->length_W){
 117                 ++r;
 118         }
 119         if(a->wordv[r] == 0){
 120                 return (uint32_t)(-1);
 121         }
 122         while((x&a->wordv[r])==0){
 123                 ++b;
 124                 x <<= 1;
 125         }
 126         return r*BIGINT_WORD_SIZE+b;
 127 }
 128
 129 /******************************************************************************/
 130
 131 void bigint_copy(bigint_t *dest, const bigint_t *src){
 132     if(dest->wordv != src->wordv){
 133             memcpy(dest->wordv, src->wordv, src->length_W * sizeof(bigint_word_t));
 134     }
 135     dest->length_W = src->length_W;
 136         dest->info = src->info;
 137 }
 138
 139 /******************************************************************************/
 140
 141 /* this should be implemented in assembly */
 142 void bigint_add_u(bigint_t *dest, const bigint_t *a, const bigint_t *b){
 143         uint16_t i;
 144         bigint_wordplus_t t = 0LL;
 145         if(a->length_W < b->length_W){
 146                 XCHG_PTR(a,b);
 147         }
 148         for(i = 0; i < b->length_W; ++i){
 149 //              t = (bigint_wordplus_t)(a->wordv[i]) + (bigint_wordplus_t)(b->wordv[i]) + t;
 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         bigint_adjust(dest);
 165 }
 166
 167 /******************************************************************************/
 168
 169 /* this should be implemented in assembly */
 170 void bigint_add_scale_u(bigint_t *dest, const bigint_t *a, uint16_t scale){
 171         if(a->length_W == 0){
 172                 return;
 173         }
 174         if(scale == 0){
 175                 bigint_add_u(dest, dest, a);
 176                 return;
 177         }
 178         bigint_t x;
 179 #if BIGINT_WORD_SIZE == 8
 180         memset(dest->wordv + dest->length_W, 0, MAX(dest->length_W, a->length_W + scale) - dest->length_W);
 181         x.wordv = dest->wordv + scale;
 182         x.length_W = dest->length_W - scale;
 183         if((int16_t)x.length_W < 0){
 184                 x.length_W = 0;
 185                 x.info = 0;
 186         } else {
 187                 x.info = dest->info;
 188         }
 189         bigint_add_u(&x, &x, a);
 190         dest->length_W = x.length_W + scale;
 191         dest->info = 0;
 192         bigint_adjust(dest);
 193 #else
 194         bigint_t s;
 195         uint16_t word_shift = scale / sizeof(bigint_word_t), byte_shift = scale % sizeof(bigint_word_t);
 196         bigint_word_t bv[a->length_W + 1];
 197         s.wordv = bv;
 198         bv[0] = bv[a->length_W] = 0;
 199         memcpy((uint8_t*)bv + byte_shift, a->wordv, a->length_W * sizeof(bigint_word_t));
 200         s.length_W = a->length_W + 1;
 201         bigint_adjust(&s);
 202         memset(dest->wordv + dest->length_W, 0, (MAX(dest->length_W, s.length_W + word_shift) - dest->length_W) * sizeof(bigint_word_t));
 203         x.wordv = dest->wordv + word_shift;
 204         x.length_W = dest->length_W - word_shift;
 205         if((int16_t)x.length_W < 0){
 206                 x.length_W = 0;
 207                 x.info = 0;
 208         }else{
 209                 x.info = dest->info;
 210         }
 211         bigint_add_u(&x, &x, &s);
 212         dest->length_W = x.length_W + word_shift;
 213         dest->info = 0;
 214         bigint_adjust(dest);
 215 #endif
 216 }
 217
 218 /******************************************************************************/
 219
 220 /* this should be implemented in assembly */
 221 void bigint_sub_u(bigint_t *dest, const bigint_t *a, const bigint_t *b){
 222         int8_t borrow=0;
 223         int8_t  r;
 224         bigint_wordplus_signed_t t=0LL;
 225         uint16_t i;
 226         if(b->length_W == 0){
 227                 bigint_copy(dest, a);
 228                 SET_POS(dest);
 229                 return;
 230         }
 231         if(a->length_W == 0){
 232                 bigint_copy(dest, b);
 233                 SET_NEG(dest);
 234                 return;
 235         }
 236     r = bigint_cmp_u(a,b);
 237     if(r == 0){
 238         bigint_set_zero(dest);
 239         return;
 240     }
 241         if(r < 0){
 242                 bigint_sub_u(dest, b, a);
 243                 SET_NEG(dest);
 244                 return;
 245         }
 246         for(i = 0; i < a->length_W; ++i){
 247                 t = a->wordv[i];
 248                 if(i < b->length_W){
 249                         t -= b->wordv[i];
 250                 }
 251                 t -= borrow;
 252                 dest->wordv[i]=(bigint_word_t)t;
 253                 if(t<0){
 254                         borrow = 1;
 255                 }else{
 256                         borrow = 0;
 257                 }
 258         }
 259         SET_POS(dest);
 260         dest->length_W = i;
 261         bigint_adjust(dest);
 262 }
 263
 264 /******************************************************************************/
 265
 266 int8_t bigint_cmp_u(const bigint_t *a, const bigint_t *b){
 267         if(a->length_W > b->length_W){
 268                 return 1;
 269         }
 270         if(a->length_W < b->length_W){
 271                 return -1;
 272         }
 273         if(a->length_W==0){
 274                 return 0;
 275         }
 276         uint16_t i;
 277         i = a->length_W-1;
 278         do{
 279                 if(a->wordv[i] != b->wordv[i]){
 280                         if(a->wordv[i] > b->wordv[i]){
 281                                 return 1;
 282                         }else{
 283                                 return -1;
 284                         }
 285                 }
 286         }while(i--);
 287         return 0;
 288 }
 289
 290 /******************************************************************************/
 291
 292 void bigint_add_s(bigint_t *dest, const bigint_t *a, const bigint_t *b){
 293         uint8_t s;
 294         s  = GET_SIGN(a)?2:0;
 295         s |= GET_SIGN(b)?1:0;
 296         switch(s){
 297                 case 0: /* both positive */
 298                         bigint_add_u(dest, a,b);
 299                         SET_POS(dest);
 300                         break;
 301                 case 1: /* a positive, b negative */
 302                         bigint_sub_u(dest, a, b);
 303                         break;
 304                 case 2: /* a negative, b positive */
 305                         bigint_sub_u(dest, b, a);
 306                         break;
 307                 case 3: /* both negative */
 308                         bigint_add_u(dest, a, b);
 309                         SET_NEG(dest);
 310                         break;
 311                 default: /* how can this happen?*/
 312                         break;
 313         }
 314 }
 315
 316 /******************************************************************************/
 317
 318 void bigint_sub_s(bigint_t *dest, const bigint_t *a, const bigint_t *b){
 319         uint8_t s;
 320         s  = GET_SIGN(a)?2:0;
 321         s |= GET_SIGN(b)?1:0;
 322         switch(s){
 323                 case 0: /* both positive */
 324                         bigint_sub_u(dest, a,b);
 325                         break;
 326                 case 1: /* a positive, b negative */
 327                         bigint_add_u(dest, a, b);
 328                         SET_POS(dest);
 329                         break;
 330                 case 2: /* a negative, b positive */
 331                         bigint_add_u(dest, a, b);
 332                         SET_NEG(dest);
 333                         break;
 334                 case 3: /* both negative */
 335                         bigint_sub_u(dest, b, a);
 336                         break;
 337                 default: /* how can this happen?*/
 338                                         break;
 339         }
 340
 341 }
 342
 343 /******************************************************************************/
 344
 345 int8_t bigint_cmp_s(const bigint_t *a, const bigint_t *b){
 346         uint8_t s;
 347         if(a->length_W==0 && b->length_W==0){
 348                 return 0;
 349         }
 350         s  = GET_SIGN(a)?2:0;
 351         s |= GET_SIGN(b)?1:0;
 352         switch(s){
 353                 case 0: /* both positive */
 354                         return bigint_cmp_u(a, b);
 355                         break;
 356                 case 1: /* a positive, b negative */
 357                         return 1;
 358                         break;
 359                 case 2: /* a negative, b positive */
 360                         return -1;
 361                         break;
 362                 case 3: /* both negative */
 363                         return bigint_cmp_u(b, a);
 364                         break;
 365                 default: /* how can this happen?*/
 366                                         break;
 367         }
 368         return 0; /* just to satisfy the compiler */
 369 }
 370
 371 /******************************************************************************/
 372
 373 void bigint_shiftleft(bigint_t *a, uint16_t shift){
 374         uint16_t byteshift, words_to_shift;
 375         int16_t i;
 376         uint8_t bitshift;
 377         bigint_word_t *p;
 378         bigint_wordplus_t t = 0;
 379         if(shift == 0){
 380                 return;
 381         }
 382         byteshift = shift / 8;
 383         bitshift = shift & 7;
 384
 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, uint16_t shift){
 412         uint16_t byteshift;
 413         uint16_t i;
 414         uint8_t bitshift;
 415         bigint_wordplus_t t = 0;
 416         byteshift = shift / 8;
 417         bitshift = shift & 7;
 418
 419         if(byteshift >= a->length_W * sizeof(bigint_word_t)){ /* we would shift out more than we have */
 420                 bigint_set_zero(a);
 421                 return;
 422         }
 423         if(byteshift == a->length_W * sizeof(bigint_word_t) - 1 && bitshift > GET_FBS(a)){
 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         }
 431
 432     byteshift /= sizeof(bigint_word_t); /* byteshift is now wordshift */
 433     a->length_W -= byteshift;
 434         if(bitshift != 0 && a->length_W){
 435          /* shift to the right */
 436                 i = a->length_W - 1;
 437                 do{
 438                         t |= ((bigint_wordplus_t)(a->wordv[i])) << (BIGINT_WORD_SIZE - bitshift);
 439                         a->wordv[i] = (bigint_word_t)(t >> BIGINT_WORD_SIZE);
 440                         t <<= BIGINT_WORD_SIZE;
 441                 }while(i--);
 442         }
 443         bigint_adjust(a);
 444 }
 445
 446 /******************************************************************************/
 447
 448 void bigint_xor(bigint_t *dest, const bigint_t *a){
 449         uint16_t i;
 450         for(i=0; i<a->length_W; ++i){
 451                 dest->wordv[i] ^= a->wordv[i];
 452         }
 453         bigint_adjust(dest);
 454 }
 455
 456 /******************************************************************************/
 457
 458 void bigint_set_zero(bigint_t *a){
 459         a->length_W=0;
 460 }
 461
 462 /******************************************************************************/
 463
 464 /* using the Karatsuba-Algorithm */
 465 /* x*y = (xh*yh)*b**2n + ((xh+xl)*(yh+yl) - xh*yh - xl*yl)*b**n + yh*yl */
 466 void bigint_mul_u(bigint_t *dest, const bigint_t *a, const bigint_t *b){
 467         if(a->length_W == 0 || b->length_W == 0){
 468                 bigint_set_zero(dest);
 469                 return;
 470         }
 471         if(dest == a || dest == b){
 472                 bigint_t d;
 473                 bigint_word_t d_b[a->length_W + b->length_W];
 474                 d.wordv = d_b;
 475                 bigint_mul_u(&d, a, b);
 476                 bigint_copy(dest, &d);
 477                 return;
 478         }
 479         if(a->length_W == 1 || b->length_W == 1){
 480                 if(a->length_W != 1){
 481                         XCHG_PTR(a,b);
 482                 }
 483                 bigint_wordplus_t t = 0;
 484                 uint16_t i;
 485                 bigint_word_t x = a->wordv[0];
 486                 for(i=0; i < b->length_W; ++i){
 487                         t += ((bigint_wordplus_t)b->wordv[i]) * ((bigint_wordplus_t)x);
 488                         dest->wordv[i] = (bigint_word_t)t;
 489                         t >>= BIGINT_WORD_SIZE;
 490                 }
 491                 dest->length_W = i;
 492                 if(t){
 493                     dest->wordv[i] = (bigint_word_t)t;
 494                     dest->length_W += 1;
 495                 }
 496                 dest->info = 0;
 497                 bigint_adjust(dest);
 498                 return;
 499         }
 500         if(a->length_W * sizeof(bigint_word_t) <= 4 && b->length_W * sizeof(bigint_word_t) <= 4){
 501                 uint32_t p=0, q=0;
 502                 uint64_t r;
 503                 memcpy(&p, a->wordv, a->length_W*sizeof(bigint_word_t));
 504                 memcpy(&q, b->wordv, b->length_W*sizeof(bigint_word_t));
 505                 r = (uint64_t)p * (uint64_t)q;
 506                 memcpy(dest->wordv, &r, (dest->length_W = a->length_W + b->length_W)*sizeof(bigint_word_t));
 507                 bigint_adjust(dest);
 508                 return;
 509         }
 510         /* split a in xh & xl; split b in yh & yl */
 511         const uint16_t n = (MAX(a->length_W, b->length_W)+1)/2;
 512         bigint_t xl, xh, yl, yh;
 513         xl.wordv = a->wordv;
 514         yl.wordv = b->wordv;
 515         if(a->length_W<=n){
 516                 bigint_set_zero(&xh);
 517                 xl.length_W = a->length_W;
 518                 xl.info = a->info;
 519         }else{
 520                 xl.length_W=n;
 521                 xl.info = 0;
 522                 bigint_adjust(&xl);
 523                 xh.wordv = &(a->wordv[n]);
 524                 xh.length_W = a->length_W-n;
 525                 xh.info = a->info;
 526         }
 527         if(b->length_W<=n){
 528                 bigint_set_zero(&yh);
 529                 yl.length_W = b->length_W;
 530                 yl.info = b->info;
 531         }else{
 532                 yl.length_W=n;
 533                 yl.info = 0;
 534                 bigint_adjust(&yl);
 535                 yh.wordv = &(b->wordv[n]);
 536                 yh.length_W = b->length_W-n;
 537                 yh.info = b->info;
 538         }
 539         /* now we have split up a and b */
 540         /* remember we want to do:
 541          * x*y = (xh * b**n + xl) * (yh * b**n + yl)
 542          *     = (xh * yh) * b**2n + xh * b**n * yl + yh * b**n * xl + xl * yl
 543          *     = (xh * yh) * b**2n + (xh * yl + yh * xl) * b**n + xl *yl
 544          *     // xh * yl + yh * xl = (xh + yh) * (xl + yl) - xh * yh - xl * yl
 545          * x*y = (xh * yh) * b**2n + ((xh+xl)*(yh+yl) - xh*yh - xl*yl)*b**n + xl*yl
 546          *          5              9     2   4   3    7   5   6   1         8   1
 547          */
 548         bigint_word_t  tmp_b[2 * n + 2], m_b[2 * (n + 1)];
 549         bigint_t tmp, tmp2, m;
 550         tmp.wordv = tmp_b;
 551         tmp2.wordv = &(tmp_b[n + 1]);
 552         m.wordv = m_b;
 553
 554         bigint_mul_u(dest, &xl, &yl);  /* 1: dest <= xl*yl     */
 555         bigint_add_u(&tmp2, &xh, &xl); /* 2: tmp2 <= xh+xl     */
 556         bigint_add_u(&tmp, &yh, &yl);  /* 3: tmp  <= yh+yl     */
 557         bigint_mul_u(&m, &tmp2, &tmp); /* 4: m    <= tmp2*tmp  */
 558         bigint_mul_u(&tmp, &xh, &yh);  /* 5: tmp  <= xh*yh     */
 559         bigint_sub_u(&m, &m, dest);    /* 6: m    <= m-dest    */
 560     bigint_sub_u(&m, &m, &tmp);    /* 7: m    <= m-tmp     */
 561         bigint_add_scale_u(dest, &m, n * sizeof(bigint_word_t));       /* 8: dest <= dest+m**n*/
 562         bigint_add_scale_u(dest, &tmp, 2 * n * sizeof(bigint_word_t)); /* 9: dest <= dest+tmp**(2*n) */
 563 }
 564
 565 /******************************************************************************/
 566
 567 void bigint_mul_s(bigint_t *dest, const bigint_t *a, const bigint_t *b){
 568         uint8_t s;
 569         s  = GET_SIGN(a)?2:0;
 570         s |= GET_SIGN(b)?1:0;
 571         switch(s){
 572                 case 0: /* both positive */
 573                         bigint_mul_u(dest, a,b);
 574                         SET_POS(dest);
 575                         break;
 576                 case 1: /* a positive, b negative */
 577                         bigint_mul_u(dest, a,b);
 578                         SET_NEG(dest);
 579                         break;
 580                 case 2: /* a negative, b positive */
 581                         bigint_mul_u(dest, a,b);
 582                         SET_NEG(dest);
 583                         break;
 584                 case 3: /* both negative */
 585                         bigint_mul_u(dest, a,b);
 586                         SET_POS(dest);
 587                         break;
 588                 default: /* how can this happen?*/
 589                         break;
 590         }
 591 }
 592
 593 /******************************************************************************/
 594
 595
 596 #if DEBUG_SQUARE
 597 unsigned square_depth = 0;
 598 #endif
 599
 600 /* square */
 601 /* (xh*b^n+xl)^2 = xh^2*b^2n + 2*xh*xl*b^n + xl^2 */
 602 void bigint_square(bigint_t *dest, const bigint_t *a){
 603         if(a->length_W * sizeof(bigint_word_t) <= 4){
 604                 uint64_t r = 0;
 605                 memcpy(&r, a->wordv, a->length_W * sizeof(bigint_word_t));
 606                 r = r * r;
 607                 memcpy(dest->wordv, &r, 2 * a->length_W * sizeof(bigint_word_t));
 608                 SET_POS(dest);
 609                 dest->length_W = 2 * a->length_W;
 610                 bigint_adjust(dest);
 611                 return;
 612         }
 613         if(dest->wordv == a->wordv){
 614                 bigint_t d;
 615                 bigint_word_t d_b[a->length_W*2];
 616                 d.wordv = d_b;
 617                 bigint_square(&d, a);
 618                 bigint_copy(dest, &d);
 619                 return;
 620         }
 621
 622 #if DEBUG_SQUARE
 623         square_depth += 1;
 624 #endif
 625
 626         uint16_t n;
 627         n=(a->length_W+1)/2;
 628         bigint_t xh, xl, tmp; /* x-high, x-low, temp */
 629         bigint_word_t buffer[2*n+1];
 630         xl.wordv = a->wordv;
 631         xl.length_W = n;
 632         xl.info = 0;
 633         xh.wordv = &(a->wordv[n]);
 634         xh.length_W = a->length_W-n;
 635         xh.info = a->info;
 636         bigint_adjust(&xl);
 637         tmp.wordv = buffer;
 638 /* (xh * b**n + xl)**2 = xh**2 * b**2n + 2 * xh * xl * b**n + xl**2 */
 639 #if DEBUG_SQUARE
 640         if(square_depth == 1){
 641         cli_putstr("\r\nDBG (a): xl: "); bigint_print_hex(&xl);
 642         cli_putstr("\r\nDBG (b): xh: "); bigint_print_hex(&xh);
 643         }
 644 #endif
 645         bigint_square(dest, &xl);
 646 #if DEBUG_SQUARE
 647         if(square_depth == 1){
 648             cli_putstr("\r\nDBG (1): xl**2: "); bigint_print_hex(dest);
 649         }
 650 #endif
 651     bigint_square(&tmp, &xh);
 652 #if DEBUG_SQUARE
 653     if(square_depth == 1){
 654         cli_putstr("\r\nDBG (2): xh**2: "); bigint_print_hex(&tmp);
 655     }
 656 #endif
 657         bigint_add_scale_u(dest, &tmp, 2 * n * sizeof(bigint_word_t));
 658 #if DEBUG_SQUARE
 659         if(square_depth == 1){
 660             cli_putstr("\r\nDBG (3): xl**2 + xh**2*n**2: "); bigint_print_hex(dest);
 661         }
 662 #endif
 663         bigint_mul_u(&tmp, &xl, &xh);
 664 #if DEBUG_SQUARE
 665         if(square_depth == 1){
 666             cli_putstr("\r\nDBG (4): xl*xh: "); bigint_print_hex(&tmp);
 667         }
 668 #endif
 669         bigint_shiftleft(&tmp, 1);
 670 #if DEBUG_SQUARE
 671         if(square_depth == 1){
 672             cli_putstr("\r\nDBG (5): xl*xh*2: "); bigint_print_hex(&tmp);
 673         }
 674 #endif
 675         bigint_add_scale_u(dest, &tmp, n * sizeof(bigint_word_t));
 676 #if DEBUG_SQUARE
 677         if(square_depth == 1){
 678             cli_putstr("\r\nDBG (6): x**2: "); bigint_print_hex(dest);
 679             cli_putstr("\r\n");
 680         }
 681         square_depth -= 1;
 682 #endif
 683 }
 684
 685 /******************************************************************************/
 686
 687 void bigint_sub_u_bitscale(bigint_t *a, const bigint_t *b, uint16_t bitscale){
 688         bigint_t tmp, x;
 689         bigint_word_t tmp_b[b->length_W + 1];
 690         const uint16_t word_shift = bitscale / BIGINT_WORD_SIZE;
 691
 692         if(a->length_W < b->length_W + word_shift){
 693 #if DEBUG
 694                 cli_putstr("\r\nDBG: *bang*\r\n");
 695 #endif
 696                 bigint_set_zero(a);
 697                 return;
 698         }
 699         tmp.wordv = tmp_b;
 700         bigint_copy(&tmp, b);
 701         bigint_shiftleft(&tmp, bitscale % BIGINT_WORD_SIZE);
 702
 703         x.info = a->info;
 704         x.wordv = &(a->wordv[word_shift]);
 705         x.length_W = a->length_W - word_shift;
 706
 707         bigint_sub_u(&x, &x, &tmp);
 708         bigint_adjust(a);
 709         return;
 710 }
 711
 712 /******************************************************************************/
 713
 714 void bigint_reduce(bigint_t *a, const bigint_t *r){
 715 //      bigint_adjust((bigint_t*)r);
 716         uint8_t rfbs = GET_FBS(r);
 717 #if DEBUG
 718         cli_putstr("\r\nDBG: (a) = "); bigint_print_hex(a);
 719 #endif
 720         if(r->length_W==0 || a->length_W==0){
 721                 return;
 722         }
 723
 724         if(bigint_length_b(a) + 3 > bigint_length_b(r)){
 725         if((r->length_W*sizeof(bigint_word_t)<=4) && (a->length_W*sizeof(bigint_word_t)<=4)){
 726             uint32_t p=0, q=0;
 727             memcpy(&p, a->wordv, a->length_W*sizeof(bigint_word_t));
 728             memcpy(&q, r->wordv, r->length_W*sizeof(bigint_word_t));
 729             p %= q;
 730             memcpy(a->wordv, &p, a->length_W*sizeof(bigint_word_t));
 731             bigint_adjust(a);
 732     //          cli_putstr("\r\nDBG: (0) = "); bigint_print_hex(a);
 733             return;
 734         }
 735         uint16_t shift;
 736         while(a->length_W > r->length_W){
 737             shift = (a->length_W - r->length_W) * 8 * sizeof(bigint_word_t) + GET_FBS(a) - rfbs - 1;
 738             /*
 739             if((a->wordv[a->length_W-1] & ((1LL<<GET_FBS(a)) - 1)) > r->wordv[r->length_W-1]){
 740                 // cli_putc('~');
 741                 cli_putstr("\r\n ~ [a] = ");
 742                 cli_hexdump_rev(&a->wordv[a->length_W-1], 4);
 743                 cli_putstr("  [r] = ");
 744                 cli_hexdump_rev(&r->wordv[r->length_W-1], 4);
 745                 shift += 1;
 746             }
 747             */
 748     //          cli_putstr("\r\nDBG: (p) shift = "); cli_hexdump_rev(&shift, 2);
 749     //          cli_putstr(" a_len = "); cli_hexdump_rev(&a->length_W, 2);
 750     //          cli_putstr(" r_len = "); cli_hexdump_rev(&r->length_W, 2);
 751     //          uart_flush(0);
 752             bigint_sub_u_bitscale(a, r, shift);
 753     //          cli_putstr("\r\nDBG: (1) = "); bigint_print_hex(a);
 754         }
 755         while((GET_FBS(a) > rfbs) && (a->length_W == r->length_W)){
 756             shift = GET_FBS(a)-rfbs-1;
 757     //          cli_putstr("\r\nDBG: (q) shift = "); cli_hexdump_rev(&shift, 2);
 758             bigint_sub_u_bitscale(a, r, shift);
 759     //          cli_putstr("\r\nDBG: (2) = "); bigint_print_hex(a);
 760         }
 761         }
 762         while(bigint_cmp_u(a,r)>=0){
 763                 bigint_sub_u(a,a,r);
 764 //              cli_putstr("\r\nDBG: (3) = "); bigint_print_hex(a);
 765         }
 766         bigint_adjust(a);
 767 //      cli_putstr("\r\nDBG: (a) = "); bigint_print_hex(a);
 768 //      cli_putstr("\r\n");
 769 }
 770
 771 /******************************************************************************/
 772
 773 /* calculate dest = a**exp % r */
 774 /* using square&multiply */
 775 void bigint_expmod_u_sam(bigint_t *dest, const bigint_t *a, const bigint_t *exp, const bigint_t *r){
 776         if(a->length_W==0 || r->length_W==0){
 777                 return;
 778         }
 779
 780         bigint_t res, base;
 781         bigint_word_t t, base_b[MAX(a->length_W,r->length_W)], res_b[r->length_W*2];
 782         uint16_t i;
 783         uint8_t j;
 784 //      uint16_t *xaddr = &i;
 785 //      cli_putstr("\r\npre-alloc (");
 786 //      cli_hexdump_rev(&xaddr, 4);
 787 //      cli_putstr(") ...");
 788         res.wordv = res_b;
 789         base.wordv = base_b;
 790         bigint_copy(&base, a);
 791 //      cli_putstr("\r\npost-copy");
 792         bigint_reduce(&base, r);
 793         res.wordv[0]=1;
 794         res.length_W=1;
 795         res.info = 0;
 796         bigint_adjust(&res);
 797         if(exp->length_W == 0){
 798                 bigint_copy(dest, &res);
 799                 return;
 800         }
 801         uint8_t flag = 0;
 802         t=exp->wordv[exp->length_W - 1];
 803         for(i=exp->length_W; i > 0; --i){
 804                 t = exp->wordv[i - 1];
 805                 for(j=BIGINT_WORD_SIZE; j > 0; --j){
 806                         if(!flag){
 807                                 if(t & (1<<(BIGINT_WORD_SIZE-1))){
 808                                         flag = 1;
 809                                 }
 810                         }
 811                         if(flag){
 812                                 bigint_square(&res, &res);
 813                                 bigint_reduce(&res, r);
 814                                 if(t & (1 << (BIGINT_WORD_SIZE - 1))){
 815                                         bigint_mul_u(&res, &res, &base);
 816                                         bigint_reduce(&res, r);
 817                                 }
 818                         }
 819                         t <<= 1;
 820                 }
 821         }
 822
 823 //      cli_putc('+');
 824         SET_POS(&res);
 825         bigint_copy(dest, &res);
 826 }
 827
 828 /******************************************************************************/
 829 #if 1
 830 #define cli_putstr(a)
 831 #define cli_putstr_P(a)
 832 #define bigint_print_hex(a)
 833 #define cli_hexdump_rev(a,b)
 834 #define uart_flush(a)
 835 #define printf_P(...)
 836 #endif
 837 /* gcd <-- gcd(x,y) a*x+b*y=gcd */
 838 void bigint_gcdext(bigint_t *gcd, bigint_t *a, bigint_t *b, const bigint_t *x, const bigint_t *y){
 839          uint16_t i = 0;
 840          printf_P(PSTR("\nDBG: gcdext( "));
 841          bigint_print_hex(x);
 842      printf_P(PSTR(", "));
 843      bigint_print_hex(y);
 844      printf_P(PSTR(")\n"));
 845          if(x->length_W == 0 || y->length_W == 0){
 846              printf_P(PSTR("\nDBG: got zero in gcd <%s %s %d>\n"), __FILE__, __func__, __LINE__);
 847              if(gcd){
 848                  bigint_set_zero(gcd);
 849              }
 850              if(a){
 851                  bigint_set_zero(a);
 852              }
 853          if(b){
 854              bigint_set_zero(b);
 855          }
 856                  return;
 857          }
 858          if(x->length_W == 1 && x->wordv[0] == 1){
 859              if(gcd){
 860              gcd->length_W = 1;
 861              gcd->wordv[0] = 1;
 862              gcd->info = 0;
 863              }
 864                  if(a){
 865                          a->length_W = 1;
 866                          a->wordv[0] = 1;
 867                          SET_POS(a);
 868                          bigint_adjust(a);
 869                  }
 870                  if(b){
 871                          bigint_set_zero(b);
 872                  }
 873                  return;
 874          }
 875          if(y->length_W == 1 && y->wordv[0] == 1){
 876                  if(gcd){
 877              gcd->length_W = 1;
 878              gcd->wordv[0] = 1;
 879              gcd->info = 0;
 880                  }
 881                  if(b){
 882                          b->length_W = 1;
 883                          b->wordv[0] = 1;
 884                          SET_POS(b);
 885                          bigint_adjust(b);
 886                  }
 887                  if(a){
 888                          bigint_set_zero(a);
 889                  }
 890                  return;
 891          }
 892
 893          while(x->wordv[i] == 0 && y->wordv[i] == 0){
 894                  ++i;
 895          }
 896          bigint_word_t g_b[i + 2], x_b[x->length_W - i], y_b[y->length_W - i];
 897          bigint_word_t u_b[x->length_W - i], v_b[y->length_W - i];
 898          bigint_word_t a_b[y->length_W + 2], c_b[y->length_W + 2];
 899          bigint_word_t b_b[x->length_W + 2], d_b[x->length_W + 2];
 900      bigint_t g, x_, y_, u, v, a_, b_, c_, d_;
 901
 902          g.wordv = g_b;
 903          x_.wordv = x_b;
 904          y_.wordv = y_b;
 905          memset(g_b, 0, i * sizeof(bigint_word_t));
 906          g_b[i] = 1;
 907          g.length_W = i + 1;
 908          g.info = 0;
 909          x_.info = y_.info = 0;
 910          x_.length_W = x->length_W - i;
 911          y_.length_W = y->length_W - i;
 912          memcpy(x_.wordv, x->wordv + i, x_.length_W * sizeof(bigint_word_t));
 913          memcpy(y_.wordv, y->wordv + i, y_.length_W * sizeof(bigint_word_t));
 914          for(i = 0; (x_.wordv[0] & (1 << i)) == 0 && (y_.wordv[0] & (1 << i)) == 0; ++i){
 915          }
 916
 917          bigint_adjust(&x_);
 918          bigint_adjust(&y_);
 919
 920          if(i){
 921                  bigint_shiftleft(&g, i);
 922                  bigint_shiftright(&x_, i);
 923                  bigint_shiftright(&y_, i);
 924          }
 925
 926          u.wordv = u_b;
 927          v.wordv = v_b;
 928          a_.wordv = a_b;
 929          b_.wordv = b_b;
 930          c_.wordv = c_b;
 931          d_.wordv = d_b;
 932
 933          bigint_copy(&u, &x_);
 934          bigint_copy(&v, &y_);
 935          a_.wordv[0] = 1;
 936          a_.length_W = 1;
 937          a_.info = 0;
 938          d_.wordv[0] = 1;
 939          d_.length_W = 1;
 940          d_.info = 0;
 941          bigint_set_zero(&b_);
 942          bigint_set_zero(&c_);
 943      printf_P(PSTR("\nloop: x_ = "));
 944      bigint_print_hex(&x_);
 945      printf_P(PSTR("; y_ = "));
 946      bigint_print_hex(&y_);
 947          do{
 948                  printf_P(PSTR("\nDBG (gcdext) 0"));
 949                  while((u.wordv[0] & 1) == 0){
 950                          printf_P(PSTR("\nDBG (gcdext) 0.1"));
 951                          bigint_shiftright(&u, 1);
 952                          if((a_.wordv[0] & 1) || (b_.wordv[0] & 1)){
 953                                  bigint_add_s(&a_, &a_, &y_);
 954                                  bigint_sub_s(&b_, &b_, &x_);
 955                          }
 956                          bigint_shiftright(&a_, 1);
 957                          bigint_shiftright(&b_, 1);
 958                          printf_P(PSTR(" a_ = "));
 959                          bigint_print_hex(&a_);
 960                          printf_P(PSTR("; b_ = "));
 961                          bigint_print_hex(&b_);
 962                  }
 963                  while((v.wordv[0] & 1) == 0){
 964                          printf_P(PSTR("\nDBG (gcdext) 0.2"));
 965                          bigint_shiftright(&v, 1);
 966                          if((c_.wordv[0] & 1) || (d_.wordv[0] & 1)){
 967                                  bigint_add_s(&c_, &c_, &y_);
 968                                  bigint_sub_s(&d_, &d_, &x_);
 969                          }
 970              printf_P(PSTR(" c* = "));
 971              bigint_print_hex(&c_);
 972                          bigint_shiftright(&c_, 1);
 973                          bigint_shiftright(&d_, 1);
 974              printf_P(PSTR(" c_ = "));
 975              bigint_print_hex(&c_);
 976              printf_P(PSTR("; d_ = "));
 977              bigint_print_hex(&d_);
 978                  }
 979                  if(bigint_cmp_u(&u, &v) >= 0){
 980              printf_P(PSTR("\nDBG (gcdext) 0.3"));
 981                         bigint_sub_u(&u, &u, &v);
 982                         bigint_sub_s(&a_, &a_, &c_);
 983                         bigint_sub_s(&b_, &b_, &d_);
 984             printf_P(PSTR(" a_ = "));
 985             bigint_print_hex(&a_);
 986             printf_P(PSTR("; b_ = "));
 987             bigint_print_hex(&b_);
 988                  }else{
 989              printf_P(PSTR("\nDBG (gcdext) 0.4"));
 990                         bigint_sub_u(&v, &v, &u);
 991                         bigint_sub_s(&c_, &c_, &a_);
 992                         bigint_sub_s(&d_, &d_, &b_);
 993             printf_P(PSTR(" c_ = "));
 994             bigint_print_hex(&c_);
 995             printf_P(PSTR("; d_ = "));
 996             bigint_print_hex(&d_);
 997                  }
 998          }while(u.length_W);
 999          if(gcd){
1000                  bigint_mul_s(gcd, &v, &g);
1001          }
1002          if(a){
1003                 bigint_copy(a, &c_);
1004          }
1005          if(b){
1006                  bigint_copy(b, &d_);
1007          }
1008 }
1009
1010 /******************************************************************************/
1011
1012 void bigint_inverse(bigint_t *dest, const bigint_t *a, const bigint_t *m){
1013         bigint_gcdext(NULL, dest, NULL, a, m);
1014         while(dest->info&BIGINT_NEG_MASK){
1015                 bigint_add_s(dest, dest, m);
1016         }
1017 }
1018
1019 /******************************************************************************/
1020
1021 void bigint_changeendianess(bigint_t *a){
1022         uint8_t t, *p, *q;
1023         p = (uint8_t*)(a->wordv);
1024         q = p + a->length_W * sizeof(bigint_word_t) - 1;
1025         while(p<q){
1026                 t = *p;
1027                 *p = *q;
1028                 *q = t;
1029                 ++p; --q;
1030         }
1031 }
1032
1033 /******************************************************************************/
1034
1035
1036
1037
1038 /******************************************************************************/
1039
1040 void bigint_expmod_u(bigint_t *dest, const bigint_t *a, const bigint_t *exp, const bigint_t *r){
1041     bigint_expmod_u_sam(dest, a, exp, r);
1042 }
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058