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