我一直在研究。net反汇编和GCC源代码,但似乎在任何地方都找不到sin()和其他数学函数的实际实现…他们似乎总是在引用别的东西。

有人能帮我找到他们吗?我觉得不太可能运行C语言的所有硬件都支持三角函数,所以一定有软件算法,对吧?


我知道有几种方法可以计算函数,并且为了好玩,我写了自己的例程来使用泰勒级数计算函数。我很好奇生产语言是如何做到这一点的,因为我的所有实现总是慢几个数量级,尽管我认为我的算法非常聪明(显然不是)。


当前回答

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

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

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

题外话:将数学函数分开的原因由来已久,据我所知,可能是在共享库可用之前,它仅仅是为了在非常古老的Unix系统中减少可执行程序的大小。

Now the compiler may optimize the standard C library function sin() (provided by libm.so) to be replaced with an call to a native instruction to your CPU/FPU's built-in sin() function, which exists as an FPU instruction (FSIN for x86/x87) on newer processors like the Core 2 series (this is correct pretty much as far back as the i486DX). This would depend on optimization flags passed to the gcc compiler. If the compiler was told to write code that would execute on any i386 or newer processor, it would not make such an optimization. The -mcpu=486 flag would inform the compiler that it was safe to make such an optimization.

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

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

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

其他回答

OK kiddies, time for the pros.... This is one of my biggest complaints with inexperienced software engineers. They come in calculating transcendental functions from scratch (using Taylor's series) as if nobody had ever done these calculations before in their lives. Not true. This is a well defined problem and has been approached thousands of times by very clever software and hardware engineers and has a well defined solution. Basically, most of the transcendental functions use Chebyshev Polynomials to calculate them. As to which polynomials are used depends on the circumstances. First, the bible on this matter is a book called "Computer Approximations" by Hart and Cheney. In that book, you can decide if you have a hardware adder, multiplier, divider, etc, and decide which operations are fastest. e.g. If you had a really fast divider, the fastest way to calculate sine might be P1(x)/P2(x) where P1, P2 are Chebyshev polynomials. Without the fast divider, it might be just P(x), where P has much more terms than P1 or P2....so it'd be slower. So, first step is to determine your hardware and what it can do. Then you choose the appropriate combination of Chebyshev polynomials (is usually of the form cos(ax) = aP(x) for cosine for example, again where P is a Chebyshev polynomial). Then you decide what decimal precision you want. e.g. if you want 7 digits precision, you look that up in the appropriate table in the book I mentioned, and it will give you (for precision = 7.33) a number N = 4 and a polynomial number 3502. N is the order of the polynomial (so it's p4.x^4 + p3.x^3 + p2.x^2 + p1.x + p0), because N=4. Then you look up the actual value of the p4,p3,p2,p1,p0 values in the back of the book under 3502 (they'll be in floating point). Then you implement your algorithm in software in the form: (((p4.x + p3).x + p2).x + p1).x + p0 ....and this is how you'd calculate cosine to 7 decimal places on that hardware.

请注意,在FPU中大多数硬件实现的超越操作通常涉及一些微码和类似的操作(取决于硬件)。 切比雪夫多项式用于大多数先验多项式,但不是全部。例:使用Newton raphson方法的两次迭代,首先使用查询表,使用平方根更快。 同样,《计算机逼近》这本书会告诉你。

If you plan on implmementing these functions, I'd recommend to anyone that they get a copy of that book. It really is the bible for these kinds of algorithms. Note that there are bunches of alternative means for calculating these values like cordics, etc, but these tend to be best for specific algorithms where you only need low precision. To guarantee the precision every time, the chebyshev polynomials are the way to go. Like I said, well defined problem. Has been solved for 50 years now.....and thats how it's done.

Now, that being said, there are techniques whereby the Chebyshev polynomials can be used to get a single precision result with a low degree polynomial (like the example for cosine above). Then, there are other techniques to interpolate between values to increase the accuracy without having to go to a much larger polynomial, such as "Gal's Accurate Tables Method". This latter technique is what the post referring to the ACM literature is referring to. But ultimately, the Chebyshev Polynomials are what are used to get 90% of the way there.

享受。

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

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

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

盲汉回答的改进版代码

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

关于sin(), cos(),tan()这样的三角函数,在5年之后,没有提到高质量三角函数的一个重要方面:极差约简。

任何这些函数的早期步骤都是将角度(以弧度为单位)减小到2*π区间。但是π是无理数,所以像x =余数(x, 2*M_PI)这样的简单简化会引入误差,因为M_PI或机器pi是π的近似值。那么,如何求x =余数(x, 2*π)呢?

早期的库使用扩展精度或精心设计的编程来提供高质量的结果,但仍然在有限的double范围内。当请求一个较大的值,如sin(pow(2,30))时,结果是无意义的或0.0,并且可能将错误标志设置为TLOSS完全损失精度或PLOSS部分损失精度。

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

一个好的报告是大论点的论据缩减:好到最后一位(1992)。它涵盖了这个问题很好:讨论了需要和事情是如何在各种平台(SPARC, PC, HP, 30+其他),并提供了一个解决方案算法,为所有双从-DBL_MAX到DBL_MAX的高质量结果。


如果原始参数以度为单位,但可能值很大,则首先使用fmod()以提高精度。一个好的fmod()将不会引入任何错误,从而提供出色的范围缩小。

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

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