Skip to content

Commit 05757d4

Browse files
tk-yoshimuratk-yoshimura
committed
+ erf pade
1 parent 8da730b commit 05757d4

5 files changed

Lines changed: 359 additions & 55 deletions

File tree

‎ErrorFunctionApproximation/ErrorFunctionMarkII.cs‎

Lines changed: 252 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,252 @@
1+
using static System.Math;
2+
3+
namespace ErrorFunctionApproximation {
4+
public static class ErrorFunctionMarkII {
5+
public static double Erf(double x) {
6+
if (double.IsNaN(x)) {
7+
return double.NaN;
8+
}
9+
if (x < 0.0) {
10+
return -Erf(Abs(x));
11+
}
12+
13+
if (x < 0.5) {
14+
return ErfNearZero(x);
15+
}
16+
if (x < 1.0) {
17+
return 1.0 - ErfcPlusP5(x);
18+
}
19+
if (x < 2.0) {
20+
return 1.0 - ErfcPlus1(x);
21+
}
22+
if (x < 4.0) {
23+
return 1.0 - ErfcPlus2(x);
24+
}
25+
if (x < 8.0) {
26+
return 1.0 - ErfcPlus4(x);
27+
}
28+
29+
return 1.0;
30+
}
31+
32+
public static double Erfc(double x) {
33+
if (double.IsNaN(x)) {
34+
return double.NaN;
35+
}
36+
37+
if (x < 0.5) {
38+
return 1.0 - Erf(x);
39+
}
40+
if (x < 1.0) {
41+
return ErfcPlusP5(x);
42+
}
43+
if (x < 2.0) {
44+
return ErfcPlus1(x);
45+
}
46+
if (x < 4.0) {
47+
return ErfcPlus2(x);
48+
}
49+
if (x < 8.0) {
50+
return ErfcPlus4(x);
51+
}
52+
if (x < 16.0) {
53+
return ErfcPlus8(x);
54+
}
55+
if (x < 27.25) {
56+
return ErfcPlus16(x);
57+
}
58+
59+
return 0.0;
60+
}
61+
62+
private static double ErfNearZero(double x) {
63+
double w = x * x;
64+
65+
double u = (
66+
1.12837916709551257390e0 + w * (
67+
1.41992244117798429234e-1 + w * (
68+
4.52195919383661078468e-2 + w * (
69+
1.81571926199746178752e-3 + w *
70+
1.92899332746423909221e-4))))
71+
/ (
72+
1.00000000000000000000e0 + w * (
73+
4.59170683275987299444e-1 + w * (
74+
9.31317143590956177199e-2 + w * (
75+
1.05454995673356468502e-2 + w * (
76+
6.75953355582711700711e-4 + w *
77+
1.99751561392946360834e-5)))));
78+
79+
double y = x * u;
80+
81+
return y;
82+
}
83+
84+
private static double ErfcPlusP5(double x) {
85+
double w = x - 0.5;
86+
87+
double u = (
88+
1.00000000000000000000e0 + w * (
89+
1.34099976999121329590e0 + w * (
90+
8.56762978144288861443e-1 + w * (
91+
3.16570587181473309350e-1 + w * (
92+
7.05798099110973092777e-2 + w * (
93+
8.92998004186820217446e-3 + w *
94+
5.00239328430236539100e-4))))))
95+
/ (
96+
1.62419308574807066971e0 + w * (
97+
3.53051729947024079596e0 + w * (
98+
3.38347446398070396153e0 + w * (
99+
1.85669816101497046144e0 + w * (
100+
6.31463529762465612180e-1 + w * (
101+
1.33441854070762465425e-1 + w * (
102+
1.62726723432911469379e-2 + w *
103+
8.86583249851894459219e-4)))))));
104+
105+
double y = Exp(-x * x) * u;
106+
107+
return y;
108+
}
109+
110+
private static double ErfcPlus1(double x) {
111+
double w = x - 1.0;
112+
113+
double u = (
114+
1.00000000000000000000e0 + w * (
115+
1.50314224721654127642e0 + w * (
116+
1.04476336774122843193e0 + w * (
117+
4.26889879050711362051e-1 + w * (
118+
1.09935724575582226688e-1 + w * (
119+
1.77735261175907078409e-2 + w * (
120+
1.66761213238888643314e-3 + w *
121+
7.00293677081427709050e-5)))))))
122+
/ (
123+
2.33872406651000647661e0 + w * (
124+
5.00980365221224034376e0 + w * (
125+
4.80015965515186399462e0 + w * (
126+
2.69016407758612841845e0 + w * (
127+
9.65783251517281689174e-1 + w * (
128+
2.27773183728755448353e-1 + w * (
129+
3.45207187382309254048e-2 + w * (
130+
3.07987987255160486906e-3 + w *
131+
1.24124178972110154556e-4))))))));
132+
133+
double y = Exp(-x * x) * u;
134+
135+
return y;
136+
}
137+
138+
private static double ErfcPlus2(double x) {
139+
double w = x - 2.0;
140+
141+
double u = (
142+
1.00000000000000000000e0 + w * (
143+
1.54789172892709681253e0 + w * (
144+
1.08732691359437197452e0 + w * (
145+
4.50760649484611848536e-1 + w * (
146+
1.20340580111242787211e-1 + w * (
147+
2.11598919778364788400e-2 + w * (
148+
2.39169506541366076969e-3 + w * (
149+
1.58871052785301256315e-4 + w *
150+
4.75003546471038792860e-6))))))))
151+
/ (
152+
3.91549306725230897235e0 + w * (
153+
7.69806507033900444546e0 + w * (
154+
6.83556849290924291721e0 + w * (
155+
3.60033851351004657099e0 + w * (
156+
1.24047612045121162177e0 + w * (
157+
2.90158593990291689591e-1 + w * (
158+
4.61156483147089573007e-2 + w * (
159+
4.80656185549192240090e-3 + w * (
160+
2.98430051654960595459e-4 + w *
161+
8.41921851559063070821e-6)))))))));
162+
163+
double y = Exp(-x * x) * u;
164+
165+
return y;
166+
}
167+
168+
private static double ErfcPlus4(double x) {
169+
double w = x - 4.0;
170+
171+
double u = (
172+
1.00000000000000000000e0 + w * (
173+
1.25746170311257204438e0 + w * (
174+
7.03808332270199051823e-1 + w * (
175+
2.28892390212618636873e-1 + w * (
176+
4.72951835938663487710e-2 + w * (
177+
6.35712058181220292220e-3 + w * (
178+
5.42839661005351321710e-4 + w * (
179+
2.69272815311508576646e-5 + w *
180+
5.94216970374386786069e-7))))))))
181+
/ (
182+
7.29929897048781762900e0 + w * (
183+
1.09039745244909594419e1 + w * (
184+
7.31699851898786526319e0 + w * (
185+
2.89571228902039492291e0 + w * (
186+
7.45055426257682045346e-1 + w * (
187+
1.29292865258882582078e-1 + w * (
188+
1.51380930848769294948e-2 + w * (
189+
1.15359431397539239090e-3 + w * (
190+
5.19402524758985255038e-5 + w *
191+
1.05322215737521862732e-6)))))))));
192+
193+
double y = Exp(-x * x) * u;
194+
195+
return y;
196+
}
197+
198+
private static double ErfcPlus8(double x) {
199+
double w = x - 8.0;
200+
201+
double u = (
202+
1.00000000000000000000e0 + w * (
203+
1.25746170311257204438e0 + w * (
204+
7.03808332270199051823e-1 + w * (
205+
2.28892390212618636873e-1 + w * (
206+
4.72951835938663487710e-2 + w * (
207+
6.35712058181220292220e-3 + w * (
208+
5.42839661005351321710e-4 + w * (
209+
2.69272815311508576646e-5 + w *
210+
5.94216970374386786069e-7))))))))
211+
/ (
212+
7.29929897048781762900e0 + w * (
213+
1.09039745244909594419e1 + w * (
214+
7.31699851898786526319e0 + w * (
215+
2.89571228902039492291e0 + w * (
216+
7.45055426257682045346e-1 + w * (
217+
1.29292865258882582078e-1 + w * (
218+
1.51380930848769294948e-2 + w * (
219+
1.15359431397539239090e-3 + w * (
220+
5.19402524758985255038e-5 + w *
221+
1.05322215737521862732e-6)))))))));
222+
223+
double y = Exp(-x * x) * u;
224+
225+
return y;
226+
}
227+
228+
private static double ErfcPlus16(double x) {
229+
double w = x - 16.0;
230+
231+
double u = (
232+
1.00000000000000000000e0 + w * (
233+
3.01131564001861798828e-1 + w * (
234+
3.63647817146354655538e-2 + w * (
235+
2.20132258260920223731e-3 + w * (
236+
6.67982674643708596687e-5 + w *
237+
8.12863261006534035097e-7)))))
238+
/ (
239+
2.84144365162813186597e1 + w * (
240+
1.03255156563650534155e1 + w * (
241+
1.56620755225060452355e0 + w * (
242+
1.26930452727791373837e-1 + w * (
243+
5.79681261235647022363e-3 + w * (
244+
1.41449048276726826659e-4 + w *
245+
1.44076261723069040035e-6))))));
246+
247+
double y = Exp(-x * x) * u;
248+
249+
return y;
250+
}
251+
}
252+
}

‎ErrorFunctionApproximation/Program.cs‎

Lines changed: 50 additions & 54 deletions
Original file line numberDiff line numberDiff line change
@@ -7,76 +7,72 @@
77
namespace ErrorFunctionApproximation {
88
class Program {
99
static void Main(string[] args) {
10-
MultiPrecision<Pow2.N32>[] xs = new MultiPrecision<Pow2.N32>[]{
11-
0.5, 1, 2, 4, 8, 16, 27.25
12-
};
13-
1410
MultiPrecision<Pow2.N32> dx = 1d / 1024;
1511

16-
using StreamWriter sw = new("../../../../results/pade_erfc_e32.txt");
17-
sw.WriteLine("pade approximant exp(-x^2)/erfc(x)");
12+
using StreamWriter sw = new("../../../../results/pade_erf_e32.txt");
13+
sw.WriteLine("pade approximant f(x^2) f(x):=erf(x)/x");
14+
15+
MultiPrecision<Pow2.N32> x0 = 0, x1 = 0.5;
1816

19-
for (int j = 0; j < xs.Length - 1; j++) {
20-
MultiPrecision<Pow2.N32> x0 = xs[j], x1 = xs[j + 1];
17+
sw.WriteLine($"\nrange x in [{x0}, {x1}]");
2118

22-
sw.WriteLine($"\nrange x in [{x0}, {x1}]");
19+
MultiPrecision<Pow2.N32> c = 2 / MultiPrecision<Pow2.N32>.Sqrt(MultiPrecision<Pow2.N32>.PI);
2320

24-
sw.WriteLine("expected");
25-
List<MultiPrecision<Pow2.N32>> expecteds = new();
26-
for (MultiPrecision<Pow2.N32> x = x0; x <= x1; x += dx) {
27-
MultiPrecision<Pow2.N32> y = MultiPrecision<Pow2.N32>.Exp(-x * x) / MultiPrecision<Pow2.N32>.Erfc(x);
21+
sw.WriteLine("expected");
22+
List<MultiPrecision<Pow2.N32>> expecteds = new();
23+
for (MultiPrecision<Pow2.N32> x = x0; x <= x1; x += dx) {
24+
MultiPrecision<Pow2.N32> y = (x > 0) ? MultiPrecision<Pow2.N32>.Erf(x) / (x * c) : 1;
2825

29-
expecteds.Add(y);
26+
expecteds.Add(y);
3027

31-
//sw.WriteLine($"{x},{y:e40}");
32-
}
28+
//sw.WriteLine($"{x},{y:e40}");
29+
}
3330

34-
sw.WriteLine($"diffs x = {x0}");
35-
List<MultiPrecision<Pow2.N32>> diffs = new();
36-
for (int d = 0; d <= 128; d++) {
37-
MultiPrecision<Pow2.N32> dy = ErfcExpScaled<Pow2.N32>.Diff(d, x0) * MultiPrecision<Pow2.N32>.TaylorSequence[d];
38-
diffs.Add(dy);
31+
sw.WriteLine($"diffs x = {x0}");
32+
List<MultiPrecision<Pow2.N32>> diffs = new();
33+
for (int d = 0; d <= 128; d++) {
34+
MultiPrecision<Pow2.N32> dy = 1 / ErfTaylorSeries.Coef<Pow2.N32>(d);
35+
diffs.Add(dy);
3936

40-
//sw.WriteLine($"{d},{dy:e40}");
41-
}
42-
sw.Flush();
43-
44-
sw.WriteLine("pade results");
45-
bool is_finished = false;
37+
//sw.WriteLine($"{d},{dy:e40}");
38+
}
39+
sw.Flush();
4640

47-
for (int m = 4; m <= 64 && !is_finished; m++) {
48-
for (int n = m - 1; n <= m + 1 && m + n < 128 && !is_finished; n++) {
49-
(MultiPrecision<Pow2.N32>[] ms, MultiPrecision<Pow2.N32>[] ns) =
50-
PadeSolver<Pow2.N32>.Solve(diffs.Take(m + n + 1).ToArray(), m, n);
41+
sw.WriteLine("pade results");
42+
bool is_finished = false;
5143

52-
MultiPrecision<Pow2.N32> err = 0;
53-
for ((int i, MultiPrecision<Pow2.N32> x) = (0, x0); i < expecteds.Count; i++, x += dx) {
54-
MultiPrecision<Pow2.N32> expected = expecteds[i];
55-
MultiPrecision<Pow2.N32> actual = PadeSolver<Pow2.N32>.Approx(x - x0, ms, ns);
44+
for (int m = 4; m <= 64 && !is_finished; m++) {
45+
for (int n = m - 1; n <= m + 1 && m + n < 128 && !is_finished; n++) {
46+
(MultiPrecision<Pow2.N32>[] ms, MultiPrecision<Pow2.N32>[] ns) =
47+
PadeSolver<Pow2.N32>.Solve(diffs.Take(m + n + 1).ToArray(), m, n);
5648

57-
err = MultiPrecision<Pow2.N32>.Max(err, MultiPrecision<Pow2.N32>.Abs(expected / actual - 1));
58-
}
49+
MultiPrecision<Pow2.N32> err = 0;
50+
for ((int i, MultiPrecision<Pow2.N32> x) = (0, x0); i < expecteds.Count; i++, x += dx) {
51+
MultiPrecision<Pow2.N32> expected = expecteds[i];
52+
MultiPrecision<Pow2.N32> actual = PadeSolver<Pow2.N32>.Approx(MultiPrecision<Pow2.N32>.Square(x - x0), ms, ns);
5953

60-
Console.WriteLine($"m={m},n={n}");
61-
Console.WriteLine($"{err:e20}");
54+
err = MultiPrecision<Pow2.N32>.Max(err, MultiPrecision<Pow2.N32>.Abs(expected / actual - 1));
55+
}
6256

63-
if (err < 1e-32 && ms.All((v) => v > 0) && ns.All((v) => v > 0)) {
64-
sw.WriteLine($"m={m},n={n}");
65-
sw.WriteLine("ms");
66-
for (int i = 0; i < ms.Length; i++) {
67-
sw.WriteLine($"{ms[i]:e40}");
68-
}
69-
sw.WriteLine("ns");
70-
for (int i = 0; i < ns.Length; i++) {
71-
sw.WriteLine($"{ns[i]:e40}");
72-
}
73-
sw.WriteLine("relative err");
74-
sw.WriteLine($"{err:e20}");
75-
sw.Flush();
57+
Console.WriteLine($"m={m},n={n}");
58+
Console.WriteLine($"{err:e20}");
7659

77-
is_finished = true;
78-
break;
60+
if (err < 1e-32 && ms.All((v) => v > 0) && ns.All((v) => v > 0)) {
61+
sw.WriteLine($"m={m},n={n}");
62+
sw.WriteLine("ms");
63+
for (int i = 0; i < ms.Length; i++) {
64+
sw.WriteLine($"{(ms[i] * c):e40}");
7965
}
66+
sw.WriteLine("ns");
67+
for (int i = 0; i < ns.Length; i++) {
68+
sw.WriteLine($"{ns[i]:e40}");
69+
}
70+
sw.WriteLine("relative err");
71+
sw.WriteLine($"{err:e20}");
72+
sw.Flush();
73+
74+
is_finished = true;
75+
break;
8076
}
8177
}
8278
}

‎README.md‎

Lines changed: 5 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -118,5 +118,9 @@ larger 27.25:
118118

119119
### Pade Approximant
120120
![pade definition](https://github.com/tk-yoshimura/ErrorFunctionApproximation/blob/main/figures/pade_definition.svg)
121+
[C# code](https://github.com/tk-yoshimura/ErrorFunctionApproximation/blob/main/ErrorFunctionApproximation/ErrorFunctionMarkII.cs)
122+
121123
[pade erfc precision16](https://github.com/tk-yoshimura/ErrorFunctionApproximation/blob/main/results/pade_erfc_e16.txt)
122-
[pade erfc precision32](https://github.com/tk-yoshimura/ErrorFunctionApproximation/blob/main/results/pade_erfc_e32.txt)
124+
[pade erfc precision32](https://github.com/tk-yoshimura/ErrorFunctionApproximation/blob/main/results/pade_erfc_e32.txt)
125+
[pade erf precision16](https://github.com/tk-yoshimura/ErrorFunctionApproximation/blob/main/results/pade_erf_e16.txt)
126+
[pade erf precision32](https://github.com/tk-yoshimura/ErrorFunctionApproximation/blob/main/results/pade_erf_e32.txt)

0 commit comments

Comments
 (0)