extend namespace Math {
    public real pi = imprecise
    (3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679,
    256);
    
    protected real e = imprecise 
    (2.7182818284590452353602874713526624977572470936999595749669676277240766303535475945713821785251664274,
     256);

    real function errbounds (real v)
    {
	int	prec;

	prec = precision (v);
	if (prec == 0)
	    prec = 256;
	return imprecise (2**(prec+3), prec + 3);
    }

    public real function sqrt (real v)
    {
	if (v < 0)
	    raise invalid_argument ("sqrt of negative number", 0, v);
	real function real_sqrt(real v)
	{
	    real	err;
	    real	prev, cur;
	    int	n, iter;

	    v = imprecise (v);
	    err = 1/(2**(precision (v)+3));
	    prev = imprecise (1 / (0.75 * 2**(floor (exponent(v)/2))), precision(v));
	    iter = precision (v) + 15;
	    while (iter--)
	    {
		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 function cbrt (real v)
    {
	real function real_cbrt(real v)
	{
	    real	err;
	    real	prev, cur;
	    int	n;
	    int	s;

	    s = sign (v);
	    v = imprecise (abs (v));
	    err = 1/(2**(precision (v)+3));
	    prev = imprecise (1 / (0.75 * 2**(floor (exponent(v)/3))), precision(v));
	    for (;;)
	    {
		cur = 1/3 * (2 * prev + v / prev**2);
		if (abs (cur - prev) < err)
		    break;
		prev = cur;
	    }
	    return s * abs (cur);
	}

	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 function exp (real v)
    {
	if (v < 0)
	    return 1/exp(-v);
	if (v == 0)
	    return 1;
	v = imprecise (v);
	int prec = precision (v);
	/*
	 * For values greater than one,
	 * use e ** floor (v) * exp (v - floor (v));
	 *
	 * Otherwise the power series takes a long time
	 * to converge
	 */
	if (v > 1)
	{
	    int	    p = floor (v);
	    static  real big_e = e;
	    static  int big_e_bits = precision (big_e);

	    if (big_e_bits < prec * 2)
	    {
		big_e_bits = prec * 2;
		big_e = exp (imprecise (1, big_e_bits));
	    }
	    real    my_e = imprecise (big_e, prec * 2);
	    return imprecise (exp (v - p) * (my_e ** p), prec);
	}
	
	real err = errbounds (v);
	real e = imprecise (0, prec*2);
	real num = imprecise (1, prec*2);
	real den = imprecise (1, prec*2);
	real loop = imprecise (0, prec*2);
	
	v = imprecise (v, prec*2);
	
	/*
	 * Traditional power series
	 *
	 *  exp(n) = 1 + n/1 + n**2/2! + n**3/3!
	 */
	int iter = prec;
	while (iter--)
	{
	    real term = num/den;
	    e = e + term;
	    if (abs (e) >= abs (term * err))
		break;
	    num *= v;
	    loop++;
	    den *= loop;
	}
	return imprecise (e, prec);
    }

    public real function log (real a)
    {
	/*
	 * y = x/(x+2)
	 *
	 * log((1+y)/(1-y)) = 2 * (y + y**3/3 + y**5/5 +...)
	 */
	real function log_near_1 (real a)
	{
	    real    err = errbounds (a);
	    real    v, term;
	    int	    i;
	    int	    iter;
	    int	    prec = precision(a);

	    a = imprecise (a / (a + 2), prec*2);
	    v = 0;
	    iter = bit_width (prec) + 15;
	    for (i = 1; --iter; i += 2)
	    {
		term = a**i/i;
		v = v + term;
		/* printf ("v %g term %g term * err %g\n", v, term,
			   term*err); */
		if (v > term * err)
		    break;
	    }
	    return imprecise (v*2, prec);
	}
	
	if (a == 1)
	    return 0;
	if (a <= 0)
	    raise invalid_argument ("log: must be positive", 0, a);
	if (abs (a-1) < 1e-4)
	    return log_near_1 (a-1);
	if (a < 1)
	    return -log(1/a);

	a = imprecise (a);

	int	prec = precision(a);
	int	maxiter = bit_width (prec * 2) + 15;
	int	iprec = prec + maxiter * 4;
        real    err = errbounds (a);
	
	/* printf ("maxiter %g\n", maxiter); */
	
	a = imprecise (a, iprec);
	
	/*
	 * estimate = lg(a) / lg(e)
	 */

	real v = imprecise ((exponent (a)-1) / 1.44269504088896, 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--)
	{
	    real term = a * exp (-v) - one;
	    v = v + term;
	    if (v > abs (term) * err)
		break;
	}

	return imprecise (v + epow, prec);
    }

    /*
     * log10(x) = log10(e) * log(x)
     *
     * log10(e) = log(e) / log(10) = 1/log(10)
     */
    public real function log10 (real a)
    {
	static real	loge = 0;

	a = imprecise (a);
	if (!loge || 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 function log2 (real a)
    {
	static real	loge = 0;

	a = imprecise (a);
	if (!loge || precision (loge) < precision (a))
	    loge = 1/log(imprecise (2, precision (a)));
	return loge * log(a);
    }
    
    /*
     * Calculate pi using the formula:
     *
     *  PI = 24*atan (1/8) + 8*atan (1/57) + 4*atan (1/239);
     */
    
    real function calculate_pi (int prec)
    {
	/*
	 * Estimate the number of digits available for
	 * the specified value (v) after a certain number of
	 * loops (p)
	 */
	real function 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 function 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 function 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--)
	    {
		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 function pi_value (int prec)
    {
	static real local_pi = pi;
    
	if (precision (local_pi) < prec)
	    local_pi = calculate_pi (prec);
	return imprecise (local_pi, prec);
    }
    
    real function 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 function sin (real a)
    {
	/*
	 * sin(x) = x - x**3/3! + x**5/5! ... 
	 */
	real function raw_sin (real a)
	{
	    real    err;
	    real    v, term;
	    real    a4, aj, ai;
	    int	    i, j;
	    int	    iter;
	    int	    prec;

	    err = errbounds (a);
	    prec = precision(a);
	    a = imprecise(a,prec*2);
	    i = 1;
	    j = 3;
	    a4 = a**4;
	    ai = a**i;
	    aj = a**j;
	    iter = prec + 8;
	    v = 0;
	    while (iter--)
	    {
		term = ai/i! - aj/j!;
		v += term;
		if (abs (v) > abs(term) * err)
		    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 function do_5x (real a)
	{
	    return 16 * a**5 - 20 * a**3 + 5 * a;
	}
	
	real function 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 function cos (real a)
    {
	/*
	 * cos(x) = 1 - x**2/2! + x**4/4! - x**6/6! ... 
	 */

	real function raw_cos (real a)
	{
	    real    err = errbounds (a);
	    real    v, term;
	    real    ai, aj, a4;
	    int	    i, j;
	    int	    iter;
	    int	    prec = precision(a);

	    a = imprecise(a, prec*2);
	    i = 0;
	    j = 2;
	    ai = 1;
	    aj = a**2;
	    a4 = a**4;
	    iter = prec + 8;
	    v = 0;
	    while (iter--)
	    {
		term = ai/i! - aj/j!;
		v += term;
		if (abs (v) > abs(term) * err)
		    break;
		ai *= a4;
		aj *= a4;
		i += 4;
		j += 4;
	    }
	    return imprecise (v + term);
	}

	/* cos(4x) = 8 * (cos**4(x) - cos**2(x)) + 1 */
	real function do_4x (real c)
	{
	    return 8 * (c**4 - c**2) + 1;
	}
	 
	real function 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 function cos_to_sin (real v)
    {
	return sqrt (1 - v**2);
    }
    
    public void function sin_cos (real a, *real sinp, *real cosp)
    {
	real	c, s;

	a = limit_angle_to_pi (a);
	c = cos (a);
	s = sign(a) * abs (cos_to_sin(c));
	if (cosp)
	    *cosp = c;
	if (sinp)
	    *sinp = s;
    }
    
    public real function tan (real a)
    {
	real	c, s;

	a = imprecise (a);
	sin_cos (a, &s, &c);
	return s/c;
    }
    
    public real function atan (real v)
    {
	/* 
	 * atan(x) = x - x**3/3 + x**5/5 - ...
	 */
	real function raw_atan (real v)
	{
	    real    err = errbounds (v);
	    real    a, term;
	    real    vi, vj, v4;
	    int	    i, j;
	    int	    iter;
	    int	    prec = precision(v);

	    v = imprecise (v, prec*2);
	    i = 1;
	    j = 3;
	    vi = v**i;
	    vj = v**j;
	    v4 = v**4;
	    a = 0;
	    iter = prec + 8;
	    while (iter--)
	    {
		term = vi/i - vj/j;
		a += term;
		if (abs (a) > abs (term) * err)
		    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 function asin (real v)
    {
	v = imprecise (v);
	if (abs (v) > 1)
	    return 1/0;
	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 function acos (real v)
    {
	v = imprecise(v);
	if (abs (v) > 1)
	    return 1/0;
	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 (x/y)
     */
    public real function atan2 (real x, real y)
    {
	x = imprecise (x);
	y = imprecise (y);
	if (y == 0)
	{
	    if (x == 0)
		return 0;
	    return sign(x) * pi_value (precision (x))/2;
	}
	if (x == 0)
	{
	    if (y > 0)
		return 0;
	    return pi_value (precision (y));
	}

	real	a;

	a = atan (x/y);
	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 function pow (real a, real b)
    {
	real    result;
	if (is_int (b))
	{
	    int flip = 0;

	    if (b < 0)
	    {
		flip = 1;
		b = -b;
	    }
	    result = 1;
	    while (b)
	    {
		if (b % 2)
		    result *= a;
		if (b //= 2)
		    a *= a;
	    }
	    if (flip)
		result = 1/result;
	}
	else
	{
	    result = exp (b * log(a));
	}
	return result;
    }

    public real function assign_pow (*real a, real b)
    {
	return *a = pow (*a, b);
    }
}

import Math;

extend namespace Math {
    
}