Skip to content

Commit a153100

Browse files
tk-yoshimuratk-yoshimura
committed
+ series inverf
1 parent e762667 commit a153100

6 files changed

Lines changed: 563 additions & 88 deletions

File tree

‎ErrorFunctionApproximation/InvErfTaylorSeries.cs‎

Lines changed: 30 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,30 @@
1+
using MultiPrecision;
2+
using System;
3+
using System.Collections.Generic;
4+
5+
namespace ErrorFunctionApproximation {
6+
public static class InvErfTaylorSeries {
7+
private static readonly List<Fraction> coefs = new() {
8+
1, 1
9+
};
10+
11+
public static MultiPrecision<N> Coef<N>(int n) where N : struct, IConstant {
12+
if (n < 0) {
13+
throw new ArgumentOutOfRangeException(nameof(n));
14+
}
15+
16+
for (int k = coefs.Count; k <= n; k++) {
17+
Fraction coef = 0;
18+
19+
for (int m = 0; m < k; m++) {
20+
coef += coefs[m] * coefs[k - m - 1] / ((m + 1) * (2 * m + 1));
21+
}
22+
23+
coefs.Add(coef);
24+
}
25+
26+
return MultiPrecision<N>.Div(coefs[n].Numer, coefs[n].Denom) / (2 * n + 1)
27+
* MultiPrecision<N>.Pow(MultiPrecision<N>.Sqrt(MultiPrecision<N>.PI) / 2, 2 * n + 1);
28+
}
29+
}
30+
}

‎ErrorFunctionApproximation/Program.cs‎

Lines changed: 59 additions & 86 deletions
Original file line numberDiff line numberDiff line change
@@ -7,101 +7,74 @@
77
namespace ErrorFunctionApproximation {
88
class Program {
99
static void Main(string[] args) {
10-
MultiPrecision<Pow2.N32>[] xs = new MultiPrecision<Pow2.N32>[]{
11-
1, 2, 4, 8, 16, 32
12-
};
10+
MultiPrecision<Pow2.N32> dx = 1d / 1024;
1311

14-
static MultiPrecision<Pow2.N32> f32(MultiPrecision<Pow2.N32> x) {
15-
return MultiPrecision<Pow2.N32>.InverseErfc(MultiPrecision<Pow2.N32>.Pow2(-x * x));
16-
};
12+
using StreamWriter sw = new("../../../../results/series_inverf_e20.txt");
13+
sw.WriteLine("pade approximant f(x^2) f(x):=inverf(x)/x");
1714

18-
static MultiPrecision<Pow2.N64> f64(MultiPrecision<Pow2.N64> x) {
19-
return MultiPrecision<Pow2.N64>.InverseErfc(MultiPrecision<Pow2.N64>.Pow2(-x * x));
20-
};
15+
MultiPrecision<Pow2.N32> x0 = 0, x1 = 0.5;
2116

22-
using StreamWriter sw = new("../../../../results/pade_inverfc_e16.txt");
23-
sw.WriteLine("pade approximant inverse_erfc(2^(-x^2))");
17+
sw.WriteLine($"\nrange x in [{x0}, {x1}]");
2418

25-
for (int j = 0; j < xs.Length - 1; j++) {
26-
MultiPrecision<Pow2.N32> x0 = xs[j], x1 = xs[j + 1];
27-
MultiPrecision<Pow2.N32> dx = (x1 - x0) / 2048;
19+
MultiPrecision<Pow2.N32> c = MultiPrecision<Pow2.N32>.Sqrt(MultiPrecision<Pow2.N32>.PI) / 2;
2820

29-
sw.WriteLine($"\nrange x in [{x0}, {x1}]");
30-
31-
//sw.WriteLine("expected");
32-
List<MultiPrecision<Pow2.N32>> expecteds = new();
33-
for (MultiPrecision<Pow2.N32> x = x0; x <= x1; x += dx) {
34-
MultiPrecision<Pow2.N32> y = f32(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>.InverseErf(x) / x : c;
3525

36-
expecteds.Add(y);
37-
38-
if ((x % 0.125) == 0) {
39-
Console.Write(".");
40-
}
26+
expecteds.Add(y);
4127

42-
//sw.WriteLine($"{x},{y:e40}");
43-
}
44-
45-
Console.Write("\n");
28+
//sw.WriteLine($"{x},{y:e40}");
29+
}
4630

47-
sw.WriteLine($"diffs x = {x0}");
48-
MultiPrecision<Pow2.N64>[] diffs = FiniteDifference<Pow2.N64>.Diff(x0.Convert<Pow2.N64>(), f64, Math.ScaleB(1, -20));
49-
for (int i = 0; i < diffs.Length; i++) {
50-
sw.WriteLine($"{i},{diffs[i]:e40}");
51-
}
52-
sw.Flush();
53-
54-
MultiPrecision<Pow2.N32>[] cs = new MultiPrecision<Pow2.N32>[diffs.Length + 1];
55-
cs[0] = f32(x0);
56-
for (int i = 0; i < diffs.Length; i++) {
57-
cs[i + 1] = diffs[i].Convert<Pow2.N32>() * MultiPrecision<Pow2.N32>.TaylorSequence[i + 1];
58-
}
59-
60-
sw.WriteLine("pade results");
61-
bool is_finished = false;
62-
63-
for (int m = 4; m < diffs.Length && !is_finished; m++) {
64-
for (int n = m - 1; n <= m + 1 && m + n < diffs.Length; n++) {
65-
(MultiPrecision<Pow2.N32>[] ms, MultiPrecision<Pow2.N32>[] ns) =
66-
PadeSolver<Pow2.N32>.Solve(cs.Take(m + n + 1).ToArray(), m, n);
67-
68-
MultiPrecision<Pow2.N32> err = 0;
69-
for ((int i, MultiPrecision<Pow2.N32> x) = (0, x0); i < expecteds.Count; i++, x += dx) {
70-
MultiPrecision<Pow2.N32> expected = expecteds[i];
71-
MultiPrecision<Pow2.N32> actual = PadeSolver<Pow2.N32>.Approx(x - x0, ms, ns);
72-
73-
err = MultiPrecision<Pow2.N32>.Max(err, MultiPrecision<Pow2.N32>.Abs(expected / actual - 1));
74-
}
75-
76-
Console.WriteLine($"m={m},n={n}");
77-
Console.WriteLine($"{err:e20}");
78-
79-
if (err < 1e-16) {
80-
sw.WriteLine($"m={m},n={n}");
81-
sw.WriteLine("ms");
82-
for (int i = 0; i < ms.Length; i++) {
83-
sw.WriteLine($"{ms[i]:e20}");
84-
}
85-
sw.WriteLine("ns");
86-
for (int i = 0; i < ns.Length; i++) {
87-
sw.WriteLine($"{ns[i]:e20}");
88-
}
89-
sw.WriteLine("relative err");
90-
sw.WriteLine($"{err:e20}");
91-
sw.Flush();
92-
93-
if (ms.All((v) => v > 0) && ns.All((v) => v > 0)) {
94-
is_finished = true;
95-
}
96-
}
97-
}
98-
}
99-
100-
if (!is_finished) {
101-
Console.WriteLine("not convergence");
102-
}
31+
sw.WriteLine($"series x = {x0}");
32+
List<MultiPrecision<Pow2.N32>> diffs = new();
33+
for (int d = 0; d <= 64; d++) {
34+
MultiPrecision<Pow2.N32> dy = InvErfTaylorSeries.Coef<Pow2.N32>(d);
35+
diffs.Add(dy);
36+
37+
sw.WriteLine($"{d},{dy:e20}");
10338
}
104-
39+
sw.Flush();
40+
41+
//sw.WriteLine("pade results");
42+
//bool is_finished = false;
43+
//
44+
//for (int m = 4; m <= 64 && !is_finished; m++) {
45+
// for (int n = m - 1; n <= m + 1 && m + n < 128; n++) {
46+
// (MultiPrecision<Pow2.N32>[] ms, MultiPrecision<Pow2.N32>[] ns) =
47+
// PadeSolver<Pow2.N32>.Solve(diffs.Take(m + n + 1).ToArray(), m, n);
48+
//
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);
53+
//
54+
// err = MultiPrecision<Pow2.N32>.Max(err, MultiPrecision<Pow2.N32>.Abs(expected / actual - 1));
55+
// }
56+
//
57+
// Console.WriteLine($"m={m},n={n}");
58+
// Console.WriteLine($"{err:e20}");
59+
//
60+
// if (err < 1e-32) {
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}");
65+
// }
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+
// }
76+
// }
77+
//}
10578

10679
Console.WriteLine("END");
10780
Console.Read();

‎README.md‎

Lines changed: 8 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -127,8 +127,14 @@ larger 27.25:
127127

128128
## Inverse Function
129129

130-
### Erfc NearZero
130+
### InverseErfc NearZero
131131
![inverse erfc nearzero](https://github.com/tk-yoshimura/ErrorFunctionApproximation/blob/main/figures/inverse_erfc_nz.svg)
132132

133133
[pade inverfc precision16](https://github.com/tk-yoshimura/ErrorFunctionApproximation/blob/main/results/pade_inverfc_e16.txt)
134-
[pade inverfc precision32](https://github.com/tk-yoshimura/ErrorFunctionApproximation/blob/main/results/pade_inverfc_e32.txt)
134+
[pade inverfc precision32](https://github.com/tk-yoshimura/ErrorFunctionApproximation/blob/main/results/pade_inverfc_e32.txt)
135+
136+
### InverseErf NearZero
137+
138+
[series inverf e20](https://github.com/tk-yoshimura/ErrorFunctionApproximation/blob/main/results/series_inverf_e20.txt)
139+
[series inverf e40](https://github.com/tk-yoshimura/ErrorFunctionApproximation/blob/main/results/series_inverf_e40.txt)
140+
[series inverf e80](https://github.com/tk-yoshimura/ErrorFunctionApproximation/blob/main/results/series_inverf_e80.txt)

‎results/series_inverf_e20.txt‎

Lines changed: 70 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,70 @@
1+
pade approximant f(x^2) f(x):=inverf(x)/x
2+
3+
range x in [0, 0.5]
4+
expected
5+
series x = 0
6+
0,8.86226925452758013649e-1
7+
1,2.32013666534654493554e-1
8+
2,1.27556175305597958254e-1
9+
3,8.65521292415475337296e-2
10+
4,6.49596177453854133820e-2
11+
5,5.17312819846163741126e-2
12+
6,4.28367206517973498447e-2
13+
7,3.64659293085316263256e-2
14+
8,3.16890050216054468096e-2
15+
9,2.79806329649952247334e-2
16+
10,2.50222758411983494572e-2
17+
11,2.26098633188975744328e-2
18+
12,2.06067803790590017188e-2
19+
13,1.89182172507788544635e-2
20+
14,1.74763705628565461904e-2
21+
15,1.62315009876852512753e-2
22+
16,1.51463150632478055204e-2
23+
17,1.41923160025099641512e-2
24+
18,1.33473641974212971503e-2
25+
19,1.25940048713320698416e-2
26+
20,1.19182959363920398739e-2
27+
21,1.13089701059225367722e-2
28+
22,1.07568253033179575778e-2
29+
23,1.02542740818534682141e-2
30+
24,9.79500577007117565181e-3
31+
25,9.37372981918208151683e-3
32+
26,8.98597850284337808804e-3
33+
27,8.62795358070943465254e-3
34+
28,8.29640592773924140441e-3
35+
29,7.98854016260335483777e-3
36+
30,7.70193843225976116637e-3
37+
31,7.43449901783153060384e-3
38+
32,7.18438651126831244412e-3
39+
33,6.94999110106471485712e-3
40+
34,6.72989508534152981828e-3
41+
35,6.52284516145006176435e-3
42+
36,6.32772936434314323002e-3
43+
37,6.14355777038415591689e-3
44+
38,5.96944626973445629437e-3
45+
39,5.80460285383430176030e-3
46+
40,5.64831597555156001751e-3
47+
41,5.49994462620396350840e-3
48+
42,5.35890984168644547029e-3
49+
43,5.22468740368749060984e-3
50+
44,5.09680154470620556637e-3
51+
45,4.97481949973904039883e-3
52+
46,4.85834677495848553682e-3
53+
47,4.74702302588497662306e-3
54+
48,4.64051845555903192536e-3
55+
49,4.53853065790708485618e-3
56+
50,4.44078184352718353227e-3
57+
51,4.34701639502151030530e-3
58+
52,4.25699870718272882770e-3
59+
53,4.17051127412617157112e-3
60+
54,4.08735299110902163849e-3
61+
55,4.00733764349814205679e-3
62+
56,3.93029255930647814844e-3
63+
57,3.85605740504821755180e-3
64+
58,3.78448310747382581560e-3
65+
59,3.71543088612604792433e-3
66+
60,3.64877138367909283721e-3
67+
61,3.58438388274456894823e-3
68+
62,3.52215559929786878265e-3
69+
63,3.46198104413766023086e-3
70+
64,3.40376144487207088729e-3

0 commit comments

Comments
 (0)