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

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


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


当前回答

关于sin(), cos(),tan()这样的三角函数,在5年之后,没有提到高质量三角函数的一个重要方面:极差约简。

任何这些函数的早期步骤都是将角度(以弧度为单位)减小到2*π区间。但是π是无理数,所以像x =余数(x, 2*M_PI)这样的简单简化会引入误差,因为M_PI或机器pi是π的近似值。那么,如何求x =余数(x, 2*π)呢?

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

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

一个好的报告是大论点的论据缩减:好到最后一位(1992)。它涵盖了这个问题很好:讨论了需要和事情是如何在各种平台(SPARC, PC, HP, 30+其他),并提供了一个解决方案算法,为所有双从-DBL_MAX到DBL_MAX的高质量结果。


如果原始参数以度为单位,但可能值很大,则首先使用fmod()以提高精度。一个好的fmod()将不会引入任何错误,从而提供出色的范围缩小。

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

各种三角恒等式和remquo()提供了更多的改进。示例:信德()

其他回答

如果你想要一个软件实现,而不是硬件实现,可以在《数值公式》的第5章中找到这个问题的明确答案。我的副本在一个盒子里,所以我不能给出细节,但简短的版本(如果我没记错的话)是你把tan(theta/2)作为你的基本操作,然后从那里计算其他的。计算是用级数近似完成的,但它比泰勒级数收敛得快得多。

抱歉,我没拿到书就想不起来了。

库函数的实际实现取决于特定的编译器和/或库提供程序。不管它是用硬件还是软件,不管它是不是泰勒展开,等等,都会有所不同。

我意识到这完全没有帮助。

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

在某些情况下,最大误差不是你感兴趣的,而是最大相对误差。例如,对于正弦函数,x = 0附近的误差应该比较大的值小得多;你想要一个小的相对误差。所以你可以计算sinx / x的切比雪夫多项式,然后把这个多项式乘以x。

Next you have to figure out how to evaluate the polynomial. You want to evaluate it in such a way that the intermediate values are small and therefore rounding errors are small. Otherwise the rounding errors might become a lot larger than errors in the polynomial. And with functions like the sine function, if you are careless then it may be possible that the result that you calculate for sin x is greater than the result for sin y even when x < y. So careful choice of the calculation order and calculation of upper bounds for the rounding error are needed.

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

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

For example for sin (x), where you need coefficients for x, x^3, x^5, x^7 etc. you do the following: Calculate the best approximation of sin x with a polynomial (ax + bx^3 + cx^5 + dx^7) with higher than double precision, then round a to double precision, giving A. The difference between a and A would be quite large. Now calculate the best approximation of (sin x - Ax) with a polynomial (b x^3 + cx^5 + dx^7). You get different coefficients, because they adapt to the difference between a and A. Round b to double precision B. Then approximate (sin x - Ax - Bx^3) with a polynomial cx^5 + dx^7 and so on. You will get a polynomial that is almost as good as the original Chebyshev polynomial, but much better than Chebyshev rounded to double precision.

Next you should take into account the rounding errors in the choice of polynomial. You found a polynomial with minimum error in the polynomial ignoring rounding error, but you want to optimise polynomial plus rounding error. Once you have the Chebyshev polynomial, you can calculate bounds for the rounding error. Say f (x) is your function, P (x) is the polynomial, and E (x) is the rounding error. You don't want to optimise | f (x) - P (x) |, you want to optimise | f (x) - P (x) +/- E (x) |. You will get a slightly different polynomial that tries to keep the polynomial errors down where the rounding error is large, and relaxes the polynomial errors a bit where the rounding error is small.

所有这些将使您轻松地获得最多0.55倍于最后一位的舍入误差,其中+,-,*,/的舍入误差最多为0.50倍于最后一位。

关于sin(), cos(),tan()这样的三角函数,在5年之后,没有提到高质量三角函数的一个重要方面:极差约简。

任何这些函数的早期步骤都是将角度(以弧度为单位)减小到2*π区间。但是π是无理数,所以像x =余数(x, 2*M_PI)这样的简单简化会引入误差,因为M_PI或机器pi是π的近似值。那么,如何求x =余数(x, 2*π)呢?

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

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

一个好的报告是大论点的论据缩减:好到最后一位(1992)。它涵盖了这个问题很好:讨论了需要和事情是如何在各种平台(SPARC, PC, HP, 30+其他),并提供了一个解决方案算法,为所有双从-DBL_MAX到DBL_MAX的高质量结果。


如果原始参数以度为单位,但可能值很大,则首先使用fmod()以提高精度。一个好的fmod()将不会引入任何错误,从而提供出色的范围缩小。

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

各种三角恒等式和remquo()提供了更多的改进。示例:信德()

在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