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) / BIGINT_WORD_SIZE))) {
 704             return ret;
 705         }
 706         q->length_W = (i + BIGINT_WORD_SIZE) / BIGINT_WORD_SIZE;
 707         memset(q->wordv, 0, q->allocated_W * sizeof(bigint_word_t));
 708     }
 709     if (r) {
 710         if ((ret = check_size(q, (lb + BIGINT_WORD_SIZE - 1) / BIGINT_WORD_SIZE))) {
 711             return ret;
 712         }
 713     }
 714     if ((ret = bigint_copy(&a_, a))) {
 715         return ret;
 716     }
 717     if ((ret = bigint_shiftleft(&x_, b, i))) {
 718         bigint_free(&a_);
 719         return ret;
 720     }
 721 #if 0
 722     printf("DBG: x' = ");
 723     bigint_print_hex(&x_);
 724     printf("; b = ");
 725     bigint_print_hex(b);
 726     printf("; la = %d; lb = %d; i = %d\n", la, lb, i);
 727 #endif
 728     do {
 729         if (bigint_cmp_u(&a_, &x_) >= 0) {
 730             bigint_sub_u(&a_, &a_, &x_);
 731             if (q) {
 732                 q->wordv[i / BIGINT_WORD_SIZE] |= 1 << (i % BIGINT_WORD_SIZE);
 733             }
 734         }
 735         bigint_shiftright_1bit(&x_);
 736     } while (i--);
 737     if (r) {
 738         bigint_copy(r, &a_);
 739     }
 740     bigint_free(&x_);
 741     bigint_free(&a_);
 742     return 0;
 743 }
 744
 745 /**
 746  * UNSAFE!!!
 747  */
 748 int bigint_sub_s(bigint_t *dest, const bigint_t *a, const bigint_t *b) {
 749     int8_t s;
 750     int r;
 751     if (a->length_W == 0) {
 752         bigint_copy(dest, b);
 753         dest->info ^= BIGINT_SIGN_MASK;
 754         return 0;
 755     }
 756     if ((a->info ^ b->info) & BIGINT_SIGN_MASK) {
 757         /* different signs */
 758         if ((r = bigint_add_auto_u(dest, a, b))) {
 759             return r;
 760         }
 761         dest->info = a->info;
 762     } else {
 763         s = bigint_cmp_u(a, b);
 764         if (s >= 0) {
 765             if (( r = bigint_sub_u(dest, a, b))) {
 766                 return r;
 767             }
 768             dest->info = a->info;
 769         } else {
 770             if (( r = bigint_sub_u(dest, b, a))) {
 771                 return r;
 772             }
 773             dest->info = (~a->info) & BIGINT_SIGN_MASK;
 774         }
 775     }
 776     return 0;
 777 }
 778
 779 /**
 780  * UNSAFE!!!
 781  */
 782 int bigint_add_s(bigint_t *dest, const bigint_t *a, const bigint_t *b) {
 783     int8_t s;
 784     int r;
 785     if ((a->info ^ b->info) & BIGINT_SIGN_MASK) {
 786         /* different signs */
 787         s = bigint_cmp_u(a, b);
 788         if (s >= 0) {
 789             if (( r = bigint_sub_u(dest, a, b))) {
 790                 return r;
 791             }
 792             dest->info = a->info;
 793         } else {
 794             if (( r = bigint_sub_u(dest, b, a))) {
 795                 return r;
 796             }
 797             dest->info = (~a->info) & BIGINT_SIGN_MASK;
 798         }
 799     } else {
 800         if ((r = bigint_add_auto_u(dest, a, b))) {
 801             return r;
 802         }
 803         dest->info = a->info;
 804     }
 805     return 0;
 806 }
 807
 808 int8_t bigint_is_even(const bigint_t *a) {
 809     if (!a) {
 810         return E_VALUE;
 811     }
 812     if (a->length_W == 0 || (a->wordv[0] & 1) == 0) {
 813         return 1;
 814     }
 815     return 0;
 816 }
 817
 818 int8_t bigint_is_odd(const bigint_t *a) {
 819     if (!a) {
 820         return E_VALUE;
 821     }
 822     if (a->length_W > 0 && (a->wordv[0] & 1) == 1) {
 823         return 1;
 824     }
 825     return 0;
 826 }
 827
 828 /**
 829  * UNSAFE!!!
 830  * (a, b) -> (x, y, v)
 831  * v = gcd(a,b)
 832  * x * a + y * b = v
 833  */
 834 int bigint_gcdext(bigint_t *gcd, bigint_t *x, bigint_t *y, const bigint_t *a, const bigint_t *b) {
 835     bigint_length_t g = 0;
 836     bigint_t u = {NULL, 0, 0, 0}, v = {NULL, 0, 0, 0};
 837     bigint_t x_ = {NULL, 0, 0, 0}, y_ = {NULL, 0, 0, 0};
 838     bigint_t A = {NULL, 0, 0, 0}, B = {NULL, 0, 0, 0};
 839     bigint_t C = {NULL, 0, 0, 0}, D = {NULL, 0, 0, 0};
 840     int r;
 841     if ((r = bigint_copy(&x_, a))) {
 842         return r;
 843     }
 844     if ((r = bigint_copy(&y_, b))) {
 845         bigint_free(&x_);
 846         return r;
 847     }
 848     bigint_adjust_length(&x_);
 849     bigint_adjust_length(&y_);
 850     if (x_.length_W == 0 || y_.length_W == 0) {
 851         bigint_free(&y_);
 852         bigint_free(&x_);
 853         return E_VALUE;
 854     }
 855     while (((x_.wordv[0] | y_.wordv[0]) & 1) == 0) {
 856         ++g;
 857         bigint_shiftright_1bit(&x_);
 858         bigint_shiftright_1bit(&y_);
 859     }
 860     if ((r = check_size(&A, 1 + (bigint_get_first_set_bit(&y_) + BIGINT_WORD_SIZE) / BIGINT_WORD_SIZE))) {
 861         bigint_free(&y_);
 862         bigint_free(&x_);
 863         return r;
 864     }
 865     if ((r = check_size(&C, 1 + (bigint_get_first_set_bit(&y_) + BIGINT_WORD_SIZE) / BIGINT_WORD_SIZE))) {
 866         bigint_free(&A);
 867         bigint_free(&y_);
 868         bigint_free(&x_);
 869         return r;
 870     }
 871     if ((r = check_size(&B, 1 + (bigint_get_first_set_bit(&x_) + BIGINT_WORD_SIZE) / BIGINT_WORD_SIZE))) {
 872         bigint_free(&C);
 873         bigint_free(&A);
 874         bigint_free(&y_);
 875         bigint_free(&x_);
 876         return r;
 877     }
 878     if ((r = check_size(&D, 1 + (bigint_get_first_set_bit(&x_) + BIGINT_WORD_SIZE) / BIGINT_WORD_SIZE))) {
 879         bigint_free(&B);
 880         bigint_free(&C);
 881         bigint_free(&A);
 882         bigint_free(&y_);
 883         bigint_free(&x_);
 884         return r;
 885     }
 886     if ((r = bigint_copy(&u, &x_))) {
 887         bigint_free(&D);
 888         bigint_free(&B);
 889         bigint_free(&C);
 890         bigint_free(&A);
 891         bigint_free(&y_);
 892         bigint_free(&x_);
 893         return r;
 894     }
 895     if ((r = bigint_copy(&v, &y_))) {
 896         bigint_free(&u);
 897         bigint_free(&D);
 898         bigint_free(&B);
 899         bigint_free(&C);
 900         bigint_free(&A);
 901         bigint_free(&y_);
 902         bigint_free(&x_);
 903         return r;
 904     }
 905     A.wordv[0] = D.wordv[0] = 1;
 906     A.length_W = D.length_W = 1;
 907     C.length_W = B.length_W = 0;
 908
 909     do {
 910         while ((u.wordv[0] & 1) == 0) {
 911             bigint_shiftright_1bit(&u);
 912             if (bigint_is_odd(&A) || bigint_is_odd(&B)) {
 913                 bigint_add_s(&A, &A, &y_);
 914                 bigint_sub_s(&B, &B, &x_);
 915             }
 916             bigint_shiftright_1bit(&A);
 917             bigint_shiftright_1bit(&B);
 918         }
 919         while ((v.wordv[0] & 1) == 0) {
 920             bigint_shiftright_1bit(&v);
 921             if (bigint_is_odd(&C) || bigint_is_odd(&D)) {
 922                 bigint_add_s(&C, &C, &y_);
 923                 bigint_sub_s(&D, &D, &x_);
 924             }
 925             bigint_shiftright_1bit(&C);
 926             bigint_shiftright_1bit(&D);
 927         }
 928         if (bigint_cmp_s(&u, &v) >= 0) {
 929             bigint_sub_s(&u, &u, &v);
 930             bigint_sub_s(&A, &A, &C);
 931             bigint_sub_s(&B, &B, &D);
 932         } else {
 933             bigint_sub_s(&v, &v, &u);
 934             bigint_sub_s(&C, &C, &A);
 935             bigint_sub_s(&D, &D, &B);
 936         }
 937         bigint_adjust_length(&u);
 938     } while (u.length_W != 0);
 939     if (gcd) {
 940         bigint_shiftleft(gcd, &v, g);
 941     }
 942     if (x) {
 943         bigint_copy(x, &C);
 944     }
 945     if (y) {
 946         bigint_copy(y, &D);
 947     }
 948     bigint_free(&v);
 949     bigint_free(&u);
 950     bigint_free(&D);
 951     bigint_free(&B);
 952     bigint_free(&C);
 953     bigint_free(&A);
 954     bigint_free(&y_);
 955     bigint_free(&x_);
 956     return 0;
 957 }