/* elliel.c * * Incomplete elliptic integral of the second kind * * * * SYNOPSIS: * * long double phi, m, y, elliel(); * * y = elliel( 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 * IEEE -10,10 50000 2.7e-18 2.3e-19 * * */ /* Cephes Math Library Release 2.3: November, 1995 Copyright 1984, 1987, 1993, 1995 by Stephen L. Moshier */ /* Incomplete elliptic integral of second kind */ #include <math.h> #ifdef ANSIPROT extern long double sqrtl ( long double ); extern long double fabsl ( long double ); extern long double logl ( long double ); extern long double sinl ( long double ); extern long double tanl ( long double ); extern long double atanl ( long double ); extern long double floorl ( long double ); extern long double ellpel ( long double ); extern long double ellpkl ( long double ); long double elliel ( long double, long double ); #else long double sqrtl(), fabsl(), logl(), sinl(), tanl(), atanl(), floorl(); long double ellpel(), ellpkl(), elliel(); #endif extern long double PIL, PIO2L, MACHEPL; long double elliel( phi, m ) long double phi, m; { long double a, b, c, e, temp, lphi, t, E; int d, mod, npio2, sign; if( m == 0.0L ) return( phi ); lphi = phi; npio2 = floorl( lphi/PIO2L ); if( npio2 & 1 ) npio2 += 1; lphi = lphi - npio2 * PIO2L; if( lphi < 0.0L ) { lphi = -lphi; sign = -1; } else { sign = 1; } a = 1.0L - m; E = ellpel( a ); if( a == 0.0L ) { temp = sinl( lphi ); goto done; } t = tanl( lphi ); b = sqrtl(a); if( fabsl(t) > 10.0L ) { /* Transform the amplitude */ e = 1.0L/(b*t); /* ... but avoid multiple recursions. */ if( fabsl(e) < 10.0L ) { e = atanl(e); temp = E + m * sinl( lphi ) * sinl( e ) - elliel( e, m ); goto done; } } c = sqrtl(m); a = 1.0L; d = 1; e = 0.0L; mod = 0; while( fabsl(c/a) > MACHEPL ) { temp = b/a; lphi = lphi + atanl(t*temp) + mod * PIL; mod = (lphi + PIO2L)/PIL; t = t * ( 1.0L + temp )/( 1.0L - temp * t * t ); c = 0.5L*( a - b ); temp = sqrtl( a * b ); a = 0.5L*( a + b ); b = temp; d += d; e += c * sinl(lphi); } temp = E / ellpkl( 1.0L - m ); temp *= (atanl(t) + mod * PIL)/(d * a); temp += e; done: if( sign < 0 ) temp = -temp; temp += npio2 * E; return( temp ); }