diff options
Diffstat (limited to 'gnu/usr.bin/dc/decimal.c')
| -rw-r--r-- | gnu/usr.bin/dc/decimal.c | 1235 |
1 files changed, 1235 insertions, 0 deletions
diff --git a/gnu/usr.bin/dc/decimal.c b/gnu/usr.bin/dc/decimal.c new file mode 100644 index 000000000000..1fb95c18aa14 --- /dev/null +++ b/gnu/usr.bin/dc/decimal.c @@ -0,0 +1,1235 @@ +/* + * 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_round_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 |