C中的逆误差函数

Question

是否可以在 C 中计算逆误差函数？

我可以在计算误差函数的<math.h> 中找到erf(x)，但我找不到任何逆向函数。

Answer 1

目前，ISO C 标准数学库不包括erfinv()，或其单精度变体erfinvf()。但是，创建自己的版本并不太难，我在下面通过具有合理准确性和性能的erfinvf() 实现来演示。

查看graph of the inverse error function，我们观察到它是高度非线性的，因此难以用多项式逼近。处理这种情况的一种策略是通过将更简单的基本函数（其本身可以以高性能和极好的精度计算）和更容易服从多项式逼近或有理数的相当线性的函数组合来“线性化”这样的函数低度的近似值。

这里有一些文献中已知的erfinv 线性化方法，所有这些方法都基于对数。通常，作者将逆误差函数的主要的、相当线性的部分从零到非常大致大约 0.9 的切换点和从切换点到统一的 tail 部分区分开来。以下，log()表示自然对数，R()表示有理逼近，P()表示多项式逼近。

A. J. Strecok，“关于误差函数倒数的计算。” 计算数学，卷。 22，第 101 期（1968 年 1 月），第 144-158 页(online)

β(x) = (-log(1-x²]))^½; erfinv(x) = x · R(x²) [main]; R(x) · β(x) [尾]

J。 M. Blair, C. A. Edwards, J. H. Johnson，“误差函数倒数的有理切比雪夫近似”。 计算数学，卷。 30, No. 136 (Oct. 1976), pp. 827-830 (online)

ξ = (-log(1-x))^-½; erfinv(x) = x · R(x²) [main]; ξ^-1 · R(ξ) [尾]

M. Giles，“近似 erfinv 函数。” GPU Computing Gems Jade 版，第 109-116 页。 2011.(online)

w = -log(1-x²); s = √w; erfinv(x) = x · P(w) [主]； x · P(s) [尾]

下面的解决方案通常遵循 Giles 的方法，但对其进行了简化，不需要尾部的平方根，即它使用 x · P(w) 类型的两个近似值。该代码最大限度地利用了融合乘加运算 FMA，它通过 C 中的标准数学函数 fma() 和 fmaf() 公开。许多常见的计算平台，例如 IBM Power、Arm64、x86-64 和 GPU 在硬件中提供此操作。在不存在硬件支持的情况下，使用fma{f}() 可能会使下面的代码慢得无法接受，因为操作需要由标准数学库模拟。此外，FMA are known to exist 的模拟功能不正确。

标准数学库的对数函数logf()的精度会影响下面my_erfinvf()的精度。只要该库提供了一个错误 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;
}

Answer 2

又快又脏，公差在 +-6e-3 以下。工作基于 Sergei Winitzki 的“误差函数及其逆函数的方便近似”。

C/C++ 代码：

float myErfInv2(float x){
   float tt1, tt2, lnx, sgn;
   sgn = (x < 0) ? -1.0f : 1.0f;

   x = (1 - x)*(1 + x);        // x = 1 - x*x;
   lnx = logf(x);

   tt1 = 2/(PI*0.147) + 0.5f * lnx;
   tt2 = 1/(0.147) * lnx;

   return(sgn*sqrtf(-tt1 + sqrtf(tt1*tt1 - tt2)));
}

MATLAB 完整性检查：

clear all,  close all, clc 

x = linspace(-1, 1,10000);
% x = 1 - logspace(-8,-15,1000);

a = 0.15449436008930206298828125;
% a = 0.147;

u = log(1-x.^2);  
u1 = 2/(pi*a) + u/2;    u2 = u/a;
y = sign(x).*sqrt(-u1+sqrt(u1.^2 - u2)); 

f = erfinv(x); axis equal

figure(1);
plot(x, [y; f]); legend('Approx. erf(x)', 'erf(x)')

figure(2);
e = f-y;
plot(x, e);

MATLAB 绘图：

Answer 3

我不认为这是<math.h> 中的standard 实现，但还有其他C math 库实现了逆错误函数erfinv(x)，您可以使用。

Answer 4

又快又脏：如果允许的精度低于我与反双曲正切共享我自己的近似值 - 参数由蒙特卡尔模拟寻求，其中所有随机值都在 0.5 和 1.5 之间：

p1 = 1.4872301551536515
p2 = 0.5739159012216655
p3 = 0.5803635928651558

atanh( p^( 1 / p3 ) ) / p2 )^( 1 / p1 )

这来自我的 erf 函数近似与双曲正切的代数重新排序，其中 RMSE 误差为 0.000367354，x 介于 1 和 4 之间：

tanh( x^p1 * p2 )^p3