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