一个更快但精度较低的Intel汇编fsin函数？

Question

自从 Pentium 时代以来，用于计算 x86 下 sin(x) 函数的 fsin 函数似乎甚至没有使用 SSE 寄存器，我想知道是否有更新且更好的指令集用于计算三角函数。 我通常编写 C++ 代码并进行一些 asm 优化，所以任何符合从 C++、C 到汇编线路的东西都可以为我所用。谢谢。目前我正在使用带有 gcc 和 clang（即使 clang 实际上不提供任何与 FPU 相关的优化）的 Linux 64 位系统。

编辑：

1. 我已经实现了一个 sin 函数，它通常比启用 sse 的 std::sin 快两倍。 2. 我的函数永远不会比 fsin 慢，尽管 fsin 通常更准确，但考虑到 fsin 从未超过我的 sin 实现性能，我现在将保留自己的 sin 函数，而且我的 sin 具有完全可移植性，而 fsin 仅适用于 x86。 3. 我需要这个函数进行实时计算，因此我将准确度换成速度，我认为 4-5 位小数精度足够了。 4. 不要采用基于表的方法，我不使用它，它会破坏缓存，使所有东西变慢，不要使用基于内存访问或查找表的算法。

- user2485710

---

这可能会很有用："使用英特尔SSE2指令的快速三角函数" - Alex Reinking

你需要在问题中说明这些目标：是否希望在牺牲精度的情况下获得更快的速度，否则没有人能够帮助你... - Marc Glisse

@OliCharlesworth，我喜欢你的伪随机数生成器方法，但我认为它不会像要求的那样准确到第4-5位小数。 - user2485710

1. 你想要4-5位小数的“绝对”精度还是4-5位小数的“相对”精度？
2. 输入有多灵活？它可以是一个可缩放的“int”吗？
3. 输出可以是一个可缩放的“int”吗？输入范围是-1到1弧度还是其他什么？

- chux - Reinstate Monica

绝对精度意味着结果应该在数学上正确答案的+/- 0.00001范围内。0.01弧度的正弦值将为0.01000，而1弧度的正弦值将为0.84147。相对精度意味着结果应该在数学上正确答案的+/- 0.00001倍范围内。0.01弧度的正弦值将为0.0099998，而1弧度的正弦值将为0.84147。由于您正在考虑使用定点或整数，因此似乎需要绝对精度。 - chux - Reinstate Monica

Answer 1

如果您需要一个在-π … π范围内优化绝对精度的正弦近似值，请使用：

x *（1 + x * x *（-0.1661251158026961831813227851437597220432 + x * x *（8.03943560729777481878247432892823524338e-3 + x * x * -1.4941402004593877749503989396238510717e-4））

可以使用以下方式实现：

float xx = x * x;
float s = x + (x * xx) * (-0.16612511580269618f + xx * (8.0394356072977748e-3f + xx * -1.49414020045938777495e-4f));

也许根据目标架构的特点进行优化。此外，在链接的博客文章中没有提到，如果你是用汇编实现的，一定要使用FMADD指令。如果是用C或C++实现，如果你使用了C99标准函数fmaf()，请确保生成FMADD。模拟版本比乘法和加法更昂贵，因为fmaf()所做的不完全等同于乘法后再加法（因此只是简单实现是错误的）。

在-π和π之间，sin(x)与上述多项式的区别如下图所示：

这个多项式被优化以减小它与正弦函数在-π到π之间的差异，而不仅仅是某人认为这是一个好主意。

如果你只需要[-1...1]的定义区间，那么可以通过忽略其余部分来使该多项式在该区间上更加精确。再次运行优化算法生成：

x * (1 + x * x * (-1.666659904470566774477504230733785739156e-1 + x * x *(8.329797530524482484880881032235130379746e-3 + x * x *(-1.928379009208489415662312713847811393721e-4)))

绝对误差图：

如果这对你来说太精确了，可以优化低阶多项式以实现相同的目标。然后绝对误差会更大，但你会省去一两个乘法。

- Pascal Cuoq

---

@user2485710，你的问题涉及到“sin”，所以我回答了关于“sin”的问题。不管怎样，使用的方法是Remez算法，并且已经在我的答案中提供了非常清晰的解释链接：http://lolengine.net/blog/2011/12/21/better-function-approximations 。虽然它如何工作并非必须理解才能使用它（我自己也不是很懂）。 - Pascal Cuoq

@user2485710 需要理解的是多项式逼近的原则（否则你最终会尝试用形如aX^2 + bX的多项式逼近sin函数，然后你必须到处调用abs()，这很荒谬，就像Xavier Holt答案中的“Nick's version”一样）。您还需要了解浮点数的基本知识，以便您知道将X的系数固定为1是有益的。我使用了我已经提供链接的LolRemez，但由于上述所有问题，正确使用它很复杂。 - Pascal Cuoq

1. 你在谈论什么样的多项式逼近？
2. 那个lolremez库是一段可怕的代码，我不认为它能通过任何形式的代码审查。
3. 我仍然没有任何关于Remez以外算法的参考，但是Remez仍然需要一种算法来将其转换为可以计算的东西。你一直在说这些关于多项式的事情，但你没有给出任何具体的参考，我知道那些是多项式，我知道我的目标是超越函数，问题是如何用前者计算后者。

- user2485710

@user2485710 1) http://en.wikipedia.org/wiki/Approximation_theory。有关此方面的书籍已经存在，我不会为你写一本书。2) 如果你不喜欢它，就不要使用它。在像Maple这样的工具中有可用的实现，但我没有访问这些工具的权限，也没有任何迹象表明它们的实现更加清晰。你明白这段代码不会被打包到最终产品中，对吧？3) 我已经给你提供了我所使用的所有链接，但如果你认为我使用的工具太“可怕”而拒绝使用它们，那么我无法再为你提供帮助。 - Pascal Cuoq

从你的回答中，我无法提取出任何关于Remez算法的名称、参考或算法。我已经给了你这个名字和一个开源实现的链接。“你正在从其他来源复制粘贴代码”实际上，我为你运行算法，因为它很难使用，而你的问题是关于更快但不太准确的fsin。不客气。LolRemez附带有教程，http://lolengine.net/wiki/doc/maths/remez，但你已经拒绝了那个实现，称其为“可怕的”，我也不知道是否有其他免费的实现或教程。 - Pascal Cuoq

Answer 2

如果您可以接受近似值（如果您试图击败硬件，我假设您可以接受），您应该查看Nick在DevMaster上的sin实现：

http://devmaster.net/posts/9648/fast-and-accurate-sine-cosine

他有两个版本：一个是“快速而粗糙”的方法，另一个是“慢速而准确”的方法。在回复中，有人估计相对误差分别为12％和0.2％。我自己实现了一下，在我的机器上运行时间分别为硬件时间的1/14和1/8。希望能帮到你！

PS：如果您自己做这件事，可以重构慢/准确的方法以避免乘法，并略微改进Nick的版本，但我不记得具体是如何...

- Xavier Holt

---

你可以重构那个慢而准确的方法，避免乘法并稍微改进一下比Nick版本更好。当霍纳形式优于某人的多项式评估方案时，应避免对所谓的“快速和准确”实现做出大胆的声明。这篇博客文章的标题应该是“快速而不准确的正弦”，因为两个版本都是如此。 - Pascal Cuoq

@user2485710 标题是“快速准确的正弦/余弦”，而不是“近似”。任何返回IEEE 754数字的函数都可以假定其精度受该格式的限制。当一个函数产生的结果与真实结果相差不超过1 ULP时，可以认为该函数是准确的。这篇文章所描述的是一个不准确但快速的正弦函数（来自一个从未听说过霍纳方案的人）。 - Pascal Cuoq

@PascalCuoq，你是在暗示说通过一些定点计算我们可以做得更好吗？有更好的算法吗？ - user2485710

@user2485710：请在您的问题中说明所需的定义间隔，我会向您展示一种比您已经编写的函数更适合绝对准确度的函数（我假设您对绝对精度感兴趣。您还应该让这更清楚明了）。 - Pascal Cuoq

@ user2485710 - 帕斯卡尔说帖子的标题应该包含“近似”一词，因为该方法在典型的浮点意义下 不 准确。霍纳法是发现计算机有效形式的多项式的一种方式。将其应用于此处原始代码可略微提高（5％？）速度。虽然只有轻微改进，但如果您想要速度，这绝对值得一试。 - Xavier Holt