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

是的，也有计算sin的软件算法。基本上，用数字计算机计算这些东西通常是用数值方法来完成的，就像近似表示函数的泰勒级数一样。

数值方法可以将函数近似到任意精度，因为浮点数的精度是有限的，所以它们非常适合这些任务。

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

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

如果你想要更大的控制权，你应该仔细评估你到底需要什么。一些典型的方法是查找表的插值、程序集调用(通常很慢)或其他近似方案，如Newton-Raphson的平方根。

特别是对于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;


}

像正弦和余弦这样的函数是在微处理器内部的微码中实现的。例如，英特尔芯片就有相应的组装指令。C编译器将生成调用这些汇编指令的代码。(相反，Java编译器不会。Java在软件而不是硬件中计算三角函数，因此运行速度要慢得多。)

芯片不使用泰勒级数来计算三角函数，至少不完全是。首先，他们使用CORDIC，但他们也可能使用一个短的泰勒级数来优化CORDIC的结果，或者用于特殊情况，例如在非常小的角度下以相对较高的精度计算正弦。有关更多解释，请参见StackOverflow回答。

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

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

诸如sin()之类的超越函数的软件实现通常使用多项式逼近，通常从泰勒级数获得。

使用泰勒级数并尝试找到级数项之间的关系，这样你就不会一次又一次地计算东西

下面是一个关于余窦的例子:

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^{2 p})来获得新的和项。

无论何时这样一个函数被求值，那么在某种程度上很可能有:

• 内插的值表(用于快速，不准确的应用程序-例如计算机图形)
• 收敛到期望值的级数的求值——可能是不一个泰勒级数，更可能是基于像Clenshaw-Curtis这样的花哨正交的东西。

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

在GNU libm中，sin的实现是系统相关的。因此，你可以在sysdeps的适当子目录中找到每个平台的实现。

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

这段代码非常复杂。没有一种软件算法在整个x值范围内尽可能快且准确，因此标准库实现了几种不同的算法，它的第一项工作是查看x并决定使用哪种算法。

• 当x非常非常接近0时，sin(x) == x是正确答案。

• 再往前一点，sin(x)使用我们熟悉的Taylor级数。然而，这只在接近0时是准确的，所以……

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

• 当|x| > 2时，上述算法都不能工作，因此代码开始计算一些接近0的值，可以将其提供给sin或cos。

• 还有另一个分支要处理x是NaN或无穷大。

这段代码使用了一些我以前从未见过的数值技巧，尽管据我所知，它们可能在浮点专家中很有名。有时几行代码需要几段文字来解释。例如，这两条线

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

(有时)用于将x降为接近0的值，该值与x相差π/2的倍数，特别是xn ×π/ 2。这种没有划分或分支的方式相当聪明。但是没有任何评论!

旧的32位版本的GCC/glibc使用fsin指令，这对于某些输入是非常不准确的。有一个精彩的博客文章用两行代码说明了这一点。

fdlibm在纯C中实现sin要比glibc简单得多，而且注释很好。源代码:fdlibm / s_sin.c和fdlibm / k_sin.c

我将尝试在一个C程序中回答sin()的情况，用GCC的C编译器在当前的x86处理器(让我们说一个英特尔酷睿2 Duo)上编译。

在C语言中，标准C库包含一些常见的数学函数，而这些函数并不包含在语言本身中(例如，pow、sin和cos分别表示幂、正弦和余弦)。它们的头文件包含在math.h中。

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

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

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

现在，如果程序执行sin()函数的软件版本，它将基于CORDIC(坐标旋转数字计算机)或BKM算法或更多的来执行，可能是现在通常用于计算此类超越函数的表格或幂级数计算。[Src: http://en.wikipedia.org/wiki/Cordic应用程序)

gcc的任何最新版本(大约2.9x以后)也提供了内置的sin __builtin_sin()版本，作为优化，它将用于取代对C库版本的标准调用。

我相信这是非常清楚的，但希望给你更多的信息比你期望的，和许多出发点，以了解更多自己。

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

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

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

如果你想要在软件而不是硬件中实现，可以在数值的食谱的第5章中找到这个问题的明确答案。我的副本在一个盒子里，所以我不能给出详细信息，但简短的版本(如果我没记错的话)是你把tan(theta/2)作为你的基本操作，并从那里计算其他操作。计算是用级数近似完成的，但它收敛多比泰勒级数更快。

抱歉，我没拿到书就想不起来了。

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

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

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

我希望这能有所帮助。

好了，孩子们，是时候上场了.... 这是我对缺乏经验的软件工程师最大的抱怨之一。他们从零开始计算超越函数(使用泰勒级数)，就好像以前从来没有人做过这些计算一样。不正确的。这是一个定义明确的问题，已经被非常聪明的软件和硬件工程师处理了数千次，并且有一个定义明确的解决方案。 基本上，大多数超越函数都使用切比雪夫多项式来计算。至于使用哪些多项式取决于具体情况。首先，关于这个问题的圣经是哈特和切尼合著的《计算机逼近》一书。在那本书中，你可以决定是否有硬件加法器、乘法器、除法器等，并决定哪些操作是最快的。例:如果你有一个非常快的除法器，计算正弦的最快方法可能是P1(x)/P2(x)，其中P1, P2是切比雪夫多项式。如果没有快速除法，它可能只是P(x)，其中P比P1或P2有更多的项....所以会慢一些。因此，第一步是确定硬件及其功能。然后选择切比雪夫多项式的适当组合(例如，对于余弦，通常是cos(ax) = aP(x)的形式，同样，P是切比雪夫多项式)。然后决定需要的十进制精度。例如，如果你想要7位数的精度，你可以在我提到的书中适当的表格中查找，它会给你一个数字N = 4和一个多项式数3502(对于精度= 7.33)。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位的方法

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

享受。

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

如果你想要sin那么

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

如果你想要cos那么

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

如果你想要sqrt那么

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

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

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

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

接下来你要弄清楚如何求多项式的值。你想用这样一种方式来计算它，中间值很小，因此舍入误差也很小。否则舍入误差可能会比多项式中的误差大得多。对于sin函数这样的函数，如果你不小心，你计算sinx的结果可能比siny的结果大，即使x <y.因此需要仔细选择计算顺序，计算舍入误差的上界。

例如，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的下一个更大的数字，那么有时sin y会小于sin x。相反，计算sinx = x - x^3 * (1/6 - x^2/ 120 + x^4/5040…)，这是不可能发生的。

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

以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)计算(sin x - Ax)的最佳近似。你会得到不同的系数，因为它们适应a和a之间的差异，四舍五入b到双精度b，然后用多项式cx^5 + dx^7近似(sin x - Ax - Bx^3)，以此类推。你会得到一个多项式几乎和原来的切比雪夫多项式一样好，但比切比雪夫四舍五入到两倍精度要好得多。

接下来，你应该考虑到舍入误差在多项式的选择。你在忽略舍入误差的多项式中找到了误差最小的多项式，但你想优化多项式加上舍入误差。一旦你有了切比雪夫多项式，你就可以计算舍入误差的边界。假设f (x)是你的函数，P (x)是多项式，E (x)是舍入误差。你不想优化|f (x) - P (x) |，你想优化|f (x) - P (x) +/- E (x) |。你会得到一个稍微不同的多项式，它试图在舍入误差大的地方减小多项式误差，在舍入误差小的地方减小多项式误差。

所有这些将使您轻松地获得最多0.55倍于最后一位的舍入误差，其中+，-，*，/的舍入误差最多为0.50倍于最后一位。

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

没有什么比点击源代码，看看人们是如何在常用的库中实际完成它的了;让我们特别看看一个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

关于像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+其他)上的需求和情况，并提供了一个解决方案算法，从-DBL_MAX到DBL_MAX给出了所有 double的高质量结果。

如果原来的参数是度数，但可能值很大，首先使用fmod()来提高精度。一个好的fmod()将引入没有错误，从而提供出色的范围缩减。

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

各种三角恒等式和remquo()提供了更多的改进。示例:信德()

盲汉回答的改进版代码

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

它如何做到这一点的本质在于Gerald Wheatley从应用数值分析中摘录的这段话:

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

.

.

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