From 7ce331c01ce6eb7b3f5c715a38a24359da9c6ee2 Mon Sep 17 00:00:00 2001 From: Eric Andersen Date: Thu, 22 Nov 2001 14:04:29 +0000 Subject: Totally rework the math library, this time based on the MacOs X math library (which is itself based on the math lib from FreeBSD). -Erik --- libm/double/ellie.c | 148 ---------------------------------------------------- 1 file changed, 148 deletions(-) delete mode 100644 libm/double/ellie.c (limited to 'libm/double/ellie.c') diff --git a/libm/double/ellie.c b/libm/double/ellie.c deleted file mode 100644 index 4f3379aa6..000000000 --- a/libm/double/ellie.c +++ /dev/null @@ -1,148 +0,0 @@ -/* ellie.c - * - * Incomplete elliptic integral of the second kind - * - * - * - * SYNOPSIS: - * - * double phi, m, y, ellie(); - * - * y = ellie( phi, m ); - * - * - * - * DESCRIPTION: - * - * Approximates the integral - * - * - * phi - * - - * | | - * | 2 - * E(phi_\m) = | sqrt( 1 - m sin t ) dt - * | - * | | - * - - * 0 - * - * of amplitude phi and modulus m, using the arithmetic - - * geometric mean algorithm. - * - * - * - * ACCURACY: - * - * Tested at random arguments with phi in [-10, 10] and m in - * [0, 1]. - * Relative error: - * arithmetic domain # trials peak rms - * DEC 0,2 2000 1.9e-16 3.4e-17 - * IEEE -10,10 150000 3.3e-15 1.4e-16 - * - * - */ - - -/* -Cephes Math Library Release 2.8: June, 2000 -Copyright 1984, 1987, 1993, 2000 by Stephen L. Moshier -*/ - -/* Incomplete elliptic integral of second kind */ -#include -extern double PI, PIO2, MACHEP; -#ifdef ANSIPROT -extern double sqrt ( double ); -extern double fabs ( double ); -extern double log ( double ); -extern double sin ( double x ); -extern double tan ( double x ); -extern double atan ( double ); -extern double floor ( double ); -extern double ellpe ( double ); -extern double ellpk ( double ); -double ellie ( double, double ); -#else -double sqrt(), fabs(), log(), sin(), tan(), atan(), floor(); -double ellpe(), ellpk(), ellie(); -#endif - -double ellie( phi, m ) -double phi, m; -{ -double a, b, c, e, temp; -double lphi, t, E; -int d, mod, npio2, sign; - -if( m == 0.0 ) - return( phi ); -lphi = phi; -npio2 = floor( lphi/PIO2 ); -if( npio2 & 1 ) - npio2 += 1; -lphi = lphi - npio2 * PIO2; -if( lphi < 0.0 ) - { - lphi = -lphi; - sign = -1; - } -else - { - sign = 1; - } -a = 1.0 - m; -E = ellpe( a ); -if( a == 0.0 ) - { - temp = sin( lphi ); - goto done; - } -t = tan( lphi ); -b = sqrt(a); -/* Thanks to Brian Fitzgerald - for pointing out an instability near odd multiples of pi/2. */ -if( fabs(t) > 10.0 ) - { - /* Transform the amplitude */ - e = 1.0/(b*t); - /* ... but avoid multiple recursions. */ - if( fabs(e) < 10.0 ) - { - e = atan(e); - temp = E + m * sin( lphi ) * sin( e ) - ellie( e, m ); - goto done; - } - } -c = sqrt(m); -a = 1.0; -d = 1; -e = 0.0; -mod = 0; - -while( fabs(c/a) > MACHEP ) - { - temp = b/a; - lphi = lphi + atan(t*temp) + mod * PI; - mod = (lphi + PIO2)/PI; - t = t * ( 1.0 + temp )/( 1.0 - temp * t * t ); - c = ( a - b )/2.0; - temp = sqrt( a * b ); - a = ( a + b )/2.0; - b = temp; - d += d; - e += c * sin(lphi); - } - -temp = E / ellpk( 1.0 - m ); -temp *= (atan(t) + mod * PI)/(d * a); -temp += e; - -done: - -if( sign < 0 ) - temp = -temp; -temp += npio2 * E; -return( temp ); -} -- cgit v1.2.3