C 如何计算 sin() 和其他数学函数？

Question

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

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

---

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

Answer 1

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

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

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

• 当 x 非常 very 接近 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 指令，这对于某些输入来说出乎意料地不准确。有一个fascinating blog post illustrating this with just 2 lines of code。

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

Answer 2

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

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

Answer 3

好孩子们，该是专业人士的时间了...... 这是我对缺乏经验的软件工程师最大的抱怨之一。他们从头开始计算超越函数（使用泰勒级数），就好像以前没有人在他们的生活中做过这些计算一样。不对。这是一个定义明确的问题，非常聪明的软件和硬件工程师已经处理了数千次，并且有一个定义明确的解决方案。 基本上，大多数超越函数都使用切比雪夫多项式来计算它们。至于使用哪些多项式取决于具体情况。首先，关于这件事的圣经是哈特和切尼的一本名为《计算机近似》的书。在那本书中，您可以决定是否有硬件加法器、乘法器、除法器等，并决定哪些运算最快。例如如果你有一个非常快的除法器，计算正弦的最快方法可能是 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。然后你在书后面的 3502 下查找 p4,p3,p2,p1,p0 值的实际值（它们将是浮点数）。然后，您以以下形式在软件中实现您的算法： (((p4.x + p3).x + p2).x + p1).x + p0 ....这就是您在该硬件上计算余弦到小数点后 7 位的方式。

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

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

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

享受吧。

Answer 4

特别是对于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;

}

Answer 5

是的，也有计算sin 的软件算法。基本上，用数字计算机计算这类东西通常是使用numerical methods 来完成的，比如近似代表函数的Taylor series。

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

Answer 6

使用Taylor series 并尝试找到系列项之间的关系，这样您就不会一次又一次地计算事情

这是余弦的一个例子：

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）

Answer 7

关于像sin()、cos()、tan() 这样的三角函数，5 年后没有提到高质量三角函数的一个重要方面：范围缩小。 p>

任何这些函数的早​​期步骤是将角度（以弧度为单位）减小到 2*π 区间的范围内。但是 π 是不合理的，所以像 x = remainder(x, 2*M_PI) 这样的简单归约会引入错误，因为 M_PI 或机器 pi 是 π 的近似值。那么x = remainder(x, 2*π)怎么办？

早期的库使用扩展精度或精心设计的编程来提供高质量的结果，但仍超出double 的有限范围。当请求像sin(pow(2,30)) 这样的大值时，结果毫无意义或0.0 并且可能将error flag 设置为TLOSS 完全丢失精度或PLOSS 部分精度丢失。

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

一个好的报告是Argument reduction for huge arguments: Good to the last bit (1992)。它很好地涵盖了这个问题：讨论了各种平台（SPARC、PC、HP、30 多个其他平台）上的需求和情况，并提供了一种解决方案算法，可以为 all double 提供高质量的结果来自 @ 987654338@转DBL_MAX。

---

如果原始参数以度为单位，但可能值很大，请先使用fmod() 以提高精度。一个好的fmod() 将引入no error 并因此提供出色的范围缩小。

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

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

Answer 8

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

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

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

Answer 9

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

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

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

例如，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...) 这不可能发生。

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

例如对于 sin (x)，您需要 x、x^3、x^5、x^7 等的系数。您可以执行以下操作：使用多项式 (ax + bx) 计算 sin x 的最佳近似值^3 + cx^5 + dx^7) 高于双精度，然后将a四舍五入为双精度，得到A。a和A之间的差异会很大。现在用多项式 (b x^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 倍。

Answer 10

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

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

Answer 11

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

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

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

Answer 12

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

很抱歉，如果没有拿到这本书，我就记不住了。

Answer 13

没有什么比点击源代码并查看某人在常用库中实际完成它的方式更好的了；让我们特别看一个 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 的原理得到一个小的修正添加

Answer 14

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

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

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

Answer 15

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

Answer 16

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

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

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

我希望这会有所帮助。

Answer 17

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

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

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

旁白：将数学函数分开的原因是历史性的，只是为了减少非常旧 Unix 系统中可执行程序的大小，可能在共享之前据我所知，图书馆是可用的。

现在编译器可以优化标准 C 库函数 sin()（由 libm.so 提供）替换为对 CPU/FPU 内置 sin() 函数的本机指令的调用，该函数作为FPU 指令（FSIN for x86/x87）在 Core 2 系列等较新的处理器上（这几乎可以追溯到 i486DX）。这将取决于传递给 gcc 编译器的优化标志。如果编译器被告知编写可以在任何 i386 或更新的处理器上执行的代码，它就不会进行这样的优化。 -mcpu=486 标志将通知编译器进行此类优化是安全的。

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

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

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

Answer 18

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

Answer 19

计算正弦/余弦/正切实际上很容易通过使用泰勒级数的代码来完成。自己写一个大概需要 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));  
}

Answer 20

如果你想要sin 那么

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

如果你想要cos 那么

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

如果你想要sqrt 那么

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

既然机器指令可以使用，为什么还要使用不准确的代码呢？

Answer 21

来自 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);
}

Answer 22

它如何做到这一点的本质在于杰拉德·惠特利 (Gerald Wheatley) 的 Applied Numerical Analysis 摘录：

当您的软件程序要求计算机获取 或 ，你有没有想过它是如何获得 如果它可以计算的最强大的函数是多项式的值？ 它不会在表格中查找这些并进行插值！而是， 计算机从某些函数中逼近除多项式以外的所有函数 为非常准确地给出值而量身定制的多项式。

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