dsa/bigint.c

   1 /* bigint.c */
   2 /*
   3     This file is part of the AVR-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 #include "cli.h"
  37 #include "bigint_io.h"
  38 */
  39 #ifndef MAX
  40  #define MAX(a,b) (((a)>(b))?(a):(b))
  41 #endif
  42
  43 #ifndef MIN
  44  #define MIN(a,b) (((a)<(b))?(a):(b))
  45 #endif
  46
  47 #define SET_FBS(a, v) do{(a)->info &=0xF8; (a)->info |= (v);}while(0)
  48 #define GET_FBS(a)   ((a)->info&BIGINT_FBS_MASK)
  49 #define SET_NEG(a)   (a)->info |= BIGINT_NEG_MASK
  50 #define SET_POS(a)   (a)->info &= ~BIGINT_NEG_MASK
  51 #define XCHG(a,b)    do{(a)^=(b); (b)^=(a); (a)^=(b);}while(0)
  52 #define XCHG_PTR(a,b)    do{ a = (void*)(((uint16_t)(a)) ^ ((uint16_t)(b))); \
  53                                  b = (void*)(((uint16_t)(a)) ^ ((uint16_t)(b))); \
  54                                  a = (void*)(((uint16_t)(a)) ^ ((uint16_t)(b)));}while(0)
  55
  56 #define GET_SIGN(a) ((a)->info&BIGINT_NEG_MASK)
  57
  58 /******************************************************************************/
  59 void bigint_adjust(bigint_t* a){
  60         while(a->length_B!=0 && a->wordv[a->length_B-1]==0){
  61                 a->length_B--;
  62         }
  63         if(a->length_B==0){
  64                 a->info=0;
  65                 return;
  66         }
  67         uint8_t t;
  68         uint8_t i = 0x07;
  69         t = a->wordv[a->length_B-1];
  70         while((t&0x80)==0 && i){
  71                 t<<=1;
  72                 i--;
  73         }
  74         SET_FBS(a, i);
  75 }
  76
  77 /******************************************************************************/
  78
  79 void bigint_copy(bigint_t* dest, const bigint_t* src){
  80         memcpy(dest->wordv, src->wordv, src->length_B);
  81         dest->length_B = src->length_B;
  82         dest->info = src->info;
  83 }
  84
  85 /******************************************************************************/
  86
  87 /* this should be implemented in assembly */
  88 /*
  89 void bigint_add_u(bigint_t* dest, const bigint_t* a, const bigint_t* b){
  90         uint16_t t=0, i;
  91         if(a->length_B < b->length_B){
  92                 XCHG_PTR(a,b);
  93         }
  94         for(i=0; i<b->length_B; ++i){
  95                 t = a->wordv[i] + b->wordv[i] + t;
  96                 dest->wordv[i] = (uint8_t)t;
  97                 t>>=8;
  98         }
  99         for(; i<a->length_B; ++i){
 100                 t = a->wordv[i] + t;
 101                 dest->wordv[i] = (uint8_t)t;
 102                 t>>=8;
 103         }
 104         dest->wordv[i++] = t;
 105         dest->length_B = i;
 106         bigint_adjust(dest);
 107 }
 108 */
 109 /******************************************************************************/
 110
 111 /* this should be implemented in assembly */
 112 void bigint_add_scale_u(bigint_t* dest, const bigint_t* a, uint16_t scale){
 113         uint16_t i,j=0;
 114         uint16_t t=0;
 115         if(scale>dest->length_B)
 116                 memset(dest->wordv+dest->length_B, 0, scale-dest->length_B);
 117         for(i=scale; i<a->length_B+scale; ++i,++j){
 118                 t = a->wordv[j] + t;
 119                 if(dest->length_B>i){
 120                         t += dest->wordv[i];
 121                 }
 122                 dest->wordv[i] = (uint8_t)t;
 123                 t>>=8;
 124         }
 125         while(t){
 126                 if(dest->length_B>i){
 127                         t = dest->wordv[i] + t;
 128                 }
 129                 dest->wordv[i] = (uint8_t)t;
 130                 t>>=8;
 131                 ++i;
 132         }
 133         if(dest->length_B < i){
 134                 dest->length_B = i;
 135         }
 136         bigint_adjust(dest);
 137 }
 138
 139 /******************************************************************************/
 140
 141 /* this should be implemented in assembly */
 142 void bigint_sub_u(bigint_t* dest, const bigint_t* a, const bigint_t* b){
 143         int8_t borrow=0;
 144         int8_t  r;
 145         int16_t t;
 146         uint16_t i, min, max;
 147         min = MIN(a->length_B, b->length_B);
 148         max = MAX(a->length_B, b->length_B);
 149         r = bigint_cmp_u(a,b);
 150         if(r==0){
 151                 dest->length_B = 0;
 152                 dest->wordv[0] = 0;
 153                 bigint_adjust(dest);
 154                 return;
 155         }
 156         if(b->length_B==0){
 157                 dest->length_B = a->length_B;
 158                 memcpy(dest->wordv, a->wordv, a->length_B);
 159                 dest->info = a->info;
 160                 SET_POS(dest);
 161                 return;
 162         }
 163         if(a->length_B==0){
 164                         dest->length_B = b->length_B;
 165                         memcpy(dest->wordv, b->wordv, b->length_B);
 166                         dest->info = b->info;
 167                         SET_NEG(dest);
 168                         return;
 169         }
 170         if(r<0){
 171                 bigint_sub_u(dest, b, a);
 172                 SET_NEG(dest);
 173         }else{
 174                 for(i=0; i<min; ++i){
 175                         t = a->wordv[i] - b->wordv[i] - borrow;
 176                         if(t<0){
 177                                 borrow = 1;
 178                                 dest->wordv[i]=(uint8_t)t;
 179                         }else{
 180                                 borrow = 0;
 181                                 dest->wordv[i]=(uint8_t)t;
 182                         }
 183                 }
 184                 for(;i<max; ++i){
 185                         t = a->wordv[i] - borrow;
 186                         if(t<0){
 187                                 borrow = 1;
 188                                 dest->wordv[i]=(uint8_t)t;
 189                         }else{
 190                                 borrow = 0;
 191                                 dest->wordv[i]=(uint8_t)t;
 192                         }
 193
 194                 }
 195                 SET_POS(dest);
 196                 dest->length_B = i;
 197                 bigint_adjust(dest);
 198         }
 199 }
 200
 201 /******************************************************************************/
 202
 203 int8_t bigint_cmp_u(const bigint_t* a, const bigint_t* b){
 204         if(a->length_B > b->length_B){
 205                 return 1;
 206         }
 207         if(a->length_B < b->length_B){
 208                 return -1;
 209         }
 210         if(a->length_B==0){
 211                 return 0;
 212         }
 213         uint16_t i;
 214         i = a->length_B-1;
 215         do{
 216                 if(a->wordv[i]!=b->wordv[i]){
 217                         if(a->wordv[i]>b->wordv[i]){
 218                                 return 1;
 219                         }else{
 220                                 return -1;
 221                         }
 222                 }
 223         }while(i--);
 224         return 0;
 225 }
 226
 227 /******************************************************************************/
 228
 229 void bigint_add_s(bigint_t* dest, const bigint_t* a, const bigint_t* b){
 230         uint8_t s;
 231         s  = GET_SIGN(a)?2:0;
 232         s |= GET_SIGN(b)?1:0;
 233         switch(s){
 234                 case 0: /* both positive */
 235                         bigint_add_u(dest, a,b);
 236                         SET_POS(dest);
 237                         break;
 238                 case 1: /* a positive, b negative */
 239                         bigint_sub_u(dest, a, b);
 240                         break;
 241                 case 2: /* a negative, b positive */
 242                         bigint_sub_u(dest, b, a);
 243                         break;
 244                 case 3: /* both negative */
 245                         bigint_add_u(dest, a, b);
 246                         SET_NEG(dest);
 247                         break;
 248                 default: /* how can this happen?*/
 249                         break;
 250         }
 251 }
 252
 253 /******************************************************************************/
 254
 255 void bigint_sub_s(bigint_t* dest, const bigint_t* a, const bigint_t* b){
 256         uint8_t s;
 257         s  = GET_SIGN(a)?2:0;
 258         s |= GET_SIGN(b)?1:0;
 259         switch(s){
 260                 case 0: /* both positive */
 261                         bigint_sub_u(dest, a,b);
 262                         break;
 263                 case 1: /* a positive, b negative */
 264                         bigint_add_u(dest, a, b);
 265                         SET_POS(dest);
 266                         break;
 267                 case 2: /* a negative, b positive */
 268                         bigint_add_u(dest, a, b);
 269                         SET_NEG(dest);
 270                         break;
 271                 case 3: /* both negative */
 272                         bigint_sub_u(dest, b, a);
 273                         break;
 274                 default: /* how can this happen?*/
 275                                         break;
 276         }
 277
 278 }
 279
 280 /******************************************************************************/
 281
 282 int8_t bigint_cmp_s(const bigint_t* a, const bigint_t* b){
 283         uint8_t s;
 284         if(a->length_B==0 && b->length_B==0){
 285                 return 0;
 286         }
 287         s  = GET_SIGN(a)?2:0;
 288         s |= GET_SIGN(b)?1:0;
 289         switch(s){
 290                 case 0: /* both positive */
 291                         return bigint_cmp_u(a, b);
 292                         break;
 293                 case 1: /* a positive, b negative */
 294                         return 1;
 295                         break;
 296                 case 2: /* a negative, b positive */
 297                         return -1;
 298                         break;
 299                 case 3: /* both negative */
 300                         return bigint_cmp_u(b, a);
 301                         break;
 302                 default: /* how can this happen?*/
 303                                         break;
 304         }
 305         return 0; /* just to satisfy the compiler */
 306 }
 307
 308 /******************************************************************************/
 309
 310 void bigint_shiftleft(bigint_t* a, uint16_t shift){
 311         uint16_t byteshift;
 312         uint16_t i;
 313         uint8_t bitshift;
 314         uint16_t t=0;
 315         byteshift = (shift+3)/8;
 316         bitshift = shift&7;
 317         memmove(a->wordv+byteshift, a->wordv, a->length_B);
 318         memset(a->wordv, 0, byteshift);
 319         if(bitshift!=0){
 320                 if(bitshift<=4){ /* shift to the left */
 321                         for(i=byteshift; i<a->length_B+byteshift; ++i){
 322                                 t |= (a->wordv[i])<<bitshift;
 323                                 a->wordv[i] = (uint8_t)t;
 324                                 t >>= 8;
 325                         }
 326                         a->wordv[i] = (uint8_t)t;
 327                         byteshift++;
 328                 }else{ /* shift to the right */
 329                         for(i=a->length_B+byteshift-1; i>byteshift-1; --i){
 330                                 t |= (a->wordv[i])<<(bitshift);
 331                                 a->wordv[i] = (uint8_t)(t>>8);
 332                                 t <<= 8;
 333                         }
 334                         t |= (a->wordv[i])<<(bitshift);
 335                         a->wordv[i] = (uint8_t)(t>>8);
 336                 }
 337         }
 338         a->length_B += byteshift;
 339         bigint_adjust(a);
 340 }
 341
 342 /******************************************************************************/
 343
 344 void bigint_shiftright(bigint_t* a, uint16_t shift){
 345         uint16_t byteshift;
 346         uint16_t i;
 347         uint8_t bitshift;
 348         uint16_t t=0;
 349         byteshift = shift/8;
 350         bitshift = shift&7;
 351         if(byteshift >= a->length_B){ /* we would shift out more than we have */
 352                 bigint_set_zero(a);
 353                 return;
 354         }
 355         if(byteshift == a->length_B-1 && bitshift>GET_FBS(a)){
 356                 bigint_set_zero(a);
 357                 return;
 358         }
 359         if(byteshift){
 360                 memmove(a->wordv, a->wordv+byteshift, a->length_B-byteshift);
 361                 memset(a->wordv+a->length_B-byteshift, 0,  byteshift);
 362         }
 363         if(bitshift!=0){
 364          /* shift to the right */
 365                 for(i=a->length_B-byteshift-1; i>0; --i){
 366                         t |= (a->wordv[i])<<(8-bitshift);
 367                         a->wordv[i] = (uint8_t)(t>>8);
 368                         t <<= 8;
 369                 }
 370                 t |= (a->wordv[0])<<(8-bitshift);
 371                 a->wordv[0] = (uint8_t)(t>>8);
 372         }
 373         a->length_B -= byteshift;
 374         bigint_adjust(a);
 375 }
 376
 377 /******************************************************************************/
 378
 379 void bigint_xor(bigint_t* dest, const bigint_t* a){
 380         uint16_t i;
 381         for(i=0; i<a->length_B; ++i){
 382                 dest->wordv[i] ^= a->wordv[i];
 383         }
 384         bigint_adjust(dest);
 385 }
 386
 387 /******************************************************************************/
 388
 389 void bigint_set_zero(bigint_t* a){
 390         a->length_B=0;
 391 }
 392
 393 /******************************************************************************/
 394
 395 /* using the Karatsuba-Algorithm */
 396 /* x*y = (xh*yh)*b**2n + ((xh+xl)*(yh+yl) - xh*yh - xl*yl)*b**n + yh*yl */
 397 void bigint_mul_u(bigint_t* dest, const bigint_t* a, const bigint_t* b){
 398         if(a->length_B==0 || b->length_B==0){
 399                 bigint_set_zero(dest);
 400                 return;
 401         }
 402         if(dest==a || dest==b){
 403                 bigint_t d;
 404                 uint8_t d_b[a->length_B+b->length_B];
 405                 d.wordv = d_b;
 406                 bigint_mul_u(&d, a, b);
 407                 bigint_copy(dest, &d);
 408                 return;
 409         }
 410         if(a->length_B==1 || b->length_B==1){
 411                 if(a->length_B!=1){
 412                         XCHG_PTR(a,b);
 413                 }
 414                 uint16_t i, t=0;
 415                 uint8_t x = a->wordv[0];
 416                 for(i=0; i<b->length_B; ++i){
 417                         t += b->wordv[i]*x;
 418                         dest->wordv[i] = (uint8_t)t;
 419                         t>>=8;
 420                 }
 421                 dest->wordv[i] = (uint8_t)t;
 422                 dest->length_B=i+1;
 423                 bigint_adjust(dest);
 424                 return;
 425         }
 426         if(a->length_B<=4 && b->length_B<=4){
 427                 uint32_t p=0, q=0;
 428                 uint64_t r;
 429                 memcpy(&p, a->wordv, a->length_B);
 430                 memcpy(&q, b->wordv, b->length_B);
 431                 r = (uint64_t)p*(uint64_t)q;
 432                 memcpy(dest->wordv, &r, a->length_B+b->length_B);
 433                 dest->length_B =  a->length_B+b->length_B;
 434                 bigint_adjust(dest);
 435                 return;
 436         }
 437         bigint_set_zero(dest);
 438         /* split a in xh & xl; split b in yh & yl */
 439         uint16_t n;
 440         n=(MAX(a->length_B, b->length_B)+1)/2;
 441         bigint_t xl, xh, yl, yh;
 442         xl.wordv = a->wordv;
 443         yl.wordv = b->wordv;
 444         if(a->length_B<=n){
 445                 xh.info=0;
 446                 xh.length_B = 0;
 447                 xl.length_B = a->length_B;
 448                 xl.info = 0;
 449         }else{
 450                 xl.length_B=n;
 451                 xl.info = 0;
 452                 bigint_adjust(&xl);
 453                 xh.wordv = a->wordv+n;
 454                 xh.length_B = a->length_B-n;
 455                 xh.info = 0;
 456         }
 457         if(b->length_B<=n){
 458                 yh.info=0;
 459                 yh.length_B = 0;
 460                 yl.length_B = b->length_B;
 461                 yl.info = b->info;
 462         }else{
 463                 yl.length_B=n;
 464                 yl.info = 0;
 465                 bigint_adjust(&yl);
 466                 yh.wordv = b->wordv+n;
 467                 yh.length_B = b->length_B-n;
 468                 yh.info = 0;
 469         }
 470         /* now we have split up a and b */
 471         uint8_t  tmp_b[2*n+2], m_b[2*(n+1)];
 472         bigint_t tmp, tmp2, m;
 473         tmp.wordv = tmp_b;
 474         tmp2.wordv = tmp_b+n+1;
 475         m.wordv = m_b;
 476
 477         bigint_mul_u(dest, &xl, &yl);  /* dest <= xl*yl     */
 478         bigint_add_u(&tmp2, &xh, &xl); /* tmp2 <= xh+xl     */
 479         bigint_add_u(&tmp, &yh, &yl);  /* tmp  <= yh+yl     */
 480         bigint_mul_u(&m, &tmp2, &tmp); /* m    <= tmp2*tmp  */
 481         bigint_mul_u(&tmp, &xh, &yh);  /* h    <= xh*yh     */
 482         bigint_sub_u(&m, &m, dest);    /* m    <= m-dest    */
 483     bigint_sub_u(&m, &m, &tmp);    /* m    <= m-h       */
 484         bigint_add_scale_u(dest, &m, n);
 485         bigint_add_scale_u(dest, &tmp, 2*n);
 486 }
 487
 488 /******************************************************************************/
 489
 490 void bigint_mul_s(bigint_t* dest, const bigint_t* a, const bigint_t* b){
 491         uint8_t s;
 492         s  = GET_SIGN(a)?2:0;
 493         s |= GET_SIGN(b)?1:0;
 494         switch(s){
 495                 case 0: /* both positive */
 496                         bigint_mul_u(dest, a,b);
 497                         SET_POS(dest);
 498                         break;
 499                 case 1: /* a positive, b negative */
 500                         bigint_mul_u(dest, a,b);
 501                         SET_NEG(dest);
 502                         break;
 503                 case 2: /* a negative, b positive */
 504                         bigint_mul_u(dest, a,b);
 505                         SET_NEG(dest);
 506                         break;
 507                 case 3: /* both negative */
 508                         bigint_mul_u(dest, a,b);
 509                         SET_POS(dest);
 510                         break;
 511                 default: /* how can this happen?*/
 512                         break;
 513         }
 514 }
 515
 516 /******************************************************************************/
 517
 518 /* square */
 519 /* (xh*b^n+xl)^2 = xh^2*b^2n + 2*xh*xl*b^n + xl^2 */
 520 void bigint_square(bigint_t* dest, const bigint_t* a){
 521         if(a->length_B<=4){
 522                 uint64_t r=0;
 523                 memcpy(&r, a->wordv, a->length_B);
 524                 r = r*r;
 525                 memcpy(dest->wordv, &r, 2*a->length_B);
 526                 SET_POS(dest);
 527                 dest->length_B=2*a->length_B;
 528                 bigint_adjust(dest);
 529                 return;
 530         }
 531         if(dest==a){
 532                 bigint_t d;
 533                 uint8_t d_b[a->length_B*2];
 534                 d.wordv = d_b;
 535                 bigint_square(&d, a);
 536                 bigint_copy(dest, &d);
 537                 return;
 538         }
 539         uint16_t n;
 540         n=(a->length_B+1)/2;
 541         bigint_t xh, xl, tmp; /* x-high, x-low, temp */
 542         uint8_t buffer[2*n+1];
 543         xl.wordv = a->wordv;
 544         xl.length_B = n;
 545         xh.wordv = a->wordv+n;
 546         xh.length_B = a->length_B-n;
 547         tmp.wordv = buffer;
 548         bigint_square(dest, &xl);
 549         bigint_square(&tmp, &xh);
 550         bigint_add_scale_u(dest, &tmp, 2*n);
 551         bigint_mul_u(&tmp, &xl, &xh);
 552         bigint_shiftleft(&tmp, 1);
 553         bigint_add_scale_u(dest, &tmp, n);
 554 }
 555
 556 /******************************************************************************/
 557
 558 void bigint_sub_u_bitscale(bigint_t* a, const bigint_t* b, uint16_t bitscale){
 559         bigint_t tmp;
 560         uint8_t tmp_b[b->length_B+1];
 561         uint16_t i,j,byteshift=bitscale/8;
 562         uint8_t borrow=0;
 563         int16_t t;
 564
 565         if(a->length_B < b->length_B+byteshift){
 566                 bigint_set_zero(a);
 567                 return;
 568         }
 569
 570         tmp.wordv = tmp_b;
 571         bigint_copy(&tmp, b);
 572         bigint_shiftleft(&tmp, bitscale&7);
 573
 574         for(j=0,i=byteshift; i<tmp.length_B+byteshift; ++i, ++j){
 575                 t = a->wordv[i] - tmp.wordv[j] - borrow;
 576                 a->wordv[i] = (uint8_t)t;
 577                 if(t<0){
 578                         borrow = 1;
 579                 }else{
 580                         borrow = 0;
 581                 }
 582         }
 583         while(borrow){
 584                 if(i+1 > a->length_B){
 585                         bigint_set_zero(a);
 586                         return;
 587                 }
 588                 a->wordv[i] -= borrow;
 589                 if(a->wordv[i]!=0xff){
 590                         borrow=0;
 591                 }
 592                 ++i;
 593         }
 594         bigint_adjust(a);
 595 }
 596
 597 /******************************************************************************/
 598
 599 void bigint_reduce(bigint_t* a, const bigint_t* r){
 600 //      bigint_adjust(r);
 601         uint8_t rfbs = GET_FBS(r);
 602
 603         if(r->length_B==0 || a->length_B==0){
 604                 return;
 605         }
 606         while(a->length_B > r->length_B){
 607                 bigint_sub_u_bitscale(a, r, (a->length_B-r->length_B)*8+GET_FBS(a)-rfbs-1);
 608         }
 609         while((GET_FBS(a) > rfbs+1) && (a->length_B == r->length_B)){
 610                 bigint_sub_u_bitscale(a, r, GET_FBS(a)-rfbs-1);
 611         }
 612         while(bigint_cmp_u(a,r)>=0){
 613                 bigint_sub_u(a,a,r);
 614         }
 615         bigint_adjust(a);
 616 }
 617
 618 /******************************************************************************/
 619
 620 /* calculate dest = a**exp % r */
 621 /* using square&multiply */
 622 void bigint_expmod_u(bigint_t* dest, const bigint_t* a, const bigint_t* exp, const bigint_t* r){
 623         if(a->length_B==0 || r->length_B==0){
 624                 return;
 625         }
 626
 627         bigint_t res, base;
 628         uint8_t base_b[MAX(a->length_B,r->length_B*2)], res_b[r->length_B*2];
 629         uint16_t i;
 630         uint8_t j, t;
 631         res.wordv = res_b;
 632         base.wordv = base_b;
 633         bigint_copy(&base, a);
 634         bigint_reduce(&base, r);
 635         res.wordv[0]=1;
 636         res.length_B=1;
 637         res.info = 0;
 638         bigint_adjust(&res);
 639         for(i=0; i+1<exp->length_B; ++i){
 640                 t=exp->wordv[i];
 641                 for(j=0; j<8; ++j){
 642                         if(t&1){
 643                                 bigint_mul_u(&res, &res, &base);
 644                                 bigint_reduce(&res, r);
 645                         }
 646                         bigint_square(&base, &base);
 647                         bigint_reduce(&base, r);
 648                         t>>=1;
 649                 }
 650         }
 651         t=exp->wordv[i];
 652         while(t){
 653                 if(t&1){
 654                         bigint_mul_u(&res, &res, &base);
 655                         bigint_reduce(&res, r);
 656                 }
 657                 bigint_square(&base, &base);
 658                 bigint_reduce(&base, r);
 659                 t>>=1;
 660         }
 661         SET_POS(&res);
 662         bigint_copy(dest, &res);
 663 }
 664
 665 /******************************************************************************/
 666 /* gcd <-- gcd(x,y) a*x+b*y=gcd */
 667 void bigint_gcdext(bigint_t* gcd, bigint_t* a, bigint_t* b, const bigint_t* x, const bigint_t* y){
 668          bigint_t g, x_, y_, u, v, a_, b_, c_, d_;
 669          volatile uint16_t i=0;
 670          if(x->length_B==0 || y->length_B==0){
 671                  return;
 672          }
 673          while(x->wordv[i]==0 && y->wordv[i]==0){
 674                  ++i;
 675          }
 676          uint8_t g_b[i+2], x_b[x->length_B-i], y_b[y->length_B-i];
 677          uint8_t u_b[x->length_B-i], v_b[y->length_B-i];
 678          uint8_t a_b[y->length_B+2], c_b[y->length_B+2];
 679          uint8_t b_b[x->length_B+2], d_b[x->length_B+2];
 680
 681          g.wordv = g_b;
 682          x_.wordv = x_b;
 683          y_.wordv = y_b;
 684          memset(g_b, 0, i);
 685          g_b[i]=1;
 686          g.length_B = i+1;
 687          g.info=0;
 688          x_.info = y_.info = 0;
 689          x_.length_B = x->length_B-i;
 690          y_.length_B = y->length_B-i;
 691          memcpy(x_.wordv, x->wordv+i, x_.length_B);
 692          memcpy(y_.wordv, y->wordv+i, y_.length_B);
 693          for(i=0; (x_.wordv[0]&(1<<i))==0 && (y_.wordv[0]&(1<<i))==0; ++i){
 694          }
 695
 696          bigint_adjust(&x_);
 697          bigint_adjust(&y_);
 698
 699          if(i){
 700                  bigint_shiftleft(&g, i);
 701                  bigint_shiftright(&x_, i);
 702                  bigint_shiftright(&y_, i);
 703          }
 704          u.wordv = u_b;
 705          v.wordv = v_b;
 706          a_.wordv = a_b;
 707          b_.wordv = b_b;
 708          c_.wordv = c_b;
 709          d_.wordv = d_b;
 710
 711          bigint_copy(&u, &x_);
 712          bigint_copy(&v, &y_);
 713          a_.wordv[0] = 1;
 714          a_.length_B = 1;
 715          a_.info = 0;
 716          d_.wordv[0] = 1;
 717          d_.length_B = 1;
 718          d_.info = 0;
 719          bigint_set_zero(&b_);
 720          bigint_set_zero(&c_);
 721          do{
 722                  while((u.wordv[0]&1)==0){
 723                          bigint_shiftright(&u, 1);
 724                          if((a_.wordv[0]&1) || (b_.wordv[0]&1)){
 725                                  bigint_add_s(&a_, &a_, &y_);
 726                                  bigint_sub_s(&b_, &b_, &x_);
 727                          }
 728                          bigint_shiftright(&a_, 1);
 729                          bigint_shiftright(&b_, 1);
 730                  }
 731                  while((v.wordv[0]&1)==0){
 732                          bigint_shiftright(&v, 1);
 733                          if((c_.wordv[0]&1) || (d_.wordv[0]&1)){
 734                                  bigint_add_s(&c_, &c_, &y_);
 735                                  bigint_sub_s(&d_, &d_, &x_);
 736                          }
 737                          bigint_shiftright(&c_, 1);
 738                          bigint_shiftright(&d_, 1);
 739
 740                  }
 741                  if(bigint_cmp_u(&u, &v)>=0){
 742                         bigint_sub_u(&u, &u, &v);
 743                         bigint_sub_s(&a_, &a_, &c_);
 744                         bigint_sub_s(&b_, &b_, &d_);
 745                  }else{
 746                         bigint_sub_u(&v, &v, &u);
 747                         bigint_sub_s(&c_, &c_, &a_);
 748                         bigint_sub_s(&d_, &d_, &b_);
 749                  }
 750          }while(u.length_B);
 751          if(gcd){
 752                  bigint_mul_s(gcd, &v, &g);
 753          }
 754          if(a){
 755                 bigint_copy(a, &c_);
 756          }
 757          if(b){
 758                  bigint_copy(b, &d_);
 759          }
 760 }
 761
 762 /******************************************************************************/
 763
 764 void bigint_inverse(bigint_t* dest, const bigint_t* a, const bigint_t* m){
 765         bigint_gcdext(NULL, dest, NULL, a, m);
 766         while(dest->info&BIGINT_NEG_MASK){
 767                 bigint_add_s(dest, dest, m);
 768         }
 769 }
 770
 771 /******************************************************************************/
 772
 773 void bigint_changeendianess(bigint_t* a){
 774         uint8_t t, *p, *q;
 775         p = a->wordv;
 776         q = p+a->length_B-1;
 777         while(p<q){
 778                 t = *p;
 779                 *p = *q;
 780                 *q = t;
 781                 ++p; --q;
 782         }
 783 }
 784
 785 /******************************************************************************/
 786
 787
 788
 789
 790
 791
 792
 793
 794
 795
 796
 797
 798
 799
 800
 801
 802
 803
 804
 805
 806
 807
 808