1 /* number.c: Implements arbitrary precision numbers. */ 2 3 /* This file is part of GNU bc. 4 Copyright (C) 1991, 1992, 1993, 1994, 1997, 1998 Free Software Foundation, Inc. 5 6 This program is free software; you can redistribute it and/or modify 7 it under the terms of the GNU General Public License as published by 8 the Free Software Foundation; either version 2 of the License , or 9 (at your option) any later version. 10 11 This program is distributed in the hope that it will be useful, 12 but WITHOUT ANY WARRANTY; without even the implied warranty of 13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 14 GNU General Public License for more details. 15 16 You should have received a copy of the GNU General Public License 17 along with this program; see the file COPYING. If not, write to 18 the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. 19 20 You may contact the author by: 21 e-mail: phil@cs.wwu.edu 22 us-mail: Philip A. Nelson 23 Computer Science Department, 9062 24 Western Washington University 25 Bellingham, WA 98226-9062 26 27 *************************************************************************/ 28 29 #include "bcdefs.h" 30 #include "proto.h" 31 #include "global.h" 32 33 /* Storage used for special numbers. */ 34 bc_num _zero_; 35 bc_num _one_; 36 bc_num _two_; 37 38 39 /* "Frees" a bc_num NUM. Actually decreases reference count and only 40 frees the storage if reference count is zero. */ 41 42 void 43 free_num (num) 44 bc_num *num; 45 { 46 if (*num == NULL) return; 47 (*num)->n_refs--; 48 if ((*num)->n_refs == 0) free(*num); 49 *num = NULL; 50 } 51 52 53 /* new_num allocates a number and sets fields to known values. */ 54 55 bc_num 56 new_num (length, scale) 57 int length, scale; 58 { 59 bc_num temp; 60 61 temp = (bc_num) malloc (sizeof(bc_struct)+length+scale); 62 if (temp == NULL) out_of_memory (); 63 temp->n_sign = PLUS; 64 temp->n_len = length; 65 temp->n_scale = scale; 66 temp->n_refs = 1; 67 temp->n_value[0] = 0; 68 return temp; 69 } 70 71 72 /* Intitialize the number package! */ 73 74 void 75 init_numbers () 76 { 77 _zero_ = new_num (1,0); 78 _one_ = new_num (1,0); 79 _one_->n_value[0] = 1; 80 _two_ = new_num (1,0); 81 _two_->n_value[0] = 2; 82 } 83 84 85 /* Make a copy of a number! Just increments the reference count! */ 86 87 bc_num 88 copy_num (num) 89 bc_num num; 90 { 91 num->n_refs++; 92 return num; 93 } 94 95 96 /* Initialize a number NUM by making it a copy of zero. */ 97 98 void 99 init_num (num) 100 bc_num *num; 101 { 102 *num = copy_num (_zero_); 103 } 104 105 106 /* Convert an integer VAL to a bc number NUM. */ 107 108 void 109 int2num (num, val) 110 bc_num *num; 111 int val; 112 { 113 char buffer[30]; 114 char *bptr, *vptr; 115 int ix = 1; 116 char neg = 0; 117 118 /* Sign. */ 119 if (val < 0) 120 { 121 neg = 1; 122 val = -val; 123 } 124 125 /* Get things going. */ 126 bptr = buffer; 127 *bptr++ = val % BASE; 128 val = val / BASE; 129 130 /* Extract remaining digits. */ 131 while (val != 0) 132 { 133 *bptr++ = val % BASE; 134 val = val / BASE; 135 ix++; /* Count the digits. */ 136 } 137 138 /* Make the number. */ 139 free_num (num); 140 *num = new_num (ix, 0); 141 if (neg) (*num)->n_sign = MINUS; 142 143 /* Assign the digits. */ 144 vptr = (*num)->n_value; 145 while (ix-- > 0) 146 *vptr++ = *--bptr; 147 } 148 149 150 /* Convert a number NUM to a long. The function returns only the integer 151 part of the number. For numbers that are too large to represent as 152 a long, this function returns a zero. This can be detected by checking 153 the NUM for zero after having a zero returned. */ 154 155 long 156 num2long (num) 157 bc_num num; 158 { 159 long val; 160 char *nptr; 161 int index; 162 163 /* Extract the int value, ignore the fraction. */ 164 val = 0; 165 nptr = num->n_value; 166 for (index=num->n_len; (index>0) && (val<=(LONG_MAX/BASE)); index--) 167 val = val*BASE + *nptr++; 168 169 /* Check for overflow. If overflow, return zero. */ 170 if (index>0) val = 0; 171 if (val < 0) val = 0; 172 173 /* Return the value. */ 174 if (num->n_sign == PLUS) 175 return (val); 176 else 177 return (-val); 178 } 179 180 181 /* The following are some math routines for numbers. */ 182 _PROTOTYPE(static int _do_compare, (bc_num n1, bc_num n2, int use_sign, 183 int ignore_last)); 184 _PROTOTYPE(static void _rm_leading_zeros, (bc_num num)); 185 _PROTOTYPE(static bc_num _do_add, (bc_num n1, bc_num n2, int scale_min)); 186 _PROTOTYPE(static bc_num _do_sub, (bc_num n1, bc_num n2, int scale_min)); 187 _PROTOTYPE(static void _one_mult, (unsigned char *num, int size, int digit, 188 unsigned char *result)); 189 190 191 192 /* Compare two bc numbers. Return value is 0 if equal, -1 if N1 is less 193 than N2 and +1 if N1 is greater than N2. If USE_SIGN is false, just 194 compare the magnitudes. */ 195 196 static int 197 _do_compare (n1, n2, use_sign, ignore_last) 198 bc_num n1, n2; 199 int use_sign; 200 int ignore_last; 201 { 202 char *n1ptr, *n2ptr; 203 int count; 204 205 /* First, compare signs. */ 206 if (use_sign && n1->n_sign != n2->n_sign) 207 { 208 if (n1->n_sign == PLUS) 209 return (1); /* Positive N1 > Negative N2 */ 210 else 211 return (-1); /* Negative N1 < Positive N1 */ 212 } 213 214 /* Now compare the magnitude. */ 215 if (n1->n_len != n2->n_len) 216 { 217 if (n1->n_len > n2->n_len) 218 { 219 /* Magnitude of n1 > n2. */ 220 if (!use_sign || n1->n_sign == PLUS) 221 return (1); 222 else 223 return (-1); 224 } 225 else 226 { 227 /* Magnitude of n1 < n2. */ 228 if (!use_sign || n1->n_sign == PLUS) 229 return (-1); 230 else 231 return (1); 232 } 233 } 234 235 /* If we get here, they have the same number of integer digits. 236 check the integer part and the equal length part of the fraction. */ 237 count = n1->n_len + MIN (n1->n_scale, n2->n_scale); 238 n1ptr = n1->n_value; 239 n2ptr = n2->n_value; 240 241 while ((count > 0) && (*n1ptr == *n2ptr)) 242 { 243 n1ptr++; 244 n2ptr++; 245 count--; 246 } 247 if (ignore_last && count == 1 && n1->n_scale == n2->n_scale) 248 return (0); 249 if (count != 0) 250 { 251 if (*n1ptr > *n2ptr) 252 { 253 /* Magnitude of n1 > n2. */ 254 if (!use_sign || n1->n_sign == PLUS) 255 return (1); 256 else 257 return (-1); 258 } 259 else 260 { 261 /* Magnitude of n1 < n2. */ 262 if (!use_sign || n1->n_sign == PLUS) 263 return (-1); 264 else 265 return (1); 266 } 267 } 268 269 /* They are equal up to the last part of the equal part of the fraction. */ 270 if (n1->n_scale != n2->n_scale) 271 if (n1->n_scale > n2->n_scale) 272 { 273 for (count = n1->n_scale-n2->n_scale; count>0; count--) 274 if (*n1ptr++ != 0) 275 { 276 /* Magnitude of n1 > n2. */ 277 if (!use_sign || n1->n_sign == PLUS) 278 return (1); 279 else 280 return (-1); 281 } 282 } 283 else 284 { 285 for (count = n2->n_scale-n1->n_scale; count>0; count--) 286 if (*n2ptr++ != 0) 287 { 288 /* Magnitude of n1 < n2. */ 289 if (!use_sign || n1->n_sign == PLUS) 290 return (-1); 291 else 292 return (1); 293 } 294 } 295 296 /* They must be equal! */ 297 return (0); 298 } 299 300 301 /* This is the "user callable" routine to compare numbers N1 and N2. */ 302 303 int 304 bc_compare (n1, n2) 305 bc_num n1, n2; 306 { 307 return _do_compare (n1, n2, TRUE, FALSE); 308 } 309 310 311 /* In some places we need to check if the number NUM is zero. */ 312 313 char 314 is_zero (num) 315 bc_num num; 316 { 317 int count; 318 char *nptr; 319 320 /* Quick check. */ 321 if (num == _zero_) return TRUE; 322 323 /* Initialize */ 324 count = num->n_len + num->n_scale; 325 nptr = num->n_value; 326 327 /* The check */ 328 while ((count > 0) && (*nptr++ == 0)) count--; 329 330 if (count != 0) 331 return FALSE; 332 else 333 return TRUE; 334 } 335 336 337 /* In some places we need to check if the number is negative. */ 338 339 char 340 is_neg (num) 341 bc_num num; 342 { 343 return num->n_sign == MINUS; 344 } 345 346 347 /* For many things, we may have leading zeros in a number NUM. 348 _rm_leading_zeros just moves the data to the correct 349 place and adjusts the length. */ 350 351 static void 352 _rm_leading_zeros (num) 353 bc_num num; 354 { 355 int bytes; 356 char *dst, *src; 357 358 /* Do a quick check to see if we need to do it. */ 359 if (*num->n_value != 0) return; 360 361 /* The first "digit" is 0, find the first non-zero digit in the second 362 or greater "digit" to the left of the decimal place. */ 363 bytes = num->n_len; 364 src = num->n_value; 365 while (bytes > 1 && *src == 0) src++, bytes--; 366 num->n_len = bytes; 367 bytes += num->n_scale; 368 dst = num->n_value; 369 while (bytes-- > 0) *dst++ = *src++; 370 371 } 372 373 374 /* Perform addition: N1 is added to N2 and the value is 375 returned. The signs of N1 and N2 are ignored. 376 SCALE_MIN is to set the minimum scale of the result. */ 377 378 static bc_num 379 _do_add (n1, n2, scale_min) 380 bc_num n1, n2; 381 int scale_min; 382 { 383 bc_num sum; 384 int sum_scale, sum_digits; 385 char *n1ptr, *n2ptr, *sumptr; 386 int carry, n1bytes, n2bytes; 387 int count; 388 389 /* Prepare sum. */ 390 sum_scale = MAX (n1->n_scale, n2->n_scale); 391 sum_digits = MAX (n1->n_len, n2->n_len) + 1; 392 sum = new_num (sum_digits, MAX(sum_scale, scale_min)); 393 394 /* Zero extra digits made by scale_min. */ 395 if (scale_min > sum_scale) 396 { 397 sumptr = (char *) (sum->n_value + sum_scale + sum_digits); 398 for (count = scale_min - sum_scale; count > 0; count--) 399 *sumptr++ = 0; 400 } 401 402 /* Start with the fraction part. Initialize the pointers. */ 403 n1bytes = n1->n_scale; 404 n2bytes = n2->n_scale; 405 n1ptr = (char *) (n1->n_value + n1->n_len + n1bytes - 1); 406 n2ptr = (char *) (n2->n_value + n2->n_len + n2bytes - 1); 407 sumptr = (char *) (sum->n_value + sum_scale + sum_digits - 1); 408 409 /* Add the fraction part. First copy the longer fraction.*/ 410 if (n1bytes != n2bytes) 411 { 412 if (n1bytes > n2bytes) 413 while (n1bytes>n2bytes) 414 { *sumptr-- = *n1ptr--; n1bytes--;} 415 else 416 while (n2bytes>n1bytes) 417 { *sumptr-- = *n2ptr--; n2bytes--;} 418 } 419 420 /* Now add the remaining fraction part and equal size integer parts. */ 421 n1bytes += n1->n_len; 422 n2bytes += n2->n_len; 423 carry = 0; 424 while ((n1bytes > 0) && (n2bytes > 0)) 425 { 426 *sumptr = *n1ptr-- + *n2ptr-- + carry; 427 if (*sumptr > (BASE-1)) 428 { 429 carry = 1; 430 *sumptr -= BASE; 431 } 432 else 433 carry = 0; 434 sumptr--; 435 n1bytes--; 436 n2bytes--; 437 } 438 439 /* Now add carry the longer integer part. */ 440 if (n1bytes == 0) 441 { n1bytes = n2bytes; n1ptr = n2ptr; } 442 while (n1bytes-- > 0) 443 { 444 *sumptr = *n1ptr-- + carry; 445 if (*sumptr > (BASE-1)) 446 { 447 carry = 1; 448 *sumptr -= BASE; 449 } 450 else 451 carry = 0; 452 sumptr--; 453 } 454 455 /* Set final carry. */ 456 if (carry == 1) 457 *sumptr += 1; 458 459 /* Adjust sum and return. */ 460 _rm_leading_zeros (sum); 461 return sum; 462 } 463 464 465 /* Perform subtraction: N2 is subtracted from N1 and the value is 466 returned. The signs of N1 and N2 are ignored. Also, N1 is 467 assumed to be larger than N2. SCALE_MIN is the minimum scale 468 of the result. */ 469 470 static bc_num 471 _do_sub (n1, n2, scale_min) 472 bc_num n1, n2; 473 int scale_min; 474 { 475 bc_num diff; 476 int diff_scale, diff_len; 477 int min_scale, min_len; 478 char *n1ptr, *n2ptr, *diffptr; 479 int borrow, count, val; 480 481 /* Allocate temporary storage. */ 482 diff_len = MAX (n1->n_len, n2->n_len); 483 diff_scale = MAX (n1->n_scale, n2->n_scale); 484 min_len = MIN (n1->n_len, n2->n_len); 485 min_scale = MIN (n1->n_scale, n2->n_scale); 486 diff = new_num (diff_len, MAX(diff_scale, scale_min)); 487 488 /* Zero extra digits made by scale_min. */ 489 if (scale_min > diff_scale) 490 { 491 diffptr = (char *) (diff->n_value + diff_len + diff_scale); 492 for (count = scale_min - diff_scale; count > 0; count--) 493 *diffptr++ = 0; 494 } 495 496 /* Initialize the subtract. */ 497 n1ptr = (char *) (n1->n_value + n1->n_len + n1->n_scale -1); 498 n2ptr = (char *) (n2->n_value + n2->n_len + n2->n_scale -1); 499 diffptr = (char *) (diff->n_value + diff_len + diff_scale -1); 500 501 /* Subtract the numbers. */ 502 borrow = 0; 503 504 /* Take care of the longer scaled number. */ 505 if (n1->n_scale != min_scale) 506 { 507 /* n1 has the longer scale */ 508 for (count = n1->n_scale - min_scale; count > 0; count--) 509 *diffptr-- = *n1ptr--; 510 } 511 else 512 { 513 /* n2 has the longer scale */ 514 for (count = n2->n_scale - min_scale; count > 0; count--) 515 { 516 val = - *n2ptr-- - borrow; 517 if (val < 0) 518 { 519 val += BASE; 520 borrow = 1; 521 } 522 else 523 borrow = 0; 524 *diffptr-- = val; 525 } 526 } 527 528 /* Now do the equal length scale and integer parts. */ 529 530 for (count = 0; count < min_len + min_scale; count++) 531 { 532 val = *n1ptr-- - *n2ptr-- - borrow; 533 if (val < 0) 534 { 535 val += BASE; 536 borrow = 1; 537 } 538 else 539 borrow = 0; 540 *diffptr-- = val; 541 } 542 543 /* If n1 has more digits then n2, we now do that subtract. */ 544 if (diff_len != min_len) 545 { 546 for (count = diff_len - min_len; count > 0; count--) 547 { 548 val = *n1ptr-- - borrow; 549 if (val < 0) 550 { 551 val += BASE; 552 borrow = 1; 553 } 554 else 555 borrow = 0; 556 *diffptr-- = val; 557 } 558 } 559 560 /* Clean up and return. */ 561 _rm_leading_zeros (diff); 562 return diff; 563 } 564 565 566 /* Here is the full add routine that takes care of negative numbers. 567 N1 is added to N2 and the result placed into RESULT. SCALE_MIN 568 is the minimum scale for the result. */ 569 570 void 571 bc_add (n1, n2, result, scale_min) 572 bc_num n1, n2, *result; 573 int scale_min; 574 { 575 bc_num sum; 576 int cmp_res; 577 int res_scale; 578 579 if (n1->n_sign == n2->n_sign) 580 { 581 sum = _do_add (n1, n2, scale_min); 582 sum->n_sign = n1->n_sign; 583 } 584 else 585 { 586 /* subtraction must be done. */ 587 cmp_res = _do_compare (n1, n2, FALSE, FALSE); /* Compare magnitudes. */ 588 switch (cmp_res) 589 { 590 case -1: 591 /* n1 is less than n2, subtract n1 from n2. */ 592 sum = _do_sub (n2, n1, scale_min); 593 sum->n_sign = n2->n_sign; 594 break; 595 case 0: 596 /* They are equal! return zero with the correct scale! */ 597 res_scale = MAX (scale_min, MAX(n1->n_scale, n2->n_scale)); 598 sum = new_num (1, res_scale); 599 memset (sum->n_value, 0, res_scale+1); 600 break; 601 case 1: 602 /* n2 is less than n1, subtract n2 from n1. */ 603 sum = _do_sub (n1, n2, scale_min); 604 sum->n_sign = n1->n_sign; 605 } 606 } 607 608 /* Clean up and return. */ 609 free_num (result); 610 *result = sum; 611 } 612 613 614 /* Here is the full subtract routine that takes care of negative numbers. 615 N2 is subtracted from N1 and the result placed in RESULT. SCALE_MIN 616 is the minimum scale for the result. */ 617 618 void 619 bc_sub (n1, n2, result, scale_min) 620 bc_num n1, n2, *result; 621 int scale_min; 622 { 623 bc_num diff; 624 int cmp_res; 625 int res_scale; 626 627 if (n1->n_sign != n2->n_sign) 628 { 629 diff = _do_add (n1, n2, scale_min); 630 diff->n_sign = n1->n_sign; 631 } 632 else 633 { 634 /* subtraction must be done. */ 635 cmp_res = _do_compare (n1, n2, FALSE, FALSE); /* Compare magnitudes. */ 636 switch (cmp_res) 637 { 638 case -1: 639 /* n1 is less than n2, subtract n1 from n2. */ 640 diff = _do_sub (n2, n1, scale_min); 641 diff->n_sign = (n2->n_sign == PLUS ? MINUS : PLUS); 642 break; 643 case 0: 644 /* They are equal! return zero! */ 645 res_scale = MAX (scale_min, MAX(n1->n_scale, n2->n_scale)); 646 diff = new_num (1, res_scale); 647 memset (diff->n_value, 0, res_scale+1); 648 break; 649 case 1: 650 /* n2 is less than n1, subtract n2 from n1. */ 651 diff = _do_sub (n1, n2, scale_min); 652 diff->n_sign = n1->n_sign; 653 break; 654 } 655 } 656 657 /* Clean up and return. */ 658 free_num (result); 659 *result = diff; 660 } 661 662 663 /* The multiply routine. N2 time N1 is put int PROD with the scale of 664 the result being MIN(N2 scale+N1 scale, MAX (SCALE, N2 scale, N1 scale)). 665 */ 666 667 void 668 bc_multiply (n1, n2, prod, scale) 669 bc_num n1, n2, *prod; 670 int scale; 671 { 672 bc_num pval; /* For the working storage. */ 673 char *n1ptr, *n2ptr, *pvptr; /* Work pointers. */ 674 char *n1end, *n2end; /* To the end of n1 and n2. */ 675 676 int indx; 677 int len1, len2, total_digits; 678 long sum; 679 int full_scale, prod_scale; 680 int toss; 681 682 /* Initialize things. */ 683 len1 = n1->n_len + n1->n_scale; 684 len2 = n2->n_len + n2->n_scale; 685 total_digits = len1 + len2; 686 full_scale = n1->n_scale + n2->n_scale; 687 prod_scale = MIN(full_scale,MAX(scale,MAX(n1->n_scale,n2->n_scale))); 688 toss = full_scale - prod_scale; 689 pval = new_num (total_digits-full_scale, prod_scale); 690 pval->n_sign = ( n1->n_sign == n2->n_sign ? PLUS : MINUS ); 691 n1end = (char *) (n1->n_value + len1 - 1); 692 n2end = (char *) (n2->n_value + len2 - 1); 693 pvptr = (char *) (pval->n_value + total_digits - toss - 1); 694 sum = 0; 695 696 /* Here are the loops... */ 697 for (indx = 0; indx < toss; indx++) 698 { 699 n1ptr = (char *) (n1end - MAX(0, indx-len2+1)); 700 n2ptr = (char *) (n2end - MIN(indx, len2-1)); 701 while ((n1ptr >= n1->n_value) && (n2ptr <= |