lib/libm/ellpjf.c

   1 /*                                                      ellpjf.c
   2  *
   3  *      Jacobian Elliptic Functions
   4  *
   5  *
   6  *
   7  * SYNOPSIS:
   8  *
   9  * float 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  *    IEEE      sn          10000       1.7e-6      2.2e-7
  44  *    IEEE      cn          10000       1.6e-6      2.2e-7
  45  *    IEEE      dn          10000       1.4e-3      1.9e-5
  46  *    IEEE      phi         10000       3.9e-7*     6.7e-8*
  47  *
  48  *  Peak error observed in consistency check using addition
  49  * theorem for sn(u+v) was 4e-16 (absolute).  Also tested by
  50  * the above relation to the incomplete elliptic integral.
  51  * Accuracy deteriorates when u is large.
  52  *
  53  */
  54 \f
  55 /*                                                      ellpj.c         */
  56
  57
  58 /*
  59 Cephes Math Library Release 2.2:  July, 1992
  60 Copyright 1984, 1987, 1992 by Stephen L. Moshier
  61 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  62 */
  63
  64 #include "mconf.h"
  65 extern double PIO2F, MACHEPF;
  66
  67 #define fabsf(x) ( (x) < 0 ? -(x) : (x) )
  68
  69 #ifdef ANSIC
  70 float sqrtf(float), sinf(float), cosf(float), asinf(float), tanhf(float);
  71 float sinhf(float), coshf(float), atanf(float), expf(float);
  72 #else
  73 float sqrtf(), sinf(), cosf(), asinf(), tanhf();
  74 float sinhf(), coshf(), atanf(), expf();
  75 #endif
  76
  77 #ifdef ANSIC
  78 int ellpjf( float uu, float mm,
  79    float *sn, float *cn, float *dn, float *ph )
  80 #else
  81 int ellpjf( uu, mm, sn, cn, dn, ph )
  82 double uu, mm;
  83 float *sn, *cn, *dn, *ph;
  84 #endif
  85 {
  86 float u, m, ai, b, phi, t, twon;
  87 float a[10], c[10];
  88 int i;
  89
  90 u = uu;
  91 m = mm;
  92 /* Check for special cases */
  93
  94 if( m < 0.0 || m > 1.0 )
  95         {
  96         mtherr( "ellpjf", DOMAIN );
  97         return(-1);
  98         }
  99 if( m < 1.0e-5 )
 100         {
 101         t = sinf(u);
 102         b = cosf(u);
 103         ai = 0.25 * m * (u - t*b);
 104         *sn = t - ai*b;
 105         *cn = b + ai*t;
 106         *ph = u - ai;
 107         *dn = 1.0 - 0.5*m*t*t;
 108         return(0);
 109         }
 110
 111 if( m >= 0.99999 )
 112         {
 113         ai = 0.25 * (1.0-m);
 114         b = coshf(u);
 115         t = tanhf(u);
 116         phi = 1.0/b;
 117         twon = b * sinhf(u);
 118         *sn = t + ai * (twon - u)/(b*b);
 119         *ph = 2.0*atanf(expf(u)) - PIO2F + ai*(twon - u)/b;
 120         ai *= t * phi;
 121         *cn = phi - ai * (twon - u);
 122         *dn = phi + ai * (twon + u);
 123         return(0);
 124         }
 125
 126
 127 /*      A. G. M. scale          */
 128 a[0] = 1.0;
 129 b = sqrtf(1.0 - m);
 130 c[0] = sqrtf(m);
 131 twon = 1.0;
 132 i = 0;
 133
 134 while( fabsf( (c[i]/a[i]) ) > MACHEPF )
 135         {
 136         if( i > 8 )
 137                 {
 138 /*              mtherr( "ellpjf", OVERFLOW );*/
 139                 break;
 140                 }
 141         ai = a[i];
 142         ++i;
 143         c[i] = 0.5 * ( ai - b );
 144         t = sqrtf( ai * b );
 145         a[i] = 0.5 * ( ai + b );
 146         b = t;
 147         twon += twon;
 148         }
 149
 150
 151 /* backward recurrence */
 152 phi = twon * a[i] * u;
 153 do
 154         {
 155         t = c[i] * sinf(phi) / a[i];
 156         b = phi;
 157         phi = 0.5 * (asinf(t) + phi);
 158         }
 159 while( --i );
 160
 161 *sn = sinf(phi);
 162 t = cosf(phi);
 163 *cn = t;
 164 *dn = t/cosf(phi-b);
 165 *ph = phi;
 166 return(0);
 167 }