Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/2.1.2.2591/ALGLIB-2.1.2.2591/dawson.cs @ 2645

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

extracted external libraries and adapted dependent plugins (ticket #837)

File size: 6.0 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 dawson
33    {
34        /*************************************************************************
35        Dawson's Integral
36
37        Approximates the integral
38
39                                    x
40                                    -
41                             2     | |        2
42         dawsn(x)  =  exp( -x  )   |    exp( t  ) dt
43                                 | |
44                                  -
45                                  0
46
47        Three different rational approximations are employed, for
48        the intervals 0 to 3.25; 3.25 to 6.25; and 6.25 up.
49
50        ACCURACY:
51
52                             Relative error:
53        arithmetic   domain     # trials      peak         rms
54           IEEE      0,10        10000       6.9e-16     1.0e-16
55
56        Cephes Math Library Release 2.8:  June, 2000
57        Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
58        *************************************************************************/
59        public static double dawsonintegral(double x)
60        {
61            double result = 0;
62            double x2 = 0;
63            double y = 0;
64            int sg = 0;
65            double an = 0;
66            double ad = 0;
67            double bn = 0;
68            double bd = 0;
69            double cn = 0;
70            double cd = 0;
71
72            sg = 1;
73            if( (double)(x)<(double)(0) )
74            {
75                sg = -1;
76                x = -x;
77            }
78            if( (double)(x)<(double)(3.25) )
79            {
80                x2 = x*x;
81                an = 1.13681498971755972054E-11;
82                an = an*x2+8.49262267667473811108E-10;
83                an = an*x2+1.94434204175553054283E-8;
84                an = an*x2+9.53151741254484363489E-7;
85                an = an*x2+3.07828309874913200438E-6;
86                an = an*x2+3.52513368520288738649E-4;
87                an = an*x2+-8.50149846724410912031E-4;
88                an = an*x2+4.22618223005546594270E-2;
89                an = an*x2+-9.17480371773452345351E-2;
90                an = an*x2+9.99999999999999994612E-1;
91                ad = 2.40372073066762605484E-11;
92                ad = ad*x2+1.48864681368493396752E-9;
93                ad = ad*x2+5.21265281010541664570E-8;
94                ad = ad*x2+1.27258478273186970203E-6;
95                ad = ad*x2+2.32490249820789513991E-5;
96                ad = ad*x2+3.25524741826057911661E-4;
97                ad = ad*x2+3.48805814657162590916E-3;
98                ad = ad*x2+2.79448531198828973716E-2;
99                ad = ad*x2+1.58874241960120565368E-1;
100                ad = ad*x2+5.74918629489320327824E-1;
101                ad = ad*x2+1.00000000000000000539E0;
102                y = x*an/ad;
103                result = sg*y;
104                return result;
105            }
106            x2 = 1.0/(x*x);
107            if( (double)(x)<(double)(6.25) )
108            {
109                bn = 5.08955156417900903354E-1;
110                bn = bn*x2-2.44754418142697847934E-1;
111                bn = bn*x2+9.41512335303534411857E-2;
112                bn = bn*x2-2.18711255142039025206E-2;
113                bn = bn*x2+3.66207612329569181322E-3;
114                bn = bn*x2-4.23209114460388756528E-4;
115                bn = bn*x2+3.59641304793896631888E-5;
116                bn = bn*x2-2.14640351719968974225E-6;
117                bn = bn*x2+9.10010780076391431042E-8;
118                bn = bn*x2-2.40274520828250956942E-9;
119                bn = bn*x2+3.59233385440928410398E-11;
120                bd = 1.00000000000000000000E0;
121                bd = bd*x2-6.31839869873368190192E-1;
122                bd = bd*x2+2.36706788228248691528E-1;
123                bd = bd*x2-5.31806367003223277662E-2;
124                bd = bd*x2+8.48041718586295374409E-3;
125                bd = bd*x2-9.47996768486665330168E-4;
126                bd = bd*x2+7.81025592944552338085E-5;
127                bd = bd*x2-4.55875153252442634831E-6;
128                bd = bd*x2+1.89100358111421846170E-7;
129                bd = bd*x2-4.91324691331920606875E-9;
130                bd = bd*x2+7.18466403235734541950E-11;
131                y = 1.0/x+x2*bn/(bd*x);
132                result = sg*0.5*y;
133                return result;
134            }
135            if( (double)(x)>(double)(1.0E9) )
136            {
137                result = sg*0.5/x;
138                return result;
139            }
140            cn = -5.90592860534773254987E-1;
141            cn = cn*x2+6.29235242724368800674E-1;
142            cn = cn*x2-1.72858975380388136411E-1;
143            cn = cn*x2+1.64837047825189632310E-2;
144            cn = cn*x2-4.86827613020462700845E-4;
145            cd = 1.00000000000000000000E0;
146            cd = cd*x2-2.69820057197544900361E0;
147            cd = cd*x2+1.73270799045947845857E0;
148            cd = cd*x2-3.93708582281939493482E-1;
149            cd = cd*x2+3.44278924041233391079E-2;
150            cd = cd*x2-9.73655226040941223894E-4;
151            y = 1.0/x+x2*cn/(cd*x);
152            result = sg*0.5*y;
153            return result;
154        }
155    }
156}
Note: See TracBrowser for help on using the repository browser.