contrib/bc/src/num.c

   1 /*
   2  * *****************************************************************************
   3  *
   4  * SPDX-License-Identifier: BSD-2-Clause
   5  *
   6  * Copyright (c) 2018-2020 Gavin D. Howard and contributors.
   7  *
   8  * Redistribution and use in source and binary forms, with or without
   9  * modification, are permitted provided that the following conditions are met:
  10  *
  11  * * Redistributions of source code must retain the above copyright notice, this
  12  *   list of conditions and the following disclaimer.
  13  *
  14  * * Redistributions in binary form must reproduce the above copyright notice,
  15  *   this list of conditions and the following disclaimer in the documentation
  16  *   and/or other materials provided with the distribution.
  17  *
  18  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
  19  * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  20  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  21  * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
  22  * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
  23  * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
  24  * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
  25  * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
  26  * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
  27  * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
  28  * POSSIBILITY OF SUCH DAMAGE.
  29  *
  30  * *****************************************************************************
  31  *
  32  * Code for the number type.
  33  *
  34  */
  35
  36 #include <assert.h>
  37 #include <ctype.h>
  38 #include <stdbool.h>
  39 #include <stdlib.h>
  40 #include <string.h>
  41 #include <setjmp.h>
  42 #include <limits.h>
  43
  44 #include <status.h>
  45 #include <num.h>
  46 #include <rand.h>
  47 #include <vm.h>
  48
  49 static void bc_num_m(BcNum *a, BcNum *b, BcNum *restrict c, size_t scale);
  50
  51 static inline ssize_t bc_num_neg(size_t n, bool neg) {
  52         return (((ssize_t) n) ^ -((ssize_t) neg)) + neg;
  53 }
  54
  55 ssize_t bc_num_cmpZero(const BcNum *n) {
  56         return bc_num_neg((n)->len != 0, (n)->neg);
  57 }
  58
  59 static inline size_t bc_num_int(const BcNum *n) {
  60         return n->len ? n->len - n->rdx : 0;
  61 }
  62
  63 static void bc_num_expand(BcNum *restrict n, size_t req) {
  64
  65         assert(n != NULL);
  66
  67         req = req >= BC_NUM_DEF_SIZE ? req : BC_NUM_DEF_SIZE;
  68
  69         if (req > n->cap) {
  70
  71                 BC_SIG_LOCK;
  72
  73                 n->num = bc_vm_realloc(n->num, BC_NUM_SIZE(req));
  74                 n->cap = req;
  75
  76                 BC_SIG_UNLOCK;
  77         }
  78 }
  79
  80 static void bc_num_setToZero(BcNum *restrict n, size_t scale) {
  81         assert(n != NULL);
  82         n->scale = scale;
  83         n->len = n->rdx = 0;
  84         n->neg = false;
  85 }
  86
  87 static inline void bc_num_zero(BcNum *restrict n) {
  88         bc_num_setToZero(n, 0);
  89 }
  90
  91 void bc_num_one(BcNum *restrict n) {
  92         bc_num_zero(n);
  93         n->len = 1;
  94         n->num[0] = 1;
  95 }
  96
  97 static void bc_num_clean(BcNum *restrict n) {
  98
  99         while (BC_NUM_NONZERO(n) && !n->num[n->len - 1]) n->len -= 1;
 100
 101         if (BC_NUM_ZERO(n)) {
 102                 n->neg = false;
 103                 n->rdx = 0;
 104         }
 105         else if (n->len < n->rdx) n->len = n->rdx;
 106 }
 107
 108 static size_t bc_num_log10(size_t i) {
 109         size_t len;
 110         for (len = 1; i; i /= BC_BASE, ++len);
 111         assert(len - 1 <= BC_BASE_DIGS + 1);
 112         return len - 1;
 113 }
 114
 115 static inline size_t bc_num_zeroDigits(const BcDig *n) {
 116         assert(*n >= 0);
 117         assert(((size_t) *n) < BC_BASE_POW);
 118         return BC_BASE_DIGS - bc_num_log10((size_t) *n);
 119 }
 120
 121 static size_t bc_num_intDigits(const BcNum *n) {
 122         size_t digits = bc_num_int(n) * BC_BASE_DIGS;
 123         if (digits > 0) digits -= bc_num_zeroDigits(n->num + n->len - 1);
 124         return digits;
 125 }
 126
 127 static size_t bc_num_nonzeroLen(const BcNum *restrict n) {
 128         size_t i, len = n->len;
 129         assert(len == n->rdx);
 130         for (i = len - 1; i < len && !n->num[i]; --i);
 131         assert(i + 1 > 0);
 132         return i + 1;
 133 }
 134
 135 static BcDig bc_num_addDigits(BcDig a, BcDig b, bool *carry) {
 136
 137         assert(((BcBigDig) BC_BASE_POW) * 2 == ((BcDig) BC_BASE_POW) * 2);
 138         assert(a < BC_BASE_POW);
 139         assert(b < BC_BASE_POW);
 140
 141         a += b + *carry;
 142         *carry = (a >= BC_BASE_POW);
 143         if (*carry) a -= BC_BASE_POW;
 144
 145         assert(a >= 0);
 146         assert(a < BC_BASE_POW);
 147
 148         return a;
 149 }
 150
 151 static BcDig bc_num_subDigits(BcDig a, BcDig b, bool *carry) {
 152
 153         assert(a < BC_BASE_POW);
 154         assert(b < BC_BASE_POW);
 155
 156         b += *carry;
 157         *carry = (a < b);
 158         if (*carry) a += BC_BASE_POW;
 159
 160         assert(a - b >= 0);
 161         assert(a - b < BC_BASE_POW);
 162
 163         return a - b;
 164 }
 165
 166 static void bc_num_addArrays(BcDig *restrict a, const BcDig *restrict b,
 167                              size_t len)
 168 {
 169         size_t i;
 170         bool carry = false;
 171
 172         for (i = 0; i < len; ++i) a[i] = bc_num_addDigits(a[i], b[i], &carry);
 173
 174         for (; carry; ++i) a[i] = bc_num_addDigits(a[i], 0, &carry);
 175 }
 176
 177 static void bc_num_subArrays(BcDig *restrict a, const BcDig *restrict b,
 178                              size_t len)
 179 {
 180         size_t i;
 181         bool carry = false;
 182
 183         for (i = 0; i < len; ++i) a[i] = bc_num_subDigits(a[i], b[i], &carry);
 184
 185         for (; carry; ++i) a[i] = bc_num_subDigits(a[i], 0, &carry);
 186 }
 187
 188 static void bc_num_mulArray(const BcNum *restrict a, BcBigDig b,
 189                             BcNum *restrict c)
 190 {
 191         size_t i;
 192         BcBigDig carry = 0;
 193
 194         assert(b <= BC_BASE_POW);
 195
 196         if (a->len + 1 > c->cap) bc_num_expand(c, a->len + 1);
 197
 198         memset(c->num, 0, BC_NUM_SIZE(c->cap));
 199
 200         for (i = 0; i < a->len; ++i) {
 201                 BcBigDig in = ((BcBigDig) a->num[i]) * b + carry;
 202                 c->num[i] = in % BC_BASE_POW;
 203                 carry = in / BC_BASE_POW;
 204         }
 205
 206         assert(carry < BC_BASE_POW);
 207         c->num[i] = (BcDig) carry;
 208         c->len = a->len;
 209         c->len += (carry != 0);
 210
 211         bc_num_clean(c);
 212
 213         assert(!c->neg || BC_NUM_NONZERO(c));
 214         assert(c->rdx <= c->len || !c->len);
 215         assert(!c->len || c->num[c->len - 1] || c->rdx == c->len);
 216 }
 217
 218 static void bc_num_divArray(const BcNum *restrict a, BcBigDig b,
 219                             BcNum *restrict c, BcBigDig *rem)
 220 {
 221         size_t i;
 222         BcBigDig carry = 0;
 223
 224         assert(c->cap >= a->len);
 225
 226         for (i = a->len - 1; i < a->len; --i) {
 227                 BcBigDig in = ((BcBigDig) a->num[i]) + carry * BC_BASE_POW;
 228                 assert(in / b < BC_BASE_POW);
 229                 c->num[i] = (BcDig) (in / b);
 230                 carry = in % b;
 231         }
 232
 233         c->len = a->len;
 234         bc_num_clean(c);
 235         *rem = carry;
 236
 237         assert(!c->neg || BC_NUM_NONZERO(c));
 238         assert(c->rdx <= c->len || !c->len);
 239         assert(!c->len || c->num[c->len - 1] || c->rdx == c->len);
 240 }
 241
 242 static ssize_t bc_num_compare(const BcDig *restrict a, const BcDig *restrict b,
 243                               size_t len)
 244 {
 245         size_t i;
 246         BcDig c = 0;
 247         for (i = len - 1; i < len && !(c = a[i] - b[i]); --i);
 248         return bc_num_neg(i + 1, c < 0);
 249 }
 250
 251 ssize_t bc_num_cmp(const BcNum *a, const BcNum *b) {
 252
 253         size_t i, min, a_int, b_int, diff;
 254         BcDig *max_num, *min_num;
 255         bool a_max, neg = false;
 256         ssize_t cmp;
 257
 258         assert(a != NULL && b != NULL);
 259
 260         if (a == b) return 0;
 261         if (BC_NUM_ZERO(a)) return bc_num_neg(b->len != 0, !b->neg);
 262         if (BC_NUM_ZERO(b)) return bc_num_cmpZero(a);
 263         if (a->neg) {
 264                 if (b->neg) neg = true;
 265                 else return -1;
 266         }
 267         else if (b->neg) return 1;
 268
 269         a_int = bc_num_int(a);
 270         b_int = bc_num_int(b);
 271         a_int -= b_int;
 272
 273         if (a_int) return neg ? -((ssize_t) a_int) : (ssize_t) a_int;
 274
 275         a_max = (a->rdx > b->rdx);
 276
 277         if (a_max) {
 278                 min = b->rdx;
 279                 diff = a->rdx - b->rdx;
 280                 max_num = a->num + diff;
 281                 min_num = b->num;
 282         }
 283         else {
 284                 min = a->rdx;
 285                 diff = b->rdx - a->rdx;
 286                 max_num = b->num + diff;
 287                 min_num = a->num;
 288         }
 289
 290         cmp = bc_num_compare(max_num, min_num, b_int + min);
 291
 292         if (cmp) return bc_num_neg((size_t) cmp, !a_max == !neg);
 293
 294         for (max_num -= diff, i = diff - 1; i < diff; --i) {
 295                 if (max_num[i]) return bc_num_neg(1, !a_max == !neg);
 296         }
 297
 298         return 0;
 299 }
 300
 301 void bc_num_truncate(BcNum *restrict n, size_t places) {
 302
 303         size_t places_rdx;
 304
 305         if (!places) return;
 306
 307         places_rdx = n->rdx ? n->rdx - BC_NUM_RDX(n->scale - places) : 0;
 308         assert(places <= n->scale && (BC_NUM_ZERO(n) || places_rdx <= n->len));
 309
 310         n->scale -= places;
 311         n->rdx -= places_rdx;
 312
 313         if (BC_NUM_NONZERO(n)) {
 314
 315                 size_t pow;
 316
 317                 pow = n->scale % BC_BASE_DIGS;
 318                 pow = pow ? BC_BASE_DIGS - pow : 0;
 319                 pow = bc_num_pow10[pow];
 320
 321                 n->len -= places_rdx;
 322                 memmove(n->num, n->num + places_rdx, BC_NUM_SIZE(n->len));
 323
 324                 // Clear the lower part of the last digit.
 325                 if (BC_NUM_NONZERO(n)) n->num[0] -= n->num[0] % (BcDig) pow;
 326
 327                 bc_num_clean(n);
 328         }
 329 }
 330
 331 static void bc_num_extend(BcNum *restrict n, size_t places) {
 332
 333         size_t places_rdx;
 334
 335         if (!places) return;
 336         if (BC_NUM_ZERO(n)) {
 337                 n->scale += places;
 338                 return;
 339         }
 340
 341         places_rdx = BC_NUM_RDX(places + n->scale) - n->rdx;
 342
 343         if (places_rdx) {
 344                 bc_num_expand(n, bc_vm_growSize(n->len, places_rdx));
 345                 memmove(n->num + places_rdx, n->num, BC_NUM_SIZE(n->len));
 346                 memset(n->num, 0, BC_NUM_SIZE(places_rdx));
 347         }
 348
 349         n->rdx += places_rdx;
 350         n->scale += places;
 351         n->len += places_rdx;
 352
 353         assert(n->rdx == BC_NUM_RDX(n->scale));
 354 }
 355
 356 static void bc_num_retireMul(BcNum *restrict n, size_t scale,
 357                              bool neg1, bool neg2)
 358 {
 359         if (n->scale < scale) bc_num_extend(n, scale - n->scale);
 360         else bc_num_truncate(n, n->scale - scale);
 361
 362         bc_num_clean(n);
 363         if (BC_NUM_NONZERO(n)) n->neg = (!neg1 != !neg2);
 364 }
 365
 366 static void bc_num_split(const BcNum *restrict n, size_t idx,
 367                          BcNum *restrict a, BcNum *restrict b)
 368 {
 369         assert(BC_NUM_ZERO(a));
 370         assert(BC_NUM_ZERO(b));
 371
 372         if (idx < n->len) {
 373
 374                 b->len = n->len - idx;
 375                 a->len = idx;
 376                 a->scale = a->rdx = b->scale = b->rdx = 0;
 377
 378                 assert(a->cap >= a->len);
 379                 assert(b->cap >= b->len);
 380
 381                 memcpy(b->num, n->num + idx, BC_NUM_SIZE(b->len));
 382                 memcpy(a->num, n->num, BC_NUM_SIZE(idx));
 383
 384                 bc_num_clean(b);
 385         }
 386         else bc_num_copy(a, n);
 387
 388         bc_num_clean(a);
 389 }
 390
 391 static size_t bc_num_shiftZero(BcNum *restrict n) {
 392
 393         size_t i;
 394
 395         assert(!n->rdx || BC_NUM_ZERO(n));
 396
 397         for (i = 0; i < n->len && !n->num[i]; ++i);
 398
 399         n->len -= i;
 400         n->num += i;
 401
 402         return i;
 403 }
 404
 405 static void bc_num_unshiftZero(BcNum *restrict n, size_t places_rdx) {
 406         n->len += places_rdx;
 407         n->num -= places_rdx;
 408 }
 409
 410 static void bc_num_shift(BcNum *restrict n, BcBigDig dig) {
 411
 412         size_t i, len = n->len;
 413         BcBigDig carry = 0, pow;
 414         BcDig *ptr = n->num;
 415
 416         assert(dig < BC_BASE_DIGS);
 417
 418         pow = bc_num_pow10[dig];
 419         dig = bc_num_pow10[BC_BASE_DIGS - dig];
 420
 421         for (i = len - 1; i < len; --i) {
 422                 BcBigDig in, temp;
 423                 in = ((BcBigDig) ptr[i]);
 424                 temp = carry * dig;
 425                 carry = in % pow;
 426                 ptr[i] = ((BcDig) (in / pow)) + (BcDig) temp;
 427         }
 428
 429         assert(!carry);
 430 }
 431
 432 static void bc_num_shiftLeft(BcNum *restrict n, size_t places) {
 433
 434         BcBigDig dig;
 435         size_t places_rdx;
 436         bool shift;
 437
 438         if (!places) return;
 439         if (places > n->scale) {
 440                 size_t size = bc_vm_growSize(BC_NUM_RDX(places - n->scale), n->len);
 441                 if (size > SIZE_MAX - 1) bc_vm_err(BC_ERROR_MATH_OVERFLOW);
 442         }
 443         if (BC_NUM_ZERO(n)) {
 444                 if (n->scale >= places) n->scale -= places;
 445                 else n->scale = 0;
 446                 return;
 447         }
 448
 449         dig = (BcBigDig) (places % BC_BASE_DIGS);
 450         shift = (dig != 0);
 451         places_rdx = BC_NUM_RDX(places);
 452
 453         if (n->scale) {
 454
 455                 if (n->rdx >= places_rdx) {
 456
 457                         size_t mod = n->scale % BC_BASE_DIGS, revdig;
 458
 459                         mod = mod ? mod : BC_BASE_DIGS;
 460                         revdig = dig ? BC_BASE_DIGS - dig : 0;
 461
 462                         if (mod + revdig > BC_BASE_DIGS) places_rdx = 1;
 463                         else places_rdx = 0;
 464                 }
 465                 else places_rdx -= n->rdx;
 466         }
 467
 468         if (places_rdx) {
 469                 bc_num_expand(n, bc_vm_growSize(n->len, places_rdx));
 470                 memmove(n->num + places_rdx, n->num, BC_NUM_SIZE(n->len));
 471                 memset(n->num, 0, BC_NUM_SIZE(places_rdx));
 472                 n->len += places_rdx;
 473         }
 474
 475         if (places > n->scale) n->scale = n->rdx = 0;
 476         else {
 477                 n->scale -= places;
 478                 n->rdx = BC_NUM_RDX(n->scale);
 479         }
 480
 481         if (shift) bc_num_shift(n, BC_BASE_DIGS - dig);
 482
 483         bc_num_clean(n);
 484 }
 485
 486 static void bc_num_shiftRight(BcNum *restrict n, size_t places) {
 487
 488         BcBigDig dig;
 489         size_t places_rdx, scale, scale_mod, int_len, expand;
 490         bool shift;
 491
 492         if (!places) return;
 493         if (BC_NUM_ZERO(n)) {
 494                 n->scale += places;
 495                 bc_num_expand(n, BC_NUM_RDX(n->scale));
 496                 return;
 497         }
 498
 499         dig = (BcBigDig) (places % BC_BASE_DIGS);
 500         shift = (dig != 0);
 501         scale = n->scale;
 502         scale_mod = scale % BC_BASE_DIGS;
 503         scale_mod = scale_mod ? scale_mod : BC_BASE_DIGS;
 504         int_len = bc_num_int(n);
 505         places_rdx = BC_NUM_RDX(places);
 506
 507         if (scale_mod + dig > BC_BASE_DIGS) {
 508                 expand = places_rdx - 1;
 509                 places_rdx = 1;
 510         }
 511         else {
 512                 expand = places_rdx;
 513                 places_rdx = 0;
 514         }
 515
 516         if (expand > int_len) expand -= int_len;
 517         else expand = 0;
 518
 519         bc_num_extend(n, places_rdx * BC_BASE_DIGS);
 520         bc_num_expand(n, bc_vm_growSize(expand, n->len));
 521         memset(n->num + n->len, 0, BC_NUM_SIZE(expand));
 522         n->len += expand;
 523         n->scale = n->rdx = 0;
 524
 525         if (shift) bc_num_shift(n, dig);
 526
 527         n->scale = scale + places;
 528         n->rdx = BC_NUM_RDX(n->scale);
 529
 530         bc_num_clean(n);
 531
 532         assert(n->rdx <= n->len && n->len <= n->cap);
 533         assert(n->rdx == BC_NUM_RDX(n->scale));
 534 }
 535
 536 static void bc_num_inv(BcNum *a, BcNum *b, size_t scale) {
 537
 538         BcNum one;
 539         BcDig num[2];
 540
 541         assert(BC_NUM_NONZERO(a));
 542
 543         bc_num_setup(&one, num, sizeof(num) / sizeof(BcDig));
 544         bc_num_one(&one);
 545
 546         bc_num_div(&one, a, b, scale);
 547 }
 548
 549 #if BC_ENABLE_EXTRA_MATH
 550 static void bc_num_intop(const BcNum *a, const BcNum *b, BcNum *restrict c,
 551                          BcBigDig *v)
 552 {
 553         if (BC_ERR(b->rdx)) bc_vm_err(BC_ERROR_MATH_NON_INTEGER);
 554         bc_num_copy(c, a);
 555         bc_num_bigdig(b, v);
 556 }
 557 #endif // BC_ENABLE_EXTRA_MATH
 558
 559 static void bc_num_as(BcNum *a, BcNum *b, BcNum *restrict c, size_t sub) {
 560
 561         BcDig *ptr_c, *ptr_l, *ptr_r;
 562         size_t i, min_rdx, max_rdx, diff, a_int, b_int, min_len, max_len, max_int;
 563         size_t len_l, len_r;
 564         bool b_neg, do_sub, do_rev_sub, carry;
 565
 566         // Because this function doesn't need to use scale (per the bc spec),
 567         // I am hijacking it to say whether it's doing an add or a subtract.
 568         // Convert substraction to addition of negative second operand.
 569
 570         if (BC_NUM_ZERO(b)) {
 571                 bc_num_copy(c, a);
 572                 return;
 573         }
 574         if (BC_NUM_ZERO(a)) {
 575                 bc_num_copy(c, b);
 576                 c->neg = (b->neg != sub);
 577                 return;
 578         }
 579
 580         // Invert sign of b if it is to be subtracted. This operation must
 581         // preced the tests for any of the operands being zero.
 582         b_neg = (b->neg != sub);
 583
 584         // Actually add the numbers if their signs are equal, else subtract.
 585         do_sub = (a->neg != b_neg);
 586
 587         a_int = bc_num_int(a);
 588         b_int = bc_num_int(b);
 589         max_int = BC_MAX(a_int, b_int);
 590
 591         min_rdx = BC_MIN(a->rdx, b->rdx);
 592         max_rdx = BC_MAX(a->rdx, b->rdx);
 593         diff = max_rdx - min_rdx;
 594
 595         max_len = max_int + max_rdx;
 596
 597         if (do_sub) {
 598
 599                 // Check whether b has to be subtracted from a or a from b.
 600                 if (a_int != b_int) do_rev_sub = (a_int < b_int);
 601                 else if (a->rdx > b->rdx)
 602                         do_rev_sub = (bc_num_compare(a->num + diff, b->num, b->len) < 0);
 603                 else
 604                         do_rev_sub = (bc_num_compare(a->num, b->num + diff, a->len) <= 0);
 605         }
 606         else {
 607
 608                 // The result array of the addition might come out one element
 609                 // longer than the bigger of the operand arrays.
 610                 max_len += 1;
 611                 do_rev_sub = (a_int < b_int);
 612         }
 613
 614         assert(max_len <= c->cap);
 615
 616         if (do_rev_sub) {
 617                 ptr_l = b->num;
 618                 ptr_r = a->num;
 619                 len_l = b->len;
 620                 len_r = a->len;
 621         }
 622         else {
 623                 ptr_l = a->num;
 624                 ptr_r = b->num;
 625                 len_l = a->len;
 626                 len_r = b->len;
 627         }
 628
 629         ptr_c = c->num;
 630         carry = false;
 631
 632         if (diff) {
 633
 634                 // If the rdx values of the operands do not match, the result will
 635                 // have low end elements that are the positive or negative trailing
 636                 // elements of the operand with higher rdx value.
 637                 if ((a->rdx > b->rdx) != do_rev_sub) {
 638
 639                         // !do_rev_sub && a->rdx > b->rdx || do_rev_sub && b->rdx > a->rdx
 640                         // The left operand has BcDig values that need to be copied,
 641                         // either from a or from b (in case of a reversed subtraction).
 642                         memcpy(ptr_c, ptr_l, BC_NUM_SIZE(diff));
 643                         ptr_l += diff;
 644                         len_l -= diff;
 645                 }
 646                 else {
 647
 648                         // The right operand has BcDig values that need to be copied
 649                         // or subtracted from zero (in case of a subtraction).
 650                         if (do_sub) {
 651
 652                                 // do_sub (do_rev_sub && a->rdx > b->rdx ||
 653                                 // !do_rev_sub && b->rdx > a->rdx)
 654                                 for (i = 0; i < diff; i++)
 655                                         ptr_c[i] = bc_num_subDigits(0, ptr_r[i], &carry);
 656                         }
 657                         else {
 658
 659                                 // !do_sub && b->rdx > a->rdx
 660                                 memcpy(ptr_c, ptr_r, BC_NUM_SIZE(diff));
 661                         }
 662
 663                         ptr_r += diff;
 664                         len_r -= diff;
 665                 }
 666
 667                 ptr_c += diff;
 668         }
 669
 670         min_len = BC_MIN(len_l, len_r);
 671
 672         // After dealing with possible low array elements that depend on only one
 673         // operand, the actual add or subtract can be performed as if the rdx of
 674         // both operands was the same.
 675         // Inlining takes care of eliminating constant zero arguments to
 676         // addDigit/subDigit (checked in disassembly of resulting bc binary
 677         // compiled with gcc and clang).
 678         if (do_sub) {
 679                 for (i = 0; i < min_len; ++i)
 680                         ptr_c[i] = bc_num_subDigits(ptr_l[i], ptr_r[i], &carry);
 681                 for (; i < len_l; ++i) ptr_c[i] = bc_num_subDigits(ptr_l[i], 0, &carry);
 682         }
 683         else {
 684                 for (i = 0; i < min_len; ++i)
 685                         ptr_c[i] = bc_num_addDigits(ptr_l[i], ptr_r[i], &carry);
 686                 for (; i < len_l; ++i) ptr_c[i] = bc_num_addDigits(ptr_l[i], 0, &carry);
 687                 ptr_c[i] = bc_num_addDigits(0, 0, &carry);
 688         }
 689
 690         assert(carry == false);
 691
 692         // The result has the same sign as a, unless the operation was a
 693         // reverse subtraction (b - a).
 694         c->neg = (a->neg != (do_sub && do_rev_sub));
 695         c->len = max_len;
 696         c->rdx = max_rdx;
 697         c->scale = BC_MAX(a->scale, b->scale);
 698
 699         bc_num_clean(c);
 700 }
 701
 702 static void bc_num_m_simp(const BcNum *a, const BcNum *b, BcNum *restrict c)
 703 {
 704         size_t i, alen = a->len, blen = b->len, clen;
 705         BcDig *ptr_a = a->num, *ptr_b = b->num, *ptr_c;
 706         BcBigDig sum = 0, carry = 0;
 707
 708         assert(sizeof(sum) >= sizeof(BcDig) * 2);
 709         assert(!a->rdx && !b->rdx);
 710
 711         clen = bc_vm_growSize(alen, blen);
 712         bc_num_expand(c, bc_vm_growSize(clen, 1));
 713
 714         ptr_c = c->num;
 715         memset(ptr_c, 0, BC_NUM_SIZE(c->cap));
 716
 717         for (i = 0; i < clen; ++i) {
 718
 719                 ssize_t sidx = (ssize_t) (i - blen + 1);
 720                 size_t j = (size_t) BC_MAX(0, sidx), k = BC_MIN(i, blen - 1);
 721
 722                 for (; j < alen && k < blen; ++j, --k) {
 723
 724                         sum += ((BcBigDig) ptr_a[j]) * ((BcBigDig) ptr_b[k]);
 725
 726                         if (sum >= ((BcBigDig) BC_BASE_POW) * BC_BASE_POW) {
 727                                 carry += sum / BC_BASE_POW;
 728                                 sum %= BC_BASE_POW;
 729                         }
 730                 }
 731
 732                 if (sum >= BC_BASE_POW) {
 733                         carry += sum / BC_BASE_POW;
 734                         sum %= BC_BASE_POW;
 735                 }
 736
 737                 ptr_c[i] = (BcDig) sum;
 738                 assert(ptr_c[i] < BC_BASE_POW);
 739                 sum = carry;
 740                 carry = 0;
 741         }
 742
 743         // This should always be true because there should be no carry on the last
 744         // digit; multiplication never goes above the sum of both lengths.
 745         assert(!sum);
 746
 747         c->len = clen;
 748 }
 749
 750 static void bc_num_shiftAddSub(BcNum *restrict n, const BcNum *restrict a,
 751                                size_t shift, BcNumShiftAddOp op)
 752 {
 753         assert(n->len >= shift + a->len);
 754         assert(!n->rdx && !a->rdx);
 755         op(n->num + shift, a->num, a->len);
 756 }
 757
 758 static void bc_num_k(BcNum *a, BcNum *b, BcNum *restrict c) {
 759
 760         size_t max, max2, total;
 761         BcNum l1, h1, l2, h2, m2, m1, z0, z1, z2, temp;
 762         BcDig *digs, *dig_ptr;
 763         BcNumShiftAddOp op;
 764         bool aone = BC_NUM_ONE(a);
 765
 766         assert(BC_NUM_ZERO(c));
 767
 768         if (BC_NUM_ZERO(a) || BC_NUM_ZERO(b)) return;
 769         if (aone || BC_NUM_ONE(b)) {
 770                 bc_num_copy(c, aone ? b : a);
 771                 if ((aone && a->neg) || b->neg) c->neg = !c->neg;
 772                 return;
 773         }
 774         if (a->len < BC_NUM_KARATSUBA_LEN || b->len < BC_NUM_KARATSUBA_LEN) {
 775                 bc_num_m_simp(a, b, c);
 776                 return;
 777         }
 778
 779         max = BC_MAX(a->len, b->len);
 780         max = BC_MAX(max, BC_NUM_DEF_SIZE);
 781         max2 = (max + 1) / 2;
 782
 783         total = bc_vm_arraySize(BC_NUM_KARATSUBA_ALLOCS, max);
 784
 785         BC_SIG_LOCK;
 786
 787         digs = dig_ptr = bc_vm_malloc(BC_NUM_SIZE(total));
 788
 789         bc_num_setup(&l1, dig_ptr, max);
 790         dig_ptr += max;
 791         bc_num_setup(&h1, dig_ptr, max);
 792         dig_ptr += max;
 793         bc_num_setup(&l2, dig_ptr, max);
 794         dig_ptr += max;
 795         bc_num_setup(&h2, dig_ptr, max);
 796         dig_ptr += max;
 797         bc_num_setup(&m1, dig_ptr, max);
 798         dig_ptr += max;
 799         bc_num_setup(&m2, dig_ptr, max);
 800         max = bc_vm_growSize(max, 1);
 801         bc_num_init(&z0, max);
 802         bc_num_init(&z1, max);
 803         bc_num_init(&z2, max);
 804         max = bc_vm_growSize(max, max) + 1;
 805         bc_num_init(&temp, max);
 806
 807         BC_SETJMP_LOCKED(err);
 808
 809         BC_SIG_UNLOCK;
 810
 811         bc_num_split(a, max2, &l1, &h1);
 812         bc_num_split(b, max2, &l2, &h2);
 813
 814         bc_num_expand(c, max);
 815         c->len = max;
 816         memset(c->num, 0, BC_NUM_SIZE(c->len));
 817
 818         bc_num_sub(&h1, &l1, &m1, 0);
 819         bc_num_sub(&l2, &h2, &m2, 0);
 820
 821         if (BC_NUM_NONZERO(&h1) && BC_NUM_NONZERO(&h2)) {
 822
 823                 bc_num_m(&h1, &h2, &z2, 0);
 824                 bc_num_clean(&z2);
 825
 826                 bc_num_shiftAddSub(c, &z2, max2 * 2, bc_num_addArrays);
 827                 bc_num_shiftAddSub(c, &z2, max2, bc_num_addArrays);
 828         }
 829
 830         if (BC_NUM_NONZERO(&l1) && BC_NUM_NONZERO(&l2)) {
 831
 832                 bc_num_m(&l1, &l2, &z0, 0);
 833                 bc_num_clean(&z0);
 834
 835                 bc_num_shiftAddSub(c, &z0, max2, bc_num_addArrays);
 836                 bc_num_shiftAddSub(c, &z0, 0, bc_num_addArrays);
 837         }
 838
 839         if (BC_NUM_NONZERO(&m1) && BC_NUM_NONZERO(&m2)) {
 840
 841                 bc_num_m(&m1, &m2, &z1, 0);
 842                 bc_num_clean(&z1);
 843
 844                 op = (m1.neg != m2.neg) ? bc_num_subArrays : bc_num_addArrays;
 845                 bc_num_shiftAddSub(c, &z1, max2, op);
 846         }
 847
 848 err:
 849         BC_SIG_MAYLOCK;
 850         free(digs);
 851         bc_num_free(&temp);
 852         bc_num_free(&z2);
 853         bc_num_free(&z1);
 854         bc_num_free(&z0);
 855         BC_LONGJMP_CONT;
 856 }
 857
 858 static void bc_num_m(BcNum *a, BcNum *b, BcNum *restrict c, size_t scale) {
 859
 860         BcNum cpa, cpb;
 861         size_t ascale, bscale, ardx, brdx, azero = 0, bzero = 0, zero, len, rscale;
 862
 863         bc_num_zero(c);
 864         ascale = a->scale;
 865         bscale = b->scale;
 866         scale = BC_MAX(scale, ascale);
 867         scale = BC_MAX(scale, bscale);
 868
 869         rscale = ascale + bscale;
 870         scale = BC_MIN(rscale, scale);
 871
 872         if ((a->len == 1 || b->len == 1) && !a->rdx && !b->rdx) {
 873
 874                 BcNum *operand;
 875                 BcBigDig dig;
 876
 877                 if (a->len == 1) {
 878                         dig = (BcBigDig) a->num[0];
 879                         operand = b;
 880                 }
 881                 else {
 882                         dig = (BcBigDig) b->num[0];
 883                         operand = a;
 884                 }
 885
 886                 bc_num_mulArray(operand, dig, c);
 887
 888                 if (BC_NUM_NONZERO(c)) c->neg = (a->neg != b->neg);
 889
 890                 return;
 891         }
 892
 893         BC_SIG_LOCK;
 894
 895         bc_num_init(&cpa, a->len + a->rdx);
 896         bc_num_init(&cpb, b->len + b->rdx);
 897
 898         BC_SETJMP_LOCKED(err);
 899
 900         BC_SIG_UNLOCK;
 901
 902         bc_num_copy(&cpa, a);
 903         bc_num_copy(&cpb, b);
 904
 905         cpa.neg = cpb.neg = false;
 906
 907         ardx = cpa.rdx * BC_BASE_DIGS;
 908         bc_num_shiftLeft(&cpa, ardx);
 909
 910         brdx = cpb.rdx * BC_BASE_DIGS;
 911         bc_num_shiftLeft(&cpb, brdx);
 912
 913         // We need to reset the jump here because azero and bzero are used in the
 914         // cleanup, and local variables are not guaranteed to be the same after a
 915         // jump.
 916         BC_SIG_LOCK;
 917
 918         BC_UNSETJMP;
 919
 920         azero = bc_num_shiftZero(&cpa);
 921         bzero = bc_num_shiftZero(&cpb);
 922
 923         BC_SETJMP_LOCKED(err);
 924
 925         BC_SIG_UNLOCK;
 926
 927         bc_num_clean(&cpa);
 928         bc_num_clean(&cpb);
 929
 930         bc_num_k(&cpa, &cpb, c);
 931
 932         zero = bc_vm_growSize(azero, bzero);
 933         len = bc_vm_growSize(c->len, zero);
 934
 935         bc_num_expand(c, len);
 936         bc_num_shiftLeft(c, (len - c->len) * BC_BASE_DIGS);
 937         bc_num_shiftRight(c, ardx + brdx);
 938
 939         bc_num_retireMul(c, scale, a->neg, b->neg);
 940
 941 err:
 942         BC_SIG_MAYLOCK;
 943         bc_num_unshiftZero(&cpb, bzero);
 944         bc_num_free(&cpb);
 945         bc_num_unshiftZero(&cpa, azero);
 946         bc_num_free(&cpa);
 947         BC_LONGJMP_CONT;
 948 }
 949
 950 static bool bc_num_nonZeroDig(BcDig *restrict a, size_t len) {
 951         size_t i;
 952         bool nonzero = false;
 953         for (i = len - 1; !nonzero && i < len; --i) nonzero = (a[i] != 0);
 954         return nonzero;
 955 }
 956
 957 static ssize_t bc_num_divCmp(const BcDig *a, const BcNum *b, size_t len) {
 958
 959         ssize_t cmp;
 960
 961         if (b->len > len && a[len]) cmp = bc_num_compare(a, b->num, len + 1);
 962         else if (b->len <= len) {
 963                 if (a[len]) cmp = 1;
 964                 else cmp = bc_num_compare(a, b->num, len);
 965         }
 966         else cmp = -1;
 967
 968         return cmp;
 969 }
 970
 971 static void bc_num_divExtend(BcNum *restrict a, BcNum *restrict b,
 972                              BcBigDig divisor)
 973 {
 974         size_t pow;
 975
 976         assert(divisor < BC_BASE_POW);
 977
 978         pow = BC_BASE_DIGS - bc_num_log10((size_t) divisor);
 979
 980         bc_num_shiftLeft(a, pow);
 981         bc_num_shiftLeft(b, pow);
 982 }
 983
 984 static void bc_num_d_long(BcNum *restrict a, BcNum *restrict b,
 985                           BcNum *restrict c, size_t scale)
 986 {
 987         BcBigDig divisor;
 988         size_t len, end, i, rdx;
 989         BcNum cpb;
 990         bool nonzero = false;
 991
 992         assert(b->len < a->len);
 993         len = b->len;
 994         end = a->len - len;
 995         assert(len >= 1);
 996
 997         bc_num_expand(c, a->len);
 998         memset(c->num, 0, c->cap * sizeof(BcDig));
 999
1000         c->rdx = a->rdx;
1001         c->scale = a->scale;
1002         c->len = a->len;
1003
1004         divisor = (BcBigDig) b->num[len - 1];
1005
1006         if (len > 1 && bc_num_nonZeroDig(b->num, len - 1)) {
1007
1008                 nonzero = (divisor > 1 << ((10 * BC_BASE_DIGS) / 6 + 1));
1009
1010                 if (!nonzero) {
1011
1012                         bc_num_divExtend(a, b, divisor);
1013
1014                         len = BC_MAX(a->len, b->len);
1015                         bc_num_expand(a, len + 1);
1016
1017                         if (len + 1 > a->len) a->len = len + 1;
1018
1019                         len = b->len;
1020                         end = a->len - len;
1021                         divisor = (BcBigDig) b->num[len - 1];
1022
1023                         nonzero = bc_num_nonZeroDig(b->num, len - 1);
1024                 }
1025         }
1026
1027         divisor += nonzero;
1028
1029         bc_num_expand(c, a->len);
1030         memset(c->num, 0, BC_NUM_SIZE(c->cap));
1031
1032         assert(c->scale >= scale);
1033         rdx = c->rdx - BC_NUM_RDX(scale);
1034
1035         BC_SIG_LOCK;
1036
1037         bc_num_init(&cpb, len + 1);
1038
1039         BC_SETJMP_LOCKED(err);
1040
1041         BC_SIG_UNLOCK;
1042
1043         i = end - 1;
1044
1045         for (; i < end && i >= rdx && BC_NUM_NONZERO(a); --i) {
1046
1047                 ssize_t cmp;
1048                 BcDig *n;
1049                 BcBigDig result;
1050
1051                 n = a->num + i;
1052                 assert(n >= a->num);
1053                 result = 0;
1054
1055                 cmp = bc_num_divCmp(n, b, len);
1056
1057                 while (cmp >= 0) {
1058
1059                         BcBigDig n1, dividend, q;
1060
1061                         n1 = (BcBigDig) n[len];
1062                         dividend = n1 * BC_BASE_POW + (BcBigDig) n[len - 1];
1063                         q = (dividend / divisor);
1064
1065                         if (q <= 1) {
1066                                 q = 1;
1067                                 bc_num_subArrays(n, b->num, len);
1068                         }
1069                         else {
1070
1071                                 assert(q <= BC_BASE_POW);
1072
1073                                 bc_num_mulArray(b, (BcBigDig) q, &cpb);
1074                                 bc_num_subArrays(n, cpb.num, cpb.len);
1075                         }
1076
1077                         result += q;
1078                         assert(result <= BC_BASE_POW);
1079
1080                         if (nonzero) cmp = bc_num_divCmp(n, b, len);
1081                         else cmp = -1;
1082                 }
1083
1084                 assert(result < BC_BASE_POW);
1085
1086                 c->num[i] = (BcDig) result;
1087         }
1088
1089 err:
1090         BC_SIG_MAYLOCK;
1091         bc_num_free(&cpb);
1092         BC_LONGJMP_CONT;
1093 }
1094
1095 static void bc_num_d(BcNum *a, BcNum *b, BcNum *restrict c, size_t scale) {
1096
1097         size_t len;
1098         BcNum cpa, cpb;
1099
1100         if (BC_NUM_ZERO(b)) bc_vm_err(BC_ERROR_MATH_DIVIDE_BY_ZERO);
1101         if (BC_NUM_ZERO(a)) {
1102                 bc_num_setToZero(c, scale);
1103                 return;
1104         }
1105         if (BC_NUM_ONE(b)) {
1106                 bc_num_copy(c, a);
1107                 bc_num_retireMul(c, scale, a->neg, b->neg);
1108                 return;
1109         }
1110         if (!a->rdx && !b->rdx && b->len == 1 && !scale) {
1111                 BcBigDig rem;
1112                 bc_num_divArray(a, (BcBigDig) b->num[0], c, &rem);
1113                 bc_num_retireMul(c, scale, a->neg, b->neg);
1114                 return;
1115         }
1116
1117         len = bc_num_mulReq(a, b, scale);
1118
1119         BC_SIG_LOCK;
1120
1121         bc_num_init(&cpa, len);
1122         bc_num_copy(&cpa, a);
1123         bc_num_createCopy(&cpb, b);
1124
1125         BC_SETJMP_LOCKED(err);
1126
1127         BC_SIG_UNLOCK;
1128
1129         len = b->len;
1130
1131         if (len > cpa.len) {
1132                 bc_num_expand(&cpa, bc_vm_growSize(len, 2));
1133                 bc_num_extend(&cpa, (len - cpa.len) * BC_BASE_DIGS);
1134         }
1135
1136         cpa.scale = cpa.rdx * BC_BASE_DIGS;
1137
1138         bc_num_extend(&cpa, b->scale);
1139         cpa.rdx -= BC_NUM_RDX(b->scale);
1140         cpa.scale = cpa.rdx * BC_BASE_DIGS;
1141
1142         if (scale > cpa.scale) {
1143                 bc_num_extend(&cpa, scale);
1144                 cpa.scale = cpa.rdx * BC_BASE_DIGS;
1145         }
1146
1147         if (cpa.cap == cpa.len) bc_num_expand(&cpa, bc_vm_growSize(cpa.len, 1));
1148
1149         // We want an extra zero in front to make things simpler.
1150         cpa.num[cpa.len++] = 0;
1151
1152         if (cpa.rdx == cpa.len) cpa.len = bc_num_nonzeroLen(&cpa);
1153         if (cpb.rdx == cpb.len) cpb.len = bc_num_nonzeroLen(&cpb);
1154         cpb.scale = cpb.rdx = 0;
1155
1156         bc_num_d_long(&cpa, &cpb, c, scale);
1157
1158         bc_num_retireMul(c, scale, a->neg, b->neg);
1159
1160 err:
1161         BC_SIG_MAYLOCK;
1162         bc_num_free(&cpb);
1163         bc_num_free(&cpa);
1164         BC_LONGJMP_CONT;
1165 }
1166
1167 static void bc_num_r(BcNum *a, BcNum *b, BcNum *restrict c,
1168                      BcNum *restrict d, size_t scale, size_t ts)
1169 {
1170         BcNum temp;
1171         bool neg;
1172
1173         if (BC_NUM_ZERO(b)) bc_vm_err(BC_ERROR_MATH_DIVIDE_BY_ZERO);
1174         if (BC_NUM_ZERO(a)) {
1175                 bc_num_setToZero(c, ts);
1176                 bc_num_setToZero(d, ts);
1177                 return;
1178         }
1179
1180         BC_SIG_LOCK;
1181
1182         bc_num_init(&temp, d->cap);
1183
1184         BC_SETJMP_LOCKED(err);
1185
1186         BC_SIG_UNLOCK;
1187
1188         bc_num_d(a, b, c, scale);
1189
1190         if (scale) scale = ts + 1;
1191
1192         bc_num_m(c, b, &temp, scale);
1193         bc_num_sub(a, &temp, d, scale);
1194
1195         if (ts > d->scale && BC_NUM_NONZERO(d)) bc_num_extend(d, ts - d->scale);
1196
1197         neg = d->neg;
1198         bc_num_retireMul(d, ts, a->neg, b->neg);
1199         d->neg = BC_NUM_NONZERO(d) ? neg : false;
1200
1201 err:
1202         BC_SIG_MAYLOCK;
1203         bc_num_free(&temp);
1204         BC_LONGJMP_CONT;
1205 }
1206
1207 static void bc_num_rem(BcNum *a, BcNum *b, BcNum *restrict c, size_t scale) {
1208
1209         BcNum c1;
1210         size_t ts;
1211
1212         ts = bc_vm_growSize(scale, b->scale);
1213         ts = BC_MAX(ts, a->scale);
1214
1215         BC_SIG_LOCK;
1216
1217         bc_num_init(&c1, bc_num_mulReq(a, b, ts));
1218
1219         BC_SETJMP_LOCKED(err);
1220
1221         BC_SIG_UNLOCK;
1222
1223         bc_num_r(a, b, &c1, c, scale, ts);
1224
1225 err:
1226         BC_SIG_MAYLOCK;
1227         bc_num_free(&c1);
1228         BC_LONGJMP_CONT;
1229 }
1230
1231 static void bc_num_p(BcNum *a, BcNum *b, BcNum *restrict c, size_t scale) {
1232
1233         BcNum copy;
1234         BcBigDig pow = 0;
1235         size_t i, powrdx, resrdx;
1236         bool neg, zero;
1237
1238         if (BC_ERR(b->rdx)) bc_vm_err(BC_ERROR_MATH_NON_INTEGER);
1239
1240         if (BC_NUM_ZERO(b)) {
1241                 bc_num_one(c);
1242                 return;
1243         }
1244         if (BC_NUM_ZERO(a)) {
1245                 if (b->neg) bc_vm_err(BC_ERROR_MATH_DIVIDE_BY_ZERO);
1246                 bc_num_setToZero(c, scale);
1247                 return;
1248         }
1249         if (BC_NUM_ONE(b)) {
1250                 if (!b->neg) bc_num_copy(c, a);
1251                 else bc_num_inv(a, c, scale);
1252                 return;
1253         }
1254
1255         BC_SIG_LOCK;
1256
1257         neg = b->neg;
1258         b->neg = false;
1259         bc_num_bigdig(b, &pow);
1260         b->neg = neg;
1261
1262         bc_num_createCopy(&copy, a);
1263
1264         BC_SETJMP_LOCKED(err);
1265
1266         BC_SIG_UNLOCK;
1267
1268         if (!neg) {
1269                 size_t max = BC_MAX(scale, a->scale), scalepow = a->scale * pow;
1270                 scale = BC_MIN(scalepow, max);
1271         }
1272
1273         for (powrdx = a->scale; !(pow & 1); pow >>= 1) {
1274                 powrdx <<= 1;
1275                 bc_num_mul(&copy, &copy, &copy, powrdx);
1276         }
1277
1278         bc_num_copy(c, &copy);
1279         resrdx = powrdx;
1280
1281         while (pow >>= 1) {
1282
1283                 powrdx <<= 1;
1284                 bc_num_mul(&copy, &copy, &copy, powrdx);
1285
1286                 if (pow & 1) {
1287                         resrdx += powrdx;
1288                         bc_num_mul(c, &copy, c, resrdx);
1289                 }
1290         }
1291
1292         if (neg) bc_num_inv(c, c, scale);
1293
1294         if (c->scale > scale) bc_num_truncate(c, c->scale - scale);
1295
1296         // We can't use bc_num_clean() here.
1297         for (zero = true, i = 0; zero && i < c->len; ++i) zero = !c->num[i];
1298         if (zero) bc_num_setToZero(c, scale);
1299
1300 err:
1301         BC_SIG_MAYLOCK;
1302         bc_num_free(&copy);
1303         BC_LONGJMP_CONT;
1304 }
1305
1306 #if BC_ENABLE_EXTRA_MATH
1307 static void bc_num_place(BcNum *a, BcNum *b, BcNum *restrict c, size_t scale) {
1308
1309         BcBigDig val = 0;
1310
1311         BC_UNUSED(scale);
1312
1313         bc_num_intop(a, b, c, &val);
1314
1315         if (val < c->scale) bc_num_truncate(c, c->scale - val);
1316         else if (val > c->scale) bc_num_extend(c, val - c->scale);
1317 }
1318
1319 static void bc_num_left(BcNum *a, BcNum *b, BcNum *restrict c, size_t scale) {
1320
1321         BcBigDig val = 0;
1322
1323         BC_UNUSED(scale);
1324
1325         bc_num_intop(a, b, c, &val);
1326
1327         bc_num_shiftLeft(c, (size_t) val);
1328 }
1329
1330 static void bc_num_right(BcNum *a, BcNum *b, BcNum *restrict c, size_t scale) {
1331
1332         BcBigDig val = 0;
1333
1334         BC_UNUSED(scale);
1335
1336         bc_num_intop(a, b, c, &val);
1337
1338         if (BC_NUM_ZERO(c)) return;
1339
1340         bc_num_shiftRight(c, (size_t) val);
1341 }
1342 #endif // BC_ENABLE_EXTRA_MATH
1343
1344 static void bc_num_binary(BcNum *a, BcNum *b, BcNum *c, size_t scale,
1345                           BcNumBinaryOp op, size_t req)
1346 {
1347         BcNum num2, *ptr_a, *ptr_b;
1348         bool init = false;
1349
1350         assert(a != NULL && b != NULL && c != NULL && op != NULL);
1351
1352         BC_SIG_LOCK;
1353
1354         if (c == a) {
1355
1356                 ptr_a = &num2;
1357
1358                 memcpy(ptr_a, c, sizeof(BcNum));
1359                 init = true;
1360         }
1361         else ptr_a = a;
1362
1363         if (c == b) {
1364
1365                 ptr_b = &num2;
1366
1367                 if (c != a) {
1368                         memcpy(ptr_b, c, sizeof(BcNum));
1369                         init = true;
1370                 }
1371         }
1372         else ptr_b = b;
1373
1374         if (init) {
1375
1376                 bc_num_init(c, req);
1377
1378                 BC_SETJMP_LOCKED(err);
1379                 BC_SIG_UNLOCK;
1380         }
1381         else {
1382                 BC_SIG_UNLOCK;
1383                 bc_num_expand(c, req);
1384         }
1385
1386         op(ptr_a, ptr_b, c, scale);
1387
1388         assert(!c->neg || BC_NUM_NONZERO(c));
1389         assert(c->rdx <= c->len || !c->len);
1390         assert(!c->len || c->num[c->len - 1] || c->rdx == c->len);
1391
1392 err:
1393         if (init) {
1394                 BC_SIG_MAYLOCK;
1395                 bc_num_free(&num2);
1396                 BC_LONGJMP_CONT;
1397         }
1398 }
1399
1400 #ifndef NDEBUG
1401 static bool bc_num_strValid(const char *val) {
1402
1403         bool radix = false;
1404         size_t i, len = strlen(val);
1405
1406         if (!len) return true;
1407
1408         for (i = 0; i < len; ++i) {
1409
1410                 BcDig c = val[i];
1411
1412                 if (c == '.') {
1413
1414                         if (radix) return false;
1415
1416                         radix = true;
1417                         continue;
1418                 }
1419
1420                 if (!(isdigit(c) || isupper(c))) return false;
1421         }
1422
1423         return true;
1424 }
1425 #endif // NDEBUG
1426
1427 static BcBigDig bc_num_parseChar(char c, size_t base_t) {
1428
1429         if (isupper(c)) {
1430                 c = BC_NUM_NUM_LETTER(c);
1431                 c = ((size_t) c) >= base_t ? (char) base_t - 1 : c;
1432         }
1433         else c -= '0';
1434
1435         return (BcBigDig) (uchar) c;
1436 }
1437
1438 static void bc_num_parseDecimal(BcNum *restrict n, const char *restrict val) {
1439
1440         size_t len, i, temp, mod;
1441         const char *ptr;
1442         bool zero = true, rdx;
1443
1444         for (i = 0; val[i] == '0'; ++i);
1445
1446         val += i;
1447         assert(!val[0] || isalnum(val[0]) || val[0] == '.');
1448
1449         // All 0's. We can just return, since this
1450         // procedure expects a virgin (already 0) BcNum.
1451         if (!val[0]) return;
1452
1453         len = strlen(val);
1454
1455         ptr = strchr(val, '.');
1456         rdx = (ptr != NULL);
1457
1458         for (i = 0; i < len && (zero = (val[i] == '0' || val[i] == '.')); ++i);
1459
1460         n->scale = (size_t) (rdx * ((val + len) - (ptr + 1)));
1461         n->rdx = BC_NUM_RDX(n->scale);
1462
1463         i = len - (ptr == val ? 0 : i) - rdx;
1464         temp = BC_NUM_ROUND_POW(i);
1465         mod = n->scale % BC_BASE_DIGS;
1466         i = mod ? BC_BASE_DIGS - mod : 0;
1467         n->len = ((temp + i) / BC_BASE_DIGS);
1468
1469         bc_num_expand(n, n->len);
1470         memset(n->num, 0, BC_NUM_SIZE(n->len));
1471
1472         if (zero) n->len = n->rdx = 0;
1473         else {
1474
1475                 BcBigDig exp, pow;
1476
1477                 assert(i <= BC_NUM_BIGDIG_MAX);
1478
1479                 exp = (BcBigDig) i;
1480                 pow = bc_num_pow10[exp];
1481
1482                 for (i = len - 1; i < len; --i, ++exp) {
1483
1484                         char c = val[i];
1485
1486                         if (c == '.') exp -= 1;
1487                         else {
1488
1489                                 size_t idx = exp / BC_BASE_DIGS;
1490
1491                                 if (isupper(c)) c = '9';
1492                                 n->num[idx] += (((BcBigDig) c) - '0') * pow;
1493
1494                                 if ((exp + 1) % BC_BASE_DIGS == 0) pow = 1;
1495                                 else pow *= BC_BASE;
1496                         }
1497                 }
1498         }
1499 }
1500
1501 static void bc_num_parseBase(BcNum *restrict n, const char *restrict val,
1502                              BcBigDig base)
1503 {
1504         BcNum temp, mult1, mult2, result1, result2, *m1, *m2, *ptr;
1505         char c = 0;
1506         bool zero = true;
1507         BcBigDig v;
1508         size_t i, digs, len = strlen(val);
1509
1510         for (i = 0; zero && i < len; ++i) zero = (val[i] == '.' || val[i] == '0');
1511         if (zero) return;
1512
1513         BC_SIG_LOCK;
1514
1515         bc_num_init(&temp, BC_NUM_BIGDIG_LOG10);
1516         bc_num_init(&mult1, BC_NUM_BIGDIG_LOG10);
1517
1518         BC_SETJMP_LOCKED(int_err);
1519
1520         BC_SIG_UNLOCK;
1521
1522         for (i = 0; i < len && (c = val[i]) && c != '.'; ++i) {
1523
1524                 v = bc_num_parseChar(c, base);
1525
1526                 bc_num_mulArray(n, base, &mult1);
1527                 bc_num_bigdig2num(&temp, v);
1528                 bc_num_add(&mult1, &temp, n, 0);
1529         }
1530
1531         if (i == len && !(c = val[i])) goto int_err;
1532
1533         assert(c == '.');
1534
1535         BC_SIG_LOCK;
1536
1537         BC_UNSETJMP;
1538
1539         bc_num_init(&mult2, BC_NUM_BIGDIG_LOG10);
1540         bc_num_init(&result1, BC_NUM_DEF_SIZE);
1541         bc_num_init(&result2, BC_NUM_DEF_SIZE);
1542         bc_num_one(&mult1);
1543
1544         BC_SETJMP_LOCKED(err);
1545
1546         BC_SIG_UNLOCK;
1547
1548         m1 = &mult1;
1549         m2 = &mult2;
1550
1551         for (i += 1, digs = 0; i < len && (c = val[i]); ++i, ++digs) {
1552
1553                 v = bc_num_parseChar(c, base);
1554
1555                 bc_num_mulArray(&result1, base, &result2);
1556
1557                 bc_num_bigdig2num(&temp, v);
1558                 bc_num_add(&result2, &temp, &result1, 0);
1559                 bc_num_mulArray(m1, base, m2);
1560
1561                 if (m2->len < m2->rdx) m2->len = m2->rdx;
1562
1563                 ptr = m1;
1564                 m1 = m2;
1565                 m2 = ptr;
1566         }
1567
1568         // This one cannot be a divide by 0 because mult starts out at 1, then is
1569         // multiplied by base, and base cannot be 0, so mult cannot be 0.
1570         bc_num_div(&result1, m1, &result2, digs * 2);
1571         bc_num_truncate(&result2, digs);
1572         bc_num_add(n, &result2, n, digs);
1573
1574         if (BC_NUM_NONZERO(n)) {
1575                 if (n->scale < digs) bc_num_extend(n, digs - n->scale);
1576         }
1577         else bc_num_zero(n);
1578
1579 err:
1580         BC_SIG_MAYLOCK;
1581         bc_num_free(&result2);
1582         bc_num_free(&result1);
1583         bc_num_free(&mult2);
1584 int_err:
1585         BC_SIG_MAYLOCK;
1586         bc_num_free(&mult1);
1587         bc_num_free(&temp);
1588         BC_LONGJMP_CONT;
1589 }
1590
1591 static void bc_num_printNewline(void) {
1592         if (vm.nchars >= vm.line_len - 1) {
1593                 bc_vm_putchar('\\');
1594                 bc_vm_putchar('\n');
1595         }
1596 }
1597
1598 static void bc_num_putchar(int c) {
1599         if (c != '\n') bc_num_printNewline();
1600         bc_vm_putchar(c);
1601 }
1602
1603 #if DC_ENABLED
1604 static void bc_num_printChar(size_t n, size_t len, bool rdx) {
1605         BC_UNUSED(rdx);
1606         BC_UNUSED(len);
1607         assert(len == 1);
1608         bc_vm_putchar((uchar) n);
1609 }
1610 #endif // DC_ENABLED
1611
1612 static void bc_num_printDigits(size_t n, size_t len, bool rdx) {
1613
1614         size_t exp, pow;
1615
1616         bc_num_putchar(rdx ? '.' : ' ');
1617
1618         for (exp = 0, pow = 1; exp < len - 1; ++exp, pow *= BC_BASE);
1619
1620         for (exp = 0; exp < len; pow /= BC_BASE, ++exp) {
1621                 size_t dig = n / pow;
1622                 n -= dig * pow;
1623                 bc_num_putchar(((uchar) dig) + '0');
1624         }
1625 }
1626
1627 static void bc_num_printHex(size_t n, size_t len, bool rdx) {
1628
1629         BC_UNUSED(len);
1630
1631         assert(len == 1);
1632
1633         if (rdx) bc_num_putchar('.');
1634
1635         bc_num_putchar(bc_num_hex_digits[n]);
1636 }
1637
1638 static void bc_num_printDecimal(const BcNum *restrict n) {
1639
1640         size_t i, j, rdx = n->rdx;
1641         bool zero = true;
1642         size_t buffer[BC_BASE_DIGS];
1643
1644         if (n->neg) bc_num_putchar('-');
1645
1646         for (i = n->len - 1; i < n->len; --i) {
1647
1648                 BcDig n9 = n->num[i];
1649                 size_t temp;
1650                 bool irdx = (i == rdx - 1);
1651
1652                 zero = (zero & !irdx);
1653                 temp = n->scale % BC_BASE_DIGS;
1654                 temp = i || !temp ? 0 : BC_BASE_DIGS - temp;
1655
1656                 memset(buffer, 0, BC_BASE_DIGS * sizeof(size_t));
1657
1658                 for (j = 0; n9 && j < BC_BASE_DIGS; ++j) {
1659                         buffer[j] = n9 % BC_BASE;
1660                         n9 /= BC_BASE;
1661                 }
1662
1663                 for (j = BC_BASE_DIGS - 1; j < BC_BASE_DIGS && j >= temp; --j) {
1664                         bool print_rdx = (irdx & (j == BC_BASE_DIGS - 1));
1665                         zero = (zero && buffer[j] == 0);
1666                         if (!zero) bc_num_printHex(buffer[j], 1, print_rdx);
1667                 }
1668         }
1669 }
1670
1671 #if BC_ENABLE_EXTRA_MATH
1672 static void bc_num_printExponent(const BcNum *restrict n, bool eng) {
1673
1674         bool neg = (n->len <= n->rdx);
1675         BcNum temp, exp;
1676         size_t places, mod;
1677         BcDig digs[BC_NUM_BIGDIG_LOG10];
1678
1679         BC_SIG_LOCK;
1680
1681         bc_num_createCopy(&temp, n);
1682
1683         BC_SETJMP_LOCKED(exit);
1684
1685         BC_SIG_UNLOCK;
1686
1687         if (neg) {
1688
1689                 size_t i, idx = bc_num_nonzeroLen(n) - 1;
1690
1691                 places = 1;
1692
1693                 for (i = BC_BASE_DIGS - 1; i < BC_BASE_DIGS; --i) {
1694                         if (bc_num_pow10[i] > (BcBigDig) n->num[idx]) places += 1;
1695                         else break;
1696                 }
1697
1698                 places += (n->rdx - (idx + 1)) * BC_BASE_DIGS;
1699                 mod = places % 3;
1700
1701                 if (eng && mod != 0) places += 3 - mod;
1702                 bc_num_shiftLeft(&temp, places);
1703         }
1704         else {
1705                 places = bc_num_intDigits(n) - 1;
1706                 mod = places % 3;
1707                 if (eng && mod != 0) places -= 3 - (3 - mod);
1708                 bc_num_shiftRight(&temp, places);
1709         }
1710
1711         bc_num_printDecimal(&temp);
1712         bc_num_putchar('e');
1713
1714         if (!places) {
1715                 bc_num_printHex(0, 1, false);
1716                 goto exit;
1717         }
1718
1719         if (neg) bc_num_putchar('-');
1720
1721         bc_num_setup(&exp, digs, BC_NUM_BIGDIG_LOG10);
1722         bc_num_bigdig2num(&exp, (BcBigDig) places);
1723
1724         bc_num_printDecimal(&exp);
1725
1726 exit:
1727         BC_SIG_MAYLOCK;
1728         bc_num_free(&temp);
1729         BC_LONGJMP_CONT;
1730 }
1731 #endif // BC_ENABLE_EXTRA_MATH
1732
1733 static void bc_num_printFixup(BcNum *restrict n, BcBigDig rem,
1734                               BcBigDig pow, size_t idx)
1735 {
1736         size_t i, len = n->len - idx;
1737         BcBigDig acc;
1738         BcDig *a = n->num + idx;
1739
1740         if (len < 2) return;
1741
1742         for (i = len - 1; i > 0; --i) {
1743
1744                 acc = ((BcBigDig) a[i]) * rem + ((BcBigDig) a[i - 1]);
1745                 a[i - 1] = (BcDig) (acc % pow);
1746                 acc /= pow;
1747                 acc += (BcBigDig) a[i];
1748
1749                 if (acc >= BC_BASE_POW) {
1750
1751                         if (i == len - 1) {
1752                                 len = bc_vm_growSize(len, 1);
1753                                 bc_num_expand(n, bc_vm_growSize(len, idx));
1754                                 a = n->num + idx;
1755                                 a[len - 1] = 0;
1756                         }
1757
1758                         a[i + 1] += acc / BC_BASE_POW;
1759                         acc %= BC_BASE_POW;
1760                 }
1761
1762                 assert(acc < BC_BASE_POW);
1763                 a[i] = (BcDig) acc;
1764         }
1765
1766         n->len = len + idx;
1767 }
1768
1769 static void bc_num_printPrepare(BcNum *restrict n, BcBigDig rem,
1770                                 BcBigDig pow)
1771 {
1772         size_t i;
1773
1774         for (i = 0; i < n->len; ++i) bc_num_printFixup(n, rem, pow, i);
1775
1776         for (i = 0; i < n->len; ++i) {
1777
1778                 assert(pow == ((BcBigDig) ((BcDig) pow)));
1779
1780                 if (n->num[i] >= (BcDig) pow) {
1781
1782                         if (i + 1 == n->len) {
1783                                 n->len = bc_vm_growSize(n->len, 1);
1784                                 bc_num_expand(n, n->len);
1785                                 n->num[i + 1] = 0;
1786                         }
1787
1788                         assert(pow < BC_BASE_POW);
1789                         n->num[i + 1] += n->num[i] / ((BcDig) pow);
1790                         n->num[i] %= (BcDig) pow;
1791                 }
1792         }
1793 }
1794
1795 static void bc_num_printNum(BcNum *restrict n, BcBigDig base,
1796                             size_t len, BcNumDigitOp print)
1797 {
1798         BcVec stack;
1799         BcNum intp, fracp1, fracp2, digit, flen1, flen2, *n1, *n2, *temp;
1800         BcBigDig dig = 0, *ptr, acc, exp;
1801         size_t i, j;
1802         bool radix;
1803         BcDig digit_digs[BC_NUM_BIGDIG_LOG10 + 1];
1804
1805         assert(base > 1);
1806
1807         if (BC_NUM_ZERO(n)) {
1808                 print(0, len, false);
1809                 return;
1810         }
1811
1812         // This function uses an algorithm that Stefan Esser <se@freebsd.org> came
1813         // up with to print the integer part of a number. What it does is convert
1814         // intp into a number of the specified base, but it does it directly,
1815         // instead of just doing a series of divisions and printing the remainders
1816         // in reverse order.
1817         //
1818         // Let me explain in a bit more detail:
1819         //
1820         // The algorithm takes the current least significant digit (after intp has
1821         // been converted to an integer) and the next to least significant digit,
1822         // and it converts the least significant digit into one of the specified
1823         // base, putting any overflow into the next to least significant digit. It
1824         // iterates through the whole number, from least significant to most
1825         // significant, doing this conversion. At the end of that iteration, the
1826         // least significant digit is converted, but the others are not, so it
1827         // iterates again, starting at the next to least significant digit. It keeps
1828         // doing that conversion, skipping one more digit than the last time, until
1829         // all digits have been converted. Then it prints them in reverse order.
1830         //
1831         // That is the gist of the algorithm. It leaves out several things, such as
1832         // the fact that digits are not always converted into the specified base,
1833         // but into something close, basically a power of the specified base. In
1834         // Stefan's words, "You could consider BcDigs to be of base 10^BC_BASE_DIGS
1835         // in the normal case and obase^N for the largest value of N that satisfies
1836         // obase^N <= 10^BC_BASE_DIGS. [This means that] the result is not in base
1837         // "obase", but in base "obase^N", which happens to be printable as a number
1838         // of base "obase" without consideration for neighbouring BcDigs." This fact
1839         // is what necessitates the existence of the loop later in this function.
1840         //
1841         // The conversion happens in bc_num_printPrepare() where the outer loop
1842         // happens and bc_num_printFixup() where the inner loop, or actual
1843         // conversion, happens.
1844
1845         BC_SIG_LOCK;
1846
1847         bc_vec_init(&stack, sizeof(BcBigDig), NULL);
1848         bc_num_init(&fracp1, n->rdx);
1849
1850         bc_num_createCopy(&intp, n);
1851
1852         BC_SETJMP_LOCKED(err);
1853
1854         BC_SIG_UNLOCK;
1855
1856         bc_num_truncate(&intp, intp.scale);
1857
1858         bc_num_sub(n, &intp, &fracp1, 0);
1859
1860         if (base != vm.last_base) {
1861
1862                 vm.last_pow = 1;
1863                 vm.last_exp = 0;
1864
1865                 while (vm.last_pow * base <= BC_BASE_POW) {
1866                         vm.last_pow *= base;
1867                         vm.last_exp += 1;
1868                 }
1869
1870                 vm.last_rem = BC_BASE_POW - vm.last_pow;
1871                 vm.last_base = base;
1872         }
1873
1874         exp = vm.last_exp;
1875
1876         if (vm.last_rem != 0) bc_num_printPrepare(&intp, vm.last_rem, vm.last_pow);
1877
1878         for (i = 0; i < intp.len; ++i) {
1879
1880                 acc = (BcBigDig) intp.num[i];
1881
1882                 for (j = 0; j < exp && (i < intp.len - 1 || acc != 0); ++j)
1883                 {
1884                         if (j != exp - 1) {
1885                                 dig = acc % base;
1886                                 acc /= base;
1887                         }
1888                         else {
1889                                 dig = acc;
1890                                 acc = 0;
1891                         }
1892
1893                         assert(dig < base);
1894
1895                         bc_vec_push(&stack, &dig);
1896                 }
1897
1898                 assert(acc == 0);
1899         }
1900
1901         for (i = 0; i < stack.len; ++i) {
1902                 ptr = bc_vec_item_rev(&stack, i);
1903                 assert(ptr != NULL);
1904                 print(*ptr, len, false);
1905         }
1906
1907         if (!n->scale) goto err;
1908
1909         BC_SIG_LOCK;
1910
1911         BC_UNSETJMP;
1912
1913         bc_num_init(&fracp2, n->rdx);
1914         bc_num_setup(&digit, digit_digs, sizeof(digit_digs) / sizeof(BcDig));
1915         bc_num_init(&flen1, BC_NUM_BIGDIG_LOG10);
1916         bc_num_init(&flen2, BC_NUM_BIGDIG_LOG10);
1917
1918         BC_SETJMP_LOCKED(frac_err);
1919
1920         BC_SIG_UNLOCK;
1921
1922         bc_num_one(&flen1);
1923
1924         radix = true;
1925         n1 = &flen1;
1926         n2 = &flen2;
1927
1928         fracp2.scale = n->scale;
1929         fracp2.rdx = BC_NUM_RDX(fracp2.scale);
1930
1931         while (bc_num_intDigits(n1) < n->scale + 1) {
1932
1933                 bc_num_expand(&fracp2, fracp1.len + 1);
1934                 bc_num_mulArray(&fracp1, base, &fracp2);
1935                 if (fracp2.len < fracp2.rdx) fracp2.len = fracp2.rdx;
1936
1937                 // fracp is guaranteed to be non-negative and small enough.
1938                 bc_num_bigdig2(&fracp2, &dig);
1939
1940                 bc_num_bigdig2num(&digit, dig);
1941                 bc_num_sub(&fracp2, &digit, &fracp1, 0);
1942
1943                 print(dig, len, radix);
1944                 bc_num_mulArray(n1, base, n2);
1945
1946                 radix = false;
1947                 temp = n1;
1948                 n1 = n2;
1949                 n2 = temp;
1950         }
1951
1952 frac_err:
1953         BC_SIG_MAYLOCK;
1954         bc_num_free(&flen2);
1955         bc_num_free(&flen1);
1956         bc_num_free(&fracp2);
1957 err:
1958         BC_SIG_MAYLOCK;
1959         bc_num_free(&fracp1);
1960         bc_num_free(&intp);
1961         bc_vec_free(&stack);
1962         BC_LONGJMP_CONT;
1963 }
1964
1965 static void bc_num_printBase(BcNum *restrict n, BcBigDig base) {
1966
1967         size_t width;
1968         BcNumDigitOp print;
1969         bool neg = n->neg;
1970
1971         if (neg) bc_num_putchar('-');
1972
1973         n->neg = false;
1974
1975         if (base <= BC_NUM_MAX_POSIX_IBASE) {
1976                 width = 1;
1977                 print = bc_num_printHex;
1978         }
1979         else {
1980                 assert(base <= BC_BASE_POW);
1981                 width = bc_num_log10(base - 1);
1982                 print = bc_num_printDigits;
1983         }
1984
1985         bc_num_printNum(n, base, width, print);
1986         n->neg = neg;
1987 }
1988
1989 #if DC_ENABLED
1990 void bc_num_stream(BcNum *restrict n, BcBigDig base) {
1991         bc_num_printNum(n, base, 1, bc_num_printChar);
1992 }
1993 #endif // DC_ENABLED
1994
1995 void bc_num_setup(BcNum *restrict n, BcDig *restrict num, size_t cap) {
1996         assert(n != NULL);
1997         n->num = num;
1998         n->cap = cap;
1999         bc_num_zero(n);
2000 }
2001
2002 void bc_num_init(BcNum *restrict n, size_t req) {
2003
2004         BcDig *num;
2005
2006         BC_SIG_ASSERT_LOCKED;
2007
2008         assert(n != NULL);
2009
2010         req = req >= BC_NUM_DEF_SIZE ? req : BC_NUM_DEF_SIZE;
2011
2012         if (req == BC_NUM_DEF_SIZE && vm.temps.len) {
2013                 BcNum *nptr = bc_vec_top(&vm.temps);
2014                 num = nptr->num;
2015                 bc_vec_pop(&vm.temps);
2016         }
2017         else num = bc_vm_malloc(BC_NUM_SIZE(req));
2018
2019         bc_num_setup(n, num, req);
2020 }
2021
2022 void bc_num_clear(BcNum *restrict n) {
2023         n->num = NULL;
2024         n->cap = 0;
2025 }
2026
2027 void bc_num_free(void *num) {
2028
2029         BcNum *n = (BcNum*) num;
2030
2031         BC_SIG_ASSERT_LOCKED;
2032
2033         assert(n != NULL);
2034
2035         if (n->cap == BC_NUM_DEF_SIZE) bc_vec_push(&vm.temps, n);
2036         else free(n->num);
2037 }
2038
2039 void bc_num_copy(BcNum *d, const BcNum *s) {
2040         assert(d != NULL && s != NULL);
2041         if (d == s) return;
2042         bc_num_expand(d, s->len);
2043         d->len = s->len;
2044         d->neg = s->neg;
2045         d->rdx = s->rdx;
2046         d->scale = s->scale;
2047         memcpy(d->num, s->num, BC_NUM_SIZE(d->len));
2048 }
2049
2050 void bc_num_createCopy(BcNum *d, const BcNum *s) {
2051         BC_SIG_ASSERT_LOCKED;
2052         bc_num_init(d, s->len);
2053         bc_num_copy(d, s);
2054 }
2055
2056 void bc_num_createFromBigdig(BcNum *n, BcBigDig val) {
2057         BC_SIG_ASSERT_LOCKED;
2058         bc_num_init(n, (BC_NUM_BIGDIG_LOG10 - 1) / BC_BASE_DIGS + 1);
2059         bc_num_bigdig2num(n, val);
2060 }
2061
2062 size_t bc_num_scale(const BcNum *restrict n) {
2063         return n->scale;
2064 }
2065
2066 size_t bc_num_len(const BcNum *restrict n) {
2067
2068         size_t len = n->len;
2069
2070         if (BC_NUM_ZERO(n)) return 0;
2071
2072         if (n->rdx == len) {
2073
2074                 size_t zero, scale;
2075
2076                 len = bc_num_nonzeroLen(n);
2077
2078                 scale = n->scale % BC_BASE_DIGS;
2079                 scale = scale ? scale : BC_BASE_DIGS;
2080
2081                 zero = bc_num_zeroDigits(n->num + len - 1);
2082
2083                 len = len * BC_BASE_DIGS - zero - (BC_BASE_DIGS - scale);
2084         }
2085         else len = bc_num_intDigits(n) + n->scale;
2086
2087         return len;
2088 }
2089
2090 void bc_num_parse(BcNum *restrict n, const char *restrict val,
2091                   BcBigDig base, bool letter)
2092 {
2093         assert(n != NULL && val != NULL && base);
2094         assert(base >= BC_NUM_MIN_BASE && base <= vm.maxes[BC_PROG_GLOBALS_IBASE]);
2095         assert(bc_num_strValid(val));
2096
2097         if (letter) {
2098                 BcBigDig dig = bc_num_parseChar(val[0], BC_NUM_MAX_LBASE);
2099                 bc_num_bigdig2num(n, dig);
2100         }
2101         else if (base == BC_BASE) bc_num_parseDecimal(n, val);
2102         else bc_num_parseBase(n, val, base);
2103 }
2104
2105 void bc_num_print(BcNum *restrict n, BcBigDig base, bool newline) {
2106
2107         assert(n != NULL);
2108         assert(BC_ENABLE_EXTRA_MATH || base >= BC_NUM_MIN_BASE);
2109
2110         bc_num_printNewline();
2111
2112         if (BC_NUM_ZERO(n)) bc_num_printHex(0, 1, false);
2113         else if (base == BC_BASE) bc_num_printDecimal(n);
2114 #if BC_ENABLE_EXTRA_MATH
2115         else if (base == 0 || base == 1)
2116                 bc_num_printExponent(n, base != 0);
2117 #endif // BC_ENABLE_EXTRA_MATH
2118         else bc_num_printBase(n, base);
2119
2120         if (newline) bc_num_putchar('\n');
2121 }
2122
2123 void bc_num_bigdig2(const BcNum *restrict n, BcBigDig *result) {
2124
2125         // This function returns no errors because it's guaranteed to succeed if
2126         // its preconditions are met. Those preconditions include both parameters
2127         // being non-NULL, n being non-negative, and n being less than vm.max. If
2128         // all of that is true, then we can just convert without worrying about
2129         // negative errors or overflow. We also don't care about signals because
2130         // this function should execute in only a few iterations, meaning that
2131         // ignoring signals here should be fine.
2132
2133         BcBigDig r = 0;
2134
2135         assert(n != NULL && result != NULL);
2136         assert(!n->neg);
2137         assert(bc_num_cmp(n, &vm.max) < 0);
2138         assert(n->len - n->rdx <= 3);
2139
2140         // There is a small speed win from unrolling the loop here, and since it
2141         // only adds 53 bytes, I decided that it was worth it.
2142         switch (n->len - n->rdx) {
2143                 case 3:
2144                         r = (BcBigDig) n->num[n->rdx + 2];
2145                         // Fallthrough.
2146                 case 2:
2147                         r = r * BC_BASE_POW + (BcBigDig) n->num[n->rdx + 1];
2148                         // Fallthrough.
2149                 case 1:
2150                         r = r * BC_BASE_POW + (BcBigDig) n->num[n->rdx];
2151         }
2152
2153         *result = r;
2154 }
2155
2156 void bc_num_bigdig(const BcNum *restrict n, BcBigDig *result) {
2157
2158         assert(n != NULL && result != NULL);
2159
2160         if (BC_ERR(n->neg)) bc_vm_err(BC_ERROR_MATH_NEGATIVE);
2161         if (BC_ERR(bc_num_cmp(n, &vm.max) >= 0))
2162                 bc_vm_err(BC_ERROR_MATH_OVERFLOW);
2163
2164         bc_num_bigdig2(n, result);
2165 }
2166
2167 void bc_num_bigdig2num(BcNum *restrict n, BcBigDig val) {
2168
2169         BcDig *ptr;
2170         size_t i;
2171
2172         assert(n != NULL);
2173
2174         bc_num_zero(n);
2175
2176         if (!val) return;
2177
2178         bc_num_expand(n, BC_NUM_BIGDIG_LOG10);
2179
2180         for (ptr = n->num, i = 0; val; ++i, val /= BC_BASE_POW)
2181                 ptr[i] = val % BC_BASE_POW;
2182
2183         n->len = i;
2184 }
2185
2186 #if BC_ENABLE_EXTRA_MATH && BC_ENABLE_RAND
2187 void bc_num_rng(const BcNum *restrict n, BcRNG *rng) {
2188
2189         BcNum pow, temp, temp2, intn, frac;
2190         BcRand state1, state2, inc1, inc2;
2191         BcDig pow_num[BC_RAND_NUM_SIZE];
2192
2193         bc_num_setup(&pow, pow_num, sizeof(pow_num) / sizeof(BcDig));
2194
2195         BC_SIG_LOCK;
2196
2197         bc_num_init(&temp, n->len);
2198         bc_num_init(&temp2, n->len);
2199         bc_num_init(&frac, n->rdx);
2200         bc_num_init(&intn, bc_num_int(n));
2201
2202         BC_SETJMP_LOCKED(err);
2203
2204         BC_SIG_UNLOCK;
2205
2206         bc_num_mul(&vm.max, &vm.max, &pow, 0);
2207
2208         memcpy(frac.num, n->num, BC_NUM_SIZE(n->rdx));
2209         frac.len = n->rdx;
2210         frac.rdx = n->rdx;
2211         frac.scale = n->scale;
2212
2213         bc_num_mul(&frac, &pow, &temp, 0);
2214
2215         bc_num_truncate(&temp, temp.scale);
2216         bc_num_copy(&frac, &temp);
2217
2218         memcpy(intn.num, n->num + n->rdx, BC_NUM_SIZE(bc_num_int(n)));
2219         intn.len = bc_num_int(n);
2220
2221         // This assert is here because it has to be true. It is also here to justify
2222         // the use of BC_ERROR_SIGNAL_ONLY() on each of the divmod's and mod's
2223         // below.
2224         assert(BC_NUM_NONZERO(&vm.max));
2225
2226         if (BC_NUM_NONZERO(&frac)) {
2227
2228                 bc_num_divmod(&frac, &vm.max, &temp, &temp2, 0);
2229
2230                 // frac is guaranteed to be smaller than vm.max * vm.max (pow).
2231                 // This means that when dividing frac by vm.max, as above, the
2232                 // quotient and remainder are both guaranteed to be less than vm.max,
2233                 // which means we can use bc_num_bigdig2() here and not worry about
2234                 // overflow.
2235                 bc_num_bigdig2(&temp2, (BcBigDig*) &state1);
2236                 bc_num_bigdig2(&temp, (BcBigDig*) &state2);
2237         }
2238         else state1 = state2 = 0;
2239
2240         if (BC_NUM_NONZERO(&intn)) {
2241
2242                 bc_num_divmod(&intn, &vm.max, &temp, &temp2, 0);
2243
2244                 // Because temp2 is the mod of vm.max, from above, it is guaranteed
2245                 // to be small enough to use bc_num_bigdig2().
2246                 bc_num_bigdig2(&temp2, (BcBigDig*) &inc1);
2247
2248                 if (bc_num_cmp(&temp, &vm.max) >= 0) {
2249                         bc_num_copy(&temp2, &temp);
2250                         bc_num_mod(&temp2, &vm.max, &temp, 0);
2251                 }
2252
2253                 // The if statement above ensures that temp is less than vm.max, which
2254                 // means that we can use bc_num_bigdig2() here.
2255                 bc_num_bigdig2(&temp, (BcBigDig*) &inc2);
2256         }
2257         else inc1 = inc2 = 0;
2258
2259         bc_rand_seed(rng, state1, state2, inc1, inc2);
2260
2261 err:
2262         BC_SIG_MAYLOCK;
2263         bc_num_free(&intn);
2264         bc_num_free(&frac);
2265         bc_num_free(&temp2);
2266         bc_num_free(&temp);
2267         BC_LONGJMP_CONT;
2268 }
2269
2270 void bc_num_createFromRNG(BcNum *restrict n, BcRNG *rng) {
2271
2272         BcRand s1, s2, i1, i2;
2273         BcNum pow, conv, temp1, temp2, temp3;
2274         BcDig pow_num[BC_RAND_NUM_SIZE];
2275         BcDig temp1_num[BC_RAND_NUM_SIZE], temp2_num[BC_RAND_NUM_SIZE];
2276         BcDig conv_num[BC_NUM_BIGDIG_LOG10];
2277
2278         BC_SIG_LOCK;
2279
2280         bc_num_init(&temp3, 2 * BC_RAND_NUM_SIZE);
2281
2282         BC_SETJMP_LOCKED(err);
2283
2284         BC_SIG_UNLOCK;
2285
2286         bc_num_setup(&pow, pow_num, sizeof(pow_num) / sizeof(BcDig));
2287         bc_num_setup(&temp1, temp1_num, sizeof(temp1_num) / sizeof(BcDig));
2288         bc_num_setup(&temp2, temp2_num, sizeof(temp2_num) / sizeof(BcDig));
2289         bc_num_setup(&conv, conv_num, sizeof(conv_num) / sizeof(BcDig));
2290
2291         // This assert is here because it has to be true. It is also here to justify
2292         // the assumption that pow is not zero.
2293         assert(BC_NUM_NONZERO(&vm.max));
2294
2295         bc_num_mul(&vm.max, &vm.max, &pow, 0);
2296
2297         // Because this is true, we can just use BC_ERROR_SIGNAL_ONLY() below when
2298         // dividing by pow.
2299         assert(BC_NUM_NONZERO(&pow));
2300
2301         bc_rand_getRands(rng, &s1, &s2, &i1, &i2);
2302
2303         bc_num_bigdig2num(&conv, (BcBigDig) s2);
2304
2305         bc_num_mul(&conv, &vm.max, &temp1, 0);
2306
2307         bc_num_bigdig2num(&conv, (BcBigDig) s1);
2308
2309         bc_num_add(&conv, &temp1, &temp2, 0);
2310
2311         bc_num_div(&temp2, &pow, &temp3, BC_RAND_STATE_BITS);
2312
2313         bc_num_bigdig2num(&conv, (BcBigDig) i2);
2314
2315         bc_num_mul(&conv, &vm.max, &temp1, 0);
2316
2317         bc_num_bigdig2num(&conv, (BcBigDig) i1);
2318
2319         bc_num_add(&conv, &temp1, &temp2, 0);
2320
2321         bc_num_add(&temp2, &temp3, n, 0);
2322
2323 err:
2324         BC_SIG_MAYLOCK;
2325         bc_num_free(&temp3);
2326         BC_LONGJMP_CONT;
2327 }
2328
2329 void bc_num_irand(const BcNum *restrict a, BcNum *restrict b,
2330                   BcRNG *restrict rng)
2331 {
2332         BcRand r;
2333         BcBigDig modl;
2334         BcNum pow, pow2, cp, cp2, mod, temp1, temp2, rand;
2335         BcNum *p1, *p2, *t1, *t2, *c1, *c2, *tmp;
2336         BcDig rand_num[BC_NUM_BIGDIG_LOG10];
2337         bool carry;
2338         ssize_t cmp;
2339
2340         assert(a != b);
2341
2342         if (BC_ERR(a->neg)) bc_vm_err(BC_ERROR_MATH_NEGATIVE);
2343         if (BC_ERR(a->rdx)) bc_vm_err(BC_ERROR_MATH_NON_INTEGER);
2344         if (BC_NUM_ZERO(a) || BC_NUM_ONE(a)) return;
2345
2346         cmp = bc_num_cmp(a, &vm.max);
2347
2348         if (cmp <= 0) {
2349
2350                 BcRand bits = 0;
2351
2352                 if (cmp < 0) bc_num_bigdig2(a, (BcBigDig*) &bits);
2353
2354                 // This condition means that bits is a power of 2. In that case, we
2355                 // can just grab a full-size int and mask out the unneeded bits.
2356                 // Also, this condition says that 0 is a power of 2, which works for
2357                 // us, since a value of 0 means a == rng->max. The bitmask will mask
2358                 // nothing in that case as well.
2359                 if (!(bits & (bits - 1))) r = bc_rand_int(rng) & (bits - 1);
2360                 else r = bc_rand_bounded(rng, bits);
2361
2362                 // We made sure that r is less than vm.max,
2363                 // so we can use bc_num_bigdig2() here.
2364                 bc_num_bigdig2num(b, r);
2365
2366                 return;
2367         }
2368
2369         // In the case where a is less than rng->max, we have to make sure we have
2370         // an exclusive bound. This ensures that it happens. (See below.)
2371         carry = (cmp < 0);
2372
2373         BC_SIG_LOCK;
2374
2375         bc_num_createCopy(&cp, a);
2376
2377         bc_num_init(&cp2, cp.len);
2378         bc_num_init(&mod, BC_NUM_BIGDIG_LOG10);
2379         bc_num_init(&temp1, BC_NUM_DEF_SIZE);
2380         bc_num_init(&temp2, BC_NUM_DEF_SIZE);
2381         bc_num_init(&pow2, BC_NUM_DEF_SIZE);
2382         bc_num_init(&pow, BC_NUM_DEF_SIZE);
2383         bc_num_one(&pow);
2384         bc_num_setup(&rand, rand_num, sizeof(rand_num) / sizeof(BcDig));
2385
2386         BC_SETJMP_LOCKED(err);
2387
2388         BC_SIG_UNLOCK;
2389
2390         p1 = &pow;
2391         p2 = &pow2;
2392         t1 = &temp1;
2393         t2 = &temp2;
2394         c1 = &cp;
2395         c2 = &cp2;
2396
2397         // This assert is here because it has to be true. It is also here to justify
2398         // the use of BC_ERROR_SIGNAL_ONLY() on each of the divmod's and mod's
2399         // below.
2400         assert(BC_NUM_NONZERO(&vm.max));
2401
2402         while (BC_NUM_NONZERO(c1)) {
2403
2404                 bc_num_divmod(c1, &vm.max, c2, &mod, 0);
2405
2406                 // Because mod is the mod of vm.max, it is guaranteed to be smaller,
2407                 // which means we can use bc_num_bigdig2() here.
2408                 bc_num_bigdig(&mod, &modl);
2409
2410                 if (bc_num_cmp(c1, &vm.max) < 0) {
2411
2412                         // In this case, if there is no carry, then we know we can generate
2413                         // an integer *equal* to modl. Thus, we add one if there is no
2414                         // carry. Otherwise, we add zero, and we are still bounded properly.
2415                         // Since the last portion is guaranteed to be greater than 1, we
2416                         // know modl isn't 0 unless there is no carry.
2417                         modl += !carry;
2418
2419                         if (modl == 1) r = 0;
2420                         else if (!modl) r = bc_rand_int(rng);
2421                         else r = bc_rand_bounded(rng, (BcRand) modl);
2422                 }
2423                 else {
2424                         if (modl) modl -= carry;
2425                         r = bc_rand_int(rng);
2426                         carry = (r >= (BcRand) modl);
2427                 }
2428
2429                 bc_num_bigdig2num(&rand, r);
2430
2431                 bc_num_mul(&rand, p1, p2, 0);
2432                 bc_num_add(p2, t1, t2, 0);
2433
2434                 if (BC_NUM_NONZERO(c2)) {
2435
2436                         bc_num_mul(&vm.max, p1, p2, 0);
2437
2438                         tmp = p1;
2439                         p1 = p2;
2440                         p2 = tmp;
2441
2442                         tmp = c1;
2443                         c1 = c2;
2444                         c2 = tmp;
2445                 }
2446                 else c1 = c2;
2447
2448                 tmp = t1;
2449                 t1 = t2;
2450                 t2 = tmp;
2451         }
2452
2453         bc_num_copy(b, t1);
2454         bc_num_clean(b);
2455
2456 err:
2457         BC_SIG_MAYLOCK;
2458         bc_num_free(&pow);
2459         bc_num_free(&pow2);
2460         bc_num_free(&temp2);
2461         bc_num_free(&temp1);
2462         bc_num_free(&mod);
2463         bc_num_free(&cp2);
2464         bc_num_free(&cp);
2465         BC_LONGJMP_CONT;
2466 }
2467 #endif // BC_ENABLE_EXTRA_MATH && BC_ENABLE_RAND
2468
2469 size_t bc_num_addReq(const BcNum *a, const BcNum *b, size_t scale) {
2470
2471         size_t aint, bint, ardx, brdx;
2472
2473         BC_UNUSED(scale);
2474
2475         ardx = a->rdx;
2476         aint = bc_num_int(a);
2477         assert(aint <= a->len && ardx <= a->len);
2478
2479         brdx = b->rdx;
2480         bint = bc_num_int(b);
2481         assert(bint <= b->len && brdx <= b->len);
2482
2483         ardx = BC_MAX(ardx, brdx);
2484         aint = BC_MAX(aint, bint);
2485
2486         return bc_vm_growSize(bc_vm_growSize(ardx, aint), 1);
2487 }
2488
2489 size_t bc_num_mulReq(const BcNum *a, const BcNum *b, size_t scale) {
2490         size_t max, rdx;
2491         rdx = bc_vm_growSize(a->rdx, b->rdx);
2492         max = BC_NUM_RDX(scale);
2493         max = bc_vm_growSize(BC_MAX(max, rdx), 1);
2494         rdx = bc_vm_growSize(bc_vm_growSize(bc_num_int(a), bc_num_int(b)), max);
2495         return rdx;
2496 }
2497
2498 size_t bc_num_powReq(const BcNum *a, const BcNum *b, size_t scale) {
2499         BC_UNUSED(scale);
2500         return bc_vm_growSize(bc_vm_growSize(a->len, b->len), 1);
2501 }
2502
2503 #if BC_ENABLE_EXTRA_MATH
2504 size_t bc_num_placesReq(const BcNum *a, const BcNum *b, size_t scale) {
2505         BC_UNUSED(scale);
2506         return a->len + b->len - a->rdx - b->rdx;
2507 }
2508 #endif // BC_ENABLE_EXTRA_MATH
2509
2510 void bc_num_add(BcNum *a, BcNum *b, BcNum *c, size_t scale) {
2511         bc_num_binary(a, b, c, false, bc_num_as, bc_num_addReq(a, b, scale));
2512 }
2513
2514 void bc_num_sub(BcNum *a, BcNum *b, BcNum *c, size_t scale) {
2515         bc_num_binary(a, b, c, true, bc_num_as, bc_num_addReq(a, b, scale));
2516 }
2517
2518 void bc_num_mul(BcNum *a, BcNum *b, BcNum *c, size_t scale) {
2519         bc_num_binary(a, b, c, scale, bc_num_m, bc_num_mulReq(a, b, scale));
2520 }
2521
2522 void bc_num_div(BcNum *a, BcNum *b, BcNum *c, size_t scale) {
2523         bc_num_binary(a, b, c, scale, bc_num_d, bc_num_mulReq(a, b, scale));
2524 }
2525
2526 void bc_num_mod(BcNum *a, BcNum *b, BcNum *c, size_t scale) {
2527         bc_num_binary(a, b, c, scale, bc_num_rem, bc_num_mulReq(a, b, scale));
2528 }
2529
2530 void bc_num_pow(BcNum *a, BcNum *b, BcNum *c, size_t scale) {
2531         bc_num_binary(a, b, c, scale, bc_num_p, bc_num_powReq(a, b, scale));
2532 }
2533
2534 #if BC_ENABLE_EXTRA_MATH
2535 void bc_num_places(BcNum *a, BcNum *b, BcNum *c, size_t scale) {
2536         bc_num_binary(a, b, c, scale, bc_num_place, bc_num_placesReq(a, b, scale));
2537 }
2538
2539 void bc_num_lshift(BcNum *a, BcNum *b, BcNum *c, size_t scale) {
2540         bc_num_binary(a, b, c, scale, bc_num_left, bc_num_placesReq(a, b, scale));
2541 }
2542
2543 void bc_num_rshift(BcNum *a, BcNum *b, BcNum *c, size_t scale) {
2544         bc_num_binary(a, b, c, scale, bc_num_right, bc_num_placesReq(a, b, scale));
2545 }
2546 #endif // BC_ENABLE_EXTRA_MATH
2547
2548 void bc_num_sqrt(BcNum *restrict a, BcNum *restrict b, size_t scale) {
2549
2550         BcNum num1, num2, half, f, fprime, *x0, *x1, *temp;
2551         size_t pow, len, rdx, req, digs, digs1, digs2, resscale;
2552         BcDig half_digs[1];
2553
2554         assert(a != NULL && b != NULL && a != b);
2555
2556         if (BC_ERR(a->neg)) bc_vm_err(BC_ERROR_MATH_NEGATIVE);
2557
2558         if (a->scale > scale) scale = a->scale;
2559
2560         len = bc_vm_growSize(bc_num_intDigits(a), 1);
2561         rdx = BC_NUM_RDX(scale);
2562         req = bc_vm_growSize(BC_MAX(rdx, a->rdx), len >> 1);
2563
2564         BC_SIG_LOCK;
2565
2566         bc_num_init(b, bc_vm_growSize(req, 1));
2567
2568         BC_SIG_UNLOCK;
2569
2570         if (BC_NUM_ZERO(a)) {
2571                 bc_num_setToZero(b, scale);
2572                 return;
2573         }
2574         if (BC_NUM_ONE(a)) {
2575                 bc_num_one(b);
2576                 bc_num_extend(b, scale);
2577                 return;
2578         }
2579
2580         rdx = BC_NUM_RDX(scale);
2581         rdx = BC_MAX(rdx, a->rdx);
2582         len = bc_vm_growSize(a->len, rdx);
2583
2584         BC_SIG_LOCK;
2585
2586         bc_num_init(&num1, len);
2587         bc_num_init(&num2, len);
2588         bc_num_setup(&half, half_digs, sizeof(half_digs) / sizeof(BcDig));
2589
2590         bc_num_one(&half);
2591         half.num[0] = BC_BASE_POW / 2;
2592         half.len = 1;
2593         half.rdx = 1;
2594         half.scale = 1;
2595
2596         bc_num_init(&f, len);
2597         bc_num_init(&fprime, len);
2598
2599         BC_SETJMP_LOCKED(err);
2600
2601         BC_SIG_UNLOCK;
2602
2603         x0 = &num1;
2604         x1 = &num2;
2605
2606         bc_num_one(x0);
2607         pow = bc_num_intDigits(a);
2608
2609         if (pow) {
2610
2611                 if (pow & 1) x0->num[0] = 2;
2612                 else x0->num[0] = 6;
2613
2614                 pow -= 2 - (pow & 1);
2615                 bc_num_shiftLeft(x0, pow / 2);
2616         }
2617
2618         x0->scale = x0->rdx = digs = digs1 = digs2 = 0;
2619         resscale = (scale + BC_BASE_DIGS) + 2;
2620
2621         while (bc_num_cmp(x1, x0)) {
2622
2623                 assert(BC_NUM_NONZERO(x0));
2624
2625                 bc_num_div(a, x0, &f, resscale);
2626                 bc_num_add(x0, &f, &fprime, resscale);
2627                 bc_num_mul(&fprime, &half, x1, resscale);
2628
2629                 temp = x0;
2630                 x0 = x1;
2631                 x1 = temp;
2632         }
2633
2634         bc_num_copy(b, x0);
2635         if (b->scale > scale) bc_num_truncate(b, b->scale - scale);
2636
2637         assert(!b->neg || BC_NUM_NONZERO(b));
2638         assert(b->rdx <= b->len || !b->len);
2639         assert(!b->len || b->num[b->len - 1] || b->rdx == b->len);
2640
2641 err:
2642         BC_SIG_MAYLOCK;
2643         bc_num_free(&fprime);
2644         bc_num_free(&f);
2645         bc_num_free(&num2);
2646         bc_num_free(&num1);
2647         BC_LONGJMP_CONT;
2648 }
2649
2650 void bc_num_divmod(BcNum *a, BcNum *b, BcNum *c, BcNum *d, size_t scale) {
2651
2652         BcNum num2, *ptr_a;
2653         bool init = false;
2654         size_t ts, len;
2655
2656         ts = BC_MAX(scale + b->scale, a->scale);
2657         len = bc_num_mulReq(a, b, ts);
2658
2659         assert(a != NULL && b != NULL && c != NULL && d != NULL);
2660         assert(c != d && a != d && b != d && b != c);
2661
2662         if (c == a) {
2663
2664                 memcpy(&num2, c, sizeof(BcNum));
2665                 ptr_a = &num2;
2666
2667                 BC_SIG_LOCK;
2668
2669                 bc_num_init(c, len);
2670
2671                 init = true;
2672
2673                 BC_SETJMP_LOCKED(err);
2674
2675                 BC_SIG_UNLOCK;
2676         }
2677         else {
2678                 ptr_a = a;
2679                 bc_num_expand(c, len);
2680         }
2681
2682         if (BC_NUM_NONZERO(a) && !a->rdx && !b->rdx && b->len == 1 && !scale) {
2683
2684                 BcBigDig rem;
2685
2686                 bc_num_divArray(ptr_a, (BcBigDig) b->num[0], c, &rem);
2687
2688                 assert(rem < BC_BASE_POW);
2689
2690                 d->num[0] = (BcDig) rem;
2691                 d->len = (rem != 0);
2692         }
2693         else bc_num_r(ptr_a, b, c, d, scale, ts);
2694
2695         assert(!c->neg || BC_NUM_NONZERO(c));
2696         assert(c->rdx <= c->len || !c->len);
2697         assert(!c->len || c->num[c->len - 1] || c->rdx == c->len);
2698         assert(!d->neg || BC_NUM_NONZERO(d));
2699         assert(d->rdx <= d->len || !d->len);
2700         assert(!d->len || d->num[d->len - 1] || d->rdx == d->len);
2701
2702 err:
2703         if (init) {
2704                 BC_SIG_MAYLOCK;
2705                 bc_num_free(&num2);
2706                 BC_LONGJMP_CONT;
2707         }
2708 }
2709
2710 #if DC_ENABLED
2711 void bc_num_modexp(BcNum *a, BcNum *b, BcNum *c, BcNum *restrict d) {
2712
2713         BcNum base, exp, two, temp;
2714         BcDig two_digs[2];
2715
2716         assert(a != NULL && b != NULL && c != NULL && d != NULL);
2717         assert(a != d && b != d && c != d);
2718
2719         if (BC_ERR(BC_NUM_ZERO(c))) bc_vm_err(BC_ERROR_MATH_DIVIDE_BY_ZERO);
2720         if (BC_ERR(b->neg)) bc_vm_err(BC_ERROR_MATH_NEGATIVE);
2721         if (BC_ERR(a->rdx || b->rdx || c->rdx))
2722                 bc_vm_err(BC_ERROR_MATH_NON_INTEGER);
2723
2724         bc_num_expand(d, c->len);
2725
2726         BC_SIG_LOCK;
2727
2728         bc_num_init(&base, c->len);
2729         bc_num_setup(&two, two_digs, sizeof(two_digs) / sizeof(BcDig));
2730         bc_num_init(&temp, b->len + 1);
2731         bc_num_createCopy(&exp, b);
2732
2733         BC_SETJMP_LOCKED(err);
2734
2735         BC_SIG_UNLOCK;
2736
2737         bc_num_one(&two);
2738         two.num[0] = 2;
2739         bc_num_one(d);
2740
2741         // We already checked for 0.
2742         bc_num_rem(a, c, &base, 0);
2743
2744         while (BC_NUM_NONZERO(&exp)) {
2745
2746                 // Num two cannot be 0, so no errors.
2747                 bc_num_divmod(&exp, &two, &exp, &temp, 0);
2748
2749                 if (BC_NUM_ONE(&temp) && !temp.neg) {
2750
2751                         bc_num_mul(d, &base, &temp, 0);
2752
2753                         // We already checked for 0.
2754                         bc_num_rem(&temp, c, d, 0);
2755                 }
2756
2757                 bc_num_mul(&base, &base, &temp, 0);
2758
2759                 // We already checked for 0.
2760                 bc_num_rem(&temp, c, &base, 0);
2761         }
2762
2763 err:
2764         BC_SIG_MAYLOCK;
2765         bc_num_free(&exp);
2766         bc_num_free(&temp);
2767         bc_num_free(&base);
2768         BC_LONGJMP_CONT;
2769         assert(!d->neg || d->len);
2770         assert(!d->len || d->num[d->len - 1] || d->rdx == d->len);
2771 }
2772 #endif // DC_ENABLED
2773
2774 #if BC_DEBUG_CODE
2775 void bc_num_printDebug(const BcNum *n, const char *name, bool emptyline) {
2776         bc_file_puts(&vm.fout, name);
2777         bc_file_puts(&vm.fout, ": ");
2778         bc_num_printDecimal(n);
2779         bc_file_putchar(&vm.fout, '\n');
2780         if (emptyline) bc_file_putchar(&vm.fout, '\n');
2781         vm.nchars = 0;
2782 }
2783
2784 void bc_num_printDigs(const BcDig *n, size_t len, bool emptyline) {
2785
2786         size_t i;
2787
2788         for (i = len - 1; i < len; --i)
2789                 bc_file_printf(&vm.fout, " %lu", (unsigned long) n[i]);
2790
2791         bc_file_putchar(&vm.fout, '\n');
2792         if (emptyline) bc_file_putchar(&vm.fout, '\n');
2793         vm.nchars = 0;
2794 }
2795
2796 void bc_num_printWithDigs(const BcNum *n, const char *name, bool emptyline) {
2797         bc_file_puts(&vm.fout, name);
2798         bc_file_printf(&vm.fout, " len: %zu, rdx: %zu, scale: %zu\n",
2799                        name, n->len, n->rdx, n->scale);
2800         bc_num_printDigs(n->num, n->len, emptyline);
2801 }
2802
2803 void bc_num_dump(const char *varname, const BcNum *n) {
2804
2805         ulong i, scale = n->scale;
2806
2807         bc_file_printf(&vm.ferr, "\n%s = %s", varname,
2808                        n->len ? (n->neg ? "-" : "+") : "0 ");
2809
2810         for (i = n->len - 1; i < n->len; --i) {
2811
2812                 if (i + 1 == n->rdx) bc_file_puts(&vm.ferr, ". ");
2813
2814                 if (scale / BC_BASE_DIGS != n->rdx - i - 1)
2815                         bc_file_printf(&vm.ferr, "%lu ", (unsigned long) n->num[i]);
2816                 else {
2817
2818                         int mod = scale % BC_BASE_DIGS;
2819                         int d = BC_BASE_DIGS - mod;
2820                         BcDig div;
2821
2822                         if (mod != 0) {
2823                                 div = n->num[i] / ((BcDig) bc_num_pow10[(ulong) d]);
2824                                 bc_file_printf(&vm.ferr, "%lu", (unsigned long) div);
2825                         }
2826
2827                         div = n->num[i] % ((BcDig) bc_num_pow10[(ulong) d]);
2828                         bc_file_printf(&vm.ferr, " ' %lu ", (unsigned long) div);
2829                 }
2830         }
2831
2832         bc_file_printf(&vm.ferr, "(%zu | %zu.%zu / %zu) %lu\n",
2833                        n->scale, n->len, n->rdx, n->cap,
2834                        (unsigned long) (void*) n->num);
2835 }
2836 #endif // BC_DEBUG_CODE