准确高效地计算反朗之万函数

计算科学 算法 效率 准确性 特殊功能
2021-12-09 10:15:52

朗之万函数L(x)=coth(x)1x发生在与弹性体和顺磁性材料相关的计算中。只要解决了原点附近的减法抵消问题,例如通过在该区域中使用极小极大多项式近似,就可以轻松准确且高性能地计算它。

为了完全消除对消问题,选择一个良好的多项式近似转换点很重要;一些实现统一切换,这是次优的。更接近的切换点2.0是可取的。忠实地实现指数函数exp,通常由现代数学库提供,最大误差为1.3ulp 可以直接实现。示例性 C 实现如下所示。

然而,许多应用程序需要朗之万函数L1(x),没有封闭形式的表示;需要数值近似。文献中使用的大多数近似值都非常不准确,相对误差在 0.05% 到 2% 的范围内。有用的概述由以下人员提供:

Benjamin C. Marchi 和 Ellen M. Arruda,“逆朗之万近似的误差最小化方法”。流变学报(2015) 54:887-902

Martin Kröger,“朗之万和布里渊逆函数的简单、可接受和准确的近似值,与强聚合物变形和流动相关。” 非牛顿流体力学杂志223 (2015) 77-87

最近的工作表明,反朗之万函数的计算缺乏准确性对其使用的科学计算的整体准确性产生负面影响:

Amine Ammar,“逆朗之万近似对非线性稀聚合物福克-普朗克方程解的影响。” 非牛顿流体力学杂志231(2016):1-5。

因此,非常需要准确计算反朗之万函数。除了准确性之外,一些用例还受益于特定的功能形式,以便于替换、区分和集成。抛开这些额外的要求,专注于最小化相对误差,反朗之万函数如何L1(x)计算准确(接近机器精度)并具有良好的性能?

#include <math.h>  // fabs, exp, copysign

/* compute Langevin function L(a) = coth(a) - 1/a. Max error typically < 1.3 ulp */
float langevinf (float a)
{
    float r, t;

    t = fabsf (a);
    if (t > 1.8125f) {
        float e;
        e = 1.0f / (expf (2.0f * t) - 1.0f);
        r = (1.0f - 1.0f / t) + 2.0f * e;
        r = copysignf (r, a);
    } else {
        float s;
        s = a * a;
        r =         7.70960469e-8f;
        r = r * s - 1.65101926e-6f;
        r = r * s + 2.03457112e-5f;
        r = r * s - 2.10521728e-4f;
        r = r * s + 2.11580913e-3f;
        r = r * s - 2.22220998e-2f;
        r = r * s + 8.33333284e-2f;
        r = r * a + 0.25f * a;
    }
    return r;
}

/* compute Langevin function L(a) = coth(a) - 1/a. Max error typically < 1.3 ulp */
double langevin (double a)
{
    double r, t;

    t = fabs (a);
    if (t > 1.8125) {
        double e;
        e = 1.0 / (exp (2.0 * t) - 1.0);
        r = (1.0 - 1.0 / t) + 2.0 * e;
        r = copysign (r, a);
    } else {
        double s, t, u;
        s = a * a;
        t = s * s;
        r = 2.6224020345921745e-16 * s - 9.0382507125284313e-15;
        u = 1.6058476918889783e-13 * s - 2.0544032928760129e-12;
        r = r * t + u;
        u = 2.2318295632954910e-11 * s - 2.2645716105162602e-10;
        r = r * t + u;
        u = 2.2484176131435534e-09 * s - 2.2212056151537571e-08;
        r = r * t + u;
        u = 2.1925751794205211e-07 * s - 2.1644032423757088e-06;
        r = r * t + u;
        u = 2.1377798795067705e-05 * s - 2.1164021156459166e-04;
        r = r * t + u;
        u = 2.1164021163937877e-03 * s - 2.2222222222221859e-02;
        r = r * t + u;
        r = r * s + 8.3333333333333329e-02;
        r = r * a + 0.25 * a;
    }
    return r;
}
1个回答

反朗之万函数L1(x)是奇函数。因此,只需要考虑区间的近似值[0,1]; 负半平面通过关于原点的对称处理。由于单位的奇异性,多项式近似不是研究的富有成果的方向,尽管一些作者已经探索了零附近的高度泰勒展开,例如

Mikhail Itskov、Roozbeh Dargazany 和 Karl Hörnes,“反函数的泰勒展开与朗之万函数的应用”。固体数学与力学17, No. 7 (2012): 693-701

由于需要大量计算,论文中的 115 项级数展开并不实用。有理逼近非常适合逼近具有奇点的函数,并且文献中的大多数工作都集中在这些方面,特别是 Pade 逼近的形式。近年来,该方法已通过增加对误差项的补偿来改进,例如

本杰明 C. 马奇和艾伦 M. 阿鲁达。“广义误差最小化,有理逆朗之万近似。” 固体数学与力学(2018):1081286517754131。

一些作者试图用三角函数来近似反朗之万函数的形状,特别是,tansin

Jörgen Stefan Bergström,“弹性材料的大应变随时间变化的行为”。博士论文,麻省理工学院,1999。

Grant Keady,“朗之万函数和截断指数分布”。arXiv 预印本 arXiv:1501.02535 (2015)。

Rafayel Petrosyan,“改进了一些聚合物扩展模型的近似值”。流变学学报56,第 56 期。1 (2017): 21-26。

虽然较新的作品将相对误差降低到大约105,我知道只有一篇论文试图达到接近机器精度。它通过使用具有数千个间隔的分段近似来做到这一点,每个间隔都适合三次样条。系数表的存储要求达到数百 KBytes。这超过了大多数处理器上一级缓存的大小,因此该方案在实践中可能不适合:

何塞·玛丽亚·贝尼特斯和弗朗西斯科·哈维尔·蒙坦斯,“一种简单而有效的数值程序,可以高精度地计算逆朗之万函数”。arXiv 预印本 arXiv:1806.08068 (2018)。

对于另一种方法,我注意到反朗之万函数和反误差函数之间的整体形状相似。正如我在另一个答案中指出的那样,后者可以通过使用变换来有效地近似t=log(1x2)然后生成形式的近似值xP(t), 在哪里P是多项式。自从L(x)11x接近统一,因此1x1是一个很好的近似值L1(x)在奇点附近。实验上,对于 IEEE-754 单精度,使用该近似值适用于[5764,1],保持最大误差小于2.6ulp。

我的实验表明,正如预期的那样,L1(x)以与逆误差函数相同的方式,在多项式逼近中只需要适度更多的项。然而,生成高质量的极小极大近似需要一些努力,因为近似问题的条件较差。从 Remez 算法开始生成初始近似值,我使用启发式细化(类似于模拟退火)来得出最终结果。我的初步实验表明,这种方法对于需要超过单精度精度的用例是可扩展的,但对于达到完全双精度精度来说计算成本太高,这需要 40 次左右的多项式。

大多数现代硬件平台都提供融合乘加运算 (FMA),有助于减少舍入误差并针对某些减法取消提供保护。它的使用还倾向于提高性能,因为执行时间通常相同或仅略高于浮点乘法。该操作通常通过库函数公开,例如fma()在 C 和 C++ 中。但是,这映射到不支持 FMA 的硬件平台上的缓慢仿真,并且某些仿真存在功能错误

IEEE-754 (2008) 单精度 ( ) 的示例 C 实现binary32如下所示。虽然达到的准确度将取决于标准数学库执行的质量log,大多数现代数学库都提供了忠实全面的实现。通过我尝试的三种不同实现,logf代码注释中说明的错误范围得以保持。

#include <math.h> // fabs, log, copysign, fma

/* 
  Compute inverse Langevin function accurate to almost machine precision

  USE_FMA == 0: max. ulp error < 4.27, max. relative error < 4.43e-7
  USE_FMA == 1: max. ulp error < 3.64, max. relative error < 3.84e-7
*/
float langevininvf (float x)
{
    float p, r, t;
    if ((fabsf (x) > 0.890625f) && (fabsf (x) <= 1.0f)) {
        r = copysignf (1.0f / (fabsf (x) - 1.0f), x);
    } else {
#if USE_FMA
        t = fmaf (x, 0.0f - x, 1.0f); // compute 1-x*x accurately
        t = logf (t);
        p =              2.18808651e-4f;  //  0x1.cae000p-13
        p = fmaf (p, t, -7.90076610e-3f); // -0x1.02e46ep-7
        p = fmaf (p, t, -7.12909698e-2f); // -0x1.240200p-4
        p = fmaf (p, t, -2.40409270e-1f); // -0x1.ec5bb2p-3
        p = fmaf (p, t, -4.14386481e-1f); // -0x1.a854eep-2
        p = fmaf (p, t, -4.05752033e-1f); // -0x1.9f7d76p-2
        p = fmaf (p, t, -2.56382942e-1f); // -0x1.068940p-2
        p = fmaf (p, t, -1.22061931e-1f); // -0x1.f3f736p-4
        p = fmaf (p, t,  5.00488468e-2f); //  0x1.9a000ap-5
        p = fmaf (p, t, -1.84208602e-1f); // -0x1.79425cp-3
        p = fmaf (p, t,  3.98338169e-1f); //  0x1.97e5f6p-2 
        p = fmaf (p, t, -9.00006115e-1f); // -0x1.cccd9ap-1 
        p = fmaf (p, t,  5.00000000e-1f); //  0x1.000000p-1
        t = x + x;
        r = fmaf (p, t, t);
#else // USE_FMA
        t = (1.0f - x) * (1.0f + x); // compute 1-x*x accurately
        t = logf (t);
        p =         2.18868256e-4f; //  0x1.cb0000p-13
        p = p * t - 7.90085457e-3f; // -0x1.02e52cp-7
        p = p * t - 7.12914243e-2f; // -0x1.24027ap-4
        p = p * t - 2.40408897e-1f; // -0x1.ec5b80p-3
        p = p * t - 4.14384902e-1f; // -0x1.a85484p-2
        p = p * t - 4.05750901e-1f; // -0x1.9f7d2ap-2
        p = p * t - 2.56382436e-1f; // -0x1.06891ep-2
        p = p * t - 1.22061789e-1f; // -0x1.f3f710p-4
        p = p * t + 5.00487871e-2f; //  0x1.99ffeap-5
        p = p * t - 1.84208602e-1f; // -0x1.79425cp-3
        p = p * t + 3.98338020e-1f; //  0x1.97e5ecp-2
        p = p * t - 9.00006115e-1f; // -0x1.cccd9ap-1
        p = p * t + 5.00000000e-1f; //  0x1.000000p-1
        t = x + x;
        r = p * t + t;
#endif // USE_FMA
    }
    return r;
}

为简单的超越函数提供快速硬件支持的计算平台有助于特别有效地实现所提出的算法。下面显示了一个使用 CUDA 的 GPU 实现,它只转换为大约 25 条指令。

/* max. ulp error < 3.925, max. relative error < 3.95e-7 */
__device__ float langevininvf (float x)
{
    float fa, p, r, t;
    fa = fabsf (x);
    if ((fa > (57.0f / 64.0f)) && (fa <= 1.0f)) {
        t = fa - 1.0f;
        asm ("rcp.approx.ftz.f32 %0,%0;" : "+f"(t));
        r = copysignf (t, x);
    } else {
        t = fmaf (x, -x, 1.0f);
        asm ("lg2.approx.ftz.f32 %0,%0;" : "+f"(t));
        p =              2.69152224e-6f;  //  0x1.694000p-19
        p = fmaf (p, t, -1.40199758e-4f); // -0x1.26052cp-13
        p = fmaf (p, t, -1.82510854e-3f); // -0x1.de70f6p-10
        p = fmaf (p, t, -8.87932349e-3f); // -0x1.22f52ap-7
        p = fmaf (p, t, -2.20804960e-2f); // -0x1.69c450p-6
        p = fmaf (p, t, -3.11916713e-2f); // -0x1.ff0b5ap-6
        p = fmaf (p, t, -2.84342300e-2f); // -0x1.d1ddcep-6
        p = fmaf (p, t, -1.95302274e-2f); // -0x1.3ffbb6p-6
        p = fmaf (p, t,  1.15530221e-2f); //  0x1.7a91c6p-7
        p = fmaf (p, t, -6.13459945e-2f); // -0x1.f68be0p-5
        p = fmaf (p, t,  1.91382647e-1f); //  0x1.87f3a0p-3
        p = fmaf (p, t, -6.23836875e-1f); // -0x1.3f678cp-1
        p = fmaf (p, t,  5.00000000e-1f); //  0x1.000000p-1
        t = x + x;
        r = fmaf (p, t, t);
    }
    return r;
}