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