C如何计算sin()和其他数学函数？

Question

我一直在研究.NET反汇编和GCC源代码,但似乎无法找到实际实现sin()和其他数学函数的任何地方......它们似乎总是引用其他东西.

谁能帮我找到它们？我觉得C运行的所有硬件都不太可能支持硬件中的触发功能,因此某处必须有软件算法,对吧？

---

我知道有几种方法可以计算函数,并编写了我自己的例程来计算函数使用泰勒系列来获得乐趣.我很好奇真正的生产语言是如何做到的,因为我的所有实现总是慢几个数量级,即使我认为我的算法非常聪明(显然它们不是).

Answer 1

在GNU libm中,实现sin依赖于系统.因此,您可以在sysdeps的相应子目录中的某个位置找到每个平台的实现.

一个目录包含由IBM提供的C实现.自2011年10月以来,这是在您调用sin()典型的x86-64 Linux系统时实际运行的代码.它显然比fsin汇编指令快.源代码:sysdeps/ieee754/dbl-64/s_sin.c,寻找__sin (double x).

这段代码非常复杂.没有一种软件算法尽可能快,并且在整个x值范围内也是准确的,因此库实现了许多不同的算法,它的第一项工作是查看x并决定使用哪种算法.在某些地区,它使用的是看似熟悉的泰勒系列.一些算法首先计算快速结果,然后如果不够准确,则丢弃它并回退到较慢的算法.

较旧的32位版本的GCC/glibc使用了该fsin指令,这对于某些输入而言令人惊讶地不准确.有一个引人入胜的博客文章只用2行代码就可以说明这一点.

fdlibm sin在纯C中的实现比glibc简单得多,并且得到很好的评论.源代码:fdlibm/s_sin.c和fdlibm/k_sin.c

Answer 2

正弦和余弦等函数在微处理器内部的微码中实现.例如,英特尔芯片有这些的组装说明.AC编译器将生成调用这些汇编指令的代码.(相比之下,Java编译器不会.Java评估软件中的trig函数而不是硬件,因此运行速度要慢得多.)

芯片不使用泰勒级数计算三角函数,至少不完全.首先,他们使用CORDIC,但他们也可能使用短的泰勒系列来完善CORDIC的结果或特殊情况,例如计算正弦,对于非常小的角度具有高相对精度.有关更多说明,请参阅此StackOverflow答案.

Answer 3

好孩子,专业人士的时间......这是我对没有经验的软件工程师的最大抱怨之一.他们从头开始计算超越函数(使用泰勒系列),好像在他们的生活中没有人曾经做过这些计算.不对.这是一个定义明确的问题,非常聪明的软件和硬件工程师已经接触过数千次,并且有一个定义明确的解决方案.基本上,大多数超越函数使用切比雪夫多项式来计算它们.至于使用哪种多项式取决于具体情况.首先,关于这件事的圣经是哈特和切尼的一本名为"计算机近似"的书.在该书中,您可以决定是否有硬件加法器,乘法器,除法器等,并确定哪些操作最快.例如,如果你有一个非常快的分频器,计算正弦的最快方法可能是P1(x)/ P2(x),其中P1,P2是切比雪夫多项式.没有快速分频器,它可能只是P(x),其中P具有比P1或P2更多的项......所以它会更慢.因此,第一步是确定您的硬件及其功能.然后选择Chebyshev多项式的适当组合(例如,对于余弦通常具有cos(ax)= aP(x)的形式,其中P是Chebyshev多项式).然后你决定你想要的小数精度.例如,如果你想要7位精度,你可以在我提到的书中的相应表中查找,它会给你(精度= 7.33)一个数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的十进制数硬件上的位置.

请注意,FPU中的超越操作的大多数硬件实现通常涉及一些微代码和这样的操作(取决于硬件).切比雪夫多项式用于大多数超越,但不是全部.例如,Square root使用查找表首先使用Newton raphson方法的双重迭代更快.同样,那本书"计算机近似"会告诉你.

如果您计划对这些功能进行推理,我建议任何人都能获得该书的副本.它确实是这类算法的圣经.请注意,有许多替代方法可用于计算这些值,如cordics等,但这些方法往往最适合您只需要低精度的特定算法.为了保证每次精度,切比雪夫多项式是最佳选择.就像我说的,明确的问题.已经解决了50年了......那就是它如何完成的.

现在,正如所说的那样,有些技术可以使用切比雪夫多项式来获得具有低次多项式的单精度结果(如上面的余弦的例子).然后,还有其他技术可以在值之间进行插值以提高精度,而无需使用更大的多项式,例如"Gal的精确表格方法".后一种技术是指ACM文献所指的内容.但最终,Chebyshev多项式是用来获得90%的方式.

请享用.

Answer 4

对于sin具体而言,利用泰勒展开会给你:

sin(x):= x - x ^ 3/3!+ x ^ 5/5! - x ^ 7/7!+ ...(1)

您将继续添加术语,直到它们之间的差异低于可接受的容差水平或仅仅是有限的步骤(更快,但不太精确).一个例子是:

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

注意:(1)因为小角度的近似sin(x)= x而起作用.对于更大的角度,您需要计算越来越多的术语以获得可接受的结果.您可以使用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;

}

Answer 5

使用泰勒系列并尝试找出系列术语之间的关系,这样你就不会一次又一次地计算事物

这是cosinus的一个例子:

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

使用这个我们可以使用已经使用过的一个得到总和的新项(我们避免使用阶乘和x ^ 2p)

解释http://img514.imageshack.us/img514/1959/82098830.jpg

Answer 6

是的,还有用于计算的软件算法sin.基本上,使用数字计算机计算这些类型的东西通常使用数值方法来完成,例如近似表示函数的泰勒级数.

数值方法可以将函数逼近任意精度,并且由于浮点数中的精度是有限的,因此它们非常适合这些任务.

Answer 7

这是一个复杂的问题.类似Intel的x86系列CPU具有该sin()功能的硬件实现,但它是x87 FPU的一部分,在64位模式下不再使用(其中使用SSE2寄存器).在该模式中,使用软件实现.

那里有几个这样的实现.一个是在fdlibm中,用于Java.据我所知,glibc实现包含部分fdlibm,以及IBM提供的其他部分.

超常函数的软件实现,例如sin()通常使用多项式的近似,通常从泰勒级数获得.

Answer 8

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

在某些情况下,最大误差不是您感兴趣的,而是最大相对误差.例如,对于正弦函数,x = 0附近的误差应该比较大的值小得多; 你想要一个小的相对错误.因此,您将计算sin x/x的Chebyshev多项式,并将该多项式乘以x.

接下来,您必须弄清楚如何评估多项式.您希望以中间值较小的方式对其进行评估,因此舍入误差很小.否则,舍入误差可能比多项式中的误差大得多.并且使用正弦函数之类的函数,如果你不小心,那么即使x <y,你为sin x计算的结果也可能大于sin y的结果.因此需要仔细选择计算顺序和计算舍入误差的上限.

例如,sin x = x - x ^ 3/6 + x ^ 5/120 - x ^ 7/5040 ......如果你计算天真sin x = x*(1 - x ^ 2/6 + x ^ 4/120 - X ^五千零四十零分之六...),则在括号中的那功能被降低,并且将发生的是,如果y是一个较大数为x,则有时罪ý将比罪x小.相反,计算sin x = x - x ^ 3*(1/6 - x ^ 2/120 + x ^ 4/5040 ......),这是不可能的.

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

例如,对于sin(x),您需要x,x ^ 3,x ^ 5,x ^ 7等系数,您可以执行以下操作:使用多项式计算sin x的最佳近似值(ax + bx ^ 3 + cx ^ 5 + dx ^ 7)高于双精度,然后将a舍入为双精度,得到A. a和A之间的差异将非常大.现在用多项式(bx ^ 3 + cx ^ 5 + dx ^ 7)计算(sin x-Ax)的最佳近似值.你得到不同的系数,因为它们适应a和A之间的差异.圆b到双精度B.然后用多项式cx ^ 5 + dx ^ 7近似(sin x - Ax - Bx ^ 3),依此类推.你将得到一个几乎与原始Chebyshev多项式一样好的多项式,但比Chebyshev更好地舍入到双精度.

接下来,您应该考虑多项式选择中的舍入误差.您在多项式中找到了一个忽略舍入误差的最小误差的多项式,但是您希望优化多项式加舍入误差.一旦有了Chebyshev多项式,就可以计算出舍入误差的界限.假设f(x)是你的函数,P(x)是多项式,E(x)是舍入误差.你不想优化| f(x) - P(x)|,你想优化| f(x) - P(x)+/- E(x)|.您将获得稍微不同的多项式,该多项式试图在舍入误差较大的地方保持多项式误差,并且在舍入误差较小的情况下放松多项式误差.

所有这些都可以让您轻松舍入最后一位最多0.55倍的错误,其中+, - ,*,/的舍入误差最多为最后一位的0.50倍.

Answer 9

关于三角函数像sin(),cos(),tan()一直没有提到,5年后,高品质的三角函数的一个重要方面:范围减少.

任何这些函数的早​​期步骤是以弧度将角度减小到2*π区间的范围.但π是不合理的,因此x = remainder(x, 2*M_PI)引入误差M_PI或机器pi 等简单的约简是π的近似值.那么,该怎么办x = remainder(x, 2*?)？

早期的图书馆使用扩展精度或精心编程来提供高质量的结果,但仍然在有限的范围内double.当请求大值时sin(pow(2,30)),结果没有意义,或者0.0可能将错误标志设置为TLOSS完全丢失精度或PLOSS部分精度损失.

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

一个好的报告是对于大论证的论证减少:好到最后一点(1992).它涵盖了问题得好:讨论的需要,事情是如何在各种平台(SPARC,PC,HP,其他30+),并提供了一个解决方案算法给出了质量的结果都 double来自-DBL_MAX于DBL_MAX.

如果原始参数以度为单位,但可能具有较大值,请fmod()首先使用以提高精度.好的fmod()将不会引入错误,因此可以提供出色的范围缩小.

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

各种触发身份并remquo()提供更多改进.示例:sind()

Answer 10

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

我意识到这绝对没有帮助.

Answer 11

它们通常以软件实现,并且在大多数情况下不会使用相应的硬件(即,aseembly)调用.但是,正如Jason指出的那样,这些都是特定于实现的.

请注意,这些软件例程不是编译器源的一部分,而是可以在诸如clib的对应库中找到,或者用于GNU编译器的glibc.请参阅http://www.gnu.org/software/libc/manual/html_mono/libc.html#Trig-Functions

如果您想要更好的控制,您应该仔细评估您的需求.一些典型的方法是查找表的插值,汇编调用(通常很慢),或其他近似方案,如平方根的Newton-Raphson.

Answer 12

我将尝试回答sin()一个C程序的情况,在当前的x86处理器上用GCC的C编译器编译(假设是Intel Core 2 Duo).

在C语言中的标准C库包括普通的数学函数,不包括在语言本身(例如pow,sin与cos用于电源,正弦,余弦及分别).其标题包含在math.h中.

现在在GNU/Linux系统上,这些库函数由glibc(GNU libc或GNU C Library)提供.但是GCC编译器希望您使用编译器标志链接到数学库(libm.so)-lm以启用这些数学函数的使用.我不确定为什么它不是标准C库的一部分.这些将是浮点函数的软件版本,或"软浮点".

旁白:将数学函数分开的原因是历史性的,并且仅仅是为了减少非常旧的Unix系统中可执行程序的大小,可能在共享库可用之前,据我所知.

现在编译器可以优化标准C库函数sin()(由提供者libm.so)替换为对CPU/FPU的内置sin()函数的本机指令的调用,该函数作为FPU指令(FSIN对于x86/x87)存在更新的处理器,如Core 2系列(这与i486DX相当正确).这将取决于传递给gcc编译器的优化标志.如果编译器被告知编写可在任何i386或更新的处理器上执行的代码,则不会进行这样的优化.该-mcpu=486标志将通知编译器进行这样的优化是安全的.

现在,如果程序执行sin()函数的软件版本,它将基于CORDIC(坐标旋转数字计算机)或BKM算法,或更可能是现在常用于计算的表格或幂级数计算这样的先验功能.[Src:http://en.wikipedia.org/wiki/Cordic#Application]

任何最近的(自2.9x左右)版本的gcc也提供了sin的内置版本,__builtin_sin()它将用于替换对C库版本的标准调用,作为优化.

我确信它和泥浆一样清晰,但希望能为您提供比您预期更多的信息,还有许多跳跃点可以让您自己了解更多信息.

Answer 13

如果你想在软件而不是硬件中实现,那么寻找这个问题的明确答案的地方就是数字食谱的第5章.我的副本放在一个盒子里,所以我不能提供详细信息,但简短的版本(如果我记得这个权利)是你tan(theta/2)作为原始操作并从那里计算其他人.计算是通过一系列近似来完成的,但它的收敛速度比泰勒级数要快得多.

对不起,如果不把手放在书上,我就无法记住更多.

Answer 14

正如许多人所指出的那样,它依赖于实现.但据我了解你的问题,你对数学函数的真实软件实现很感兴趣,但只是没有设法找到它.如果是这种情况,那么你在这里:

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

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

我希望这有帮助.

Answer 15

没有什么比击中源代码并看到某个人在一个常用的库中实际完成它的方式; 让我们看一下特别是一个C库实现.我选择了uLibC.

这是sin函数:

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

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

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

然后对这两部分进行操作.如果没有尾部,则使用13次多项式生成近似答案.如果有尾部,则基于以下原则得到小的校正加法:sin(x+y) = sin(x) + sin'(x')y