Skip to content

Commit f43da78

Browse files
tk-yoshimuratk-yoshimura
committed
erfc pade
1 parent 516d527 commit f43da78

10 files changed

Lines changed: 231 additions & 83 deletions

‎ErrorFunctionApproximation/ErfConvergenseTester.cs‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -11,7 +11,7 @@ public static (MultiPrecision<N> y, int convergense_n) Erf(MultiPrecision<N> z)
1111
MultiPrecision<Plus4<N>> s = MultiPrecision<Plus4<N>>.One;
1212

1313
MultiPrecision<N> s_prev = null;
14-
14+
1515
int n = 1;
1616
while (s_prev != s.Convert<N>()) {
1717
s_prev = s.Convert<N>();

‎ErrorFunctionApproximation/ErfTaylorSeries.cs‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,5 @@
1-
using System;
2-
using MultiPrecision;
1+
using MultiPrecision;
2+
using System;
33
using System.Collections.Generic;
44
using System.Numerics;
55

@@ -22,7 +22,7 @@ public static MultiPrecision<N> Coef<N>(int n) where N : struct, IConstant {
2222
return coefs[n];
2323
}
2424

25-
for(int k = factorials.Count; k <= n; k++){
25+
for (int k = factorials.Count; k <= n; k++) {
2626
BigInteger fact = factorials[k - 1] * k;
2727
BigInteger coef = fact * checked(2 * k + 1);
2828

‎ErrorFunctionApproximation/ErfcContinuedFraction.cs‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -14,7 +14,7 @@ public static MultiPrecision<N> Frac(MultiPrecision<N> z, long n) {
1414
MultiPrecision<Plus4<N>> z_ex = z.Convert<Plus4<N>>();
1515
MultiPrecision<Plus4<N>> w = z_ex * z_ex;
1616

17-
MultiPrecision<Plus4<N>> f =
17+
MultiPrecision<Plus4<N>> f =
1818
(MultiPrecision<Plus4<N>>.Sqrt(25 + w * (440 + w * (488 + w * 16 * (10 + w))))
1919
- 5 + w * 4 * (1 + w))
2020
/ (20 + w * 8);

‎ErrorFunctionApproximation/ErfcExpScaled.cs‎

Lines changed: 102 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,102 @@
1+
using MultiPrecision;
2+
using System;
3+
using System.Collections.Generic;
4+
using System.Collections.ObjectModel;
5+
using System.Linq;
6+
7+
namespace ErrorFunctionApproximation {
8+
public static class ErfcExpScaled<N> where N : struct, IConstant {
9+
private static readonly Dictionary<int, ReadOnlyCollection<(MultiPrecision<Double<N>> c, int n, int m)>> table = new() {
10+
{ 0, new ReadOnlyCollection<(MultiPrecision<Double<N>> c, int n, int m)>(new (MultiPrecision<Double<N>> c, int n, int m)[]{ (1, 0, 1) })},
11+
};
12+
13+
public static MultiPrecision<N> Diff(int d, MultiPrecision<N> x) {
14+
MultiPrecision<Double<N>> x_ex = x.Convert<Double<N>>();
15+
MultiPrecision<Double<N>> y = MultiPrecision<Double<N>>.Exp(-x_ex * x_ex) / MultiPrecision<Double<N>>.Erfc(x_ex);
16+
17+
MultiPrecision<Double<N>> s = 0, s0, s1;
18+
19+
ReadOnlyCollection<(MultiPrecision<Double<N>> c, int n, int m)> coef = Coef(d);
20+
21+
int i = 0;
22+
for (i = 0; i < coef.Count - 1; i += 2) {
23+
(MultiPrecision<Double<N>> c0, int n0, int m0) = coef[i];
24+
(MultiPrecision<Double<N>> c1, int n1, int m1) = coef[i + 1];
25+
26+
(int n, int m) = (Math.Min(n0, n1), Math.Min(m0, m1));
27+
28+
if (n0 > n) {
29+
s0 = c0 * MultiPrecision<Double<N>>.Pow(x_ex, n0 - n) * MultiPrecision<Double<N>>.Pow(y, m0 - m);
30+
}
31+
else {
32+
s0 = c0 * MultiPrecision<Double<N>>.Pow(y, m0 - m);
33+
}
34+
35+
if (n1 > n) {
36+
s1 = c1 * MultiPrecision<Double<N>>.Pow(x_ex, n1 - n) * MultiPrecision<Double<N>>.Pow(y, m1 - m);
37+
}
38+
else {
39+
s1 = c1 * MultiPrecision<Double<N>>.Pow(y, m1 - m);
40+
}
41+
42+
s += MultiPrecision<Double<N>>.Pow(x_ex, n) * MultiPrecision<Double<N>>.Pow(y, m) * (s0 + s1);
43+
}
44+
45+
for (; i < coef.Count; i++) {
46+
(MultiPrecision<Double<N>> c, int n, int m) = coef[i];
47+
48+
if (n > 0) {
49+
s += c * MultiPrecision<Double<N>>.Pow(x_ex, n) * MultiPrecision<Double<N>>.Pow(y, m);
50+
}
51+
else {
52+
s += c * MultiPrecision<Double<N>>.Pow(y, m);
53+
}
54+
}
55+
56+
return s.Convert<N>();
57+
}
58+
59+
private static ReadOnlyCollection<(MultiPrecision<Double<N>> c, int n, int m)> Coef(int d) {
60+
if (d < 0) {
61+
throw new ArgumentOutOfRangeException(nameof(d));
62+
}
63+
64+
if (table.ContainsKey(d)) {
65+
return table[d];
66+
}
67+
68+
MultiPrecision<Double<N>> inv_sqrtpi = 1 / MultiPrecision<Double<N>>.Sqrt(MultiPrecision<Double<N>>.PI);
69+
70+
Dictionary<(int n, int m), MultiPrecision<Double<N>>> t = new();
71+
72+
foreach ((MultiPrecision<Double<N>> c, int n, int m) in Coef(d - 1)) {
73+
if (!t.ContainsKey((n, m + 1))) {
74+
t.Add((n, m + 1), 0);
75+
}
76+
t[(n, m + 1)] += 2 * m * inv_sqrtpi * c;
77+
78+
if (!t.ContainsKey((n + 1, m))) {
79+
t.Add((n + 1, m), 0);
80+
}
81+
t[(n + 1, m)] -= 2 * m * c;
82+
83+
if (n > 0) {
84+
if (!t.ContainsKey((n - 1, m))) {
85+
t.Add((n - 1, m), 0);
86+
}
87+
t[(n - 1, m)] += n * c;
88+
}
89+
}
90+
91+
ReadOnlyCollection<(MultiPrecision<Double<N>> c, int n, int m)> coef = new(
92+
t.
93+
OrderByDescending((item) => item.Key.m).ThenByDescending((item) => item.Key.n).
94+
Select((item) => (item.Value, item.Key.n, item.Key.m)).ToArray()
95+
);
96+
97+
table.Add(d, coef);
98+
99+
return coef;
100+
}
101+
}
102+
}

‎ErrorFunctionApproximation/ErrorFunction.cs‎

Lines changed: 7 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -141,11 +141,11 @@ public static double InverseErf(double z) {
141141
if (z == 1) {
142142
return double.PositiveInfinity;
143143
}
144-
if (z < 0) {
144+
if (z < 0) {
145145
return -InverseErf(Abs(z));
146146
}
147147

148-
if (z < 1.220703125e-4) {
148+
if (z < 1.220703125e-4) {
149149
double w = PI * z * z;
150150
double t = Sqrt(PI) * ((40320 + w * (3360 + w * (588 + w * 127))) / 80640);
151151

@@ -171,7 +171,7 @@ public static double InverseErfc(double z) {
171171
if (z == 2) {
172172
return double.NegativeInfinity;
173173
}
174-
if (z >= 0.5) {
174+
if (z >= 0.5) {
175175
return InverseErf(1 - z);
176176
}
177177

@@ -256,14 +256,14 @@ private static double ErfcLarger4(double z) {
256256
return v * f;
257257
}
258258

259-
private static double InverseErfRootFinding(double x) {
259+
private static double InverseErfRootFinding(double x) {
260260
const double a = 0.147;
261261

262262
double s = 2 * c;
263263
double lg = Log(1 - x * x), lga = 2 / (PI * a) + lg / 2;
264264
double z = Sqrt(Sqrt(lga * lga - lg / a) - lga);
265265

266-
for (int i = 0; i < 8; i++) {
266+
for (int i = 0; i < 8; i++) {
267267
double y = Erf(z) - x;
268268
double df1 = Exp(-z * z) * s;
269269
double df2 = -2 * z * df1;
@@ -279,14 +279,14 @@ private static double InverseErfRootFinding(double x) {
279279
return z;
280280
}
281281

282-
private static double InverseErfcRootFinding(double x) {
282+
private static double InverseErfcRootFinding(double x) {
283283
const double a = 0.147;
284284

285285
double s = 2 * c;
286286
double lg = Log((2 - x) * x), lga = 2 / (PI * a) + lg / 2;
287287
double z = Sqrt(Sqrt(lga * lga - lg / a) - lga);
288288

289-
for (int i = 0; i < 8; i++) {
289+
for (int i = 0; i < 8; i++) {
290290
double y = Erfc(z) - x;
291291
double df1 = -Exp(-z * z) * s;
292292
double df2 = -2 * z * df1;

‎ErrorFunctionApproximation/ErrorFunctionApproximation.csproj‎

Lines changed: 3 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -2,13 +2,12 @@
22

33
<PropertyGroup>
44
<OutputType>Exe</OutputType>
5-
<TargetFramework>net5.0</TargetFramework>
5+
<TargetFramework>net6.0</TargetFramework>
66
</PropertyGroup>
77

88
<ItemGroup>
9-
<Reference Include="MultiPrecision">
10-
<HintPath>..\import\MultiPrecision.dll</HintPath>
11-
</Reference>
9+
<PackageReference Include="TYoshimura.MultiPrecision" Version="5.1.0" />
10+
<PackageReference Include="TYoshimura.MultiPrecision.Algebra" Version="1.1.0" />
1211
</ItemGroup>
1312

1413
</Project>

‎ErrorFunctionApproximation/Expand.cs‎

Lines changed: 5 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
using MultiPrecision;
22

3-
namespace ErrorFunctionApproximation {
3+
namespace ErrorFunctionApproximation {
44
internal struct Plus4<N> : IConstant where N : struct, IConstant {
55
public int Value => checked(default(N).Value + 4);
66
}
@@ -16,4 +16,8 @@ internal struct Expand50<N> : IConstant where N : struct, IConstant {
1616
internal struct Expand75<N> : IConstant where N : struct, IConstant {
1717
public int Value => checked(default(N).Value * 7 / 4);
1818
}
19+
20+
internal struct Double<N> : IConstant where N : struct, IConstant {
21+
public int Value => checked(default(N).Value * 2);
22+
}
1923
}

‎ErrorFunctionApproximation/PadeSolver.cs‎

Lines changed: 56 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,56 @@
1+
using MultiPrecision;
2+
using MultiPrecisionAlgebra;
3+
using System;
4+
using System.Linq;
5+
6+
namespace ErrorFunctionApproximation {
7+
internal static class PadeSolver<N> where N : struct, IConstant {
8+
public static (MultiPrecision<N>[] ms, MultiPrecision<N>[] ns) Solve(MultiPrecision<N>[] cs, int m, int n) {
9+
if (m < 0) {
10+
throw new ArgumentOutOfRangeException(nameof(m));
11+
}
12+
if (n < 0) {
13+
throw new ArgumentOutOfRangeException(nameof(n));
14+
}
15+
if (cs.Length != checked(m + n + 1)) {
16+
throw new ArgumentException("Illegal length.", nameof(cs));
17+
}
18+
19+
int k = m + n;
20+
21+
Matrix<N> a = new(k, k);
22+
Vector<N> c = cs[1..];
23+
24+
for (int i = 0; i < m; i++) {
25+
a[i, i] = 1;
26+
}
27+
28+
for (int i = m; i < k; i++) {
29+
for (int j = i - m, r = 0; j < k; j++, r++) {
30+
a[j, i] = -cs[r];
31+
}
32+
}
33+
34+
Vector<N> v = a.Inverse * c;
35+
MultiPrecision<N>[] ms = new MultiPrecision<N>[] { cs[0] }.Concat(((MultiPrecision<N>[])v)[..m]).ToArray();
36+
MultiPrecision<N>[] ns = new MultiPrecision<N>[] { 1 }.Concat(((MultiPrecision<N>[])v)[m..]).ToArray();
37+
38+
return (ms, ns);
39+
}
40+
41+
public static MultiPrecision<N> Approx(MultiPrecision<N> a, MultiPrecision<N>[] ms, MultiPrecision<N>[] ns) {
42+
MultiPrecision<N> p = ms[^1], q = ns[^1];
43+
44+
for (int i = ms.Length - 2; i >= 0; i--) {
45+
p = p * a + ms[i];
46+
}
47+
for (int i = ns.Length - 2; i >= 0; i--) {
48+
q = q * a + ns[i];
49+
}
50+
51+
MultiPrecision<N> y = p / q;
52+
53+
return y;
54+
}
55+
}
56+
}

0 commit comments

Comments
 (0)