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

有人能帮我找到他们吗?我觉得不太可能运行C语言的所有硬件都支持三角函数,所以一定有软件算法,对吧?


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


当前回答

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

芯片不使用泰勒级数来计算三角函数,至少不完全是这样。首先,他们使用CORDIC,但他们也可能使用一个短的泰勒级数来优化CORDIC的结果,或者用于特殊情况,例如在非常小的角度下以相对较高的精度计算正弦。有关更多解释,请参阅StackOverflow的回答。

其他回答

是的,也有计算罪恶的软件算法。基本上,用数字计算机计算这些东西通常是用数值方法来完成的,比如近似表示函数的泰勒级数。

数值方法可以将函数近似到任意精度,因为浮点数的精度是有限的,所以它们非常适合这些任务。

没有什么比点击源代码,看看人们是如何在常用的库中实际完成它的了;让我们特别看看一个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的原理,你会得到一个小的修正

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

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

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

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

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

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

这段代码非常复杂。没有一种软件算法在整个x值范围内尽可能快且准确,因此库实现了几种不同的算法,它的第一项工作是查看x并决定使用哪种算法。

When x is very very close to 0, sin(x) == x is the right answer. A bit further out, sin(x) uses the familiar Taylor series. However, this is only accurate near 0, so... When the angle is more than about 7°, a different algorithm is used, computing Taylor-series approximations for both sin(x) and cos(x), then using values from a precomputed table to refine the approximation. When |x| > 2, none of the above algorithms would work, so the code starts by computing some value closer to 0 that can be fed to sin or cos instead. There's yet another branch to deal with x being a NaN or infinity.

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

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

(有时)用于将x减小到接近0的值,该值与x相差π/2的倍数,特别是xn × π/2。这种没有划分或分支的方式相当聪明。但是没有任何评论!


旧的32位版本的GCC/glibc使用fsin指令,这对于某些输入是非常不准确的。有一篇精彩的博客文章用两行代码说明了这一点。

fdlibm在纯C中实现sin要比glibc简单得多,而且注释很好。源代码:fdlibm/s_sin.c和fdlibm/k_sin.c