extend namespace Math { public real pi = imprecise (3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679, 256); public real π = pi; protected real e = imprecise (2.7182818284590452353602874713526624977572470936999595749669676277240766303535475945713821785251664274, 256); public real sqrt (real v) /* * Returns square root of v to the same precision as 'v'. * If v is precise and has a precise square root, returns that. */ { if (v < 0) raise invalid_argument ("sqrt of negative number", 0, v); real real_sqrt(real v) { real err; real prev, cur; int n, iter; v = imprecise (v); err = 1/(2**(precision (v)+3)); prev = imprecise (1 / (2**(floor (exponent(v)/2))), precision(v)); iter = precision (v) + 15; while (iter-- > 0) { cur = 0.5 * prev * (3 - v * prev**2); if (abs (cur - prev) < err) break; prev = cur; } return abs (1/cur); } if (is_rational (v)) { int num, den; real num_s, den_s; num = numerator (v); den = denominator (v); num_s = real_sqrt (imprecise (num, bit_width(num) + 128)); den_s = real_sqrt (imprecise (den, bit_width(den) + 128)); num = floor (num_s + 0.5); den = floor (den_s + 0.5); if (num * num == numerator (v) && den * den == denominator (v)) { return num/den; } } return real_sqrt (v); } public real cbrt (real v) /* * Returns cube root of v to the same precision as 'v'. * If v is precise and has a precise cube root, returns that. */ { real real_cbrt(real v) { real prev, cur; int s = sign (v); v = imprecise (abs (v)); int result_bits = precision (v); int intermediate_bits = result_bits + 10; v = imprecise (v, intermediate_bits); real err = imprecise (1/(2**(result_bits+3)), intermediate_bits); cur = imprecise (1 / (0.75 * 2**(floor (exponent(v)/3))), intermediate_bits); do { prev = cur; cur = 1/3 * (2 * prev + v / prev**2); } while (abs (cur - prev) > err); return s * imprecise (abs (cur), result_bits); } if (is_rational (v)) { int num, den; real num_s, den_s; num = numerator (v); den = denominator (v); num_s = real_cbrt (imprecise (num, bit_width(num) + 128)); den_s = real_cbrt (imprecise (den, bit_width(den) + 128)); num = floor (num_s + 0.5); den = floor (den_s + 0.5); /* printf ("num %g den %g\n", num, den); */ if (num ** 3 == numerator (v) && den ** 3 == denominator (v)) { return num/den; } } return real_cbrt (v); } public real exp (real v) /* * Return e ** v; */ { if (v < 0) return 1/exp(-v); if (v == 0) return 1; v = imprecise (v); /* * Emperically determined scale factor. This * reduces the computation to working on values * near zero so that the power series converges * rapidly. Increasing this further makes the * power series converge more rapidly, but * makes the expansion step more expensive. */ int prec = precision (v); int scale; if (prec > 50) scale = 27; else scale = 12; int div = (1 << scale); int iter = prec + 1; int iprec = prec + iter; real mant = imprecise (mantissa(v), prec + iter) / div; int expo = exponent(v) + scale; real e = imprecise (0, iprec); real num = imprecise (1, iprec); real den = imprecise (1, iprec); real loop = imprecise (0, iprec); /* * Traditional power series * * exp(n) = 1 + n/1 + n**2/2! + n**3/3! */ while (iter-- > 0) { real term = num/den; e = e + term; if (exponent (e) > exponent(term) + iprec) break; num *= mant; loop++; den *= loop; } e = imprecise (e, prec + expo); e = e ** (1 << expo); return imprecise (e, prec); } public real log (real a) /* * Return natural logarithm of 'a' */ { /* * Copyright (c) 1985 Regents of the University of California. * All rights reserved. * * Redistribution and use in source and binary forms, with or without * modification, are permitted provided that the following conditions * are met: * 1. Redistributions of source code must retain the above copyright * notice, this list of conditions and the following disclaimer. * 2. Redistributions in binary form must reproduce the above copyright * notice, this list of conditions and the following disclaimer in the * documentation and/or other materials provided with the distribution. * 3. All advertising materials mentioning features or use of this software * must display the following acknowledgement: * This product includes software developed by the University of * California, Berkeley and its contributors. * 4. Neither the name of the University nor the names of its contributors * may be used to endorse or promote products derived from this software * without specific prior written permission. * * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF * SUCH DAMAGE. */ /* log__L(Z) * LOG(1+X) - 2S X * RETURN --------------- WHERE Z = S*S, S = ------- , 0 <= Z <= .0294... * S 2 + X * * DOUBLE PRECISION (VAX D FORMAT 56 bits or IEEE DOUBLE 53 BITS) * KERNEL FUNCTION FOR LOG; TO BE USED IN LOG1P, LOG, AND POW FUNCTIONS * CODED IN C BY K.C. NG, 1/19/85; * REVISED BY K.C. Ng, 2/3/85, 4/16/85. * * Method : * 1. Polynomial approximation: let s = x/(2+x). * Based on log(1+x) = log(1+s) - log(1-s) * = 2s + 2/3 s**3 + 2/5 s**5 + ....., * * (log(1+x) - 2s)/s is computed by * * z*(L1 + z*(L2 + z*(... (L7 + z*L8)...))) * * where z=s*s. (See the listing below for Lk's values.) The * coefficients are obtained by a special Remez algorithm. * * Accuracy: * Assuming no rounding error, the maximum magnitude of the approximation * error (absolute) is 2**(-58.49) for IEEE double, and 2**(-63.63) * for VAX D format. * * Constants: * The hexadecimal values are the intended ones for the following constants. * The decimal values may be used, provided that the compiler will convert * from decimal to binary accurately enough to produce the hexadecimal values * shown. */ real log__L (real z) { global real L1 = imprecise (6.6666666666667340202E-1, 64); global real L2 = imprecise (3.9999999999416702146E-1, 64); global real L3 = imprecise (2.8571428742008753154E-1, 64); global real L4 = imprecise (2.2222198607186277597E-1, 64); global real L5 = imprecise (1.8183562745289935658E-1, 64); global real L6 = imprecise (1.5314087275331442206E-1, 64); global real L7 = imprecise (1.4795612545334174692E-1, 64); return(z*(L1+z*(L2+z*(L3+z*(L4+z*(L5+z*(L6+z*L7))))))); } /* LOG(X) * RETURN THE LOGARITHM OF x * DOUBLE PRECISION (VAX D FORMAT 56 bits or IEEE DOUBLE 53 BITS) * CODED IN C BY K.C. NG, 1/19/85; * REVISED BY K.C. NG on 2/7/85, 3/7/85, 3/24/85, 4/16/85. * * Required system supported functions: * scalb(x,n) * copysign(x,y) * logb(x) * finite(x) * * Required kernel function: * log__L(z) * * Method : * 1. Argument Reduction: find k and f such that * x = 2^k * (1+f), * where sqrt(2)/2 < 1+f < sqrt(2) . * * 2. Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s) * = 2s + 2/3 s**3 + 2/5 s**5 + ....., * log(1+f) is computed by * * log(1+f) = 2s + s*log__L(s*s) * where * log__L(z) = z*(L1 + z*(L2 + z*(... (L6 + z*L7)...))) * * See log__L() for the values of the coefficients. * * 3. Finally, log(x) = k*ln2 + log(1+f). (Here n*ln2 will be stored * in two floating point number: n*ln2hi + n*ln2lo, n*ln2hi is exact * since the last 20 bits of ln2hi is 0.) * * Special cases: * log(x) is NaN with signal if x < 0 (including -INF) ; * log(+INF) is +INF; log(0) is -INF with signal; * log(NaN) is that NaN with no signal. * * Accuracy: * log(x) returns the exact log(x) nearly rounded. In a test run with * 1,536,000 random arguments on a VAX, the maximum observed error was * .826 ulps (units in the last place). * * Constants: * The hexadecimal values are the intended ones for the following constants. * The decimal values may be used, provided that the compiler will convert * from decimal to binary accurately enough to produce the hexadecimal values * shown. */ real bsd_log (real x) { global real ln2hi = imprecise (6.9314718036912381649E-1, 64); global real ln2lo = imprecise (1.9082149292705877000E-10, 64); global real sqrt2 = imprecise (1.4142135623730951455E0, 64); global real negone = imprecise (-1.0, 64); global real half = imprecise (0.5, 64); global real two = imprecise (2, 64); real s,z,t; int k,n; /* argument reduction */ k=exponent(x); x=mantissa(x); if(x >= sqrt2 ) {k += 1; x *= half;} x += negone ; /* compute log(1+x) */ s = x/(two + x); t = x*x*half; z = k*ln2lo + s*(t+log__L(s*s)); x += (z - t); return (k*ln2hi+x); } /* * Bounds checking */ if (a <= 0) raise invalid_argument ("log: must be positive", 0, a); /* * Checks to bring a into range */ if (a == 1) return 0; if (a < 1) return -log(1/a); a = imprecise (a); int prec = precision(a); /* * estimate = bsd_log (a). This gives 53 bits */ real v = bsd_log (imprecise (a, 64)); /* * Precision doubles every time around, start * with 50 bits and compute how many doublings are * needed to get the desired precision */ int maxiter = 0; int rprec = 50; while (rprec < prec) { rprec *= 2; maxiter++; } /* printf ("maxiter %g\n", maxiter); */ if (maxiter > 0) { int iprec = prec + maxiter * 16; v = imprecise (v, iprec); a = imprecise (a, iprec); int epow = floor (v); /* * compute log(v) = log(v/(e**epow)) + epow; */ v = v - epow; a /= exp(imprecise (epow, iprec)); /* * Newton's method * * v' = v - 1 + a * exp(-v); */ real one = imprecise (1, iprec); while (maxiter-- > 0) { real term = a * exp (-v) - one; /* printf ("v: %g term: %g err: %g term*err: %g\n", v, term, err, term * err); */ v = v + term; } v = v + epow; } return imprecise (v, prec); } /* * log10(x) = log10(e) * log(x) * * log10(e) = log(e) / log(10) = 1/log(10) */ public real log10 (real a) /* * Return base-10 log of 'a' */ { static real loge = 0; a = imprecise (a); if (loge == 0 || precision (loge) < precision (a)) loge = 1/log(imprecise (10, precision (a))); return loge * log(a); } /* * log2(x) = log2(e) * log(x) * * log2(e) = log(e) / log(2) = 1/log(2) */ public real log2 (real a) /* * Return base-2 log of 'a' */ { static real loge = 0; a = imprecise (a); if (loge == 0 || precision (loge) < precision (a)) loge = 1/log(imprecise (2, precision (a))); return loge * log(a); } real calculate_pi (int prec) /* * Calculate pi using the formula: * * PI = 24*atan (1/8) + 8*atan (1/57) + 4*atan (1/239); */ { /* * Estimate the number of digits available for * the specified value (v) after a certain number of * loops (p) */ real avail_prec (real v, int p) { real ret; ret = bit_width (p) - p * exponent (imprecise (v)); /* printf ("v %g p %g avail %g\n", v, p, ret); */ return ret; } /* * Compute the number of loops needed to get * the desired precision */ int loops (real v, int prec) { int p, low, high; for (high = 1; ; high *= 2) { if (avail_prec (v, high) > prec) break; } low = 1; while (high - low > 1) { p = (high + low) // 2; if (avail_prec (v, p) > prec) high = p; else low = p; } return high; } /* * Compute atan near zero * * atan(x) = x - x**3/3 + x**5/5 - ... */ real atan (rational den, int digits) { int p, q; int l; int prec, mult; real partial, result; real pv, qv, mden; p = 3; q = 5; mden = imprecise (den, digits * 4) ** 4; l = loops (1 / den, digits) // 2; /* * Need at least digits + log10(loops) for all intermediate * computations */ /* printf ("loops %d\n", l); File.flush (stdout); */ result = 1 / den; pv = 1 / (den ** p); qv = 1 / (den ** q); while (l-- > 0) { partial = pv / p - qv / q; if (partial == 0) break; result = result - partial; /* if (l % 10 == 0) { printf ("."); File.flush (stdout); } */ p += 4; q += 4; pv = pv / mden; qv = qv / mden; } /* printf ("\n"); */ return result; } real value; real part1, part2, part3; part1 = 24 *atan (8, prec + 30); part2 = 8 * atan (57, prec + 30); part3 = 4 * atan (239, prec + 30); value = part1 + part2 + part3; return imprecise (value, prec); } public real pi_value (int prec) /* * Return pi at least as precise as 'prec' */ { static real local_pi = pi; if (precision (local_pi) < prec) local_pi = calculate_pi (prec); return imprecise (local_pi, prec); } real limit_angle_to_pi (real aa) { real my_pi; aa = imprecise (aa); my_pi = pi_value (precision (aa)); if (aa > my_pi) aa = aa - 2 * my_pi; return aa % (2 * my_pi); } public real sin (real a) /* * return sine (a) */ { /* * sin(x) = x - x**3/3! + x**5/5! ... */ real raw_sin (real a) { real err; real v, term; real a4, aj, ai; int i, j; int iter; int prec; prec = precision(a); int iprec = prec * 2; a = imprecise(a,iprec); i = 1; j = 3; a4 = a**4; ai = a**i; aj = a**j; iter = prec + 8; v = 0; while (iter-- > 0) { term = ai/i! - aj/j!; /* printf ("sin iter %d term %d\n", iter, term);*/ v += term; if (exponent (v) > exponent (term) + iprec) break; ai *= a4; aj *= a4; i += 4; j += 4; } return imprecise (v + term, prec); } /* sin(5x) = 16 * sin**5(x) - 20 * sin**3(x) + 5 * sin(x) */ real do_5x (real a) { return 16 * a**5 - 20 * a**3 + 5 * a; } real big_sin (real a) { if (a > 0.01) return do_5x (big_sin (a/5)); return raw_sin (a); } a = limit_angle_to_pi (a); if (a == 0) return 0; return big_sin (a); } public real cos (real a) /* * return cosine (a) */ { /* * cos(x) = 1 - x**2/2! + x**4/4! - x**6/6! ... */ real raw_cos (real a) { real v, term; real ai, aj, a4; int i, j; int iter; int prec = precision(a); int iprec = prec * 2; a = imprecise(a, iprec); i = 0; j = 2; ai = 1; aj = a**2; a4 = a**4; iter = prec + 8; v = 0; while (iter-- > 0) { term = ai/i! - aj/j!; v += term; if (exponent (v) > exponent (term) + iprec) break; ai *= a4; aj *= a4; i += 4; j += 4; } return imprecise (v + term); } /* cos(4x) = 8 * (cos**4(x) - cos**2(x)) + 1 */ real do_4x (real c) { return 8 * (c**4 - c**2) + 1; } real big_cos (real a) { if (a > .01) return do_4x (big_cos (a/4)); return raw_cos (a); } a = limit_angle_to_pi (a); if (a == 0) return 1; return big_cos (limit_angle_to_pi (a)); } real cos_to_sin (real v) { return sqrt (1 - v**2); } public void sin_cos (real a, *real sinp, *real cosp) /* * Compute sine and cosine of 'a' simultaneously */ { real c, s; a = limit_angle_to_pi (a); c = cos (a); s = sign(a) * abs (cos_to_sin(c)); *cosp = c; *sinp = s; } public real tan (real a) /* * return tangent (a) */ { real c, s; a = imprecise (a); sin_cos (a, &s, &c); return s/c; } public real atan (real v) /* * return arctangent (v) */ { /* * atan(x) = x - x**3/3 + x**5/5 - ... */ real raw_atan (real v) { real a, term; real vi, vj, v4; int i, j; int iter; int prec = precision(v); int iprec = prec * 2; v = imprecise (v, iprec); i = 1; j = 3; vi = v**i; vj = v**j; v4 = v**4; a = 0; iter = prec + 8; while (iter-- > 0) { term = vi/i - vj/j; a += term; if (exponent (a) > exponent (term) + iprec) break; vi *= v4; vj *= v4; i += 4; j += 4; } return imprecise (a, prec); } real sqrt3; v = imprecise (v); /* * atan(v) = -atan(-v) */ if (v < 0) return -atan (-v); /* * atan(v) = pi/2 - atan(1/v) */ if (v > 1) return pi_value (precision(v))/2 - atan (1/v); /* * atan(v) = pi/6 + atan((v*sqrt(3) - 1) / (sqrt(3) + v)) */ if (v > .268) { sqrt3 = sqrt (imprecise (3,precision(v))); return (pi_value (precision(v)) / 6 + raw_atan ((v * sqrt3 - 1) / (sqrt3 + v))); } return raw_atan (v); } /* * atan(v) = asin(v/sqrt(1+v**2)) * * q = v/sqrt(1+v**2) * q*sqrt(1+v**2) = v * q**2*(1+v**2) = v**2 * q**2 + q**2v**2 = v**2 * q**2 = v**2 - q**2v**2 * q**2 = v**2 * (1 - q**2) * v**2 = q**2/(1-q**2) * v = q/sqrt(1-q**2) * * asin(q) = atan(q/sqrt(1-q**2)) */ public real asin (real v) /* * return arcsine (v) */ { v = imprecise (v); if (abs (v) > 1) raise invalid_argument ("asin argument out of range", 0, v); if (v == 1) return pi_value (precision (v))/2; if (v == -1) return -pi_value (precision (v))/2; return atan (v/sqrt(1-v**2)); } /* * acos(v) = asin (sqrt (1 - v**2)) * = atan (sqrt(1-v**2) / sqrt (1-(sqrt (1-v**2))**2)) * = atan (sqrt(1-v**2) / sqrt (1-(1-v**2))) * = atan (sqrt(1-v**2) / v) */ public real acos (real v) /* * return arccosine (v) */ { v = imprecise(v); if (abs (v) > 1) raise invalid_argument ("acos argument out of range", 0, v); if (v == 1) return 0; if (v == -1) return pi_value(precision(v)); if (v == 0) return pi_value(precision(v))/2; return atan (sqrt (1-v**2)/v); } /* * atan (y/x) */ public real atan2 (real y, real x) /* * return atan (y/x), but adjust for quadrant correctly */ { y = imprecise (y); x = imprecise (x); if (y == 0) { if (x >= 0) return 0; return pi_value (precision (y)); } if (x == 0) return sign(y) * pi_value (precision (y))/2; real a; a = atan (y/x); if (x < 0) { if (y >= 0) a += pi_value (precision (a)); else a -= pi_value (precision (a)); } return a; } /* * These two are used for the '**' and '**=' operators */ public real pow (real a, real b) /* * return a ** b; */ { real result; if (is_int (b)) { if (!is_int (a) && is_rational (a)) return pow (numerator(a), b) / pow (denominator (a), b); int flip = 0; if (b < 0) { flip = 1; b = -b; } result = 1; while (b > 0) { if (b % 2 != 0) result *= a; if ((b //= 2) != 0) a *= a; } if (flip != 0) result = 1/result; } else switch (b) { case .5: result = sqrt (a); break; case .{3}: result = cbrt (a); break; default: result = exp (b * log(a)); break; } return result; } public real assign_pow (*real a, real b) /* * return *a = *a ** b; */ { return *a = pow (*a, b); } public real max(real arg, real args ...) /* * Return maximum of all arguments */ { for (int i = 0; i < dim(args); i++) if (arg < args[i]) arg = args[i]; return arg; } public real min(real arg, real args ...) /* * Return minimum of all arguments */ { for (int i = 0; i < dim(args); i++) if (arg > args[i]) arg = args[i]; return arg; } /* * Fast integer logarithm via binary search from below (no division). * Returns floor(log(n)/log(base)) with no rounding error */ public int ilog(int base, int n) /* * Fast integer logarithm via binary search from below (no division). * Returns floor(log(n)/log(base)) with no rounding error */ { if (base <= 1) raise invalid_argument("ilog of bad base", 0, base); if (n <= 0) raise invalid_argument("ilog of bad value", 1, n); int below = 0; int above = 1; int k = base; while (k <= n) { k *= k; below = above; above *= 2; } while (true) { int q = base ** below; k = base; int nbelow = 0; int nabove = 1; while (q * k <= n) { k *= k; nbelow = nabove; nabove *= 2; } if (nbelow == 0) break; below += nbelow; } return below; } public exception lsb_0(); public int lsb(int b) /* * return the bit position of * the least significant bit of the int argument * via binary search */ { global bool mask(int b, int ul) { return (b & ((1 << (ul + 1)) - 1)) != 0; } if (b == 0) raise lsb_0(); if (b == -1) return 0; /* doubling phase */ int ul = 1; for (!mask(b, ul); ul *= 2) /* do nothing */; /* binary search phase */ int ll = 0; while (ul > ll + 1) { int step = (ul - ll) // 2; if (mask(b, ul - step)) { ul -= step; continue; } if (!mask(b, ll + step)) { ll += step; continue; } abort("error in binary search"); } if (mask(b, ll)) return ll; return ul; } } /* XXX these shouldn't be here, but it was *convenient* */ &int(string, int ...) atoi = &string_to_integer; &rational(string) atof = &string_to_real;