lib/libm/ellpef.c

   1 /*                                                      ellpef.c
   2  *
   3  *      Complete elliptic integral of the second kind
   4  *
   5  *
   6  *
   7  * SYNOPSIS:
   8  *
   9  * float m1, y, ellpef();
  10  *
  11  * y = ellpef( m1 );
  12  *
  13  *
  14  *
  15  * DESCRIPTION:
  16  *
  17  * Approximates the integral
  18  *
  19  *
  20  *            pi/2
  21  *             -
  22  *            | |                 2
  23  * E(m)  =    |    sqrt( 1 - m sin t ) dt
  24  *          | |
  25  *           -
  26  *            0
  27  *
  28  * Where m = 1 - m1, using the approximation
  29  *
  30  *      P(x)  -  x log x Q(x).
  31  *
  32  * Though there are no singularities, the argument m1 is used
  33  * rather than m for compatibility with ellpk().
  34  *
  35  * E(1) = 1; E(0) = pi/2.
  36  *
  37  *
  38  * ACCURACY:
  39  *
  40  *                      Relative error:
  41  * arithmetic   domain     # trials      peak         rms
  42  *    IEEE       0, 1       30000       1.1e-7      3.9e-8
  43  *
  44  *
  45  * ERROR MESSAGES:
  46  *
  47  *   message         condition      value returned
  48  * ellpef domain     x<0, x>1            0.0
  49  *
  50  */
  51 \f
  52 /*                                                      ellpe.c         */
  53
  54 /* Elliptic integral of second kind */
  55
  56 /*
  57 Cephes Math Library, Release 2.1:  February, 1989
  58 Copyright 1984, 1987, 1989 by Stephen L. Moshier
  59 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  60 */
  61
  62 #include "mconf.h"
  63
  64
  65 static float P[] = {
  66   1.53552577301013293365E-4,
  67   2.50888492163602060990E-3,
  68   8.68786816565889628429E-3,
  69   1.07350949056076193403E-2,
  70   7.77395492516787092951E-3,
  71   7.58395289413514708519E-3,
  72   1.15688436810574127319E-2,
  73   2.18317996015557253103E-2,
  74   5.68051945617860553470E-2,
  75   4.43147180560990850618E-1,
  76   1.00000000000000000299E0
  77 };
  78 static float Q[] = {
  79   3.27954898576485872656E-5,
  80   1.00962792679356715133E-3,
  81   6.50609489976927491433E-3,
  82   1.68862163993311317300E-2,
  83   2.61769742454493659583E-2,
  84   3.34833904888224918614E-2,
  85   4.27180926518931511717E-2,
  86   5.85936634471101055642E-2,
  87   9.37499997197644278445E-2,
  88   2.49999999999888314361E-1
  89 };
  90
  91 #ifdef ANSIC
  92 float polevlf(float, float *, int), logf(float);
  93 float ellpef( float xx)
  94 #else
  95 float polevlf(), logf();
  96 float ellpef(xx)
  97 double xx;
  98 #endif
  99 {
 100 float x;
 101
 102 x = xx;
 103 if( (x <= 0.0) || (x > 1.0) )
 104         {
 105         if( x == 0.0 )
 106                 return( 1.0 );
 107         mtherr( "ellpef", DOMAIN );
 108         return( 0.0 );
 109         }
 110 return( polevlf(x,P,10) - logf(x) * (x * polevlf(x,Q,9)) );
 111 }