precision - 计算 (x-sin(x))/x^3 的数值稳定方法

Question

我正在寻找一种使用 IEEE 浮点算法和标准三角函数为所有 x 计算 (x-sin(x))/(x^3) 的非常快速的方法。在 0 时，它应该返回 1/6。

对于 sin(x)/x，检查 x=0 并返回 1 就足够了，否则只需使用标准浮点 sin 和除法来计算它。对于 (1-cos(x))/x^2，如果 cos(x) <= 0，则此表达式可以按原样使用，否则表示为 (sin(x)/x)^2/(1+cos(x ))

但我不知道如何表达 (x-sin(x))/x^3。

到目前为止，我所拥有的最好的方法是计算无限和直到它收敛:$\sum_0^{\infty}{1/4^(n+1)sin(x/2^n)/(x/2^n)(1-cos(x/2^n))/(x/2^n)^2}$ 但我更喜欢封闭形式

Answer 1

(1 - cos x)/x² 与 (x - sin x)/x³ 根本不同/em>，因为单位可以通过三角函数构造为 sin² x + cos² x = 1，而 sin² x = 1 x 为真。这意味着我们无法将后一个函数转换为数值上有利的封闭式三角函数公式。我想了很久，也尝试用我知道的所有三角恒等式来操纵公式。我很想被证明是错误的；这似乎是 Math Stack Exchange 的一个问题。实现前一个功能最简单最准确的方法是

// (1-cos(x))/x**2
double cosm1_over_xsquared (double x)
{
    if (fabs (x) < sqrt (DBL_EPSILON)) {
        return 0.5;
    } else {
        double s = sin (x * 0.5) / x;
        return 2.0 * s * s;
    }
}

如果标准数学库计算 sin() 时误差略高于半个 ulp，则此实现计算 (1 - cos x)/x² 误差不超过 4 ulp。作为旁注，此函数也适用于 Kahan 的 self 校正技术的使用，他首先展示了用于计算 (e^x - 1)/x 在

William M. Kahan，“提议的 IEEE 浮点算术标准中的区间算术”。载于 Karl L. E. Nickel(编)，Interval Arithmetic 1980，Academic Press 1980，第 99-128 页。

// (1-cos(x))/x**2 on [-3, 3] using Kahan's self-compensation technique
double cosm1_over_xsquared_kahan (double x)
{
    double u = cos (x);
    double n = 1.0 - u;
    if (n == 0.0) {
        return 0.5;
    }
    double d = acos (u);
    return n / (d * d);
}

如果 cos() 和 acos() 的最大误差略高于 ulp 的一半，则此函数返回的结果误差小于 5 ulp。因为 cos，而不是 e^x 是周期函数，所以这种方法仅适用于代码注释中注明的受限区间。

以上建议我们应该以最大误差约为 4 ulp 的方式实现 (x - sin x)/x³。表征朴素计算，我们发现它足以满足 |x| > 1 根据此规定。尽管替代计算的输入域很窄，但 Kahan 的 self 补偿技术不对该函数不起作用。然而，数学函数实现者的旧备用，多项式极小极大逼近，工作得很好。这导致以下代码:

// (x-sin(x))/x**3
double sinmx_over_xcubed (double x)
{
    if (fabs(x) < 1.0) { // minimax approximation
        double x2 = x * x;
        double p =   7.5475867852548673E-13;
        p = p * x2 - 1.6057658525730946E-10;
        p = p * x2 + 2.5052098906959416E-8;
        p = p * x2 - 2.7557319191306421E-6;
        p = p * x2 + 1.9841269841218293E-4;
        p = p * x2 - 8.3333333333333055E-3;
        p = p * x2 + 1.6666666666666666E-1;
        return p;
    } else {
        return (x - sin (x)) / x / x / x;
    }
}