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

Posted

技术标签:

【中文标题】C 如何计算 sin() 和其他数学函数?【英文标题】:How does C compute sin() and other math functions? 【发布时间】:2010-02-17 22:22:49 【问题描述】:

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

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


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

【问题讨论】:

请注意这个实现依赖。您应该指定您最感兴趣的实现。 我标记了 .NET 和 C,因为我查看了这两个地方,但也无法弄清楚。尽管查看 .NET 反汇编,但它看起来可能正在调用非托管 C,据我所知,它们具有相同的实现。 另见:What algorithm is used by computers to calculate logarithms? 【参考方案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 的值,然后将其提供给 sincos

还有另一个分支可以处理 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

【讨论】:

要查看这确实是在 x86 上运行的代码:编译一个调用 sin() 的程序;输入gdb a.out,然后输入break sin,然后输入run,然后输入disassemble @Henry:不要误以为这是好代码。真的是可怕,别这样学写代码! @Andreas 嗯,你是对的,与 fdlibm 相比,IBM 代码看起来确实很糟糕。我编辑了答案以添加指向 fdlibm 正弦例程的链接。 @Henry: __kernel_sin 是在 k_sin.c 中定义的,但它是纯 C。再次单击它——我第一次搞砸了 URL。 链接的 sysdeps 代码特别有趣,因为它是正确四舍五入的。也就是说,它显然为所有输入值提供了最好的答案,而这直到最近才成为可能。在某些情况下,这可能会很慢,因为可能需要计算许多额外的数字以确保正确舍入。在其他情况下,它非常快——对于足够小的数字,答案就是角度。【参考方案2】:

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

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

【讨论】:

超越数学函数(如正弦和余弦)可以在当前 32 位桌面和服务器处理器中以微码或硬件指令的形式实现。情况并非总是如此,直到 i486(DX) 所有浮点计算都在 x86 系列的软件(“软浮点”)中完成,没有单独的协处理器。并非所有这些 (FPU) 都包含超越函数(例如 Weitek 3167)。 你能说得更具体点吗?如何使用泰勒级数“完善”近似值? 就“完善”答案而言,假设您正在计算正弦和余弦。假设您知道两者在某一点的确切值(例如来自 CORDIC),但想要附近点的值。然后对于小的差异 h,您可以应用泰勒近似 f(x + h) = f(x) + h f'(x) 或 f(x + h) = f(x) + h f'(x) + h^2 f''(x)/2. x86/x64 芯片有一个用于计算正弦 (fsin) 的汇编指令,但该指令有时非常不准确,因此很少使用。有关详细信息,请参阅randomascii.wordpress.com/2014/10/09/…。大多数其他处理器没有有正弦和余弦指令,因为在软件中计算它们提供了更大的灵活性,甚至可能更快。 intel芯片里面的cordic东西一般不用。首先,操作的准确性和分辨率对于许多应用来说都极为重要。当你到达第 7 位左右时,Cordic 是出了名的不准确,而且无法预测。其次,我听说他们的实现中有一个错误,这会导致更多问题。我看了一下linux gcc的sin函数,果然,它使用了chebyshev。不使用内置的东西。哦,还有,芯片里面的cordic算法比软件解决的慢。【参考方案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% 的方法的。

享受吧。

【讨论】:

我完全同意前几句话。此外,值得回顾的是,计算具有保证精度的特殊函数是一个难题。你提到的聪明人大部分时间都在做这件事。此外,从技术角度讲,最小-最大多项式是最受欢迎的,而切比雪夫多项式是它们的更简单代理。 -1 表示不专业和漫无边际(并且有点粗鲁)的语气,以及这个答案的实际非冗余 content 去除了漫无边际和屈尊俯就的事实,基本上归结为“他们经常使用切比雪夫多项式;更多详细信息请参阅这本书,真的很棒!”您知道,这很可能是绝对正确的,但这并不是我们在 SO 上想要的那种独立的 answer。不过,像这样浓缩起来,它会对这个问题做出不错的评论。 在早期的游戏开发年代,它通常使用对速度至关重要的查找表来完成)。我们通常不会为这些事情使用标准的 lib 函数。 我经常在嵌入式系统中使用查找表和 bittians(而不是弧度),但这是针对专门的应用程序(例如您的游戏)。我认为这家伙对 c 编译器如何计算浮点数的 sin 很感兴趣.... 啊,50 年前。我开始在迈凯轮系列的 Burroughs B220 上玩这种游戏。后来的 CDC 硬件,然后是摩托罗拉 68000。Arcsin 很乱——我选择了两个多项式的商并开发了代码来找到最佳系数。【参考方案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;


【讨论】:

如果你稍微调整一下系数(并将它们硬编码成多项式),你可以更快地停止大约 2 次迭代。 你可以用 DBL_EPSILON 替换这个神奇的 .000…01 吗?【参考方案5】:

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

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

【讨论】:

真正的实现可能不会使用泰勒级数,因为有更有效的方法。您只需要在域 [0...pi/2] 中正确逼近,并且有些函数可以比泰勒级数更有效地提供良好的逼近。 @David:我同意。我很小心,在我的回答中提到了“喜欢”这个词。但是泰勒展开是一个简单的解释近似函数的方法背后的想法。也就是说,我见过使用泰勒级数的软件实现(不确定它们是否经过优化)。 实际上,多项式逼近是计算三角函数最有效的方法之一。【参考方案6】:

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

这里是 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;

使用它,我们可以使用已经使用的一项来获得总和的新项(我们避免了阶乘和 x2p

【讨论】:

您知道您可以使用 Google Chart API 使用 TeX 制作这样的公式吗? code.google.com/apis/chart/docs/gallery/formulas.html【参考方案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 from @ 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()

【讨论】:

【参考方案8】:

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

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

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

【讨论】:

SSE2 寄存器用于计算 sin(),既不在 x86 也不在 x64 模式下,当然,无论模式如何,sin 都是在硬件中计算的。嘿,我们生活在 2010 年 :) @Igor:这取决于您正在查看的数学库。事实证明,x86 上最优化的数学库使用 sincos 的 SSE 软件实现,这比 FPU 上的硬件指令更快。更简单、更天真的库倾向于使用fsinfcos 指令。 @Stephen Canon:那些快速库是否像 FPU 寄存器那样具有 80 位精度?我非常隐秘地怀疑他们更喜欢速度而不是精度,这在许多情况下当然是合理的,例如在游戏中。而且我确实相信,使用 SSE 和预先计算的中间表计算 32 位精度的正弦值可能比使用完全精度的 FSIN 更快。如果您能告诉我那些快速库的名称,我将不胜感激,看看会很有趣。 @Igor:在 64 位模式下的 x86 上,至少在我所知道的所有类 Unix 系统上,精度限制为 64 位,而不是 x87 FPU 的 79 位。 sin() 的软件实现恰好比 fsin 计算的速度快大约两倍(正是因为它的精度较低)。请注意,众所周知,x87 的实际精度比其宣布的 79 位要低。 确实,msvc 运行时库中 sin() 的 32 位和 64 位实现都使用 FSIN 指令。事实上,它们给出了不同的结果,例如 sin(0.70444454416678126)。这将在 32 位程序中导致 0.64761068800896837(正确,公差为 0.5*(eps/2)),在 64 位程序中将导致 0.64761068800896848(错误)。【参考方案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 倍。

【讨论】:

这是对一个可能如何有效地计算 sin(x) 的一个很好的解释,但它似乎并没有真正回答 OP 的问题,特别是关于 C 的常见程度库/编译器计算它。 切比雪夫多项式最小化一个区间内的最大绝对值,但它们不会最小化目标函数和多项式之间的最大差异。 Minimax 多项式就是这样做的。【参考方案10】:

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

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

【讨论】:

【参考方案11】:

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

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

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

【讨论】:

【参考方案12】:

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

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

【讨论】:

【参考方案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 的原理得到一个小的修正添加

【讨论】:

【参考方案14】:

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

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

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

【讨论】:

【参考方案15】:

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

【讨论】:

【参考方案16】:

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

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

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

我希望这会有所帮助。

【讨论】:

【参考方案17】:

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

在 C 语言中,标准 C 库包含语言本身不包含的常用数学函数(例如,powsincos 分别表示幂、正弦和余弦)。其标头包含在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 库版本的标准调用,作为优化。

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

【讨论】:

【参考方案18】:

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

【讨论】:

这似乎并没有真正回答所问的问题。 OP 询问的是普通 C 编译器/库如何计算三角函数 (我很确定 ZX Spectrum 不符合条件),而不是如何计算它们应该。不过,这可能是对一些早期答案的有用评论 啊,你是对的。它应该是评论而不是答案。我有一段时间没有使用 SO,忘记了系统是如何工作的。无论如何,我认为 Spectrum 的实现是相关的,因为它的 CPU 非常慢,而且速度至关重要。那么最好的算法肯定还是很不错的,所以对于 C 库来说,使用切比雪夫多项式实现三角函数是一个好主意。【参考方案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));  

【讨论】:

这是一个相当糟糕的实现,因为它没有使用正弦和余弦系列的连续项具有非常简单的商。这意味着可以将乘法和除法的数量从这里的 O(n^2) 减少到 O(n)。通过减半和平方来进一步减少,例如在 bc(POSIX 多精度计算器)数学库中完成。 它似乎也没有按要求回答问题; OP 询问普通 C 编译器/库如何计算三角函数,而不是自定义重新实现。 我认为这是一个很好的答案,因为它回答了问题的精神(当然我只能猜测)对像 sin() 这样的“黑匣子”函数的好奇心。这是这里唯一的答案,它让人们有机会通过在几秒钟内掩盖它而不是阅读一些优化的 C 源代码来快速了解正在发生的事情。 事实上库使用更优化的版本,通过意识到一旦你有了一个术语,你可以通过乘以一些值来获得下一个术语。请参阅Blindy's answer 中的示例。您一次又一次地计算功率和阶乘,这要慢得多 例如Intel's ICC math library is able to do math functions in math.h much faster than Intel's hardware instructions themselves【参考方案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));

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

【讨论】:

可能是因为the machine instructions are also notoriously inaccurate.【参考方案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);

【讨论】:

难道不能只使用除法的余数而不是循环吗?类似于(对于积极的部分):x = x / PI - floor(x / PI)【参考方案22】:

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

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

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

【讨论】:

以上是关于C 如何计算 sin() 和其他数学函数?的主要内容,如果未能解决你的问题,请参考以下文章

初中数学三角函数表

数学实验:通过Matlab软件操作,对函数性质有进一步了解 f(x)=sin^2(1.006)x

与数学库 sin 相比,sin 函数不准确

事件命令信号函数方法

C库数学函数

第84天:NumPy 数学函数