C如何计算sin（）和其他数学函数？

问题描述 投票：222回答：20

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

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

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

c math sine
20个回答
192
投票

在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

6
投票

关于像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()

5
投票

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

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

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

5
投票

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

对不起，如果不把手放在书上，我记不住了。

5
投票

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

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

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

我希望这有帮助。

5
投票

我将尝试在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库版本的标准调用，作为优化。

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

5
投票

没有什么比击中源代码并看到某个人在一个常用的库中实际完成它的方式;让我们看一下特别是一个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的原则得到一个小的校正加法

4
投票

无论什么时候评估这样的函数，那么在某种程度上，最有可能：

  • 插值的值表（用于快速，不准确的应用程序 - 例如计算机图形）
  • 对一系列收敛到期望值的评估 - 可能不是泰勒系列，更可能是基于像Clenshaw-Curtis这样的花式正交的东西。

如果没有硬件支持，则编译器可能使用后一种方法，仅发出汇编代码（没有调试符号），而不是使用c库---这使得在调试器中跟踪实际代码变得棘手。

3
投票

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

2
投票

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

1
投票

通过使用泰勒级数的代码，计算正弦/余弦/正切实际上非常容易。自己写一个需要5秒钟。

整个过程可以用这个等式来总结：

sin and cost expansion

以下是我为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));  
}
70
投票

好孩子，专业人士的时间......这是我对没有经验的软件工程师的最大抱怨之一。他们从头开始计算超越函数（使用泰勒系列），好像在他们的生活中没有人曾经做过这些计算。不对。这是一个定义明确的问题，非常聪明的软件和硬件工程师已经接触过数千次，并且有一个定义明确的解决方案。基本上，大多数超越函数使用切比雪夫多项式来计算它们。至于使用哪种多项式取决于具体情况。首先，关于这件事的圣经是哈特和切尼的一本名为“计算机近似”的书。在该书中，您可以决定是否有硬件加法器，乘法器，除法器等，并确定哪些操作最快。例如如果你有一个非常快的分频器，计算正弦的最快方法可能是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中的超越操作的大多数硬件实现通常涉及一些微代码和这样的操作（取决于硬件）。切比雪夫多项式用于大多数超越，但不是全部。例如使用查找表首先使用牛顿raphson方法的双重迭代，平方根更快。同样，那本书“计算机近似”会告诉你。

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

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

请享用。

-1
投票

如果你想要sin那么

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

如果你想要cos那么

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

如果你想要sqrt那么

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

那么为什么在机器指令时会使用不准确的代码呢？

63
投票

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

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

14
投票

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

}
12
投票

是的，有用于计算sin的软件算法。基本上，使用数字计算机计算这类东西通常使用numerical methods来完成，例如近似代表函数的Taylor series。

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

11
投票

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

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

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

11
投票

使用Taylor series并尝试查找系列术语之间的关系，这样您就不会一次又一次地计算事物

以下是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;
}

使用这个，我们可以使用已经使用的一个来获得总和的新项（我们避免使用阶乘和x2p）

explanation

10
投票

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

在某些情况下，最大误差不是您感兴趣的，而是最大相对误差。例如，对于正弦函数，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 ^ 6/5040 ...），那么括号中的函数正在减少，并且如果y是x的下一个更大的数字，那么有时sin y将小于sin 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之间的差异将非常大。现在用多项式（b x ^ 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倍。

6
投票

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

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

最新问题
© www.soinside.com 2019 - 2024. All rights reserved.