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 internally 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 2.5e-16 6.8e-17
46 *
47 * ERROR MESSAGES:
48 *
49 * message condition value returned
50 * ellpk domain x<0, x>1 0.0
51 *
52 */
53
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 double P[] = {
66 1.37982864606273237150E-4,
67 2.28025724005875567385E-3,
68 7.97404013220415179367E-3,
69 9.85821379021226008714E-3,
70 6.87489687449949877925E-3,
71 6.18901033637687613229E-3,
72 8.79078273952743772254E-3,
73 1.49380448916805252718E-2,
74 3.08851465246711995998E-2,
75 9.65735902811690126535E-2,
76 1.38629436111989062502E0
77 };
78
79 static double Q[] = {
80 2.94078955048598507511E-5,
81 9.14184723865917226571E-4,
82 5.94058303753167793257E-3,
83 1.54850516649762399335E-2,
84 2.39089602715924892727E-2,
85 3.01204715227604046988E-2,
86 3.73774314173823228969E-2,
87 4.88280347570998239232E-2,
88 7.03124996963957469739E-2,
89 1.24999999999870820058E-1,
90 4.99999999999999999821E-1
91 };
92
93 static double C1 = 1.3862943611198906188E0; /* log(4) */
94
95 extern double MACHEP;
96
ellpk(x)97 double ellpk(x)
98 double x;
99 {
100
101 if (x < 0.0) {
102 sf_error("ellpk", SF_ERROR_DOMAIN, NULL);
103 return (NPY_NAN);
104 }
105
106 if (x > 1.0) {
107 if (cephes_isinf(x)) {
108 return 0.0;
109 }
110 return ellpk(1/x)/sqrt(x);
111 }
112
113 if (x > MACHEP) {
114 return (polevl(x, P, 10) - log(x) * polevl(x, Q, 10));
115 }
116 else {
117 if (x == 0.0) {
118 sf_error("ellpk", SF_ERROR_SINGULAR, NULL);
119 return (NPY_INFINITY);
120 }
121 else {
122 return (C1 - 0.5 * log(x));
123 }
124 }
125 }
126