diff options
Diffstat (limited to 'libm/float/ellpjf.c')
-rw-r--r-- | libm/float/ellpjf.c | 161 |
1 files changed, 0 insertions, 161 deletions
diff --git a/libm/float/ellpjf.c b/libm/float/ellpjf.c deleted file mode 100644 index 552f5ffe4..000000000 --- a/libm/float/ellpjf.c +++ /dev/null @@ -1,161 +0,0 @@ -/* ellpjf.c - * - * Jacobian Elliptic Functions - * - * - * - * SYNOPSIS: - * - * float 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 - * IEEE sn 10000 1.7e-6 2.2e-7 - * IEEE cn 10000 1.6e-6 2.2e-7 - * IEEE dn 10000 1.4e-3 1.9e-5 - * IEEE phi 10000 3.9e-7* 6.7e-8* - * - * 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.2: July, 1992 -Copyright 1984, 1987, 1992 by Stephen L. Moshier -Direct inquiries to 30 Frost Street, Cambridge, MA 02140 -*/ - -#include <math.h> -extern float PIO2F, MACHEPF; - -#define fabsf(x) ( (x) < 0 ? -(x) : (x) ) - -#ifdef ANSIC -float sqrtf(float), sinf(float), cosf(float), asinf(float), tanhf(float); -float sinhf(float), coshf(float), atanf(float), expf(float); -#else -float sqrtf(), sinf(), cosf(), asinf(), tanhf(); -float sinhf(), coshf(), atanf(), expf(); -#endif - -int ellpjf( float uu, float mm, - float *sn, float *cn, float *dn, float *ph ) -{ -float u, m, ai, b, phi, t, twon; -float a[10], c[10]; -int i; - -u = uu; -m = mm; -/* Check for special cases */ - -if( m < 0.0 || m > 1.0 ) - { - mtherr( "ellpjf", DOMAIN ); - return(-1); - } -if( m < 1.0e-5 ) - { - t = sinf(u); - b = cosf(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.99999 ) - { - ai = 0.25 * (1.0-m); - b = coshf(u); - t = tanhf(u); - phi = 1.0/b; - twon = b * sinhf(u); - *sn = t + ai * (twon - u)/(b*b); - *ph = 2.0*atanf(expf(u)) - PIO2F + 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 = sqrtf(1.0 - m); -c[0] = sqrtf(m); -twon = 1.0; -i = 0; - -while( fabsf( (c[i]/a[i]) ) > MACHEPF ) - { - if( i > 8 ) - { -/* mtherr( "ellpjf", OVERFLOW );*/ - break; - } - ai = a[i]; - ++i; - c[i] = 0.5 * ( ai - b ); - t = sqrtf( ai * b ); - a[i] = 0.5 * ( ai + b ); - b = t; - twon += twon; - } - - -/* backward recurrence */ -phi = twon * a[i] * u; -do - { - t = c[i] * sinf(phi) / a[i]; - b = phi; - phi = 0.5 * (asinf(t) + phi); - } -while( --i ); - -*sn = sinf(phi); -t = cosf(phi); -*cn = t; -*dn = t/cosf(phi-b); -*ph = phi; -return(0); -} |