bigint2/bigint2.c

   1 /*
   2  * bigint2.c
   3  *
   4  *  Created on: 09.04.2014
   5  *      Author: bg
   6  */
   7
   8 #include <bigint2.h>
   9 #include <bigint2_io.h>
  10 #include <stdio.h>
  11 #include <string.h>
  12
  13 #define CHECKS 1
  14
  15 #define E_PARAM -1;
  16 #define E_OOM   -2;
  17 #define E_VALUE -3;
  18
  19
  20 #define MAX(a,b) ((a) > (b) ? (a) : (b))
  21 #define MIN(a,b) ((a) < (b) ? (a) : (b))
  22
  23 void *(*int_realloc)(void *ptr, size_t size);
  24 void (*int_free)(void *ptr);
  25
  26 /**
  27  * \brief used to check if bigint is large enough
  28  * Checks if the memory of a is large enough to hold min_size
  29  * words, and if this is not the case tries to increase the allocated memory
  30  * amount (conserving the memory content).
  31  * If the allocated memory insufficient and the buffer could not be enlarged
  32  * E_OOM is returned, otherwise 0 is returned.
  33  */
  34 static
  35 int check_size(bigint_t *a, bigint_length_t min_size) {
  36     bigint_word_t *tmp;
  37     if (a->allocated_W >= min_size) {
  38         return 0;
  39     }
  40     tmp = int_realloc(a->wordv, min_size * sizeof(bigint_word_t));
  41     if (!tmp) {
  42         return E_OOM;
  43     }
  44     a->wordv = tmp;
  45     a->allocated_W = min_size;
  46     return 0;
  47 }
  48
  49 int bigint_copy(bigint_t *dest, const bigint_t *a) {
  50     int r;
  51     if (dest == a) {
  52         return 0;
  53     }
  54     r = check_size(dest, a->length_W);
  55     if (r) {
  56         return r;
  57     }
  58     dest->info = a->info;
  59     dest->length_W = a->length_W;
  60     memcpy(dest->wordv, a->wordv, a->length_W * sizeof(bigint_word_t));
  61
  62     return 0;
  63 }
  64
  65 int bigint_free(bigint_t *a) {
  66     if(a->allocated_W) {
  67         int_free(a->wordv);
  68     }
  69     memset(a, 0, sizeof(bigint_t));
  70     return 0;
  71 }
  72
  73 /**
  74  * \brief dest = |a| + |b|
  75  * Adds the bigints a and b and stores the result in dest.
  76  * Signs are ignored.
  77  */
  78 int bigint_add_u(bigint_t *dest, const bigint_t *a, const bigint_t *b) {
  79     bigint_word_t t, c1, c2, c = 0;
  80     bigint_length_t i, j;
  81     int r;
  82 #if CHECKS
  83     if (!dest || !a || !b) {
  84         return E_PARAM;
  85     }
  86 #endif
  87 #if CHECKS
  88     if ((r = check_size(dest, a->length_W))) {
  89         return r;
  90     }
  91 #endif
  92     j = MIN(a->length_W, b->length_W);
  93     for (i = 0; i < j; ++i) {
  94         t = a->wordv[i] + b->wordv[i];
  95         c1 = t < a->wordv[i];
  96         dest->wordv[i] = t + c;
  97         c2 = dest->wordv[i] < t;
  98         c = c1 | c2;
  99     }
 100     for (; i < a->length_W; ++i) {
 101         t = a->wordv[i];
 102         dest->wordv[i] = t + c;
 103         c = dest->wordv[i] < t;
 104     }
 105     dest->length_W = i;
 106     return 0;
 107 }
 108
 109 /** UNSAFE!!!
 110  * \brief dest = |a| + |b|
 111  * Adds the bigints a and b and stores the result in dest.
 112  * Signs are ignored.
 113  */
 114 int bigint_add_auto_u(bigint_t *dest, const bigint_t *a, const bigint_t *b) {
 115     bigint_word_t t, c1, c2, c = 0;
 116     bigint_length_t i, j;
 117     int r;
 118 #if CHECKS
 119     if (!dest || !a || !b) {
 120         return E_PARAM;
 121     }
 122 #endif
 123     if (a->length_W < b->length_W) {
 124         const bigint_t *t;
 125         t = a;
 126         a = b;
 127         b = t;
 128     }
 129 #if CHECKS
 130     if ((r = check_size(dest, a->length_W + 1))) {
 131         return r;
 132     }
 133 #endif
 134     j = MIN(a->length_W, b->length_W);
 135     for (i = 0; i < j; ++i) {
 136         t = a->wordv[i] + b->wordv[i];
 137         c1 = t < a->wordv[i];
 138         dest->wordv[i] = t + c;
 139         c2 = dest->wordv[i] < t;
 140         c = c1 | c2;
 141     }
 142     for (; i < a->length_W; ++i) {
 143         t = a->wordv[i];
 144         dest->wordv[i] = t + c;
 145         c = dest->wordv[i] < t;
 146     }
 147     if (c) {
 148         dest->wordv[i++] = c;
 149     }
 150     dest->length_W = i;
 151     return 0;
 152 }
 153
 154
 155 /**
 156  * \brief dest = ~a
 157  */
 158 int bigint_invert(bigint_t *dest, const bigint_t *a) {
 159     bigint_length_t i;
 160     int r;
 161 #if CHECKS
 162     if ((r = check_size(dest, a->length_W))) {
 163         return r;
 164     }
 165 #endif
 166     i = a->length_W;
 167     while(i--) {
 168         dest->wordv[i] = ~a->wordv[i];
 169     }
 170     dest->length_W = a->length_W;
 171     return 0;
 172 }
 173
 174 /**
 175  * \brief dest = a + 1
 176  */
 177 int bigint_add_one(bigint_t *dest, const bigint_t *a) {
 178     bigint_t one = {
 179         .wordv = (bigint_word_t*)"\x01",
 180         .length_W = 1,
 181         .allocated_W = 1,
 182         .info = 0
 183     };
 184     return bigint_add_u(dest, a, &one);
 185 }
 186
 187 /**
 188  * \brief dest = -a
 189  */
 190 int bigint_negate(bigint_t *dest, const bigint_t *a) {
 191     int r = 0;
 192 #if CHECKS
 193     if ((r = check_size(dest, a->length_W))) {
 194         return r;
 195     }
 196 #endif
 197     r |= bigint_invert(dest, a);
 198     r |= bigint_add_one(dest, dest);
 199     return r;
 200 }
 201
 202 /**
 203  * \brief dest = |a| - |b|
 204  * Subtracts b from a and stores the result in dest.
 205  * Signs are ignored
 206  */
 207 int bigint_sub_u(bigint_t *dest, const bigint_t *a, const bigint_t *b) {
 208     bigint_length_t i, j;
 209     bigint_word_t t, c1, c2, c = 0;
 210     int r;
 211 #if CHECKS
 212     if (!dest || !a || !b) {
 213         return E_PARAM;
 214     }
 215 #endif
 216 #if CHECKS
 217     if ((r = check_size(dest, MAX(a->length_W, b->length_W)))) {
 218         return r;
 219     }
 220
 221 #endif
 222     j = MIN(a->length_W, b->length_W);
 223     for (i = 0; i < j; ++i) {
 224         t = a->wordv[i] - b->wordv[i];
 225         c1 = t > a->wordv[i];
 226         dest->wordv[i] = t - c;
 227         c2 = dest->wordv[i] > t;
 228         c = c1 | c2;
 229     }
 230     for (; i < a->length_W; ++i) {
 231         t = a->wordv[i];
 232         dest->wordv[i] = t - c;
 233         c = dest->wordv[i] > t;
 234     }
 235     dest->length_W = a->length_W;
 236     return 0;
 237 }
 238
 239 /**
 240  * \brief a <<= 1
 241  */
 242 int bigint_shiftleft_words(bigint_t *dest, const bigint_t *a, bigint_length_t s) {
 243 #if CHECKS
 244     if (!a) {
 245         return E_PARAM;
 246     }
 247 #endif
 248     if (s == 0) {
 249         bigint_copy(dest, a);
 250         return 0;
 251     }
 252     int r;
 253     if ((r = check_size(dest, a->length_W + s))) {
 254         return r;
 255     }
 256     memmove(&dest->wordv[s], &a->wordv[0], a->length_W * sizeof(bigint_word_t));
 257     memset(&dest->wordv[0], 0, s * sizeof(bigint_word_t));
 258     dest->length_W = a->length_W + s;
 259     return 0;
 260 }
 261
 262
 263 /**
 264  * \brief a <<= 1
 265  */
 266 int bigint_shiftleft_1bit(bigint_t *a) {
 267     uint8_t c1 = 0, c2;
 268     bigint_word_t t;
 269     bigint_length_t i;
 270 #if CHECKS
 271     if (!a) {
 272         return E_PARAM;
 273     }
 274 #endif
 275
 276     for(i = 0; i < a->length_W; ++i) {
 277         t = a->wordv[i];
 278         c2 = t >> (BIGINT_WORD_SIZE - 1);
 279         t <<= 1;
 280         t |= c1;
 281         a->wordv[i] = t;
 282         c1 = c2;
 283     }
 284
 285     return 0;
 286 }
 287
 288 /**
 289  * \brief a <<= s
 290  */
 291 static
 292 int bigint_shiftleft_small(bigint_t *dest, const bigint_t *a, uint_fast8_t s) {
 293     bigint_wordplus_t t = 0;
 294     bigint_length_t i, l;
 295 #if CHECKS
 296     if (!a) {
 297         return E_PARAM;
 298     }
 299 #endif
 300     l = bigint_get_first_set_bit(a) + BIGINT_WORD_SIZE;
 301 #if CHECKS
 302     int r;
 303     if ((r = check_size(dest, (l + s ) / BIGINT_WORD_SIZE))) {
 304         return r;
 305     }
 306
 307 #endif
 308     dest->length_W = (l + s ) / BIGINT_WORD_SIZE;
 309     l /= BIGINT_WORD_SIZE;
 310     for(i = 0; i < l; ++i) {
 311         t |= (bigint_wordplus_t)a->wordv[i] << s;
 312         dest->wordv[i] = t;
 313         t >>= BIGINT_WORD_SIZE;
 314     }
 315     if (t) {
 316         dest->wordv[i] = t;
 317     }
 318     return 0;
 319 }
 320
 321 /**
 322  * \brief dest = a << s
 323  */
 324 int bigint_shiftleft(bigint_t *dest, const bigint_t *a, bigint_length_t s) {
 325     int r;
 326     if((r = bigint_shiftleft_small(dest, a, s % BIGINT_WORD_SIZE))) {
 327         return r;
 328     }
 329     if ((r = bigint_shiftleft_words(dest, dest, s / BIGINT_WORD_SIZE))) {
 330         return r;
 331     }
 332     return 0;
 333 }
 334
 335 /**
 336  * \brief a >>= 1
 337  */
 338 int bigint_shiftright_1bit(bigint_t *a) {
 339     uint8_t c1 = 0, c2;
 340     bigint_word_t t;
 341     bigint_length_t i;
 342 #if CHECKS
 343     if (!a) {
 344         return E_PARAM;
 345     }
 346 #endif
 347
 348     for(i = a->length_W; i != 0; --i) {
 349         t = a->wordv[i - 1];
 350         c2 = t & 1;
 351         t >>= 1;
 352         t |= c1 << (BIGINT_WORD_SIZE - 1);
 353         a->wordv[i - 1] = t;
 354         c1 = c2;
 355     }
 356
 357     return 0;
 358 }
 359
 360 /**
 361  * \brief dest = a ** 2
 362  */
 363 int bigint_square(bigint_t *dest, const bigint_t *a) {
 364     bigint_word_t c1, c2;
 365     bigint_wordplus_t uv;
 366     bigint_word_t q;
 367     bigint_length_t i, j;
 368 #if CHECKS
 369     int r;
 370     if (!dest || !a) {
 371         return E_PARAM;
 372     }
 373     if ((r = check_size(dest, a->length_W * 2))) {
 374         return r;
 375     }
 376 #endif
 377     if (dest == a) {
 378         bigint_t t = {NULL, 0, 0, 0};
 379         bigint_copy(&t, a);
 380         r = bigint_square(dest, &t);
 381         bigint_free(&t);
 382         return r;
 383     }
 384     memset(dest->wordv, 0, 2 * a->length_W * sizeof(bigint_word_t));
 385     for(i = 0; i < a->length_W; ++i) {
 386         uv =   (bigint_wordplus_t)dest->wordv[2 * i]
 387              + (bigint_wordplus_t)a->wordv[i]
 388              * (bigint_wordplus_t)a->wordv[i];
 389         dest->wordv[2 * i] = uv;
 390         c1 = uv >> BIGINT_WORD_SIZE;
 391         c2 = 0;
 392         for (j = i + 1; j < a->length_W; ++j) {
 393             uv =   (bigint_wordplus_t)a->wordv[i]
 394                  * (bigint_wordplus_t)a->wordv[j];
 395             q = uv >> (2 * BIGINT_WORD_SIZE - 1);
 396             uv <<= 1;
 397             uv += dest->wordv[i + j];
 398             q += (uv < dest->wordv[i + j]); /* this might be dangerous! XXX */
 399             uv += c1;
 400             q += (uv < c1); /* this might be dangerous! XXX */
 401             dest->wordv[i + j] = uv;
 402             c1 = c2 + (uv >> BIGINT_WORD_SIZE);
 403             c2 = q + (c1 < c2); /* this might be dangerous! XXX */
 404         }
 405         dest->wordv[i + a->length_W] += c1;
 406         dest->wordv[i + a->length_W + 1] = c2 + (dest->wordv[i + a->length_W] < c1);  /* this might be dangerous! XXX */
 407     }
 408     dest->length_W = a->length_W * 2;
 409     return 0;
 410 }
 411
 412
 413 /**
 414  * \brief dest = |a * b|
 415  * unsigned multiply a bigint (a) by a word (b)
 416  */
 417 int bigint_mul_word(bigint_t *dest, const bigint_t *a, const bigint_word_t b) {
 418     bigint_wordplus_t c1 = 0;
 419     bigint_length_t i;
 420 #if CHECKS
 421     int r;
 422     if (!dest || !a || !b) {
 423         return E_PARAM;
 424     }
 425     if ((r = check_size(dest, a->length_W + 1))) {
 426         return r;
 427     }
 428 #endif
 429     for(i = 0; i < a->length_W; ++i) {
 430         c1 += (bigint_wordplus_t)b * (bigint_wordplus_t)a->wordv[i];
 431         dest->wordv[i] = c1;
 432         c1 >>= BIGINT_WORD_SIZE;
 433     }
 434     dest->wordv[i++] = c1;
 435     dest->length_W = i;
 436     BIGINT_SET_POS(dest);
 437     return 0;
 438 }
 439
 440 /**
 441  * \brief dest = a * b
 442  */
 443 int bigint_mul_schoolbook(bigint_t *dest, const bigint_t *a, const bigint_t *b) {
 444     bigint_wordplus_t v;
 445     bigint_word_t c;
 446     bigint_length_t i, j;
 447     int r;
 448 #if CHECKS
 449     if (!dest || !a || !b) {
 450         return E_PARAM;
 451     }
 452     if ((r = check_size(dest, a->length_W + b->length_W))) {
 453         return r;
 454     }
 455 #endif
 456     if (dest == a) {
 457         bigint_t t = {NULL, 0, 0, 0};
 458         bigint_copy(&t, a);
 459         r = bigint_mul_schoolbook(dest, &t, b);
 460         bigint_free(&t);
 461         return r;
 462     }
 463     if (dest == b) {
 464         bigint_t t = {NULL, 0, 0, 0};
 465         bigint_copy(&t, b);
 466         r = bigint_mul_schoolbook(dest, a, &t);
 467         bigint_free(&t);
 468         return r;
 469     }
 470     memset(dest->wordv, 0, (a->length_W + b->length_W) * sizeof(bigint_word_t));
 471     for(i = 0; i < b->length_W; ++i) {
 472         c = 0;
 473         for(j = 0; j < a->length_W; ++j) {
 474             v =   (bigint_wordplus_t)a->wordv[j]
 475                 * (bigint_wordplus_t)b->wordv[i];
 476             v += (bigint_wordplus_t)dest->wordv[i + j];
 477             v += (bigint_wordplus_t)c;
 478             dest->wordv[i + j] = v;
 479             c = v >> BIGINT_WORD_SIZE;
 480         }
 481         dest->wordv[i + j] = c;
 482     }
 483
 484     dest->length_W = a->length_W + b->length_W;
 485     dest->info &= ~BIGINT_SIGN_MASK;
 486     dest->info |= (a->info ^ b->info) & BIGINT_SIGN_MASK;
 487     return 0;
 488 }
 489
 490 /*
 491  * UNSAFE!!!
 492  */
 493 int8_t bigint_cmp_u(const bigint_t *a, const bigint_t *b) {
 494     int8_t r = 1;
 495     bigint_length_t i;
 496     if (a->length_W == 0 && b->length_W == 0) {
 497         return 0;
 498     }
 499     if (b->length_W == 0) {
 500         return 1;
 501     }
 502     if (a->length_W == 0) {
 503         return -1;
 504     }
 505     if (a->length_W < b->length_W) {
 506         const bigint_t *t;
 507         t = a;
 508         a = b;
 509         b = t;
 510         r = -1;
 511     }
 512     for (i = a->length_W - 1; i > b->length_W - 1; --i) {
 513         if (a->wordv[i]) {
 514             return r;
 515         }
 516     }
 517     ++i;
 518     do {
 519         --i;
 520         if (a->wordv[i] != b->wordv[i]) {
 521             return a->wordv[i] > b->wordv[i] ? r : -r;
 522         }
 523     } while (i);
 524     return 0;
 525 }
 526
 527 int8_t bigint_cmp_s(const bigint_t *a, const bigint_t *b) {
 528     if ( bigint_get_first_set_bit(a) == (bigint_length_t)-1 &&
 529          bigint_get_first_set_bit(b) == (bigint_length_t)-1 ) {
 530         return 0;
 531     }
 532     if ((a->info & BIGINT_SIGN_MASK) ^ (b->info & BIGINT_SIGN_MASK)) {
 533         return (a->info & BIGINT_SIGN_MASK) ? -1 : 1;
 534     } else {
 535         int r;
 536         r = bigint_cmp_u(a, b);
 537         return (a->info & BIGINT_SIGN_MASK) ? -r : r;
 538     }
 539     return 0;
 540 }
 541
 542 /*
 543  * UNSAFE!!!
 544  * a <=> b * (B ** s)
 545  */
 546 int8_t bigint_cmp_scale(const bigint_t *a, const bigint_t *b, bigint_length_t s) {
 547     bigint_length_t i;
 548     if (s == 0) {
 549         return bigint_cmp_u(a, b);
 550     }
 551     if (a->length_W == 0 && b->length_W == 0) {
 552         return 0;
 553     }
 554     if (b->length_W == 0) {
 555         return 1;
 556     }
 557     if (a->length_W == 0) {
 558         return -1;
 559     }
 560     if (s > a->length_W) {
 561         return -1;
 562     }
 563     for (i = a->length_W - 1; i > b->length_W + s - 1; --i) {
 564         if (a->wordv[i]) {
 565             return 1;
 566         }
 567     }
 568     for (; i > s - 1; --i) {
 569         if (a->wordv[i] != b->wordv[i + s]){
 570             return a->wordv[i] > b->wordv[i + s] ? 1 : -1;
 571         }
 572     }
 573         ++i;
 574     do {
 575         --i;
 576         if (a->wordv[i]) {
 577             return 1;
 578         }
 579     } while (i);
 580     return 0;
 581 }
 582
 583 /**
 584  * XXXXXXXXXXXX
 585  * \brief dest = |a| - |b| * B ** s
 586  * Subtracts b from a and stores the result in dest.
 587  * Signs are ignored
 588  */
 589 int bigint_sub_scale(bigint_t *dest, const bigint_t *a, const bigint_t *b, bigint_length_t s) {
 590     bigint_length_t i, j;
 591     bigint_word_t t, c1, c2, c = 0;
 592     int r;
 593 #if CHECKS
 594     if (!dest || !a || !b) {
 595         return E_PARAM;
 596     }
 597 #endif
 598 #if CHECKS
 599     if ((r = check_size(dest, MAX(a->length_W, b->length_W)))) {
 600         return r;
 601     }
 602
 603 #endif
 604     j = MIN(a->length_W, b->length_W);
 605     for (i = 0; i < j; ++i) {
 606         t = a->wordv[i + s] - b->wordv[i];
 607         c1 = t > a->wordv[i + s];
 608         dest->wordv[i + s] = t - c;
 609         c2 = dest->wordv[i + s] > t;
 610         c = c1 | c2;
 611     }
 612     for (; i < a->length_W; ++i) {
 613         t = a->wordv[i + s];
 614         dest->wordv[i + s] = t - c;
 615         c = dest->wordv[i + s] > t;
 616     }
 617     dest->length_W = a->length_W;
 618     return 0;
 619 }
 620
 621 /*
 622  * UNSAFE!!!
 623  */
 624 int bigint_adjust_length(bigint_t *a) {
 625 #if CHECKS
 626     if (!a) {
 627         return E_PARAM;
 628     }
 629 #endif
 630     while (a->length_W) {
 631         if (a->wordv[a->length_W - 1]) {
 632             return 0;
 633         }
 634         a->length_W--;
 635     }
 636     return 0;
 637 }
 638
 639 bigint_length_t bigint_get_first_set_bit(const bigint_t *a) {
 640     bigint_length_t r;
 641     bigint_word_t t;
 642     if (!a) {
 643         return (bigint_length_t)-1;
 644     }
 645     if (a->length_W == 0) {
 646         return (bigint_length_t)-1;
 647     }
 648     r = a->length_W - 1;
 649     while (r && a->wordv[r] == 0) {
 650         --r;
 651     }
 652     t = a->wordv[r];
 653     if (!t) {
 654         return (bigint_length_t)-1;
 655     }
 656     r *= BIGINT_WORD_SIZE;
 657     t >>= 1;
 658     while (t) {
 659         t >>= 1;
 660         ++r;
 661     }
 662     return r;
 663 }
 664
 665 /*
 666  * UNSAFE!!!
 667  * a = q * b + r
 668  */
 669 int bigint_divide(bigint_t *q, bigint_t *r, const bigint_t *a, const bigint_t *b) {
 670     bigint_t a_ = {NULL, 0, 0, 0}, x_ = {NULL, 0, 0, 0};
 671     bigint_length_t i, la, lb;
 672     int ret;
 673 #if CHECKS
 674     if (!a || !b || (!q && !r)) {
 675         return E_PARAM;
 676     }
 677 #endif
 678     la = bigint_get_first_set_bit(a);
 679     lb = bigint_get_first_set_bit(b);
 680     if (lb == (bigint_length_t)-1) {
 681         return E_VALUE;
 682     }
 683     if (la == (bigint_length_t)-1) {
 684         if (q) {
 685             q->length_W = 0;
 686         }
 687         if (r) {
 688             r->length_W = 0;
 689         }
 690         return 0;
 691     }
 692     if (la < lb) {
 693         if (q) {
 694             q->length_W = 0;
 695         }
 696         if (r) {
 697             bigint_copy(r, a);
 698         }
 699         return 0;
 700     }
 701     i = la - lb;
 702     if (q) {
 703         if ((ret = check_size(q, (i + BIGINT_WORD_SIZE - 1) / BIGINT_WORD_SIZE))) {
 704             return ret;
 705         }
 706         memset(q->wordv, 0, q->allocated_W * sizeof(bigint_word_t));
 707     }
 708     if (r) {
 709         if ((ret = check_size(q, (lb + BIGINT_WORD_SIZE - 1) / BIGINT_WORD_SIZE))) {
 710             return ret;
 711         }
 712     }
 713     if ((ret = bigint_copy(&a_, a))) {
 714         return ret;
 715     }
 716     if ((ret = bigint_shiftleft(&x_, b, i))) {
 717         bigint_free(&a_);
 718         return ret;
 719     }
 720 #if 0
 721     printf("DBG: x' = ");
 722     bigint_print_hex(&x_);
 723     printf("; b = ");
 724     bigint_print_hex(b);
 725     printf("; la = %d; lb = %d; i = %d\n", la, lb, i);
 726 #endif
 727     do {
 728         if (bigint_cmp_u(&a_, &x_) >= 0) {
 729             bigint_sub_u(&a_, &a_, &x_);
 730             if (q) {
 731                 q->wordv[i / BIGINT_WORD_SIZE] |= 1 << (i % BIGINT_WORD_SIZE);
 732             }
 733         }
 734         bigint_shiftright_1bit(&x_);
 735     } while (i--);
 736     if (r) {
 737         bigint_copy(r, &a_);
 738     }
 739     bigint_free(&x_);
 740     bigint_free(&a_);
 741     return 0;
 742 }
 743
 744 /**
 745  * UNSAFE!!!
 746  */
 747 int bigint_sub_s(bigint_t *dest, const bigint_t *a, const bigint_t *b) {
 748     int8_t s;
 749     int r;
 750     if (a->length_W == 0) {
 751         bigint_copy(dest, b);
 752         dest->info ^= BIGINT_SIGN_MASK;
 753         return 0;
 754     }
 755     if ((a->info ^ b->info) & BIGINT_SIGN_MASK) {
 756         /* different signs */
 757         if ((r = bigint_add_auto_u(dest, a, b))) {
 758             return r;
 759         }
 760         dest->info = a->info;
 761     } else {
 762         s = bigint_cmp_u(a, b);
 763         if (s >= 0) {
 764             if (( r = bigint_sub_u(dest, a, b))) {
 765                 return r;
 766             }
 767             dest->info = a->info;
 768         } else {
 769             if (( r = bigint_sub_u(dest, b, a))) {
 770                 return r;
 771             }
 772             dest->info = (~a->info) & BIGINT_SIGN_MASK;
 773         }
 774     }
 775     return 0;
 776 }
 777
 778 /**
 779  * UNSAFE!!!
 780  */
 781 int bigint_add_s(bigint_t *dest, const bigint_t *a, const bigint_t *b) {
 782     int8_t s;
 783     int r;
 784     if ((a->info ^ b->info) & BIGINT_SIGN_MASK) {
 785         /* different signs */
 786         s = bigint_cmp_u(a, b);
 787         if (s >= 0) {
 788             if (( r = bigint_sub_u(dest, a, b))) {
 789                 return r;
 790             }
 791             dest->info = a->info;
 792         } else {
 793             if (( r = bigint_sub_u(dest, b, a))) {
 794                 return r;
 795             }
 796             dest->info = (~a->info) & BIGINT_SIGN_MASK;
 797         }
 798     } else {
 799         if ((r = bigint_add_auto_u(dest, a, b))) {
 800             return r;
 801         }
 802         dest->info = a->info;
 803     }
 804     return 0;
 805 }
 806
 807 int8_t bigint_is_even(const bigint_t *a) {
 808     if (!a) {
 809         return E_VALUE;
 810     }
 811     if (a->length_W == 0 || (a->wordv[0] & 1) == 0) {
 812         return 1;
 813     }
 814     return 0;
 815 }
 816
 817 int8_t bigint_is_odd(const bigint_t *a) {
 818     if (!a) {
 819         return E_VALUE;
 820     }
 821     if (a->length_W > 0 && (a->wordv[0] & 1) == 1) {
 822         return 1;
 823     }
 824     return 0;
 825 }
 826
 827 /**
 828  * UNSAFE!!!
 829  * (a, b) -> (x, y, v)
 830  * v = gcd(a,b)
 831  * x * a + y * b = v
 832  */
 833 int bigint_gcdext(bigint_t *gcd, bigint_t *x, bigint_t *y, const bigint_t *a, const bigint_t *b) {
 834     bigint_length_t g = 0;
 835     bigint_t u = {NULL, 0, 0, 0}, v = {NULL, 0, 0, 0};
 836     bigint_t x_ = {NULL, 0, 0, 0}, y_ = {NULL, 0, 0, 0};
 837     bigint_t A = {NULL, 0, 0, 0}, B = {NULL, 0, 0, 0};
 838     bigint_t C = {NULL, 0, 0, 0}, D = {NULL, 0, 0, 0};
 839     int r;
 840     if ((r = bigint_copy(&x_, a))) {
 841         return r;
 842     }
 843     if ((r = bigint_copy(&y_, b))) {
 844         bigint_free(&x_);
 845         return r;
 846     }
 847     bigint_adjust_length(&x_);
 848     bigint_adjust_length(&y_);
 849     if (x_.length_W == 0 || y_.length_W == 0) {
 850         bigint_free(&y_);
 851         bigint_free(&x_);
 852         return E_VALUE;
 853     }
 854     while (((x_.wordv[0] | y_.wordv[0]) & 1) == 0) {
 855         ++g;
 856         bigint_shiftright_1bit(&x_);
 857         bigint_shiftright_1bit(&y_);
 858     }
 859     if ((r = check_size(&A, 1 + (bigint_get_first_set_bit(&y_) + BIGINT_WORD_SIZE) / BIGINT_WORD_SIZE))) {
 860         bigint_free(&y_);
 861         bigint_free(&x_);
 862         return r;
 863     }
 864     if ((r = check_size(&C, 1 + (bigint_get_first_set_bit(&y_) + BIGINT_WORD_SIZE) / BIGINT_WORD_SIZE))) {
 865         bigint_free(&A);
 866         bigint_free(&y_);
 867         bigint_free(&x_);
 868         return r;
 869     }
 870     if ((r = check_size(&B, 1 + (bigint_get_first_set_bit(&x_) + BIGINT_WORD_SIZE) / BIGINT_WORD_SIZE))) {
 871         bigint_free(&C);
 872         bigint_free(&A);
 873         bigint_free(&y_);
 874         bigint_free(&x_);
 875         return r;
 876     }
 877     if ((r = check_size(&D, 1 + (bigint_get_first_set_bit(&x_) + BIGINT_WORD_SIZE) / BIGINT_WORD_SIZE))) {
 878         bigint_free(&B);
 879         bigint_free(&C);
 880         bigint_free(&A);
 881         bigint_free(&y_);
 882         bigint_free(&x_);
 883         return r;
 884     }
 885     if ((r = bigint_copy(&u, &x_))) {
 886         bigint_free(&D);
 887         bigint_free(&B);
 888         bigint_free(&C);
 889         bigint_free(&A);
 890         bigint_free(&y_);
 891         bigint_free(&x_);
 892         return r;
 893     }
 894     if ((r = bigint_copy(&v, &y_))) {
 895         bigint_free(&u);
 896         bigint_free(&D);
 897         bigint_free(&B);
 898         bigint_free(&C);
 899         bigint_free(&A);
 900         bigint_free(&y_);
 901         bigint_free(&x_);
 902         return r;
 903     }
 904     A.wordv[0] = D.wordv[0] = 1;
 905     A.length_W = D.length_W = 1;
 906     C.length_W = B.length_W = 0;
 907
 908     do {
 909         while ((u.wordv[0] & 1) == 0) {
 910             bigint_shiftright_1bit(&u);
 911             if (bigint_is_odd(&A) || bigint_is_odd(&B)) {
 912                 bigint_add_s(&A, &A, &y_);
 913                 bigint_sub_s(&B, &B, &x_);
 914             }
 915             bigint_shiftright_1bit(&A);
 916             bigint_shiftright_1bit(&B);
 917         }
 918         while ((v.wordv[0] & 1) == 0) {
 919             bigint_shiftright_1bit(&v);
 920             if (bigint_is_odd(&C) || bigint_is_odd(&D)) {
 921                 bigint_add_s(&C, &C, &y_);
 922                 bigint_sub_s(&D, &D, &x_);
 923             }
 924             bigint_shiftright_1bit(&C);
 925             bigint_shiftright_1bit(&D);
 926         }
 927         if (bigint_cmp_s(&u, &v) >= 0) {
 928             bigint_sub_s(&u, &u, &v);
 929             bigint_sub_s(&A, &A, &C);
 930             bigint_sub_s(&B, &B, &D);
 931         } else {
 932             bigint_sub_s(&v, &v, &u);
 933             bigint_sub_s(&C, &C, &A);
 934             bigint_sub_s(&D, &D, &B);
 935         }
 936         bigint_adjust_length(&u);
 937     } while (u.length_W != 0);
 938     if (gcd) {
 939         bigint_shiftleft(gcd, &v, g);
 940     }
 941     if (x) {
 942         bigint_copy(x, &C);
 943     }
 944     if (y) {
 945         bigint_copy(y, &D);
 946     }
 947     bigint_free(&v);
 948     bigint_free(&u);
 949     bigint_free(&D);
 950     bigint_free(&B);
 951     bigint_free(&C);
 952     bigint_free(&A);
 953     bigint_free(&y_);
 954     bigint_free(&x_);
 955     return 0;
 956 }