diff options
Diffstat (limited to 'libm/ldouble/ellpjl.c')
| -rw-r--r-- | libm/ldouble/ellpjl.c | 164 |
1 files changed, 0 insertions, 164 deletions
diff --git a/libm/ldouble/ellpjl.c b/libm/ldouble/ellpjl.c deleted file mode 100644 index bb57fe6a1..000000000 --- a/libm/ldouble/ellpjl.c +++ /dev/null @@ -1,164 +0,0 @@ -/* ellpjl.c - * - * Jacobian Elliptic Functions - * - * - * - * SYNOPSIS: - * - * long double u, m, sn, cn, dn, phi; - * int ellpjl(); - * - * ellpjl( u, m, _&sn, _&cn, _&dn, _&phi ); - * - * - * - * DESCRIPTION: - * - * - * Evaluates the Jacobian elliptic functions sn(u|m), cn(u|m), - * and dn(u|m) of parameter m between 0 and 1, and real - * argument u. - * - * These functions are periodic, with quarter-period on the - * real axis equal to the complete elliptic integral - * ellpk(1.0-m). - * - * Relation to incomplete elliptic integral: - * If u = ellik(phi,m), then sn(u|m) = sin(phi), - * and cn(u|m) = cos(phi). Phi is called the amplitude of u. - * - * Computation is by means of the arithmetic-geometric mean - * algorithm, except when m is within 1e-12 of 0 or 1. In the - * latter case with m close to 1, the approximation applies - * only for phi < pi/2. - * - * ACCURACY: - * - * Tested at random points with u between 0 and 10, m between - * 0 and 1. - * - * Absolute error (* = relative error): - * arithmetic function # trials peak rms - * IEEE sn 10000 1.7e-18 2.3e-19 - * IEEE cn 20000 1.6e-18 2.2e-19 - * IEEE dn 10000 4.7e-15 2.7e-17 - * IEEE phi 10000 4.0e-19* 6.6e-20* - * - * Accuracy deteriorates when u is large. - * - */ - -/* -Cephes Math Library Release 2.3: November, 1995 -Copyright 1984, 1987, 1995 by Stephen L. Moshier -*/ - -#include <math.h> -#ifdef ANSIPROT -extern long double sqrtl ( long double ); -extern long double fabsl ( long double ); -extern long double sinl ( long double ); -extern long double cosl ( long double ); -extern long double asinl ( long double ); -extern long double tanhl ( long double ); -extern long double sinhl ( long double ); -extern long double coshl ( long double ); -extern long double atanl ( long double ); -extern long double expl ( long double ); -#else -long double sqrtl(), fabsl(), sinl(), cosl(), asinl(), tanhl(); -long double sinhl(), coshl(), atanl(), expl(); -#endif -extern long double PIO2L, MACHEPL; - -int ellpjl( u, m, sn, cn, dn, ph ) -long double u, m; -long double *sn, *cn, *dn, *ph; -{ -long double ai, b, phi, t, twon; -long double a[9], c[9]; -int i; - - -/* Check for special cases */ - -if( m < 0.0L || m > 1.0L ) - { - mtherr( "ellpjl", DOMAIN ); - *sn = 0.0L; - *cn = 0.0L; - *ph = 0.0L; - *dn = 0.0L; - return(-1); - } -if( m < 1.0e-12L ) - { - t = sinl(u); - b = cosl(u); - ai = 0.25L * m * (u - t*b); - *sn = t - ai*b; - *cn = b + ai*t; - *ph = u - ai; - *dn = 1.0L - 0.5L*m*t*t; - return(0); - } - -if( m >= 0.999999999999L ) - { - ai = 0.25L * (1.0L-m); - b = coshl(u); - t = tanhl(u); - phi = 1.0L/b; - twon = b * sinhl(u); - *sn = t + ai * (twon - u)/(b*b); - *ph = 2.0L*atanl(expl(u)) - PIO2L + ai*(twon - u)/b; - ai *= t * phi; - *cn = phi - ai * (twon - u); - *dn = phi + ai * (twon + u); - return(0); - } - - -/* A. G. M. scale */ -a[0] = 1.0L; -b = sqrtl(1.0L - m); -c[0] = sqrtl(m); -twon = 1.0L; -i = 0; - -while( fabsl(c[i]/a[i]) > MACHEPL ) - { - if( i > 7 ) - { - mtherr( "ellpjl", OVERFLOW ); - goto done; - } - ai = a[i]; - ++i; - c[i] = 0.5L * ( ai - b ); - t = sqrtl( ai * b ); - a[i] = 0.5L * ( ai + b ); - b = t; - twon *= 2.0L; - } - -done: - -/* backward recurrence */ -phi = twon * a[i] * u; -do - { - t = c[i] * sinl(phi) / a[i]; - b = phi; - phi = 0.5L * (asinl(t) + phi); - } -while( --i ); - -*sn = sinl(phi); -t = cosl(phi); -*cn = t; -*dn = t/cosl(phi-b); -*ph = phi; -return(0); -} |