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