1 /* ellik.c
2 *
3 * Incomplete elliptic integral of the first kind
4 *
5 *
6 *
7 * SYNOPSIS:
8 *
9 * double phi, m, y, ellik();
10 *
11 * y = ellik( phi, m );
12 *
13 *
14 *
15 * DESCRIPTION:
16 *
17 * Approximates the integral
18 *
19 *
20 *
21 * phi
22 * -
23 * | |
24 * | dt
25 * F(phi | m) = | ------------------
26 * | 2
27 * | | sqrt( 1 - m sin t )
28 * -
29 * 0
30 *
31 * of amplitude phi and modulus m, using the arithmetic -
32 * geometric mean algorithm.
33 *
34 *
35 *
36 *
37 * ACCURACY:
38 *
39 * Tested at random points with m in [0, 1] and phi as indicated.
40 *
41 * Relative error:
42 * arithmetic domain # trials peak rms
43 * IEEE -10,10 200000 7.4e-16 1.0e-16
44 *
45 *
46 */
47
48
49 /*
50 * Cephes Math Library Release 2.0: April, 1987
51 * Copyright 1984, 1987 by Stephen L. Moshier
52 * Direct inquiries to 30 Frost Street, Cambridge, MA 02140
53 */
54 /* Copyright 2014, Eric W. Moore */
55
56 /* Incomplete elliptic integral of first kind */
57
58 #include "mconf.h"
59 extern double MACHEP;
60
61 #include <numpy/npy_math.h>
62
63 static double ellik_neg_m(double phi, double m);
64
ellik(double phi,double m)65 double ellik(double phi, double m)
66 {
67 double a, b, c, e, temp, t, K, denom, npio2;
68 int d, mod, sign;
69
70 if (cephes_isnan(phi) || cephes_isnan(m))
71 return NPY_NAN;
72 if (m > 1.0)
73 return NPY_NAN;
74 if (cephes_isinf(phi) || cephes_isinf(m))
75 {
76 if (cephes_isinf(m) && cephes_isfinite(phi))
77 return 0.0;
78 else if (cephes_isinf(phi) && cephes_isfinite(m))
79 return phi;
80 else
81 return NPY_NAN;
82 }
83 if (m == 0.0)
84 return (phi);
85 a = 1.0 - m;
86 if (a == 0.0) {
87 if (fabs(phi) >= (double)NPY_PI_2) {
88 sf_error("ellik", SF_ERROR_SINGULAR, NULL);
89 return (NPY_INFINITY);
90 }
91 /* DLMF 19.6.8, and 4.23.42 */
92 return npy_asinh(tan(phi));
93 }
94 npio2 = floor(phi / NPY_PI_2);
95 if (fmod(fabs(npio2), 2.0) == 1.0)
96 npio2 += 1;
97 if (npio2 != 0.0) {
98 K = ellpk(a);
99 phi = phi - npio2 * NPY_PI_2;
100 }
101 else
102 K = 0.0;
103 if (phi < 0.0) {
104 phi = -phi;
105 sign = -1;
106 }
107 else
108 sign = 0;
109 if (a > 1.0) {
110 temp = ellik_neg_m(phi, m);
111 goto done;
112 }
113 b = sqrt(a);
114 t = tan(phi);
115 if (fabs(t) > 10.0) {
116 /* Transform the amplitude */
117 e = 1.0 / (b * t);
118 /* ... but avoid multiple recursions. */
119 if (fabs(e) < 10.0) {
120 e = atan(e);
121 if (npio2 == 0)
122 K = ellpk(a);
123 temp = K - ellik(e, m);
124 goto done;
125 }
126 }
127 a = 1.0;
128 c = sqrt(m);
129 d = 1;
130 mod = 0;
131
132 while (fabs(c / a) > MACHEP) {
133 temp = b / a;
134 phi = phi + atan(t * temp) + mod * NPY_PI;
135 denom = 1.0 - temp * t * t;
136 if (fabs(denom) > 10*MACHEP) {
137 t = t * (1.0 + temp) / denom;
138 mod = (phi + NPY_PI_2) / NPY_PI;
139 }
140 else {
141 t = tan(phi);
142 mod = (int)floor((phi - atan(t))/NPY_PI);
143 }
144 c = (a - b) / 2.0;
145 temp = sqrt(a * b);
146 a = (a + b) / 2.0;
147 b = temp;
148 d += d;
149 }
150
151 temp = (atan(t) + mod * NPY_PI) / (d * a);
152
153 done:
154 if (sign < 0)
155 temp = -temp;
156 temp += npio2 * K;
157 return (temp);
158 }
159
160 /* N.B. This will evaluate its arguments multiple times. */
161 #define MAX3(a, b, c) (a > b ? (a > c ? a : c) : (b > c ? b : c))
162
163 /* To calculate legendre's incomplete elliptical integral of the first kind for
164 * negative m, we use a power series in phi for small m*phi*phi, an asymptotic
165 * series in m for large m*phi*phi* and the relation to Carlson's symmetric
166 * integral of the first kind.
167 *
168 * F(phi, m) = sin(phi) * R_F(cos(phi)^2, 1 - m * sin(phi)^2, 1.0)
169 * = R_F(c-1, c-m, c)
170 *
171 * where c = csc(phi)^2. We use the second form of this for (approximately)
172 * phi > 1/(sqrt(DBL_MAX) ~ 1e-154, where csc(phi)^2 overflows. Elsewhere we
173 * use the first form, accounting for the smallness of phi.
174 *
175 * The algorithm used is described in Carlson, B. C. Numerical computation of
176 * real or complex elliptic integrals. (1994) https://arxiv.org/abs/math/9409227
177 * Most variable names reflect Carlson's usage.
178 *
179 * In this routine, we assume m < 0 and 0 > phi > pi/2.
180 */
ellik_neg_m(double phi,double m)181 double ellik_neg_m(double phi, double m)
182 {
183 double x, y, z, x1, y1, z1, A0, A, Q, X, Y, Z, E2, E3, scale;
184 int n = 0;
185 double mpp = (m*phi)*phi;
186
187 if (-mpp < 1e-6 && phi < -m) {
188 return phi + (-mpp*phi*phi/30.0 + 3.0*mpp*mpp/40.0 + mpp/6.0)*phi;
189 }
190
191 if (-mpp > 4e7) {
192 double sm = sqrt(-m);
193 double sp = sin(phi);
194 double cp = cos(phi);
195
196 double a = log(4*sp*sm/(1+cp));
197 double b = -(1 + cp/sp/sp - a) / 4 / m;
198 return (a + b) / sm;
199 }
200
201 if (phi > 1e-153 && m > -1e305) {
202 double s = sin(phi);
203 double csc2 = 1.0 / (s*s);
204 scale = 1.0;
205 x = 1.0 / (tan(phi) * tan(phi));
206 y = csc2 - m;
207 z = csc2;
208 }
209 else {
210 scale = phi;
211 x = 1.0;
212 y = 1 - m*scale*scale;
213 z = 1.0;
214 }
215
216 if (x == y && x == z) {
217 return scale / sqrt(x);
218 }
219
220 A0 = (x + y + z) / 3.0;
221 A = A0;
222 x1 = x; y1 = y; z1 = z;
223 /* Carlson gives 1/pow(3*r, 1.0/6.0) for this constant. if r == eps,
224 * it is ~338.38. */
225 Q = 400.0 * MAX3(fabs(A0-x), fabs(A0-y), fabs(A0-z));
226
227 while (Q > fabs(A) && n <= 100) {
228 double sx = sqrt(x1);
229 double sy = sqrt(y1);
230 double sz = sqrt(z1);
231 double lam = sx*sy + sx*sz + sy*sz;
232 x1 = (x1 + lam) / 4.0;
233 y1 = (y1 + lam) / 4.0;
234 z1 = (z1 + lam) / 4.0;
235 A = (x1 + y1 + z1) / 3.0;
236 n += 1;
237 Q /= 4;
238 }
239 X = (A0 - x) / A / (1 << 2*n);
240 Y = (A0 - y) / A / (1 << 2*n);
241 Z = -(X + Y);
242
243 E2 = X*Y - Z*Z;
244 E3 = X*Y*Z;
245
246 return scale * (1.0 - E2/10.0 + E3/14.0 + E2*E2/24.0
247 - 3.0*E2*E3/44.0) / sqrt(A);
248 }
249