是否可以在 C 中计算反误差函数?
我可以找到 erf(x)
在 <math.h>
它计算误差函数,但我找不到任何东西来做反函数。
最佳答案
目前,ISO C 标准数学库不包括 erfinv()
或其单精度变体 erfinvf()
。然而,创建自己的版本并不太难,我在下面通过具有合理准确性和性能的 erfinvf()
实现来演示。
查看 graph of the inverse error function我们观察到它是高度非线性的,因此很难用多项式近似。处理这种情况的一种策略是通过将更简单的初等函数(它们本身可以以高性能和出色的精度计算)和更容易服从多项式逼近或有理数的相当线性的函数组合来“线性化”这样的函数低阶近似值。
这里是文献中已知的一些erfinv
线性化方法,所有这些方法都基于对数。通常,作者会区分逆误差函数的一个主要的、相当线性的部分,从零到大约 0.9 的切换点,以及从切换点到 1 的 tail 部分。下面log()表示自然对数,R()表示有理逼近,P()表示多项式逼近。
一个。 J. Strecok,“关于误差函数的反函数的计算”。 计算数学,卷。 22,第 101 期(1968 年 1 月),第 144-158 页 (online)
β(x) = (-log(1-x2]))½; erfinv(x) = x · R(x2) [main]; R(x) · β(x) [tail]
J. M. Blair、C. A. Edwards、J. H. Johnson,“误差函数反函数的有理切比雪夫近似”。 计算数学,卷。 30,第 136 期(1976 年 10 月),第 827-830 页 (online)
ξ = (-log(1-x))-½; erfinv(x) = x · R(x2) [main]; ξ-1 · R(ξ) [tail]
米。贾尔斯,“逼近 erfinv 函数。” GPU Computing Gems Jade Edition,第 109-116 页。 2011. (online)
w = -log(1-x2); s = √w; erfinv(x) = x · P(w) [main]; x · P(s) [tail]
下面的解决方案通常遵循 Giles 的方法,但对其进行了简化,不需要尾部的平方根,即它使用 x · P(w) 类型的两个近似值。该代码最大限度地利用了融合乘加运算 FMA,它通过 C 中的标准数学函数 fma()
和 fmaf()
公开。许多常见的计算平台, 如
IBM Power、Arm64、x86-64 和 GPU 在硬件中提供此操作。在不存在硬件支持的情况下,使用 fma{f}()
可能会使下面的代码慢得无法接受,因为该操作需要由标准数学库模拟。此外,功能不正确的 FMA 仿真 are known to exist .
标准数学库的对数函数logf()
的精度会对下面my_erfinvf()
的精度产生一定的影响。只要该库提供了错误 < 1 ulp 的忠实全面的实现,规定的错误范围就应该成立,我尝试过的几个库也是如此。为了提高可重现性,我包含了我自己的可移植的忠实实现,my_logf()
。
#include <math.h>
float my_logf (float);
/* compute inverse error functions with maximum error of 2.35793 ulp */
float my_erfinvf (float a)
{
float p, r, t;
t = fmaf (a, 0.0f - a, 1.0f);
t = my_logf (t);
if (fabsf(t) > 6.125f) { // maximum ulp error = 2.35793
p = 3.03697567e-10f; // 0x1.4deb44p-32
p = fmaf (p, t, 2.93243101e-8f); // 0x1.f7c9aep-26
p = fmaf (p, t, 1.22150334e-6f); // 0x1.47e512p-20
p = fmaf (p, t, 2.84108955e-5f); // 0x1.dca7dep-16
p = fmaf (p, t, 3.93552968e-4f); // 0x1.9cab92p-12
p = fmaf (p, t, 3.02698812e-3f); // 0x1.8cc0dep-9
p = fmaf (p, t, 4.83185798e-3f); // 0x1.3ca920p-8
p = fmaf (p, t, -2.64646143e-1f); // -0x1.0eff66p-2
p = fmaf (p, t, 8.40016484e-1f); // 0x1.ae16a4p-1
} else { // maximum ulp error = 2.35002
p = 5.43877832e-9f; // 0x1.75c000p-28
p = fmaf (p, t, 1.43285448e-7f); // 0x1.33b402p-23
p = fmaf (p, t, 1.22774793e-6f); // 0x1.499232p-20
p = fmaf (p, t, 1.12963626e-7f); // 0x1.e52cd2p-24
p = fmaf (p, t, -5.61530760e-5f); // -0x1.d70bd0p-15
p = fmaf (p, t, -1.47697632e-4f); // -0x1.35be90p-13
p = fmaf (p, t, 2.31468678e-3f); // 0x1.2f6400p-9
p = fmaf (p, t, 1.15392581e-2f); // 0x1.7a1e50p-7
p = fmaf (p, t, -2.32015476e-1f); // -0x1.db2aeep-3
p = fmaf (p, t, 8.86226892e-1f); // 0x1.c5bf88p-1
}
r = a * p;
return r;
}
/* compute natural logarithm with a maximum error of 0.85089 ulp */
float my_logf (float a)
{
float i, m, r, s, t;
int e;
m = frexpf (a, &e);
if (m < 0.666666667f) { // 0x1.555556p-1
m = m + m;
e = e - 1;
}
i = (float)e;
/* m in [2/3, 4/3] */
m = m - 1.0f;
s = m * m;
/* Compute log1p(m) for m in [-1/3, 1/3] */
r = -0.130310059f; // -0x1.0ae000p-3
t = 0.140869141f; // 0x1.208000p-3
r = fmaf (r, s, -0.121484190f); // -0x1.f19968p-4
t = fmaf (t, s, 0.139814854f); // 0x1.1e5740p-3
r = fmaf (r, s, -0.166846052f); // -0x1.55b362p-3
t = fmaf (t, s, 0.200120345f); // 0x1.99d8b2p-3
r = fmaf (r, s, -0.249996200f); // -0x1.fffe02p-3
r = fmaf (t, m, r);
r = fmaf (r, m, 0.333331972f); // 0x1.5554fap-2
r = fmaf (r, m, -0.500000000f); // -0x1.000000p-1
r = fmaf (r, s, m);
r = fmaf (i, 0.693147182f, r); // 0x1.62e430p-1 // log(2)
if (!((a > 0.0f) && (a <= 3.40282346e+38f))) { // 0x1.fffffep+127
r = a + a; // silence NaNs if necessary
if (a < 0.0f) r = ( 0.0f / 0.0f); // NaN
if (a == 0.0f) r = (-1.0f / 0.0f); // -Inf
}
return r;
}
关于c - C 中的逆误差函数,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/27229371/