boost :: math :: sinc_pi是否不必要地复杂？

Question

这不是关于模板黑客或处理编译器怪癖的问题。我理解为什么Boost库就是这样。这是关于sinc_piBoost数学库中该函数使用的实际算法的。

---

该功能sinc(x)等效于sin(x)/x。

在Boost数学库的文档中sinc_pi()，它说“ Taylor系列最初用于确保准确性”。这似乎是荒谬的，因为浮点数的除法运算不会比乘法运算造成更多的精度损失。除非的特定实现存在错误，否则sin的天真的方法

double sinc(double x) {if(x == 0) return 1; else return sin(x)/x;}

似乎会没事的。

我已经测试了这一点，朴素版本与Boost数学工具包中的版本之间的最大相对差异仅为所用类型的epsilon的一半左右（对于float和double而言），这使其与离散化的规模相同错误。此外，此最大差异不会出现在0附近，而是出现在Boost版本使用部分泰勒级数（即abs(x) < epsilon**(1/4)）的时间间隔的结尾处。这看起来实际上是泰勒级数逼近，实际上是（非常轻微的）错误，这是由于区间末尾附近的精度损失或由于多次运算引起的反复舍入。

这是我编写的用于测试此结果的程序的结果，该程序对float0到1之间的每一个进行迭代，并计算Boost结果与朴素的结果之间的相对差：

Test for type float:

Max deviation from Boost result is 5.96081e-08 relative difference
equals 0.500029 * epsilon
at x = 0.0185723
which is epsilon ** 0.25003

这是程序的代码。它可以用于对任何浮点类型执行相同的测试，并且运行大约一分钟。

#include <cmath>
#include <iostream>

#include "boost/math/special_functions/sinc.hpp"

template <class T>
T sinc_naive(T x) { using namespace std; if (x == 0) return 1; else return sin(x) / x; }



template <class T>
void run_sinc_test()
{
    using namespace std;

    T eps = std::numeric_limits<T>::epsilon();

    T max_rel_err = 0;
    T x_at_max_rel_err = 0;

    for (T x = 0; x < 1; x = nextafter(static_cast<float>(x), 1.0f))
    {
        T boost_result = boost::math::sinc_pi(x);
        T naive_result = sinc_naive(x);

        if (boost_result != naive_result)
        {
            T rel_err = abs(boost_result - naive_result) / boost_result;

            if (rel_err > max_rel_err)
            {
                max_rel_err = rel_err;
                x_at_max_rel_err = x;
            }
        }
    }

    cout << "Max deviation from Boost result is " << max_rel_err << " relative difference" << endl;
    cout << "equals " << max_rel_err / eps << " * epsilon" << endl;
    cout << "at x = " << x_at_max_rel_err << endl;
    cout << "which is epsilon ** " << log(x_at_max_rel_err) / log(eps) << endl;
    cout << endl;
}



int main()
{
    using namespace std;

    cout << "Test for type float:" << endl << endl;
    run_sinc_test<float>();
    cout << endl;
    cin.ignore();
}

Answer 1

经过一番调查，我从原作者那里挖出了一个讨论。

\n

\n

[sin(x)] 在 x=0 时表现良好，sinc(x) 也是如此。[\xe2\x80\xa6] 我的解决方案\n将具有更好的性能或较小的参数，即|x| < pow(x, 1/6)，\n因为大多数处理器需要比\n1- (1/6) * x *x 更多的时间来计算 sin(x)。

\n

\n

来自https://lists.boost.org/Archives/boost/2001/05/12421.php。

\n

我发现的最早的关于使用泰勒展开来确保准确性的参考文献是很久以后的，并且是由另一个人提出的。所以看起来这与性能有关，而不是准确性。如果您想确定，您可能需要与相关人员联系。

\n\n

具体来说sinc_pi，我发现了以下交流。请注意，它们用于sinc_a指代以下形式的函数族sin(x*a)/(x*a)。

\n

\n

\n

sinc_a(x) 的优点是什么？解决非常大的 x 的舍入问题？那么对于非常大的参数改进 sin(x) 就更重要了。

\n

\n

这个家族的这个特定成员的主要兴趣在于它需要更少的计算，而且它本身就是一个特殊的函数，因为它比它的兄弟更常见。

\n

\n

来自https://lists.boost.org/Archives/boost/2001/05/12485.php。

\n