278

我一直在研究 .NET 反汇编和 GCC 源代码，但似乎无法在任何地方找到实际实现sin()和其他数学函数……它们似乎总是在引用其他东西。

谁能帮我找到他们？我觉得 C 将运行的所有硬件都不太可能支持硬件中的三角函数，所以某处必须有软件算法，对吧？

我知道可以计算函数的几种方法，并且为了好玩，我已经编写了自己的例程来使用泰勒级数来计算函数。我很好奇真实的生产语言是如何做到的，因为我所有的实现总是慢几个数量级，尽管我认为我的算法非常聪明（显然他们不是）。

4

22 回答 22

248

在 GNU libm 中， 的实现sin依赖于系统。因此，您可以在sysdeps的相应子目录中找到每个平台的实现。

一个目录包括 IBM 提供的 C 语言实现。sin()自 2011 年 10 月以来，这是您在典型的 x86-64 Linux 系统上调用时实际运行的代码。它显然比fsin汇编指令快。源代码：sysdeps/ieee754/dbl-64/s_sin.c，查找__sin (double x).

这段代码非常复杂。没有一种软件算法在整个x值范围内尽可能快且准确，因此该库实现了几种不同的算法，其首要任务是查看x并决定使用哪种算法。

  • 当x非常非常接近 0 时，sin(x) == x是正确答案。

  • 稍远一点，sin(x)使用熟悉的泰勒级数。但是，这仅在 0 附近准确，所以...

  • 当角度大于约 7° 时，使用不同的算法，计算 sin(x) 和 cos(x) 的泰勒级数近似值，然后使用预先计算的表中的值来改进近似值。

  • 什么时候 | x | > 2，上述算法都不起作用，因此代码首先计算一些接近 0 的值，该值可以输入sin或cos代替。

  • 还有另一个分支来处理x是 NaN 或无穷大。

这段代码使用了一些我以前从未见过的数字技巧，尽管据我所知，它们在浮点专家中可能是众所周知的。有时几行代码需要几段来解释。例如，这两行

double t = (x * hpinv + toint);
double xn = t - toint;

（有时）用于将x减少到接近 0 的值，该值与x相差π/2 的倍数，特别是xn× π/2。不进行划分或分支的方式非常聪明。但是一点评论都没有！

较旧的 32 位版本的 GCC/glibc 使用了该fsin指令，这对于某些输入来说出人意料地不准确。有一篇引人入胜的博客文章仅用 2 行代码说明了这一点。

fdlibmsin在纯 C 中的实现比 glibc 简单得多，并且得到了很好的评论。源代码：fdlibm/s_sin.c和fdlibm/k_sin.c

于 2010-02-17T23:34:26.990 回答
72

正弦和余弦等函数在微处理器内部的微码中实现。例如，英特尔芯片有这些的组装说明。AC 编译器将生成调用这些汇编指令的代码。（相比之下，Java 编译器不会。Java 在软件而不是硬件中评估三角函数，因此运行速度要慢得多。）

芯片不使用泰勒级数来计算三角函数，至少不完全如此。首先，他们使用CORDIC，但他们也可能使用短泰勒级数来完善 CORDIC 的结果，或者用于特殊情况，例如计算非常小角度的高相对精度的正弦。有关更多说明，请参阅此StackOverflow 答案。

于 2010-02-17T22:33:45.817 回答
66

好的孩子们，该是专业人士的时间了……这是我对缺乏经验的软件工程师最大的抱怨之一。他们从头开始计算超越函数（使用泰勒级数），就好像以前没有人在他们的生活中做过这些计算一样。不对。这是一个定义明确的问题，非常聪明的软件和硬件工程师已经处理了数千次，并且有一个定义明确的解决方案。基本上，大多数超越函数都使用切比雪夫多项式来计算它们。至于使用哪些多项式取决于具体情况。首先，关于这件事的圣经是哈特和切尼的一本名为《计算机近似》的书。在那本书中，您可以决定是否有硬件加法器、乘法器、除法器等，并决定哪些运算最快。例如，如果你有一个非常快的分频器，计算正弦的最快方法可能是 P1(x)/P2(x)，其中 P1、P2 是切比雪夫多项式。如果没有快速除法器，它可能只是 P(x)，其中 P 的项比 P1 或 P2 多得多......所以它会更慢。因此，第一步是确定您的硬件及其功能。然后选择切比雪夫多项式的适当组合（例如，余弦的形式通常为 cos(ax) = aP(x)，其中 P 是切比雪夫多项式）。然后你决定你想要什么小数精度。例如，如果您想要 7 位精度，您可以在我提到的书中的相应表格中查找，它会给您（精度 = 7.33）一个数字 N = 4 和一个多项式数字 3502。N 是多项式（所以它是 p4.x^4 + p3.x^3 + p2.x^2 + p1.x + p0），因为 N=4。然后查找 p4,p3,p2,p1 的实际值，书后面的 p0 值低于 3502（它们将是浮点数）。然后你在软件中以如下形式实现你的算法： (((p4.x + p3).x + p2).x + p1).x + p0 ....这就是你如何计算余弦到十进制的 7放在那个硬件上。

请注意，FPU 中超越运算的大多数硬件实现通常涉及一些微码和类似这样的操作（取决于硬件）。切比雪夫多项式用于大多数超越数，但不是全部。例如，平方根使用牛顿拉夫森方法的双重迭代更快，首先使用查找表。同样，《计算机近似》这本书会告诉你这一点。

如果您计划实现这些功能，我会向任何人推荐他们获得该书的副本。它确实是这类算法的圣经。请注意，有许多替代方法可以计算这些值，例如 cordics 等，但这些方法往往最适合您只需要低精度的特定算法。为了保证每次的精度，切比雪夫多项式是要走的路。就像我说的，定义明确的问题。现在已经解决了 50 年......这就是它的完成方式。

现在，话虽如此，有一些技术可以使用切比雪夫多项式来获得具有低次多项式的单精度结果（如上面的余弦示例）。然后，还有其他技术可以在值之间进行插值以提高准确性，而不必使用更大的多项式，例如“Gal's Accurate Tables Method”。后一种技术就是引用 ACM 文献的帖子所指的内容。但最终，切比雪夫多项式是用来实现 90% 的方法的。

享受。

于 2013-02-14T06:55:34.380 回答
16

具体来说sin，使用泰勒展开会给你：

sin(x) := x - x^3/3！+ x^5/5！- x^7/7！+ ... (1)

您将继续添加术语，直到它们之间的差异低于可接受的公差水平或仅用于有限数量的步骤（更快，但不太精确）。一个例子是这样的：

float sin(float x)
{
  float res=0, pow=x, fact=1;
  for(int i=0; i<5; ++i)
  {
    res+=pow/fact;
    pow*=-1*x*x;
    fact*=(2*(i+1))*(2*(i+1)+1);
  }

  return res;
}

注意：(1) 之所以有效，是因为小角度的近似值 sin(x)=x。对于更大的角度，您需要计算越来越多的项才能获得可接受的结果。您可以使用 while 参数并继续以达到一定的准确性：

double sin (double x){
    int i = 1;
    double cur = x;
    double acc = 1;
    double fact= 1;
    double pow = x;
    while (fabs(acc) > .00000001 &&   i < 100){
        fact *= ((2*i)*(2*i+1));
        pow *= -1 * x*x; 
        acc =  pow / fact;
        cur += acc;
        i++;
    }
    return cur;

}
于 2010-02-17T22:33:32.593 回答
13

是的，也有用于计算的软件算法sin。基本上，用数字计算机计算这类东西通常是使用数值方法完成的，比如逼近代表函数的泰勒级数。

数值方法可以将函数逼近到任意精度，并且由于浮点数的精度是有限的，因此它们非常适合这些任务。

于 2010-02-17T22:25:16.123 回答
12

使用泰勒级数并尝试找出级数项之间的关系，这样您就不会一次又一次地计算事物

这是 cosinus 的示例：

double cosinus(double x, double prec)
{
    double t, s ;
    int p;
    p = 0;
    s = 1.0;
    t = 1.0;
    while(fabs(t/s) > prec)
    {
        p++;
        t = (-t * x * x) / ((2 * p - 1) * (2 * p));
        s += t;
    }
    return s;
}

使用它，我们可以使用已经使用的一项来获得总和的新项（我们避免了阶乘和 x 2p）

解释

于 2010-02-17T22:38:53.470 回答
12

关于三角函数，如sin(), 5cos()年tan()后没有提及高质量三角函数的一个重要方面：范围缩减。

这些函数中的任何一个的早期步骤是以弧度为单位的角度减小到 2*π 区间的范围。但是 π 是非理性的，所以简单的归约（如x = remainder(x, 2*M_PI)引入误差M_PI）或机器 pi 是 π 的近似值。那么，怎么办x = remainder(x, 2*π)？

早期的库使用扩展精度或精心设计的编程来提供高质量的结果，但仍然在有限的double. 当请求一个较大的值时sin(pow(2,30))，结果是没有意义的，或者0.0可能会设置一个错误标志，例如TLOSS完全丢失精度或PLOSS部分丢失精度。

sin()将大值的范围缩小到像 -π 到 π 这样的区间是一个具有挑战性的问题，它可以与基本 trig 函数（如本身）的挑战相媲美。

一个好的报告是Argument reduction for large arguments: Good to the last bit (1992)。它很好地涵盖了这个问题：讨论了各种平台（SPARC、PC、HP、30+ 其他）上的需求和情况，并提供了一种解决方案算法，可以为从 double.-DBL_MAX到DBL_MAX.

如果原始参数以度为单位，但可能具有很大的值，fmod()请先使用以提高精度。好的fmod()将不会引入错误，因此可以提供出色的范围缩小。

// sin(degrees2radians(x))
sin(degrees2radians(fmod(x, 360.0))); // -360.0 < fmod(x,360) < +360.0

各种触发标识并remquo()提供更多改进。示例：sind()

于 2015-07-22T18:20:50.620 回答
11

这是一个复杂的问题。x86 系列的类似 Intel 的 CPU 具有该sin()功能的硬件实现，但它是 x87 FPU 的一部分，不再用于 64 位模式（使用 SSE2 寄存器代替）。在该模式下，使用软件实现。

有几种这样的实现。一个在fdlibm中，用于 Java。据我所知，glibc 实现包含 fdlibm 的一部分，以及 IBM 贡献的其他部分。

超越函数的软件实现，例如sin()通常使用多项式的近似值，通常从泰勒级数获得。

于 2010-02-17T22:36:42.743 回答
10

正如另一个答案中提到的，切比雪夫多项式是函数与多项式之间的最大差异尽可能小的多项式。这是一个很好的开始。

在某些情况下，最大误差不是您感兴趣的，而是最大相对误差。例如，对于正弦函数，x = 0 附近的误差应该比较大的值小得多；你想要一个小的相对误差。因此，您将计算 sin x / x 的切比雪夫多项式，并将该多项式乘以 x。

接下来，您必须弄清楚如何评估多项式。您希望以这样一种方式评估它，即中间值很小，因此舍入误差很小。否则，舍入误差可能会比多项式中的误差大得多。对于像 sine 函数这样的函数，如果你不小心，那么即使 x < y，你计算 sin x 的结果也可能大于 sin y 的结果。因此需要仔细选择计算顺序和计算舍入误差的上限。

例如，sin x = x - x^3/6 + x^5 / 120 - x^7 / 5040... 如果你天真地计算 sin x = x * (1 - x^2/6 + x^4/ 120 - x^6/5040...)，那么括号中的函数是递减的，如果 y 是 x 的下一个更大的数，那么有时 sin y会小于 sin x。相反，计算 sin x = x - x^3 * (1/6 - x^2 / 120 + x^4/5040...) 这不可能发生。

例如，在计算切比雪夫多项式时，通常需要将系数四舍五入到双精度。但是，虽然切比雪夫多项式是最优的，但系数四舍五入到双精度的切比雪夫多项式并不是具有双精度系数的最优多项式！

例如对于 sin (x)，您需要 x、x^3、x^5、x^7 等的系数。您可以执行以下操作： 使用多项式 (ax + bx^3 +) 计算 sin x 的最佳近似值cx^5 + dx^7) 高于双精度，然后将a四舍五入为双精度，得到A。a和A之间的差异会很大。现在用多项式 (bx^3 + cx^5 + dx^7) 计算 (sin x - Ax) 的最佳近似值。您会得到不同的系数，因为它们适应 a 和 A 之间的差异。将 b 舍入为双精度 B。然后用多项式 cx^5 + dx^7 逼近 (sin x - Ax - Bx^3)，依此类推。您将得到一个几乎与原始切比雪夫多项式一样好的多项式，但比舍比雪夫四舍五入到双精度要好得多。

接下来，您应该考虑多项式选择中的舍入误差。您在忽略舍入误差的多项式中找到了一个误差最小的多项式，但您想优化多项式加舍入误差。一旦有了切比雪夫多项式，就可以计算舍入误差的界限。假设 f (x) 是您的函数，P (x) 是多项式，E (x) 是舍入误差。你不想优化 | f(x) - P(x)|，你要优化| f (x) - P (x) +/- E (x) |。您将得到一个略有不同的多项式，它试图在舍入误差较大的地方保持多项式误差，并在舍入误差较小的地方稍微放宽多项式误差。

所有这一切都会让您轻松地舍入误差至多为最后一位的 0.55 倍，其中 +、-、*、/ 的舍入误差至多为最后一位的 0.50 倍。

于 2015-03-20T16:42:27.140 回答
6

库函数的实际实现取决于特定的编译器和/或库提供者。无论是在硬件还是软件中完成，是否是泰勒展开式等等，都会有所不同。

我意识到这绝对没有帮助。

于 2010-02-17T23:51:55.703 回答
5

它们通常在软件中实现，并且在大多数情况下不会使用相应的硬件（即 asembly）调用。然而，正如 Jason 所指出的，这些是特定于实现的。

请注意，这些软件例程不是编译器源代码的一部分，而是可以在相应的库中找到，例如 GNU 编译器的 clib 或 glibc。见http://www.gnu.org/software/libc/manual/html_mono/libc.html#Trig-Functions

如果你想要更大的控制权，你应该仔细评估你到底需要什么。一些典型的方法是查找表的插值、汇编调用（通常很慢）或其他近似方案，例如用于平方根的 Newton-Raphson。

于 2010-02-17T22:32:27.283 回答
5

如果你想在软件而不是硬件中实现，寻找这个问题的明确答案的地方是数字食谱的第 5 章。我的副本在一个盒子里，所以我不能提供细节，但简短的版本（如果我没记错的话）是你把tan(theta/2)它作为你的原始操作并从那里计算其他操作。计算是通过级数近似完成的，但它的收敛速度比泰勒级数快得多。

抱歉，如果没有拿到这本书，我就无法记住更多。

于 2010-02-18T00:34:11.253 回答
5

没有什么比点击源代码并查看某人在常用库中实际完成它的方式更好的了。让我们特别看一个 C 库实现。我选择了uLibC。

这是 sin 函数：

http://git.uclibc.org/uClibc/tree/libm/s_sin.c

看起来它处理了一些特殊情况，然后执行一些参数缩减以将输入映射到范围 [-pi/4,pi/4]，（将参数分成两部分，大部分和尾部）打电话之前

http://git.uclibc.org/uClibc/tree/libm/k_sin.c

然后对这两个部分进行操作。如果没有尾，则使用 13 次多项式生成近似答案。如果有尾，则根据以下原理得到一个小的修正加法：sin(x+y) = sin(x) + sin'(x')y

于 2015-07-18T10:02:21.487 回答
4

每当评估这样的函数时，在某种程度上，最有可能的是：

  • 插值表（用于快速、不准确的应用程序 - 例如计算机图形）
  • 对收敛到期望值的级数的评估——可能不是泰勒级数，更可能是基于像 Clenshaw-Curtis 这样的奇特正交的东西。

如果没有硬件支持，那么编译器可能会使用后一种方法，只发出汇编代码（没有调试符号），而不是使用 ac 库——这让您很难在调试器中跟踪实际代码。

于 2010-02-17T22:41:01.703 回答
4

我将尝试回答sin()在当前 x86 处理器（假设是 Intel Core 2 Duo）上使用 GCC 的 C 编译器编译的 C 程序中的情况。

在 C 语言中，标准 C 库包括语言本身不包括的通用数学函数（例如pow，sin分别cos用于幂、正弦和余弦）。其中的标头包含在math.h中。

现在在 GNU/Linux 系统上，这些库函数由 glibc（GNU libc 或 GNU C 库）提供。但是 GCC 编译器希望您使用编译器标志链接到数学库( libm.so)-lm以启用这些数学函数的使用。我不确定为什么它不是标准 C 库的一部分。这些将是浮点函数的软件版本，或“软浮点”。

顺便说一句：据我所知，将数学函数分开的原因是历史性的，并且仅仅是为了减少非常旧的 Unix 系统中可执行程序的大小，可能在共享库可用之前。

现在编译器可以优化标准 C 库函数（sin()由较新的处理器，例如 Core 2 系列（这几乎可以追溯到 i486DX）。这将取决于传递给 gcc 编译器的优化标志。如果告诉编译器编写可以在任何 i386 或更新的处理器上执行的代码，它就不会进行这样的优化。该标志将通知编译器进行这样的优化是安全的。libm.soFSIN-mcpu=486

现在，如果程序执行 sin() 函数的软件版本，它将基于CORDIC（坐标旋转数字计算机）或BKM 算法，或者更可能是现在常用的表格或幂级数计算来计算这样的超越功能。[来源：http://en.wikipedia.org/wiki/Cordic#Application]

gcc 的任何最新版本（大约自 2.9x 以来）还提供了一个内置版本的 sin，__builtin_sin()它将用于替换对 C 库版本的标准调用，作为优化。

我敢肯定，这就像泥巴一样清晰，但希望能给你提供比你预期更多的信息，以及许多让你自己了解更多信息的起点。

于 2010-02-17T23:50:09.223 回答
4

如果您想查看这些函数在 C 中的实际 GNU 实现，请查看最新的 glibc 主干。请参阅GNU C 库。

于 2010-02-17T23:56:49.493 回答
4

正如许多人指出的那样，它依赖于实现。但据我了解您的问题，您对数学函数的真正软件实现感兴趣，但只是没有找到一个。如果是这种情况，那么您在这里：

  • 从http://ftp.gnu.org/gnu/glibc/下载 glibc 源代码
  • 查看dosincos.c位于解压缩的 glibc 根\sysdeps\ieee754\dbl-64 文件夹中的文件
  • 同样，您可以找到数学库其余部分的实现，只需查找具有适当名称的文件

您还可以查看带有.tbl扩展名的文件，它们的内容只不过是二进制形式的不同函数的预计算值的巨大表。这就是实现如此之快的原因：他们无需计算他们使用的任何系列的所有系数，而是进行快速查找，这要快得多。顺便说一句，他们确实使用 Tailor 系列来计算正弦和余弦。

我希望这有帮助。

于 2010-02-18T03:26:19.630 回答
3

不要使用泰勒级数。正如上面的几个人所指出的，切比雪夫多项式既更快又更准确。这是一个实现（最初来自 ZX Spectrum ROM）：https ://albertveli.wordpress.com/2015/01/10/zx-sine/

于 2015-07-18T08:39:57.363 回答
2

通过使用泰勒级数的代码计算正弦/余弦/正切实际上非常容易。自己写一个大概需要 5 秒。

整个过程可以用这个方程来概括：

罪恶和成本膨胀

以下是我为 C 编写的一些例程：

double _pow(double a, double b) {
    double c = 1;
    for (int i=0; i<b; i++)
        c *= a;
    return c;
}

double _fact(double x) {
    double ret = 1;
    for (int i=1; i<=x; i++) 
        ret *= i;
    return ret;
}

double _sin(double x) {
    double y = x;
    double s = -1;
    for (int i=3; i<=100; i+=2) {
        y+=s*(_pow(x,i)/_fact(i));
        s *= -1;
    }  
    return y;
}
double _cos(double x) {
    double y = 1;
    double s = -1;
    for (int i=2; i<=100; i+=2) {
        y+=s*(_pow(x,i)/_fact(i));
        s *= -1;
    }  
    return y;
}
double _tan(double x) {
     return (_sin(x)/_cos(x));  
}
于 2014-04-02T22:07:06.540 回答
0

如果你想sin那么

 __asm__ __volatile__("fsin" : "=t"(vsin) : "0"(xrads));

如果你想cos那么

 __asm__ __volatile__("fcos" : "=t"(vcos) : "0"(xrads));

如果你想sqrt那么

 __asm__ __volatile__("fsqrt" : "=t"(vsqrt) : "0"(value));

那么为什么在机器指令可以使用时使用不准确的代码呢？

于 2015-03-20T15:44:41.080 回答
0

来自 Blindy 答案的改进版代码

#define EPSILON .0000000000001
// this is smallest effective threshold, at least on my OS (WSL ubuntu 18)
// possibly because factorial part turns 0 at some point
// and it happens faster then series element turns 0;
// validation was made against sin() from <math.h>
double ft_sin(double x)
{
    int k = 2;
    double r = x;
    double acc = 1;
    double den = 1;
    double num = x;

//  precision drops rapidly when x is not close to 0
//  so move x to 0 as close as possible
    while (x > PI)
        x -= PI;
    while (x < -PI)
        x += PI;
    if (x > PI / 2)
        return (ft_sin(PI - x));
    if (x < -PI / 2)
        return (ft_sin(-PI - x));
//  not using fabs for performance reasons
    while (acc > EPSILON || acc < -EPSILON)
    {
        num *= -x * x;
        den *= k * (k + 1);
        acc = num / den;
        r += acc;
        k += 2;
    }
    return (r);
}
于 2019-10-02T13:54:02.803 回答
-1

它如何做到这一点的本质在于Gerald Wheatley的应用数值分析摘录：

当您的软件程序要求计算机获得 在此处输入图像描述or的值时在此处输入图像描述，您是否想知道如果它可以计算的最强大的函数是多项式，它如何获得这些值？它不会在表格中查找这些并进行插值！相反，计算机从一些经过剪裁以非常准确地给出值的多项式中逼近除多项式之外的所有函数。

上面要提到的几点是，一些算法实际上是从表中插值的，尽管只针对前几次迭代。还要注意它是如何提到计算机使用近似多项式而不指定哪种类型的近似多项式的。正如线程中的其他人所指出的那样，在这种情况下，切比雪夫多项式比泰勒多项式更有效。

于 2020-02-16T09:13:19.263 回答