在C++中最快的正弦,余弦和平方根的实现(不需要太精确)

Question

我正在谷歌搜索过去一小时的问题,但只有泰勒系列或一些示例代码的要点太慢或根本不编译.好吧,我发现谷歌的答案大多是"Google it,它已经被问到了",但遗憾的是它不是 ......

我在低端Pentium 4上分析我的游戏,发现大约85%的执行时间浪费在计算窦,cosinus和平方根(来自Visual Studio中的标准C++库)上,这似乎与CPU密切相关(在我的I7上,相同的功能只有5%的执行时间,并且游戏更快了waaaaaaaaaa).我不能优化这三个函数,也不能在一次传递中计算正弦和余弦(相互依赖),但我不需要太精确的模拟结果,所以我可以使用更快的逼近.

那么,问题是:在C++中计算float的正弦,余弦和平方根的最快方法是什么？

编辑 查找表更加痛苦,因为在现代CPU上产生的Cache Miss比Taylor系列更昂贵.这些天CPU很快,而缓存则不然.

我犯了一个错误,我虽然需要为Taylor系列计算几个阶乘,我现在看到它们可以实现为常量.

所以更新的问题是:对于平方根还有任何快速优化吗？

EDIT2

我使用平方根计算距离,而不是规范化 - 不能使用快速反平方根算法(如评论中所指出:http://en.wikipedia.org/wiki/Fast_inverse_square_root

EDIT3

我也无法在平方距离上操作,我需要精确的距离进行计算

Answer 1

首先,泰勒系列不是实现正弦/余弦的最佳/最快方式.它也不是专业库实现这些三角函数的方式,并且知道最佳数值实现允许您调整精度以更有效地获得速度.此外,StackOverflow中已经广泛讨论了这个问题.这只是一个例子.

其次,您在旧/新PCS之间看到的巨大差异是由于现代英特尔架构具有用于计算元素三角函数的显式汇编代码.在执行速度上击败他们是相当困难的.

最后,我们来谈谈旧PC上的代码.检查gsl gnu科学库 (或数值配方)的实现,你会发现他们基本上使用了Chebyshev近似公式.

切比雪夫逼近收敛速度更快,因此您需要评估更少的术语.我不会在这里编写实现细节,因为StackOverflow上已经发布了非常好的答案.例如,检查一下.只需调整此系列中的术语数量即可更改精度/速度之间的平衡.

顺便说一下:针对这类问题的规则0:如果你想要一些特殊函数或数值方法的实现细节,你应该在任何进一步的行动之前看一下GSL代码 - GSL是标准的数值库.

编辑:您可以通过在gcc/icc中包含积极的浮点优化标志来缩短执行时间.这会降低精度,但似乎这正是你想要的.

编辑2:您可以尝试制作粗网格并使用gsl例程(gsl_interp_cspline_periodic用于具有周期条件的样条曲线)来对该表进行样条化(与线性插值相比,样条线将减少误差=>您需要更少的点在您的表格上= >减少缓存未命中率)!

Answer 2

这是C++中保证最快的正弦函数:

double FastSin(double x)
{
    return 0;
}

哦,你想要比1.0 |更准确的准确度？好好读一读.

20世纪70年代的工程师在这个领域取得了一些奇妙的发现,但新的程序员根本不知道这些方法存在,因为它们不是作为标准计算机科学课程的一部分进行教学的.

您需要首先了解所有应用程序都没有"完美"实现这些功能.因此,对"哪一个最快"等问题的表面答案保证是错误的.

大多数提出这个问题的人都不理解在性能和准确性之间进行权衡的重要性.特别是,在进行任何其他操作之前,您必须对计算的准确性做出一些选择.您可以在结果中容忍多少错误？10 ^ -4？10 ^ -16？

除非您可以在任何方法中量化错误,否则请不要使用它.

没有人单独使用泰勒系列来近似软件中的超越.除了某些高度特殊的情况,泰勒系列通常在共同的输入范围内缓慢接近目标.

祖父母用来有效计算超验数的算法统称为CORDIC,并且很简单,可以在硬件中实现.以下是C中记录良好的CORDIC实现.CORDIC实现通常需要非常小的查找表,但大多数实现甚至不需要硬件乘法器可用.大多数CORDIC实现允许您在准确性方面进行权衡,包括我链接的那个.

多年来,对原始CORDIC算法进行了大量的渐进式改进.例如,去年日本的一些研究人员发表了一篇关于改进的CORDIC 的文章,该文章具有更好的旋转角度,从而减少了所需的操作.

如果你有硬件乘数(你几乎可以肯定),或者如果你买不起CORDIC要求的查找表,你总是可以使用Chebyshev多项式来做同样的事情.切比雪夫多项式需要乘法,但这在现代硬件上很少成为问题.我们喜欢切比雪夫多项式,因为它们对于给定的近似具有高度可预测的最大误差.在输入范围内,切比雪夫多项式中最后一项的最大值限制了结果中的误差.随着术语数量的增加,这个错误会变小. 这是切比雪夫多项式的一个例子,它给出了一个巨大范围内的正弦近似,忽略了正弦函数的自然对称性,只是通过在其上抛出更多系数来解决近似问题.

我们也喜欢切比雪夫多项式,因为近似误差在输出范围内均匀分布.如果你正在写的音频插件或进行数字信号处理,切比雪夫多项式给你一个便宜的和可预见的抖动效果"是免费的."

如果你想在特定范围内找到自己的切比雪夫多项式系数,许多数学库都会调用找到那些系数" Chebyshev fit "或类似的系数.

现在,平方根通常使用Newton-Raphson算法的一些变体计算,通常具有固定数量的迭代.通常情况下,当有人为了做平方根而采用"惊人的新"算法时,它只是伪装的Newton-Raphson.

Newton-Raphson,CORDIC和Chebyshev多项式让您可以权衡速度以获得准确性,因此答案可能与您想要的一样不精确.

最后,当您完成所有花哨的基准测试和微优化后,请确保您的"快速"版本实际上比库版本更快. 这是fsin()的典型库实现, 在域中从-pi/4到pi/4.它只是不是那么该死.

有些人一生致力于有效地解决这些问题,并产生了一些令人着迷的结果.当您准备加入旧学校时,请拿一份数字食谱.

TL:博士; go google"sine approximation"或"cosine approximation"或"square root approximation"或" approximation theory".

Answer 3

对于平方根,有一种称为位移的方法.

由IEEE-754定义的浮点数使用某些位表示基于多个的2的描述时间.一些位用于表示基值.

float squareRoot(float x)
{
  unsigned int i = *(unsigned int*) &x;

  // adjust bias
  i  += 127 << 23;
  // approximation of square root
  i >>= 1;

  return *(float*) &i;
}

这是计算方根的恒定时间

Answer 4

最快的方法是预先计算值,使用像这个示例中的表:

用C++创建正弦查找表

但如果你坚持在运行时计算,你可以使用正弦或余弦的泰勒级数展开...

有关泰勒系列的更多信息,请访问http://en.wikipedia.org/wiki/Taylor_series

让它运行良好的关键之一是预先计算因子并截断合理数量的项.阶乘在分母中迅速增长,因此您不需要携带超过几个术语.

另外......每次都不要将x ^ n从开头乘以...例如,将x ^ 3乘以x再乘两次,然后乘以另外两次来计算指数.

Answer 5

基于http://forum.devmaster.net/t/fast-and-accurate-sine-cosine/9648的想法和一些手动重写以提高微基准测试的性能,我最终得到了以下余弦实现,用于HPC物理模拟,该模拟由大量空间上的重复cos调用引起的瓶颈.它比查询表更准确,更快,最明显的是不需要除法.

template<typename T>
inline T cos(T x) noexcept
{
    constexpr T tp = 1./(2.*M_PI);
    x *= tp;
    x -= T(.25) + std::floor(x + T(.25));
    x *= T(16.) * (std::abs(x) - T(.5));
    #if EXTRA_PRECISION
    x += T(.225) * x * (std::abs(x) - T(1.));
    #endif
    return x;
}

在循环中使用时,至少英特尔编译器在向量化这个函数方面也足够聪明.

如果EXTRA_PRECISION被定义,最大误差为-π〜π范围在约0.00109,假定T是double因为它在大多数C++实现通常定义.否则,相同范围的最大误差约为0.056.

Answer 6

我尝试了 millianw 的答案，它给了我 4.5 倍的加速，所以这太棒了。

然而，millianw 链接到的原始文章计算正弦，而不是余弦，并且它的做法有些不同。（看起来更简单。）

可以预见的是，15 年后，该文章的 URL ( http://forum.devmaster.net/t/fast-and-accurate-sine-cosine/9648 ) 今天给出了 404，所以我通过 archive.org 获取了它，然后我我将其添加到这里供后代使用。

不幸的是，尽管该文章包含多个图像，但 archive.org 只存储了前 2 个图像。另外，作者的个人资料页面（http://forum.devmaster.net/users/Nick）没有被存储，所以我想我们永远不会知道尼克是谁。

===================================================

快速准确的正弦/余弦

尼克 06 年 4 月

大家好，

在某些情况下，您需要以非常高性能运行的正弦和余弦的良好近似值。一个例子是实现圆形表面的动态细分，与 Quake 3 中的类似。或者实现波浪运动，以防没有可用的顶点着色器 2.0。

标准 C sinf() 和 cosf() 函数非常慢，并且提供的精度比我们实际需要的要高得多。我们真正想要的是一种能够在精度和性能之间提供最佳折衷的近似值。最著名的近似方法是使用约 0 的泰勒级数（也称为麦克劳林级数），对于正弦，它变为：

x - 1/6 x^3 + 1/120 x^5 - 1/5040 x^7 + ...

当我们绘制它时，我们得到：taylor.gif。

绿线是实正弦，红线是泰勒级数的前四项。这似乎是一个可以接受的近似值，但让我们仔细看看：taylor_zoom.gif。

在 pi/2 之前它的表现非常好，但之后它很快就会偏离。在 pi 处，它的计算结果为 -0.075 而不是 0。使用它进行波浪模拟将导致不可接受的急速运动。

我们可以添加另一个项，这实际上可以显着减少误差，但这使得公式相当冗长。对于 4 项版本，我们已经需要 7 次乘法和 3 次加法。泰勒级数无法为我们提供所需的精度和性能。

然而，我们确实注意到我们需要 sine(pi) = 0。从 taylor_zoom.gif 中我们还可以看到另一件事：这看起来非常像抛物线！因此，让我们尝试找到尽可能与其匹配的抛物线公式。抛物线的通用公式为 A + B x + C x\^2。所以这给了我们三个自由度。明显的选择是我们想要 sine(0) = 0、sine(pi/2) = 1 和 sine(pi) = 0。这给了我们以下三个方程：

A + B 0 + C 0^2 = 0
A + B pi/2 + C (pi/2)^2 = 1
A + B pi + C pi^2 = 0

其解为 A = 0，B = 4/pi，C = -4/pi\^2。所以我们的抛物线近似变为 4/pi x - 4/pi\^2 x\^2。绘制此图我们得到：parabola.gif。这看起来比 4 项泰勒级数更糟糕，对吗？错误的！最大绝对误差为 0.056。此外，这种近似将为我们提供平滑的波动运动，并且只需 3 次乘法和 1 次加法即可计算！

不幸的是它还不太实用。这就是我们在 [-pi, pi] 范围内得到的：负数.gif。很明显我们至少需要一个完整的周期。但同样清楚的是，这只是另一条抛物线，围绕原点镜像。其公式为 4/pi x + 4/pi\^2 x\^2。所以直接的（伪 C）解决方案是：

if(x > 0)
{
    y = 4/pi x - 4/pi^2 x^2;
}
else
{
    y = 4/pi x + 4/pi^2 x^2;
}

不过，添加分支并不是一个好主意。它使代码明显变慢。但看看这两个部分到底有多相似。根据 x 的符号，减法变为加法。在首次尝试消除分支时，我们可以使用 x / abs(x)“提取”x 的符号。除法的开销非常大，但看看结果公式：4/pi x - x / abs(x) 4/pi\^2 x\^2。通过反转除法，我们可以将其简化为非常漂亮且干净的 4/pi x - 4/pi\^2 x abs(x)。因此，只需进行一项额外操作，我们就可以获得正弦近似值的两半！下面是该公式的图表，确认了结果：abs.gif。

现在让我们看看余弦。基本三角学告诉我们 cosine(x) = sine(pi/2 + x)。就是这样，将 pi/2 加到 x 上吗？不，我们实际上再次得到了抛物线不需要的部分：shift_sine.gif。我们需要做的是当 x > pi/2 时“环绕”。这可以通过减去 2 pi 来完成。所以代码就变成了：

x += pi/2;

if(x > pi)   // Original x > pi/2
{
    x -= 2 * pi;   // Wrap: cos(x) = cos(x - 2 pi)
}

y = sine(x);

又一个分支。为了消除它，我们可以使用二进制逻辑，如下所示：

x -= (x > pi) & (2 * pi);

请注意，这根本不是有效的 C 代码。但它应该澄清这是如何运作的。当 x > pi 为 false 时，& 运算会将右侧部分归零，因此减法不会执行任何操作，这是完全等价的。我将把它作为练习留给读者来为此创建工作 C 代码（或者只是继续阅读）。显然，余弦比正弦需要更多的操作，但似乎没有其他方法，而且它仍然非常快。

现在，最大误差为 0.056 已经不错了，但显然 4 项泰勒级数的平均误差仍然较小。回想一下我们的正弦是什么样子的：abs.gif。那么我们能做些什么来以最小的成本进一步提高精度呢？当然，当前版本已经适用于许多情况，其中看起来像正弦的东西与真正的正弦一样好。但对于其他情况，这还不够好。

查看这些图表，您会发现我们的近似值总是高估真实的正弦值，除了 0、pi/2 和 pi 以外。因此，我们需要的是在不触及这些重要点的情况下“缩小规模”。解决方案是使用平方抛物线，如下所示：squared.gif。请注意它如何保留这些重要点，但它总是低于真实的正弦值。因此我们可以使用两者的加权平均值来获得更好的近似值：

Q (4/pi x - 4/pi^2 x^2) + P (4/pi x - 4/pi^2 x^2)^2

当 Q + P = 1 时。您可以使用绝对误差或相对误差的精确最小化，但我会为您节省数学运算。绝对误差的最佳权重为 Q = 0.775，P = 0.225，相对误差的最佳权重为 Q = 0.782，P = 0.218。我会使用前者。生成的图表是：average.gif。红线去哪儿了？它几乎完全被绿线覆盖，这立即显示了这个近似值的真实效果。最大误差约为 0.001，提高了 50 倍！公式看起来很长，但括号之间的部分与抛物线的值相同，只需计算一次。事实上，只需要 2 次额外的乘法和 2 次加法即可实现这种精度提升。

为了使其也适用于负 x，我们需要第二个 abs() 操作，这应该不足为奇。正弦的最终 C 代码变为：

float sine(float x)
{
    const float B = 4/pi;
    const float C = -4/(pi*pi);

    float y = B * x + C * x * abs(x);

    #ifdef EXTRA_PRECISION
    //  const float Q = 0.775;
        const float P = 0.225;

        y = P * (y * abs(y) - y) + y;   // Q * y + P * y * abs(y)
    #endif
}

所以我们只需要 5 次乘法和 3 次加法；如果我们忽略 abs()，仍然比 4 项泰勒更快，而且更精确！余弦版本只需要对 x 进行额外的移位和换行操作。

最后但并非最不重要的一点是，如果我不包含 SIMD 优化汇编版本，我就不会成为 Nick。它允许非常有效地执行包装操作，所以我会给你余弦：

// cos(x) = sin(x + pi/2)
addps xmm0, PI_2
movaps xmm1, xmm0
cmpnltps xmm1, PI
andps xmm1, PIx2
subps xmm0, xmm1

// Parabola
movaps xmm1, xmm0
andps xmm1, abs
mulps xmm1, xmm0
mulps xmm0, B
mulps xmm1, C
addps xmm0, xmm1

// Extra precision
movaps xmm1, xmm0
andps xmm1, abs
mulps xmm1, xmm0
subps xmm1, xmm0
mulps xmm1, P
addps xmm0, xmm1

该代码并行计算四个余弦，对于大多数 CPU 架构来说，每个余弦的峰值性能约为 9 个时钟周期。理想情况下，正弦波只需要 6 个时钟周期。较低精度的版本甚至可以在每个正弦 3 个时钟周期运行...并且不要忘记 -pi 和 pi 之间的所有输入都是有效的，并且公式在 -pi、-pi/2、0、pi/2 和圆周率。

因此，结论是不要再使用泰勒级数来近似正弦或余弦！为了在本文中添加有用的讨论，我很想听听是否有人知道其他超越函数（如指数、对数和幂函数）的良好近似值。

干杯，

缺口

===================================================

通过访问网络存档页面，您可能还会发现本文后面的评论很有趣：

http://web.archive.org/web/20141220225551/http://forum.devmaster.net/t/fast-and-accurate-sine-cosine/9648

Answer 7

我使用以下CORDIC代码以四倍精度计算三角函数。常数 N 确定所需精度的位数（例如 N=26 将给出单精度精度）。根据所需的精度，预先计算的存储空间可以很小并且适合缓存。它只需要加法和乘法运算，也很容易矢量化。

该算法预先计算 0.5^i, i=1,...,N 的 sin 和 cos 值。然后，我们可以组合这些预先计算的值，计算任何角度的 sin 和 cos，分辨率高达 0.5^N

template <class QuadReal_t>
QuadReal_t sin(const QuadReal_t a){
  const int N=128;
  static std::vector<QuadReal_t> theta;
  static std::vector<QuadReal_t> sinval;
  static std::vector<QuadReal_t> cosval;
  if(theta.size()==0){
    #pragma omp critical (QUAD_SIN)
    if(theta.size()==0){
      theta.resize(N);
      sinval.resize(N);
      cosval.resize(N);

      QuadReal_t t=1.0;
      for(int i=0;i<N;i++){
        theta[i]=t;
        t=t*0.5;
      }

      sinval[N-1]=theta[N-1];
      cosval[N-1]=1.0-sinval[N-1]*sinval[N-1]/2;
      for(int i=N-2;i>=0;i--){
        sinval[i]=2.0*sinval[i+1]*cosval[i+1];
        cosval[i]=sqrt(1.0-sinval[i]*sinval[i]);
      }
    }
  }

  QuadReal_t t=(a<0.0?-a:a);
  QuadReal_t sval=0.0;
  QuadReal_t cval=1.0;
  for(int i=0;i<N;i++){
    while(theta[i]<=t){
      QuadReal_t sval_=sval*cosval[i]+cval*sinval[i];
      QuadReal_t cval_=cval*cosval[i]-sval*sinval[i];
      sval=sval_;
      cval=cval_;
      t=t-theta[i];
    }
  }
  return (a<0.0?-sval:sval);
}

Answer 8

对于x86,硬件FP平方根指令很快(sqrtss是sqrt Scalar单精度).单精度比双精度快,所以绝对使用float而不是double代码,你可以负担得起使用较少的精度.

对于32位代码,通常需要编译器选项才能使用SSE指令进行FP数学运算,而不是x87.(例如-mfpmath=sse)

要使C sqrt()或sqrtf()函数内联为sqrtsd或sqrtss,您需要编译-fno-math-errno.在errnoNaN上设置数学函数通常被认为是设计错误,但标准需要它.如果没有该选项,gcc会内联它,但会执行compare +分支以查看结果是否为NaN,如果是,则调用库函数以便进行设置errno.如果您的程序errno在数学函数后没有检查,则使用时没有危险-fno-math-errno.

你不需要任何"不安全"的部分-ffast-math来获取sqrt和其他一些函数更好地内联,但-ffast-math可以产生很大的不同(例如,允许编译器在更改结果的情况下自动向量化,因为FP数学不是关联的.

例如用gcc6.3编译float foo(float a){ return sqrtf(a); }

foo:    # with -O3 -fno-math-errno.
    sqrtss  xmm0, xmm0
    ret

foo:   # with just -O3
    pxor    xmm2, xmm2   # clang just checks for NaN, instead of comparing against zero.
    sqrtss  xmm1, xmm0
    ucomiss xmm2, xmm0
    ja      .L8          # take the slow path if 0.0 > a
    movaps  xmm0, xmm1
    ret

.L8:                     # errno-setting path
    sub     rsp, 24
    movss   DWORD PTR [rsp+12], xmm1   # store the sqrtss result because the x86-64 SysV ABI has no call-preserved xmm regs.
    call    sqrtf                      # call sqrtf just to set errno
    movss   xmm1, DWORD PTR [rsp+12]
    add     rsp, 24
    movaps  xmm0, xmm1    # extra mov because gcc reloaded into the wrong register.
    ret

gcc的NaN案例代码似乎过于复杂; 它甚至没有使用sqrtf返回值!无论如何,-fno-math-errno对于sqrtf()你的程序中的每个呼叫站点,这是你实际上没有的混乱.大多数情况下,这只是代码膨胀,并且.L8当采用> = 0.0的数字的sqrt时,没有一个块会运行,但在快速路径中仍然有几个额外的指令.

如果您知道输入为sqrt非零,则可以使用快速但非常近似的倒数sqrt指令rsqrtps(或rsqrtss标量版本).一次Newton-Raphson迭代使其达到与硬件单精度sqrt指令几乎相同的精度,但并不完全相同.

sqrt(x) = x * 1/sqrt(x),for x!=0,因此您可以使用rsqrt和乘法计算sqrt.这些都很快,即使在P4上(2013年仍然相关)？

在P4上,可能值得使用rsqrt+ Newton迭代来替换单个sqrt,即使您不需要通过它来划分任何东西.

又见回答我最近写了一篇关于计算处理时零sqrt(x)为x*rsqrt(x),用牛顿迭代法.如果您想将FP值转换为整数,我会包含一些关于舍入错误的讨论,并链接到其他相关问题.

P4:

• sqrtss:23c延迟,不流水线
• sqrtsd:38c延迟,不流水线
• fsqrt (x87):43c延迟,不流水线

• rsqrtss/ mulss:4c + 6c延迟.可能是每3c吞吐量一个,因为它们显然不需要相同的执行单元(mmx与fp).

• SIMD打包版本有点慢

SKYLAKE微架构:

• sqrtss/ sqrtps:12c延迟,每3c吞吐量一个
• sqrtsd/ sqrtpd:15-16c延迟,每4-6c吞吐量一个
• fsqrt (x87):14-21cc延迟,每4-7c吞吐量一个
• rsqrtss/ mulss:4c + 4c延迟.每1c吞吐量一个.
• SIMD 128b矢量版本速度相同.256b矢量版本的延迟稍高,几乎是吞吐量的一半.该rsqrtss版本具有256b向量的完整性能.

使用牛顿迭代,rsqrt如果更快,版本并不多.

来自Agner Fog实验测试的数字.查看他的微指南,了解使代码快速或慢速运行的原因.另请参阅x86标签wiki上的链接.

IDK如何最好地计算sin/cos.我已经读过硬件fsin/ fcos(并且fsincos两者一次只能稍微慢一点)不是最快的方式,但IDK是什么.

Answer 9

QT 具有使用带插值的查找表的正弦 (qFastSin) 和余弦 (qFastCos) 的快速实现。我在我的代码中使用它，它们比 std:sin/cos 更快，并且足够精确，满足我的需要：

https://code.woboq.org/qt5/qtbase/src/corelib/kernel/qmath.h.html#_Z8qFastSind

#define QT_SINE_TABLE_SIZE 256


inline qreal qFastSin(qreal x)
{
   int si = int(x * (0.5 * QT_SINE_TABLE_SIZE / M_PI)); // Would be more accurate with qRound, but slower.
   qreal d = x - si * (2.0 * M_PI / QT_SINE_TABLE_SIZE);
   int ci = si + QT_SINE_TABLE_SIZE / 4;
   si &= QT_SINE_TABLE_SIZE - 1;
   ci &= QT_SINE_TABLE_SIZE - 1;
   return qt_sine_table[si] + (qt_sine_table[ci] - 0.5 * qt_sine_table[si] * d) * d;
}

inline qreal qFastCos(qreal x)
{
   int ci = int(x * (0.5 * QT_SINE_TABLE_SIZE / M_PI)); // Would be more accurate with qRound, but slower.
   qreal d = x - ci * (2.0 * M_PI / QT_SINE_TABLE_SIZE);
   int si = ci + QT_SINE_TABLE_SIZE / 4;
   si &= QT_SINE_TABLE_SIZE - 1;
   ci &= QT_SINE_TABLE_SIZE - 1;
   return qt_sine_table[si] - (qt_sine_table[ci] + 0.5 * qt_sine_table[si] * d) * d;
}

LUT 和许可证可以在这里找到：https : //code.woboq.org/qt5/qtbase/src/corelib/kernel/qmath.cpp.html#qt_sine_table

这对函数采用弧度输入。LUT 涵盖了整个 2? 输入范围。该函数使用差值在值之间进行插值d，使用余弦（再次使用正弦进行类似的插值）作为导数。

Answer 10

分享我的代码，它是一个六次多项式，没什么特别的，但重新排列以避免pows. 在 Core i7 上，这比标准实现慢 2.3 倍，尽管 [0..2*PI] 范围要快一些。对于旧处理器，这可能是标准 sin/cos 的替代方案。

/*
    On [-1000..+1000] range with 0.001 step average error is: +/- 0.000011, max error: +/- 0.000060
    On [-100..+100] range with 0.001 step average error is:   +/- 0.000009, max error: +/- 0.000034
    On [-10..+10] range with 0.001 step average error is:     +/- 0.000009, max error: +/- 0.000030
    Error distribution ensures there's no discontinuity.
*/

const double PI          = 3.141592653589793;
const double HALF_PI     = 1.570796326794897;
const double DOUBLE_PI   = 6.283185307179586;
const double SIN_CURVE_A = 0.0415896;
const double SIN_CURVE_B = 0.00129810625032;

double cos1(double x) {
    if (x < 0) {
        int q = -x / DOUBLE_PI;
        q += 1;
        double y = q * DOUBLE_PI;
        x = -(x - y);
    }
    if (x >= DOUBLE_PI) {
        int q = x / DOUBLE_PI;
        double y = q * DOUBLE_PI;
        x = x - y;
    }
    int s = 1;
    if (x >= PI) {
        s = -1;
        x -= PI;
    }
    if (x > HALF_PI) {
        x = PI - x;
        s = -s;
    }
    double z = x * x;
    double r = z * (z * (SIN_CURVE_A - SIN_CURVE_B * z) - 0.5) + 1.0;
    if (r > 1.0) r = r - 2.0;
    if (s > 0) return r;
    else return -r;
}

double sin1(double x) {
    return cos1(x - HALF_PI);
}