C如何计算sin()和其他数学函数?
Posted
tags:
篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了C如何计算sin()和其他数学函数?相关的知识,希望对你有一定的参考价值。
我一直在研究.NET反汇编和GCC源代码,但似乎无法找到sin()
和其他数学函数的实际实现...它们似乎总是引用其他东西。
谁能帮我找到它们?我觉得C运行的所有硬件都不太可能支持硬件中的触发功能,因此某处必须有软件算法,对吧?
我知道有几种方法可以计算函数,并编写了我自己的例程来计算函数使用泰勒系列来获得乐趣。我很好奇真正的生产语言是如何做到的,因为我的所有实现总是慢几个数量级,即使我认为我的算法非常聪明(显然它们不是)。
在GNU libm中,sin
的实现依赖于系统。因此,您可以在sysdeps的相应子目录中找到每个平台的实现。
一个目录包含由IBM提供的C实现。自2011年10月以来,这是在典型的x86-64 Linux系统上调用sin()
时实际运行的代码。它显然比fsin
汇编指令快。源代码:sysdeps/ieee754/dbl-64/s_sin.c,寻找__sin (double x)
。
这段代码非常复杂。没有一种软件算法尽可能快,并且在整个x值范围内也是准确的,因此库实现了许多不同的算法,它的第一项工作是查看x并决定使用哪种算法。在某些地区,它使用的是看似熟悉的泰勒系列。一些算法首先计算快速结果,然后如果不够准确,则丢弃它并回退到较慢的算法。
较旧的32位版本的GCC / glibc使用了fsin
指令,这对于某些输入而言令人惊讶地不准确。有一个fascinating blog post illustrating this with just 2 lines of code。
fdlibm在纯C中实现sin
比glibc更简单,并且评论很好。源代码:fdlibm/s_sin.c和fdlibm/k_sin.c
关于像sin()
,cos()
,tan()
这样的三角函数,在5年后,没有提到高质量trig函数的一个重要方面:范围缩小。
任何这些函数的早期步骤是以弧度将角度减小到2 *π间隔的范围。但是π是不合理的,因此像x = remainder(x, 2*M_PI)
那样简单的减少会引入误差,因为M_PI
或机器pi是π的近似值。那么,如何做x = remainder(x, 2*π)
?
早期的图书馆使用扩展精度或精心编程来提供高质量的结果,但仍然在有限的double
范围内。当像sin(pow(2,30))
那样请求一个大的值时,结果是无意义的或0.0
,并且可能与error flag设置为像TLOSS
完全丧失精度或PLOSS
部分精度损失。
将大值的大范围减小到像-π到π这样的区间是一个具有挑战性的问题,可以与sin()
本身的基本触发函数相媲美。
一个好的报告是Argument reduction for huge arguments: Good to the last bit (1992)。它很好地解决了这个问题:讨论了各种平台(SPARC,PC,HP,30+其他)的需求和方式,并提供了一种解决方案算法,可以为double
和-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()
它们通常以软件实现,并且在大多数情况下不会使用相应的硬件(即,aseembly)调用。但是,正如Jason指出的那样,这些都是特定于实现的。
请注意,这些软件例程不是编译器源代码的一部分,而是可以在对应的库中找到,例如clib,或GNU编译器的glibc。见http://www.gnu.org/software/libc/manual/html_mono/libc.html#Trig-Functions
如果您想要更好的控制,您应该仔细评估您的需求。一些典型的方法是查找表的插值,汇编调用(通常很慢),或其他近似方案,如平方根的Newton-Raphson。
如果你想在软件而不是硬件中实现,那么寻找这个问题的明确答案的地方是Numerical Recipes的第5章。我的副本放在一个盒子里,所以我不能提供详细信息,但是短版本(如果我还记得这个)是你将tan(theta/2)
作为原始操作并从那里计算其他的。计算是通过一系列近似来完成的,但它的收敛速度比泰勒级数要快得多。
对不起,如果不把手放在书上,我记不住了。
正如许多人所指出的那样,它依赖于实现。但据我了解你的问题,你对数学函数的真实软件实现很感兴趣,但只是没有设法找到它。如果是这种情况,那么你在这里:
- 从http://ftp.gnu.org/gnu/glibc/下载glibc源代码
- 查看位于解压缩的glibc root \ sysdeps \ ieee754 \ dbl-64文件夹中的文件
dosincos.c
- 类似地,您可以找到数学库其余部分的实现,只需查找具有适当名称的文件即可
您还可以查看具有.tbl
扩展名的文件,它们的内容只不过是二进制形式的不同函数的预计算值的巨大表。这就是为什么实现如此之快:不是计算他们使用的任何系列的所有系数,而是快速查找,这要快得多。顺便说一句,他们确实使用Tailor系列来计算正弦和余弦。
我希望这有帮助。
我将尝试在C程序中回答sin()
的情况,在当前的x86处理器上使用GCC的C编译器进行编译(比如英特尔酷睿2双核处理器)。
在C语言中,标准C库包括常见的数学函数,不包括在语言本身中(例如分别用于幂,正弦和余弦的pow
,sin
和cos
)。其标题包含在math.h中。
现在在GNU / Linux系统上,这些库函数由glibc(GNU libc或GNU C Library)提供。但GCC编译器希望您使用math library编译器标志链接到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 algorithm,或者更可能是现在常用于计算这样的表格或幂级数计算超越功能。 [Src:http://en.wikipedia.org/wiki/Cordic#Application]
任何最近的(自2.9x左右)版本的gcc也提供了一个内置版本的sin,__builtin_sin()
,它将用来取代对C库版本的标准调用,作为优化。
我确信它和泥浆一样清晰,但希望能为您提供比您预期更多的信息,还有许多跳跃点可以让您自己了解更多信息。
没有什么比击中源代码并看到某个人在一个常用的库中实际完成它的方式;让我们看一下特别是一个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
的原则得到一个小的校正加法
无论什么时候评估这样的函数,那么在某种程度上,最有可能:
- 插值的值表(用于快速,不准确的应用程序 - 例如计算机图形)
- 对一系列收敛到期望值的评估 - 可能不是泰勒系列,更可能是基于像Clenshaw-Curtis这样的花式正交的东西。
如果没有硬件支持,则编译器可能使用后一种方法,仅发出汇编代码(没有调试符号),而不是使用c库---这使得在调试器中跟踪实际代码变得棘手。
如果你想在C中查看这些函数的实际GNU实现,请查看最新的glibc主干。见GNU C Library。
不要使用泰勒系列。正如上面几个人所指出的那样,切比雪夫多项式更快更准确。这是一个实现(最初来自ZX Spectrum ROM):https://albertveli.wordpress.com/2015/01/10/zx-sine/
通过使用泰勒级数的代码,计算正弦/余弦/正切实际上非常容易。自己写一个需要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&l以上是关于C如何计算sin()和其他数学函数?的主要内容,如果未能解决你的问题,请参考以下文章