cephes/ellpk.c

   1 /*                                                      ellpk.c
   2  *
   3  *      Complete elliptic integral of the first kind
   4  *
   5  *
   6  *
   7  * SYNOPSIS:
   8  *
   9  * double m1, y, ellpk();
  10  *
  11  * y = ellpk( 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  *    DEC        0,1        16000       3.5e-17     1.1e-17
  46  *    IEEE       0,1        30000       2.5e-16     6.8e-17
  47  *
  48  * ERROR MESSAGES:
  49  *
  50  *   message         condition      value returned
  51  * ellpk domain       x<0, x>1           0.0
  52  *
  53  */
  54 \f
  55 /*                                                      ellpk.c */
  56
  57
  58 /*
  59 Cephes Math Library, Release 2.0:  April, 1987
  60 Copyright 1984, 1987 by Stephen L. Moshier
  61 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  62 */
  63
  64 #include "mconf.h"
  65 #include "cephes.h"
  66
  67 #ifdef DEC
  68 static unsigned short P[] =
  69 {
  70 0035020,0127576,0040430,0051544,
  71 0036025,0070136,0042703,0153716,
  72 0036402,0122614,0062555,0077777,
  73 0036441,0102130,0072334,0025172,
  74 0036341,0043320,0117242,0172076,
  75 0036312,0146456,0077242,0154141,
  76 0036420,0003467,0013727,0035407,
  77 0036564,0137263,0110651,0020237,
  78 0036775,0001330,0144056,0020305,
  79 0037305,0144137,0157521,0141734,
  80 0040261,0071027,0173721,0147572
  81 };
  82 static unsigned short Q[] =
  83 {
  84 0034366,0130371,0103453,0077633,
  85 0035557,0122745,0173515,0113016,
  86 0036302,0124470,0167304,0074473,
  87 0036575,0132403,0117226,0117576,
  88 0036703,0156271,0047124,0147733,
  89 0036766,0137465,0002053,0157312,
  90 0037031,0014423,0154274,0176515,
  91 0037107,0177747,0143216,0016145,
  92 0037217,0177777,0172621,0074000,
  93 0037377,0177777,0177776,0156435,
  94 0040000,0000000,0000000,0000000
  95 };
  96 static unsigned short ac1[] = {0040261,0071027,0173721,0147572};
  97 #define C1 (*(double *)ac1)
  98 #endif
  99
 100 #ifdef IBMPC
 101 static unsigned short P[] =
 102 {
 103 0x0a6d,0xc823,0x15ef,0x3f22,
 104 0x7afa,0xc8b8,0xae0b,0x3f62,
 105 0xb000,0x8cad,0x54b1,0x3f80,
 106 0x854f,0x0e9b,0x308b,0x3f84,
 107 0x5e88,0x13d4,0x28da,0x3f7c,
 108 0x5b0c,0xcfd4,0x59a5,0x3f79,
 109 0xe761,0xe2fa,0x00e6,0x3f82,
 110 0x2414,0x7235,0x97d6,0x3f8e,
 111 0xc419,0x1905,0xa05b,0x3f9f,
 112 0x387c,0xfbea,0xb90b,0x3fb8,
 113 0x39ef,0xfefa,0x2e42,0x3ff6
 114 };
 115 static unsigned short Q[] =
 116 {
 117 0x6ff3,0x30e5,0xd61f,0x3efe,
 118 0xb2c2,0xbee9,0xf4bc,0x3f4d,
 119 0x8f27,0x1dd8,0x5527,0x3f78,
 120 0xd3f0,0x73d2,0xb6a0,0x3f8f,
 121 0x99fb,0x29ca,0x7b97,0x3f98,
 122 0x7bd9,0xa085,0xd7e6,0x3f9e,
 123 0x9faa,0x7b17,0x2322,0x3fa3,
 124 0xc38d,0xf8d1,0xfffc,0x3fa8,
 125 0x2f00,0xfeb2,0xffff,0x3fb1,
 126 0xdba4,0xffff,0xffff,0x3fbf,
 127 0x0000,0x0000,0x0000,0x3fe0
 128 };
 129 static unsigned short ac1[] = {0x39ef,0xfefa,0x2e42,0x3ff6};
 130 #define C1 (*(double *)ac1)
 131 #endif
 132
 133 #ifdef MIEEE
 134 static unsigned short P[] =
 135 {
 136 0x3f22,0x15ef,0xc823,0x0a6d,
 137 0x3f62,0xae0b,0xc8b8,0x7afa,
 138 0x3f80,0x54b1,0x8cad,0xb000,
 139 0x3f84,0x308b,0x0e9b,0x854f,
 140 0x3f7c,0x28da,0x13d4,0x5e88,
 141 0x3f79,0x59a5,0xcfd4,0x5b0c,
 142 0x3f82,0x00e6,0xe2fa,0xe761,
 143 0x3f8e,0x97d6,0x7235,0x2414,
 144 0x3f9f,0xa05b,0x1905,0xc419,
 145 0x3fb8,0xb90b,0xfbea,0x387c,
 146 0x3ff6,0x2e42,0xfefa,0x39ef
 147 };
 148 static unsigned short Q[] =
 149 {
 150 0x3efe,0xd61f,0x30e5,0x6ff3,
 151 0x3f4d,0xf4bc,0xbee9,0xb2c2,
 152 0x3f78,0x5527,0x1dd8,0x8f27,
 153 0x3f8f,0xb6a0,0x73d2,0xd3f0,
 154 0x3f98,0x7b97,0x29ca,0x99fb,
 155 0x3f9e,0xd7e6,0xa085,0x7bd9,
 156 0x3fa3,0x2322,0x7b17,0x9faa,
 157 0x3fa8,0xfffc,0xf8d1,0xc38d,
 158 0x3fb1,0xffff,0xfeb2,0x2f00,
 159 0x3fbf,0xffff,0xffff,0xdba4,
 160 0x3fe0,0x0000,0x0000,0x0000
 161 };
 162 static unsigned short ac1[] = {
 163 0x3ff6,0x2e42,0xfefa,0x39ef
 164 };
 165 #define C1 (*(double *)ac1)
 166 #endif
 167
 168 #ifdef UNK
 169 static double P[] =
 170 {
 171  1.37982864606273237150E-4,
 172  2.28025724005875567385E-3,
 173  7.97404013220415179367E-3,
 174  9.85821379021226008714E-3,
 175  6.87489687449949877925E-3,
 176  6.18901033637687613229E-3,
 177  8.79078273952743772254E-3,
 178  1.49380448916805252718E-2,
 179  3.08851465246711995998E-2,
 180  9.65735902811690126535E-2,
 181  1.38629436111989062502E0
 182 };
 183
 184 static double Q[] =
 185 {
 186  2.94078955048598507511E-5,
 187  9.14184723865917226571E-4,
 188  5.94058303753167793257E-3,
 189  1.54850516649762399335E-2,
 190  2.39089602715924892727E-2,
 191  3.01204715227604046988E-2,
 192  3.73774314173823228969E-2,
 193  4.88280347570998239232E-2,
 194  7.03124996963957469739E-2,
 195  1.24999999999870820058E-1,
 196  4.99999999999999999821E-1
 197 };
 198 static double C1 = 1.3862943611198906188E0; /* log(4) */
 199 #endif
 200
 201 extern double MACHEP, MAXNUM;
 202
 203 double ellpk(x)
 204 double x;
 205 {
 206
 207 if( (x < 0.0) || (x > 1.0) )
 208         {
 209         mtherr( "ellpk", DOMAIN );
 210         return( 0.0 );
 211         }
 212
 213 if( x > MACHEP )
 214         {
 215         return( polevl(x,P,10) - log(x) * polevl(x,Q,10) );
 216         }
 217 else
 218         {
 219         if( x == 0.0 )
 220                 {
 221                 mtherr( "ellpk", SING );
 222                 return( MAXNUM );
 223                 }
 224         else
 225                 {
 226                 return( C1 - 0.5 * log(x) );
 227                 }
 228         }
 229 }