home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD2.mdf / c / library / dos / math / cephes / ellf / ellpk.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-11-17  |  4.7 KB  |  231 lines
  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.  
  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.  
  66. #ifdef DEC
  67. static short P[] =
  68. {
  69. 0035020,0127576,0040430,0051544,
  70. 0036025,0070136,0042703,0153716,
  71. 0036402,0122614,0062555,0077777,
  72. 0036441,0102130,0072334,0025172,
  73. 0036341,0043320,0117242,0172076,
  74. 0036312,0146456,0077242,0154141,
  75. 0036420,0003467,0013727,0035407,
  76. 0036564,0137263,0110651,0020237,
  77. 0036775,0001330,0144056,0020305,
  78. 0037305,0144137,0157521,0141734,
  79. 0040261,0071027,0173721,0147572
  80. };
  81. static short Q[] =
  82. {
  83. 0034366,0130371,0103453,0077633,
  84. 0035557,0122745,0173515,0113016,
  85. 0036302,0124470,0167304,0074473,
  86. 0036575,0132403,0117226,0117576,
  87. 0036703,0156271,0047124,0147733,
  88. 0036766,0137465,0002053,0157312,
  89. 0037031,0014423,0154274,0176515,
  90. 0037107,0177747,0143216,0016145,
  91. 0037217,0177777,0172621,0074000,
  92. 0037377,0177777,0177776,0156435,
  93. 0040000,0000000,0000000,0000000
  94. };
  95. static short ac1[] = {0040261,0071027,0173721,0147572};
  96. #define C1 (*(double *)ac1)
  97. #endif
  98.  
  99. #ifdef IBMPC
  100. static short P[] =
  101. {
  102. 0x0a6d,0xc823,0x15ef,0x3f22,
  103. 0x7afa,0xc8b8,0xae0b,0x3f62,
  104. 0xb000,0x8cad,0x54b1,0x3f80,
  105. 0x854f,0x0e9b,0x308b,0x3f84,
  106. 0x5e88,0x13d4,0x28da,0x3f7c,
  107. 0x5b0c,0xcfd4,0x59a5,0x3f79,
  108. 0xe761,0xe2fa,0x00e6,0x3f82,
  109. 0x2414,0x7235,0x97d6,0x3f8e,
  110. 0xc419,0x1905,0xa05b,0x3f9f,
  111. 0x387c,0xfbea,0xb90b,0x3fb8,
  112. 0x39ef,0xfefa,0x2e42,0x3ff6
  113. };
  114. static short Q[] =
  115. {
  116. 0x6ff3,0x30e5,0xd61f,0x3efe,
  117. 0xb2c2,0xbee9,0xf4bc,0x3f4d,
  118. 0x8f27,0x1dd8,0x5527,0x3f78,
  119. 0xd3f0,0x73d2,0xb6a0,0x3f8f,
  120. 0x99fb,0x29ca,0x7b97,0x3f98,
  121. 0x7bd9,0xa085,0xd7e6,0x3f9e,
  122. 0x9faa,0x7b17,0x2322,0x3fa3,
  123. 0xc38d,0xf8d1,0xfffc,0x3fa8,
  124. 0x2f00,0xfeb2,0xffff,0x3fb1,
  125. 0xdba4,0xffff,0xffff,0x3fbf,
  126. 0x0000,0x0000,0x0000,0x3fe0
  127. };
  128. static short ac1[] = {0x39ef,0xfefa,0x2e42,0x3ff6};
  129. #define C1 (*(double *)ac1)
  130. #endif
  131.  
  132. #ifdef MIEEE
  133. static short P[] =
  134. {
  135. 0x3f22,0x15ef,0xc823,0x0a6d,
  136. 0x3f62,0xae0b,0xc8b8,0x7afa,
  137. 0x3f80,0x54b1,0x8cad,0xb000,
  138. 0x3f84,0x308b,0x0e9b,0x854f,
  139. 0x3f7c,0x28da,0x13d4,0x5e88,
  140. 0x3f79,0x59a5,0xcfd4,0x5b0c,
  141. 0x3f82,0x00e6,0xe2fa,0xe761,
  142. 0x3f8e,0x97d6,0x7235,0x2414,
  143. 0x3f9f,0xa05b,0x1905,0xc419,
  144. 0x3fb8,0xb90b,0xfbea,0x387c,
  145. 0x3ff6,0x2e42,0xfefa,0x39ef
  146. };
  147. static short Q[] =
  148. {
  149. 0x3efe,0xd61f,0x30e5,0x6ff3,
  150. 0x3f4d,0xf4bc,0xbee9,0xb2c2,
  151. 0x3f78,0x5527,0x1dd8,0x8f27,
  152. 0x3f8f,0xb6a0,0x73d2,0xd3f0,
  153. 0x3f98,0x7b97,0x29ca,0x99fb,
  154. 0x3f9e,0xd7e6,0xa085,0x7bd9,
  155. 0x3fa3,0x2322,0x7b17,0x9faa,
  156. 0x3fa8,0xfffc,0xf8d1,0xc38d,
  157. 0x3fb1,0xffff,0xfeb2,0x2f00,
  158. 0x3fbf,0xffff,0xffff,0xdba4,
  159. 0x3fe0,0x0000,0x0000,0x0000
  160. };
  161. static short ac1[] = {
  162. 0x3ff6,0x2e42,0xfefa,0x39ef
  163. };
  164. #define C1 (*(double *)ac1)
  165. #endif
  166.  
  167. #ifdef UNK
  168. static double P[] =
  169. {
  170.  1.37982864606273237150E-4,
  171.  2.28025724005875567385E-3,
  172.  7.97404013220415179367E-3,
  173.  9.85821379021226008714E-3,
  174.  6.87489687449949877925E-3,
  175.  6.18901033637687613229E-3,
  176.  8.79078273952743772254E-3,
  177.  1.49380448916805252718E-2,
  178.  3.08851465246711995998E-2,
  179.  9.65735902811690126535E-2,
  180.  1.38629436111989062502E0
  181. };
  182.  
  183. static double Q[] =
  184. {
  185.  2.94078955048598507511E-5,
  186.  9.14184723865917226571E-4,
  187.  5.94058303753167793257E-3,
  188.  1.54850516649762399335E-2,
  189.  2.39089602715924892727E-2,
  190.  3.01204715227604046988E-2,
  191.  3.73774314173823228969E-2,
  192.  4.88280347570998239232E-2,
  193.  7.03124996963957469739E-2,
  194.  1.24999999999870820058E-1,
  195.  4.99999999999999999821E-1
  196. };
  197. static double C1 = 1.3862943611198906188E0; /* log(4) */
  198. #endif
  199.  
  200. extern double MACHEP, MAXNUM;
  201.  
  202. double ellpk(x)
  203. double x;
  204. {
  205. double polevl(), p1evl(), log();
  206.  
  207.  
  208. if( (x < 0.0) || (x > 1.0) )
  209.     {
  210.     mtherr( "ellpk", DOMAIN );
  211.     return( 0.0 );
  212.     }
  213.  
  214. if( x > MACHEP )
  215.     {
  216.     return( polevl(x,P,10) - log(x) * polevl(x,Q,10) );
  217.     }
  218. else
  219.     {
  220.     if( x == 0.0 )
  221.         {
  222.         mtherr( "ellpk", SING );
  223.         return( MAXNUM );
  224.         }
  225.     else
  226.         {
  227.         return( C1 - 0.5 * log(x) );
  228.         }
  229.     }
  230. }
  231.