Mercurial > projects > ldc

comparison tango/tango/math/Elliptic.d @ 132:1700239cab2e trunk

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[svn r136] MAJOR UNSTABLE UPDATE!!! Initial commit after moving to Tango instead of Phobos. Lots of bugfixes... This build is not suitable for most things.
author lindquist
date Fri, 11 Jan 2008 17:57:40 +0100
parents
children
comparison
equal deleted inserted replaced
131:5825d48b27d1 132:1700239cab2e
1 /**
2 * Elliptic integrals.
3 * The functions are named similarly to the names used in Mathematica.
4 *
5 * License: BSD style: $(LICENSE)
6 * Copyright: Based on the CEPHES math library, which is
7 * Copyright (C) 1994 Stephen L. Moshier (moshier@world.std.com).
8 * Authors: Stephen L. Moshier (original C code). Conversion to D by Don Clugston
9 *
10 * References:
11 * $(LINK http://en.wikipedia.org/wiki/Elliptic_integral)
12 *
13 * Eric W. Weisstein. "Elliptic Integral of the First Kind." From MathWorld--A Wolfram Web Resource. $(LINK http://mathworld.wolfram.com/EllipticIntegraloftheFirstKind.html)
14 *
15 * $(LINK http://www.netlib.org/cephes/ldoubdoc.html)
16 *
17 * Macros:
18 * TABLE_SV = <table border=1 cellpadding=4 cellspacing=0>
19 * <caption>Special Values</caption>
20 * $0</table>
21 * SVH = $(TR $(TH $1) $(TH $2))
22 * SV = $(TR $(TD $1) $(TD $2))
23 * GAMMA = &#915;
24 * INTEGRATE = $(BIG &#8747;<sub>$(SMALL $1)</sub><sup>$2</sup>)
25 * POWER = $1<sup>$2</sup>
26 * NAN = $(RED NAN)
27 */
28 /**
29 * Macros:
30 * TABLE_SV = <table border=1 cellpadding=4 cellspacing=0>
31 * <caption>Special Values</caption>
32 * $0</table>
33 * SVH = $(TR $(TH $1) $(TH $2))
34 * SV = $(TR $(TD $1) $(TD $2))
35 *
36 * NAN = $(RED NAN)
37 * SUP = <span style="vertical-align:super;font-size:smaller">$0</span>
38 * GAMMA = &#915;
39 * INTEGRAL = &#8747;
40 * INTEGRATE = $(BIG &#8747;<sub>$(SMALL $1)</sub><sup>$2</sup>)
41 * POWER = $1<sup>$2</sup>
42 * BIGSUM = $(BIG &Sigma; <sup>$2</sup><sub>$(SMALL $1)</sub>)
43 * CHOOSE = $(BIG &#40;) <sup>$(SMALL $1)</sup><sub>$(SMALL $2)</sub> $(BIG &#41;)
44 */
45
46 module tango.math.Elliptic;
47
48 import tango.math.Math;
49 import tango.math.IEEE;
50
51 /* These functions are based on code from:
52 Cephes Math Library, Release 2.3: October, 1995
53 Copyright 1984, 1987, 1995 by Stephen L. Moshier
54 */
55
56
57 /**
58 * Incomplete elliptic integral of the first kind
59 *
60 * Approximates the integral
61 * F(phi | m) = $(INTEGRATE 0, phi) dt/ (sqrt( 1- m $(POWER sin, 2) t))
62 *
63 * of amplitude phi and modulus m, using the arithmetic -
64 * geometric mean algorithm.
65 */
66
67 real ellipticF(real phi, real m )
68 {
69 real a, b, c, e, temp, t, K;
70 int d, mod, sign, npio2;
71
72 if( m == 0.0L )
73 return phi;
74 a = 1.0L - m;
75 if( a == 0.0L ) {
76 if ( fabs(phi) >= PI_2 ) return real.infinity;
77 return log( tan( 0.5L*(PI_2 + phi) ) );
78 }
79 npio2 = cast(int)floor( phi/PI_2 );
80 if ( npio2 & 1 )
81 npio2 += 1;
82 if ( npio2 ) {
83 K = ellipticKComplete( a );
84 phi = phi - npio2 * PI_2;
85 } else
86 K = 0.0L;
87 if( phi < 0.0L ){
88 phi = -phi;
89 sign = -1;
90 } else
91 sign = 0;
92 b = sqrt(a);
93 t = tan( phi );
94 if( fabs(t) > 10.0L ) {
95 /* Transform the amplitude */
96 e = 1.0L/(b*t);
97 /* ... but avoid multiple recursions. */
98 if( fabs(e) < 10.0L ){
99 e = atan(e);
100 if( npio2 == 0 )
101 K = ellipticKComplete( a );
102 temp = K - ellipticF( e, m );
103 goto done;
104 }
105 }
106 a = 1.0L;
107 c = sqrt(m);
108 d = 1;
109 mod = 0;
110
111 while( fabs(c/a) > real.epsilon ) {
112 temp = b/a;
113 phi = phi + atan(t*temp) + mod * PI;
114 mod = cast(int)((phi + PI_2)/PI);
115 t = t * ( 1.0L + temp )/( 1.0L - temp * t * t );
116 c = 0.5L * ( a - b );
117 temp = sqrt( a * b );
118 a = 0.5L * ( a + b );
119 b = temp;
120 d += d;
121 }
122 temp = (atan(t) + mod * PI)/(d * a);
123
124 done:
125 if ( sign < 0 )
126 temp = -temp;
127 temp += npio2 * K;
128 return temp;
129 }
130
131
132 /**
133 * Incomplete elliptic integral of the second kind
134 *
135 * Approximates the integral
136 *
137 * E(phi | m) = $(INTEGRATE 0, phi) sqrt( 1- m $(POWER sin, 2) t) dt
138 *
139 * of amplitude phi and modulus m, using the arithmetic -
140 * geometric mean algorithm.
141 */
142
143 real ellipticE(real phi, real m )
144 {
145 real a, b, c, e, temp, t, E;
146 int d, mod, npio2, sign;
147
148 if ( m == 0.0L ) return phi;
149 real lphi = phi;
150 npio2 = cast(int)floor( lphi/PI_2 );
151 if( npio2 & 1 )
152 npio2 += 1;
153 lphi = lphi - npio2 * PI_2;
154 if( lphi < 0.0L ){
155 lphi = -lphi;
156 sign = -1;
157 } else {
158 sign = 1;
159 }
160 a = 1.0L - m;
161 E = ellipticEComplete( a );
162 if( a == 0.0L ) {
163 temp = sin( lphi );
164 goto done;
165 }
166 t = tan( lphi );
167 b = sqrt(a);
168 if ( fabs(t) > 10.0L ) {
169 /* Transform the amplitude */
170 e = 1.0L/(b*t);
171 /* ... but avoid multiple recursions. */
172 if( fabs(e) < 10.0L ){
173 e = atan(e);
174 temp = E + m * sin( lphi ) * sin( e ) - ellipticE( e, m );
175 goto done;
176 }
177 }
178 c = sqrt(m);
179 a = 1.0L;
180 d = 1;
181 e = 0.0L;
182 mod = 0;
183
184 while( fabs(c/a) > real.epsilon ) {
185 temp = b/a;
186 lphi = lphi + atan(t*temp) + mod * PI;
187 mod = cast(int)((lphi + PI_2)/PI);
188 t = t * ( 1.0L + temp )/( 1.0L - temp * t * t );
189 c = 0.5L*( a - b );
190 temp = sqrt( a * b );
191 a = 0.5L*( a + b );
192 b = temp;
193 d += d;
194 e += c * sin(lphi);
195 }
196
197 temp = E / ellipticKComplete( 1.0L - m );
198 temp *= (atan(t) + mod * PI)/(d * a);
199 temp += e;
200
201 done:
202 if( sign < 0 )
203 temp = -temp;
204 temp += npio2 * E;
205 return temp;
206 }
207
208
209 /**
210 * Complete elliptic integral of the first kind
211 *
212 * Approximates the integral
213 *
214 * K(m) = $(INTEGRATE 0, &pi/2) dt/ (sqrt( 1- m $(POWER sin, 2) t))
215 *
216 * where m = 1 - x, using the approximation
217 *
218 * P(x) - log x Q(x).
219 *
220 * The argument x is used rather than m so that the logarithmic
221 * singularity at x = 1 will be shifted to the origin; this
222 * preserves maximum accuracy.
223 *
224 * x must be in the range
225 * 0 <= x <= 1
226 *
227 * This is equivalent to ellipticF(PI_2, 1-x).
228 *
229 * K(0) = &pi/2.
230 */
231
232 real ellipticKComplete(real x)
233 in {
234 // assert(x>=0.0L && x<=1.0L);
235 }
236 body{
237
238 const real [] P = [
239 0x1.62e42fefa39ef35ap+0, // 1.3862943611198906189
240 0x1.8b90bfbe8ed811fcp-4, // 0.096573590279993142323
241 0x1.fa05af797624c586p-6, // 0.030885144578720423267
242 0x1.e979cdfac7249746p-7, // 0.01493761594388688915
243 0x1.1f4cc8890cff803cp-7, // 0.0087676982094322259125
244 0x1.7befb3bb1fa978acp-8, // 0.0057973684116620276454
245 0x1.2c2566aa1d5fe6b8p-8, // 0.0045798659940508010431
246 0x1.7333514e7fe57c98p-8, // 0.0056640695097481470287
247 0x1.09292d1c8621348cp-7, // 0.0080920667906392630755
248 0x1.b89ab5fe793a6062p-8, // 0.0067230886765842542487
249 0x1.28e9c44dc5e26e66p-9, // 0.002265267575136470585
250 0x1.c2c43245d445addap-13, // 0.00021494216542320112406
251 0x1.4ee247035a03e13p-20 // 1.2475397291548388385e-06
252 ];
253
254 const real [] Q = [
255 0x1p-1, // 0.5
256 0x1.fffffffffff635eap-4, // 0.12499999999999782631
257 0x1.1fffffff8a2bea1p-4, // 0.070312499993302227507
258 0x1.8ffffe6f40ec2078p-5, // 0.04882812208418620146
259 0x1.323f4dbf7f4d0c2ap-5, // 0.037383701182969303058
260 0x1.efe8a028541b50bp-6, // 0.030267864612427881354
261 0x1.9d58c49718d6617cp-6, // 0.025228683455123323041
262 0x1.4d1a8d2292ff6e2ep-6, // 0.020331037356569904872
263 0x1.b637687027d664aap-7, // 0.013373304362459048444
264 0x1.687a640ae5c71332p-8, // 0.0055004591221382442135
265 0x1.0f9c30a94a1dcb4ep-10, // 0.001036110372590318803
266 0x1.d321746708e92d48p-15 // 5.568631677757315399e-05
267 ];
268
269 const real LOG4 = 0x1.62e42fefa39ef358p+0; // log(4)
270
271 if( x > real.epsilon )
272 return poly(x,P) - log(x) * poly(x,Q);
273 if ( x == 0.0L )
274 return real.infinity;
275 return LOG4 - 0.5L * log(x);
276 }
277
278 /**
279 * Complete elliptic integral of the second kind
280 *
281 * Approximates the integral
282 *
283 * E(m) = $(INTEGRATE 0, &pi/2) sqrt( 1- m $(POWER sin, 2) t) dt
284 *
285 * Where m = 1 - m1, using the approximation
286 *
287 * P(x) - x log x Q(x).
288 *
289 * Though there are no singularities, the argument m1 is used
290 * rather than m for compatibility with ellipticKComplete().
291 *
292 * E(1) = 1; E(0) = &pi/2.
293 * m must be in the range 0 <= m <= 1.
294 */
295
296 real ellipticEComplete(real x)
297 in {
298 assert(x>=0 && x<=1.0);
299 }
300 body {
301 const real [] P = [
302 0x1.c5c85fdf473f78f2p-2, // 0.44314718055994670505
303 0x1.d1591f9e9a66477p-5, // 0.056805192715569305834
304 0x1.65af6a7a61f587cp-6, // 0.021831373198011179718
305 0x1.7a4d48ed00d5745ap-7, // 0.011544857605264509506
306 0x1.d4f5fe4f93b60688p-8, // 0.0071557756305783152481
307 0x1.4cb71c73bac8656ap-8, // 0.0050768322432573952962
308 0x1.4a9167859a1d0312p-8, // 0.0050440671671840438539
309 0x1.dd296daa7b1f5b7ap-8, // 0.0072809117068399675418
310 0x1.04f2c29224ba99b6p-7, // 0.0079635095646944542686
311 0x1.0f5820e2d80194d8p-8, // 0.0041403847015715420009
312 0x1.95ee634752ca69b6p-11, // 0.00077425232385887751162
313 0x1.0c58aa9ab404f4fp-15 // 3.1989378120323412946e-05
314 ];
315
316 const real [] Q = [
317 0x1.ffffffffffffb1cep-3, // 0.24999999999999986434
318 0x1.7ffffffff29eaa0cp-4, // 0.093749999999239422678
319 0x1.dfffffbd51eb098p-5, // 0.058593749514839092674
320 0x1.5dffd791cb834c92p-5, // 0.04272453406734691973
321 0x1.1397b63c2f09a8ep-5, // 0.033641677787700181541
322 0x1.c567cde5931e75bcp-6, // 0.02767367465121309044
323 0x1.75e0cae852be9ddcp-6, // 0.022819708015315777007
324 0x1.12bb968236d4e434p-6, // 0.016768357258894633433
325 0x1.1f6572c1c402d07cp-7, // 0.0087706384979640787504
326 0x1.452c6909f88b8306p-9, // 0.0024808767529843311337
327 0x1.1f7504e72d664054p-12, // 0.00027414045912208516032
328 0x1.ad17054dc46913e2p-18 // 6.3939381343012054851e-06
329 ];
330 if (x==0)
331 return 1.0L;
332 return 1.0L + x * poly(x,P) - log(x) * (x * poly(x,Q) );
333 }
334
335 unittest {
336 assert( ellipticF(1, 0)==1);
337 assert(ellipticEComplete(0)==1);
338 assert(ellipticEComplete(1)==PI_2);
339 assert(feqrel(ellipticKComplete(1),PI_2)>= real.mant_dig-1);
340 assert(ellipticKComplete(0)==real.infinity);
341 // assert(ellipticKComplete(1)==0); //-real.infinity);
342
343 real x=0.5653L;
344 assert(ellipticKComplete(1-x) == ellipticF(PI_2, x) );
345 assert(ellipticEComplete(1-x) == ellipticE(PI_2, x) );
346 }