我一直在研究。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倍于最后一位。

其他回答

如果你想犯罪

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

如果你想的话,因为

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

如果你想要根号方根

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

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

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

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

先验函数的软件实现,如sin(),通常使用多项式逼近,通常从泰勒级数获得。

正如许多人指出的那样,它依赖于实现。但就我对你的问题的理解而言,你对数学函数的真正软件实现感兴趣,但只是没有找到一个。如果是这样的话,那么你是这样的:

从http://ftp.gnu.org/gnu/glibc/下载glibc源代码 查看位于解包的glibc根\sysdeps\ieee754\dbl-64文件夹中的文件dosincosc 类似地,您可以找到其余数学库的实现,只需查找具有适当名称的文件

您也可以看看扩展名为.tbl的文件,它们的内容只不过是以二进制形式的不同函数的预计算值的巨大表格。这就是为什么实现如此之快:而不是计算他们使用的任何级数的所有系数,他们只是做一个快速查找,这要快得多。顺便说一下,他们确实用裁缝级数来计算正弦和余弦。

我希望这能有所帮助。

计算正弦/余弦/正切其实很容易通过代码使用泰勒级数来实现。自己写一个只需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));  
}

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