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

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


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


当前回答

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

在某些情况下,最大误差不是你感兴趣的,而是最大相对误差。例如,对于正弦函数,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倍于最后一位。

其他回答

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

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

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

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

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

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

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

如果你想犯罪

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

如果你想的话,因为

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

如果你想要根号方根

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

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