scid.functions source code

1 /** Special functions.
2 
3     The functions in this module are ported (with permission) to D from
4     $(LINK2 http://www.moshier.net/,Stephen L. Moshier)'s
5     $(LINK2 http://www.netlib.org/cephes/,Cephes) library.
6 
7     Authors:    Lars Tandle Kyllingstad
8     Copyright:  Copyright (c) 2010–2011, Lars T. Kyllingstad. All rights reserved.
9     License:    Boost License 1.0
10 */
11 module scid.functions;
12 
13 import std.math;
14 
15 
16 
17 version (LittleEndian)
18 {
19     static if (real.sizeof == 12)
20     {
21         version = LEReal80;
22         private enum short PAD = 0;
23     }
24 }
25 
26 
27 
28 
29 /*  Evaluate Chebyshev series.
30 
31     Coefficients are stored in reverse order, i.e. the
32     zeroth-order term is the last in the array.
33 */
34 private real chebEval(real x, const real[] coefficients) @safe pure nothrow
35 {
36     real
37         b0 = coefficients[0],
38         b1 = 0,
39         b2;
40 
41     foreach (c; coefficients[1 .. $])
42     {
43         b2 = b1;
44         b1 = b0;
45         b0 = x * b1 - b2 + c;
46     }
47 
48     return 0.5 * (b0 - b2);
49 }
50 
51 
52 
53 
54 /** Modified Bessel function of the first kind, order 0. */
55 real besselI0(real x) pure
56 {
57     // Chebyshev coefficients for exp(-x) I0(x)
58     // in the interval [0,8].
59     //
60     // lim(x->0){ exp(-x) I0(x) } = 1.
61     immutable real[30] chebCoeffsA = [
62        -4.41534164647933937950e-18L,
63         3.33079451882223809783e-17L,
64        -2.43127984654795469359e-16L,
65         1.71539128555513303061e-15L,
66        -1.16853328779934516808e-14L,
67         7.67618549860493561688e-14L,
68        -4.85644678311192946090e-13L,
69         2.95505266312963983461e-12L,
70        -1.72682629144155570723e-11L,
71         9.67580903537323691224e-11L,
72        -5.18979560163526290666e-10L,
73         2.65982372468238665035e-9L,
74        -1.30002500998624804212e-8L,
75         6.04699502254191894932e-8L,
76        -2.67079385394061173391e-7L,
77         1.11738753912010371815e-6L,
78        -4.41673835845875056359e-6L,
79         1.64484480707288970893e-5L,
80        -5.75419501008210370398e-5L,
81         1.88502885095841655729e-4L,
82        -5.76375574538582365885e-4L,
83         1.63947561694133579842e-3L,
84        -4.32430999505057594430e-3L,
85         1.05464603945949983183e-2L,
86        -2.37374148058994688156e-2L,
87         4.93052842396707084878e-2L,
88        -9.49010970480476444210e-2L,
89         1.71620901522208775349e-1L,
90        -3.04682672343198398683e-1L,
91         6.76795274409476084995e-1L  ];
92 
93     // Chebyshev coefficients for exp(-x) sqrt(x) I0(x)
94     // in the inverted interval [8,infinity].
95     //
96     // lim(x->inf){ exp(-x) sqrt(x) I0(x) } = 1/sqrt(2pi).
97     immutable real[25] chebCoeffsB = [
98        -7.23318048787475395456e-18L,
99        -4.83050448594418207126e-18L,
100         4.46562142029675999901e-17L,
101         3.46122286769746109310e-17L,
102        -2.82762398051658348494e-16L,
103        -3.42548561967721913462e-16L,
104         1.77256013305652638360e-15L,
105         3.81168066935262242075e-15L,
106        -9.55484669882830764870e-15L,
107        -4.15056934728722208663e-14L,
108         1.54008621752140982691e-14L,
109         3.85277838274214270114e-13L,
110         7.18012445138366623367e-13L,
111        -1.79417853150680611778e-12L,
112        -1.32158118404477131188e-11L,
113        -3.14991652796324136454e-11L,
114         1.18891471078464383424e-11L,
115         4.94060238822496958910e-10L,
116         3.39623202570838634515e-9L,
117         2.26666899049817806459e-8L,
118         2.04891858946906374183e-7L,
119         2.89137052083475648297e-6L,
120         6.88975834691682398426e-5L,
121         3.36911647825569408990e-3L,
122         8.04490411014108831608e-1L  ];
123 
124     auto y = fabs(x);
125 
126     if (y < 8)
127         return exp(y) * chebEval(y/2 - 2, chebCoeffsA);
128     else
129         return exp(y) * chebEval(32/y - 2, chebCoeffsB) / sqrt(y);
130 }
131 
132 
133 unittest
134 {
135     assert (approxEqual(besselI0( 1), 1.2660658777520083356e0, 8.2e-17));
136     assert (approxEqual(besselI0(10), 2.8157166284662544715e3, 8.2e-17));
137     assert (besselI0(- 1) == besselI0( 1));
138     assert (besselI0(-10) == besselI0(10));
139 }
140 
141 
142 
143 
144 /** Modified Bessel function of the first kind, order 1. */
145 real besselI1(real x) pure
146 {
147     /* Chebyshev coefficients for exp(-x) I1(x) / x
148      * in the interval [0,8].
149      *
150      * lim(x->0){ exp(-x) I1(x) / x } = 1/2.
151      */
152     immutable real[29] chebCoeffsA = [
153         2.77791411276104639959e-18L,
154        -2.11142121435816608115e-17L,
155         1.55363195773620046921e-16L,
156        -1.10559694773538630805e-15L,
157         7.60068429473540693410e-15L,
158        -5.04218550472791168711e-14L,
159         3.22379336594557470981e-13L,
160        -1.98397439776494371520e-12L,
161         1.17361862988909016308e-11L,
162        -6.66348972350202774223e-11L,
163         3.62559028155211703701e-10L,
164        -1.88724975172282928790e-9L,
165         9.38153738649577178388e-9L,
166        -4.44505912879632808065e-8L,
167         2.00329475355213526229e-7L,
168        -8.56872026469545474066e-7L,
169         3.47025130813767847674e-6L,
170        -1.32731636560394358279e-5L,
171         4.78156510755005422638e-5L,
172        -1.61760815825896745588e-4L,
173         5.12285956168575772895e-4L,
174        -1.51357245063125314899e-3L,
175         4.15642294431288815669e-3L,
176        -1.05640848946261981558e-2L,
177         2.47264490306265168283e-2L,
178        -5.29459812080949914269e-2L,
179         1.02643658689847095384e-1L,
180        -1.76416518357834055153e-1L,
181         2.52587186443633654823e-1L  ];
182 
183     /* Chebyshev coefficients for exp(-x) sqrt(x) I1(x)
184      * in the inverted interval [8,infinity].
185      *
186      * lim(x->inf){ exp(-x) sqrt(x) I1(x) } = 1/sqrt(2pi).
187      */
188     immutable real[25] chebCoeffsB = [
189         7.51729631084210481353e-18L,
190         4.41434832307170791151e-18L,
191        -4.65030536848935832153e-17L,
192        -3.20952592199342395980e-17L,
193         2.96262899764595013876e-16L,
194         3.30820231092092828324e-16L,
195        -1.88035477551078244854e-15L,
196        -3.81440307243700780478e-15L,
197         1.04202769841288027642e-14L,
198         4.27244001671195135429e-14L,
199        -2.10154184277266431302e-14L,
200        -4.08355111109219731823e-13L,
201        -7.19855177624590851209e-13L,
202         2.03562854414708950722e-12L,
203         1.41258074366137813316e-11L,
204         3.25260358301548823856e-11L,
205        -1.89749581235054123450e-11L,
206        -5.58974346219658380687e-10L,
207        -3.83538038596423702205e-9L,
208        -2.63146884688951950684e-8L,
209        -2.51223623787020892529e-7L,
210        -3.88256480887769039346e-6L,
211        -1.10588938762623716291e-4L,
212        -9.76109749136146840777e-3L,
213         7.78576235018280120474e-1L  ];
214 
215     auto z = fabs(x);
216     if (z <= 8)
217     {
218         immutable y = z/2 - 2;
219         z = chebEval(y, chebCoeffsA[]) * z * exp(z);
220     }
221     else
222     {
223         z = exp(z) * chebEval(32/z - 2, chebCoeffsB[]) / sqrt(z);
224     }
225     if (x < 0) z = -z;
226 
227     return z;
228 }
229 
230 
231 unittest
232 {
233     assert (approxEqual(besselI1( 1), 5.6515910399248502721e-1, 1.2e-16));
234     assert (approxEqual(besselI1(10), 2.6709883037012546543e3,  1.2e-16));
235     assert (besselI1(- 1) == -besselI1( 1));
236     assert (besselI1(-10) == -besselI1(10));
237 }
238 
239 
240 
241 
242 /** Modified Bessel function of the second kind, order 0. */
243 real besselK0(real x) pure
244 in { assert(x > 0, "Can only calculate K0 for positive argument"); }
245 body {
246     /* Chebyshev coefficients for K0(x) + log(x/2) I0(x)
247      * in the interval [0,2].  The odd order coefficients are all
248      * zero; only the even order coefficients are listed.
249      * 
250      * lim(x->0){ K0(x) + log(x/2) I0(x) } = -EUL.
251      */
252     immutable real[10] chebCoeffsA = [
253         1.37446543561352307156e-16L,
254         4.25981614279661018399e-14L,
255         1.03496952576338420167e-11L,
256         1.90451637722020886025e-9L,
257         2.53479107902614945675e-7L,
258         2.28621210311945178607e-5L,
259         1.26461541144692592338e-3L,
260         3.59799365153615016266e-2L,
261         3.44289899924628486886e-1L,
262        -5.35327393233902768720e-1L  ];
263 
264     /* Chebyshev coefficients for exp(x) sqrt(x) K0(x)
265      * in the inverted interval [2,infinity].
266      * 
267      * lim(x->inf){ exp(x) sqrt(x) K0(x) } = sqrt(pi/2).
268      */
269     immutable real[25] chebCoeffsB = [
270         5.30043377268626276149e-18L,
271        -1.64758043015242134646e-17L,
272         5.21039150503902756861e-17L,
273        -1.67823109680541210385e-16L,
274         5.51205597852431940784e-16L,
275        -1.84859337734377901440e-15L,
276         6.34007647740507060557e-15L,
277        -2.22751332699166985548e-14L,
278         8.03289077536357521100e-14L,
279        -2.98009692317273043925e-13L,
280         1.14034058820847496303e-12L,
281        -4.51459788337394416547e-12L,
282         1.85594911495471785253e-11L,
283        -7.95748924447710747776e-11L,
284         3.57739728140030116597e-10L,
285        -1.69753450938905987466e-9L,
286         8.57403401741422608519e-9L,
287        -4.66048989768794782956e-8L,
288         2.76681363944501510342e-7L,
289        -1.83175552271911948767e-6L,
290         1.39498137188764993662e-5L,
291        -1.28495495816278026384e-4L,
292         1.56988388573005337491e-3L,
293        -3.14481013119645005427e-2L,
294         2.44030308206595545468e0L   ];
295 
296     if (x <= 2)
297         return chebEval(x^^2 - 2, chebCoeffsA) - log(x/2) * besselI0(x);
298     else
299         return exp(-x) * chebEval(8/x - 2, chebCoeffsB) / sqrt(x);
300 }
301 
302 
303 unittest
304 {
305     assert (approxEqual(besselK0( 1), 4.2102443824070833334e-1, 1.3e-16));
306     assert (approxEqual(besselK0(10), 1.7780062316167651811e-5, 1.3e-16));
307 }
308 
309 
310 
311 
312 /** Modified Bessel function of the second kind, order 1. */
313 real besselK1(real x) pure
314 in { assert(x > 0, "Can only calculate K1 for positive argument"); }
315 body {
316     /* Chebyshev coefficients for x(K1(x) - log(x/2) I1(x))
317      * in the interval [0,2].
318      * 
319      * lim(x->0){ x(K1(x) - log(x/2) I1(x)) } = 1.
320      */
321     immutable real[11] chebCoeffsA = [
322        -7.02386347938628759343e-18L,
323        -2.42744985051936593393e-15L,
324        -6.66690169419932900609e-13L,
325        -1.41148839263352776110e-10L,
326        -2.21338763073472585583e-8L,
327        -2.43340614156596823496e-6L,
328        -1.73028895751305206302e-4L,
329        -6.97572385963986435018e-3L,
330        -1.22611180822657148235e-1L,
331        -3.53155960776544875667e-1L,
332         1.52530022733894777053e0L  ];
333 
334     /* Chebyshev coefficients for exp(x) sqrt(x) K1(x)
335      * in the interval [2,infinity].
336      *
337      * lim(x->inf){ exp(x) sqrt(x) K1(x) } = sqrt(pi/2).
338      */
339     immutable real[25] chebCoeffsB = [
340        -5.75674448366501715755e-18L,
341         1.79405087314755922667e-17L,
342        -5.68946255844285935196e-17L,
343         1.83809354436663880070e-16L,
344        -6.05704724837331885336e-16L,
345         2.03870316562433424052e-15L,
346        -7.01983709041831346144e-15L,
347         2.47715442448130437068e-14L,
348        -8.97670518232499435011e-14L,
349         3.34841966607842919884e-13L,
350        -1.28917396095102890680e-12L,
351         5.13963967348173025100e-12L,
352        -2.12996783842756842877e-11L,
353         9.21831518760500529508e-11L,
354        -4.19035475934189648750e-10L,
355         2.01504975519703286596e-9L,
356        -1.03457624656780970260e-8L,
357         5.74108412545004946722e-8L,
358        -3.50196060308781257119e-7L,
359         2.40648494783721712015e-6L,
360        -1.93619797416608296024e-5L,
361         1.95215518471351631108e-4L,
362        -2.85781685962277938680e-3L,
363         1.03923736576817238437e-1L,
364         2.72062619048444266945e0L   ];
365 
366     real z = 0.5 * x;
367     if (x <= 2.0)
368         return log(z) * besselI1(x) + chebEval(x^^2 - 2, chebCoeffsA) / x;
369     else
370         return exp(-x) * chebEval(8.0/x - 2.0, chebCoeffsB) / sqrt(x);
371 }
372 
373 
374 unittest
375 {
376     assert (approxEqual(besselK1( 1), 6.0190723019723457474e-1, 8.9e-17));
377     assert (approxEqual(besselK1(10), 1.8648773453825584597e-5, 8.9e-17));
378 }
379