59

我在谷歌上搜索了过去一个小时的问题,但只有泰勒级数或一些示例代码太慢或根本无法编译。好吧,我在谷歌上找到的大多数答案是“谷歌它,它已经被问过了”,但遗憾的是它不是......

我在低端 Pentium 4 上分析我的游戏,发现大约 85% 的执行时间浪费在计算正弦、余弦和平方根(来自 Visual Studio 中的标准 C++ 库),这似乎严重依赖 CPU(在我的 I7 上,相同的功能只获得了 5% 的执行时间,而且游戏速度更快)。我无法优化这三个函数,也无法一次性计算正弦和余弦(相互依赖),但我的模拟不需要太准确的结果,因此我可以接受更快的近似。

所以,问题是:在 C++ 中计算浮点数的正弦、余弦和平方根的最快方法是什么?

编辑 查找表更痛苦,因为在现代 CPU 上产生的 Cache Miss 比泰勒级数更昂贵。这些天 CPU 的速度实在是太快了,而缓存却不是。

我犯了一个错误,虽然我需要计算泰勒级数的几个阶乘,现在我看到它们可以实现为常数。

所以更新的问题:平方根是否也有任何快速优化?

编辑2

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

编辑3

我也不能对平方距离进行操作,我需要精确的距离进行计算

4

18 回答 18

98

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

double FastSin(double x)
{
    return 0;
}

哦,您想要比 |1.0| 更高的准确度?好吧,这是一个同样快的正弦函数:

double FastSin(double x)
{
    return x;
}

当 x 接近于零时, 这个答案实际上并不糟糕。对于较小的 x,sin(x) 约等于 x,因为 x 是 sin(x) 的泰勒展开式的第一项。

什么,对你来说还不够准确?好好读下去。

1970 年代的工程师在该领域取得了一些奇妙的发现,但新程序员根本不知道这些方法的存在,因为它们没有作为标准计算机科学课程的一部分进行教授。

您需要首先了解对于所有应用程序都没有这些功能的“完美”实现。因此,对诸如“哪个最快”之类的问题的肤浅回答肯定是错误的。

大多数提出这个问题的人都不了解性能和准确性之间权衡的重要性。特别是,在您做任何其他事情之前,您将不得不对计算的准确性做出一些选择。您可以在结果中容忍多少错误?10^-4?10^-16?

除非你可以用任何方法量化误差,否则不要使用它。 查看我下面的所有随机答案,发布一堆随机未注释的源代码,没有清楚地记录使用的算法及其在输入范围内的确切最大误差?“我猜这个错误大约是一种咕哝咕哝。” 那是严格的丛林联盟。如果您不知道如何计算PRECISE最大误差,以FULL精度,在您的近似函数中,在整个输入范围内......那么您不知道如何编写近似函数!

没有人单独使用泰勒级数来近似软件中的超越数。除了某些高度特殊的情况外,泰勒级数通常会在常见的输入范围内缓慢接近目标。

您的祖父母用来有效计算超越数的算法统称为CORDIC,并且非常简单,可以在硬件中实现。这是C 中一个有据可查的 CORDIC 实现。CORDIC 实现通常需要一个非常小的查找表,但大多数实现甚至不需要硬件乘法器可用。大多数 CORDIC 实现都允许您在性能与准确性之间进行权衡,包括我链接的那个。

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

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

我们也喜欢切比雪夫多项式,因为近似误差在输出范围内均匀分布。如果您正在编写音频插件或进行数字信号处理,切比雪夫多项式“免费”为您提供廉价且可预测的抖动效果。

如果您想在特定范围内找到自己的切比雪夫多项式系数,许多数学库将查找这些系数的过程称为“切比雪夫拟合”或类似的东西。

和现在一样,平方根通常使用牛顿-拉夫森算法的一些变体来计算,通常使用固定次数的迭代。通常,当有人想出一个“惊人的新”算法来计算平方根时,它只是变相的 Newton-Raphson。

Newton-Raphson、CORDIC 和 Chebyshev 多项式让您可以在速度与准确性之间进行权衡,因此答案可以随心所欲地不精确。

最后,当您完成所有花哨的基准测试和微优化后,请确保您的“快速”版本实际上比库版本快。 这是一个典型的 fsin() 库实现,定义 在从 -pi/4 到 pi/4 的域上。而且它并没有那么慢。

最后提醒您:您几乎肯定会使用 IEEE-754 数学来执行您的估计,并且每当您使用一堆乘法执行 IEEE-754 数学时,几十年前做出的一些晦涩的工程决策将再次困扰您你,以舍入误差的形式。这些错误一开始很小,但它们变得越来越大,越来越大,越来越大!在你生命中的某个时刻,请阅读“每个计算机科学家都应该知道的关于浮点数的知识”并有适量的恐惧。请记住,如果您开始编写自己的超越函数,则需要对浮点舍入导致的实际误差进行基准测试和测量,而不仅仅是最大理论误差。这不是理论上的问题。在多个项目中,“快速数学”编译设置让我很头疼。

tl:博士;去谷歌“正弦近似”或“余弦近似”或“平方根近似”或“近似理论”。

于 2016-06-21T14:43:44.807 回答
50

首先,泰勒级数不是实现正弦/余弦的最佳/最快方法。这也不是专业库实现这些三角函数的方式,并且了解最佳数值实现可以让您调整精度以更有效地获得速度。此外,这个问题已经在 StackOverflow 中广泛讨论过。这里只是一个例子

其次,您看到的新旧 PCS 之间的巨大差异是由于现代英特尔架构具有明确的汇编代码来计算基本三角函数。在执行速度上很难击败他们。

最后,让我们谈谈您的旧 PC 上的代码。检查gsl gnu 科学库 (或数值配方)实现,您会发现它们基本上使用 Chebyshev Approximation Formula。

Chebyshev 近似收敛速度更快,因此您需要评估的项更少。我不会在这里写实现细节,因为 StackOverflow 上已经发布了非常好的答案。例如检查这个。只需调整该系列的术语数量即可改变准确性/速度之间的平衡。

对于此类问题,如果您想了解一些特殊函数或数值方法的实现细节,您应该在采取任何进一步行动之前查看 GSL 代码 - GSL 是标准数值库。

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

EDIT2:您可以尝试制作粗糙的 sin 网格并使用 gsl 例程(gsl_interp_cspline_periodic 用于具有周期性条件的样条)来样条该表(样条将减少与线性插值相比的错误 => 你需要更少的点在你的表 = > 更少的缓存未命中)!

于 2013-09-07T00:00:00.057 回答
30

基于http://forum.devmaster.net/t/fast-and-accurate-sine-cosine/9648的想法和一些手动重写以提高微基准测试中的性能,我最终得到了以下余弦实现,即用于因在大量空间上重复 cos 调用而成为瓶颈的 HPC 物理模拟。它比查找表足够准确且速度更快,最值得注意的是不需要除法。

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,假设Tdouble通常在大多数 C++ 实现中定义。否则,相同范围的最大误差约为 0.056。

于 2015-01-20T16:24:40.007 回答
24

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

由 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;
}

这是计算平方根的常数时间

于 2013-09-06T16:46:13.590 回答
23

最快的方法是预先计算值并使用此示例中的表:

在 C++ 中创建正弦查找表

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

泰勒级数的正弦

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

使其正常工作的关键之一是预先计算阶乘并截断合理数量的项。阶乘在分母中的增长非常快,因此您不需要携带多个项。

另外...不要每次都从一开始就将您的 x^n 相乘...例如,将 x^3 乘以 x 再乘以两次,然后再乘以两次以计算指数。

于 2013-09-06T16:28:15.110 回答
8

对于 x86,硬件 FP 平方根指令很快(sqrtss是 sqrt 标量单精度)。单精度比双精度快,所以在你可以负担得起使用较少精度的代码时,一定要使用float而不是。double

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

要使 Csqrt()sqrtf()函数像sqrtsdor一样内联sqrtss您需要使用-fno-math-errno. 在 NaN 上设置数学函数errno通常被认为是设计错误,但标准要求这样做。如果没有该选项,gcc 会将其内联,然后执行比较+分支以查看结果是否为 NaN,如果是,则调用库函数以便它可以设置errno. 如果您的程序不检查errno数学函数后,使用-fno-math-errno.

您不需要-ffast-mathget的任何“不安全”部分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

用于 NaN 案例的 gcc 代码似乎过于复杂;它甚至不使用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)as时处理零的答案。x*rsqrt(x)如果您想将 FP 值转换为整数,我包含了一些关于舍入误差的讨论,以及指向其他相关问题的链接。


P4:

  • sqrtss: 23c 延迟,未流水线化
  • sqrtsd: 38c 延迟,未流水线化
  • fsqrt(x87):43c 延迟,未流水线化
  • rsqrtss/ mulss: 4c + 6c 延迟。可能每 3c 吞吐量一个,因为它们显然不需要相同的执行单元(mmx 与 fp)。

  • SIMD 打包版本有点慢


天湖:

  • 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 的实验测试的数字。请参阅他的 microarch 指南,了解是什么让代码运行得快或慢。另请参阅标签 wiki 上的链接。

IDK 如何最好地计算 sin/cos。我已经读过硬件fsin/ fcos(以及同时执行两者的唯一稍慢fsincos的)不是最快的方式,但 IDK 是什么。

于 2016-04-16T05:39:51.097 回答
6

QT 具有正弦 (qFastSin) 和余弦 (qFastCos) 的快速实现,它们使用带有插值的查找表并覆盖任何输入值(即使在 0-2PI 范围之外)。我在我的代码中使用它,它们比 std:sin/cos 快(快约 5 倍)并且足够精确以满足我的需要(与 std::sin/cos 的最大差异为 ~0.00000246408):

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π 输入范围。& 运算符提供了对周期性的修复,将 si(sin 索引)和 ci(cos 索引)值带回到 LUT 范围内。该函数使用差值在值之间进行插值d,使用余弦(再次使用正弦进行类似的插值)作为导数。

于 2018-10-16T17:37:39.787 回答
5

分享我的代码,它是一个 6 次多项式,没什么特别的,只是为了避免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);
}
于 2019-04-08T16:38:00.140 回答
5

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

但是,millianw 链接到的原始文章计算的是正弦,而不是余弦,并且它的计算方式有所不同。(看起来更简单。)

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

不幸的是,即使这篇文章包含多张图片,archive.org 也只存储了前两张。另外,作者的个人资料页面(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。

泰勒.gif

绿线是真正的正弦,红线是泰勒级数的前四项。这似乎是一个可以接受的近似值,但让我们仔细看看:taylor_zoom.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。绘制这个我们得到:抛物线.gif。这看起来比 4 项泰勒级数更糟糕,对吧?错误的!最大绝对误差为 0.056。此外,这种近似将使我们得到平滑的波动,并且只需 3 次乘法和 1 次加法即可计算!

不幸的是,它还不是很实用。这就是我们在 [-pi, pi] 范围内得到的:negative.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。

现在让我们看看余弦。基本三角学告诉我们余弦(x) = 正弦(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 为假时,& 操作将右手部分归零,因此减法不执行任何操作,这是完全等价的。我将把它作为练习留给读者为此创建工作 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 项 Taylor 更快,而且更精确!余弦版本只需要对 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 和圆周率。

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

干杯,

缺口

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

通过访问 Web 存档页面,您可能还会发现本文后面的评论很有趣:

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

于 2021-03-30T09:42:24.750 回答
4

我使用以下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);
}
于 2015-01-15T15:16:00.840 回答
4

这是之前在akellehe 的回答中给出的泰勒级数方法的实现。

unsigned int Math::SIN_LOOP = 15;
unsigned int Math::COS_LOOP = 15;

// sin(x) = x - x^3/3! + x^5/5! - x^7/7! + ...
template <class T>
T Math::sin(T x)
{
    T Sum       = 0;
    T Power     = x;
    T Sign      = 1;
    const T x2  = x * x;
    T Fact      = 1.0;
    for (unsigned int i=1; i<SIN_LOOP; i+=2)
    {
        Sum     += Sign * Power / Fact;
        Power   *= x2;
        Fact    *= (i + 1) * (i + 2);
        Sign    *= -1.0;
    }
    return Sum;
}

// cos(x) = 1 - x^2/2! + x^4/4! - x^6/6! + ...
template <class T>
T Math::cos(T x)
{
    T Sum       = x;
    T Power     = x;
    T Sign      = 1.0;
    const T x2  = x * x;
    T Fact      = 1.0;
    for (unsigned int i=3; i<COS_LOOP; i+=2)
    {
        Power   *= x2;
        Fact    *= i * (i - 1);
        Sign    *= -1.0;
        Sum     += Sign * Power / Fact;
    }
    return Sum;
}
于 2015-06-18T07:22:31.370 回答
2

所以让我换个说法,这个想法来自使用Remez 算法在区间 [-pi/4,+pi/4] 上逼近余弦和正弦函数,误差有界。然后使用范围缩小的浮点余数和输出余弦和整数商的正弦的 LUT,可以将近似值移动到任何角度参数。

它只是独一无二的,我认为可以对其进行扩展,以在有界误差方面做出更有效的算法。

void sincos_fast(float x, float *pS, float *pC){
    float cosOff4LUT[] = { 0x1.000000p+00,  0x1.6A09E6p-01,  0x0.000000p+00, -0x1.6A09E6p-01, -0x1.000000p+00, -0x1.6A09E6p-01,  0x0.000000p+00,  0x1.6A09E6p-01 };

    int     m, ms, mc;
    float   xI, xR, xR2;
    float   c, s, cy, sy;

    // Cody & Waite's range reduction Algorithm, [-pi/4, pi/4]
    xI  = floorf(x * 0x1.45F306p+00 + 0.5);              // This is 4/pi.
    xR  = (x - xI * 0x1.920000p-01) - xI*0x1.FB5444p-13; // This is pi/4 in two parts per C&W.
    m   = (int) xI;
    xR2 = xR*xR;

    // Find cosine & sine index for angle offsets indices
    mc = (  m  ) & 0x7;     // two's complement permits upper modulus for negative numbers =P
    ms = (m + 6) & 0x7;     // phase correction for sine.

    // Find cosine & sine
    cy = cosOff4LUT[mc];     // Load angle offset neighborhood cosine value 
    sy = cosOff4LUT[ms];     // Load angle offset neighborhood sine value 

    c = 0xf.ff79fp-4 + xR2 * (-0x7.e58e9p-4);               // TOL = 1.2786e-4
    // c = 0xf.ffffdp-4 + xR2 * (-0x7.ffebep-4 + xR2 * 0xa.956a9p-8);  // TOL = 1.7882e-7

    s = xR * (0xf.ffbf7p-4 + xR2 * (-0x2.a41d0cp-4));   // TOL = 4.835251e-6
    // s = xR * (0xf.fffffp-4 + xR2 * (-0x2.aaa65cp-4 + xR2 * 0x2.1ea25p-8));  // TOL = 1.1841e-8

    *pC = c*cy - s*sy;      
    *pS = c*sy + s*cy;
}

float sqrt_fast(float x){
    union {float f; int i; } X, Y;
    float ScOff;
    uint8_t e;

    X.f = x;
    e = (X.i >> 23);           // f.SFPbits.e;

    if(x <= 0) return(0.0f);

    ScOff = ((e & 1) != 0) ? 1.0f : 0x1.6a09e6p0;  // NOTE: If exp=EVEN, b/c (exp-127) a (EVEN - ODD) := ODD; but a (ODD - ODD) := EVEN!!

    e = ((e + 127) >> 1);                            // NOTE: If exp=ODD,  b/c (exp-127) then flr((exp-127)/2)
    X.i = (X.i & ((1uL << 23) - 1)) | (0x7F << 23);  // Mask mantissa, force exponent to zero.
    Y.i = (((uint32_t) e) << 23);

    // Error grows with square root of the exponent. Unfortunately no work around like inverse square root... :(
    // Y.f *= ScOff * (0x9.5f61ap-4 + X.f*(0x6.a09e68p-4));        // Error = +-1.78e-2 * 2^(flr(log2(x)/2))
    // Y.f *= ScOff * (0x7.2181d8p-4 + X.f*(0xa.05406p-4 + X.f*(-0x1.23a14cp-4)));      // Error = +-7.64e-5 * 2^(flr(log2(x)/2))
    // Y.f *= ScOff * (0x5.f10e7p-4 + X.f*(0xc.8f2p-4 +X.f*(-0x2.e41a4cp-4 + X.f*(0x6.441e6p-8))));     // Error =  8.21e-5 * 2^(flr(log2(x)/2))
    // Y.f *= ScOff * (0x5.32eb88p-4 + X.f*(0xe.abbf5p-4 + X.f*(-0x5.18ee2p-4 + X.f*(0x1.655efp-4 + X.f*(-0x2.b11518p-8)))));   // Error = +-9.92e-6 * 2^(flr(log2(x)/2))
    // Y.f *= ScOff * (0x4.adde5p-4 + X.f*(0x1.08448cp0 + X.f*(-0x7.ae1248p-4 + X.f*(0x3.2cf7a8p-4 + X.f*(-0xc.5c1e2p-8 + X.f*(0x1.4b6dp-8))))));   // Error = +-1.38e-6 * 2^(flr(log2(x)/2))
    // Y.f *= ScOff * (0x4.4a17fp-4 + X.f*(0x1.22d44p0 + X.f*(-0xa.972e8p-4 + X.f*(0x5.dd53fp-4 + X.f*(-0x2.273c08p-4 + X.f*(0x7.466cb8p-8 + X.f*(-0xa.ac00ep-12)))))));    // Error = +-2.9e-7 * 2^(flr(log2(x)/2))
    Y.f *= ScOff * (0x3.fbb3e8p-4 + X.f*(0x1.3b2a3cp0 + X.f*(-0xd.cbb39p-4 + X.f*(0x9.9444ep-4 + X.f*(-0x4.b5ea38p-4 + X.f*(0x1.802f9ep-4 + X.f*(-0x4.6f0adp-8 + X.f*(0x5.c24a28p-12 ))))))));   // Error = +-2.7e-6 * 2^(flr(log2(x)/2))

    return(Y.f);
}

较长的表达式更长、更慢,但更精确。多项式是根据霍纳规则编写的

于 2017-04-03T04:29:35.333 回答
2

这是一个非常快的 sinus 实现,它的工作原理如下:

它具有平方根复数的算术实现

从复数的分析数学中,您知道当复数是平方根时,角度是一半

您可以取一个您已经知道其角度的复数(例如 i,角度为 90 度或 PI / 2 弧度)

比通过平方根可以得到复数形式 cos (90 / 2^n) + i sin (90 / 2^n)

从复数的分析数学中,您知道当两个数字相乘时,它们的角度相加

您可以将数字 k(作为 sin 或 cos 中的参数)显示为角度 90 / 2^n 的总和,然后通过将您预先计算的复数相乘得到结果值

结果将采用 cos k + i sin k 的形式

#define PI 3.14159265
#define complex pair <float, float>

/* this is square root function, uses binary search and halves mantisa */

float sqrt(float a) {

    float b = a;

    int *x = (int*) (&b); // here I get integer pointer to float b which allows me to directly change bits from float reperesentation

    int c = ((*x >> 23) & 255) - 127; // here I get mantisa, that is exponent of 2 (floats are like scientific notation 1.111010101... * 2^n)

    if(c < 0) c = -((-c) >> 1); // ---
                                //   |--> This is for halfing the mantisa
    else c >>= 1;               // ---

    *x &= ~(255 << 23); // here space reserved for mantisa is filled with 0s

    *x |= (c + 127) << 23; // here new mantisa is put in place

    for(int i = 0; i < 5; i++) b = (b + a / b) / 2; // here normal square root approximation runs 5 times (I assume even 2 or 3 would be enough)

    return b;
}

/* this is a square root for complex numbers (I derived it in paper), you'll need it later */

complex croot(complex x) {

    float c = x.first, d = x.second;

    return make_pair(sqrt((c + sqrt(c * c + d * d)) / 2), sqrt((-c + sqrt(c * c + d * d)) / 2) * (d < 0 ? -1 : 1));
}

/* this is for multiplying complex numbers, you'll also need it later */

complex mul(complex x, complex y) {

    float a = x.first, b = x.second, c = y.first, d = y.second;

    return make_pair(a * c - b * d, a * d + b * c);
}

/* this function calculates both sinus and cosinus */

complex roots[24];

float angles[24];

void init() {

    complex c = make_pair(-1, 0); // first number is going to be -1

    float alpha = PI; // angle of -1 is PI

    for(int i = 0; i < 24; i++) {

        roots[i] = c; // save current c

        angles[i] = alpha; // save current angle

        c = croot(c); // root c

        alpha *= 0.5; // halve alpha
    }
}

complex cosin(float k) {

    complex r = make_pair(1, 0); // at start 1

    for(int i = 0; i < 24; i++) {

        if(k >= angles[i]) { // if current k is bigger than angle of c

            k -= angles[i]; // reduce k by that number

            r = mul(r, roots[i]); // multiply the result by c
        }
    }

    return r; // here you'll have a complex number equal to cos k + i sin k.
}

float sin(float k) {

    return cosin(k).second;
}

float cos(float k) {

    return cosin(k).first;
}

现在如果你仍然觉得它很慢,你可以减少函数中的迭代次数cosin(注意精度会降低)

于 2020-02-21T16:56:27.333 回答
2

此公式给出了保留 90 度倍数的导数的正弦函数的近似值。推导类似于Bhaskara I 的正弦逼近公式,但约束是将 0、90 和 180 度的值和导数设置为正弦函数的值和导数。如果您需要该功能在任何地方都平滑,您可以使用它。

#define PI 3.141592653589793

double fast_sin(double x) {
    x /= 2 * PI;
    x -= (int) x;

    if (x <= 0.5) {
        double t = 2 * x * (2 * x - 1);
        return (PI * t) / ((PI - 4) * t - 1);
    }
    else {
        double t = 2 * (1 - x) * (1 - 2 * x);
        return -(PI * t) / ((PI - 4) * t - 1);
    }
}

double fast_cos(double x) {
    return fast_sin(x + 0.5 * PI);
}

至于它的速度,它std::sin()每次调用至少比函数平均快 0.3 微秒。最大绝对误差为 0.0051。

于 2020-06-17T09:04:03.820 回答
1

我没有测试过这有多快,但我觉得它仍然很有趣。

对于xin [0,1],将 的零度泰勒近似sqrt((2-x)/3)Hardy's (1959)的余弦近似相结合给出了一个好的近似:

cos(x*pi/2) = 1 - x^2 / (0.81649658 + 0.183503 x)

其中分母中的系数来自sqrt(2/3)1 - sqrt(2/3)

如果你愿意拉伸到平方根的一级近似值,那么你会得到一个非常好的近似值

cos(x*pi/2) = 1 - x^2 / (0.816497 - 0.0206207 x + 0.204124 x^2)

sqrt(2/3)其中分母中的系数是1-sqrt(2/3)-1/(2sqrt(6))1/(2sqrt(6))

等等...

这是上面二次版本的图像:

在此处输入图像描述

于 2021-07-05T10:15:19.960 回答
0

只需将 FPU 与内联 x86 一起用于 Wintel 应用程序。据报道,直接 CPU sqrt 函数在速度上仍优于任何其他算法。我的自定义 x86 数学库代码适用于标准 MSVC++ 2005 及更高版本。如果您想要我介绍的更高精度,则需要单独的浮点/双精度版本。有时编译器的“__inline”策略会出错,所以为了安全起见,您可以将其删除。有了经验,您可以切换到宏来完全避免每次调用函数。

extern __inline float  __fastcall fs_sin(float x);
extern __inline double __fastcall fs_Sin(double x);
extern __inline float  __fastcall fs_cos(float x);
extern __inline double __fastcall fs_Cos(double x);
extern __inline float  __fastcall fs_atan(float x);
extern __inline double __fastcall fs_Atan(double x);
extern __inline float  __fastcall fs_sqrt(float x);
extern __inline double __fastcall fs_Sqrt(double x);
extern __inline float  __fastcall fs_log(float x);
extern __inline double __fastcall fs_Log(double x);

extern __inline float __fastcall fs_sqrt(float x) { __asm {
FLD x  ;// Load/Push input value
FSQRT
}}

extern __inline double __fastcall fs_Sqrt(double x) { __asm {
FLD x  ;// Load/Push input value
FSQRT
}}

extern __inline float __fastcall fs_sin(float x) { __asm {
FLD x  ;// Load/Push input value
FSIN
}}

extern __inline double __fastcall fs_Sin(double x) { __asm {
FLD x  ;// Load/Push input value
FSIN
}}    

extern __inline float __fastcall fs_cos(float x) { __asm {
FLD x  ;// Load/Push input value
FCOS
}}

extern __inline double __fastcall fs_Cos(double x) { __asm {
FLD x  ;// Load/Push input value
FCOS
}}

extern __inline float __fastcall fs_tan(float x) { __asm {
FLD x  ;// Load/Push input value
FPTAN
}}

extern __inline double __fastcall fs_Tan(double x) { __asm {
FLD x  ;// Load/Push input value
FPTAN
}}

extern __inline float __fastcall fs_log(float x) { __asm {
FLDLN2
FLD x
FYL2X
FSTP ST(1) ;// Pop1, Pop2 occurs on return
}}

extern __inline double __fastcall fs_Log(double x) { __asm {
FLDLN2
FLD x
FYL2X
FSTP ST(1) ;// Pop1, Pop2 occurs on return
}}
于 2019-07-25T19:08:36.850 回答
0

这是一个可能的加速,这在很大程度上取决于您的应用程序。(您的程序可能根本无法使用它,但我将其发布在这里,因为它可能。)我也只是在这里发布数学,代码由您决定。

对于我的应用程序,我需要计算围绕一个完整圆的每个角度步长 (dA) 的正弦和余弦。

这样我就可以利用一些三角身份:

cos(-A) = cos(A)

sin(-A) = -sin(A)

这样就可以了,所以我只需要计算半圈的 sin 和 cos。

我还设置了指向我的输出数组的指针,这也加快了我的计算速度。我不确定这一点,但我相信我的编译器矢量化了我的计算。

第三个是使用:

sin(A+dA) = sin(a)*cos(dA) + cos(a)*sin(dA)

cos(a+dA) = cos(a)*cos(dA) - sin(a)*sin(dA)

这使得我只需要实际计算一个角度的正弦和余弦 - 其余的都是用两个乘法和一个加法计算的。(这需要注意的是,在计算 sin(dA) 和 cos(dA) 时,舍入误差可能会在您绕到圆的一半时累积。再一次,如果您使用它,您的应用程序就是一切。)

于 2020-12-08T17:42:12.443 回答
-1

超过 100000000 次测试,milianw 的答案比 std::cos 实现慢 2 倍。但是,您可以通过执行以下步骤来更快地运行它:

-> 使用浮动

-> 不要使用地板,而是使用 static_cast

-> 不要使用 abs 而是三元条件

-> 使用#define 常量进行除法

-> 使用宏来避免函数调用

// 1 / (2 * PI)
#define FPII 0.159154943091895
//PI / 2
#define PI2 1.570796326794896619

#define _cos(x)         x *= FPII;\
                        x -= .25f + static_cast<int>(x + .25f) - 1;\
                        x *= 16.f * ((x >= 0 ? x : -x) - .5f);
#define _sin(x)         x -= PI2; _cos(x);

对 std::cos 和 _cos(x) 的调用超过 100000000 次,std::cos 运行时间约为 14 秒,而 _cos(x) 则运行时间约为 3 秒(_sin(x) 运行时间要多一点)

于 2016-05-11T07:26:11.927 回答