d025cfd6c5
Change a round to a truncate. Problem reported from Christoph Kukulies: 9.8 2 / p did an IOT trap. There is one remaining problem.: 9.8 2 % p shows on other machines 1.8 but does here in the moment 1.
1236 lines
29 KiB
C
1236 lines
29 KiB
C
/*
|
||
* Arbitrary precision decimal arithmetic.
|
||
*
|
||
* Copyright (C) 1984 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, 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; if not, you can either send email to this
|
||
* program's author (see below) or write to: The Free Software Foundation,
|
||
* Inc.; 675 Mass Ave. Cambridge, MA 02139, USA.
|
||
*/
|
||
|
||
/* Some known problems:
|
||
|
||
Another problem with decimal_div is found when you try to
|
||
divide a number with > scale fraction digits by 1. The
|
||
expected result is simply truncation, but all sorts of things
|
||
happen instead. An example is that the result of .99999998/1
|
||
with scale set to 6 is .000001
|
||
|
||
There are some problems in the behavior of the decimal package
|
||
related to printing and parsing. The
|
||
printer is weird about very large output radices, tending to want
|
||
to output single ASCII characters for any and all digits (even
|
||
in radices > 127). The UNIX bc approach is to print digit groups
|
||
separated by spaces. There is a rather overwrought workaround in
|
||
the function decputc() in bcmisc.c, but it would be better if
|
||
decimal.c got a fix for this. */
|
||
|
||
/* For stand-alone testing, compile with -DTEST.
|
||
This DTESTable feature defines a `main' function
|
||
which is a simple loop that accepts input of the form
|
||
number space op space number newline
|
||
where op is +, -, *, /, %, p or r,
|
||
and performs the operation and prints the operands and result.
|
||
`p' means print the first number in the radix spec'd by the second.
|
||
`r' means read the first one in the radix specified by the second
|
||
(and print the result in decimal).
|
||
Divide in this test keeps three fraction digits. */
|
||
|
||
#include "decimal.h"
|
||
|
||
#define MAX(a, b) (((a) > (b) ? (a) : (b)))
|
||
|
||
/* Some constant decimal numbers */
|
||
|
||
struct decimal decimal_zero = {0, 0, 0, 0, 0};
|
||
|
||
struct decimal decimal_one = {0, 0, 1, 0, 1};
|
||
|
||
/*** Assumes RADIX is even ***/
|
||
struct decimal decimal_half = {0, 1, 0, 0, RADIX / 2};
|
||
|
||
decimal static decimal_add1 (), decimal_sub1 ();
|
||
static void add_scaled ();
|
||
static int subtract_scaled ();
|
||
|
||
/* Create and return a decimal number that has `before' digits before
|
||
the decimal point and `after' digits after. The digits themselves are
|
||
initialized to zero. */
|
||
|
||
decimal
|
||
make_decimal (before, after)
|
||
int before, after;
|
||
{
|
||
decimal result;
|
||
if (before >= 1<<16)
|
||
{
|
||
decimal_error ("%d too many decimal digits", before);
|
||
return 0;
|
||
}
|
||
if (after >= 1<<15)
|
||
{
|
||
decimal_error ("%d too many decimal digits", after);
|
||
return 0;
|
||
}
|
||
result = (decimal) malloc (sizeof (struct decimal) + before + after - 1);
|
||
result->sign = 0;
|
||
result->before = before;
|
||
result->after = after;
|
||
result->refcnt = 0;
|
||
bzero (result->contents, before + after);
|
||
return result;
|
||
}
|
||
|
||
/* Create a copy of the decimal number `b' and return it. */
|
||
|
||
decimal
|
||
decimal_copy (b)
|
||
decimal b;
|
||
{
|
||
decimal result = make_decimal (b->before, b->after);
|
||
bcopy (b->contents, result->contents, LENGTH(b));
|
||
result->sign = b->sign;
|
||
return result;
|
||
}
|
||
|
||
/* Copy a decimal number `b' but extend or truncate to exactly
|
||
`digits' fraction digits. */
|
||
|
||
static decimal
|
||
decimal_copy_1 (b, digits)
|
||
decimal b;
|
||
int digits;
|
||
{
|
||
if (digits > b->after)
|
||
{
|
||
decimal result = make_decimal (b->before, digits);
|
||
bcopy (b->contents, result->contents + (digits - (int) b->after), LENGTH(b));
|
||
return result;
|
||
}
|
||
else
|
||
return decimal_trunc_digits (b, digits);
|
||
}
|
||
|
||
/* flush specified number `digits' of trailing fraction digits,
|
||
and flush any trailing fraction zero digits exposed after they are gone.
|
||
The number `b' is actually modified; no new storage is allocated.
|
||
That is why this is not global. */
|
||
|
||
static void
|
||
flush_trailing_digits (b, digits)
|
||
decimal b;
|
||
int digits;
|
||
{
|
||
int flush = digits;
|
||
int maxdig = b->after;
|
||
|
||
while (flush < maxdig && !b->contents [flush])
|
||
flush++;
|
||
|
||
if (flush)
|
||
{
|
||
int i;
|
||
|
||
b->after -= flush;
|
||
for (i = 0; i < LENGTH (b); i++)
|
||
b->contents[i] = b->contents[flush + i];
|
||
}
|
||
|
||
}
|
||
|
||
/* Return nonzero integer if the value of decimal number `b' is zero. */
|
||
|
||
int
|
||
decimal_zerop (b)
|
||
decimal b;
|
||
{
|
||
return !LENGTH(b);
|
||
}
|
||
|
||
/* Compare two decimal numbers arithmetically.
|
||
The value is < 0 if b1 < b2, > 0 if b1 > b2, 0 if b1 = b2.
|
||
This is the same way that `strcmp' reports the result of comparing
|
||
strings. */
|
||
|
||
int
|
||
decimal_compare (b1, b2)
|
||
decimal b1, b2;
|
||
{
|
||
int l1, l2;
|
||
char *p1, *p2, *s1, *s2;
|
||
int i;
|
||
|
||
/* If signs differ, deduce result from the signs */
|
||
|
||
if (b2->sign && !b1->sign) return 1;
|
||
if (b1->sign && !b2->sign) return -1;
|
||
|
||
/* If same sign but number of nonfraction digits differs,
|
||
the one with more of them is farther from zero. */
|
||
|
||
if (b1->before != b2->before)
|
||
if (b1->sign)
|
||
return (int) (b2->before - b1->before);
|
||
else
|
||
return (int) (b1->before - b2->before);
|
||
|
||
/* Else compare the numbers digit by digit from high end */
|
||
l1 = LENGTH(b1);
|
||
l2 = LENGTH(b2);
|
||
s1 = b1->contents; /* Start of number -- don't back up digit pointer past here */
|
||
s2 = b2->contents;
|
||
p1 = b1->contents + l1; /* Scanning pointer, for fetching digits. */
|
||
p2 = b2->contents + l2;
|
||
for (i = MAX(l1, l2); i >= 0; i--)
|
||
{
|
||
int r = ((p1 != s1) ? *--p1 : 0) - ((p2 != s2) ? *--p2 : 0);
|
||
if (r)
|
||
return b1->sign ? -r : r;
|
||
}
|
||
return 0;
|
||
}
|
||
|
||
/* Return the number of digits stored in decimal number `b' */
|
||
|
||
int
|
||
decimal_length (b)
|
||
decimal b;
|
||
{
|
||
return LENGTH(b);
|
||
}
|
||
|
||
/* Return the number of fraction digits stored in decimal number `b'. */
|
||
|
||
int
|
||
decimal_after (b)
|
||
decimal b;
|
||
{
|
||
return b->after;
|
||
}
|
||
|
||
/* Round decimal number `b' to have only `digits' fraction digits.
|
||
Result is rounded to nearest unit in the last remaining digit.
|
||
Return the result, another decimal number. */
|
||
|
||
decimal
|
||
decimal_round_digits (b, digits)
|
||
decimal b;
|
||
int digits;
|
||
{
|
||
decimal result;
|
||
int old;
|
||
|
||
if (b->after <= digits) return decimal_copy (b);
|
||
|
||
if (digits < 0)
|
||
{
|
||
decimal_error ("request to keep negative number of digits %d", digits);
|
||
return decimal_copy (b);
|
||
}
|
||
|
||
result = make_decimal (b->before + 1, b->after);
|
||
result->sign = b->sign;
|
||
bcopy (b->contents, result->contents, LENGTH(b));
|
||
|
||
old = result->after;
|
||
|
||
/* Add .5 * last place to keep, so that we round rather than truncate */
|
||
/* Note this ignores sign of result, so if result is negative
|
||
it is subtracting */
|
||
|
||
add_scaled (result, DECIMAL_HALF, 1, old - digits - 1);
|
||
|
||
/* Flush desired digits, and any trailing zeros exposed by them. */
|
||
|
||
flush_trailing_digits (result, old - digits);
|
||
|
||
/* Flush leading digits -- always is one, unless was a carry into it */
|
||
|
||
while (result->before > 0
|
||
&& result->contents[LENGTH(result) - 1] == 0)
|
||
result->before--;
|
||
|
||
return result;
|
||
}
|
||
|
||
/* Truncate decimal number `b' to have only `digits' fraction digits.
|
||
Any fraction digits in `b' beyond that are dropped and ignored.
|
||
Truncation is toward zero.
|
||
Return the result, another decimal number. */
|
||
|
||
decimal
|
||
decimal_trunc_digits (b, digits)
|
||
decimal b;
|
||
int digits;
|
||
{
|
||
decimal result = decimal_copy (b);
|
||
int old = result->after;
|
||
|
||
if (old <= digits) return result;
|
||
|
||
if (digits < 0)
|
||
{
|
||
decimal_error ("request to keep negative number of digits %d", digits);
|
||
return result;
|
||
}
|
||
|
||
flush_trailing_digits (result, old - digits);
|
||
|
||
return result;
|
||
}
|
||
|
||
/* Return the fractional part of decimal number `b':
|
||
that is, `b' - decimal_trunc_digits (`b') */
|
||
|
||
decimal
|
||
decimal_fraction (b)
|
||
decimal b;
|
||
{
|
||
decimal result = make_decimal (0, b->after);
|
||
bcopy (b->contents, result->contents, b->after);
|
||
return result;
|
||
}
|
||
|
||
/* return an integer whose value is that of decimal `b', sans its fraction. */
|
||
|
||
int
|
||
decimal_to_int (b)
|
||
decimal b;
|
||
{
|
||
int result = 0;
|
||
int i;
|
||
int end = b->after;
|
||
|
||
for (i = LENGTH(b) - 1; i >= end; i--)
|
||
{
|
||
result *= RADIX;
|
||
result += b->contents[i];
|
||
}
|
||
return result;
|
||
}
|
||
|
||
/* return a decimal whose value is the integer i. */
|
||
|
||
decimal
|
||
decimal_from_int (i)
|
||
int i;
|
||
{
|
||
int log, tem;
|
||
decimal result;
|
||
|
||
for (log = 0, tem = (i > 0 ? i : - i); tem; log++, tem /= RADIX);
|
||
|
||
result = make_decimal (log, 0);
|
||
|
||
for (log = 0, tem = (i > 0 ? i : - i); tem; log++, tem /= RADIX)
|
||
result->contents[log] = tem % RADIX;
|
||
|
||
if (i < 0) result->sign = 1;
|
||
return result;
|
||
}
|
||
|
||
/* Return (as an integer) the result of dividing decimal number `b' by
|
||
integer `divisor'.
|
||
This is used in printing decimal numbers in other radices. */
|
||
|
||
int
|
||
decimal_int_rem (b, divisor)
|
||
decimal b;
|
||
int divisor;
|
||
{
|
||
int len = LENGTH(b);
|
||
int end = b->after;
|
||
int accum = 0;
|
||
int i;
|
||
|
||
for (i = len - 1; i >= end; i--)
|
||
{
|
||
accum %= divisor;
|
||
accum *= RADIX;
|
||
accum += b->contents[i];
|
||
}
|
||
return accum % divisor;
|
||
}
|
||
|
||
/* Convert digit `digit' to a character and output it by calling
|
||
`charout' with it as arg. */
|
||
|
||
static void
|
||
print_digit (digit, charout)
|
||
int digit;
|
||
void (*charout) ();
|
||
{
|
||
if (digit < 10)
|
||
charout ('0' + digit);
|
||
else
|
||
charout ('A' + digit - 10);
|
||
}
|
||
|
||
/* print decimal number `b' in radix `radix', assuming it is an integer.
|
||
`r' is `radix' expressed as a decimal number. */
|
||
|
||
static
|
||
decimal_print_1 (b, r, radix, charout)
|
||
decimal b, r;
|
||
int radix;
|
||
void (*charout) ();
|
||
{
|
||
int digit = decimal_int_rem (b, radix);
|
||
decimal rest = decimal_div (b, r, 0);
|
||
|
||
if (!decimal_zerop (rest))
|
||
decimal_print_1 (rest, r, radix, charout);
|
||
|
||
print_digit (digit, charout);
|
||
|
||
free (rest);
|
||
}
|
||
|
||
/* User entry: print decimal number `b' in radix `radix' (an integer),
|
||
outputting characters by calling `charout'. */
|
||
|
||
void
|
||
decimal_print (b, charout, radix)
|
||
decimal b;
|
||
void (*charout) ();
|
||
int radix;
|
||
{
|
||
if (b->sign) charout ('-');
|
||
|
||
if (radix == RADIX)
|
||
{
|
||
/* decimal output => just print the digits, inserting a point in
|
||
the proper place. */
|
||
int i;
|
||
int before = b->before;
|
||
int len = before + b->after;
|
||
for (i = 0; i < len; i++)
|
||
{
|
||
if (i == before) charout ('.');
|
||
/* Broken if RADIX /= 10
|
||
charout ('0' + b->contents [len - 1 - i]); */
|
||
print_digit (b->contents [len - 1 - i], charout);
|
||
}
|
||
if (!len)
|
||
charout ('0');
|
||
}
|
||
else
|
||
{
|
||
/* nonstandard radix: must use multiply and divide to determine the
|
||
digits of the number in that radix. */
|
||
|
||
int i;
|
||
extern double log10 ();
|
||
/* Compute the number of fraction digits we want to have in the
|
||
new radix. They should contain the same amount of
|
||
information as the decimal digits we have. */
|
||
int nfrac = (b->after / log10 ((double) radix) + .99);
|
||
decimal r = decimal_from_int (radix);
|
||
decimal intpart = decimal_trunc_digits (b, 0);
|
||
|
||
/* print integer part */
|
||
decimal_print_1 (intpart, r, radix, charout);
|
||
free (intpart);
|
||
|
||
/* print fraction part */
|
||
if (nfrac)
|
||
{
|
||
decimal tem1, tem2;
|
||
tem1 = decimal_fraction (b);
|
||
charout ('.');
|
||
/* repeatedly multiply by `radix', print integer part as one digit,
|
||
and flush the integer part. */
|
||
for (i = 0; i < nfrac; i++)
|
||
{
|
||
tem2 = decimal_mul (tem1, r);
|
||
free (tem1);
|
||
print_digit (decimal_to_int (tem2), charout);
|
||
tem1 = decimal_fraction (tem2);
|
||
free (tem2);
|
||
}
|
||
free (tem1);
|
||
}
|
||
free (r);
|
||
}
|
||
}
|
||
|
||
static int
|
||
decode_digit (digitchar)
|
||
char digitchar;
|
||
{
|
||
if ('0' <= digitchar && digitchar <= '9')
|
||
return digitchar - '0';
|
||
if ('a' <= digitchar && digitchar <= 'z')
|
||
return digitchar - 'a' + 10;
|
||
if ('A' <= digitchar && digitchar <= 'Z')
|
||
return digitchar - 'A' + 10;
|
||
return -1;
|
||
}
|
||
|
||
/* Parse string `s' into a number using radix `radix'
|
||
and return result as a decimal number. */
|
||
|
||
decimal
|
||
decimal_parse (s, radix)
|
||
char *s;
|
||
int radix;
|
||
{
|
||
int i, len, before = -1;
|
||
char *p;
|
||
char c;
|
||
decimal result;
|
||
int negative = 0;
|
||
int excess_digit = 0;
|
||
|
||
if (*s == '-')
|
||
{
|
||
s++;
|
||
negative = 1;
|
||
}
|
||
|
||
/* First scan for valid characters.
|
||
Count total num digits, and count num before the decimal point. */
|
||
|
||
p = s;
|
||
i = 0;
|
||
while (c = *p++)
|
||
{
|
||
if (c == '.')
|
||
{
|
||
if (before >= 0)
|
||
decimal_error ("two decimal points in %s", s);
|
||
before = i;
|
||
}
|
||
else if (c == '0' && !i && before < 0)
|
||
s++; /* Discard leading zeros */
|
||
else if (decode_digit (c) >= 0)
|
||
{
|
||
i++;
|
||
if (decode_digit (c) > RADIX)
|
||
excess_digit = 1;
|
||
}
|
||
else
|
||
decimal_error ("invalid number %s", s);
|
||
}
|
||
|
||
len = i;
|
||
if (before < 0) before = i;
|
||
|
||
p = s;
|
||
|
||
/* Now parse those digits */
|
||
|
||
if (radix != RADIX || excess_digit)
|
||
{
|
||
decimal r = decimal_from_int (radix);
|
||
extern double log10 ();
|
||
int digits = (len - before) * log10 ((double) radix) + .99;
|
||
result = decimal_copy (DECIMAL_ZERO);
|
||
|
||
/* Parse all the digits into an integer, ignoring decimal point,
|
||
by multiplying by `radix'. */
|
||
|
||
while (i > 0 && (c = *p++))
|
||
{
|
||
if (c != '.')
|
||
{
|
||
decimal newdig = decimal_from_int (decode_digit (c));
|
||
decimal prod = decimal_mul (result, r);
|
||
decimal newresult = decimal_add (newdig, prod);
|
||
|
||
free (newdig); free (prod); free (result);
|
||
result = newresult;
|
||
i--;
|
||
}
|
||
}
|
||
|
||
/* Now put decimal point in right place
|
||
by dividing by `radix' once for each digit
|
||
that really should have followed the decimal point. */
|
||
|
||
for (i = before; i < len; i++)
|
||
{
|
||
decimal newresult = decimal_div (result, r, digits);
|
||
free (result);
|
||
result = newresult;
|
||
}
|
||
free (r);
|
||
}
|
||
else
|
||
{
|
||
/* radix is standard - just copy the digits into a decimal number. */
|
||
|
||
int tem;
|
||
result = make_decimal (before, len - before);
|
||
|
||
while (i > 0 && (c = *p++))
|
||
{
|
||
if ((c != '.') &&
|
||
((tem = decode_digit (c)) >= 0))
|
||
result->contents [--i] = tem;
|
||
}
|
||
}
|
||
|
||
if (negative) result->sign = 1;
|
||
flush_trailing_digits (result, 0);
|
||
return result;
|
||
}
|
||
|
||
/* Add b1 and b2, considering their signs */
|
||
|
||
decimal
|
||
decimal_add (b1, b2)
|
||
decimal b1, b2;
|
||
{
|
||
decimal v;
|
||
|
||
if (b1->sign != b2->sign)
|
||
v = decimal_sub1 (b1, b2);
|
||
else
|
||
v = decimal_add1 (b1, b2);
|
||
if (b1->sign && !decimal_zerop (v))
|
||
v->sign = !v->sign;
|
||
return v;
|
||
}
|
||
|
||
/* Add b1 and minus b2, considering their signs */
|
||
|
||
decimal
|
||
decimal_sub (b1, b2)
|
||
decimal b1, b2;
|
||
{
|
||
decimal v;
|
||
|
||
if (b1->sign != b2->sign)
|
||
v = decimal_add1 (b1, b2);
|
||
else
|
||
v = decimal_sub1 (b1, b2);
|
||
if (b1->sign && !decimal_zerop (v))
|
||
v->sign = !v->sign;
|
||
return v;
|
||
}
|
||
|
||
/* Return the negation of b2. */
|
||
|
||
decimal
|
||
decimal_neg (b2)
|
||
decimal b2;
|
||
{
|
||
decimal v = decimal_copy (b2);
|
||
|
||
if (!decimal_zerop (v))
|
||
v->sign = !v->sign;
|
||
return v;
|
||
}
|
||
|
||
/* add magnitudes of b1 and b2, ignoring their signs. */
|
||
|
||
static decimal
|
||
decimal_add1 (b1, b2)
|
||
decimal b1, b2;
|
||
{
|
||
int before = MAX (b1->before, b2->before);
|
||
int after = MAX (b1->after, b2->after);
|
||
|
||
int len = before+after+1;
|
||
decimal result = make_decimal (before+1, after);
|
||
|
||
int i;
|
||
char *s1 = b1->contents;
|
||
char *s2 = b2->contents;
|
||
char *p1 = s1 + b1->after - after;
|
||
char *p2 = s2 + b2->after - after;
|
||
char *e1 = s1 + b1->before + b1->after;
|
||
char *e2 = s2 + b2->before + b2->after;
|
||
char *pr = result->contents;
|
||
int accum = 0;
|
||
|
||
for (i = 0; i < len; i++, p1++, p2++)
|
||
{
|
||
accum /= RADIX;
|
||
if (p1 >= s1 && p1 < e1) accum += *p1;
|
||
if (p2 >= s2 && p2 < e2) accum += *p2;
|
||
*pr++ = accum % RADIX;
|
||
}
|
||
if (!accum)
|
||
(result->before)--;
|
||
|
||
flush_trailing_digits (result, 0);
|
||
|
||
return result;
|
||
}
|
||
|
||
/* subtract magnitude of b2 from that or b1, returning signed decimal
|
||
number. */
|
||
|
||
static decimal
|
||
decimal_sub1 (b1, b2)
|
||
decimal b1, b2;
|
||
{
|
||
int before = MAX (b1->before, b2->before);
|
||
int after = MAX (b1->after, b2->after);
|
||
|
||
int len = before+after;
|
||
decimal result = make_decimal (before, after);
|
||
|
||
int i;
|
||
char *s1 = b1->contents;
|
||
char *s2 = b2->contents;
|
||
char *p1 = s1 + b1->after - after;
|
||
char *p2 = s2 + b2->after - after;
|
||
char *e1 = s1 + b1->before + b1->after;
|
||
char *e2 = s2 + b2->before + b2->after;
|
||
char *pr = result->contents;
|
||
int accum = 0;
|
||
|
||
for (i = 0; i < len; i++, p1++, p2++)
|
||
{
|
||
if (p1 >= s1 && p1 < e1) accum += *p1;
|
||
if (p2 >= s2 && p2 < e2) accum -= *p2;
|
||
if (accum < 0 && accum % RADIX)
|
||
*pr = RADIX - (- accum) % RADIX;
|
||
else
|
||
*pr = accum % RADIX;
|
||
accum -= *pr++;
|
||
accum /= RADIX;
|
||
}
|
||
|
||
/* If result is negative, subtract it from RADIX**length
|
||
so that we get the right digits for sign-magnitude
|
||
rather than RADIX-complement */
|
||
|
||
if (accum)
|
||
{
|
||
result->sign = 1;
|
||
pr = result->contents;
|
||
accum = 0;
|
||
for (i = 0; i < len; i++)
|
||
{
|
||
accum -= *pr;
|
||
if (accum)
|
||
*pr = accum + RADIX;
|
||
else
|
||
*pr = 0;
|
||
accum -= *pr++;
|
||
accum /= RADIX;
|
||
}
|
||
}
|
||
|
||
/* flush leading nonfraction zero digits */
|
||
|
||
while (result->before && *--pr == 0)
|
||
(result->before)--;
|
||
|
||
flush_trailing_digits (result, 0);
|
||
|
||
return result;
|
||
}
|
||
|
||
/* multiply b1 and b2 keeping `digits' fraction digits */
|
||
|
||
decimal
|
||
decimal_mul_rounded (b1, b2, digits)
|
||
decimal b1, b2;
|
||
int digits;
|
||
{
|
||
decimal tem = decimal_mul (b1, b2);
|
||
decimal result = decimal_round_digits (tem, digits);
|
||
free (tem);
|
||
return result;
|
||
}
|
||
|
||
/* multiply b1 and b2 keeping the right number of fraction digits
|
||
for the `dc' program with precision = `digits'. */
|
||
|
||
decimal
|
||
decimal_mul_dc (b1, b2, digits)
|
||
decimal b1, b2;
|
||
int digits;
|
||
{
|
||
decimal tem = decimal_mul (b1, b2);
|
||
decimal result
|
||
= decimal_round_digits (tem, MAX (digits, MAX (b1->after, b2->after)));
|
||
free (tem);
|
||
return result;
|
||
}
|
||
|
||
/* multiply b1 and b2 as decimal error-free values;
|
||
keep LENGTH(b1) plus LENGTH(b2) significant figures. */
|
||
|
||
decimal
|
||
decimal_mul (b1, b2)
|
||
decimal b1, b2;
|
||
{
|
||
decimal result = make_decimal (b1->before + b2->before, b1->after + b2->after);
|
||
int i;
|
||
int length2 = LENGTH(b2);
|
||
char *pr;
|
||
|
||
for (i = 0; i < length2; i++)
|
||
add_scaled (result, b1, b2->contents[i], i);
|
||
|
||
/* flush leading nonfraction zero digits */
|
||
|
||
pr = result->contents + LENGTH(result);
|
||
while (result->before && *--pr == 0)
|
||
(result->before)--;
|
||
|
||
flush_trailing_digits (result, 0); /* flush trailing zeros */
|
||
|
||
/* Set sign properly */
|
||
|
||
if (b1->sign != b2->sign && LENGTH(result))
|
||
result->sign = 1;
|
||
|
||
return result;
|
||
}
|
||
|
||
/* Modify decimal number `into' by adding `from',
|
||
multiplied by `factor' (which should be nonnegative and less than RADIX)
|
||
and shifted left `scale' digits at the least significant end. */
|
||
|
||
static void
|
||
add_scaled (into, from, factor, scale)
|
||
decimal into, from;
|
||
int factor, scale;
|
||
{
|
||
char *pf = from->contents;
|
||
char *pi = into->contents + scale;
|
||
int lengthf = LENGTH(from);
|
||
int lengthi = LENGTH(into) - scale;
|
||
|
||
int accum = 0;
|
||
int i;
|
||
|
||
for (i = 0; i < lengthi; i++)
|
||
{
|
||
accum /= RADIX;
|
||
if (i < lengthf)
|
||
accum += *pf++ * factor;
|
||
accum += *pi;
|
||
*pi++ = accum % RADIX;
|
||
}
|
||
}
|
||
|
||
/* Divide decimal number `b1' by `b2', keeping at most `digits'
|
||
fraction digits.
|
||
Returns the result as a decimal number.
|
||
|
||
When division is not exact, the quotient is truncated toward zero. */
|
||
|
||
decimal
|
||
decimal_div (b1, b2, digits)
|
||
decimal b1, b2;
|
||
int digits;
|
||
{
|
||
decimal result = make_decimal (MAX(1, (int) (1 + b1->before - b2->before)), digits);
|
||
|
||
/* b1copy holds what is left of the dividend,
|
||
that is not accounted for by the quotient digits already known */
|
||
|
||
decimal b1copy = decimal_copy_1 (b1, b2->after + digits);
|
||
int length1 = LENGTH(b1copy);
|
||
int length2 = LENGTH(b2);
|
||
int lengthr = LENGTH(result);
|
||
int i;
|
||
|
||
/* leading_divisor_digits contains the first two divisor digits, as
|
||
an integer */
|
||
|
||
int leading_divisor_digits = b2->contents[length2-1]*RADIX;
|
||
if (length2 > 1)
|
||
leading_divisor_digits += b2->contents[length2-2];
|
||
|
||
if (decimal_zerop (b2))
|
||
{
|
||
decimal_error ("divisor is zero", 0);
|
||
return decimal_copy (DECIMAL_ZERO);
|
||
}
|
||
|
||
/* if (lengthr <= (length1 - length2))
|
||
abort(); */ /* My reasoning says this cannot happen, I hope */
|
||
|
||
for (i = length1 - length2; i >= 0; i--)
|
||
{
|
||
/* Guess the next quotient digit (in order of decreasing significance)
|
||
using integer division */
|
||
|
||
int guess;
|
||
int trial_dividend = b1copy->contents[length2+i-1]*RADIX;
|
||
if (i != length1 - length2)
|
||
trial_dividend += b1copy->contents[length2+i]*RADIX*RADIX;
|
||
if (length2 + i > 1)
|
||
trial_dividend += b1copy->contents[length2+i-2];
|
||
|
||
guess = trial_dividend / leading_divisor_digits;
|
||
|
||
/* Remove the quotient times this digit from the dividend left */
|
||
/* We may find that the quotient digit is too large,
|
||
when we consider the entire divisor.
|
||
Then we decrement the quotient digit and add the divisor back in */
|
||
|
||
if (guess && 0 > subtract_scaled (b1copy, b2, guess, i))
|
||
{
|
||
guess--;
|
||
add_scaled (b1copy, b2, 1, i);
|
||
}
|
||
|
||
if (guess >= RADIX)
|
||
{
|
||
result->contents[i + 1] += guess / RADIX;
|
||
guess %= RADIX;
|
||
}
|
||
result->contents[i] = guess;
|
||
}
|
||
|
||
free (b1copy);
|
||
|
||
result->sign = (b1->sign != b2->sign);
|
||
|
||
/* flush leading nonfraction zero digits */
|
||
|
||
{
|
||
char *pr = result->contents + lengthr;
|
||
while (result->before && *--pr == 0)
|
||
(result->before)--;
|
||
}
|
||
|
||
flush_trailing_digits (result, 0); /* Flush trailing zero fraction digits */
|
||
|
||
return result;
|
||
}
|
||
|
||
/* The remainder for the above division.
|
||
Same as `b1' - (`b1' / `b2') * 'b2'.
|
||
Note that the value depends on the number of fraction digits
|
||
that were kept in computing `b1' / `b2';
|
||
the argument `digits' specifies this.
|
||
|
||
The remainder has the same sign as the dividend.
|
||
The divisor's sign is ignored. */
|
||
|
||
decimal
|
||
decimal_rem (b1, b2, digits)
|
||
decimal b1, b2;
|
||
int digits;
|
||
{
|
||
decimal b1copy = decimal_copy_1 (b1, b2->after + digits);
|
||
int length1 = LENGTH(b1copy);
|
||
int length2 = LENGTH(b2);
|
||
int i;
|
||
|
||
int leading_divisor_digits = b2->contents[length2-1]*RADIX;
|
||
|
||
if (length2 > 1)
|
||
leading_divisor_digits += b2->contents[length2-2];
|
||
|
||
if (decimal_zerop (b2))
|
||
{
|
||
decimal_error ("divisor is zero", 0);
|
||
return decimal_copy (DECIMAL_ZERO);
|
||
}
|
||
|
||
/* Do like division, above, but throw away the quotient.
|
||
Keep only the final `rest of dividend', which becomes the remainder. */
|
||
|
||
for (i = length1 - length2; i >= 0; i--)
|
||
{
|
||
int guess;
|
||
int trial_dividend = b1copy->contents[length2+i-1]*RADIX;
|
||
if (i != length1 - length2)
|
||
trial_dividend += b1copy->contents[length2+i]*RADIX*RADIX;
|
||
if (length2 + i > 1)
|
||
trial_dividend += b1copy->contents[length2+i-2];
|
||
|
||
guess = trial_dividend / leading_divisor_digits;
|
||
|
||
if (guess && 0 > subtract_scaled (b1copy, b2, guess, i))
|
||
{
|
||
guess--;
|
||
add_scaled (b1copy, b2, 1, i);
|
||
}
|
||
/* No need to check whether guess exceeds RADIX
|
||
since we are not saving guess. */
|
||
}
|
||
|
||
/* flush leading nonfraction zero digits */
|
||
|
||
{
|
||
char *pr = b1copy->contents + length1;
|
||
while (b1copy->before && *--pr == 0)
|
||
(b1copy->before)--;
|
||
}
|
||
|
||
flush_trailing_digits (b1copy, 0);
|
||
return b1copy;
|
||
}
|
||
|
||
/* returns negative number if we chose factor too large */
|
||
|
||
static int
|
||
subtract_scaled (into, from, factor, scale)
|
||
decimal into, from;
|
||
int factor, scale;
|
||
{
|
||
char *pf = from->contents;
|
||
char *pi = into->contents + scale;
|
||
int lengthf = LENGTH(from);
|
||
int lengthi = LENGTH(into) - scale;
|
||
int accum = 0;
|
||
int i;
|
||
|
||
for (i = 0; i < lengthi && i <= lengthf; i++)
|
||
{
|
||
if (i < lengthf)
|
||
accum -= *pf++ * factor;
|
||
accum += *pi;
|
||
if (accum < 0 && accum % RADIX)
|
||
*pi = RADIX - (- accum) % RADIX;
|
||
else
|
||
*pi = accum % RADIX;
|
||
accum -= *pi++;
|
||
accum /= RADIX;
|
||
}
|
||
return accum;
|
||
}
|
||
|
||
/* Return the square root of decimal number D, using Newton's method.
|
||
Number of fraction digits returned is max of FRAC_DIGITS
|
||
and D's number of fraction digits. */
|
||
|
||
decimal
|
||
decimal_sqrt (d, frac_digits)
|
||
decimal d;
|
||
int frac_digits;
|
||
{
|
||
decimal guess;
|
||
int notdone = 1;
|
||
|
||
if (decimal_zerop (d)) return d;
|
||
if (d->sign)
|
||
{
|
||
decimal_error ("square root argument negative", 0);
|
||
return decimal_copy (DECIMAL_ZERO);
|
||
}
|
||
|
||
frac_digits = MAX (frac_digits, d->after);
|
||
|
||
/* Compute an initial guess by taking the square root
|
||
of a nearby power of RADIX. */
|
||
|
||
if (d->before)
|
||
{
|
||
guess = make_decimal ((d->before + 1) / 2, 0);
|
||
guess->contents[guess->before - 1] = 1;
|
||
}
|
||
else
|
||
{
|
||
/* Arg is less than 1; compute nearest power of RADIX */
|
||
char *p = d->contents + LENGTH(d);
|
||
char *sp = p;
|
||
|
||
while (!*--p); /* Find most significant nonzero digit */
|
||
if (sp - p == 1)
|
||
{
|
||
/* Arg is bigger than 1/RADIX; use 1 as a guess */
|
||
guess = decimal_copy (DECIMAL_ONE);
|
||
}
|
||
else
|
||
{
|
||
guess = make_decimal (0, (sp - p) / 2);
|
||
guess->contents[0] = 1;
|
||
}
|
||
}
|
||
|
||
/* Iterate doing guess = (guess + d/guess) / 2 */
|
||
|
||
while (notdone)
|
||
{
|
||
decimal tem1 = decimal_div (d, guess, frac_digits + 1);
|
||
decimal tem2 = decimal_add (guess, tem1);
|
||
decimal tem3 = decimal_mul_rounded (tem2, DECIMAL_HALF, frac_digits);
|
||
notdone = decimal_compare (guess, tem3);
|
||
free (tem1);
|
||
free (tem2);
|
||
free (guess);
|
||
guess = tem3;
|
||
if (decimal_zerop (guess)) return guess; /* Avoid divide-by-zero */
|
||
}
|
||
|
||
return guess;
|
||
}
|
||
|
||
/* Raise decimal number `base' to power of integer part of decimal
|
||
number `expt'.
|
||
This function depends on using radix 10.
|
||
It is too hard to write it to work for any value of RADIX,
|
||
so instead it is simply not available if RADIX is not ten. */
|
||
|
||
#if !(RADIX - 10)
|
||
|
||
decimal
|
||
decimal_expt (base, expt, frac_digits)
|
||
decimal base, expt;
|
||
int frac_digits;
|
||
{
|
||
decimal accum = decimal_copy (DECIMAL_ONE);
|
||
decimal basis1 = base;
|
||
int digits = expt->before;
|
||
int dig = 0; /* Expt digit being processed */
|
||
|
||
if (expt->sign)
|
||
/* If negative power, take reciprocal first thing
|
||
so that fraction digit truncation won't destroy
|
||
what will ultimately be nonfraction digits. */
|
||
basis1 = decimal_div (DECIMAL_ONE, base, frac_digits);
|
||
while (dig < digits)
|
||
{
|
||
decimal basis2, basis4, basis8, basis10;
|
||
int thisdigit = expt->contents[expt->after + dig];
|
||
|
||
/* Compute factors to multiply in for each bit of this digit */
|
||
|
||
basis2 = decimal_mul_rounded (basis1, basis1, frac_digits);
|
||
basis4 = decimal_mul_rounded (basis2, basis2, frac_digits);
|
||
basis8 = decimal_mul_rounded (basis4, basis4, frac_digits);
|
||
|
||
/* Now accumulate the factors this digit value selects */
|
||
|
||
if (thisdigit & 1)
|
||
{
|
||
decimal accum1 = decimal_mul_rounded (accum, basis1, frac_digits);
|
||
free (accum);
|
||
accum = accum1;
|
||
}
|
||
|
||
if (thisdigit & 2)
|
||
{
|
||
decimal accum1 = decimal_mul_rounded (accum, basis2, frac_digits);
|
||
free (accum);
|
||
accum = accum1;
|
||
}
|
||
|
||
if (thisdigit & 4)
|
||
{
|
||
decimal accum1 = decimal_mul_rounded (accum, basis4, frac_digits);
|
||
free (accum);
|
||
accum = accum1;
|
||
}
|
||
|
||
if (thisdigit & 8)
|
||
{
|
||
decimal accum1 = decimal_mul_rounded (accum, basis8, frac_digits);
|
||
free (accum);
|
||
accum = accum1;
|
||
}
|
||
|
||
/* If there are further digits, compute the basis1 for the next digit */
|
||
|
||
if (++dig < digits)
|
||
basis10 = decimal_mul_rounded (basis2, basis8, frac_digits);
|
||
|
||
/* Free intermediate results */
|
||
|
||
if (basis1 != base) free (basis1);
|
||
free (basis2);
|
||
free (basis4);
|
||
free (basis8);
|
||
basis1 = basis10;
|
||
}
|
||
return accum;
|
||
}
|
||
#endif
|
||
|
||
#ifdef TEST
|
||
|
||
fputchar (c)
|
||
char c;
|
||
{
|
||
putchar (c);
|
||
}
|
||
|
||
/* Top level that can be used to test the arithmetic functions */
|
||
|
||
main ()
|
||
{
|
||
char s1[40], s2[40];
|
||
decimal b1, b2, b3;
|
||
char c;
|
||
|
||
while (1)
|
||
{
|
||
scanf ("%s %c %s", s1, &c, s2);
|
||
b1 = decimal_parse (s1, RADIX);
|
||
b2 = decimal_parse (s2, RADIX);
|
||
switch (c)
|
||
{
|
||
default:
|
||
c = '+';
|
||
case '+':
|
||
b3 = decimal_add (b1, b2);
|
||
break;
|
||
case '*':
|
||
b3 = decimal_mul (b1, b2);
|
||
break;
|
||
case '/':
|
||
b3 = decimal_div (b1, b2, 3);
|
||
break;
|
||
case '%':
|
||
b3 = decimal_rem (b1, b2, 3);
|
||
break;
|
||
case 'p':
|
||
decimal_print (b1, fputchar, RADIX);
|
||
printf (" printed in base %d is ", decimal_to_int (b2));
|
||
decimal_print (b1, fputchar, decimal_to_int (b2));
|
||
printf ("\n");
|
||
continue;
|
||
case 'r':
|
||
printf ("%s read in base %d is ", s1, decimal_to_int (b2));
|
||
decimal_print (decimal_parse (s1, decimal_to_int (b2)), fputchar, RADIX);
|
||
printf ("\n");
|
||
continue;
|
||
}
|
||
decimal_print (b1, fputchar, RADIX);
|
||
printf (" %c ", c);
|
||
decimal_print (b2, fputchar, RADIX);
|
||
printf (" = ");
|
||
decimal_print (b3, fputchar, RADIX);
|
||
printf ("\n");
|
||
}
|
||
}
|
||
|
||
decimal_error (s1, s2)
|
||
char *s1, *s2;
|
||
{
|
||
printf ("\n");
|
||
printf (s1, s2);
|
||
printf ("\n");
|
||
}
|
||
|
||
static void
|
||
pbi (b)
|
||
int b;
|
||
{
|
||
decimal_print ((decimal) b, fputchar, RADIX);
|
||
}
|
||
|
||
static void
|
||
pb (b)
|
||
decimal b;
|
||
{
|
||
decimal_print (b, fputchar, RADIX);
|
||
}
|
||
|
||
#endif
|