/* number.c: Implements arbitrary precision numbers. */
/* This file is part of bc written for MINIX.
Copyright (C) 1991 Free Software Foundation, Inc.
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License , or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
You may contact the author by:
us-mail: Philip A. Nelson
Computer Science Department, 9062
Western Washington University
Bellingham, WA 98226-9062
*************************************************************************/
/* Storage used for special numbers. */
/* "Frees" a bc_num NUM. Actually decreases reference count and only
frees the storage if reference count is zero. */
if (*num
== NULL
) return;
if ((*num
)->n_refs
== 0) free(*num
);
/* new_num allocates a number and sets fields to known values. */
temp
= (bc_num
) malloc (sizeof(bc_struct
)+length
+scale
);
if (temp
== NULL
) out_of_memory ();
/* Intitialize the number package! */
/* Make a copy of a number! Just increments the reference count! */
/* Initialize a number NUM by making it a copy of zero. */
*num
= copy_num (_zero_
);
/* Convert an integer VAL to a bc number NUM. */
/* Extract remaining digits. */
ix
++; /* Count the digits. */
if (neg
) (*num
)->n_sign
= MINUS
;
/* Convert a number NUM to a long. The function returns only the integer
part of the number. For numbers that are too large to represent as
a long, this function returns a zero. This can be detected by checking
the NUM for zero after having a zero returned. */
/* Extract the int value, ignore the fraction. */
for (index
=num
->n_len
; (index
>0) && (val
<=(LONG_MAX
/10)); index
--)
/* Check for overflow. If overflow, return zero. */
/* The following are some math routines for numbers. */
_PROTOTYPE(static int _do_compare
, (bc_num n1
, bc_num n2
, int use_sign
,
_PROTOTYPE(static void _rm_leading_zeros
, (bc_num num
));
_PROTOTYPE(static bc_num _do_add
, (bc_num n1
, bc_num n2
));
_PROTOTYPE(static bc_num _do_sub
, (bc_num n1
, bc_num n2
));
_PROTOTYPE(static void _one_mult
, (unsigned char *num
, int size
, int digit
,
/* Compare two bc numbers. Return value is 0 if equal, -1 if N1 is less
than N2 and +1 if N1 is greater than N2. If USE_SIGN is false, just
compare the magnitudes. */
_do_compare (n1
, n2
, use_sign
, ignore_last
)
/* First, compare signs. */
if (use_sign
&& n1
->n_sign
!= n2
->n_sign
)
return (1); /* Positive N1 > Negative N2 */
return (-1); /* Negative N1 < Positive N1 */
/* Now compare the magnitude. */
if (n1
->n_len
!= n2
->n_len
)
if (n1
->n_len
> n2
->n_len
)
/* Magnitude of n1 > n2. */
if (!use_sign
|| n1
->n_sign
== PLUS
)
/* Magnitude of n1 < n2. */
if (!use_sign
|| n1
->n_sign
== PLUS
)
/* If we get here, they have the same number of integer digits.
check the integer part and the equal length part of the fraction. */
count
= n1
->n_len
+ MIN (n1
->n_scale
, n2
->n_scale
);
while ((count
> 0) && (*n1ptr
== *n2ptr
))
if (ignore_last
&& count
== 1 && n1
->n_scale
== n2
->n_scale
)
/* Magnitude of n1 > n2. */
if (!use_sign
|| n1
->n_sign
== PLUS
)
/* Magnitude of n1 < n2. */
if (!use_sign
|| n1
->n_sign
== PLUS
)
/* They are equal up to the last part of the equal part of the fraction. */
if (n1
->n_scale
!= n2
->n_scale
)
if (n1
->n_scale
> n2
->n_scale
)
for (count
= n1
->n_scale
-n2
->n_scale
; count
>0; count
--)
/* Magnitude of n1 > n2. */
if (!use_sign
|| n1
->n_sign
== PLUS
)
for (count
= n2
->n_scale
-n1
->n_scale
; count
>0; count
--)
/* Magnitude of n1 < n2. */
if (!use_sign
|| n1
->n_sign
== PLUS
)
/* They must be equal! */
/* This is the "user callable" routine to compare numbers N1 and N2. */
return _do_compare (n1
, n2
, TRUE
, FALSE
);
/* In some places we need to check if the number NUM is zero. */
if (num
== _zero_
) return TRUE
;
count
= num
->n_len
+ num
->n_scale
;
while ((count
> 0) && (*nptr
++ == 0)) count
--;
/* In some places we need to check if the number is negative. */
return num
->n_sign
== MINUS
;
/* For many things, we may have leading zeros in a number NUM.
_rm_leading_zeros just moves the data to the correct
place and adjusts the length. */
/* Do a quick check to see if we need to do it. */
if (*num
->n_value
!= 0) return;
/* The first digit is 0, find the first non-zero digit in the 10's or
while (bytes
> 1 && *src
== 0) src
++, bytes
--;
while (bytes
-- > 0) *dst
++ = *src
++;
/* Perform addition: N1 is added to N2 and the value is
returned. The signs of N1 and N2 are ignored. */
int sum_scale
, sum_digits
;
char *n1ptr
, *n2ptr
, *sumptr
;
int carry
, n1bytes
, n2bytes
;
sum_scale
= MAX (n1
->n_scale
, n2
->n_scale
);
sum_digits
= MAX (n1
->n_len
, n2
->n_len
) + 1;
sum
= new_num (sum_digits
,sum_scale
);
/* Start with the fraction part. Initialize the pointers. */
n1ptr
= (char *) (n1
->n_value
+ n1
->n_len
+ n1bytes
- 1);
n2ptr
= (char *) (n2
->n_value
+ n2
->n_len
+ n2bytes
- 1);
sumptr
= (char *) (sum
->n_value
+ sum_scale
+ sum_digits
- 1);
/* Add the fraction part. First copy the longer fraction.*/
{ *sumptr
-- = *n1ptr
--; n1bytes
--;}
{ *sumptr
-- = *n2ptr
--; n2bytes
--;}
/* Now add the remaining fraction part and equal size integer parts. */
while ((n1bytes
> 0) && (n2bytes
> 0))
*sumptr
= *n1ptr
-- + *n2ptr
-- + carry
;
/* Now add carry the longer integer part. */
{ n1bytes
= n2bytes
; n1ptr
= n2ptr
; }
*sumptr
= *n1ptr
-- + carry
;
/* Adjust sum and return. */
/* Perform subtraction: N2 is subtracted from N1 and the value is
returned. The signs of N1 and N2 are ignored. Also, N1 is
assumed to be larger than N2. */
int diff_scale
, diff_len
;
char *n1ptr
, *n2ptr
, *diffptr
;
/* Allocate temporary storage. */
diff_len
= MAX (n1
->n_len
, n2
->n_len
);
diff_scale
= MAX (n1
->n_scale
, n2
->n_scale
);
min_len
= MIN (n1
->n_len
, n2
->n_len
);
min_scale
= MIN (n1
->n_scale
, n2
->n_scale
);
diff
= new_num (diff_len
, diff_scale
);
/* Initialize the subtract. */
n1ptr
= (char *) (n1
->n_value
+ n1
->n_len
+ n1
->n_scale
-1);
n2ptr
= (char *) (n2
->n_value
+ n2
->n_len
+ n2
->n_scale
-1);
diffptr
= (char *) (diff
->n_value
+ diff_len
+ diff_scale
-1);
/* Subtract the numbers. */
/* Take care of the longer scaled number. */
if (n1
->n_scale
!= min_scale
)
/* n1 has the longer scale */
for (count
= n1
->n_scale
- min_scale
; count
> 0; count
--)
/* n2 has the longer scale */
for (count
= n2
->n_scale
- min_scale
; count
> 0; count
--)
val
= - *n2ptr
-- - borrow
;
/* Now do the equal length scale and integer parts. */
for (count
= 0; count
< min_len
+ min_scale
; count
++)
val
= *n1ptr
-- - *n2ptr
-- - borrow
;
/* If n1 has more digits then n2, we now do that subtract. */
for (count
= diff_len
- min_len
; count
> 0; count
--)
/* Clean up and return. */
_rm_leading_zeros (diff
);
/* Here is the full add routine that takes care of negative numbers.
N1 is added to N2 and the result placed into RESULT. */
if (n1
->n_sign
== n2
->n_sign
)
sum
->n_sign
= n1
->n_sign
;
/* subtraction must be done. */
cmp_res
= _do_compare (n1
, n2
, FALSE
, FALSE
); /* Compare magnitudes. */
/* n1 is less than n2, subtract n1 from n2. */
sum
->n_sign
= n2
->n_sign
;
/* They are equal! return zero! */
/* n2 is less than n1, subtract n2 from n1. */
sum
->n_sign
= n1
->n_sign
;
/* Clean up and return. */
/* Here is the full subtract routine that takes care of negative numbers.
N2 is subtracted from N1 and the result placed in RESULT. */
if (n1
->n_sign
!= n2
->n_sign
)
diff
->n_sign
= n1
->n_sign
;
/* subtraction must be done. */
cmp_res
= _do_compare (n1
, n2
, FALSE
, FALSE
); /* Compare magnitudes. */
/* n1 is less than n2, subtract n1 from n2. */
diff
->n_sign
= (n2
->n_sign
== PLUS
? MINUS
: PLUS
);
/* They are equal! return zero! */
diff
= copy_num (_zero_
);
/* n2 is less than n1, subtract n2 from n1. */
diff
->n_sign
= n1
->n_sign
;
/* Clean up and return. */
/* The multiply routine. N2 time N1 is put int PROD with the scale of
the result being MIN(N2 scale+N1 scale, MAX (SCALE, N2 scale, N1 scale)).
bc_multiply (n1
, n2
, prod
, scale
)
bc_num pval
; /* For the working storage. */
char *n1ptr
, *n2ptr
, *pvptr
; /* Work pointers. */
char *n1end
, *n2end
; /* To the end of n1 and n2. */
int len1
, len2
, total_digits
;
int full_scale
, prod_scale
;
len1
= n1
->n_len
+ n1
->n_scale
;
len2
= n2
->n_len
+ n2
->n_scale
;
total_digits
= len1
+ len2
;
full_scale
= n1
->n_scale
+ n2
->n_scale
;
prod_scale
= MIN(full_scale
,MAX(scale
,MAX(n1
->n_scale
,n2
->n_scale
)));
toss
= full_scale
- prod_scale
;
pval
= new_num (total_digits
-full_scale
, prod_scale
);
pval
->n_sign
= ( n1
->n_sign
== n2
->n_sign
? PLUS
: MINUS
);
n1end
= (char *) (n1
->n_value
+ len1
- 1);
n2end
= (char *) (n2
->n_value
+ len2
- 1);
pvptr
= (char *) (pval
->n_value
+ total_digits
- toss
- 1);
/* Here are the loops... */
for (indx
= 0; indx
< toss
; indx
++)
n1ptr
= (char *) (n1end
- MAX(0, indx
-len2
+1));
n2ptr
= (char *) (n2end
- MIN(indx
, len2
-1));
while ((n1ptr
>= n1
->n_value
) && (n2ptr
<= n2end
))
sum
+= *n1ptr
-- * *n2ptr
++;
for ( ; indx
< total_digits
-1; indx
++)
n1ptr
= (char *) (n1end
- MAX(0, indx
-len2
+1));
n2ptr
= (char *) (n2end
- MIN(indx
, len2
-1));
while ((n1ptr
>= n1
->n_value
) && (n2ptr
<= n2end
))
sum
+= *n1ptr
-- * *n2ptr
++;
/* Assign to prod and clean up the number. */
_rm_leading_zeros (*prod
);
/* Some utility routines for the divide: First a one digit multiply.
NUM (with SIZE digits) is multiplied by DIGIT and the result is
placed into RESULT. It is written so that NUM and RESULT can be
_one_mult (num
, size
, digit
, result
)
unsigned char *nptr
, *rptr
;
memset (result
, 0, size
);
memcpy (result
, num
, size
);
nptr
= (unsigned char *) (num
+size
-1);
rptr
= (unsigned char *) (result
+size
-1);
value
= *nptr
-- * digit
+ carry
;
if (carry
!= 0) *rptr
= carry
;
/* The full division routine. This computes N1 / N2. It returns
0 if the division is ok and the result is in QUOT. The number of
digits after the decimal point is SCALE. It returns -1 if division
by zero is tried. The algorithm is found in Knuth Vol 2. p237. */
bc_divide (n1
, n2
, quot
, scale
)
unsigned char *num1
, *num2
;
unsigned char *ptr1
, *ptr2
, *n2ptr
, *qptr
;
unsigned int len1
, len2
, scale2
, qdigits
, extra
, count
;
unsigned int qdig
, qguess
, borrow
, carry
;
/* Test for divide by zero. */
if (is_zero (n2
)) return -1;
/* Test for divide by 1. If it is we must truncate. */
if (n2
->n_len
== 1 && *n2
->n_value
== 1)
qval
= new_num (n1
->n_len
, scale
);
qval
->n_sign
= (n1
->n_sign
== n2
->n_sign
? PLUS
: MINUS
);
memset (&qval
->n_value
[n1
->n_len
],0,scale
);
memcpy (qval
->n_value
, n1
->n_value
,
n1
->n_len
+ MIN(n1
->n_scale
,scale
));
/* Set up the divide. Move the decimal point on n1 by n2's scale.
Remember, zeros on the end of num2 are wasted effort for dividing. */
n2ptr
= (unsigned char *) n2
->n_value
+n2
->n_len
+scale2
-1;
while ((scale2
> 0) && (*n2ptr
-- == 0)) scale2
--;
len1
= n1
->n_len
+ scale2
;
scale1
= n1
->n_scale
- scale2
;
num1
= (unsigned char *) malloc (n1
->n_len
+n1
->n_scale
+extra
+2);
if (num1
== NULL
) out_of_memory();
memset (num1
, 0, n1
->n_len
+n1
->n_scale
+extra
+2);
memcpy (num1
+1, n1
->n_value
, n1
->n_len
+n1
->n_scale
);
len2
= n2
->n_len
+ scale2
;
num2
= (unsigned char *) malloc (len2
+1);
if (num2
== NULL
) out_of_memory();
memcpy (num2
, n2
->n_value
, len2
);
/* Calculate the number of quotient digits. */
qdigits
= scale
+1; /* One for the zero integer part. */
qdigits
= len1
-len2
+scale
+1;
/* Allocate and zero the storage for the quotient. */
qval
= new_num (qdigits
-scale
,scale
);
memset (qval
->n_value
, 0, qdigits
);
/* Allocate storage for the temporary storage mval. */
mval
= (unsigned char *) malloc (len2
+1);
if (mval
== NULL
) out_of_memory ();
/* Now for the full divide algorithm. */
norm
= 10 / ((int)*n2ptr
+ 1);
_one_mult (num1
, len1
+scale1
+extra
+1, norm
, num1
);
_one_mult (n2ptr
, len2
, norm
, n2ptr
);
/* Initialize divide loop. */
qptr
= (unsigned char *) qval
->n_value
+len2
-len1
;
qptr
= (unsigned char *) qval
->n_value
;
while (qdig
<= len1
+scale
-len2
)
/* Calculate the quotient digit guess. */
if (*n2ptr
== num1
[qdig
])
qguess
= (num1
[qdig
]*10 + num1
[qdig
+1]) / *n2ptr
;
(num1
[qdig
]*10 + num1
[qdig
+1] - *n2ptr
*qguess
)*10
(num1
[qdig
]*10 + num1
[qdig
+1] - *n2ptr
*qguess
)*10
/* Multiply and subtract. */
_one_mult (n2ptr
, len2
, qguess
, mval
+1);
ptr1
= (unsigned char *) num1
+qdig
+len2
;
ptr2
= (unsigned char *) mval
+len2
;
for (count
= 0; count
< len2
+1; count
++)
val
= (int) *ptr1
- (int) *ptr2
-- - borrow
;
/* Test for negative result. */
ptr1
= (unsigned char *) num1
+qdig
+len2
;
ptr2
= (unsigned char *) n2ptr
+len2
-1;
for (count
= 0; count
< len2
; count
++)
val
= (int) *ptr1
+ (int) *ptr2
-- + carry
;
if (carry
== 1) *ptr1
= (*ptr1
+ 1) % 10;
/* We now know the quotient digit. */
/* Clean up and return the number. */
qval
->n_sign
= ( n1
->n_sign
== n2
->n_sign
? PLUS
: MINUS
);
if (is_zero (qval
)) qval
->n_sign
= PLUS
;
_rm_leading_zeros (qval
);
/* Clean up temporary storage. */
return 0; /* Everything is OK. */
/* Modulo for numbers. This computes NUM1 % NUM2 and puts the
bc_modulo (num1
, num2
, result
, scale
)
bc_num num1
, num2
, *result
;
/* Check for correct numbers. */
if (is_zero (num2
)) return -1;
/* Calculate final scale. */
rscale
= MAX (num1
->n_scale
, num2
->n_scale
+scale
);
bc_divide (num1
, num2
, &temp
, scale
);
bc_multiply (temp
, num2
, &temp
, rscale
);
bc_sub (num1
, temp
, result
);
return 0; /* Everything is OK. */
/* Raise NUM1 to the NUM2 power. The result is placed in RESULT.
Maximum exponent is LONG_MAX. If a NUM2 is not an integer,
only the integer part is used. */
bc_raise (num1
, num2
, result
, scale
)
bc_num num1
, num2
, *result
;
/* Check the exponent for scale digits and convert to a long. */
rt_warn ("non-zero scale in exponent");
exponent
= num2long (num2
);
if (exponent
== 0 && (num2
->n_len
> 1 || num2
->n_value
[0] != 0))
rt_error ("exponent too large in raise");
/* Special case if exponent is a zero. */
*result
= copy_num (_one_
);
/* Other initializations. */
rscale
= MIN (num1
->n_scale
*exponent
, MAX(scale
, num1
->n_scale
));
/* Do the calculation. */
bc_multiply (temp
, power
, &temp
, rscale
);
bc_multiply (power
, power
, &power
, rscale
);
exponent
= exponent
>> 1;
bc_divide (_one_
, temp
, result
, rscale
);
/* Take the square root NUM and return it in NUM with SCALE digits
after the decimal place. */
int rscale
, cmp_res
, done
;
bc_num guess
, guess1
, point5
;
cmp_res
= bc_compare (*num
, _zero_
);
*num
= copy_num (_zero_
);
cmp_res
= bc_compare (*num
, _one_
);
/* Initialize the variables. */
rscale
= MAX (scale
, (*num
)->n_scale
);
/* Calculate the initial guess. */
/* The number is between 0 and 1. Guess should start at 1. */
guess
= copy_num (_one_
);
/* The number is greater than 1. Guess should start at 10^(exp/2). */
int2num (&guess1
,(*num
)->n_len
);
bc_multiply (guess1
, point5
, &guess1
, rscale
);
bc_raise (guess
, guess1
, &guess
, rscale
);
/* Find the square root using Newton's algorithm. */
guess1
= copy_num (guess
);
bc_divide (*num
,guess
,&guess
,cscale
);
bc_add (guess
,guess1
,&guess
);
bc_multiply (guess
,point5
,&guess
,cscale
);
cmp_res
= _do_compare (guess
,guess1
,FALSE
,TRUE
);
if (cmp_res
== 0) done
= TRUE
;
/* Assign the number and clean up. */
bc_divide (guess
,_one_
,num
,rscale
);
/* The following routines provide output for bcd numbers package
using the rules of POSIX bc for output. */
/* This structure is used for saving digits in the conversion process. */
/* The reference string for digits. */
char ref_str
[] = "0123456789ABCDEF";
/* A special output routine for "multi-character digits." Exactly
SIZE characters must be output for the value VAL. If SPACE is
non-zero, we must output one space before the number. OUT_CHAR
is the actual routine for writing the characters. */
out_long (val
, size
, space
, out_char
)
if (space
) (*out_char
) (' ');
sprintf (digits
, "%ld", val
);
for (ix
=0; ix
< len
; ix
++)
(*out_char
) (digits
[ix
]);
/* Output of a bcd number. NUM is written in base O_BASE using OUT_CHAR
as the routine to do the actual output of the characters. */
out_num (num
, o_base
, out_char
)
int index
, fdigit
, pre_space
;
bc_num int_part
, frac_part
, base
, cur_dig
, t_num
;
/* The negative sign if needed. */
if (num
->n_sign
== MINUS
) (*out_char
) ('-');
/* The number is in base 10, do it the fast way. */
if (num
->n_len
> 1 || *nptr
!= 0)
for (index
=num
->n_len
; index
>0; index
--)
(*out_char
) (BCD_CHAR(*nptr
++));
for (index
=0; index
<num
->n_scale
; index
++)
(*out_char
) (BCD_CHAR(*nptr
++));
/* The number is some other base. */
bc_divide (num
, _one_
, &int_part
, 0);
bc_sub (num
, int_part
, &frac_part
);
/* Get the digits of the integer part and push them on a stack. */
while (!is_zero (int_part
))
bc_modulo (int_part
, base
, &cur_dig
, 0);
temp
= (stk_rec
*) malloc (sizeof(stk_rec
));
if (temp
== NULL
) out_of_memory();
temp
->digit
= num2long (cur_dig
);
bc_divide (int_part
, base
, &int_part
, 0);
/* Print the digits on the stack. */
(*out_char
) (ref_str
[ (int) temp
->digit
]);
out_long (temp
->digit
, base
->n_len
, 1, out_char
);
/* Get and print the digits of the fraction part. */
t_num
= copy_num (_one_
);
while (t_num
->n_len
<= num
->n_scale
) {
bc_multiply (frac_part
, base
, &frac_part
, num
->n_scale
);
fdigit
= num2long (frac_part
);
int2num (&int_part
, fdigit
);
bc_sub (frac_part
, int_part
, &frac_part
);
(*out_char
) (ref_str
[fdigit
]);
out_long (fdigit
, base
->n_len
, pre_space
, out_char
);
bc_multiply (t_num
, base
, &t_num
, 0);
/* Debugging procedures. Some are just so one can call them from the
/* p_n prints the number NUM in base 10. */
out_num (num
, 10, out_char
);
/* p_b prints a character array as if it was a string of bcd digits. */
for (i
=0; i
<len
; i
++) printf ("%c",BCD_CHAR(num
[i
]));
/* Convert strings to bc numbers. Base 10 only.*/
str2num (num
, str
, scale
)
/* Check for valid number and count digits. */
if ( (*ptr
== '+') || (*ptr
== '-')) ptr
++; /* Sign */
while (*ptr
== '0') ptr
++; /* Skip leading zeros. */
while (isdigit(*ptr
)) ptr
++, digits
++; /* digits */
if (*ptr
== '.') ptr
++; /* decimal point */
while (isdigit(*ptr
)) ptr
++, strscale
++; /* digits */
if ((*ptr
!= '\0') || (digits
+strscale
== 0))
*num
= copy_num (_zero_
);
/* Adjust numbers and allocate storage and initialize fields. */
strscale
= MIN(strscale
, scale
);
*num
= new_num (digits
, strscale
);
/* Build the whole number. */
while (*ptr
== '0') ptr
++; /* Skip leading zeros. */
for (;digits
> 0; digits
--)
*nptr
++ = CH_VAL(*ptr
++);
/* Build the fractional part. */
ptr
++; /* skip the decimal point! */
for (;strscale
> 0; strscale
--)
*nptr
++ = CH_VAL(*ptr
++);
/* Convert a numbers to a string. Base 10 only.*/
/* Allocate the string memory. */
signch
= ( num
->n_sign
== PLUS
? 0 : 1 ); /* Number of sign chars. */
str
= (char *) malloc (num
->n_len
+ num
->n_scale
+ 2 + signch
);
str
= (char *) malloc (num
->n_len
+ 1 + signch
);
if (str
== NULL
) out_of_memory();
/* The negative sign if needed. */
if (signch
) *sptr
++ = '-';
/* Load the whole number. */
for (index
=num
->n_len
; index
>0; index
--)
*sptr
++ = BCD_CHAR(*nptr
++);
for (index
=0; index
<num
->n_scale
; index
++)
*sptr
++ = BCD_CHAR(*nptr
++);
/* Terminate the string and return it! */