1/* 2 * ***************************************************************************** 3 * 4 * SPDX-License-Identifier: BSD-2-Clause 5 * 6 * Copyright (c) 2018-2021 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 48static void bc_num_m(BcNum *a, BcNum *b, BcNum *restrict c, size_t scale); 49 50static inline ssize_t bc_num_neg(size_t n, bool neg) { 51 return (((ssize_t) n) ^ -((ssize_t) neg)) + neg; 52} 53 54ssize_t bc_num_cmpZero(const BcNum *n) { 55 return bc_num_neg((n)->len != 0, BC_NUM_NEG(n)); 56} 57 58static inline size_t bc_num_int(const BcNum *n) { 59 return n->len ? n->len - BC_NUM_RDX_VAL(n) : 0; 60} 61 62static 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 79static 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 85void bc_num_zero(BcNum *restrict n) { 86 bc_num_setToZero(n, 0); 87} 88 89void bc_num_one(BcNum *restrict n) { 90 bc_num_zero(n); 91 n->len = 1; 92 n->num[0] = 1; 93} 94 95static 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 106static 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 113static 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 119static 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 125static 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 133static 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 149static 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 164static 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 175static 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 186static 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 216static 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 240static 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 249ssize_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 301void 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 332void 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 358static 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 368static 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 395static 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 409static void bc_num_unshiftZero(BcNum *restrict n, size_t places_rdx) { 410 n->len += places_rdx; 411 n->num -= places_rdx; 412} 413 414static 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 436static 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 495void 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 546static 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 560static 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 569static 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 714static 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 762static 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 770static 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 870err: 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 880static 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 977err: 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 986static 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 993static 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 1007static 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 1020static 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 1125err: 1126 BC_SIG_MAYLOCK; 1127 bc_num_free(&cpb); 1128 BC_LONGJMP_CONT; 1129} 1130 1131static 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 1200err: 1201 BC_SIG_MAYLOCK; 1202 bc_num_free(&cpb); 1203 bc_num_free(&cpa); 1204 BC_LONGJMP_CONT; 1205} 1206 1207static 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 1244err: 1245 BC_SIG_MAYLOCK; 1246 bc_num_free(&temp); 1247 BC_LONGJMP_CONT; 1248} 1249 1250static 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 1268err: 1269 BC_SIG_MAYLOCK; 1270 bc_num_free(&c1); 1271 BC_LONGJMP_CONT; 1272} 1273 1274static 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(©, 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(©, ©, ©, powrdx); 1320 } 1321 1322 bc_num_copy(c, ©); 1323 resrdx = powrdx; 1324 1325 while (pow >>= 1) { 1326 1327 powrdx <<= 1; 1328 assert(BC_NUM_RDX_VALID_NP(copy)); 1329 bc_num_mul(©, ©, ©, 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, ©, 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 1347err: 1348 BC_SIG_MAYLOCK; 1349 bc_num_free(©); 1350 BC_LONGJMP_CONT; 1351} 1352 1353#if BC_ENABLE_EXTRA_MATH 1354static 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 1366static 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 1377static 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 1391static void bc_num_binary(BcNum *a, BcNum *b, BcNum *c, size_t scale, 1392 BcNumBinOp 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 1447err: 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 1456bool 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 1482static 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 1493static 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 1561static 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 && !val[i]) goto int_err; 1592 1593 assert(val[i] == '.'); 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 1643err: 1644 BC_SIG_MAYLOCK; 1645 bc_num_free(&result2); 1646 bc_num_free(&result1); 1647 bc_num_free(&mult2); 1648int_err: 1649 BC_SIG_MAYLOCK; 1650 bc_num_free(&mult1); 1651 bc_num_free(&temp); 1652 BC_LONGJMP_CONT; 1653} 1654 1655static inline void bc_num_printNewline(void) { 1656#if !BC_ENABLE_LIBRARY 1657 if (vm.nchars >= vm.line_len - 1) { 1658 bc_vm_putchar('\\', bc_flush_none); 1659 bc_vm_putchar('\n', bc_flush_err); 1660 } 1661#endif // !BC_ENABLE_LIBRARY 1662} 1663 1664static void bc_num_putchar(int c) { 1665 if (c != '\n') bc_num_printNewline(); 1666 bc_vm_putchar(c, bc_flush_save); 1667} 1668 1669#if DC_ENABLED && !BC_ENABLE_LIBRARY 1670static 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, bc_flush_save); 1675} 1676#endif // DC_ENABLED && !BC_ENABLE_LIBRARY 1677 1678static 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 1693static 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 1704static 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 1738static 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 1792exit: 1793 BC_SIG_MAYLOCK; 1794 bc_num_free(&temp); 1795 BC_LONGJMP_CONT; 1796} 1797#endif // BC_ENABLE_EXTRA_MATH 1798 1799static 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 1835static 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 1861static 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 2023frac_err: 2024 BC_SIG_MAYLOCK; 2025 bc_num_free(&flen2); 2026 bc_num_free(&flen1); 2027 bc_num_free(&fracp2); 2028err: 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 2036static 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 2061void 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 2066void 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 2073void 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 2093void bc_num_clear(BcNum *restrict n) { 2094 n->num = NULL; 2095 n->cap = 0; 2096} 2097 2098void 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 2110void 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 2121void 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 2127void 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 2133size_t bc_num_scale(const BcNum *restrict n) { 2134 return n->scale; 2135} 2136 2137size_t bc_num_len(const BcNum *restrict n) { 2138 2139 size_t len = n->len; 2140 2141 if (BC_NUM_ZERO(n)) return n->scale ? n->scale : 1; 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 2161void 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 2177void 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 2194void 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 2237void 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 2248void 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 2268void 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 2342err: 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 2351void 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 2406err: 2407 BC_SIG_MAYLOCK; 2408 bc_num_free(&temp3); 2409 BC_LONGJMP_CONT; 2410} 2411 2412void 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 2546err: 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 2559size_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 2579size_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 2588size_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 2597size_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 2603size_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 2609void 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 2615void 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 2621void 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 2627void 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 2633void 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 2639void 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 2646void 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 2652void 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 2658void 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 2665void 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, 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 = 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 2767err: 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 2776void 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 2831err: 2832 if (init) { 2833 BC_SIG_MAYLOCK; 2834 bc_num_free(&num2); 2835 BC_LONGJMP_CONT; 2836 } 2837} 2838 2839#if DC_ENABLED 2840void 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 2897err: 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 2910void bc_num_printDebug(const BcNum *n, const char *name, bool emptyline) { 2911 bc_file_puts(&vm.fout, bc_flush_none, name); 2912 bc_file_puts(&vm.fout, bc_flush_none, ": "); 2913 bc_num_printDecimal(n); 2914 bc_file_putchar(&vm.fout, bc_flush_err, '\n'); 2915 if (emptyline) bc_file_putchar(&vm.fout, bc_flush_err, '\n'); 2916 vm.nchars = 0; 2917} 2918 2919void 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, bc_flush_err, '\n'); 2927 if (emptyline) bc_file_putchar(&vm.fout, bc_flush_err, '\n'); 2928 vm.nchars = 0; 2929} 2930 2931void bc_num_printWithDigs(const BcNum *n, const char *name, bool emptyline) { 2932 bc_file_puts(&vm.fout, bc_flush_none, 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 2938void 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)) 2948 bc_file_puts(&vm.ferr, bc_flush_none, ". "); 2949 2950 if (scale / BC_BASE_DIGS != BC_NUM_RDX_VAL(n) - i - 1) 2951 bc_file_printf(&vm.ferr, "%lu ", (unsigned long) n->num[i]); 2952 else { 2953 2954 int mod = scale % BC_BASE_DIGS; 2955 int d = BC_BASE_DIGS - mod; 2956 BcDig div; 2957 2958 if (mod != 0) { 2959 div = n->num[i] / ((BcDig) bc_num_pow10[(ulong) d]); 2960 bc_file_printf(&vm.ferr, "%lu", (unsigned long) div); 2961 } 2962 2963 div = n->num[i] % ((BcDig) bc_num_pow10[(ulong) d]); 2964 bc_file_printf(&vm.ferr, " ' %lu ", (unsigned long) div); 2965 } 2966 } 2967 2968 bc_file_printf(&vm.ferr, "(%zu | %zu.%zu / %zu) %lu\n", 2969 n->scale, n->len, BC_NUM_RDX_VAL(n), n->cap, 2970 (unsigned long) (void*) n->num); 2971 2972 bc_file_flush(&vm.ferr, bc_flush_err); 2973} 2974#endif // BC_DEBUG_CODE 2975