Skip to content

Commit c3a8ff2

Browse files
tk-yoshimuratk-yoshimura
committed
+ pade
1 parent f43da78 commit c3a8ff2

5 files changed

Lines changed: 488 additions & 43 deletions

File tree

‎ErrorFunctionApproximation/Program.cs‎

Lines changed: 60 additions & 42 deletions
Original file line numberDiff line numberDiff line change
@@ -7,59 +7,77 @@
77
namespace ErrorFunctionApproximation {
88
class Program {
99
static void Main(string[] args) {
10-
using StreamWriter sw = new("../../../../results_disused/pade_erfc_n32.txt");
11-
sw.WriteLine("pade approximant exp(-x^2)/erfc(x)");
10+
MultiPrecision<Pow2.N32>[] xs = new MultiPrecision<Pow2.N32>[]{
11+
0.5, 1, 2, 4, 8, 16, 27.25
12+
};
1213

13-
sw.WriteLine("expected");
14-
List<MultiPrecision<Pow2.N32>> expecteds = new();
15-
for (MultiPrecision<Pow2.N32> x = 0.5; x <= 27.25; x += 1d / 64) {
16-
MultiPrecision<Pow2.N32> y = MultiPrecision<Pow2.N32>.Exp(-x * x) / MultiPrecision<Pow2.N32>.Erfc(x);
14+
MultiPrecision<Pow2.N32> dx = 1d / 1024;
1715

18-
expecteds.Add(y);
16+
using StreamWriter sw = new("../../../../results/pade_erfc_e32.txt");
17+
sw.WriteLine("pade approximant exp(-x^2)/erfc(x)");
1918

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

23-
sw.WriteLine("diffs x = 0.5");
24-
List<MultiPrecision<Pow2.N32>> diffs = new();
25-
for (int d = 0; d <= 128; d++) {
26-
MultiPrecision<Pow2.N32> dy = ErfcExpScaled<Pow2.N32>.Diff(d, 0.5) * MultiPrecision<Pow2.N32>.TaylorSequence[d];
27-
diffs.Add(dy);
22+
sw.WriteLine($"\nrange x in [{x0}, {x1}]");
2823

29-
sw.WriteLine($"{d},{dy}");
30-
}
31-
sw.Flush();
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);
28+
29+
expecteds.Add(y);
30+
31+
//sw.WriteLine($"{x},{y:e40}");
32+
}
33+
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);
39+
40+
//sw.WriteLine($"{d},{dy:e40}");
41+
}
42+
sw.Flush();
3243

33-
sw.WriteLine("pade results");
34-
for (int m = 4; m <= 64; m++) {
35-
for (int n = m / 2; n <= m * 2 && m + n < 128; n++) {
36-
(MultiPrecision<Pow2.N32>[] ms, MultiPrecision<Pow2.N32>[] ns) =
37-
PadeSolver<Pow2.N32>.Solve(diffs.Take(m + n + 1).ToArray(), m, n);
44+
sw.WriteLine("pade results");
45+
bool is_finished = false;
3846

39-
MultiPrecision<Pow2.N32> err = 0;
40-
for ((int i, MultiPrecision<Pow2.N32> x) = (0, 0.5); x <= 27.25; i++, x += 1d / 64) {
41-
MultiPrecision<Pow2.N32> expected = expecteds[i];
42-
MultiPrecision<Pow2.N32> actual = PadeSolver<Pow2.N32>.Approx(x - 0.5, ms, ns);
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);
4351

44-
err = MultiPrecision<Pow2.N32>.Max(err, MultiPrecision<Pow2.N32>.Abs(expected / actual - 1));
45-
}
46-
if (err < 1e-15) {
47-
sw.WriteLine($"\nm={m},n={n}");
48-
sw.WriteLine("ms");
49-
for (int i = 0; i < ms.Length; i++) {
50-
sw.WriteLine(ms[i]);
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);
56+
57+
err = MultiPrecision<Pow2.N32>.Max(err, MultiPrecision<Pow2.N32>.Abs(expected / actual - 1));
5158
}
52-
sw.WriteLine("ns");
53-
for (int i = 0; i < ns.Length; i++) {
54-
sw.WriteLine(ns[i]);
59+
60+
Console.WriteLine($"m={m},n={n}");
61+
Console.WriteLine($"{err:e20}");
62+
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();
76+
77+
is_finished = true;
78+
break;
5579
}
56-
sw.WriteLine("relative err");
57-
sw.WriteLine($"{err:e20}");
58-
sw.Flush();
5980
}
60-
61-
Console.WriteLine($"m={m},n={n}");
62-
Console.WriteLine($"{err:e20}");
6381
}
6482
}
6583

‎README.md‎

Lines changed: 8 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -81,6 +81,8 @@ Let N be the recursion times (Fn+4(z) to Fn(z)), and N until convergence is as f
8181
![erfc convergence bitwise](https://github.com/tk-yoshimura/ErrorFunctionApproximation/blob/main/figures/erfc_convergence_bitswise.svg)
8282

8383
## Double Precision (IEEE 754) Approx
84+
### Taylor Series
85+
8486
[C# code](https://github.com/tk-yoshimura/ErrorFunctionApproximation/blob/main/ErrorFunctionApproximation/ErrorFunction.cs)
8587
[erf result](https://github.com/tk-yoshimura/ErrorFunctionApproximation/blob/main/results/erf_approx.csv)
8688
[erfc result](https://github.com/tk-yoshimura/ErrorFunctionApproximation/blob/main/results/erfc_approx.csv)
@@ -112,4 +114,9 @@ A(z) taylor series (z0=3.5)
112114
![erfc x20](https://github.com/tk-yoshimura/ErrorFunctionApproximation/blob/main/figures/erfc_approx_x20.svg)
113115

114116
larger 27.25:
115-
![erfc tolarge](https://github.com/tk-yoshimura/ErrorFunctionApproximation/blob/main/figures/erfc_approx_tolarge.svg)
117+
![erfc tolarge](https://github.com/tk-yoshimura/ErrorFunctionApproximation/blob/main/figures/erfc_approx_tolarge.svg)
118+
119+
### Pade Approximant
120+
![pade definition](https://github.com/tk-yoshimura/ErrorFunctionApproximation/blob/main/figures/pade_definition.svg)
121+
[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)

0 commit comments

Comments
 (0)