2 * Karinto Library Project
\r
4 * This software is distributed under a zlib-style license.
\r
5 * See license.txt for more information.
\r
9 using System.Collections.Generic;
\r
14 public static class MathEx
\r
17 private static readonly Polynomial[] erfcTable =
\r
19 new Polynomial( //[0]: 1 <= x < 2
\r
20 4.7019003296642450457550391704e-1,
\r
21 -5.4129916560282193120142282350e-4,
\r
22 3.2466045488989207057247716494e-2,
\r
23 -1.5551574967900995569013129661e-3,
\r
24 1.2991371652846980287617056420e-3,
\r
25 -1.2476738559919909068871115854e-4,
\r
26 4.2164948133638407169423314936e-5,
\r
27 -5.5944555649744218134732620248e-6,
\r
28 1.2391258283306184550290288637e-6,
\r
29 -1.8629332352031266041799875695e-7,
\r
30 3.3508489330180889238527074083e-8,
\r
31 -5.1453759994172323811411458769e-9,
\r
32 8.2778054007509566713316121080e-10,
\r
33 -1.2419584612081402943997565209e-10,
\r
34 1.8628225410158766603724913483e-11,
\r
35 -2.6916642541339090011279256698e-12,
\r
36 3.8321798025710167662264220379e-13,
\r
37 -5.2985573947392393836199793226e-14,
\r
38 7.2214593891835651145391699616e-15,
\r
39 -1.0552316968012789176241975783e-15,
\r
40 1.3848365749828243613160210837e-16
\r
42 new Polynomial( //[1]: 2 <= x < 4
\r
43 0.28246919208993532,//2.8246919208993525100243945507e-1,
\r
44 1.0295429295263581287408466969e-4,
\r
45 2.6352174925955979906965931429e-2,
\r
46 -2.0840262584872572824603082922e-3,
\r
47 1.6107317263933041407504226538e-3,
\r
48 -2.6413257091217055241174636048e-4,
\r
49 9.3720096524725604112989062675e-5,
\r
50 -2.0394452777501772979594808178e-5,
\r
51 5.4818843984699154429750547408e-6,
\r
52 -1.2769242863171188612896799307e-6,
\r
53 3.0855108190100308040276886416e-7,
\r
54 -7.1110449954162265412516889674e-8,
\r
55 1.6277854940656666077194365644e-8,
\r
56 -3.6388350444865150940102551512e-9,
\r
57 8.0125451384130244729739875906e-10,
\r
58 -1.7333139256038089278659882792e-10,
\r
59 3.6913251414732819339098560977e-11,
\r
60 -7.7391112861863152629855204763e-12,
\r
61 1.5985771117554783386097643078e-12,
\r
62 -3.2558812235239520370430316160e-13,
\r
63 6.5474397002885039401566856717e-14,
\r
64 -1.2908539064863451579569117003e-14,
\r
65 2.4663418649675516649805424542e-15,
\r
66 -5.0102905606542145531985288605e-16,
\r
67 1.1637629874200535139854357149e-16,
\r
68 -1.7702372520803693617903759643e-17
\r
70 new Polynomial( //[2]: 4 <= x < 8
\r
71 2.6609916698112647641230096605e-1,
\r
72 1.0052997692961233689565319255e-1,
\r
73 0.063819586227166714,//6.3819586227166696717161903119e-2,
\r
74 1.6510111747337704484097091109e-2,
\r
75 6.7745775475928297665975539362e-3,
\r
76 1.2613799574673550665257600345e-3,
\r
77 4.6314050254419785813179921962e-4,
\r
78 5.4816124629655721757719475954e-5,
\r
79 2.5069740734847835425277961803e-5,
\r
80 7.6494811710229433994592919775e-7,
\r
81 1.3069891449561428965827894398e-6,
\r
82 -1.1211722505095090109921762128e-7,
\r
83 7.7627072370605473720667209070e-8,
\r
84 -1.5067764445724103557487307330e-8,
\r
85 5.3921138750029687859330065848e-9,
\r
86 -1.3591247804716958040748207366e-9,
\r
87 4.0341965879492537779635899830e-10,
\r
88 -1.0884398300756496354533131958e-10,
\r
89 3.0363217123340332343882922235e-11,
\r
90 -8.2439513886349055414081918841e-12,
\r
91 2.2756305305902268859393373330e-12,
\r
92 -6.1246789459726400033402401708e-13,
\r
93 1.4333003386843162197015125496e-13,
\r
94 -3.7815369812690037058307445308e-14,
\r
95 1.7401564599128446023009328203e-14,
\r
96 -4.5875278594765951588791904401e-15
\r
98 new Polynomial( //[3]: 8 <= x < 16
\r
99 4.6854221014893762619745618301e-2,
\r
100 -1.5511450952249084104775598329e-2,
\r
101 5.1178905303441648577722059649e-3,
\r
102 -1.6829798529769531049332840624e-3,
\r
103 5.5160777130644620308102757353e-4,
\r
104 -1.8020184996880086091090661904e-4,
\r
105 5.8678514133519807272309470851e-5,
\r
106 -1.9045977453678276537004114523e-5,
\r
107 6.1623270905278656151589976968e-6,
\r
108 -1.9875419915464967477768116933e-6,
\r
109 6.3904356643246431611609250453e-7,
\r
110 -2.0483278667408340516487849749e-7,
\r
111 6.5453904816333940501972558207e-8,
\r
112 -2.0852118199579889358624923805e-8,
\r
113 6.6229050456699429842283832512e-9,
\r
114 -2.0972627311983894339808366077e-9,
\r
115 6.6237821648144864340143708663e-10,
\r
116 -2.0853485442305791192334927397e-10,
\r
117 6.5162153190849185531320932146e-11,
\r
118 -2.0380108335416826794297357698e-11,
\r
119 6.6463403803777970165860627602e-12,
\r
120 -2.0777605497301746376639943675e-12,
\r
121 4.6614566309518191650253075708e-13,
\r
122 -1.4074713019480904125485292343e-13,
\r
123 1.0831110710204089408511450325e-13,
\r
124 -3.3906934201035073659517665140e-14
\r
126 new Polynomial( //[4]: 16 <= x < 28
\r
127 2.8262089312268317148835811995e-2,
\r
128 -5.4152623609484073274229170763e-3,
\r
129 1.6241399412435775154107409724e-3,
\r
130 -5.0744105591774430405597196162e-3,
\r
131 -1.7616157737971578130497741968e-2,
\r
132 2.3083427166129394399962687197e-2,
\r
133 1.1164742403248709989628167597e-1,
\r
134 -3.5273994341997794105095275190e-2,
\r
135 -3.0869428535741202544195312887e-1,
\r
136 -1.3020536214165453885942308553e-2,
\r
137 4.3332671745834241123676787595e-1,
\r
138 9.4385604975718760529994282754e-2,
\r
139 -3.0238469418692319645741714650e-1,
\r
140 -9.5995136295717931771565803096e-2,
\r
141 8.3313274670729987166299100438e-2,
\r
142 3.1593997744278794785645602547e-2
\r
147 #region public methods
\r
148 public static double Erf(double x)
\r
150 double a = Math.Abs(x);
\r
153 return 1.0 - Erfc(x);
\r
159 double c1 = a * aa;
\r
160 double c2 = c1 * aa;
\r
161 double c3 = c2 * aa;
\r
162 double c4 = c3 * aa;
\r
164 double s1 = c1 * -315.0;// (3*5*7*9 / -3 / 1!);
\r
165 double s2 = s1 * 0 + c2 * 189.0 / 2.0;// (3*5*7*9 / +5 / 2!)
\r
166 double s3 = s2 * 0 + c3 * -45.0 / 2.0;// (3*5*7*9 / -7 / 3!)
\r
167 double s4 = s3 * 0 + c4 * 35.0 / 8.0;// (3*5*7*9 / 9 / 4!)
\r
168 y = a + (s1 + s2 + s3 + s4) / 945.0; //(3*5*7*9);
\r
175 for (double i = -1.0; i > -20.0; i-=1.0)
\r
183 y *= 1.1283791670955126; // 2 / pi^0.5
\r
184 return Math.Sign(x) > 0 ? y : -y;
\r
187 public static double Erfc(double x)
\r
189 double a = Math.Abs(x);
\r
192 return 1.0 - Erf(x);
\r
194 if (a > 27.226017111108501)
\r
196 return x > 0 ? 0 : 2.0;
\r
202 y = Math.Exp(-1.1688327773406599 * a * a);
\r
203 y *= erfcTable[0].GetValueAt(a + a - 3.0);
\r
207 decimal b = -1.05068636145441m * (decimal)a * (decimal)a;
\r
208 y = Math.Exp((double)b);
\r
209 y *= erfcTable[1].GetValueAt(a - 3.0);
\r
213 decimal b = -1.0292687483818601m * (decimal)a * (decimal)a;
\r
214 y = Math.Exp((double)b);
\r
215 y *= erfcTable[2].GetValueAt(a * 0.5 - 3);
\r
219 y = Math.Exp(-a * a);
\r
220 y *= erfcTable[3].GetValueAt(a * 0.25 - 3);
\r
224 y = Math.Exp(-a * a);
\r
225 y *= erfcTable[4].GetValueAt(a * 0.125 - 3);
\r
227 return x > 0 ? y : 2.0 - y;
\r