Free cookie consent management tool by TermsFeed Policy Generator

source: tags/3.3.2/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/2.5.0/ALGLIB-2.5.0/airyf.cs @ 13398

Last change on this file since 13398 was 3839, checked in by mkommend, 14 years ago

implemented first version of LR (ticket #1012)
File size: 14.7 KB
Line 
1/*************************************************************************
2Cephes Math Library Release 2.8:  June, 2000
3Copyright by Stephen L. Moshier
4
5Contributors:
6    * Sergey Bochkanov (ALGLIB project). Translation from C to
7      pseudocode.
8
9See subroutines comments for additional copyrights.
10
11>>> SOURCE LICENSE >>>
12This program is free software; you can redistribute it and/or modify
13it under the terms of the GNU General Public License as published by
14the Free Software Foundation (www.fsf.org); either version 2 of the
15License, or (at your option) any later version.
16
17This program is distributed in the hope that it will be useful,
18but WITHOUT ANY WARRANTY; without even the implied warranty of
19MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
20GNU General Public License for more details.
21
22A copy of the GNU General Public License is available at
23http://www.fsf.org/licensing/licenses
24
25>>> END OF LICENSE >>>
26*************************************************************************/
27
28using System;
29
30namespace alglib
31{
32    public class airyf
33    {
34        /*************************************************************************
35        Airy function
36
37        Solution of the differential equation
38
39        y"(x) = xy.
40
41        The function returns the two independent solutions Ai, Bi
42        and their first derivatives Ai'(x), Bi'(x).
43
44        Evaluation is by power series summation for small x,
45        by rational minimax approximations for large x.
46
47
48
49        ACCURACY:
50        Error criterion is absolute when function <= 1, relative
51        when function > 1, except * denotes relative error criterion.
52        For large negative x, the absolute error increases as x^1.5.
53        For large positive x, the relative error increases as x^1.5.
54
55        Arithmetic  domain   function  # trials      peak         rms
56        IEEE        -10, 0     Ai        10000       1.6e-15     2.7e-16
57        IEEE          0, 10    Ai        10000       2.3e-14*    1.8e-15*
58        IEEE        -10, 0     Ai'       10000       4.6e-15     7.6e-16
59        IEEE          0, 10    Ai'       10000       1.8e-14*    1.5e-15*
60        IEEE        -10, 10    Bi        30000       4.2e-15     5.3e-16
61        IEEE        -10, 10    Bi'       30000       4.9e-15     7.3e-16
62
63        Cephes Math Library Release 2.8:  June, 2000
64        Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
65        *************************************************************************/
66        public static void airy(double x,
67            ref double ai,
68            ref double aip,
69            ref double bi,
70            ref double bip)
71        {
72            double z = 0;
73            double zz = 0;
74            double t = 0;
75            double f = 0;
76            double g = 0;
77            double uf = 0;
78            double ug = 0;
79            double k = 0;
80            double zeta = 0;
81            double theta = 0;
82            int domflg = 0;
83            double c1 = 0;
84            double c2 = 0;
85            double sqrt3 = 0;
86            double sqpii = 0;
87            double afn = 0;
88            double afd = 0;
89            double agn = 0;
90            double agd = 0;
91            double apfn = 0;
92            double apfd = 0;
93            double apgn = 0;
94            double apgd = 0;
95            double an = 0;
96            double ad = 0;
97            double apn = 0;
98            double apd = 0;
99            double bn16 = 0;
100            double bd16 = 0;
101            double bppn = 0;
102            double bppd = 0;
103
104            sqpii = 5.64189583547756286948E-1;
105            c1 = 0.35502805388781723926;
106            c2 = 0.258819403792806798405;
107            sqrt3 = 1.732050807568877293527;
108            domflg = 0;
109            if( (double)(x)>(double)(25.77) )
110            {
111                ai = 0;
112                aip = 0;
113                bi = AP.Math.MaxRealNumber;
114                bip = AP.Math.MaxRealNumber;
115                return;
116            }
117            if( (double)(x)<(double)(-2.09) )
118            {
119                domflg = 15;
120                t = Math.Sqrt(-x);
121                zeta = -(2.0*x*t/3.0);
122                t = Math.Sqrt(t);
123                k = sqpii/t;
124                z = 1.0/zeta;
125                zz = z*z;
126                afn = -1.31696323418331795333E-1;
127                afn = afn*zz-6.26456544431912369773E-1;
128                afn = afn*zz-6.93158036036933542233E-1;
129                afn = afn*zz-2.79779981545119124951E-1;
130                afn = afn*zz-4.91900132609500318020E-2;
131                afn = afn*zz-4.06265923594885404393E-3;
132                afn = afn*zz-1.59276496239262096340E-4;
133                afn = afn*zz-2.77649108155232920844E-6;
134                afn = afn*zz-1.67787698489114633780E-8;
135                afd = 1.00000000000000000000E0;
136                afd = afd*zz+1.33560420706553243746E1;
137                afd = afd*zz+3.26825032795224613948E1;
138                afd = afd*zz+2.67367040941499554804E1;
139                afd = afd*zz+9.18707402907259625840E0;
140                afd = afd*zz+1.47529146771666414581E0;
141                afd = afd*zz+1.15687173795188044134E-1;
142                afd = afd*zz+4.40291641615211203805E-3;
143                afd = afd*zz+7.54720348287414296618E-5;
144                afd = afd*zz+4.51850092970580378464E-7;
145                uf = 1.0+zz*afn/afd;
146                agn = 1.97339932091685679179E-2;
147                agn = agn*zz+3.91103029615688277255E-1;
148                agn = agn*zz+1.06579897599595591108E0;
149                agn = agn*zz+9.39169229816650230044E-1;
150                agn = agn*zz+3.51465656105547619242E-1;
151                agn = agn*zz+6.33888919628925490927E-2;
152                agn = agn*zz+5.85804113048388458567E-3;
153                agn = agn*zz+2.82851600836737019778E-4;
154                agn = agn*zz+6.98793669997260967291E-6;
155                agn = agn*zz+8.11789239554389293311E-8;
156                agn = agn*zz+3.41551784765923618484E-10;
157                agd = 1.00000000000000000000E0;
158                agd = agd*zz+9.30892908077441974853E0;
159                agd = agd*zz+1.98352928718312140417E1;
160                agd = agd*zz+1.55646628932864612953E1;
161                agd = agd*zz+5.47686069422975497931E0;
162                agd = agd*zz+9.54293611618961883998E-1;
163                agd = agd*zz+8.64580826352392193095E-2;
164                agd = agd*zz+4.12656523824222607191E-3;
165                agd = agd*zz+1.01259085116509135510E-4;
166                agd = agd*zz+1.17166733214413521882E-6;
167                agd = agd*zz+4.91834570062930015649E-9;
168                ug = z*agn/agd;
169                theta = zeta+0.25*Math.PI;
170                f = Math.Sin(theta);
171                g = Math.Cos(theta);
172                ai = k*(f*uf-g*ug);
173                bi = k*(g*uf+f*ug);
174                apfn = 1.85365624022535566142E-1;
175                apfn = apfn*zz+8.86712188052584095637E-1;
176                apfn = apfn*zz+9.87391981747398547272E-1;
177                apfn = apfn*zz+4.01241082318003734092E-1;
178                apfn = apfn*zz+7.10304926289631174579E-2;
179                apfn = apfn*zz+5.90618657995661810071E-3;
180                apfn = apfn*zz+2.33051409401776799569E-4;
181                apfn = apfn*zz+4.08718778289035454598E-6;
182                apfn = apfn*zz+2.48379932900442457853E-8;
183                apfd = 1.00000000000000000000E0;
184                apfd = apfd*zz+1.47345854687502542552E1;
185                apfd = apfd*zz+3.75423933435489594466E1;
186                apfd = apfd*zz+3.14657751203046424330E1;
187                apfd = apfd*zz+1.09969125207298778536E1;
188                apfd = apfd*zz+1.78885054766999417817E0;
189                apfd = apfd*zz+1.41733275753662636873E-1;
190                apfd = apfd*zz+5.44066067017226003627E-3;
191                apfd = apfd*zz+9.39421290654511171663E-5;
192                apfd = apfd*zz+5.65978713036027009243E-7;
193                uf = 1.0+zz*apfn/apfd;
194                apgn = -3.55615429033082288335E-2;
195                apgn = apgn*zz-6.37311518129435504426E-1;
196                apgn = apgn*zz-1.70856738884312371053E0;
197                apgn = apgn*zz-1.50221872117316635393E0;
198                apgn = apgn*zz-5.63606665822102676611E-1;
199                apgn = apgn*zz-1.02101031120216891789E-1;
200                apgn = apgn*zz-9.48396695961445269093E-3;
201                apgn = apgn*zz-4.60325307486780994357E-4;
202                apgn = apgn*zz-1.14300836484517375919E-5;
203                apgn = apgn*zz-1.33415518685547420648E-7;
204                apgn = apgn*zz-5.63803833958893494476E-10;
205                apgd = 1.00000000000000000000E0;
206                apgd = apgd*zz+9.85865801696130355144E0;
207                apgd = apgd*zz+2.16401867356585941885E1;
208                apgd = apgd*zz+1.73130776389749389525E1;
209                apgd = apgd*zz+6.17872175280828766327E0;
210                apgd = apgd*zz+1.08848694396321495475E0;
211                apgd = apgd*zz+9.95005543440888479402E-2;
212                apgd = apgd*zz+4.78468199683886610842E-3;
213                apgd = apgd*zz+1.18159633322838625562E-4;
214                apgd = apgd*zz+1.37480673554219441465E-6;
215                apgd = apgd*zz+5.79912514929147598821E-9;
216                ug = z*apgn/apgd;
217                k = sqpii*t;
218                aip = -(k*(g*uf+f*ug));
219                bip = k*(f*uf-g*ug);
220                return;
221            }
222            if( (double)(x)>=(double)(2.09) )
223            {
224                domflg = 5;
225                t = Math.Sqrt(x);
226                zeta = 2.0*x*t/3.0;
227                g = Math.Exp(zeta);
228                t = Math.Sqrt(t);
229                k = 2.0*t*g;
230                z = 1.0/zeta;
231                an = 3.46538101525629032477E-1;
232                an = an*z+1.20075952739645805542E1;
233                an = an*z+7.62796053615234516538E1;
234                an = an*z+1.68089224934630576269E2;
235                an = an*z+1.59756391350164413639E2;
236                an = an*z+7.05360906840444183113E1;
237                an = an*z+1.40264691163389668864E1;
238                an = an*z+9.99999999999999995305E-1;
239                ad = 5.67594532638770212846E-1;
240                ad = ad*z+1.47562562584847203173E1;
241                ad = ad*z+8.45138970141474626562E1;
242                ad = ad*z+1.77318088145400459522E2;
243                ad = ad*z+1.64234692871529701831E2;
244                ad = ad*z+7.14778400825575695274E1;
245                ad = ad*z+1.40959135607834029598E1;
246                ad = ad*z+1.00000000000000000470E0;
247                f = an/ad;
248                ai = sqpii*f/k;
249                k = -(0.5*sqpii*t/g);
250                apn = 6.13759184814035759225E-1;
251                apn = apn*z+1.47454670787755323881E1;
252                apn = apn*z+8.20584123476060982430E1;
253                apn = apn*z+1.71184781360976385540E2;
254                apn = apn*z+1.59317847137141783523E2;
255                apn = apn*z+6.99778599330103016170E1;
256                apn = apn*z+1.39470856980481566958E1;
257                apn = apn*z+1.00000000000000000550E0;
258                apd = 3.34203677749736953049E-1;
259                apd = apd*z+1.11810297306158156705E1;
260                apd = apd*z+7.11727352147859965283E1;
261                apd = apd*z+1.58778084372838313640E2;
262                apd = apd*z+1.53206427475809220834E2;
263                apd = apd*z+6.86752304592780337944E1;
264                apd = apd*z+1.38498634758259442477E1;
265                apd = apd*z+9.99999999999999994502E-1;
266                f = apn/apd;
267                aip = f*k;
268                if( (double)(x)>(double)(8.3203353) )
269                {
270                    bn16 = -2.53240795869364152689E-1;
271                    bn16 = bn16*z+5.75285167332467384228E-1;
272                    bn16 = bn16*z-3.29907036873225371650E-1;
273                    bn16 = bn16*z+6.44404068948199951727E-2;
274                    bn16 = bn16*z-3.82519546641336734394E-3;
275                    bd16 = 1.00000000000000000000E0;
276                    bd16 = bd16*z-7.15685095054035237902E0;
277                    bd16 = bd16*z+1.06039580715664694291E1;
278                    bd16 = bd16*z-5.23246636471251500874E0;
279                    bd16 = bd16*z+9.57395864378383833152E-1;
280                    bd16 = bd16*z-5.50828147163549611107E-2;
281                    f = z*bn16/bd16;
282                    k = sqpii*g;
283                    bi = k*(1.0+f)/t;
284                    bppn = 4.65461162774651610328E-1;
285                    bppn = bppn*z-1.08992173800493920734E0;
286                    bppn = bppn*z+6.38800117371827987759E-1;
287                    bppn = bppn*z-1.26844349553102907034E-1;
288                    bppn = bppn*z+7.62487844342109852105E-3;
289                    bppd = 1.00000000000000000000E0;
290                    bppd = bppd*z-8.70622787633159124240E0;
291                    bppd = bppd*z+1.38993162704553213172E1;
292                    bppd = bppd*z-7.14116144616431159572E0;
293                    bppd = bppd*z+1.34008595960680518666E0;
294                    bppd = bppd*z-7.84273211323341930448E-2;
295                    f = z*bppn/bppd;
296                    bip = k*t*(1.0+f);
297                    return;
298                }
299            }
300            f = 1.0;
301            g = x;
302            t = 1.0;
303            uf = 1.0;
304            ug = x;
305            k = 1.0;
306            z = x*x*x;
307            while( (double)(t)>(double)(AP.Math.MachineEpsilon) )
308            {
309                uf = uf*z;
310                k = k+1.0;
311                uf = uf/k;
312                ug = ug*z;
313                k = k+1.0;
314                ug = ug/k;
315                uf = uf/k;
316                f = f+uf;
317                k = k+1.0;
318                ug = ug/k;
319                g = g+ug;
320                t = Math.Abs(uf/f);
321            }
322            uf = c1*f;
323            ug = c2*g;
324            if( domflg%2==0 )
325            {
326                ai = uf-ug;
327            }
328            if( domflg/2%2==0 )
329            {
330                bi = sqrt3*(uf+ug);
331            }
332            k = 4.0;
333            uf = x*x/2.0;
334            ug = z/3.0;
335            f = uf;
336            g = 1.0+ug;
337            uf = uf/3.0;
338            t = 1.0;
339            while( (double)(t)>(double)(AP.Math.MachineEpsilon) )
340            {
341                uf = uf*z;
342                ug = ug/k;
343                k = k+1.0;
344                ug = ug*z;
345                uf = uf/k;
346                f = f+uf;
347                k = k+1.0;
348                ug = ug/k;
349                uf = uf/k;
350                g = g+ug;
351                k = k+1.0;
352                t = Math.Abs(ug/g);
353            }
354            uf = c1*f;
355            ug = c2*g;
356            if( domflg/4%2==0 )
357            {
358                aip = uf-ug;
359            }
360            if( domflg/8%2==0 )
361            {
362                bip = sqrt3*(uf+ug);
363            }
364        }
365    }
366}
Note: See TracBrowser for help on using the repository browser.