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

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


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


当前回答

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

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

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

其他回答

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

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

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

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

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

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

它如何做到这一点的本质在于杰拉德·惠特利的《应用数值分析》节选:

当你的软件程序要求计算机获取一个值时 或者,你有没有想过它是如何得到 如果它能计算的最强大的函数是多项式? 它不会在表中查找这些并进行插值!相反, 计算机逼近除多项式以外的所有函数 一个多项式,可以精确地给出值。

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

盲汉回答的改进版代码

#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);
}