1 /* number.c: Implements arbitrary precision numbers. */ 2 3 /* This file is part of GNU bc. 4 Copyright (C) 1991, 1992, 1993, 1994, 1997, 1998 Free Software Foundation, Inc. 5 6 This program is free software; you can redistribute it and/or modify 7 it under the terms of the GNU General Public License as published by 8 the Free Software Foundation; either version 2 of the License , or 9 (at your option) any later version. 10 11 This program is distributed in the hope that it will be useful, 12 but WITHOUT ANY WARRANTY; without even the implied warranty of 13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 14 GNU General Public License for more details. 15 16 You should have received a copy of the GNU General Public License 17 along with this program; see the file COPYING. If not, write to 18 the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. 19 20 You may contact the author by: 21 e-mail: phil@cs.wwu.edu 22 us-mail: Philip A. Nelson 23 Computer Science Department, 9062 24 Western Washington University 25 Bellingham, WA 98226-9062 26 27 *************************************************************************/ 28 29 #include "bcdefs.h" 30 #include "proto.h" 31 #include "global.h" 32 33 /* Storage used for special numbers. */ 34 bc_num _zero_; 35 bc_num _one_; 36 bc_num _two_; 37 38 39 /* "Frees" a bc_num NUM. Actually decreases reference count and only 40 frees the storage if reference count is zero. */ 41 42 void 43 free_num (num) 44 bc_num *num; 45 { 46 if (*num == NULL) return; 47 (*num)->n_refs--; 48 if ((*num)->n_refs == 0) free(*num); 49 *num = NULL; 50 } 51 52 53 /* new_num allocates a number and sets fields to known values. */ 54 55 bc_num 56 new_num (length, scale) 57 int length, scale; 58 { 59 bc_num temp; 60 61 temp = (bc_num) malloc (sizeof(bc_struct)+length+scale); 62 if (temp == NULL) out_of_memory (); 63 temp->n_sign = PLUS; 64 temp->n_len = length; 65 temp->n_scale = scale; 66 temp->n_refs = 1; 67 temp->n_value[0] = 0; 68 return temp; 69 } 70 71 72 /* Intitialize the number package! */ 73 74 void 75 init_numbers () 76 { 77 _zero_ = new_num (1,0); 78 _one_ = new_num (1,0); 79 _one_->n_value[0] = 1; 80 _two_ = new_num (1,0); 81 _two_->n_value[0] = 2; 82 } 83 84 85 /* Make a copy of a number! Just increments the reference count! */ 86 87 bc_num 88 copy_num (num) 89 bc_num num; 90 { 91 num->n_refs++; 92 return num; 93 } 94 95 96 /* Initialize a number NUM by making it a copy of zero. */ 97 98 void 99 init_num (num) 100 bc_num *num; 101 { 102 *num = copy_num (_zero_); 103 } 104 105 106 /* Convert an integer VAL to a bc number NUM. */ 107 108 void 109 int2num (num, val) 110 bc_num *num; 111 int val; 112 { 113 char buffer[30]; 114 char *bptr, *vptr; 115 int ix = 1; 116 char neg = 0; 117 118 /* Sign. */ 119 if (val < 0) 120 { 121 neg = 1; 122 val = -val; 123 } 124 125 /* Get things going. */ 126 bptr = buffer; 127 *bptr++ = val % BASE; 128 val = val / BASE; 129 130 /* Extract remaining digits. */ 131 while (val != 0) 132 { 133 *bptr++ = val % BASE; 134 val = val / BASE; 135 ix++; /* Count the digits. */ 136 } 137 138 /* Make the number. */ 139 free_num (num); 140 *num = new_num (ix, 0); 141 if (neg) (*num)->n_sign = MINUS; 142 143 /* Assign the digits. */ 144 vptr = (*num)->n_value; 145 while (ix-- > 0) 146 *vptr++ = *--bptr; 147 } 148 149 150 /* Convert a number NUM to a long. The function returns only the integer 151 part of the number. For numbers that are too large to represent as 152 a long, this function returns a zero. This can be detected by checking 153 the NUM for zero after having a zero returned. */ 154 155 long 156 num2long (num) 157 bc_num num; 158 { 159 long val; 160 char *nptr; 161 int index; 162 163 /* Extract the int value, ignore the fraction. */ 164 val = 0; 165 nptr = num->n_value; 166 for (index=num->n_len; (index>0) && (val<=(LONG_MAX/BASE)); index--) 167 val = val*BASE + *nptr++; 168 169 /* Check for overflow. If overflow, return zero. */ 170 if (index>0) val = 0; 171 if (val < 0) val = 0; 172 173 /* Return the value. */ 174 if (num->n_sign == PLUS) 175 return (val); 176 else 177 return (-val); 178 } 179 180 181 /* The following are some math routines for numbers. */ 182 _PROTOTYPE(static int _do_compare, (bc_num n1, bc_num n2, int use_sign, 183 int ignore_last)); 184 _PROTOTYPE(static void _rm_leading_zeros, (bc_num num)); 185 _PROTOTYPE(static bc_num _do_add, (bc_num n1, bc_num n2, int scale_min)); 186 _PROTOTYPE(static bc_num _do_sub, (bc_num n1, bc_num n2, int scale_min)); 187 _PROTOTYPE(static void _one_mult, (unsigned char *num, int size, int digit, 188 unsigned char *result)); 189 190 191 192 /* Compare two bc numbers. Return value is 0 if equal, -1 if N1 is less 193 than N2 and +1 if N1 is greater than N2. If USE_SIGN is false, just 194 compare the magnitudes. */ 195 196 static int 197 _do_compare (n1, n2, use_sign, ignore_last) 198 bc_num n1, n2; 199 int use_sign; 200 int ignore_last; 201 { 202 char *n1ptr, *n2ptr; 203 int count; 204 205 /* First, compare signs. */ 206 if (use_sign && n1->n_sign != n2->n_sign) 207 { 208 if (n1->n_sign == PLUS) 209 return (1); /* Positive N1 > Negative N2 */ 210 else 211 return (-1); /* Negative N1 < Positive N1 */ 212 } 213 214 /* Now compare the magnitude. */ 215 if (n1->n_len != n2->n_len) 216 { 217 if (n1->n_len > n2->n_len) 218 { 219 /* Magnitude of n1 > n2. */ 220 if (!use_sign || n1->n_sign == PLUS) 221 return (1); 222 else 223 return (-1); 224 } 225 else 226 { 227 /* Magnitude of n1 < n2. */ 228 if (!use_sign || n1->n_sign == PLUS) 229 return (-1); 230 else 231 return (1); 232 } 233 } 234 235 /* If we get here, they have the same number of integer digits. 236 check the integer part and the equal length part of the fraction. */ 237 count = n1->n_len + MIN (n1->n_scale, n2->n_scale); 238 n1ptr = n1->n_value; 239 n2ptr = n2->n_value; 240 241 while ((count > 0) && (*n1ptr == *n2ptr)) 242 { 243 n1ptr++; 244 n2ptr++; 245 count--; 246 } 247 if (ignore_last && count == 1 && n1->n_scale == n2->n_scale) 248 return (0); 249 if (count != 0) 250 { 251 if (*n1ptr > *n2ptr) 252 { 253 /* Magnitude of n1 > n2. */ 254 if (!use_sign || n1->n_sign == PLUS) 255 return (1); 256 else 257 return (-1); 258 } 259 else 260 { 261 /* Magnitude of n1 < n2. */ 262 if (!use_sign || n1->n_sign == PLUS) 263 return (-1); 264 else 265 return (1); 266 } 267 } 268 269 /* They are equal up to the last part of the equal part of the fraction. */ 270 if (n1->n_scale != n2->n_scale) 271 if (n1->n_scale > n2->n_scale) 272 { 273 for (count = n1->n_scale-n2->n_scale; count>0; count--) 274 if (*n1ptr++ != 0) 275 { 276 /* Magnitude of n1 > n2. */ 277 if (!use_sign || n1->n_sign == PLUS) 278 return (1); 279 else 280 return (-1); 281 } 282 } 283 else 284 { 285 for (count = n2->n_scale-n1->n_scale; count>0; count--) 286 if (*n2ptr++ != 0) 287 { 288 /* Magnitude of n1 < n2. */ 289 if (!use_sign || n1->n_sign == PLUS) 290 return (-1); 291 else 292 return (1); 293 } 294 } 295 296 /* They must be equal! */ 297 return (0); 298 } 299 300 301 /* This is the "user callable" routine to compare numbers N1 and N2. */ 302 303 int 304 bc_compare (n1, n2) 305 bc_num n1, n2; 306 { 307 return _do_compare (n1, n2, TRUE, FALSE); 308 } 309 310 311 /* In some places we need to check if the number NUM is zero. */ 312 313 char 314 is_zero (num) 315 bc_num num; 316 { 317 int count; 318 char *nptr; 319 320 /* Quick check. */ 321 if (num == _zero_) return TRUE; 322 323 /* Initialize */ 324 count = num->n_len + num->n_scale; 325 nptr = num->n_value; 326 327 /* The check */ 328 while ((count > 0) && (*nptr++ == 0)) count--; 329 330 if (count != 0) 331 return FALSE; 332 else 333 return TRUE; 334 } 335 336 337 /* In some places we need to check if the number is negative. */ 338 339 char 340 is_neg (num) 341 bc_num num; 342 { 343 return num->n_sign == MINUS; 344 } 345 346 347 /* For many things, we may have leading zeros in a number NUM. 348 _rm_leading_zeros just moves the data to the correct 349 place and adjusts the length. */ 350 351 static void 352 _rm_leading_zeros (num) 353 bc_num num; 354 { 355 int bytes; 356 char *dst, *src; 357 358 /* Do a quick check to see if we need to do it. */ 359 if (*num->n_value != 0) return; 360 361 /* The first "digit" is 0, find the first non-zero digit in the second 362 or greater "digit" to the left of the decimal place. */ 363 bytes = num->n_len; 364 src = num->n_value; 365 while (bytes > 1 && *src == 0) src++, bytes--; 366 num->n_len = bytes; 367 bytes += num->n_scale; 368 dst = num->n_value; 369 while (bytes-- > 0) *dst++ = *src++; 370 371 } 372 373 374 /* Perform addition: N1 is added to N2 and the value is 375 returned. The signs of N1 and N2 are ignored. 376 SCALE_MIN is to set the minimum scale of the result. */ 377 378 static bc_num 379 _do_add (n1, n2, scale_min) 380 bc_num n1, n2; 381 int scale_min; 382 { 383 bc_num sum; 384 int sum_scale, sum_digits; 385 char *n1ptr, *n2ptr, *sumptr; 386 int carry, n1bytes, n2bytes; 387 int count; 388 389 /* Prepare sum. */ 390 sum_scale = MAX (n1->n_scale, n2->n_scale); 391 sum_digits = MAX (n1->n_len, n2->n_len) + 1; 392 sum = new_num (sum_digits, MAX(sum_scale, scale_min)); 393 394 /* Zero extra digits made by scale_min. */ 395 if (scale_min > sum_scale) 396 { 397 sumptr = (char *) (sum->n_value + sum_scale + sum_digits); 398 for (count = scale_min - sum_scale; count > 0; count--) 399 *sumptr++ = 0; 400 } 401 402 /* Start with the fraction part. Initialize the pointers. */ 403 n1bytes = n1->n_scale; 404 n2bytes = n2->n_scale; 405 n1ptr = (char *) (n1->n_value + n1->n_len + n1bytes - 1); 406 n2ptr = (char *) (n2->n_value + n2->n_len + n2bytes - 1); 407 sumptr = (char *) (sum->n_value + sum_scale + sum_digits - 1); 408 409 /* Add the fraction part. First copy the longer fraction.*/ 410 if (n1bytes != n2bytes) 411 { 412 if (n1bytes > n2bytes) 413 while (n1bytes>n2bytes) 414 { *sumptr-- = *n1ptr--; n1bytes--;} 415 else 416 while (n2bytes>n1bytes) 417 { *sumptr-- = *n2ptr--; n2bytes--;} 418 } 419 420 /* Now add the remaining fraction part and equal size integer parts. */ 421 n1bytes += n1->n_len; 422 n2bytes += n2->n_len; 423 carry = 0; 424 while ((n1bytes > 0) && (n2bytes > 0)) 425 { 426 *sumptr = *n1ptr-- + *n2ptr-- + carry; 427 if (*sumptr > (BASE-1)) 428 { 429 carry = 1; 430 *sumptr -= BASE; 431 } 432 else 433 carry = 0; 434 sumptr--; 435 n1bytes--; 436 n2bytes--; 437 } 438 439 /* Now add carry the longer integer part. */ 440 if (n1bytes == 0) 441 { n1bytes = n2bytes; n1ptr = n2ptr; } 442 while (n1bytes-- > 0) 443 { 444 *sumptr = *n1ptr-- + carry; 445 if (*sumptr > (BASE-1)) 446 { 447 carry = 1; 448 *sumptr -= BASE; 449 } 450 else 451 carry = 0; 452 sumptr--; 453 } 454 455 /* Set final carry. */ 456 if (carry == 1) 457 *sumptr += 1; 458 459 /* Adjust sum and return. */ 460 _rm_leading_zeros (sum); 461 return sum; 462 } 463 464 465 /* Perform subtraction: N2 is subtracted from N1 and the value is 466 returned. The signs of N1 and N2 are ignored. Also, N1 is 467 assumed to be larger than N2. SCALE_MIN is the minimum scale 468 of the result. */ 469 470 static bc_num 471 _do_sub (n1, n2, scale_min) 472 bc_num n1, n2; 473 int scale_min; 474 { 475 bc_num diff; 476 int diff_scale, diff_len; 477 int min_scale, min_len; 478 char *n1ptr, *n2ptr, *diffptr; 479 int borrow, count, val; 480 481 /* Allocate temporary storage. */ 482 diff_len = MAX (n1->n_len, n2->n_len); 483 diff_scale = MAX (n1->n_scale, n2->n_scale); 484 min_len = MIN (n1->n_len, n2->n_len); 485 min_scale = MIN (n1->n_scale, n2->n_scale); 486 diff = new_num (diff_len, MAX(diff_scale, scale_min)); 487 488 /* Zero extra digits made by scale_min. */ 489 if (scale_min > diff_scale) 490 { 491 diffptr = (char *) (diff->n_value + diff_len + diff_scale); 492 for (count = scale_min - diff_scale; count > 0; count--) 493 *diffptr++ = 0; 494 } 495 496 /* Initialize the subtract. */ 497 n1ptr = (char *) (n1->n_value + n1->n_len + n1->n_scale -1); 498 n2ptr = (char *) (n2->n_value + n2->n_len + n2->n_scale -1); 499 diffptr = (char *) (diff->n_value + diff_len + diff_scale -1); 500 501 /* Subtract the numbers. */ 502 borrow = 0; 503 504 /* Take care of the longer scaled number. */ 505 if (n1->n_scale != min_scale) 506 { 507 /* n1 has the longer scale */ 508 for (count = n1->n_scale - min_scale; count > 0; count--) 509 *diffptr-- = *n1ptr--; 510 } 511 else 512 { 513 /* n2 has the longer scale */ 514 for (count = n2->n_scale - min_scale; count > 0; count--) 515 { 516 val = - *n2ptr-- - borrow; 517 if (val < 0) 518 { 519 val += BASE; 520 borrow = 1; 521 } 522 else 523 borrow = 0; 524 *diffptr-- = val; 525 } 526 } 527 528 /* Now do the equal length scale and integer parts. */ 529 530 for (count = 0; count < min_len + min_scale; count++) 531 { 532 val = *n1ptr-- - *n2ptr-- - borrow; 533 if (val < 0) 534 { 535 val += BASE; 536 borrow = 1; 537 } 538 else 539 borrow = 0; 540 *diffptr-- = val; 541 } 542 543 /* If n1 has more digits then n2, we now do that subtract. */ 544 if (diff_len != min_len) 545 { 546 for (count = diff_len - min_len; count > 0; count--) 547 { 548 val = *n1ptr-- - borrow; 549 if (val < 0) 550 { 551 val += BASE; 552 borrow = 1; 553 } 554 else 555 borrow = 0; 556 *diffptr-- = val; 557 } 558 } 559 560 /* Clean up and return. */ 561 _rm_leading_zeros (diff); 562 return diff; 563 } 564 565 566 /* Here is the full add routine that takes care of negative numbers. 567 N1 is added to N2 and the result placed into RESULT. SCALE_MIN 568 is the minimum scale for the result. */ 569 570 void 571 bc_add (n1, n2, result, scale_min) 572 bc_num n1, n2, *result; 573 int scale_min; 574 { 575 bc_num sum; 576 int cmp_res; 577 int res_scale; 578 579 if (n1->n_sign == n2->n_sign) 580 { 581 sum = _do_add (n1, n2, scale_min); 582 sum->n_sign = n1->n_sign; 583 } 584 else 585 { 586 /* subtraction must be done. */ 587 cmp_res = _do_compare (n1, n2, FALSE, FALSE); /* Compare magnitudes. */ 588 switch (cmp_res) 589 { 590 case -1: 591 /* n1 is less than n2, subtract n1 from n2. */ 592 sum = _do_sub (n2, n1, scale_min); 593 sum->n_sign = n2->n_sign; 594 break; 595 case 0: 596 /* They are equal! return zero with the correct scale! */ 597 res_scale = MAX (scale_min, MAX(n1->n_scale, n2->n_scale)); 598 sum = new_num (1, res_scale); 599 memset (sum->n_value, 0, res_scale+1); 600 break; 601 case 1: 602 /* n2 is less than n1, subtract n2 from n1. */ 603 sum = _do_sub (n1, n2, scale_min); 604 sum->n_sign = n1->n_sign; 605 } 606 } 607 608 /* Clean up and return. */ 609 free_num (result); 610 *result = sum; 611 } 612 613 614 /* Here is the full subtract routine that takes care of negative numbers. 615 N2 is subtracted from N1 and the result placed in RESULT. SCALE_MIN 616 is the minimum scale for the result. */ 617 618 void 619 bc_sub (n1, n2, result, scale_min) 620 bc_num n1, n2, *result; 621 int scale_min; 622 { 623 bc_num diff; 624 int cmp_res; 625 int res_scale; 626 627 if (n1->n_sign != n2->n_sign) 628 { 629 diff = _do_add (n1, n2, scale_min); 630 diff->n_sign = n1->n_sign; 631 } 632 else 633 { 634 /* subtraction must be done. */ 635 cmp_res = _do_compare (n1, n2, FALSE, FALSE); /* Compare magnitudes. */ 636 switch (cmp_res) 637 { 638 case -1: 639 /* n1 is less than n2, subtract n1 from n2. */ 640 diff = _do_sub (n2, n1, scale_min); 641 diff->n_sign = (n2->n_sign == PLUS ? MINUS : PLUS); 642 break; 643 case 0: 644 /* They are equal! return zero! */ 645 res_scale = MAX (scale_min, MAX(n1->n_scale, n2->n_scale)); 646 diff = new_num (1, res_scale); 647 memset (diff->n_value, 0, res_scale+1); 648 break; 649 case 1: 650 /* n2 is less than n1, subtract n2 from n1. */ 651 diff = _do_sub (n1, n2, scale_min); 652 diff->n_sign = n1->n_sign; 653 break; 654 } 655 } 656 657 /* Clean up and return. */ 658 free_num (result); 659 *result = diff; 660 } 661 662 663 /* The multiply routine. N2 time N1 is put int PROD with the scale of 664 the result being MIN(N2 scale+N1 scale, MAX (SCALE, N2 scale, N1 scale)). 665 */ 666 667 void 668 bc_multiply (n1, n2, prod, scale) 669 bc_num n1, n2, *prod; 670 int scale; 671 { 672 bc_num pval; /* For the working storage. */ 673 char *n1ptr, *n2ptr, *pvptr; /* Work pointers. */ 674 char *n1end, *n2end; /* To the end of n1 and n2. */ 675 676 int indx; 677 int len1, len2, total_digits; 678 long sum; 679 int full_scale, prod_scale; 680 int toss; 681 682 /* Initialize things. */ 683 len1 = n1->n_len + n1->n_scale; 684 len2 = n2->n_len + n2->n_scale; 685 total_digits = len1 + len2; 686 full_scale = n1->n_scale + n2->n_scale; 687 prod_scale = MIN(full_scale,MAX(scale,MAX(n1->n_scale,n2->n_scale))); 688 toss = full_scale - prod_scale; 689 pval = new_num (total_digits-full_scale, prod_scale); 690 pval->n_sign = ( n1->n_sign == n2->n_sign ? PLUS : MINUS ); 691 n1end = (char *) (n1->n_value + len1 - 1); 692 n2end = (char *) (n2->n_value + len2 - 1); 693 pvptr = (char *) (pval->n_value + total_digits - toss - 1); 694 sum = 0; 695 696 /* Here are the loops... */ 697 for (indx = 0; indx < toss; indx++) 698 { 699 n1ptr = (char *) (n1end - MAX(0, indx-len2+1)); 700 n2ptr = (char *) (n2end - MIN(indx, len2-1)); 701 while ((n1ptr >= n1->n_value) && (n2ptr <= n2end)) 702 sum += *n1ptr-- * *n2ptr++; 703 sum = sum / BASE; 704 } 705 for ( ; indx < total_digits-1; indx++) 706 { 707 n1ptr = (char *) (n1end - MAX(0, indx-len2+1)); 708 n2ptr = (char *) (n2end - MIN(indx, len2-1)); 709 while ((n1ptr >= n1->n_value) && (n2ptr <= n2end)) 710 sum += *n1ptr-- * *n2ptr++; 711 *pvptr-- = sum % BASE; 712 sum = sum / BASE; 713 } 714 *pvptr-- = sum; 715 716 /* Assign to prod and clean up the number. */ 717 free_num (prod); 718 *prod = pval; 719 _rm_leading_zeros (*prod); 720 if (is_zero (*prod)) 721 (*prod)->n_sign = PLUS; 722 } 723 724 725 /* Some utility routines for the divide: First a one digit multiply. 726 NUM (with SIZE digits) is multiplied by DIGIT and the result is 727 placed into RESULT. It is written so that NUM and RESULT can be 728 the same pointers. */ 729 730 static void 731 _one_mult (num, size, digit, result) 732 unsigned char *num; 733 int size, digit; 734 unsigned char *result; 735 { 736 int carry, value; 737 unsigned char *nptr, *rptr; 738 739 if (digit == 0) 740 memset (result, 0, size); 741 else 742 { 743 if (digit == 1) 744 memcpy (result, num, size); 745 else 746 { 747 /* Initialize */ 748 nptr = (unsigned char *) (num+size-1); 749 rptr = (unsigned char *) (result+size-1); 750 carry = 0; 751 752 while (size-- > 0) 753 { 754 value = *nptr-- * digit + carry; 755 *rptr-- = value % BASE; 756 carry = value / BASE; 757 } 758 759 if (carry != 0) *rptr = carry; 760 } 761 } 762 } 763 764 765 /* The full division routine. This computes N1 / N2. It returns 766 0 if the division is ok and the result is in QUOT. The number of 767 digits after the decimal point is SCALE. It returns -1 if division 768 by zero is tried. The algorithm is found in Knuth Vol 2. p237. */ 769 770 int 771 bc_divide (n1, n2, quot, scale) 772 bc_num n1, n2, *quot; 773 int scale; 774 { 775 bc_num qval; 776 unsigned char *num1, *num2; 777 unsigned char *ptr1, *ptr2, *n2ptr, *qptr; 778 int scale1, val; 779 unsigned int len1, len2, scale2, qdigits, extra, count; 780 unsigned int qdig, qguess, borrow, carry; 781 unsigned char *mval; 782 char zero; 783 unsigned int norm; 784 785 /* Test for divide by zero. */ 786 if (is_zero (n2)) return -1; 787 788 /* Test for divide by 1. If it is we must truncate. */ 789 if (n2->n_scale == 0) 790 { 791 if (n2->n_len == 1 && *n2->n_value == 1) 792 { 793 qval = new_num (n1->n_len, scale); 794 qval->n_sign = (n1->n_sign == n2->n_sign ? PLUS : MINUS); 795 memset (&qval->n_value[n1->n_len],0,scale); 796 memcpy (qval->n_value, n1->n_value, 797 n1->n_len + MIN(n1->n_scale,scale)); 798 free_num (quot); 799 *quot = qval; 800 } 801 } 802 803 /* Set up the divide. Move the decimal point on n1 by n2's scale. 804 Remember, zeros on the end of num2 are wasted effort for dividing. */ 805 scale2 = n2->n_scale; 806 n2ptr = (unsigned char *) n2->n_value+n2->n_len+scale2-1; 807 while ((scale2 > 0) && (*n2ptr-- == 0)) scale2--; 808 809 len1 = n1->n_len + scale2; 810 scale1 = n1->n_scale - scale2; 811 if (scale1 < scale) 812 extra = scale - scale1; 813 else 814 extra = 0; 815 num1 = (unsigned char *) malloc (n1->n_len+n1->n_scale+extra+2); 816 if (num1 == NULL) out_of_memory(); 817 memset (num1, 0, n1->n_len+n1->n_scale+extra+2); 818 memcpy (num1+1, n1->n_value, n1->n_len+n1->n_scale); 819 820 len2 = n2->n_len + scale2; 821 num2 = (unsigned char *) malloc (len2+1); 822 if (num2 == NULL) out_of_memory(); 823 memcpy (num2, n2->n_value, len2); 824 *(num2+len2) = 0; 825 n2ptr = num2; 826 while (*n2ptr == 0) 827 { 828 n2ptr++; 829 len2--; 830 } 831 832 /* Calculate the number of quotient digits. */ 833 if (len2 > len1+scale) 834 { 835 qdigits = scale+1; 836 zero = TRUE; 837 } 838 else 839 { 840 zero = FALSE; 841 if (len2>len1) 842 qdigits = scale+1; /* One for the zero integer part. */ 843 else 844 qdigits = len1-len2+scale+1; 845 } 846 847 /* Allocate and zero the storage for the quotient. */ 848 qval = new_num (qdigits-scale,scale); 849 memset (qval->n_value, 0, qdigits); 850 851 /* Allocate storage for the temporary storage mval. */ 852 mval = (unsigned char *) malloc (len2+1); 853 if (mval == NULL) out_of_memory (); 854 855 /* Now for the full divide algorithm. */ 856 if (!zero) 857 { 858 /* Normalize */ 859 norm = 10 / ((int)*n2ptr + 1); 860 if (norm != 1) 861 { 862 _one_mult (num1, len1+scale1+extra+1, norm, num1); 863 _one_mult (n2ptr, len2, norm, n2ptr); 864 } 865 866 /* Initialize divide loop. */ 867 qdig = 0; 868 if (len2 > len1) 869 qptr = (unsigned char *) qval->n_value+len2-len1; 870 else 871 qptr = (unsigned char *) qval->n_value; 872 873 /* Loop */ 874 while (qdig <= len1+scale-len2) 875 { 876 /* Calculate the quotient digit guess. */ 877 if (*n2ptr == num1[qdig]) 878 qguess = 9; 879 else 880 qguess = (num1[qdig]*10 + num1[qdig+1]) / *n2ptr; 881 882 /* Test qguess. */ 883 if (n2ptr[1]*qguess > 884 (num1[qdig]*10 + num1[qdig+1] - *n2ptr*qguess)*10 885 + num1[qdig+2]) 886 { 887 qguess--; 888 /* And again. */ 889 if (n2ptr[1]*qguess > 890 (num1[qdig]*10 + num1[qdig+1] - *n2ptr*qguess)*10 891 + num1[qdig+2]) 892 qguess--; 893 } 894 895 /* Multiply and subtract. */ 896 borrow = 0; 897 if (qguess != 0) 898 { 899 *mval = 0; 900 _one_mult (n2ptr, len2, qguess, mval+1); 901 ptr1 = (unsigned char *) num1+qdig+len2; 902 ptr2 = (unsigned char *) mval+len2; 903 for (count = 0; count < len2+1; count++) 904 { 905 val = (int) *ptr1 - (int) *ptr2-- - borrow; 906 if (val < 0) 907 { 908 val += 10; 909 borrow = 1; 910 } 911 else 912 borrow = 0; 913 *ptr1-- = val; 914 } 915 } 916 917 /* Test for negative result. */ 918 if (borrow == 1) 919 { 920 qguess--; 921 ptr1 = (unsigned char *) num1+qdig+len2; 922 ptr2 = (unsigned char *) n2ptr+len2-1; 923 carry = 0; 924 for (count = 0; count < len2; count++) 925 { 926 val = (int) *ptr1 + (int) *ptr2-- + carry; 927 if (val > 9) 928 { 929 val -= 10; 930 carry = 1; 931 } 932 else 933 carry = 0; 934 *ptr1-- = val; 935 } 936 if (carry == 1) *ptr1 = (*ptr1 + 1) % 10; 937 } 938 939 /* We now know the quotient digit. */ 940 *qptr++ = qguess; 941 qdig++; 942 } 943 } 944 945 /* Clean up and return the number. */ 946 qval->n_sign = ( n1->n_sign == n2->n_sign ? PLUS : MINUS ); 947 if (is_zero (qval)) qval->n_sign = PLUS; 948 _rm_leading_zeros (qval); 949 free_num (quot); 950 *quot = qval; 951 952 /* Clean up temporary storage. */ 953 free (mval); 954 free (num1); 955 free (num2); 956 957 return 0; /* Everything is OK. */ 958 } 959 960 961 /* Division *and* modulo for numbers. This computes both NUM1 / NUM2 and 962 NUM1 % NUM2 and puts the results in QUOT and REM, except that if QUOT 963 is NULL then that store will be omitted. 964 */ 965 966 int 967 bc_divmod (num1, num2, quot, rem, scale) 968 bc_num num1, num2, *quot, *rem; 969 int scale; 970 { 971 bc_num quotient; 972 bc_num temp; 973 int rscale; 974 975 /* Check for correct numbers. */ 976 if (is_zero (num2)) return -1; 977 978 /* Calculate final scale. */ 979 rscale = MAX (num1->n_scale, num2->n_scale+scale); 980 init_num (&temp); 981 982 /* Calculate it. */ 983 bc_divide (num1, num2, &temp, scale); 984 if (quot) 985 quotient = copy_num(temp); 986 bc_multiply (temp, num2, &temp, rscale); 987 bc_sub (num1, temp, rem, rscale); 988 free_num (&temp); 989 990 if (quot) 991 { 992 free_num (quot); 993 *quot = quotient; 994 } 995 996 return 0; /* Everything is OK. */ 997 } 998 999 1000 /* Modulo for numbers. This computes NUM1 % NUM2 and puts the 1001 result in RESULT. */ 1002 1003 int 1004 bc_modulo (num1, num2, result, scale) 1005 bc_num num1, num2, *result; 1006 int scale; 1007 { 1008 return bc_divmod (num1, num2, NULL, result, scale); 1009 } 1010 1011 1012 /* Raise BASE to the EXPO power, reduced modulo MOD. The result is 1013 placed in RESULT. If a EXPO is not an integer, 1014 only the integer part is used. */ 1015 1016 int 1017 bc_raisemod (base, expo, mod, result, scale) 1018 bc_num base, expo, mod, *result; 1019 int scale; 1020 { 1021 bc_num power, exponent, parity, temp; 1022 int rscale; 1023 1024 /* Check for correct numbers. */ 1025 if (is_zero(mod)) return -1; 1026 if (is_neg(expo)) return -1; 1027 1028 /* Set initial values. */ 1029 power = copy_num (base); 1030 exponent = copy_num (expo); 1031 temp = copy_num (_one_); 1032 init_num (&parity); 1033 1034 /* Check the base for scale digits. */ 1035 if (base->n_scale != 0) 1036 rt_warn ("non-zero scale in base"); 1037 1038 /* Check the exponent for scale digits. */ 1039 if (exponent->n_scale != 0) 1040 { 1041 rt_warn ("non-zero scale in exponent"); 1042 bc_divide (exponent, _one_, &exponent, 0); /*truncate */ 1043 } 1044 1045 /* Check the modulus for scale digits. */ 1046 if (mod->n_scale != 0) 1047 rt_warn ("non-zero scale in modulus"); 1048 1049 /* Do the calculation. */ 1050 rscale = MAX(scale, base->n_scale); 1051 while ( !is_zero(exponent) ) 1052 { 1053 (void) bc_divmod (exponent, _two_, &exponent, &parity, 0); 1054 if ( !is_zero(parity) ) 1055 { 1056 bc_multiply (temp, power, &temp, rscale); 1057 (void) bc_modulo (temp, mod, &temp, scale); 1058 } 1059 1060 bc_multiply (power, power, &power, rscale); 1061 (void) bc_modulo (power, mod, &power, scale); 1062 } 1063 1064 /* Assign the value. */ 1065 free_num (&power); 1066 free_num (&exponent); 1067 free_num (result); 1068 *result = temp; 1069 return 0; /* Everything is OK. */ 1070 } 1071 1072 1073 /* Raise NUM1 to the NUM2 power. The result is placed in RESULT. 1074 Maximum exponent is LONG_MAX. If a NUM2 is not an integer, 1075 only the integer part is used. */ 1076 1077 void 1078 bc_raise (num1, num2, result, scale) 1079 bc_num num1, num2, *result; 1080 int scale; 1081 { 1082 bc_num temp, power; 1083 long exponent; 1084 int rscale; 1085 char neg; 1086 1087 /* Check the exponent for scale digits and convert to a long. */ 1088 if (num2->n_scale != 0) 1089 rt_warn ("non-zero scale in exponent"); 1090 exponent = num2long (num2); 1091 if (exponent == 0 && (num2->n_len > 1 || num2->n_value[0] != 0)) 1092 rt_error ("exponent too large in raise"); 1093 1094 /* Special case if exponent is a zero. */ 1095 if (exponent == 0) 1096 { 1097 free_num (result); 1098 *result = copy_num (_one_); 1099 return; 1100 } 1101 1102 /* Other initializations. */ 1103 if (exponent < 0) 1104 { 1105 neg = TRUE; 1106 exponent = -exponent; 1107 rscale = scale; 1108 } 1109 else 1110 { 1111 neg = FALSE; 1112 rscale = MIN (num1->n_scale*exponent, MAX(scale, num1->n_scale)); 1113 } 1114 1115 /* Set initial value of temp. */ 1116 power = copy_num (num1); 1117 while ((exponent & 1) == 0) 1118 { 1119 bc_multiply (power, power, &power, rscale); 1120 exponent = exponent >> 1; 1121 } 1122 temp = copy_num (power); 1123 exponent = exponent >> 1; 1124 1125 1126 /* Do the calculation. */ 1127 while (exponent > 0) 1128 { 1129 bc_multiply (power, power, &power, rscale); 1130 if ((exponent & 1) == 1) 1131 bc_multiply (temp, power, &temp, rscale); 1132 exponent = exponent >> 1; 1133 } 1134 1135 /* Assign the value. */ 1136 if (neg) 1137 { 1138 bc_divide (_one_, temp, result, rscale); 1139 free_num (&temp); 1140 } 1141 else 1142 { 1143 free_num (result); 1144 *result = temp; 1145 } 1146 free_num (&power); 1147 } 1148 1149 /* In some places we need to check if the number NUM is zero. */ 1150 1151 char 1152 is_near_zero (num, scale) 1153 bc_num num; 1154 int scale; 1155 { 1156 int count; 1157 char *nptr; 1158 1159 /* Initialize */ 1160 count = num->n_len + scale; 1161 nptr = num->n_value; 1162 1163 /* The check */ 1164 while ((count > 0) && (*nptr++ == 0)) count--; 1165 1166 if (count != 0 && (count != 1 || *--nptr != 1)) 1167 return FALSE; 1168 else 1169 return TRUE; 1170 } 1171 1172 /* Take the square root NUM and return it in NUM with SCALE digits 1173 after the decimal place. */ 1174 1175 int 1176 bc_sqrt (num, scale) 1177 bc_num *num; 1178 int scale; 1179 { 1180 int rscale, cmp_res, done; 1181 int cscale; 1182 bc_num guess, guess1, point5, diff; 1183 1184 /* Initial checks. */ 1185 cmp_res = bc_compare (*num, _zero_); 1186 if (cmp_res < 0) 1187 return 0; /* error */ 1188 else 1189 { 1190 if (cmp_res == 0) 1191 { 1192 free_num (num); 1193 *num = copy_num (_zero_); 1194 return 1; 1195 } 1196 } 1197 cmp_res = bc_compare (*num, _one_); 1198 if (cmp_res == 0) 1199 { 1200 free_num (num); 1201 *num = copy_num (_one_); 1202 return 1; 1203 } 1204 1205 /* Initialize the variables. */ 1206 rscale = MAX (scale, (*num)->n_scale); 1207 init_num (&guess); 1208 init_num (&guess1); 1209 init_num (&diff); 1210 point5 = new_num (1,1); 1211 point5->n_value[1] = 5; 1212 1213 1214 /* Calculate the initial guess. */ 1215 if (cmp_res < 0) 1216 { 1217 /* The number is between 0 and 1. Guess should start at 1. */ 1218 guess = copy_num (_one_); 1219 cscale = (*num)->n_scale; 1220 } 1221 else 1222 { 1223 /* The number is greater than 1. Guess should start at 10^(exp/2). */ 1224 int2num (&guess,10); 1225 1226 int2num (&guess1,(*num)->n_len); 1227 bc_multiply (guess1, point5, &guess1, 0); 1228 guess1->n_scale = 0; 1229 bc_raise (guess, guess1, &guess, 0); 1230 free_num (&guess1); 1231 cscale = 3; 1232 } 1233 1234 /* Find the square root using Newton's algorithm. */ 1235 done = FALSE; 1236 while (!done) 1237 { 1238 free_num (&guess1); 1239 guess1 = copy_num (guess); 1240 bc_divide (*num, guess, &guess, cscale); 1241 bc_add (guess, guess1, &guess, 0); 1242 bc_multiply (guess, point5, &guess, cscale); 1243 bc_sub (guess, guess1, &diff, cscale+1); 1244 if (is_near_zero (diff, cscale)) 1245 if (cscale < rscale+1) 1246 cscale = MIN (cscale*3, rscale+1); 1247 else 1248 done = TRUE; 1249 } 1250 1251 /* Assign the number and clean up. */ 1252 free_num (num); 1253 bc_divide (guess,_one_,num,rscale); 1254 free_num (&guess); 1255 free_num (&guess1); 1256 free_num (&point5); 1257 free_num (&diff); 1258 return 1; 1259 } 1260 1261 1262 /* The following routines provide output for bcd numbers package 1263 using the rules of POSIX bc for output. */ 1264 1265 /* This structure is used for saving digits in the conversion process. */ 1266 typedef struct stk_rec { 1267 long digit; 1268 struct stk_rec *next; 1269 } stk_rec; 1270 1271 /* The reference string for digits. */ 1272 char ref_str[] = "0123456789ABCDEF"; 1273 1274 1275 /* A special output routine for "multi-character digits." Exactly 1276 SIZE characters must be output for the value VAL. If SPACE is 1277 non-zero, we must output one space before the number. OUT_CHAR 1278 is the actual routine for writing the characters. */ 1279 1280 void 1281 out_long (val, size, space, out_char) 1282 long val; 1283 int size, space; 1284 #ifdef __STDC__ 1285 void (*out_char)(int); 1286 #else 1287 void (*out_char)(); 1288 #endif 1289 { 1290 char digits[40]; 1291 int len, ix; 1292 1293 if (space) (*out_char) (' '); 1294 sprintf (digits, "%ld", val); 1295 len = strlen (digits); 1296 while (size > len) 1297 { 1298 (*out_char) ('0'); 1299 size--; 1300 } 1301 for (ix=0; ix < len; ix++) 1302 (*out_char) (digits[ix]); 1303 } 1304 1305 /* Output of a bcd number. NUM is written in base O_BASE using OUT_CHAR 1306 as the routine to do the actual output of the characters. */ 1307 1308 void 1309 out_num (num, o_base, out_char) 1310 bc_num num; 1311 int o_base; 1312 #ifdef __STDC__ 1313 void (*out_char)(int); 1314 #else 1315 void (*out_char)(); 1316 #endif 1317 { 1318 char *nptr; 1319 int index, fdigit, pre_space; 1320 stk_rec *digits, *temp; 1321 bc_num int_part, frac_part, base, cur_dig, t_num, max_o_digit; 1322 1323 /* The negative sign if needed. */ 1324 if (num->n_sign == MINUS) (*out_char) ('-'); 1325 1326 /* Output the number. */ 1327 if (is_zero (num)) 1328 (*out_char) ('0'); 1329 else 1330 if (o_base == 10) 1331 { 1332 /* The number is in base 10, do it the fast way. */ 1333 nptr = num->n_value; 1334 if (num->n_len > 1 || *nptr != 0) 1335 for (index=num->n_len; index>0; index--) 1336 (*out_char) (BCD_CHAR(*nptr++)); 1337 else 1338 nptr++; 1339 1340 if (std_only && is_zero (num)) 1341 (*out_char) ('0'); 1342 1343 /* Now the fraction. */ 1344 if (num->n_scale > 0) 1345 { 1346 (*out_char) ('.'); 1347 for (index=0; index<num->n_scale; index++) 1348 (*out_char) (BCD_CHAR(*nptr++)); 1349 } 1350 } 1351 else 1352 { 1353 /* special case ... */ 1354 if (std_only && is_zero (num)) 1355 (*out_char) ('0'); 1356 1357 /* The number is some other base. */ 1358 digits = NULL; 1359 init_num (&int_part); 1360 bc_divide (num, _one_, &int_part, 0); 1361 init_num (&frac_part); 1362 init_num (&cur_dig); 1363 init_num (&base); 1364 bc_sub (num, int_part, &frac_part, 0); 1365 /* Make the INT_PART and FRAC_PART positive. */ 1366 int_part->n_sign = PLUS; 1367 frac_part->n_sign = PLUS; 1368 int2num (&base, o_base); 1369 init_num (&max_o_digit); 1370 int2num (&max_o_digit, o_base-1); 1371 1372 1373 /* Get the digits of the integer part and push them on a stack. */ 1374 while (!is_zero (int_part)) 1375 { 1376 bc_modulo (int_part, base, &cur_dig, 0); 1377 temp = (stk_rec *) malloc (sizeof(stk_rec)); 1378 if (temp == NULL) out_of_memory(); 1379 temp->digit = num2long (cur_dig); 1380 temp->next = digits; 1381 digits = temp; 1382 bc_divide (int_part, base, &int_part, 0); 1383 } 1384 1385 /* Print the digits on the stack. */ 1386 if (digits != NULL) 1387 { 1388 /* Output the digits. */ 1389 while (digits != NULL) 1390 { 1391 temp = digits; 1392 digits = digits->next; 1393 if (o_base <= 16) 1394 (*out_char) (ref_str[ (int) temp->digit]); 1395 else 1396 out_long (temp->digit, max_o_digit->n_len, 1, out_char); 1397 free (temp); 1398 } 1399 } 1400 1401 /* Get and print the digits of the fraction part. */ 1402 if (num->n_scale > 0) 1403 { 1404 (*out_char) ('.'); 1405 pre_space = 0; 1406 t_num = copy_num (_one_); 1407 while (t_num->n_len <= num->n_scale) { 1408 bc_multiply (frac_part, base, &frac_part, num->n_scale); 1409 fdigit = num2long (frac_part); 1410 int2num (&int_part, fdigit); 1411 bc_sub (frac_part, int_part, &frac_part, 0); 1412 if (o_base <= 16) 1413 (*out_char) (ref_str[fdigit]); 1414 else { 1415 out_long (fdigit, max_o_digit->n_len, pre_space, out_char); 1416 pre_space = 1; 1417 } 1418 bc_multiply (t_num, base, &t_num, 0); 1419 } 1420 free_num (&t_num); 1421 } 1422 1423 /* Clean up. */ 1424 free_num (&int_part); 1425 free_num (&frac_part); 1426 free_num (&base); 1427 free_num (&cur_dig); 1428 free_num (&max_o_digit); 1429 } 1430 } 1431 1432 #if DEBUG > 0 1433 1434 /* Debugging procedures. Some are just so one can call them from the 1435 debugger. */ 1436 1437 /* p_n prints the number NUM in base 10. */ 1438 1439 void 1440 p_n (num) 1441 bc_num num; 1442 { 1443 out_num (num, 10, out_char); 1444 } 1445 1446 1447 /* p_b prints a character array as if it was a string of bcd digits. */ 1448 void 1449 p_v (name, num, len) 1450 char *name; 1451 unsigned char *num; 1452 int len; 1453 { 1454 int i; 1455 printf ("%s=", name); 1456 for (i=0; i<len; i++) printf ("%c",BCD_CHAR(num[i])); 1457 printf ("\n"); 1458 } 1459 1460 1461 /* Convert strings to bc numbers. Base 10 only.*/ 1462 1463 void 1464 str2num (num, str, scale) 1465 bc_num *num; 1466 char *str; 1467 int scale; 1468 { 1469 int digits, strscale; 1470 char *ptr, *nptr; 1471 char zero_int; 1472 1473 /* Prepare num. */ 1474 free_num (num); 1475 1476 /* Check for valid number and count digits. */ 1477 ptr = str; 1478 digits = 0; 1479 strscale = 0; 1480 zero_int = FALSE; 1481 if ( (*ptr == '+') || (*ptr == '-')) ptr++; /* Sign */ 1482 while (*ptr == '0') ptr++; /* Skip leading zeros. */ 1483 while (isdigit(*ptr)) ptr++, digits++; /* digits */ 1484 if (*ptr == '.') ptr++; /* decimal point */ 1485 while (isdigit(*ptr)) ptr++, strscale++; /* digits */ 1486 if ((*ptr != '\0') || (digits+strscale == 0)) 1487 { 1488 *num = copy_num (_zero_); 1489 return; 1490 } 1491 1492 /* Adjust numbers and allocate storage and initialize fields. */ 1493 strscale = MIN(strscale, scale); 1494 if (digits == 0) 1495 { 1496 zero_int = TRUE; 1497 digits = 1; 1498 } 1499 *num = new_num (digits, strscale); 1500 1501 /* Build the whole number. */ 1502 ptr = str; 1503 if (*ptr == '-') 1504 { 1505 (*num)->n_sign = MINUS; 1506 ptr++; 1507 } 1508 else 1509 { 1510 (*num)->n_sign = PLUS; 1511 if (*ptr == '+') ptr++; 1512 } 1513 while (*ptr == '0') ptr++; /* Skip leading zeros. */ 1514 nptr = (*num)->n_value; 1515 if (zero_int) 1516 { 1517 *nptr++ = 0; 1518 digits = 0; 1519 } 1520 for (;digits > 0; digits--) 1521 *nptr++ = CH_VAL(*ptr++); 1522 1523 1524 /* Build the fractional part. */ 1525 if (strscale > 0) 1526 { 1527 ptr++; /* skip the decimal point! */ 1528 for (;strscale > 0; strscale--) 1529 *nptr++ = CH_VAL(*ptr++); 1530 } 1531 } 1532 1533 /* Convert a numbers to a string. Base 10 only.*/ 1534 1535 char 1536 *num2str (num) 1537 bc_num num; 1538 { 1539 char *str, *sptr; 1540 char *nptr; 1541 int index, signch; 1542 1543 /* Allocate the string memory. */ 1544 signch = ( num->n_sign == PLUS ? 0 : 1 ); /* Number of sign chars. */ 1545 if (num->n_scale > 0) 1546 str = (char *) malloc (num->n_len + num->n_scale + 2 + signch); 1547 else 1548 str = (char *) malloc (num->n_len + 1 + signch); 1549 if (str == NULL) out_of_memory(); 1550 1551 /* The negative sign if needed. */ 1552 sptr = str; 1553 if (signch) *sptr++ = '-'; 1554 1555 /* Load the whole number. */ 1556 nptr = num->n_value; 1557 for (index=num->n_len; index>0; index--) 1558 *sptr++ = BCD_CHAR(*nptr++); 1559 1560 /* Now the fraction. */ 1561 if (num->n_scale > 0) 1562 { 1563 *sptr++ = '.'; 1564 for (index=0; index<num->n_scale; index++) 1565 *sptr++ = BCD_CHAR(*nptr++); 1566 } 1567 1568 /* Terminate the string and return it! */ 1569 *sptr = '\0'; 1570 return (str); 1571 } 1572 #endif |