diff options
| -rw-r--r-- | big.c | 250 |
1 files changed, 124 insertions, 126 deletions
@@ -40,8 +40,6 @@ struct big_div { word big_one; -static struct big_div big__longdiv(word, word); -static struct big_div shortdiv(word, word); static word shift(word, word); static word stimes(word, word); @@ -79,6 +77,128 @@ word big__ms2d(word x) return (digit(x) * IBASE + d); } +/* divide big x by single digit n returning big quotient + and setting external b_rem as side effect */ +/* may assume - x>=0,n>0 */ +struct big_div big__shortdiv(word x, word n) +{ + word d = digit(x); + word q = 0; + + /* reverse rest(x) into q */ + while ((x = rest(x)) != 0) { + q = make(INT, d, q); + d = digit(x); /* leaving most sig. digit in d */ + } + + x = q; + word s_rem = d % n; + d = d / n; + + /* put back first digit (if not leading 0) */ + if (d || !q) { + q = make(INT, d, 0); + } else { + q = 0; + } + + /* in situ division of q by n AND destructive reversal */ + while (x) { + d = s_rem * IBASE + digit(x); + digit(x) = d / n; + s_rem = d % n; + word tmp = x; + x = rest(x); + rest(tmp) = q; + q = tmp; + } + + struct big_div bd = { + .quo = q, + .rem = make(INT, s_rem, 0), + }; + return bd; +} + +/* divide big x by big y returning quotient, leaving + remainder in extern variable b_rem */ +/* may assume - x>=0,y>0 */ +struct big_div big__longdiv(word x, word y) +{ + if (bigcmp(x, y) < 0) { + struct big_div bd = { + .quo = make(INT, 0, 0), + .rem = x, + }; + return bd; + } + + word y1 = big__msd(y); + + /* rescale if necessary */ + word scale = IBASE / (y1 + 1); + if (scale > 1) { + x = stimes(x, scale); + y = stimes(y, scale); + y1 = big__msd(y); + } + + word n = 0; + word q = 0; + word ly = big__len(y); + while (bigcmp(x, y = make(INT, 0, y)) >= 0) { + n++; + } + + y = rest(y); /* want largest y not exceeding x */ + ly += n; + for (;;) { + word d, lx = big__len(x); + if (lx < ly) { + d = 0; + + } else if (lx == ly) { + if (bigcmp(x, y) >= 0) { + x = bigsub(x, y); + d = 1; + } else { + d = 0; + } + + } else { + d = big__ms2d(x) / y1; + if (d > MAXDIGIT) { + d = MAXDIGIT; + } + + if ((d -= 2) > 0) { + x = bigsub(x, stimes(y, d)); + } else { + d = 0; + } + + if (bigcmp(x, y) >= 0) { + x = bigsub(x, y), d++; + if (bigcmp(x, y) >= 0) { + x = bigsub(x, y); + d++; + } + } + } + + q = make(INT, d, q); + if (n-- == 0) { + struct big_div bd = { + .quo = q, + .rem = scale == 1 ? x : big__shortdiv(x, scale).rem, + }; + return bd; + } + ly--; + y = rest(y); + } +} /* see Bird & Wadler p82 for explanation */ + /* ignore input signs, treat x,y as positive */ static word big__plus(word x, word y, int signbit) @@ -421,7 +541,7 @@ word big_divide(word x, word y) s2 = !s1; } - struct big_div bd = rest(y) ? big__longdiv(x, y) : shortdiv(x, digit(y)); + struct big_div bd = rest(y) ? big__longdiv(x, y) : big__shortdiv(x, digit(y)); word q = bd.quo; if (s2) { @@ -461,7 +581,7 @@ word big_remainder(word x, word y) s2 = !s1; } - struct big_div bd = rest(y) ? big__longdiv(x, y) : shortdiv(x, digit(y)); + struct big_div bd = rest(y) ? big__longdiv(x, y) : big__shortdiv(x, digit(y)); if (s2) { if (!bigzero(bd.rem)) { @@ -481,128 +601,6 @@ word big_remainder(word x, word y) magnitudes invariant under change of sign, remainder has sign of dividend, quotient negative if signs of divi(sor/dend) mixed */ -/* divide big x by single digit n returning big quotient - and setting external b_rem as side effect */ -/* may assume - x>=0,n>0 */ -struct big_div shortdiv(word x, word n) -{ - word d = digit(x); - word q = 0; - - /* reverse rest(x) into q */ - while ((x = rest(x)) != 0) { - q = make(INT, d, q); - d = digit(x); /* leaving most sig. digit in d */ - } - - x = q; - word s_rem = d % n; - d = d / n; - - /* put back first digit (if not leading 0) */ - if (d || !q) { - q = make(INT, d, 0); - } else { - q = 0; - } - - /* in situ division of q by n AND destructive reversal */ - while (x) { - d = s_rem * IBASE + digit(x); - digit(x) = d / n; - s_rem = d % n; - word tmp = x; - x = rest(x); - rest(tmp) = q; - q = tmp; - } - - struct big_div bd = { - .quo = q, - .rem = make(INT, s_rem, 0), - }; - return bd; -} - -/* divide big x by big y returning quotient, leaving - remainder in extern variable b_rem */ -/* may assume - x>=0,y>0 */ -struct big_div big__longdiv(word x, word y) -{ - if (bigcmp(x, y) < 0) { - struct big_div bd = { - .quo = make(INT, 0, 0), - .rem = x, - }; - return bd; - } - - word y1 = big__msd(y); - - /* rescale if necessary */ - word scale = IBASE / (y1 + 1); - if (scale > 1) { - x = stimes(x, scale); - y = stimes(y, scale); - y1 = big__msd(y); - } - - word n = 0; - word q = 0; - word ly = big__len(y); - while (bigcmp(x, y = make(INT, 0, y)) >= 0) { - n++; - } - - y = rest(y); /* want largest y not exceeding x */ - ly += n; - for (;;) { - word d, lx = big__len(x); - if (lx < ly) { - d = 0; - - } else if (lx == ly) { - if (bigcmp(x, y) >= 0) { - x = bigsub(x, y); - d = 1; - } else { - d = 0; - } - - } else { - d = big__ms2d(x) / y1; - if (d > MAXDIGIT) { - d = MAXDIGIT; - } - - if ((d -= 2) > 0) { - x = bigsub(x, stimes(y, d)); - } else { - d = 0; - } - - if (bigcmp(x, y) >= 0) { - x = bigsub(x, y), d++; - if (bigcmp(x, y) >= 0) { - x = bigsub(x, y); - d++; - } - } - } - - q = make(INT, d, q); - if (n-- == 0) { - struct big_div bd = { - .quo = q, - .rem = scale == 1 ? x : shortdiv(x, scale).rem, - }; - return bd; - } - ly--; - y = rest(y); - } -} /* see Bird & Wadler p82 for explanation */ - /* assumes y big_is_positive */ word bigpow(word x, word y) { |