File indexing completed on 2024-05-12 17:22:40
0001 /* number.c: Implements arbitrary precision numbers. */ 0002 /* 0003 Copyright (C) 1991, 1992, 1993, 1994, 1997, 2000 Free Software Foundation, Inc. 0004 0005 This program is free software; you can redistribute it and/or modify 0006 it under the terms of the GNU General Public License as published by 0007 the Free Software Foundation; either version 2 of the License , or 0008 (at your option) any later version. 0009 0010 This program is distributed in the hope that it will be useful, 0011 but WITHOUT ANY WARRANTY; without even the implied warranty of 0012 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0013 GNU General Public License for more details. 0014 0015 You should have received a copy of the GNU General Public License 0016 along with this program; see the file COPYING. If not, write to: 0017 0018 The Free Software Foundation, Inc. 0019 59 Temple Place, Suite 330 0020 Boston, MA 02111-1307 USA. 0021 0022 0023 You may contact the author by: 0024 e-mail: philnelson@acm.org 0025 us-mail: Philip A. Nelson 0026 Computer Science Department, 9062 0027 Western Washington University 0028 Bellingham, WA 98226-9062 0029 0030 !!!This is a patched file, the original file from bc 1.06 contains bugs 0031 in (a) bc_divide and (b) bc_int2num. A patch is applied here by 0032 Wolf Lammen, Oertzweg 45, 22307 Hamburg 0033 email ookami1 <at> gmx <dot> de 0034 One patched line fixes a nasty bug, where a division by 1 may fail 0035 occasionly when an operand is overwritten by the result. 0036 The other one lets a conversion of an integer succeed, even if 0037 the most negative integer is passed as argument 0038 0039 *************************************************************************/ 0040 0041 #include "number.h" 0042 0043 #include <stdio.h> 0044 #include <assert.h> 0045 #include <stdlib.h> 0046 #include <string.h> 0047 #include <ctype.h>/* Prototypes needed for external utility routines. */ 0048 0049 #define bc_rt_warn rt_warn 0050 #define bc_rt_error rt_error 0051 #define bc_out_of_memory out_of_memory 0052 0053 _PROTOTYPE(void rt_warn, (char *mesg ,...)); 0054 _PROTOTYPE(void rt_error, (char *mesg ,...)); 0055 _PROTOTYPE(void out_of_memory, (void)); 0056 0057 0058 void out_of_memory(void){ 0059 return; 0060 } 0061 0062 void rt_warn(char *mesg ,...){ 0063 (void)mesg; 0064 return; 0065 } 0066 0067 void rt_error(char *mesg ,...){ 0068 (void)mesg; 0069 return; 0070 } 0071 0072 /* Storage used for special numbers. */ 0073 bc_num _zero_; 0074 bc_num _one_; 0075 bc_num _two_; 0076 0077 static bc_num _bc_Free_list = NULL; 0078 0079 /* new_num allocates a number and sets fields to known values. */ 0080 0081 bc_num 0082 bc_new_num (length, scale) 0083 int length, scale; 0084 { 0085 bc_num temp; 0086 0087 if (_bc_Free_list != NULL) { 0088 temp = _bc_Free_list; 0089 _bc_Free_list = temp->n_next; 0090 } else { 0091 temp = (bc_num) malloc (sizeof(bc_struct)); 0092 if (temp == NULL) bc_out_of_memory (); 0093 } 0094 temp->n_sign = PLUS; 0095 temp->n_len = length; 0096 temp->n_scale = scale; 0097 temp->n_refs = 1; 0098 temp->n_ptr = (char *) malloc (length+scale+1); 0099 if (temp->n_ptr == NULL) bc_out_of_memory(); 0100 temp->n_value = temp->n_ptr; 0101 memset (temp->n_ptr, 0, length+scale); 0102 return temp; 0103 } 0104 0105 /* "Frees" a bc_num NUM. Actually decreases reference count and only 0106 frees the storage if reference count is zero. */ 0107 0108 void 0109 bc_free_num (num) 0110 bc_num *num; 0111 { 0112 if (*num == NULL) return; 0113 (*num)->n_refs--; 0114 if ((*num)->n_refs == 0) { 0115 if ((*num)->n_ptr) 0116 free ((*num)->n_ptr); 0117 (*num)->n_next = _bc_Free_list; 0118 _bc_Free_list = *num; 0119 } 0120 *num = NULL; 0121 } 0122 0123 0124 /* Intitialize the number package! */ 0125 0126 void 0127 bc_init_numbers () 0128 { 0129 _zero_ = bc_new_num (1,0); 0130 _one_ = bc_new_num (1,0); 0131 _one_->n_value[0] = 1; 0132 _two_ = bc_new_num (1,0); 0133 _two_->n_value[0] = 2; 0134 } 0135 0136 0137 /* Make a copy of a number! Just increments the reference count! */ 0138 0139 bc_num 0140 bc_copy_num (num) 0141 bc_num num; 0142 { 0143 num->n_refs++; 0144 return num; 0145 } 0146 0147 0148 /* Initialize a number NUM by making it a copy of zero. */ 0149 0150 void 0151 bc_init_num (num) 0152 bc_num *num; 0153 { 0154 *num = bc_copy_num (_zero_); 0155 } 0156 0157 /* For many things, we may have leading zeros in a number NUM. 0158 _bc_rm_leading_zeros just moves the data "value" pointer to the 0159 correct place and adjusts the length. */ 0160 0161 static void 0162 _bc_rm_leading_zeros (num) 0163 bc_num num; 0164 { 0165 /* We can move n_value to point to the first non zero digit! */ 0166 while (*num->n_value == 0 && num->n_len > 1) { 0167 num->n_value++; 0168 num->n_len--; 0169 } 0170 } 0171 0172 0173 /* Compare two bc numbers. Return value is 0 if equal, -1 if N1 is less 0174 than N2 and +1 if N1 is greater than N2. If USE_SIGN is false, just 0175 compare the magnitudes. */ 0176 0177 static int 0178 _bc_do_compare (n1, n2, use_sign, ignore_last) 0179 bc_num n1, n2; 0180 int use_sign; 0181 int ignore_last; 0182 { 0183 char *n1ptr, *n2ptr; 0184 int count; 0185 0186 /* First, compare signs. */ 0187 if (use_sign && n1->n_sign != n2->n_sign) 0188 { 0189 if (n1->n_sign == PLUS) 0190 return (1); /* Positive N1 > Negative N2 */ 0191 else 0192 return (-1); /* Negative N1 < Positive N1 */ 0193 } 0194 0195 /* Now compare the magnitude. */ 0196 if (n1->n_len != n2->n_len) 0197 { 0198 if (n1->n_len > n2->n_len) 0199 { 0200 /* Magnitude of n1 > n2. */ 0201 if (!use_sign || n1->n_sign == PLUS) 0202 return (1); 0203 else 0204 return (-1); 0205 } 0206 else 0207 { 0208 /* Magnitude of n1 < n2. */ 0209 if (!use_sign || n1->n_sign == PLUS) 0210 return (-1); 0211 else 0212 return (1); 0213 } 0214 } 0215 0216 /* If we get here, they have the same number of integer digits. 0217 check the integer part and the equal length part of the fraction. */ 0218 count = n1->n_len + MIN (n1->n_scale, n2->n_scale); 0219 n1ptr = n1->n_value; 0220 n2ptr = n2->n_value; 0221 0222 while ((count > 0) && (*n1ptr == *n2ptr)) 0223 { 0224 n1ptr++; 0225 n2ptr++; 0226 count--; 0227 } 0228 if (ignore_last && count == 1 && n1->n_scale == n2->n_scale) 0229 return (0); 0230 if (count != 0) 0231 { 0232 if (*n1ptr > *n2ptr) 0233 { 0234 /* Magnitude of n1 > n2. */ 0235 if (!use_sign || n1->n_sign == PLUS) 0236 return (1); 0237 else 0238 return (-1); 0239 } 0240 else 0241 { 0242 /* Magnitude of n1 < n2. */ 0243 if (!use_sign || n1->n_sign == PLUS) 0244 return (-1); 0245 else 0246 return (1); 0247 } 0248 } 0249 0250 /* They are equal up to the last part of the equal part of the fraction. */ 0251 if (n1->n_scale != n2->n_scale) 0252 { 0253 if (n1->n_scale > n2->n_scale) 0254 { 0255 for (count = n1->n_scale-n2->n_scale; count>0; count--) 0256 if (*n1ptr++ != 0) 0257 { 0258 /* Magnitude of n1 > n2. */ 0259 if (!use_sign || n1->n_sign == PLUS) 0260 return (1); 0261 else 0262 return (-1); 0263 } 0264 } 0265 else 0266 { 0267 for (count = n2->n_scale-n1->n_scale; count>0; count--) 0268 if (*n2ptr++ != 0) 0269 { 0270 /* Magnitude of n1 < n2. */ 0271 if (!use_sign || n1->n_sign == PLUS) 0272 return (-1); 0273 else 0274 return (1); 0275 } 0276 } 0277 } 0278 0279 /* They must be equal! */ 0280 return (0); 0281 } 0282 0283 0284 /* This is the "user callable" routine to compare numbers N1 and N2. */ 0285 0286 int 0287 bc_compare (n1, n2) 0288 bc_num n1, n2; 0289 { 0290 return _bc_do_compare (n1, n2, TRUE, FALSE); 0291 } 0292 0293 /* In some places we need to check if the number is negative. */ 0294 0295 char 0296 bc_is_neg (num) 0297 bc_num num; 0298 { 0299 return num->n_sign == MINUS; 0300 } 0301 0302 /* In some places we need to check if the number NUM is zero. */ 0303 0304 char 0305 bc_is_zero (num) 0306 bc_num num; 0307 { 0308 int count; 0309 char *nptr; 0310 0311 /* Quick check. */ 0312 if (num == _zero_) return TRUE; 0313 0314 /* Initialize */ 0315 count = num->n_len + num->n_scale; 0316 nptr = num->n_value; 0317 0318 /* The check */ 0319 while ((count > 0) && (*nptr++ == 0)) count--; 0320 0321 if (count != 0) 0322 return FALSE; 0323 else 0324 return TRUE; 0325 } 0326 0327 /* In some places we need to check if the number NUM is almost zero. 0328 Specifically, all but the last digit is 0 and the last digit is 1. 0329 Last digit is defined by scale. */ 0330 0331 char 0332 bc_is_near_zero (num, scale) 0333 bc_num num; 0334 int scale; 0335 { 0336 int count; 0337 char *nptr; 0338 0339 /* Error checking */ 0340 if (scale > num->n_scale) 0341 scale = num->n_scale; 0342 0343 /* Initialize */ 0344 count = num->n_len + scale; 0345 nptr = num->n_value; 0346 0347 /* The check */ 0348 while ((count > 0) && (*nptr++ == 0)) count--; 0349 0350 if (count != 0 && (count != 1 || *--nptr != 1)) 0351 return FALSE; 0352 else 0353 return TRUE; 0354 } 0355 0356 0357 /* Perform addition: N1 is added to N2 and the value is 0358 returned. The signs of N1 and N2 are ignored. 0359 SCALE_MIN is to set the minimum scale of the result. */ 0360 0361 static bc_num 0362 _bc_do_add (n1, n2, scale_min) 0363 bc_num n1, n2; 0364 int scale_min; 0365 { 0366 bc_num sum; 0367 int sum_scale, sum_digits; 0368 char *n1ptr, *n2ptr, *sumptr; 0369 int carry, n1bytes, n2bytes; 0370 int count; 0371 0372 /* Prepare sum. */ 0373 sum_scale = MAX (n1->n_scale, n2->n_scale); 0374 sum_digits = MAX (n1->n_len, n2->n_len) + 1; 0375 sum = bc_new_num (sum_digits, MAX(sum_scale, scale_min)); 0376 0377 /* Zero extra digits made by scale_min. */ 0378 if (scale_min > sum_scale) 0379 { 0380 sumptr = (char *) (sum->n_value + sum_scale + sum_digits); 0381 for (count = scale_min - sum_scale; count > 0; count--) 0382 *sumptr++ = 0; 0383 } 0384 0385 /* Start with the fraction part. Initialize the pointers. */ 0386 n1bytes = n1->n_scale; 0387 n2bytes = n2->n_scale; 0388 n1ptr = (char *) (n1->n_value + n1->n_len + n1bytes - 1); 0389 n2ptr = (char *) (n2->n_value + n2->n_len + n2bytes - 1); 0390 sumptr = (char *) (sum->n_value + sum_scale + sum_digits - 1); 0391 0392 /* Add the fraction part. First copy the longer fraction.*/ 0393 if (n1bytes != n2bytes) 0394 { 0395 if (n1bytes > n2bytes) 0396 while (n1bytes>n2bytes) 0397 { *sumptr-- = *n1ptr--; n1bytes--;} 0398 else 0399 while (n2bytes>n1bytes) 0400 { *sumptr-- = *n2ptr--; n2bytes--;} 0401 } 0402 0403 /* Now add the remaining fraction part and equal size integer parts. */ 0404 n1bytes += n1->n_len; 0405 n2bytes += n2->n_len; 0406 carry = 0; 0407 while ((n1bytes > 0) && (n2bytes > 0)) 0408 { 0409 *sumptr = *n1ptr-- + *n2ptr-- + carry; 0410 if (*sumptr > (BASE-1)) 0411 { 0412 carry = 1; 0413 *sumptr -= BASE; 0414 } 0415 else 0416 carry = 0; 0417 sumptr--; 0418 n1bytes--; 0419 n2bytes--; 0420 } 0421 0422 /* Now add carry the longer integer part. */ 0423 if (n1bytes == 0) 0424 { n1bytes = n2bytes; n1ptr = n2ptr; } 0425 while (n1bytes-- > 0) 0426 { 0427 *sumptr = *n1ptr-- + carry; 0428 if (*sumptr > (BASE-1)) 0429 { 0430 carry = 1; 0431 *sumptr -= BASE; 0432 } 0433 else 0434 carry = 0; 0435 sumptr--; 0436 } 0437 0438 /* Set final carry. */ 0439 if (carry == 1) 0440 *sumptr += 1; 0441 0442 /* Adjust sum and return. */ 0443 _bc_rm_leading_zeros (sum); 0444 return sum; 0445 } 0446 0447 0448 /* Perform subtraction: N2 is subtracted from N1 and the value is 0449 returned. The signs of N1 and N2 are ignored. Also, N1 is 0450 assumed to be larger than N2. SCALE_MIN is the minimum scale 0451 of the result. */ 0452 0453 static bc_num 0454 _bc_do_sub (n1, n2, scale_min) 0455 bc_num n1, n2; 0456 int scale_min; 0457 { 0458 bc_num diff; 0459 int diff_scale, diff_len; 0460 int min_scale, min_len; 0461 char *n1ptr, *n2ptr, *diffptr; 0462 int borrow, count, val; 0463 0464 /* Allocate temporary storage. */ 0465 diff_len = MAX (n1->n_len, n2->n_len); 0466 diff_scale = MAX (n1->n_scale, n2->n_scale); 0467 min_len = MIN (n1->n_len, n2->n_len); 0468 min_scale = MIN (n1->n_scale, n2->n_scale); 0469 diff = bc_new_num (diff_len, MAX(diff_scale, scale_min)); 0470 0471 /* Zero extra digits made by scale_min. */ 0472 if (scale_min > diff_scale) 0473 { 0474 diffptr = (char *) (diff->n_value + diff_len + diff_scale); 0475 for (count = scale_min - diff_scale; count > 0; count--) 0476 *diffptr++ = 0; 0477 } 0478 0479 /* Initialize the subtract. */ 0480 n1ptr = (char *) (n1->n_value + n1->n_len + n1->n_scale -1); 0481 n2ptr = (char *) (n2->n_value + n2->n_len + n2->n_scale -1); 0482 diffptr = (char *) (diff->n_value + diff_len + diff_scale -1); 0483 0484 /* Subtract the numbers. */ 0485 borrow = 0; 0486 0487 /* Take care of the longer scaled number. */ 0488 if (n1->n_scale != min_scale) 0489 { 0490 /* n1 has the longer scale */ 0491 for (count = n1->n_scale - min_scale; count > 0; count--) 0492 *diffptr-- = *n1ptr--; 0493 } 0494 else 0495 { 0496 /* n2 has the longer scale */ 0497 for (count = n2->n_scale - min_scale; count > 0; count--) 0498 { 0499 val = - *n2ptr-- - borrow; 0500 if (val < 0) 0501 { 0502 val += BASE; 0503 borrow = 1; 0504 } 0505 else 0506 borrow = 0; 0507 *diffptr-- = val; 0508 } 0509 } 0510 0511 /* Now do the equal length scale and integer parts. */ 0512 0513 for (count = 0; count < min_len + min_scale; count++) 0514 { 0515 val = *n1ptr-- - *n2ptr-- - borrow; 0516 if (val < 0) 0517 { 0518 val += BASE; 0519 borrow = 1; 0520 } 0521 else 0522 borrow = 0; 0523 *diffptr-- = val; 0524 } 0525 0526 /* If n1 has more digits then n2, we now do that subtract. */ 0527 if (diff_len != min_len) 0528 { 0529 for (count = diff_len - min_len; count > 0; count--) 0530 { 0531 val = *n1ptr-- - borrow; 0532 if (val < 0) 0533 { 0534 val += BASE; 0535 borrow = 1; 0536 } 0537 else 0538 borrow = 0; 0539 *diffptr-- = val; 0540 } 0541 } 0542 0543 /* Clean up and return. */ 0544 _bc_rm_leading_zeros (diff); 0545 return diff; 0546 } 0547 0548 0549 /* Here is the full subtract routine that takes care of negative numbers. 0550 N2 is subtracted from N1 and the result placed in RESULT. SCALE_MIN 0551 is the minimum scale for the result. */ 0552 0553 void 0554 bc_sub (n1, n2, result, scale_min) 0555 bc_num n1, n2, *result; 0556 int scale_min; 0557 { 0558 bc_num diff = NULL; 0559 int cmp_res; 0560 int res_scale; 0561 0562 if (n1->n_sign != n2->n_sign) 0563 { 0564 diff = _bc_do_add (n1, n2, scale_min); 0565 diff->n_sign = n1->n_sign; 0566 } 0567 else 0568 { 0569 /* subtraction must be done. */ 0570 /* Compare magnitudes. */ 0571 cmp_res = _bc_do_compare (n1, n2, FALSE, FALSE); 0572 switch (cmp_res) 0573 { 0574 case -1: 0575 /* n1 is less than n2, subtract n1 from n2. */ 0576 diff = _bc_do_sub (n2, n1, scale_min); 0577 diff->n_sign = (n2->n_sign == PLUS ? MINUS : PLUS); 0578 break; 0579 case 0: 0580 /* They are equal! return zero! */ 0581 res_scale = MAX (scale_min, MAX(n1->n_scale, n2->n_scale)); 0582 diff = bc_new_num (1, res_scale); 0583 memset (diff->n_value, 0, res_scale+1); 0584 break; 0585 case 1: 0586 /* n2 is less than n1, subtract n2 from n1. */ 0587 diff = _bc_do_sub (n1, n2, scale_min); 0588 diff->n_sign = n1->n_sign; 0589 break; 0590 } 0591 } 0592 0593 /* Clean up and return. */ 0594 bc_free_num (result); 0595 *result = diff; 0596 } 0597 0598 0599 /* Here is the full add routine that takes care of negative numbers. 0600 N1 is added to N2 and the result placed into RESULT. SCALE_MIN 0601 is the minimum scale for the result. */ 0602 0603 void 0604 bc_add (n1, n2, result, scale_min) 0605 bc_num n1, n2, *result; 0606 int scale_min; 0607 { 0608 bc_num sum = NULL; 0609 int cmp_res; 0610 int res_scale; 0611 0612 if (n1->n_sign == n2->n_sign) 0613 { 0614 sum = _bc_do_add (n1, n2, scale_min); 0615 sum->n_sign = n1->n_sign; 0616 } 0617 else 0618 { 0619 /* subtraction must be done. */ 0620 cmp_res = _bc_do_compare (n1, n2, FALSE, FALSE); /* Compare magnitudes. */ 0621 switch (cmp_res) 0622 { 0623 case -1: 0624 /* n1 is less than n2, subtract n1 from n2. */ 0625 sum = _bc_do_sub (n2, n1, scale_min); 0626 sum->n_sign = n2->n_sign; 0627 break; 0628 case 0: 0629 /* They are equal! return zero with the correct scale! */ 0630 res_scale = MAX (scale_min, MAX(n1->n_scale, n2->n_scale)); 0631 sum = bc_new_num (1, res_scale); 0632 memset (sum->n_value, 0, res_scale+1); 0633 break; 0634 case 1: 0635 /* n2 is less than n1, subtract n2 from n1. */ 0636 sum = _bc_do_sub (n1, n2, scale_min); 0637 sum->n_sign = n1->n_sign; 0638 } 0639 } 0640 0641 /* Clean up and return. */ 0642 bc_free_num (result); 0643 *result = sum; 0644 } 0645 0646 /* Recursive vs non-recursive multiply crossover ranges. */ 0647 #if defined(MULDIGITS) 0648 #include "muldigits.h" 0649 #else 0650 #define MUL_BASE_DIGITS 80 0651 #endif 0652 0653 int mul_base_digits = MUL_BASE_DIGITS; 0654 #define MUL_SMALL_DIGITS mul_base_digits/4 0655 0656 /* Multiply utility routines */ 0657 0658 static bc_num 0659 new_sub_num (length, scale, value) 0660 int length, scale; 0661 char *value; 0662 { 0663 bc_num temp; 0664 0665 if (_bc_Free_list != NULL) { 0666 temp = _bc_Free_list; 0667 _bc_Free_list = temp->n_next; 0668 } else { 0669 temp = (bc_num) malloc (sizeof(bc_struct)); 0670 if (temp == NULL) bc_out_of_memory (); 0671 } 0672 temp->n_sign = PLUS; 0673 temp->n_len = length; 0674 temp->n_scale = scale; 0675 temp->n_refs = 1; 0676 temp->n_ptr = NULL; 0677 temp->n_value = value; 0678 return temp; 0679 } 0680 0681 static void 0682 _bc_simp_mul (bc_num n1, int n1len, bc_num n2, int n2len, bc_num *prod, 0683 int full_scale) 0684 { 0685 char *n1ptr, *n2ptr, *pvptr; 0686 char *n1end, *n2end; /* To the end of n1 and n2. */ 0687 int indx, sum; 0688 const int prodlen = n1len + n2len + 1; 0689 0690 (void)full_scale; 0691 0692 *prod = bc_new_num (prodlen, 0); 0693 0694 n1end = (char *) (n1->n_value + n1len - 1); 0695 n2end = (char *) (n2->n_value + n2len - 1); 0696 pvptr = (char *) ((*prod)->n_value + prodlen - 1); 0697 sum = 0; 0698 0699 /* Here is the loop... */ 0700 for (indx = 0; indx < prodlen-1; indx++) 0701 { 0702 n1ptr = (char *) (n1end - MAX(0, indx-n2len+1)); 0703 n2ptr = (char *) (n2end - MIN(indx, n2len-1)); 0704 while ((n1ptr >= n1->n_value) && (n2ptr <= n2end)) 0705 sum += *n1ptr-- * *n2ptr++; 0706 *pvptr-- = sum % BASE; 0707 sum = sum / BASE; 0708 } 0709 *pvptr = sum; 0710 } 0711 0712 0713 /* A special adder/subtractor for the recursive divide and conquer 0714 multiply algorithm. Note: if sub is called, accum must 0715 be larger that what is being subtracted. Also, accum and val 0716 must have n_scale = 0. (e.g. they must look like integers. *) */ 0717 static void 0718 _bc_shift_addsub (bc_num accum, bc_num val, int shift, int sub) 0719 { 0720 signed char *accp, *valp; 0721 int count, carry; 0722 0723 count = val->n_len; 0724 if (val->n_value[0] == 0) 0725 count--; 0726 assert (accum->n_len+accum->n_scale >= shift+count); 0727 0728 /* Set up pointers and others */ 0729 accp = (signed char *)(accum->n_value + 0730 accum->n_len + accum->n_scale - shift - 1); 0731 valp = (signed char *)(val->n_value + val->n_len - 1); 0732 carry = 0; 0733 0734 if (sub) { 0735 /* Subtraction, carry is really borrow. */ 0736 while (count--) { 0737 *accp -= *valp-- + carry; 0738 if (*accp < 0) { 0739 carry = 1; 0740 *accp-- += BASE; 0741 } else { 0742 carry = 0; 0743 accp--; 0744 } 0745 } 0746 while (carry) { 0747 *accp -= carry; 0748 if (*accp < 0) 0749 *accp-- += BASE; 0750 else 0751 carry = 0; 0752 } 0753 } else { 0754 /* Addition */ 0755 while (count--) { 0756 *accp += *valp-- + carry; 0757 if (*accp > (BASE-1)) { 0758 carry = 1; 0759 *accp-- -= BASE; 0760 } else { 0761 carry = 0; 0762 accp--; 0763 } 0764 } 0765 while (carry) { 0766 *accp += carry; 0767 if (*accp > (BASE-1)) 0768 *accp-- -= BASE; 0769 else 0770 carry = 0; 0771 } 0772 } 0773 } 0774 0775 /* Recursive divide and conquer multiply algorithm. 0776 Based on 0777 Let u = u0 + u1*(b^n) 0778 Let v = v0 + v1*(b^n) 0779 Then uv = (B^2n+B^n)*u1*v1 + B^n*(u1-u0)*(v0-v1) + (B^n+1)*u0*v0 0780 0781 B is the base of storage, number of digits in u1,u0 close to equal. 0782 */ 0783 static void 0784 _bc_rec_mul (bc_num u, int ulen, bc_num v, int vlen, bc_num *prod, 0785 int full_scale) 0786 { 0787 bc_num u0, u1, v0, v1; 0788 /*int u0len, v0len;*/ 0789 bc_num m1, m2, m3, d1, d2; 0790 int n, prodlen, m1zero; 0791 int d1len, d2len; 0792 0793 /* Base case? */ 0794 if ((ulen+vlen) < mul_base_digits 0795 || ulen < MUL_SMALL_DIGITS 0796 || vlen < MUL_SMALL_DIGITS ) { 0797 _bc_simp_mul (u, ulen, v, vlen, prod, full_scale); 0798 return; 0799 } 0800 0801 /* Calculate n -- the u and v split point in digits. */ 0802 n = (MAX(ulen, vlen)+1) / 2; 0803 0804 /* Split u and v. */ 0805 if (ulen < n) { 0806 u1 = bc_copy_num (_zero_); 0807 u0 = new_sub_num (ulen,0, u->n_value); 0808 } else { 0809 u1 = new_sub_num (ulen-n, 0, u->n_value); 0810 u0 = new_sub_num (n, 0, u->n_value+ulen-n); 0811 } 0812 if (vlen < n) { 0813 v1 = bc_copy_num (_zero_); 0814 v0 = new_sub_num (vlen,0, v->n_value); 0815 } else { 0816 v1 = new_sub_num (vlen-n, 0, v->n_value); 0817 v0 = new_sub_num (n, 0, v->n_value+vlen-n); 0818 } 0819 _bc_rm_leading_zeros (u1); 0820 _bc_rm_leading_zeros (u0); 0821 /*u0len = u0->n_len;*/ 0822 _bc_rm_leading_zeros (v1); 0823 _bc_rm_leading_zeros (v0); 0824 /*v0len = v0->n_len;*/ 0825 0826 m1zero = bc_is_zero(u1) || bc_is_zero(v1); 0827 0828 /* Calculate sub results ... */ 0829 0830 bc_init_num(&d1); 0831 bc_init_num(&d2); 0832 bc_sub (u1, u0, &d1, 0); 0833 d1len = d1->n_len; 0834 bc_sub (v0, v1, &d2, 0); 0835 d2len = d2->n_len; 0836 0837 0838 /* Do recursive multiplies and shifted adds. */ 0839 if (m1zero) 0840 m1 = bc_copy_num (_zero_); 0841 else 0842 _bc_rec_mul (u1, u1->n_len, v1, v1->n_len, &m1, 0); 0843 0844 if (bc_is_zero(d1) || bc_is_zero(d2)) 0845 m2 = bc_copy_num (_zero_); 0846 else 0847 _bc_rec_mul (d1, d1len, d2, d2len, &m2, 0); 0848 0849 if (bc_is_zero(u0) || bc_is_zero(v0)) 0850 m3 = bc_copy_num (_zero_); 0851 else 0852 _bc_rec_mul (u0, u0->n_len, v0, v0->n_len, &m3, 0); 0853 0854 /* Initialize product */ 0855 prodlen = ulen+vlen+1; 0856 *prod = bc_new_num(prodlen, 0); 0857 0858 if (!m1zero) { 0859 _bc_shift_addsub (*prod, m1, 2*n, 0); 0860 _bc_shift_addsub (*prod, m1, n, 0); 0861 } 0862 _bc_shift_addsub (*prod, m3, n, 0); 0863 _bc_shift_addsub (*prod, m3, 0, 0); 0864 _bc_shift_addsub (*prod, m2, n, d1->n_sign != d2->n_sign); 0865 0866 /* Now clean up! */ 0867 bc_free_num (&u1); 0868 bc_free_num (&u0); 0869 bc_free_num (&v1); 0870 bc_free_num (&m1); 0871 bc_free_num (&v0); 0872 bc_free_num (&m2); 0873 bc_free_num (&m3); 0874 bc_free_num (&d1); 0875 bc_free_num (&d2); 0876 } 0877 0878 /* The multiply routine. N2 times N1 is put int PROD with the scale of 0879 the result being MIN(N2 scale+N1 scale, MAX (SCALE, N2 scale, N1 scale)). 0880 */ 0881 0882 void 0883 bc_multiply (n1, n2, prod, scale) 0884 bc_num n1, n2, *prod; 0885 int scale; 0886 { 0887 bc_num pval; 0888 int len1, len2; 0889 int full_scale, prod_scale; 0890 0891 /* Initialize things. */ 0892 len1 = n1->n_len + n1->n_scale; 0893 len2 = n2->n_len + n2->n_scale; 0894 full_scale = n1->n_scale + n2->n_scale; 0895 prod_scale = MIN(full_scale,MAX(scale,MAX(n1->n_scale,n2->n_scale))); 0896 0897 /* Do the multiply */ 0898 _bc_rec_mul (n1, len1, n2, len2, &pval, full_scale); 0899 0900 /* Assign to prod and clean up the number. */ 0901 pval->n_sign = ( n1->n_sign == n2->n_sign ? PLUS : MINUS ); 0902 pval->n_value = pval->n_ptr; 0903 pval->n_len = len2 + len1 + 1 - full_scale; 0904 pval->n_scale = prod_scale; 0905 _bc_rm_leading_zeros (pval); 0906 if (bc_is_zero (pval)) 0907 pval->n_sign = PLUS; 0908 bc_free_num (prod); 0909 *prod = pval; 0910 } 0911 0912 /* Some utility routines for the divide: First a one digit multiply. 0913 NUM (with SIZE digits) is multiplied by DIGIT and the result is 0914 placed into RESULT. It is written so that NUM and RESULT can be 0915 the same pointers. */ 0916 0917 static void 0918 _one_mult (num, size, digit, result) 0919 unsigned char *num; 0920 int size, digit; 0921 unsigned char *result; 0922 { 0923 int carry, value; 0924 unsigned char *nptr, *rptr; 0925 0926 if (digit == 0) 0927 memset (result, 0, size); 0928 else 0929 { 0930 if (digit == 1) 0931 memcpy (result, num, size); 0932 else 0933 { 0934 /* Initialize */ 0935 nptr = (unsigned char *) (num+size-1); 0936 rptr = (unsigned char *) (result+size-1); 0937 carry = 0; 0938 0939 while (size-- > 0) 0940 { 0941 value = *nptr-- * digit + carry; 0942 *rptr-- = value % BASE; 0943 carry = value / BASE; 0944 } 0945 0946 if (carry != 0) *rptr = carry; 0947 } 0948 } 0949 } 0950 0951 0952 /* The full division routine. This computes N1 / N2. It returns 0953 0 if the division is ok and the result is in QUOT. The number of 0954 digits after the decimal point is SCALE. It returns -1 if division 0955 by zero is tried. The algorithm is found in Knuth Vol 2. p237. */ 0956 0957 int 0958 bc_divide (n1, n2, quot, scale) 0959 bc_num n1, n2, *quot; 0960 int scale; 0961 { 0962 bc_num qval; 0963 unsigned char *num1, *num2; 0964 unsigned char *ptr1, *ptr2, *n2ptr, *qptr; 0965 int scale1, val; 0966 unsigned int len1, len2, scale2, qdigits, extra, count; 0967 unsigned int qdig, qguess, borrow, carry; 0968 unsigned char *mval; 0969 char zero; 0970 unsigned int norm; 0971 0972 /* Test for divide by zero. */ 0973 if (bc_is_zero (n2)) return -1; 0974 0975 /* Test for divide by 1. If it is we must truncate. */ 0976 if (n2->n_scale == 0) 0977 { 0978 if (n2->n_len == 1 && *n2->n_value == 1) 0979 { 0980 qval = bc_new_num (n1->n_len, scale); 0981 qval->n_sign = (n1->n_sign == n2->n_sign ? PLUS : MINUS); 0982 memset (&qval->n_value[n1->n_len],0,scale); 0983 memcpy (qval->n_value, n1->n_value, 0984 n1->n_len + MIN(n1->n_scale,scale)); 0985 bc_free_num (quot); 0986 *quot = qval; 0987 return 0; /* bug fix Wolf Lammen */ 0988 } 0989 } 0990 0991 /* Set up the divide. Move the decimal point on n1 by n2's scale. 0992 Remember, zeros on the end of num2 are wasted effort for dividing. */ 0993 scale2 = n2->n_scale; 0994 n2ptr = (unsigned char *) n2->n_value+n2->n_len+scale2-1; 0995 while ((scale2 > 0) && (*n2ptr-- == 0)) scale2--; 0996 0997 len1 = n1->n_len + scale2; 0998 scale1 = n1->n_scale - scale2; 0999 if (scale1 < scale) 1000 extra = scale - scale1; 1001 else 1002 extra = 0; 1003 num1 = (unsigned char *) malloc (n1->n_len+n1->n_scale+extra+2); 1004 if (num1 == NULL) bc_out_of_memory(); 1005 memset (num1, 0, n1->n_len+n1->n_scale+extra+2); 1006 memcpy (num1+1, n1->n_value, n1->n_len+n1->n_scale); 1007 1008 len2 = n2->n_len + scale2; 1009 num2 = (unsigned char *) malloc (len2+1); 1010 if (num2 == NULL) bc_out_of_memory(); 1011 memcpy (num2, n2->n_value, len2); 1012 *(num2+len2) = 0; 1013 n2ptr = num2; 1014 while (*n2ptr == 0) 1015 { 1016 n2ptr++; 1017 len2--; 1018 } 1019 1020 /* Calculate the number of quotient digits. */ 1021 if (len2 > len1+scale) 1022 { 1023 qdigits = scale+1; 1024 zero = TRUE; 1025 } 1026 else 1027 { 1028 zero = FALSE; 1029 if (len2>len1) 1030 qdigits = scale+1; /* One for the zero integer part. */ 1031 else 1032 qdigits = len1-len2+scale+1; 1033 } 1034 1035 /* Allocate and zero the storage for the quotient. */ 1036 qval = bc_new_num (qdigits-scale,scale); 1037 memset (qval->n_value, 0, qdigits); 1038 1039 /* Allocate storage for the temporary storage mval. */ 1040 mval = (unsigned char *) malloc (len2+1); 1041 if (mval == NULL) bc_out_of_memory (); 1042 1043 /* Now for the full divide algorithm. */ 1044 if (!zero) 1045 { 1046 /* Normalize */ 1047 norm = 10 / ((int)*n2ptr + 1); 1048 if (norm != 1) 1049 { 1050 _one_mult (num1, len1+scale1+extra+1, norm, num1); 1051 _one_mult (n2ptr, len2, norm, n2ptr); 1052 } 1053 1054 /* Initialize divide loop. */ 1055 qdig = 0; 1056 if (len2 > len1) 1057 qptr = (unsigned char *) qval->n_value+len2-len1; 1058 else 1059 qptr = (unsigned char *) qval->n_value; 1060 1061 /* Loop */ 1062 while (qdig <= len1+scale-len2) 1063 { 1064 /* Calculate the quotient digit guess. */ 1065 if (*n2ptr == num1[qdig]) 1066 qguess = 9; 1067 else 1068 qguess = (num1[qdig]*10 + num1[qdig+1]) / *n2ptr; 1069 1070 /* Test qguess. */ 1071 if (n2ptr[1]*qguess > 1072 (num1[qdig]*10 + num1[qdig+1] - *n2ptr*qguess)*10 1073 + num1[qdig+2]) 1074 { 1075 qguess--; 1076 /* And again. */ 1077 if (n2ptr[1]*qguess > 1078 (num1[qdig]*10 + num1[qdig+1] - *n2ptr*qguess)*10 1079 + num1[qdig+2]) 1080 qguess--; 1081 } 1082 1083 /* Multiply and subtract. */ 1084 borrow = 0; 1085 if (qguess != 0) 1086 { 1087 *mval = 0; 1088 _one_mult (n2ptr, len2, qguess, mval+1); 1089 ptr1 = (unsigned char *) num1+qdig+len2; 1090 ptr2 = (unsigned char *) mval+len2; 1091 for (count = 0; count < len2+1; count++) 1092 { 1093 val = (int) *ptr1 - (int) *ptr2-- - borrow; 1094 if (val < 0) 1095 { 1096 val += 10; 1097 borrow = 1; 1098 } 1099 else 1100 borrow = 0; 1101 *ptr1-- = val; 1102 } 1103 } 1104 1105 /* Test for negative result. */ 1106 if (borrow == 1) 1107 { 1108 qguess--; 1109 ptr1 = (unsigned char *) num1+qdig+len2; 1110 ptr2 = (unsigned char *) n2ptr+len2-1; 1111 carry = 0; 1112 for (count = 0; count < len2; count++) 1113 { 1114 val = (int) *ptr1 + (int) *ptr2-- + carry; 1115 if (val > 9) 1116 { 1117 val -= 10; 1118 carry = 1; 1119 } 1120 else 1121 carry = 0; 1122 *ptr1-- = val; 1123 } 1124 if (carry == 1) *ptr1 = (*ptr1 + 1) % 10; 1125 } 1126 1127 /* We now know the quotient digit. */ 1128 *qptr++ = qguess; 1129 qdig++; 1130 } 1131 } 1132 1133 /* Clean up and return the number. */ 1134 qval->n_sign = ( n1->n_sign == n2->n_sign ? PLUS : MINUS ); 1135 if (bc_is_zero (qval)) qval->n_sign = PLUS; 1136 _bc_rm_leading_zeros (qval); 1137 bc_free_num (quot); 1138 *quot = qval; 1139 1140 /* Clean up temporary storage. */ 1141 free (mval); 1142 free (num1); 1143 free (num2); 1144 1145 return 0; /* Everything is OK. */ 1146 } 1147 1148 1149 /* Division *and* modulo for numbers. This computes both NUM1 / NUM2 and 1150 NUM1 % NUM2 and puts the results in QUOT and REM, except that if QUOT 1151 is NULL then that store will be omitted. 1152 */ 1153 1154 int 1155 bc_divmod (num1, num2, quot, rem, scale) 1156 bc_num num1, num2, *quot, *rem; 1157 int scale; 1158 { 1159 bc_num quotient = NULL; 1160 bc_num temp; 1161 int rscale; 1162 1163 /* Check for correct numbers. */ 1164 if (bc_is_zero (num2)) return -1; 1165 1166 /* Calculate final scale. */ 1167 rscale = MAX (num1->n_scale, num2->n_scale+scale); 1168 bc_init_num(&temp); 1169 1170 /* Calculate it. */ 1171 bc_divide (num1, num2, &temp, scale); 1172 if (quot) 1173 quotient = bc_copy_num (temp); 1174 bc_multiply (temp, num2, &temp, rscale); 1175 bc_sub (num1, temp, rem, rscale); 1176 bc_free_num (&temp); 1177 1178 if (quot) 1179 { 1180 bc_free_num (quot); 1181 *quot = quotient; 1182 } 1183 1184 return 0; /* Everything is OK. */ 1185 } 1186 1187 1188 /* Modulo for numbers. This computes NUM1 % NUM2 and puts the 1189 result in RESULT. */ 1190 1191 int 1192 bc_modulo (num1, num2, result, scale) 1193 bc_num num1, num2, *result; 1194 int scale; 1195 { 1196 return bc_divmod (num1, num2, NULL, result, scale); 1197 } 1198 1199 /* Raise BASE to the EXPO power, reduced modulo MOD. The result is 1200 placed in RESULT. If a EXPO is not an integer, 1201 only the integer part is used. */ 1202 1203 int 1204 bc_raisemod (base, expo, mod, result, scale) 1205 bc_num base, expo, mod, *result; 1206 int scale; 1207 { 1208 bc_num power, exponent, parity, temp; 1209 int rscale; 1210 1211 /* Check for correct numbers. */ 1212 if (bc_is_zero(mod)) return -1; 1213 if (bc_is_neg(expo)) return -1; 1214 1215 /* Set initial values. */ 1216 power = bc_copy_num (base); 1217 exponent = bc_copy_num (expo); 1218 temp = bc_copy_num (_one_); 1219 bc_init_num(&parity); 1220 1221 /* Check the base for scale digits. */ 1222 if (base->n_scale != 0) 1223 bc_rt_warn ("non-zero scale in base"); 1224 1225 /* Check the exponent for scale digits. */ 1226 if (exponent->n_scale != 0) 1227 { 1228 bc_rt_warn ("non-zero scale in exponent"); 1229 bc_divide (exponent, _one_, &exponent, 0); /*truncate */ 1230 } 1231 1232 /* Check the modulus for scale digits. */ 1233 if (mod->n_scale != 0) 1234 bc_rt_warn ("non-zero scale in modulus"); 1235 1236 /* Do the calculation. */ 1237 rscale = MAX(scale, base->n_scale); 1238 while ( !bc_is_zero(exponent) ) 1239 { 1240 (void) bc_divmod (exponent, _two_, &exponent, &parity, 0); 1241 if ( !bc_is_zero(parity) ) 1242 { 1243 bc_multiply (temp, power, &temp, rscale); 1244 (void) bc_modulo (temp, mod, &temp, scale); 1245 } 1246 1247 bc_multiply (power, power, &power, rscale); 1248 (void) bc_modulo (power, mod, &power, scale); 1249 } 1250 1251 /* Assign the value. */ 1252 bc_free_num (&power); 1253 bc_free_num (&exponent); 1254 bc_free_num (result); 1255 *result = temp; 1256 return 0; /* Everything is OK. */ 1257 } 1258 1259 /* Raise NUM1 to the NUM2 power. The result is placed in RESULT. 1260 Maximum exponent is LONG_MAX. If a NUM2 is not an integer, 1261 only the integer part is used. */ 1262 1263 void 1264 bc_raise (num1, num2, result, scale) 1265 bc_num num1, num2, *result; 1266 int scale; 1267 { 1268 bc_num temp, power; 1269 long exponent; 1270 int rscale; 1271 int pwrscale; 1272 int calcscale; 1273 char neg; 1274 1275 /* Check the exponent for scale digits and convert to a long. */ 1276 if (num2->n_scale != 0) 1277 bc_rt_warn ("non-zero scale in exponent"); 1278 exponent = bc_num2long (num2); 1279 if (exponent == 0 && (num2->n_len > 1 || num2->n_value[0] != 0)) 1280 bc_rt_error ("exponent too large in raise"); 1281 1282 /* Special case if exponent is a zero. */ 1283 if (exponent == 0) 1284 { 1285 bc_free_num (result); 1286 *result = bc_copy_num (_one_); 1287 return; 1288 } 1289 1290 /* Other initializations. */ 1291 if (exponent < 0) 1292 { 1293 neg = TRUE; 1294 exponent = -exponent; 1295 rscale = scale; 1296 } 1297 else 1298 { 1299 neg = FALSE; 1300 rscale = MIN (num1->n_scale*exponent, MAX(scale, num1->n_scale)); 1301 } 1302 1303 /* Set initial value of temp. */ 1304 power = bc_copy_num (num1); 1305 pwrscale = num1->n_scale; 1306 while ((exponent & 1) == 0) 1307 { 1308 pwrscale = 2*pwrscale; 1309 bc_multiply (power, power, &power, pwrscale); 1310 exponent = exponent >> 1; 1311 } 1312 temp = bc_copy_num (power); 1313 calcscale = pwrscale; 1314 exponent = exponent >> 1; 1315 1316 /* Do the calculation. */ 1317 while (exponent > 0) 1318 { 1319 pwrscale = 2*pwrscale; 1320 bc_multiply (power, power, &power, pwrscale); 1321 if ((exponent & 1) == 1) { 1322 calcscale = pwrscale + calcscale; 1323 bc_multiply (temp, power, &temp, calcscale); 1324 } 1325 exponent = exponent >> 1; 1326 } 1327 1328 /* Assign the value. */ 1329 if (neg) 1330 { 1331 bc_divide (_one_, temp, result, rscale); 1332 bc_free_num (&temp); 1333 } 1334 else 1335 { 1336 bc_free_num (result); 1337 *result = temp; 1338 if ((*result)->n_scale > rscale) 1339 (*result)->n_scale = rscale; 1340 } 1341 bc_free_num (&power); 1342 } 1343 1344 /* Take the square root NUM and return it in NUM with SCALE digits 1345 after the decimal place. */ 1346 1347 int 1348 bc_sqrt (num, scale) 1349 bc_num *num; 1350 int scale; 1351 { 1352 int rscale, cmp_res, done; 1353 int cscale; 1354 bc_num guess, guess1, point5, diff; 1355 1356 /* Initial checks. */ 1357 cmp_res = bc_compare (*num, _zero_); 1358 if (cmp_res < 0) 1359 return 0; /* error */ 1360 else 1361 { 1362 if (cmp_res == 0) 1363 { 1364 bc_free_num (num); 1365 *num = bc_copy_num (_zero_); 1366 return 1; 1367 } 1368 } 1369 cmp_res = bc_compare (*num, _one_); 1370 if (cmp_res == 0) 1371 { 1372 bc_free_num (num); 1373 *num = bc_copy_num (_one_); 1374 return 1; 1375 } 1376 1377 /* Initialize the variables. */ 1378 rscale = MAX (scale, (*num)->n_scale); 1379 bc_init_num(&guess); 1380 bc_init_num(&guess1); 1381 bc_init_num(&diff); 1382 point5 = bc_new_num (1,1); 1383 point5->n_value[1] = 5; 1384 1385 1386 /* Calculate the initial guess. */ 1387 if (cmp_res < 0) 1388 { 1389 /* The number is between 0 and 1. Guess should start at 1. */ 1390 guess = bc_copy_num (_one_); 1391 cscale = (*num)->n_scale; 1392 } 1393 else 1394 { 1395 /* The number is greater than 1. Guess should start at 10^(exp/2). */ 1396 bc_int2num (&guess,10); 1397 1398 bc_int2num (&guess1,(*num)->n_len); 1399 bc_multiply (guess1, point5, &guess1, 0); 1400 guess1->n_scale = 0; 1401 bc_raise (guess, guess1, &guess, 0); 1402 bc_free_num (&guess1); 1403 cscale = 3; 1404 } 1405 1406 /* Find the square root using Newton's algorithm. */ 1407 done = FALSE; 1408 while (!done) 1409 { 1410 bc_free_num (&guess1); 1411 guess1 = bc_copy_num (guess); 1412 bc_divide (*num, guess, &guess, cscale); 1413 bc_add (guess, guess1, &guess, 0); 1414 bc_multiply (guess, point5, &guess, cscale); 1415 bc_sub (guess, guess1, &diff, cscale+1); 1416 if (bc_is_near_zero (diff, cscale)) 1417 { 1418 if (cscale < rscale+1) 1419 cscale = MIN (cscale*3, rscale+1); 1420 else 1421 done = TRUE; 1422 } 1423 } 1424 1425 /* Assign the number and clean up. */ 1426 bc_free_num (num); 1427 bc_divide (guess,_one_,num,rscale); 1428 bc_free_num (&guess); 1429 bc_free_num (&guess1); 1430 bc_free_num (&point5); 1431 bc_free_num (&diff); 1432 return 1; 1433 } 1434 1435 1436 /* The following routines provide output for bcd numbers package 1437 using the rules of POSIX bc for output. */ 1438 1439 /* This structure is used for saving digits in the conversion process. */ 1440 typedef struct stk_rec { 1441 long digit; 1442 struct stk_rec *next; 1443 } stk_rec; 1444 1445 /* The reference string for digits. */ 1446 static char ref_str[] = "0123456789ABCDEF"; 1447 1448 1449 /* A special output routine for "multi-character digits." Exactly 1450 SIZE characters must be output for the value VAL. If SPACE is 1451 non-zero, we must output one space before the number. OUT_CHAR 1452 is the actual routine for writing the characters. */ 1453 1454 void 1455 bc_out_long (val, size, space, out_char) 1456 long val; 1457 int size, space; 1458 #ifdef NUMBER__STDC__ 1459 void (*out_char)(int); 1460 #else 1461 void (*out_char)(); 1462 #endif 1463 { 1464 char digits[40]; 1465 int len, ix; 1466 1467 if (space) (*out_char) (' '); 1468 sprintf (digits, "%ld", val); 1469 len = strlen (digits); 1470 while (size > len) 1471 { 1472 (*out_char) ('0'); 1473 size--; 1474 } 1475 for (ix=0; ix < len; ix++) 1476 (*out_char) (digits[ix]); 1477 } 1478 1479 /* Output of a bcd number. NUM is written in base O_BASE using OUT_CHAR 1480 as the routine to do the actual output of the characters. */ 1481 1482 void 1483 bc_out_num (num, o_base, out_char, leading_zero) 1484 bc_num num; 1485 int o_base; 1486 #ifdef NUMBER__STDC__ 1487 void (*out_char)(int); 1488 #else 1489 void (*out_char)(); 1490 #endif 1491 int leading_zero; 1492 { 1493 char *nptr; 1494 int index, fdigit, pre_space; 1495 stk_rec *digits, *temp; 1496 bc_num int_part, frac_part, base, cur_dig, t_num, max_o_digit; 1497 1498 /* The negative sign if needed. */ 1499 if (num->n_sign == MINUS) (*out_char) ('-'); 1500 1501 /* Output the number. */ 1502 if (bc_is_zero (num)) 1503 (*out_char) ('0'); 1504 else 1505 if (o_base == 10) 1506 { 1507 /* The number is in base 10, do it the fast way. */ 1508 nptr = num->n_value; 1509 if (num->n_len > 1 || *nptr != 0) 1510 for (index=num->n_len; index>0; index--) 1511 (*out_char) (BCD_CHAR(*nptr++)); 1512 else 1513 nptr++; 1514 1515 if (leading_zero && bc_is_zero (num)) 1516 (*out_char) ('0'); 1517 1518 /* Now the fraction. */ 1519 if (num->n_scale > 0) 1520 { 1521 (*out_char) ('.'); 1522 for (index=0; index<num->n_scale; index++) 1523 (*out_char) (BCD_CHAR(*nptr++)); 1524 } 1525 } 1526 else 1527 { 1528 /* special case ... */ 1529 if (leading_zero && bc_is_zero (num)) 1530 (*out_char) ('0'); 1531 1532 /* The number is some other base. */ 1533 digits = NULL; 1534 bc_init_num (&int_part); 1535 bc_divide (num, _one_, &int_part, 0); 1536 bc_init_num (&frac_part); 1537 bc_init_num (&cur_dig); 1538 bc_init_num (&base); 1539 bc_sub (num, int_part, &frac_part, 0); 1540 /* Make the INT_PART and FRAC_PART positive. */ 1541 int_part->n_sign = PLUS; 1542 frac_part->n_sign = PLUS; 1543 bc_int2num (&base, o_base); 1544 bc_init_num (&max_o_digit); 1545 bc_int2num (&max_o_digit, o_base-1); 1546 1547 1548 /* Get the digits of the integer part and push them on a stack. */ 1549 while (!bc_is_zero (int_part)) 1550 { 1551 bc_modulo (int_part, base, &cur_dig, 0); 1552 temp = (stk_rec *) malloc (sizeof(stk_rec)); 1553 if (temp == NULL) bc_out_of_memory(); 1554 temp->digit = bc_num2long (cur_dig); 1555 temp->next = digits; 1556 digits = temp; 1557 bc_divide (int_part, base, &int_part, 0); 1558 } 1559 1560 /* Print the digits on the stack. */ 1561 if (digits != NULL) 1562 { 1563 /* Output the digits. */ 1564 while (digits != NULL) 1565 { 1566 temp = digits; 1567 digits = digits->next; 1568 if (o_base <= 16) 1569 (*out_char) (ref_str[ (int) temp->digit]); 1570 else 1571 bc_out_long (temp->digit, max_o_digit->n_len, 1, out_char); 1572 free (temp); 1573 } 1574 } 1575 1576 /* Get and print the digits of the fraction part. */ 1577 if (num->n_scale > 0) 1578 { 1579 (*out_char) ('.'); 1580 pre_space = 0; 1581 t_num = bc_copy_num (_one_); 1582 while (t_num->n_len <= num->n_scale) { 1583 bc_multiply (frac_part, base, &frac_part, num->n_scale); 1584 fdigit = bc_num2long (frac_part); 1585 bc_int2num (&int_part, fdigit); 1586 bc_sub (frac_part, int_part, &frac_part, 0); 1587 if (o_base <= 16) 1588 (*out_char) (ref_str[fdigit]); 1589 else { 1590 bc_out_long (fdigit, max_o_digit->n_len, pre_space, out_char); 1591 pre_space = 1; 1592 } 1593 bc_multiply (t_num, base, &t_num, 0); 1594 } 1595 bc_free_num (&t_num); 1596 } 1597 1598 /* Clean up. */ 1599 bc_free_num (&int_part); 1600 bc_free_num (&frac_part); 1601 bc_free_num (&base); 1602 bc_free_num (&cur_dig); 1603 bc_free_num (&max_o_digit); 1604 } 1605 } 1606 /* Convert a number NUM to a long. The function returns only the integer 1607 part of the number. For numbers that are too large to represent as 1608 a long, this function returns a zero. This can be detected by checking 1609 the NUM for zero after having a zero returned. */ 1610 1611 long 1612 bc_num2long (num) 1613 bc_num num; 1614 { 1615 long val; 1616 char *nptr; 1617 int index; 1618 1619 /* Extract the int value, ignore the fraction. */ 1620 val = 0; 1621 nptr = num->n_value; 1622 for (index=num->n_len; (index>0) && (val<=(LONG_MAX/BASE)); index--) 1623 val = val*BASE + *nptr++; 1624 1625 /* Check for overflow. If overflow, return zero. */ 1626 if (index>0) val = 0; 1627 if (val < 0) val = 0; 1628 1629 /* Return the value. */ 1630 if (num->n_sign == PLUS) 1631 return (val); 1632 else 1633 return (-val); 1634 } 1635 1636 1637 /* Convert an integer VAL to a bc number NUM. */ 1638 1639 void 1640 bc_int2num (num, val) 1641 bc_num *num; 1642 int val; 1643 { 1644 char buffer[30]; 1645 char *bptr, *vptr; 1646 int ix = 1; 1647 char neg = 0; 1648 1649 /* Sign. */ 1650 if (val < 0) 1651 { 1652 neg = 1; 1653 val = -val; 1654 } 1655 1656 /* Get things going. */ 1657 bptr = buffer; 1658 *bptr++ = (unsigned)val % BASE; /* type cast to unsigned, bug fix Wolf Lammen */ 1659 val = (unsigned)val / BASE; /* type cast to unsigned, bug fix Wolf Lammen */ 1660 1661 /* Extract remaining digits. */ 1662 while (val != 0) 1663 { 1664 *bptr++ = val % BASE; 1665 val = val / BASE; 1666 ix++; /* Count the digits. */ 1667 } 1668 1669 /* Make the number. */ 1670 bc_free_num (num); 1671 *num = bc_new_num (ix, 0); 1672 if (neg) (*num)->n_sign = MINUS; 1673 1674 /* Assign the digits. */ 1675 vptr = (*num)->n_value; 1676 while (ix-- > 0) 1677 *vptr++ = *--bptr; 1678 } 1679 1680 /* Convert a numbers to a string. Base 10 only.*/ 1681 1682 char 1683 *bc_num2str (num) 1684 bc_num num; 1685 { 1686 char *str, *sptr; 1687 char *nptr; 1688 int index, signch; 1689 1690 /* Allocate the string memory. */ 1691 signch = ( num->n_sign == PLUS ? 0 : 1 ); /* Number of sign chars. */ 1692 if (num->n_scale > 0) 1693 str = (char *) malloc (num->n_len + num->n_scale + 2 + signch); 1694 else 1695 str = (char *) malloc (num->n_len + 1 + signch); 1696 if (str == NULL) bc_out_of_memory(); 1697 1698 /* The negative sign if needed. */ 1699 sptr = str; 1700 if (signch) *sptr++ = '-'; 1701 1702 /* Load the whole number. */ 1703 nptr = num->n_value; 1704 for (index=num->n_len; index>0; index--) 1705 *sptr++ = BCD_CHAR(*nptr++); 1706 1707 /* Now the fraction. */ 1708 if (num->n_scale > 0) 1709 { 1710 *sptr++ = '.'; 1711 for (index=0; index<num->n_scale; index++) 1712 *sptr++ = BCD_CHAR(*nptr++); 1713 } 1714 1715 /* Terminate the string and return it! */ 1716 *sptr = '\0'; 1717 return (str); 1718 } 1719 /* Convert strings to bc numbers. Base 10 only.*/ 1720 1721 void 1722 bc_str2num (num, str, scale) 1723 bc_num *num; 1724 char *str; 1725 int scale; 1726 { 1727 int digits, strscale; 1728 char *ptr, *nptr; 1729 char zero_int; 1730 1731 /* Prepare num. */ 1732 bc_free_num (num); 1733 1734 /* Check for valid number and count digits. */ 1735 ptr = str; 1736 digits = 0; 1737 strscale = 0; 1738 zero_int = FALSE; 1739 if ( (*ptr == '+') || (*ptr == '-')) ptr++; /* Sign */ 1740 while (*ptr == '0') ptr++; /* Skip leading zeros. */ 1741 while (isdigit((int)*ptr)) ptr++, digits++; /* digits */ 1742 if (*ptr == '.') ptr++; /* decimal point */ 1743 while (isdigit((int)*ptr)) ptr++, strscale++; /* digits */ 1744 if ((*ptr != '\0') || (digits+strscale == 0)) 1745 { 1746 *num = bc_copy_num (_zero_); 1747 return; 1748 } 1749 1750 /* Adjust numbers and allocate storage and initialize fields. */ 1751 strscale = MIN(strscale, scale); 1752 if (digits == 0) 1753 { 1754 zero_int = TRUE; 1755 digits = 1; 1756 } 1757 *num = bc_new_num (digits, strscale); 1758 1759 /* Build the whole number. */ 1760 ptr = str; 1761 if (*ptr == '-') 1762 { 1763 (*num)->n_sign = MINUS; 1764 ptr++; 1765 } 1766 else 1767 { 1768 (*num)->n_sign = PLUS; 1769 if (*ptr == '+') ptr++; 1770 } 1771 while (*ptr == '0') ptr++; /* Skip leading zeros. */ 1772 nptr = (*num)->n_value; 1773 if (zero_int) 1774 { 1775 *nptr++ = 0; 1776 digits = 0; 1777 } 1778 for (;digits > 0; digits--) 1779 *nptr++ = CH_VAL(*ptr++); 1780 1781 1782 /* Build the fractional part. */ 1783 if (strscale > 0) 1784 { 1785 ptr++; /* skip the decimal point! */ 1786 for (;strscale > 0; strscale--) 1787 *nptr++ = CH_VAL(*ptr++); 1788 } 1789 } 1790 1791 /* pn prints the number NUM in base 10. */ 1792 1793 static void 1794 out_char (int c) 1795 { 1796 putchar(c); 1797 } 1798 1799 1800 void 1801 pn (num) 1802 bc_num num; 1803 { 1804 bc_out_num (num, 10, out_char, 0); 1805 out_char ('\n'); 1806 } 1807 1808 1809 /* pv prints a character array as if it was a string of bcd digits. */ 1810 void 1811 pv (name, num, len) 1812 char *name; 1813 unsigned char *num; 1814 int len; 1815 { 1816 int i; 1817 printf ("%s=", name); 1818 for (i=0; i<len; i++) printf ("%c",BCD_CHAR(num[i])); 1819 printf ("\n"); 1820 }