summaryrefslogtreecommitdiff
path: root/libm/ldouble/ellpjl.c
diff options
context:
space:
mode:
authorEric Andersen <andersen@codepoet.org>2001-11-22 14:04:29 +0000
committerEric Andersen <andersen@codepoet.org>2001-11-22 14:04:29 +0000
commit7ce331c01ce6eb7b3f5c715a38a24359da9c6ee2 (patch)
tree3a7e8476e868ae15f4da1b7ce26b2db6f434468c /libm/ldouble/ellpjl.c
parentc117dd5fb183afb1a4790a6f6110d88704be6bf8 (diff)
Totally rework the math library, this time based on the MacOs X
math library (which is itself based on the math lib from FreeBSD). -Erik
Diffstat (limited to 'libm/ldouble/ellpjl.c')
-rw-r--r--libm/ldouble/ellpjl.c164
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);
-}