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