1 /**
2  * Elliptic integrals.
3  * The functions are named similarly to the names used in Mathematica.
4  *
5  * Eric W. Weisstein. "Elliptic Integral of the First Kind." From MathWorld--A Wolfram Web Resource. $(LINK http://mathworld.wolfram.com/EllipticIntegraloftheFirstKind.html)
6  *
7  * $(LINK http://www.netlib.org/cephes/ldoubdoc.html)
8  *
9  * Copyright:
10  *     Based on the CEPHES math library, which is
11  *                Copyright (C) 1994 Stephen L. Moshier (moshier@world.std.com).
12  *     Some parts copyright (c) 2009-2016 dunnhumby Germany GmbH.
13  *     All rights reserved.
14  *
15  * License:
16  *     Tango Dual License: 3-Clause BSD License / Academic Free License v3.0.
17  *     See LICENSE_TANGO.txt for details.
18  *
19  * Authors: Stephen L. Moshier (original C code). Conversion to D by Don Clugston
20  *
21  * References: $(LINK http://en.wikipedia.org/wiki/Elliptic_integral)
22  *
23  * Macros:
24  *  TABLE_SV = <table border=1 cellpadding=4 cellspacing=0>
25  *      <caption>Special Values</caption>
26  *      $0</table>
27  *  SVH = $(TR $(TH $1) $(TH $2))
28  *  SV  = $(TR $(TD $1) $(TD $2))
29  *  GAMMA =  &#915;
30  *  INTEGRATE = $(BIG &#8747;<sub>$(SMALL $1)</sub><sup>$2</sup>)
31  *  POWER = $1<sup>$2</sup>
32  *  NAN = $(RED NAN)
33  */
34 /**
35  * Macros:
36  *  TABLE_SV = <table border=1 cellpadding=4 cellspacing=0>
37  *      <caption>Special Values</caption>
38  *      $0</table>
39  *  SVH = $(TR $(TH $1) $(TH $2))
40  *  SV  = $(TR $(TD $1) $(TD $2))
41  *
42  *  NAN = $(RED NAN)
43  *  SUP = <span style="vertical-align:super;font-size:smaller">$0</span>
44  *  GAMMA =  &#915;
45  *  INTEGRAL = &#8747;
46  *  INTEGRATE = $(BIG &#8747;<sub>$(SMALL $1)</sub><sup>$2</sup>)
47  *  POWER = $1<sup>$2</sup>
48  *  BIGSUM = $(BIG &Sigma; <sup>$2</sup><sub>$(SMALL $1)</sub>)
49  *  CHOOSE = $(BIG &#40;) <sup>$(SMALL $1)</sup><sub>$(SMALL $2)</sub> $(BIG &#41;)
50  */
51 
52 module ocean.math.Elliptic;
53 
54 import ocean.math.Math;
55 import ocean.math.IEEE;
56 import ocean.core.Verify;
57 
58 version (unittest) import ocean.core.Test;
59 
60 /* These functions are based on code from:
61 Cephes Math Library, Release 2.3:  October, 1995
62 Copyright 1984, 1987, 1995 by Stephen L. Moshier
63 */
64 
65 
66 /**
67  *  Incomplete elliptic integral of the first kind
68  *
69  * Approximates the integral
70  *   F(phi | m) = $(INTEGRATE 0, phi) dt/ (sqrt( 1- m $(POWER sin, 2) t))
71  *
72  * of amplitude phi and modulus m, using the arithmetic -
73  * geometric mean algorithm.
74  */
75 
76 real ellipticF(real phi, real m )
77 {
78     real a, b, c, e, temp, t, K;
79     int d, mod, sign, npio2;
80 
81     if( m == 0.0L )
82         return phi;
83     a = 1.0L - m;
84     if( a == 0.0L ) {
85         if ( fabs(phi) >= PI_2 )  return real.infinity;
86         return  log(  tan( 0.5L*(PI_2 + phi) )  );
87     }
88     npio2 = cast(int)floor( phi/PI_2 );
89     if ( npio2 & 1 )
90         npio2 += 1;
91     if ( npio2 ) {
92         K = ellipticKComplete( a );
93         phi = phi - npio2 * PI_2;
94     } else
95         K = 0.0L;
96     if( phi < 0.0L ){
97         phi = -phi;
98         sign = -1;
99     } else sign = 0;
100     b = sqrt(a);
101     t = tan( phi );
102     if( fabs(t) > 10.0L ) {
103     /* Transform the amplitude */
104     e = 1.0L/(b*t);
105     /* ... but avoid multiple recursions.  */
106     if( fabs(e) < 10.0L ){
107         e = atan(e);
108         if( npio2 == 0 )
109             K = ellipticKComplete( a );
110             temp = K - ellipticF( e, m );
111             goto done;
112         }
113     }
114     a = 1.0L;
115     c = sqrt(m);
116     d = 1;
117     mod = 0;
118 
119     while( fabs(c/a) > real.epsilon ) {
120         temp = b/a;
121         phi = phi + atan(t*temp) + mod * PI;
122         mod = cast(int)((phi + PI_2)/PI);
123         t = t * ( 1.0L + temp )/( 1.0L - temp * t * t );
124         c = 0.5L * ( a - b );
125         temp = sqrt( a * b );
126         a = 0.5L * ( a + b );
127         b = temp;
128         d += d;
129     }
130     temp = (atan(t) + mod * PI)/(d * a);
131 
132 done:
133     if ( sign < 0 )
134         temp = -temp;
135     temp += npio2 * K;
136     return temp;
137 }
138 
139 
140 /**
141  *  Incomplete elliptic integral of the second kind
142  *
143  * Approximates the integral
144  *
145  * E(phi | m) = $(INTEGRATE 0, phi) sqrt( 1- m $(POWER sin, 2) t) dt
146  *
147  * of amplitude phi and modulus m, using the arithmetic -
148  * geometric mean algorithm.
149  */
150 
151 real ellipticE(real phi, real m)
152 {
153     real a, b, c, e, temp, t, E;
154     int d, mod, npio2, sign;
155 
156     if ( m == 0.0L ) return phi;
157     real lphi = phi;
158     npio2 = cast(int)floor( lphi/PI_2 );
159     if( npio2 & 1 )
160         npio2 += 1;
161     lphi = lphi - npio2 * PI_2;
162     if( lphi < 0.0L ){
163         lphi = -lphi;
164         sign = -1;
165     } else  {
166         sign = 1;
167     }
168     a = 1.0L - m;
169     E = ellipticEComplete( a );
170     if( a == 0.0L ) {
171         temp = sin( lphi );
172         goto done;
173     }
174     t = tan( lphi );
175     b = sqrt(a);
176     if ( fabs(t) > 10.0L ) {
177         /* Transform the amplitude */
178         e = 1.0L/(b*t);
179         /* ... but avoid multiple recursions.  */
180         if( fabs(e) < 10.0L ){
181             e = atan(e);
182             temp = E + m * sin( lphi ) * sin( e ) - ellipticE( e, m );
183             goto done;
184         }
185     }
186     c = sqrt(m);
187     a = 1.0L;
188     d = 1;
189     e = 0.0L;
190     mod = 0;
191 
192     while( fabs(c/a) > real.epsilon ) {
193         temp = b/a;
194         lphi = lphi + atan(t*temp) + mod * PI;
195         mod = cast(int)((lphi + PI_2)/PI);
196         t = t * ( 1.0L + temp )/( 1.0L - temp * t * t );
197         c = 0.5L*( a - b );
198         temp = sqrt( a * b );
199         a = 0.5L*( a + b );
200         b = temp;
201         d += d;
202         e += c * sin(lphi);
203     }
204 
205     temp = E / ellipticKComplete( 1.0L - m );
206     temp *= (atan(t) + mod * PI)/(d * a);
207     temp += e;
208 
209 done:
210     if( sign < 0 )
211         temp = -temp;
212     temp += npio2 * E;
213     return temp;
214 }
215 
216 
217 /**
218  *  Complete elliptic integral of the first kind
219  *
220  * Approximates the integral
221  *
222  *   K(m) = $(INTEGRATE 0, &pi;/2) dt/ (sqrt( 1- m $(POWER sin, 2) t))
223  *
224  * where m = 1 - x, using the approximation
225  *
226  *     P(x)  -  log x Q(x).
227  *
228  * The argument x is used rather than m so that the logarithmic
229  * singularity at x = 1 will be shifted to the origin; this
230  * preserves maximum accuracy.
231  *
232  * x must be in the range
233  *  0 <= x <= 1
234  *
235  * This is equivalent to ellipticF(PI_2, 1-x).
236  *
237  * K(0) = &pi;/2.
238  */
239 
240 real ellipticKComplete(real x)
241 {
242 //    verify(x>=0.0L && x<=1.0L);
243 
244 static immutable real [] P = [
245    0x1.62e42fefa39ef35ap+0, // 1.3862943611198906189
246    0x1.8b90bfbe8ed811fcp-4, // 0.096573590279993142323
247    0x1.fa05af797624c586p-6, // 0.030885144578720423267
248    0x1.e979cdfac7249746p-7, // 0.01493761594388688915
249    0x1.1f4cc8890cff803cp-7, // 0.0087676982094322259125
250    0x1.7befb3bb1fa978acp-8, // 0.0057973684116620276454
251    0x1.2c2566aa1d5fe6b8p-8, // 0.0045798659940508010431
252    0x1.7333514e7fe57c98p-8, // 0.0056640695097481470287
253    0x1.09292d1c8621348cp-7, // 0.0080920667906392630755
254    0x1.b89ab5fe793a6062p-8, // 0.0067230886765842542487
255    0x1.28e9c44dc5e26e66p-9, // 0.002265267575136470585
256    0x1.c2c43245d445addap-13,    // 0.00021494216542320112406
257    0x1.4ee247035a03e13p-20  // 1.2475397291548388385e-06
258 ];
259 
260 static immutable real [] Q = [
261    0x1p-1,  // 0.5
262    0x1.fffffffffff635eap-4, // 0.12499999999999782631
263    0x1.1fffffff8a2bea1p-4,  // 0.070312499993302227507
264    0x1.8ffffe6f40ec2078p-5, // 0.04882812208418620146
265    0x1.323f4dbf7f4d0c2ap-5, // 0.037383701182969303058
266    0x1.efe8a028541b50bp-6,  // 0.030267864612427881354
267    0x1.9d58c49718d6617cp-6, // 0.025228683455123323041
268    0x1.4d1a8d2292ff6e2ep-6, // 0.020331037356569904872
269    0x1.b637687027d664aap-7, // 0.013373304362459048444
270    0x1.687a640ae5c71332p-8, // 0.0055004591221382442135
271    0x1.0f9c30a94a1dcb4ep-10,    // 0.001036110372590318803
272    0x1.d321746708e92d48p-15     // 5.568631677757315399e-05
273 ];
274 
275     static immutable real LOG4 = 0x1.62e42fefa39ef358p+0;  // log(4)
276 
277     if( x > real.epsilon )
278         return poly(x,P) - log(x) * poly(x,Q);
279     if ( x == 0.0L )
280         return real.infinity;
281     return LOG4 - 0.5L * log(x);
282 }
283 
284 /**
285  *  Complete elliptic integral of the second kind
286  *
287  * Approximates the integral
288  *
289  * E(m) = $(INTEGRATE 0, &pi;/2) sqrt( 1- m $(POWER sin, 2) t) dt
290  *
291  * where m = 1 - x, using the approximation
292  *
293  *      P(x)  -  x log x Q(x).
294  *
295  * Though there are no singularities, the argument m1 is used
296  * rather than m for compatibility with ellipticKComplete().
297  *
298  * E(1) = 1; E(0) = &pi;/2.
299  * m must be in the range 0 <= m <= 1.
300  */
301 
302 real ellipticEComplete(real x)
303 {
304 verify(x>=0 && x<=1.0);
305 static immutable real [] P = [
306    0x1.c5c85fdf473f78f2p-2, // 0.44314718055994670505
307    0x1.d1591f9e9a66477p-5,  // 0.056805192715569305834
308    0x1.65af6a7a61f587cp-6,  // 0.021831373198011179718
309    0x1.7a4d48ed00d5745ap-7, // 0.011544857605264509506
310    0x1.d4f5fe4f93b60688p-8, // 0.0071557756305783152481
311    0x1.4cb71c73bac8656ap-8, // 0.0050768322432573952962
312    0x1.4a9167859a1d0312p-8, // 0.0050440671671840438539
313    0x1.dd296daa7b1f5b7ap-8, // 0.0072809117068399675418
314    0x1.04f2c29224ba99b6p-7, // 0.0079635095646944542686
315    0x1.0f5820e2d80194d8p-8, // 0.0041403847015715420009
316    0x1.95ee634752ca69b6p-11,    // 0.00077425232385887751162
317    0x1.0c58aa9ab404f4fp-15  // 3.1989378120323412946e-05
318 ];
319 
320 static immutable real [] Q = [
321    0x1.ffffffffffffb1cep-3, // 0.24999999999999986434
322    0x1.7ffffffff29eaa0cp-4, // 0.093749999999239422678
323    0x1.dfffffbd51eb098p-5,  // 0.058593749514839092674
324    0x1.5dffd791cb834c92p-5, // 0.04272453406734691973
325    0x1.1397b63c2f09a8ep-5,  // 0.033641677787700181541
326    0x1.c567cde5931e75bcp-6, // 0.02767367465121309044
327    0x1.75e0cae852be9ddcp-6, // 0.022819708015315777007
328    0x1.12bb968236d4e434p-6, // 0.016768357258894633433
329    0x1.1f6572c1c402d07cp-7, // 0.0087706384979640787504
330    0x1.452c6909f88b8306p-9, // 0.0024808767529843311337
331    0x1.1f7504e72d664054p-12,    // 0.00027414045912208516032
332    0x1.ad17054dc46913e2p-18     // 6.3939381343012054851e-06
333 ];
334     if (x==0)
335         return 1.0L;
336     return 1.0L + x * poly(x,P) - log(x) * (x * poly(x,Q) );
337 }
338 
339 unittest {
340     test( ellipticF(1, 0)==1);
341     test(ellipticEComplete(0)==1);
342     test(ellipticEComplete(1)==PI_2);
343     test(feqrel(ellipticKComplete(1),PI_2)>= real.mant_dig-1);
344     test(ellipticKComplete(0)==real.infinity);
345 //    assert(ellipticKComplete(1)==0); //-real.infinity);
346 
347     real x=0.5653L;
348     test(ellipticKComplete(1-x) == ellipticF(PI_2, x) );
349     test(ellipticEComplete(1-x) == ellipticE(PI_2, x) );
350 }
351 
352 /**
353  *  Incomplete elliptic integral of the third kind
354  *
355  * Approximates the integral
356  *
357  * PI(n; phi | m) = $(INTEGRATE t=0, phi) dt/((1 - n $(POWER sin,2)t) * sqrt( 1- m $(POWER sin, 2) t))
358  *
359  * of amplitude phi, modulus m, and characteristic n using Gauss-Legendre
360  * quadrature.
361  *
362  * Note that ellipticPi(PI_2, m, 1) is infinite for any m.
363  */
364 real ellipticPi(real phi, real m, real n)
365 {
366     // BUGS: This implementation suffers from poor precision.
367     static immutable double [] t = [
368         0.9931285991850949, 0.9639719272779138,
369         0.9122344282513259, 0.8391169718222188,
370         0.7463319064601508, 0.6360536807265150,
371         0.5108670019508271, 0.3737060887154195,
372         0.2277858511416451, 0.7652652113349734e-1
373     ];
374     static immutable double [] w =[
375         0.1761400713915212e-1, 0.4060142980038694e-1,
376         0.6267204833410907e-1, 0.8327674157670475e-1,
377         0.1019301198172404, 0.1181945319615184,
378         0.1316886384491766, 0.1420961093183820,
379         0.1491729864726037, 0.1527533871307258
380     ];
381         bool b1 = (m==1) && abs(phi-90)<=1e-8;
382         bool b2 = (n==1) && abs(phi-90)<=1e-8;
383         if (b1 || b2) return real.infinity;
384         real c1 = 0.87266462599716e-2 * phi;
385         real c2 = c1;
386         double x = 0;
387         for (int i=0; i< t.length; ++i) {
388             real c0 = c2 * t[i];
389             real t1 = c1 + c0;
390             real t2 = c1 - c0;
391             real s1 = sin(t1);  //  sin(c1 * (1 + t[i]))
392             real s2 = sin(t2);  //  sin(c1 * (1 - t[i]))
393             real f1 = 1.0 / ((1.0 - n * s1 * s1) * sqrt(1.0 - m * s1 * s1));
394             real f2 = 1.0 / ((1.0 - n * s2 * s2) * sqrt(1.0 - m * s2 * s2));
395             x+= w[i]*(f1+f2);
396         }
397         return c1 * x;
398 }
399 
400 /**
401  *  Complete elliptic integral of the third kind
402  */
403 real ellipticPiComplete(real m, real n)
404 {
405     verify(m>=-1.0 && m<=1.0);
406     return ellipticPi(PI_2, m, n);
407 }