ellpe.c 4.7 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141
  1. /* ellpe.c
  2. *
  3. * Complete elliptic integral of the second kind
  4. *
  5. *
  6. *
  7. * SYNOPSIS:
  8. *
  9. * double m1, y, ellpe();
  10. *
  11. * y = ellpe( 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. * DEC 0, 1 13000 3.1e-17 9.4e-18
  43. * IEEE 0, 1 10000 2.1e-16 7.3e-17
  44. *
  45. *
  46. * ERROR MESSAGES:
  47. *
  48. * message condition value returned
  49. * ellpe domain x<0, x>1 0.0
  50. *
  51. */
  52. /* ellpe.c */
  53. /* Elliptic integral of second kind */
  54. /*
  55. Cephes Math Library, Release 2.8: June, 2000
  56. Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
  57. */
  58. #include "mconf.h"
  59. #ifdef UNK
  60. static double P[] = {1.53552577301013293365E-4, 2.50888492163602060990E-3,
  61. 8.68786816565889628429E-3, 1.07350949056076193403E-2,
  62. 7.77395492516787092951E-3, 7.58395289413514708519E-3,
  63. 1.15688436810574127319E-2, 2.18317996015557253103E-2,
  64. 5.68051945617860553470E-2, 4.43147180560990850618E-1,
  65. 1.00000000000000000299E0};
  66. static double Q[] = {3.27954898576485872656E-5, 1.00962792679356715133E-3,
  67. 6.50609489976927491433E-3, 1.68862163993311317300E-2,
  68. 2.61769742454493659583E-2, 3.34833904888224918614E-2,
  69. 4.27180926518931511717E-2, 5.85936634471101055642E-2,
  70. 9.37499997197644278445E-2, 2.49999999999888314361E-1};
  71. #endif
  72. #ifdef DEC
  73. static unsigned short P[] = {
  74. 0035041, 0001364, 0141572, 0117555, 0036044, 0066032, 0130027, 0033404,
  75. 0036416, 0053617, 0064456, 0102632, 0036457, 0161100, 0061177, 0122612,
  76. 0036376, 0136251, 0012403, 0124162, 0036370, 0101316, 0151715, 0131613,
  77. 0036475, 0105477, 0050317, 0133272, 0036662, 0154232, 0024645, 0171552,
  78. 0037150, 0126220, 0047054, 0030064, 0037742, 0162057, 0167645, 0165612,
  79. 0040200, 0000000, 0000000, 0000000};
  80. static unsigned short Q[] = {
  81. 0034411, 0106743, 0115771, 0055462, 0035604, 0052575, 0155171, 0045540,
  82. 0036325, 0030424, 0064332, 0167756, 0036612, 0052366, 0063006, 0115175,
  83. 0036726, 0070430, 0004533, 0124654, 0037011, 0022741, 0030675, 0030711,
  84. 0037056, 0174452, 0127062, 0132122, 0037157, 0177750, 0142041, 0072523,
  85. 0037277, 0177777, 0173137, 0002627, 0037577, 0177777, 0177777, 0101101};
  86. #endif
  87. #ifdef IBMPC
  88. static unsigned short P[] = {
  89. 0x53ee, 0x986f, 0x205e, 0x3f24, 0xe6e0, 0x5602, 0x8d83, 0x3f64, 0xd0b3,
  90. 0xed25, 0xcaf1, 0x3f81, 0xf4b1, 0x0c4f, 0xfc48, 0x3f85, 0x750e, 0x22a0,
  91. 0xd795, 0x3f7f, 0xb671, 0xda79, 0x1059, 0x3f7f, 0xf6d7, 0xea19, 0xb167,
  92. 0x3f87, 0xbe6d, 0x4534, 0x5b13, 0x3f96, 0x8607, 0x09c5, 0x1592, 0x3fad,
  93. 0xbd71, 0xfdf4, 0x5c85, 0x3fdc, 0x0000, 0x0000, 0x0000, 0x3ff0};
  94. static unsigned short Q[] = {
  95. 0x2b66, 0x737f, 0x31bc, 0x3f01, 0x296c, 0xbb4f, 0x8aaf, 0x3f50,
  96. 0x5dfe, 0x8d1b, 0xa622, 0x3f7a, 0xd350, 0xccc0, 0x4a9e, 0x3f91,
  97. 0x7535, 0x012b, 0xce23, 0x3f9a, 0xa639, 0x2637, 0x24bc, 0x3fa1,
  98. 0x568a, 0x55c6, 0xdf25, 0x3fa5, 0x2eaa, 0x1884, 0xfffd, 0x3fad,
  99. 0xe0b3, 0xfecb, 0xffff, 0x3fb7, 0xf048, 0xffff, 0xffff, 0x3fcf};
  100. #endif
  101. #ifdef MIEEE
  102. static unsigned short P[] = {
  103. 0x3f24, 0x205e, 0x986f, 0x53ee, 0x3f64, 0x8d83, 0x5602, 0xe6e0, 0x3f81,
  104. 0xcaf1, 0xed25, 0xd0b3, 0x3f85, 0xfc48, 0x0c4f, 0xf4b1, 0x3f7f, 0xd795,
  105. 0x22a0, 0x750e, 0x3f7f, 0x1059, 0xda79, 0xb671, 0x3f87, 0xb167, 0xea19,
  106. 0xf6d7, 0x3f96, 0x5b13, 0x4534, 0xbe6d, 0x3fad, 0x1592, 0x09c5, 0x8607,
  107. 0x3fdc, 0x5c85, 0xfdf4, 0xbd71, 0x3ff0, 0x0000, 0x0000, 0x0000};
  108. static unsigned short Q[] = {
  109. 0x3f01, 0x31bc, 0x737f, 0x2b66, 0x3f50, 0x8aaf, 0xbb4f, 0x296c,
  110. 0x3f7a, 0xa622, 0x8d1b, 0x5dfe, 0x3f91, 0x4a9e, 0xccc0, 0xd350,
  111. 0x3f9a, 0xce23, 0x012b, 0x7535, 0x3fa1, 0x24bc, 0x2637, 0xa639,
  112. 0x3fa5, 0xdf25, 0x55c6, 0x568a, 0x3fad, 0xfffd, 0x1884, 0x2eaa,
  113. 0x3fb7, 0xffff, 0xfecb, 0xe0b3, 0x3fcf, 0xffff, 0xffff, 0xf048};
  114. #endif
  115. #ifdef ANSIPROT
  116. extern double polevl(double, void *, int);
  117. extern double log(double);
  118. #else
  119. double polevl(), log();
  120. #endif
  121. double ellpe(x) double x;
  122. {
  123. if ((x <= 0.0) || (x > 1.0)) {
  124. if (x == 0.0)
  125. return (1.0);
  126. mtherr("ellpe", DOMAIN);
  127. return (0.0);
  128. }
  129. return (polevl(x, P, 10) - log(x) * (x * polevl(x, Q, 9)));
  130. }