lib/libm/ellpkf.c

   1 /*                                                      ellpkf.c
   2  *
   3  *      Complete elliptic integral of the first kind
   4  *
   5  *
   6  *
   7  * SYNOPSIS:
   8  *
   9  * float m1, y, ellpkf();
  10  *
  11  * y = ellpkf( m1 );
  12  *
  13  *
  14  *
  15  * DESCRIPTION:
  16  *
  17  * Approximates the integral
  18  *
  19  *
  20  *
  21  *            pi/2
  22  *             -
  23  *            | |
  24  *            |           dt
  25  * K(m)  =    |    ------------------
  26  *            |                   2
  27  *          | |    sqrt( 1 - m sin t )
  28  *           -
  29  *            0
  30  *
  31  * where m = 1 - m1, using the approximation
  32  *
  33  *     P(x)  -  log x Q(x).
  34  *
  35  * The argument m1 is used rather than m so that the logarithmic
  36  * singularity at m = 1 will be shifted to the origin; this
  37  * preserves maximum accuracy.
  38  *
  39  * K(0) = pi/2.
  40  *
  41  * ACCURACY:
  42  *
  43  *                      Relative error:
  44  * arithmetic   domain     # trials      peak         rms
  45  *    IEEE       0,1        30000       1.3e-7      3.4e-8
  46  *
  47  * ERROR MESSAGES:
  48  *
  49  *   message         condition      value returned
  50  * ellpkf domain      x<0, x>1           0.0
  51  *
  52  */
  53 \f
  54 /*                                                      ellpk.c */
  55
  56
  57 /*
  58 Cephes Math Library, Release 2.0:  April, 1987
  59 Copyright 1984, 1987 by Stephen L. Moshier
  60 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  61 */
  62
  63 #include "mconf.h"
  64
  65 static float P[] =
  66 {
  67  1.37982864606273237150E-4,
  68  2.28025724005875567385E-3,
  69  7.97404013220415179367E-3,
  70  9.85821379021226008714E-3,
  71  6.87489687449949877925E-3,
  72  6.18901033637687613229E-3,
  73  8.79078273952743772254E-3,
  74  1.49380448916805252718E-2,
  75  3.08851465246711995998E-2,
  76  9.65735902811690126535E-2,
  77  1.38629436111989062502E0
  78 };
  79
  80 static float Q[] =
  81 {
  82  2.94078955048598507511E-5,
  83  9.14184723865917226571E-4,
  84  5.94058303753167793257E-3,
  85  1.54850516649762399335E-2,
  86  2.39089602715924892727E-2,
  87  3.01204715227604046988E-2,
  88  3.73774314173823228969E-2,
  89  4.88280347570998239232E-2,
  90  7.03124996963957469739E-2,
  91  1.24999999999870820058E-1,
  92  4.99999999999999999821E-1
  93 };
  94 static float C1 = 1.3862943611198906188E0; /* log(4) */
  95
  96 extern float MACHEPF, MAXNUMF;
  97
  98 #ifdef ANSIC
  99 float polevlf(float, float *, int);
 100 float p1evlf(float, float *, int);
 101 float logf(float);
 102 float ellpkf(float xx)
 103 #else
 104 float polevlf(), p1evlf(), logf();
 105 float ellpkf(xx)
 106 double xx;
 107 #endif
 108 {
 109 float x;
 110
 111 x = xx;
 112 if( (x < 0.0) || (x > 1.0) )
 113         {
 114         mtherr( "ellpkf", DOMAIN );
 115         return( 0.0 );
 116         }
 117
 118 if( x > MACHEPF )
 119         {
 120         return( polevlf(x,P,10) - logf(x) * polevlf(x,Q,10) );
 121         }
 122 else
 123         {
 124         if( x == 0.0 )
 125                 {
 126                 mtherr( "ellpkf", SING );
 127                 return( MAXNUMF );
 128                 }
 129         else
 130                 {
 131                 return( C1 - 0.5 * logf(x) );
 132                 }
 133         }
 134 }