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;

	    v = imprecise (v);
	    int prec = precision(v);
	    int iprec = prec + 5 * bit_width(prec) + 10;
	    real vi = imprecise(v, iprec);
	    /* Initial estimate is 1/2 ** log2(v)/2 */
	    cur = imprecise(0.5 * (2**(exponent(v) ⫽ 2)), iprec);
	    int iter = iprec + 100;
	    do {
		/*
		 * Newton's method iteration
		 *
		 * f(y)  = y² - x
		 * f'(y) = 2y
		 *
		 * y'    = y - f(y)/f'(y)
		 *       = y - (y² - x) / 2y
		 *       = y - y/2 + x/2y
		 *       = y/2 + x/2y
		 */
		prev = cur;
		cur = prev/2 + vi/(2*prev);
		/* bail if we take too long */
		if (--iter <= 0)
		    break;
	    } while (imprecise(cur, prec+2) != imprecise(prev, prec+2));
	    return abs(imprecise(cur, prec));
	}

	if (v == 0)
		return 0;

	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 prec = precision (v);
	    /*
	     * Estimate about log2(prec) iterations, add 5 bits
	     * of precision for each one (for the 5 operations)
	     */
	    int iprec = prec + 5 * bit_width(prec) + 10;
	    v = imprecise (v, iprec);
	    /* Initial estimate is 3/4 ** log2(v)/3 */
	    cur = imprecise (0.75 * 2**(exponent(v) ⫽ 3), iprec);
	    int iter = iprec + 100;
	    do {
		prev = cur;
		/*
		 * Newton's method iteration
		 *
		 * f(y)  = y³ - x
		 * f'(y) = 3y²
		 *
		 * y'    = y - f(y)/f'(y)
		 *       = y - (y³ - x) / 3y²
		 *       = y - y/3 + x/3y²
		 *       = 2y/3 + x/3y²
		 */
		cur = 0.{6} * prev + 0.{3} * v / (prev*prev);
		/* bail if we take too long */
		if (--iter <= 0)
		    break;
	    } while (imprecise(cur, prec+2) != imprecise(prev, prec+2));
	    return s * imprecise (abs (cur), prec);
	}

        if (v == 0)
	    return 0;

	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);
	    if (num ** 3 == numerator (v) && den ** 3 == denominator (v))
		return num/den;
	}
	return real_cbrt (v);
    }

    /*
     * 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;
    }

    real calculate_e (int prec)
    /*
     * Calculate e recursively
     */
    {
	typedef struct {
	    int	p;
	    int	q;
	} split_t;

	split_t make_split(int p, int q) = (split_t) { .p = p, .q = q };

	/* Compute e-1 recursively */

	split_t binary_splitting_e(int n0, int n1)
	{
	    if (n1 - n0 == 1)
		return make_split(1, n1);

	    int nmid = (n0 + n1) >> 1;

	    split_t r0 = binary_splitting_e(n0, nmid);

	    split_t r1 = binary_splitting_e(nmid, n1);

	    return make_split(r0.p * r1.q + r1.p, r0.q * r1.q);
	}

	int log_factorial = 0;
	int log_max = prec;
	int series_size = 1;

	while (log_factorial < log_max) {
	    log_factorial += ilog(2, series_size);
	    series_size++;
	}

	split_t pq = binary_splitting_e(0, series_size);

	return imprecise(1 + pq.p / pq.q, prec);
    }

    public real e_value (int prec)
	/*
	 * return e at least as precise as 'prec'
	 */
    {
	static real local_e = e;

	if (precision (local_e) < prec)
	    local_e = calculate_e (prec);
	return imprecise(local_e, prec);
    }

    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 expo = exponent (v);
	int scale;
	if (prec > 50)
	    scale = 27;
	else
	    scale = 12;
	if (expo + scale < 0)
	    scale = -expo;
	expo += scale;
	int div = (1 << scale);

	int iter = prec + 1;
	int iprec = prec + iter;

	real mant = imprecise (mantissa(v), iprec) / div;

	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 = 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);

	int	e = 0;

	/*
	 * Bring all values within the range of 1 <= a <= 2.
	 */
	if (a > 2) {
	    e = exponent(a) - 1;
	    a = mantissa(a) * 2;
	}

	/*
	 * 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++;
	}

	if (maxiter > 0)
	{
	    int	iprec = prec + maxiter * 16;

	    v = imprecise (v, iprec);
	    a = imprecise (a, iprec);

	    /*
	     * Newton's method
	     *
	     * v = log(a)
	     *
	     * exp(v) = a
	     *
	     * exp(v) - a = 0
	     *
	     * v' = v - (exp(v) - a) / exp(v)
	     *
	     * v' = v - 1 + a/exp(v)
	     *
	     *  v' = v - 1 + a * exp(-v);
	     */
	    real one = imprecise (1, iprec);

	    while (maxiter-- > 0)
		v = (v - one) + a / exp(v);
	}

	/* Mix in the log of the exponent */
	if (e != 0) {
	    int eprec = prec + 16;
	    static real log2 = 0;
	    if (log2 == 0 || precision(log2) < eprec)
		log2 = log(imprecise(2, eprec));
	    v += imprecise(e, eprec) * log2;
	}

	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);
    }

    /* Normalize angle to -π < aa <= π */
    real limit_angle_to_pi (real aa)
    {
	real	my_pi;

	aa = imprecise (aa);
	my_pi = pi_value (precision (aa));
	if (aa + my_pi == aa)
		raise invalid_argument ("argument not precise enough", 0, aa);
	aa %= 2 * my_pi;
	if (aa > my_pi)
		aa -= 2 * my_pi;
	return aa;
    }

    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) * 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 (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 (x == 0) {
		if (y == 0)
			return 0;
		if (y >= 0)
			return pi_value(precision(y))/2;
		else
			return -pi_value(precision(y))/2;
	}
	real a = atan(y/x);
	if (x < 0) {
		real p = pi_value(precision(y));
		if (y >= 0)
			a += p;
		else
			a -= p;
	}
	return a;
    }

    /*
     *	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) = atan2(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 atan2(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)))
     *		= atan2 (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 atan2 (sqrt (1-v**2), v);
    }

    /*
     * These two are used for the '**' and '**=' operators
     */
    public real pow (real a, real b)
	/*
	 * return a ** b;
	 */
    {
	real    result;
	if (a == 0) {
	    if (b == 0)
		return 1;
	    return 0;
	}

	if (is_int (b))
	{
	    if (!is_int (a) && is_rational (a))
		return pow (numerator(a), b) / pow (denominator (a), b);
	    bool flip = false;

	    if (b < 0)
	    {
		flip = true;
		b = -b;
	    }
	    result = 1;

	    int prec = precision(a);
	    /* Increase precision to avoid dropping bits */
	    if (prec != 0 && b > 0)
		a = imprecise(a, prec + (ilog(2, b) + 1) * 4 + 5);

	    while (b > 0)
	    {
		if ((b & 1) != 0)
		    result *= a;
		b >>= 1;
		if (b != 0)
		    a *= a;
	    }
	    if (flip)
		result = 1/result;
	    if (prec != 0)
		result = imprecise(result, prec);
	}
	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;
    }


    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;
    }

    public int choose(int n, int k)
        /* Number of ways of choosing k items from a set of
           n distinct items. */
    {
        /* This keeps the size of the intermediate terms
           smaller below, and also makes the bounds check
           slightly easier. */
        if (k > (n + 1) ⫽ 2)
            k = n - k;
        if (n < 0 || k < 0)
            return 0;
        int c = 1;
        /* This keeps the size of the intermediate terms
           down a bit compared to the traditional
           computation.  It probably doesn't matter, but oh
           well. */
        for (int i = n - k + 1; i <= n; i++)
            c *= i;
        return c / k!;
    }

    real(real) _abs = abs;
    public real(real) abs = _abs;
    int(real) _precision = precision;
    public int(real) precision = _precision;
}

/* XXX these shouldn't be here, but it was *convenient* */
&int(string, int ...) atoi = &string_to_integer;
&rational(string) atof = &string_to_real;