lib/libm/ellpj.c

   1 /*                                                      ellpj.c
   2  *
   3  *      Jacobian Elliptic Functions
   4  *
   5  *
   6  *
   7  * SYNOPSIS:
   8  *
   9  * double u, m, sn, cn, dn, phi;
  10  * int ellpj();
  11  *
  12  * ellpj( u, m, _&sn, _&cn, _&dn, _&phi );
  13  *
  14  *
  15  *
  16  * DESCRIPTION:
  17  *
  18  *
  19  * Evaluates the Jacobian elliptic functions sn(u|m), cn(u|m),
  20  * and dn(u|m) of parameter m between 0 and 1, and real
  21  * argument u.
  22  *
  23  * These functions are periodic, with quarter-period on the
  24  * real axis equal to the complete elliptic integral
  25  * ellpk(1.0-m).
  26  *
  27  * Relation to incomplete elliptic integral:
  28  * If u = ellik(phi,m), then sn(u|m) = sin(phi),
  29  * and cn(u|m) = cos(phi).  Phi is called the amplitude of u.
  30  *
  31  * Computation is by means of the arithmetic-geometric mean
  32  * algorithm, except when m is within 1e-9 of 0 or 1.  In the
  33  * latter case with m close to 1, the approximation applies
  34  * only for phi < pi/2.
  35  *
  36  * ACCURACY:
  37  *
  38  * Tested at random points with u between 0 and 10, m between
  39  * 0 and 1.
  40  *
  41  *            Absolute error (* = relative error):
  42  * arithmetic   function   # trials      peak         rms
  43  *    DEC       sn           1800       4.5e-16     8.7e-17
  44  *    IEEE      phi         10000       9.2e-16*    1.4e-16*
  45  *    IEEE      sn          50000       4.1e-15     4.6e-16
  46  *    IEEE      cn          40000       3.6e-15     4.4e-16
  47  *    IEEE      dn          10000       1.3e-12     1.8e-14
  48  *
  49  *  Peak error observed in consistency check using addition
  50  * theorem for sn(u+v) was 4e-16 (absolute).  Also tested by
  51  * the above relation to the incomplete elliptic integral.
  52  * Accuracy deteriorates when u is large.
  53  *
  54  */
  55 \f
  56 /*                                                      ellpj.c         */
  57
  58
  59 /*
  60 Cephes Math Library Release 2.8:  June, 2000
  61 Copyright 1984, 1987, 2000 by Stephen L. Moshier
  62 */
  63
  64 #include "mconf.h"
  65 #ifdef ANSIPROT
  66 extern double sqrt ( double );
  67 extern double fabs ( double );
  68 extern double sin ( double );
  69 extern double cos ( double );
  70 extern double asin ( double );
  71 extern double tanh ( double );
  72 extern double sinh ( double );
  73 extern double cosh ( double );
  74 extern double atan ( double );
  75 extern double exp ( double );
  76 #else
  77 double sqrt(), fabs(), sin(), cos(), asin(), tanh();
  78 double sinh(), cosh(), atan(), exp();
  79 #endif
  80 extern double PIO2, MACHEP;
  81
  82 int ellpj( u, m, sn, cn, dn, ph )
  83 double u, m;
  84 double *sn, *cn, *dn, *ph;
  85 {
  86 double ai, b, phi, t, twon;
  87 double a[9], c[9];
  88 int i;
  89
  90
  91 /* Check for special cases */
  92
  93 if( m < 0.0 || m > 1.0 )
  94         {
  95         mtherr( "ellpj", DOMAIN );
  96         *sn = 0.0;
  97         *cn = 0.0;
  98         *ph = 0.0;
  99         *dn = 0.0;
 100         return(-1);
 101         }
 102 if( m < 1.0e-9 )
 103         {
 104         t = sin(u);
 105         b = cos(u);
 106         ai = 0.25 * m * (u - t*b);
 107         *sn = t - ai*b;
 108         *cn = b + ai*t;
 109         *ph = u - ai;
 110         *dn = 1.0 - 0.5*m*t*t;
 111         return(0);
 112         }
 113
 114 if( m >= 0.9999999999 )
 115         {
 116         ai = 0.25 * (1.0-m);
 117         b = cosh(u);
 118         t = tanh(u);
 119         phi = 1.0/b;
 120         twon = b * sinh(u);
 121         *sn = t + ai * (twon - u)/(b*b);
 122         *ph = 2.0*atan(exp(u)) - PIO2 + ai*(twon - u)/b;
 123         ai *= t * phi;
 124         *cn = phi - ai * (twon - u);
 125         *dn = phi + ai * (twon + u);
 126         return(0);
 127         }
 128
 129
 130 /*      A. G. M. scale          */
 131 a[0] = 1.0;
 132 b = sqrt(1.0 - m);
 133 c[0] = sqrt(m);
 134 twon = 1.0;
 135 i = 0;
 136
 137 while( fabs(c[i]/a[i]) > MACHEP )
 138         {
 139         if( i > 7 )
 140                 {
 141                 mtherr( "ellpj", OVERFLOW );
 142                 goto done;
 143                 }
 144         ai = a[i];
 145         ++i;
 146         c[i] = ( ai - b )/2.0;
 147         t = sqrt( ai * b );
 148         a[i] = ( ai + b )/2.0;
 149         b = t;
 150         twon *= 2.0;
 151         }
 152
 153 done:
 154
 155 /* backward recurrence */
 156 phi = twon * a[i] * u;
 157 do
 158         {
 159         t = c[i] * sin(phi) / a[i];
 160         b = phi;
 161         phi = (asin(t) + phi)/2.0;
 162         }
 163 while( --i );
 164
 165 *sn = sin(phi);
 166 t = cos(phi);
 167 *cn = t;
 168 *dn = t/cos(phi-b);
 169 *ph = phi;
 170 return(0);
 171 }