Skip to content

Commit e762667

Browse files
tk-yoshimuratk-yoshimura
committed
+ pade inverfc
1 parent 75de991 commit e762667

8 files changed

Lines changed: 2343 additions & 58 deletions

File tree

‎ErrorFunctionApproximation/FiniteDifference.cs‎

Lines changed: 130 additions & 0 deletions
Large diffs are not rendered by default.

‎ErrorFunctionApproximation/Fraction.cs‎

Lines changed: 118 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,118 @@
1+
using System;
2+
using System.Numerics;
3+
4+
namespace ErrorFunctionApproximation {
5+
internal class Fraction {
6+
public BigInteger Numer { private set; get; }
7+
public BigInteger Denom { private set; get; }
8+
9+
public Fraction(long n) : this(new BigInteger(n)) { }
10+
11+
public Fraction(BigInteger n) {
12+
this.Numer = n;
13+
this.Denom = 1;
14+
}
15+
16+
public Fraction(long numer, long denom) : this(new BigInteger(numer), new BigInteger(denom)) { }
17+
18+
public Fraction(BigInteger numer, BigInteger denom) {
19+
this.Numer = numer;
20+
this.Denom = denom;
21+
22+
Reduce();
23+
}
24+
25+
public static implicit operator Fraction(long n) {
26+
return new Fraction(n);
27+
}
28+
29+
public static implicit operator Fraction(BigInteger n) {
30+
return new Fraction(n);
31+
}
32+
33+
public static implicit operator Fraction(string str) {
34+
if (!str.Contains('/')) {
35+
return BigInteger.Parse(str);
36+
}
37+
38+
string[] str_item = str.Split('/');
39+
if (str_item.Length != 2) {
40+
throw new FormatException();
41+
}
42+
43+
return new(BigInteger.Parse(str_item[0]), BigInteger.Parse(str_item[1]));
44+
}
45+
46+
public static Fraction operator +(Fraction v1, Fraction v2) {
47+
return new Fraction(v1.Numer * v2.Denom + v2.Numer * v1.Denom, v1.Denom * v2.Denom);
48+
}
49+
50+
public static Fraction operator -(Fraction v1, Fraction v2) {
51+
return new Fraction(v1.Numer * v2.Denom - v2.Numer * v1.Denom, v1.Denom * v2.Denom);
52+
}
53+
54+
public static Fraction operator *(Fraction v1, Fraction v2) {
55+
return new Fraction(v1.Numer * v2.Numer, v1.Denom * v2.Denom);
56+
}
57+
58+
public static Fraction operator *(Fraction v, BigInteger n) {
59+
return new Fraction(v.Numer * n, v.Denom);
60+
}
61+
62+
public static Fraction operator /(Fraction v1, Fraction v2) {
63+
return new Fraction(v1.Numer * v2.Denom, v1.Denom * v2.Numer);
64+
}
65+
66+
public static Fraction operator /(Fraction v, BigInteger n) {
67+
return new Fraction(v.Numer, v.Denom * n);
68+
}
69+
70+
public static Fraction Abs(Fraction v) {
71+
return new Fraction(BigInteger.Abs(v.Numer), v.Denom);
72+
}
73+
74+
public void Reduce() {
75+
if (Denom == 0) {
76+
throw new DivideByZeroException();
77+
}
78+
79+
if (Numer == 0) {
80+
Denom = 1;
81+
return;
82+
}
83+
84+
int sign = ((Numer < 0) ^ (Denom < 0)) ? -1 : +1;
85+
BigInteger n = Numer, d = Denom;
86+
87+
d = BigInteger.Abs(d);
88+
n = BigInteger.Abs(n);
89+
90+
if (n > d) {
91+
BigInteger temp = n;
92+
n = d;
93+
d = temp;
94+
}
95+
96+
BigInteger rem;
97+
while ((rem = d % n) != 0) {
98+
d = n;
99+
n = rem;
100+
}
101+
102+
Numer = (sign > 0 ? BigInteger.Abs(Numer) : -BigInteger.Abs(Numer)) / n;
103+
Denom = BigInteger.Abs(Denom) / n;
104+
}
105+
106+
public override string ToString() {
107+
return Denom != 1 ? $"{Numer}/{Denom}" : $"{Numer}";
108+
}
109+
110+
public override int GetHashCode() {
111+
return Numer.GetHashCode() ^ Denom.GetHashCode();
112+
}
113+
114+
public override bool Equals(object obj) {
115+
return obj is Fraction v && Numer == v.Numer && Denom == v.Denom;
116+
}
117+
}
118+
}

‎ErrorFunctionApproximation/Program.cs‎

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

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

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

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

19-
MultiPrecision<Pow2.N32> c = 2 / MultiPrecision<Pow2.N32>.Sqrt(MultiPrecision<Pow2.N32>.PI);
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;
2028

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;
25-
26-
expecteds.Add(y);
27-
28-
//sw.WriteLine($"{x},{y:e40}");
29-
}
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);
3035

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);
36+
expecteds.Add(y);
37+
38+
if ((x % 0.125) == 0) {
39+
Console.Write(".");
40+
}
3641

37-
//sw.WriteLine($"{d},{dy:e40}");
38-
}
39-
sw.Flush();
42+
//sw.WriteLine($"{x},{y:e40}");
43+
}
4044

41-
sw.WriteLine("pade results");
42-
bool is_finished = false;
45+
Console.Write("\n");
46+
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();
4353

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);
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+
}
4859

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);
60+
sw.WriteLine("pade results");
61+
bool is_finished = false;
5362

54-
err = MultiPrecision<Pow2.N32>.Max(err, MultiPrecision<Pow2.N32>.Abs(expected / actual - 1));
55-
}
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);
5667

57-
Console.WriteLine($"m={m},n={n}");
58-
Console.WriteLine($"{err:e20}");
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);
5972

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}");
65-
}
66-
sw.WriteLine("ns");
67-
for (int i = 0; i < ns.Length; i++) {
68-
sw.WriteLine($"{ns[i]:e40}");
73+
err = MultiPrecision<Pow2.N32>.Max(err, MultiPrecision<Pow2.N32>.Abs(expected / actual - 1));
6974
}
70-
sw.WriteLine("relative err");
71-
sw.WriteLine($"{err:e20}");
72-
sw.Flush();
7375

74-
is_finished = true;
75-
break;
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+
}
7697
}
7798
}
99+
100+
if (!is_finished) {
101+
Console.WriteLine("not convergence");
102+
}
78103
}
79104

105+
80106
Console.WriteLine("END");
81107
Console.Read();
82108
}

‎README.md‎

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -128,4 +128,7 @@ larger 27.25:
128128
## Inverse Function
129129

130130
### Erfc NearZero
131-
![inverse erfc nearzero](https://github.com/tk-yoshimura/ErrorFunctionApproximation/blob/main/figures/inverse_erfc_nz.svg)
131+
![inverse erfc nearzero](https://github.com/tk-yoshimura/ErrorFunctionApproximation/blob/main/figures/inverse_erfc_nz.svg)
132+
133+
[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)

‎results/pade_erf_e16.txt‎

Lines changed: 0 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1,9 +1,6 @@
11
pade approximant f(x^2) f(x):=erf(x)/x
22

33
range x in [0, 0.5]
4-
expected
5-
diffs x = 0
6-
pade results
74
m=4,n=5
85
ms
96
1.12837916709551257390e0

‎results/pade_erf_e32.txt‎

Lines changed: 0 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1,9 +1,6 @@
11
pade approximant f(x^2) f(x):=erf(x)/x
22

33
range x in [0, 0.5]
4-
expected
5-
diffs x = 0
6-
pade results
74
m=8,n=9
85
ms
96
1.1283791670955125738961589031215451716881e0

0 commit comments

Comments
 (0)