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