编译器对sin()使用哪种近似算法?[复制]

在gnulibm中,实现 sin sysdeps

一个目录包含一个由IBM提供的C语言实现。从2011年10月开始,这是您调用时实际运行的代码 sin() fsin 装配说明。源代码: sysdeps/ieee754/dbl-64/s_sin.c ,寻找 __sin (double x)

这个代码非常复杂。没有一种软件算法是最快的,而且在整个范围内都是精确的 值,所以库实现了几种不同的算法,它的第一个任务是 十 并决定使用哪种算法。

• 什么时候? 非常 非常 接近0, sin(x) == x

• sin(x) 使用熟悉的泰勒级数。但是,这只精确到接近0,所以。。。

• 当角度大于约7°时,使用不同的算法,计算sin(x)和cos(x)的泰勒级数近似值,然后使用预计算表中的值来优化近似值。

• 十 |>2,以上算法都不起作用,因此代码首先计算接近0的值,然后将其输入到 罪 或 cos

• 还有一个部门要处理 十

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

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

十 接近0的值与 十 xn ×π/2。这种不使用分割或分支的方法是相当聪明的。但没有任何评论!

GCC/glibc的旧32位版本使用 fsin公司 fascinating blog post illustrating this with just 2 lines of code .

fdlibm/s_sin.c 和 fdlibm/k_sin.c

芯片 不要 用泰勒级数来计算三角函数,至少不是全部。首先他们使用 CORDIC ,但他们也可以使用一个短泰勒级数来修饰CORDIC的结果,或者用于特殊情况,例如计算非常小角度的正弦相对精度。更多的解释,请看这个 StackOverflow answer .

好了,孩子们,是时候开始专业比赛了。。。。 基本上,大多数超越函数使用切比雪夫多项式来计算它们。至于用什么样的多项式要视情况而定。首先,关于这个问题的圣经是哈特和切尼写的一本书,叫做“计算机近似法”。在那本书中,你可以决定你是否有硬件加法器、乘法器、除法器等,并决定哪些操作最快。e、 如果你有一个非常快的除法器,计算正弦的最快方法可能是P1(x)/P2(x),其中P1,P2是切比雪夫多项式。如果没有快速除法器,它可能只是P(x),其中P的项比P1或P2多得多……所以它会慢一些。所以,第一步是确定您的硬件及其功能。然后选择适当的Chebyshev多项式组合(通常是cos(ax)=aP(x)的形式来表示余弦,同样,其中P是Chebyshev多项式)。然后你决定你想要什么样的小数精度。e、 如果你想要7位数的精度,你可以在我提到的那本书的适当的表格里查一下,它会给你一个数N=4和一个多项式数3502。N是多项式的阶数(因此它是p4.x^4+p3.x^3+p2.x^2+p1.x+p0),因为N=4。然后在书的后面3502下查找p4,p3,p2,p1,p0值的实际值(它们将是浮点数)。然后在软件中以以下形式实现算法: ((p4.x+p3.x+p2.x+p1.x+p0 ……这就是你在硬件上计算余弦到小数点后7位的方法。

切比雪夫多项式用于大多数超越,但不是全部。e、 首先使用查找表的牛顿-拉斐逊法的二次迭代法的平方根更快。

如果你打算实现这些功能,我建议任何人都买一本那本书。它真的是这类算法的圣经。 请注意,有许多其他方法可以用来计算这些值,如cordics等,但这些方法往往最适合于只需要较低精度的特定算法。为了保证每次计算的精度,切比雪夫多项式是一种可行的方法。就像我说的,明确的问题。已经解决了50年了…这就是为什么。

现在,也就是说,有一些技术可以使用Chebyshev多项式来获得一个低次多项式的单精度结果(就像上面的余弦例子)。然后,还有其他技术可以在值之间插值以提高精度,而不必使用更大的多项式,例如“Gal的精确表方法”。后一种技术是指的是后一篇关于ACM文献的文章。但最终,切比雪夫多项式是用来得到90%的方法。

sin 具体地说,使用泰勒展开可以得到:

您将不断添加术语,直到它们之间的差异低于可接受的公差级别,或仅限于有限的步骤(更快,但不太精确)。例如:

float sin(float x)
{
  float res=0, pow=x, fact=1;
  for(int i=0; i<5; ++i)
  {
    res+=pow/fact;
    pow*=-1*x*x;
    fact*=(2*(i+1))*(2*(i+1)+1);
  }

  return res;
}

您可以使用while参数并继续以获得一定的精度:

double sin (double x){
    int i = 1;
    double cur = x;
    double acc = 1;
    double fact= 1;
    double pow = x;
    while (fabs(acc) > .00000001 &&   i < 100){
        fact *= ((2*i)*(2*i+1));
        pow *= -1 * x*x; 
        acc =  pow / fact;
        cur += acc;
        i++;
    }
    return cur;

}

是的,有计算的软件算法 sin numerical methods 比如说 Taylor series 表示函数。

使用 Taylor series 试着找出级数项之间的关系,这样你就不会反复计算了

double cosinus(double x, double prec)
{
    double t, s ;
    int p;
    p = 0;
    s = 1.0;
    t = 1.0;
    while(fabs(t/s) > prec)
    {
        p++;
        t = (-t * x * x) / ((2 * p - 1) * (2 * p));
        s += t;
    }
    return s;
}

^2便士 )

sin() 函数,但它是x87 FPU的一部分,不再在64位模式下使用(此处使用SSE2寄存器)。在这种模式下,使用软件实现。

有几种这样的实现。一个在里面 fdlibm

超越功能的软件实现,如 通常使用多项式近似,通常从泰勒级数中获得。

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

下一步你要知道如何计算多项式。您需要以这样一种方式来计算它:中间值很小,因此舍入误差很小。否则,舍入误差可能会比多项式中的误差大得多。对于正弦函数这样的函数,如果你不小心的话,你计算sin x的结果可能会比sin y的结果大,即使是在x<y的情况下。因此需要谨慎地选择计算顺序和计算舍入误差的上限。

将 如果y是x的下一个大数,那么sin y有时会比sin x小。取而代之,计算sinx=x-x^3*(1/6-x^2/120+x^4/5040…),而这是不可能发生的。

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

例如,对于sin(x),您需要x,x^3,x^5,x^7等的系数。您可以执行以下操作:使用一个精度高于两倍的多项式(ax+bx^3+cx^5+dx^7)计算sin x的最佳近似值,然后将a四舍五入到两倍精度,得到a。a和a之间的差会很大。现在用多项式(bx^3+cx^5+dx^7)计算(sinx-Ax)的最佳近似值。你得到了不同的系数,因为它们适应了a和a之间的差异。将b取整为双精度b,然后用多项式cx^5+dx^7近似(sinx-Ax-Bx^3),依此类推。你会得到一个多项式,它几乎和原来的切比雪夫多项式一样好,但比舍比雪夫四舍五入到两倍精度要好得多。

所有这些都可以很容易地舍入误差,最多是最后一位的0.55倍,其中+、-、*、/的舍入误差最多是最后一位的0.50倍。

sin() , cos() tan() 5年后,没有提到高质量trig功能的一个重要方面: 射程缩小 .

这些函数的早期步骤是将角度(以弧度为单位)减小到2*间隔的范围内。但是,非理性的简单削减是不是像 x = remainder(x, 2*M_PI) M_PI ,或machine pi,是的近似值。那么,怎么办 x = remainder(x, 2*π) ?

早期的库使用扩展精度或精心编制的程序来提供高质量的结果,但仍然在有限的范围内 double . 当一个大值被请求时 sin(pow(2,30)) ,结果毫无意义或 0.0 也许还有一个 error flag TLOSS 完全丧失精度或 PLOSS 部分精度损失。

罪恶() ,本身。

好报告是 Argument reduction for huge arguments: Good to the last bit 全部的 双重的 从 -DBL_MAX DBL_MAX .

如果原始参数以度为单位 ,但可能有很大的价值,使用 fmod() fmod() 将介绍 no error 因此提供了出色的距离缩小。

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

各种三角身份和 remquo() sind()

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

我知道那绝对没用。

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

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

如果你想要更大的控制力,你应该仔细评估你到底需要什么。一些典型的方法是查找表的插值、汇编调用(通常很慢)或其他近似方法,如牛顿-拉斐逊平方根法。

如果您希望在软件中而不是硬件中实现,可以在 Numerical Recipes tan(theta/2) 作为你的原始操作,从那里计算其他的。计算是用级数近似法来完成的,但它是收敛的 许多的

没有什么比访问源代码并了解人们在常用库中是如何做到这一点的;让我们具体看看一个C库实现。我选择了尤利布。

http://git.uclibc.org/uClibc/tree/libm/s_sin.c

它看起来像是处理一些特殊情况,然后执行一些参数缩减,在调用之前将输入映射到范围[-pi/4,pi/4],(将参数分成两部分,一个大的部分和一个尾部)

http://git.uclibc.org/uClibc/tree/libm/k_sin.c

然后对这两部分进行操作。 如果有一个尾部,你会得到一个基于以下原则的小修正 sin(x+y) = sin(x) + sin'(x')y

无论何时对此类函数进行评估,在某种程度上很可能存在:

• 对收敛到期望值的级数的求值,可能是 不

如果没有硬件支持,那么编译器可能会使用后一种方法,只发出汇编程序代码(没有调试符号),而不是使用c库——这使得您很难在调试器中跟踪实际代码。

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

正如许多人指出的,它依赖于实现。但据我所知,你对一个真正的 软件

• 从下载glibc源代码 http://ftp.gnu.org/gnu/glibc/
• dosincos.c 解压glibc根
• 类似地,您可以找到数学库其余部分的实现,只需查找具有适当名称的文件即可

你也可以看看文件与 .tbl 扩展,它们的内容只不过是 预计算 许多的

我希望这有帮助。

我会尽力为 sin()

pow , sin cos 功率、正弦和余弦)。其标题包含在 math.h .

math library libm.so )使用 -lm 我不知道为什么它不是标准C库的一部分。 这些将是浮点函数的软件版本,或“软浮点”。

据我所知,旧的Unix系统,可能在共享库可用之前。

现在编译器可以优化标准C库函数 (提供单位 )替换为对CPU/FPU内置sin()函数的本机指令的调用,该函数作为FPU指令存在( FSIN -mcpu=486 标志将通知编译器进行这样的优化是安全的。

现在如果程序执行sin()函数的软件版本,它将基于 CORDIC BKM algorithm ,或 更多 很可能是现在常用来计算这种超越函数的表或幂级数计算。[来源: http://en.wikipedia.org/wiki/Cordic#Application]

__builtin_sin() 它将用来代替对C库版本的标准调用,作为一种优化。

我相信这一点很清楚,但希望能给你比你期望的更多的信息,还有很多起点,让你自己学到更多。

不要用泰勒级数。切比雪夫多项式既快又准确,正如上面几个人指出的那样。以下是一个实现(最初来自ZX Spectrum ROM): https://albertveli.wordpress.com/2015/01/10/zx-sine/

整个过程可以用这个等式来概括:

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

如果你想的话 sin 然后

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

如果你想的话 cos

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

如果你想的话 sqrt 然后

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

那么为什么要使用不准确的代码呢?

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

它是如何做到这一点的精髓在于 应用数值分析 杰拉尔德·惠特利:

或 ,你想知道它是如何得到 它不会在表格中查找这些数据并进行插值!更确切地说