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 void bigint_sub_u_bitscale(bigint_t* a, const bigint_t* b, uint16_t bitscale){
 687         bigint_t tmp, x;
 688         bigint_word_t tmp_b[b->length_W + 1];
 689         const uint16_t word_shift = bitscale / BIGINT_WORD_SIZE;
 690
 691         if(a->length_W < b->length_W + word_shift){
 692 #if DEBUG
 693                 cli_putstr("\r\nDBG: *bang*\r\n");
 694 #endif
 695                 bigint_set_zero(a);
 696                 return;
 697         }
 698         tmp.wordv = tmp_b;
 699         bigint_copy(&tmp, b);
 700         bigint_shiftleft(&tmp, bitscale % BIGINT_WORD_SIZE);
 701
 702         x.info = a->info;
 703         x.wordv = &(a->wordv[word_shift]);
 704         x.length_W = a->length_W - word_shift;
 705
 706         bigint_sub_u(&x, &x, &tmp);
 707         bigint_adjust(a);
 708         return;
 709 }
 710
 711 /******************************************************************************/
 712
 713 void bigint_reduce(bigint_t* a, const bigint_t* r){
 714 //      bigint_adjust((bigint_t*)r);
 715         uint8_t rfbs = GET_FBS(r);
 716 #if DEBUG
 717         cli_putstr("\r\nDBG: (a) = "); bigint_print_hex(a);
 718 #endif
 719         if(r->length_W==0 || a->length_W==0){
 720                 return;
 721         }
 722         if((r->length_W*sizeof(bigint_word_t)<=4) && (a->length_W*sizeof(bigint_word_t)<=4)){
 723                 uint32_t p=0, q=0;
 724                 memcpy(&p, a->wordv, a->length_W*sizeof(bigint_word_t));
 725                 memcpy(&q, r->wordv, r->length_W*sizeof(bigint_word_t));
 726                 p %= q;
 727                 memcpy(a->wordv, &p, a->length_W*sizeof(bigint_word_t));
 728                 bigint_adjust(a);
 729 //              cli_putstr("\r\nDBG: (0) = "); bigint_print_hex(a);
 730                 return;
 731         }
 732         uint16_t shift;
 733         while(a->length_W > r->length_W){
 734                 shift = (a->length_W - r->length_W) * 8 * sizeof(bigint_word_t) + GET_FBS(a) - rfbs - 1;
 735                 /*
 736                 if((a->wordv[a->length_W-1] & ((1LL<<GET_FBS(a)) - 1)) > r->wordv[r->length_W-1]){
 737                         // cli_putc('~');
 738                         cli_putstr("\r\n ~ [a] = ");
 739                         cli_hexdump_rev(&a->wordv[a->length_W-1], 4);
 740                         cli_putstr("  [r] = ");
 741                         cli_hexdump_rev(&r->wordv[r->length_W-1], 4);
 742                         shift += 1;
 743                 }
 744                 */
 745 //              cli_putstr("\r\nDBG: (p) shift = "); cli_hexdump_rev(&shift, 2);
 746 //              cli_putstr(" a_len = "); cli_hexdump_rev(&a->length_W, 2);
 747 //              cli_putstr(" r_len = "); cli_hexdump_rev(&r->length_W, 2);
 748 //              uart_flush(0);
 749                 bigint_sub_u_bitscale(a, r, shift);
 750 //              cli_putstr("\r\nDBG: (1) = "); bigint_print_hex(a);
 751         }
 752         while((GET_FBS(a) > rfbs) && (a->length_W == r->length_W)){
 753                 shift = GET_FBS(a)-rfbs-1;
 754 //              cli_putstr("\r\nDBG: (q) shift = "); cli_hexdump_rev(&shift, 2);
 755                 bigint_sub_u_bitscale(a, r, shift);
 756 //              cli_putstr("\r\nDBG: (2) = "); bigint_print_hex(a);
 757         }
 758         while(bigint_cmp_u(a,r)>=0){
 759                 bigint_sub_u(a,a,r);
 760 //              cli_putstr("\r\nDBG: (3) = "); bigint_print_hex(a);
 761         }
 762         bigint_adjust(a);
 763 //      cli_putstr("\r\nDBG: (a) = "); bigint_print_hex(a);
 764 //      cli_putstr("\r\n");
 765 }
 766
 767 /******************************************************************************/
 768
 769 /* calculate dest = a**exp % r */
 770 /* using square&multiply */
 771 void bigint_expmod_u(bigint_t* dest, const bigint_t* a, const bigint_t* exp, const bigint_t* r){
 772         if(a->length_W==0 || r->length_W==0){
 773                 return;
 774         }
 775
 776         bigint_t res, base;
 777         bigint_word_t t, base_b[MAX(a->length_W,r->length_W)], res_b[r->length_W*2];
 778         uint16_t i;
 779         uint8_t j;
 780 //      uint16_t *xaddr = &i;
 781 //      cli_putstr("\r\npre-alloc (");
 782 //      cli_hexdump_rev(&xaddr, 4);
 783 //      cli_putstr(") ...");
 784         res.wordv = res_b;
 785         base.wordv = base_b;
 786         bigint_copy(&base, a);
 787 //      cli_putstr("\r\npost-copy");
 788         bigint_reduce(&base, r);
 789         res.wordv[0]=1;
 790         res.length_W=1;
 791         res.info = 0;
 792         bigint_adjust(&res);
 793         if(exp->length_W == 0){
 794                 bigint_copy(dest, &res);
 795                 return;
 796         }
 797         uint8_t flag = 0;
 798         t=exp->wordv[exp->length_W - 1];
 799         for(i=exp->length_W; i > 0; --i){
 800                 t = exp->wordv[i - 1];
 801                 for(j=BIGINT_WORD_SIZE; j > 0; --j){
 802                         if(!flag){
 803                                 if(t & (1<<(BIGINT_WORD_SIZE-1))){
 804                                         flag = 1;
 805                                 }
 806                         }
 807                         if(flag){
 808                                 bigint_square(&res, &res);
 809                                 bigint_reduce(&res, r);
 810                                 if(t & (1 << (BIGINT_WORD_SIZE - 1))){
 811                                         bigint_mul_u(&res, &res, &base);
 812                                         bigint_reduce(&res, r);
 813                                 }
 814                         }
 815                         t <<= 1;
 816                 }
 817         }
 818
 819 //      cli_putc('+');
 820         SET_POS(&res);
 821         bigint_copy(dest, &res);
 822 }
 823
 824 /******************************************************************************/
 825 #if 1
 826 #define cli_putstr(a)
 827 #define cli_putstr_P(a)
 828 #define bigint_print_hex(a)
 829 #define cli_hexdump_rev(a,b)
 830 #define uart_flush(a)
 831 #define printf_P(...)
 832 #endif
 833 /* gcd <-- gcd(x,y) a*x+b*y=gcd */
 834 void bigint_gcdext(bigint_t* gcd, bigint_t* a, bigint_t* b, const bigint_t* x, const bigint_t* y){
 835          uint16_t i = 0;
 836          printf_P(PSTR("\nDBG: gcdext( "));
 837          bigint_print_hex(x);
 838      printf_P(PSTR(", "));
 839      bigint_print_hex(y);
 840      printf_P(PSTR(")\n"));
 841          if(x->length_W == 0 || y->length_W == 0){
 842              printf_P(PSTR("\nDBG: got zero in gcd <%s %s %d>\n"), __FILE__, __func__, __LINE__);
 843              if(gcd){
 844                  bigint_set_zero(gcd);
 845              }
 846              if(a){
 847                  bigint_set_zero(a);
 848              }
 849          if(b){
 850              bigint_set_zero(b);
 851          }
 852                  return;
 853          }
 854          if(x->length_W == 1 && x->wordv[0] == 1){
 855              if(gcd){
 856              gcd->length_W = 1;
 857              gcd->wordv[0] = 1;
 858              gcd->info = 0;
 859              }
 860                  if(a){
 861                          a->length_W = 1;
 862                          a->wordv[0] = 1;
 863                          SET_POS(a);
 864                          bigint_adjust(a);
 865                  }
 866                  if(b){
 867                          bigint_set_zero(b);
 868                  }
 869                  return;
 870          }
 871          if(y->length_W == 1 && y->wordv[0] == 1){
 872                  if(gcd){
 873              gcd->length_W = 1;
 874              gcd->wordv[0] = 1;
 875              gcd->info = 0;
 876                  }
 877                  if(b){
 878                          b->length_W = 1;
 879                          b->wordv[0] = 1;
 880                          SET_POS(b);
 881                          bigint_adjust(b);
 882                  }
 883                  if(a){
 884                          bigint_set_zero(a);
 885                  }
 886                  return;
 887          }
 888
 889          while(x->wordv[i] == 0 && y->wordv[i] == 0){
 890                  ++i;
 891          }
 892          bigint_word_t g_b[i + 2], x_b[x->length_W - i], y_b[y->length_W - i];
 893          bigint_word_t u_b[x->length_W - i], v_b[y->length_W - i];
 894          bigint_word_t a_b[y->length_W + 2], c_b[y->length_W + 2];
 895          bigint_word_t b_b[x->length_W + 2], d_b[x->length_W + 2];
 896      bigint_t g, x_, y_, u, v, a_, b_, c_, d_;
 897
 898          g.wordv = g_b;
 899          x_.wordv = x_b;
 900          y_.wordv = y_b;
 901          memset(g_b, 0, i * sizeof(bigint_word_t));
 902          g_b[i] = 1;
 903          g.length_W = i + 1;
 904          g.info = 0;
 905          x_.info = y_.info = 0;
 906          x_.length_W = x->length_W - i;
 907          y_.length_W = y->length_W - i;
 908          memcpy(x_.wordv, x->wordv + i, x_.length_W * sizeof(bigint_word_t));
 909          memcpy(y_.wordv, y->wordv + i, y_.length_W * sizeof(bigint_word_t));
 910          for(i = 0; (x_.wordv[0] & (1 << i)) == 0 && (y_.wordv[0] & (1 << i)) == 0; ++i){
 911          }
 912
 913          bigint_adjust(&x_);
 914          bigint_adjust(&y_);
 915
 916          if(i){
 917                  bigint_shiftleft(&g, i);
 918                  bigint_shiftright(&x_, i);
 919                  bigint_shiftright(&y_, i);
 920          }
 921
 922          u.wordv = u_b;
 923          v.wordv = v_b;
 924          a_.wordv = a_b;
 925          b_.wordv = b_b;
 926          c_.wordv = c_b;
 927          d_.wordv = d_b;
 928
 929          bigint_copy(&u, &x_);
 930          bigint_copy(&v, &y_);
 931          a_.wordv[0] = 1;
 932          a_.length_W = 1;
 933          a_.info = 0;
 934          d_.wordv[0] = 1;
 935          d_.length_W = 1;
 936          d_.info = 0;
 937          bigint_set_zero(&b_);
 938          bigint_set_zero(&c_);
 939      printf_P(PSTR("\nloop: x_ = "));
 940      bigint_print_hex(&x_);
 941      printf_P(PSTR("; y_ = "));
 942      bigint_print_hex(&y_);
 943          do{
 944                  printf_P(PSTR("\nDBG (gcdext) 0"));
 945                  while((u.wordv[0] & 1) == 0){
 946                          printf_P(PSTR("\nDBG (gcdext) 0.1"));
 947                          bigint_shiftright(&u, 1);
 948                          if((a_.wordv[0] & 1) || (b_.wordv[0] & 1)){
 949                                  bigint_add_s(&a_, &a_, &y_);
 950                                  bigint_sub_s(&b_, &b_, &x_);
 951                          }
 952                          bigint_shiftright(&a_, 1);
 953                          bigint_shiftright(&b_, 1);
 954                          printf_P(PSTR(" a_ = "));
 955                          bigint_print_hex(&a_);
 956                          printf_P(PSTR("; b_ = "));
 957                          bigint_print_hex(&b_);
 958                  }
 959                  while((v.wordv[0] & 1) == 0){
 960                          printf_P(PSTR("\nDBG (gcdext) 0.2"));
 961                          bigint_shiftright(&v, 1);
 962                          if((c_.wordv[0] & 1) || (d_.wordv[0] & 1)){
 963                                  bigint_add_s(&c_, &c_, &y_);
 964                                  bigint_sub_s(&d_, &d_, &x_);
 965                          }
 966              printf_P(PSTR(" c* = "));
 967              bigint_print_hex(&c_);
 968                          bigint_shiftright(&c_, 1);
 969                          bigint_shiftright(&d_, 1);
 970              printf_P(PSTR(" c_ = "));
 971              bigint_print_hex(&c_);
 972              printf_P(PSTR("; d_ = "));
 973              bigint_print_hex(&d_);
 974                  }
 975                  if(bigint_cmp_u(&u, &v) >= 0){
 976              printf_P(PSTR("\nDBG (gcdext) 0.3"));
 977                         bigint_sub_u(&u, &u, &v);
 978                         bigint_sub_s(&a_, &a_, &c_);
 979                         bigint_sub_s(&b_, &b_, &d_);
 980             printf_P(PSTR(" a_ = "));
 981             bigint_print_hex(&a_);
 982             printf_P(PSTR("; b_ = "));
 983             bigint_print_hex(&b_);
 984                  }else{
 985              printf_P(PSTR("\nDBG (gcdext) 0.4"));
 986                         bigint_sub_u(&v, &v, &u);
 987                         bigint_sub_s(&c_, &c_, &a_);
 988                         bigint_sub_s(&d_, &d_, &b_);
 989             printf_P(PSTR(" c_ = "));
 990             bigint_print_hex(&c_);
 991             printf_P(PSTR("; d_ = "));
 992             bigint_print_hex(&d_);
 993                  }
 994          }while(u.length_W);
 995          if(gcd){
 996                  bigint_mul_s(gcd, &v, &g);
 997          }
 998          if(a){
 999                 bigint_copy(a, &c_);
1000          }
1001          if(b){
1002                  bigint_copy(b, &d_);
1003          }
1004 }
1005
1006 /******************************************************************************/
1007
1008 void bigint_inverse(bigint_t* dest, const bigint_t* a, const bigint_t* m){
1009         bigint_gcdext(NULL, dest, NULL, a, m);
1010         while(dest->info&BIGINT_NEG_MASK){
1011                 bigint_add_s(dest, dest, m);
1012         }
1013 }
1014
1015 /******************************************************************************/
1016
1017 void bigint_changeendianess(bigint_t* a){
1018         uint8_t t, *p, *q;
1019         p = (uint8_t*)(a->wordv);
1020         q = p + a->length_W * sizeof(bigint_word_t) - 1;
1021         while(p<q){
1022                 t = *p;
1023                 *p = *q;
1024                 *q = t;
1025                 ++p; --q;
1026         }
1027 }
1028
1029 /******************************************************************************/
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052