[2563] | 1 | /*************************************************************************
|
---|
| 2 | Cephes Math Library Release 2.8: June, 2000
|
---|
| 3 | Copyright by Stephen L. Moshier
|
---|
| 4 |
|
---|
| 5 | Contributors:
|
---|
| 6 | * Sergey Bochkanov (ALGLIB project). Translation from C to
|
---|
| 7 | pseudocode.
|
---|
| 8 |
|
---|
| 9 | See subroutines comments for additional copyrights.
|
---|
| 10 |
|
---|
| 11 | >>> SOURCE LICENSE >>>
|
---|
| 12 | This program is free software; you can redistribute it and/or modify
|
---|
| 13 | it under the terms of the GNU General Public License as published by
|
---|
| 14 | the Free Software Foundation (www.fsf.org); either version 2 of the
|
---|
| 15 | License, or (at your option) any later version.
|
---|
| 16 |
|
---|
| 17 | This program is distributed in the hope that it will be useful,
|
---|
| 18 | but WITHOUT ANY WARRANTY; without even the implied warranty of
|
---|
| 19 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
---|
| 20 | GNU General Public License for more details.
|
---|
| 21 |
|
---|
| 22 | A copy of the GNU General Public License is available at
|
---|
| 23 | http://www.fsf.org/licensing/licenses
|
---|
| 24 |
|
---|
| 25 | >>> END OF LICENSE >>>
|
---|
| 26 | *************************************************************************/
|
---|
| 27 |
|
---|
| 28 | using System;
|
---|
| 29 |
|
---|
| 30 | namespace 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 | }
|
---|