diff options
Diffstat (limited to 'libm/double/ellik.c')
-rw-r--r-- | libm/double/ellik.c | 148 |
1 files changed, 0 insertions, 148 deletions
diff --git a/libm/double/ellik.c b/libm/double/ellik.c deleted file mode 100644 index 1c9053676..000000000 --- a/libm/double/ellik.c +++ /dev/null @@ -1,148 +0,0 @@ -/* ellik.c - * - * Incomplete elliptic integral of the first kind - * - * - * - * SYNOPSIS: - * - * double phi, m, y, ellik(); - * - * y = ellik( phi, m ); - * - * - * - * DESCRIPTION: - * - * Approximates the integral - * - * - * - * phi - * - - * | | - * | dt - * F(phi_\m) = | ------------------ - * | 2 - * | | sqrt( 1 - m sin t ) - * - - * 0 - * - * of amplitude phi and modulus m, using the arithmetic - - * geometric mean algorithm. - * - * - * - * - * ACCURACY: - * - * Tested at random points with m in [0, 1] and phi as indicated. - * - * Relative error: - * arithmetic domain # trials peak rms - * IEEE -10,10 200000 7.4e-16 1.0e-16 - * - * - */ - - -/* -Cephes Math Library Release 2.8: June, 2000 -Copyright 1984, 1987, 2000 by Stephen L. Moshier -*/ - -/* Incomplete elliptic integral of first kind */ - -#include <math.h> -#ifdef ANSIPROT -extern double sqrt ( double ); -extern double fabs ( double ); -extern double log ( double ); -extern double tan ( double ); -extern double atan ( double ); -extern double floor ( double ); -extern double ellpk ( double ); -double ellik ( double, double ); -#else -double sqrt(), fabs(), log(), tan(), atan(), floor(), ellpk(); -double ellik(); -#endif -extern double PI, PIO2, MACHEP, MAXNUM; - -double ellik( phi, m ) -double phi, m; -{ -double a, b, c, e, temp, t, K; -int d, mod, sign, npio2; - -if( m == 0.0 ) - return( phi ); -a = 1.0 - m; -if( a == 0.0 ) - { - if( fabs(phi) >= PIO2 ) - { - mtherr( "ellik", SING ); - return( MAXNUM ); - } - return( log( tan( (PIO2 + phi)/2.0 ) ) ); - } -npio2 = floor( phi/PIO2 ); -if( npio2 & 1 ) - npio2 += 1; -if( npio2 ) - { - K = ellpk( a ); - phi = phi - npio2 * PIO2; - } -else - K = 0.0; -if( phi < 0.0 ) - { - phi = -phi; - sign = -1; - } -else - sign = 0; -b = sqrt(a); -t = tan( phi ); -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); - if( npio2 == 0 ) - K = ellpk( a ); - temp = K - ellik( e, m ); - goto done; - } - } -a = 1.0; -c = sqrt(m); -d = 1; -mod = 0; - -while( fabs(c/a) > MACHEP ) - { - temp = b/a; - phi = phi + atan(t*temp) + mod * PI; - mod = (phi + 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; - } - -temp = (atan(t) + mod * PI)/(d * a); - -done: -if( sign < 0 ) - temp = -temp; -temp += npio2 * K; -return( temp ); -} |