diff options
Diffstat (limited to 'libm/double/ellpj.c')
| -rw-r--r-- | libm/double/ellpj.c | 171 |
1 files changed, 0 insertions, 171 deletions
diff --git a/libm/double/ellpj.c b/libm/double/ellpj.c deleted file mode 100644 index 327fc56e8..000000000 --- a/libm/double/ellpj.c +++ /dev/null @@ -1,171 +0,0 @@ -/* ellpj.c - * - * Jacobian Elliptic Functions - * - * - * - * SYNOPSIS: - * - * double u, m, sn, cn, dn, phi; - * int ellpj(); - * - * ellpj( 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-9 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 - * DEC sn 1800 4.5e-16 8.7e-17 - * IEEE phi 10000 9.2e-16* 1.4e-16* - * IEEE sn 50000 4.1e-15 4.6e-16 - * IEEE cn 40000 3.6e-15 4.4e-16 - * IEEE dn 10000 1.3e-12 1.8e-14 - * - * Peak error observed in consistency check using addition - * theorem for sn(u+v) was 4e-16 (absolute). Also tested by - * the above relation to the incomplete elliptic integral. - * Accuracy deteriorates when u is large. - * - */ - -/* ellpj.c */ - - -/* -Cephes Math Library Release 2.8: June, 2000 -Copyright 1984, 1987, 2000 by Stephen L. Moshier -*/ - -#include <math.h> -#ifdef ANSIPROT -extern double sqrt ( double ); -extern double fabs ( double ); -extern double sin ( double ); -extern double cos ( double ); -extern double asin ( double ); -extern double tanh ( double ); -extern double sinh ( double ); -extern double cosh ( double ); -extern double atan ( double ); -extern double exp ( double ); -#else -double sqrt(), fabs(), sin(), cos(), asin(), tanh(); -double sinh(), cosh(), atan(), exp(); -#endif -extern double PIO2, MACHEP; - -int ellpj( u, m, sn, cn, dn, ph ) -double u, m; -double *sn, *cn, *dn, *ph; -{ -double ai, b, phi, t, twon; -double a[9], c[9]; -int i; - - -/* Check for special cases */ - -if( m < 0.0 || m > 1.0 ) - { - mtherr( "ellpj", DOMAIN ); - *sn = 0.0; - *cn = 0.0; - *ph = 0.0; - *dn = 0.0; - return(-1); - } -if( m < 1.0e-9 ) - { - t = sin(u); - b = cos(u); - ai = 0.25 * m * (u - t*b); - *sn = t - ai*b; - *cn = b + ai*t; - *ph = u - ai; - *dn = 1.0 - 0.5*m*t*t; - return(0); - } - -if( m >= 0.9999999999 ) - { - ai = 0.25 * (1.0-m); - b = cosh(u); - t = tanh(u); - phi = 1.0/b; - twon = b * sinh(u); - *sn = t + ai * (twon - u)/(b*b); - *ph = 2.0*atan(exp(u)) - PIO2 + 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.0; -b = sqrt(1.0 - m); -c[0] = sqrt(m); -twon = 1.0; -i = 0; - -while( fabs(c[i]/a[i]) > MACHEP ) - { - if( i > 7 ) - { - mtherr( "ellpj", OVERFLOW ); - goto done; - } - ai = a[i]; - ++i; - c[i] = ( ai - b )/2.0; - t = sqrt( ai * b ); - a[i] = ( ai + b )/2.0; - b = t; - twon *= 2.0; - } - -done: - -/* backward recurrence */ -phi = twon * a[i] * u; -do - { - t = c[i] * sin(phi) / a[i]; - b = phi; - phi = (asin(t) + phi)/2.0; - } -while( --i ); - -*sn = sin(phi); -t = cos(phi); -*cn = t; -*dn = t/cos(phi-b); -*ph = phi; -return(0); -} |