lib/libm/ellief.c

   1 /*                                                      ellief.c
   2  *
   3  *      Incomplete elliptic integral of the second kind
   4  *
   5  *
   6  *
   7  * SYNOPSIS:
   8  *
   9  * float phi, m, y, ellief();
  10  *
  11  * y = ellief( phi, m );
  12  *
  13  *
  14  *
  15  * DESCRIPTION:
  16  *
  17  * Approximates the integral
  18  *
  19  *
  20  *                phi
  21  *                 -
  22  *                | |
  23  *                |                   2
  24  * E(phi\m)  =    |    sqrt( 1 - m sin t ) dt
  25  *                |
  26  *              | |
  27  *               -
  28  *                0
  29  *
  30  * of amplitude phi and modulus m, using the arithmetic -
  31  * geometric mean algorithm.
  32  *
  33  *
  34  *
  35  * ACCURACY:
  36  *
  37  * Tested at random arguments with phi in [0, 2] and m in
  38  * [0, 1].
  39  *                      Relative error:
  40  * arithmetic   domain     # trials      peak         rms
  41  *    IEEE       0,2        10000       4.5e-7      7.4e-8
  42  *
  43  *
  44  */
  45 \f
  46
  47 /*
  48 Cephes Math Library Release 2.2:  July, 1992
  49 Copyright 1984, 1987, 1992 by Stephen L. Moshier
  50 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  51 */
  52
  53 /*      Incomplete elliptic integral of second kind     */
  54
  55 #include "mconf.h"
  56
  57 extern float PIF, PIO2F, MACHEPF;
  58
  59 #define fabsf(x) ( (x) < 0 ? -(x) : (x) )
  60
  61 #ifdef ANSIC
  62 float sqrtf(float), logf(float), sinf(float), tanf(float), atanf(float);
  63 float ellpef(float), ellpkf(float);
  64 #else
  65 float sqrtf(), logf(), sinf(), tanf(), atanf();
  66 float ellpef(), ellpkf();
  67 #endif
  68
  69
  70 #ifdef ANSIC
  71 float ellief( float phia, float ma )
  72 #else
  73 float ellief( phia, ma )
  74 double phia, ma;
  75 #endif
  76 {
  77 float phi, m, a, b, c, e, temp;
  78 float lphi, t;
  79 int d, mod;
  80
  81 phi = phia;
  82 m = ma;
  83 if( m == 0.0 )
  84         return( phi );
  85 if( m == 1.0 )
  86         return( sinf(phi) );
  87 lphi = phi;
  88 if( lphi < 0.0 )
  89         lphi = -lphi;
  90 a = 1.0;
  91 b = 1.0 - m;
  92 b = sqrtf(b);
  93 c = sqrtf(m);
  94 d = 1;
  95 e = 0.0;
  96 t = tanf( lphi );
  97 mod = (lphi + PIO2F)/PIF;
  98
  99 while( fabsf(c/a) > MACHEPF )
 100         {
 101         temp = b/a;
 102         lphi = lphi + atanf(t*temp) + mod * PIF;
 103         mod = (lphi + PIO2F)/PIF;
 104         t = t * ( 1.0 + temp )/( 1.0 - temp * t * t );
 105         c = 0.5 * ( a - b );
 106         temp = sqrtf( a * b );
 107         a = 0.5 * ( a + b );
 108         b = temp;
 109         d += d;
 110         e += c * sinf(lphi);
 111         }
 112
 113 b = 1.0 - m;
 114 temp = ellpef(b)/ellpkf(b);
 115 temp *= (atanf(t) + mod * PIF)/(d * a);
 116 temp += e;
 117 if( phi < 0.0 )
 118         temp = -temp;
 119 return( temp );
 120 }