C는 sin()과 다른 함수들을 어떻게 계산합니까?

C는 sin()과 다른 함수들을 어떻게 계산합니까?

자세히 살펴봤어요.NET 디스어셈블리와 GCC 소스 코드를 사용하지만 실제 구현은 어디서도 찾을 수 없는 것 같습니다.sin()그리고 다른 수학 함수들...항상 다른 걸 언급하는 것 같아요

찾는 거 도와줄 사람 있어요?C가 동작하는 모든 하드웨어가 하드웨어에서 trig 함수를 지원할 가능성은 낮다고 생각하기 때문에 소프트웨어 알고리즘이 어딘가에 있을 것입니다.

---

함수를 계산하는 방법은 몇 가지 알고 있으며, Taylor 시리즈를 사용하여 함수를 계산하는 루틴을 재미로 작성했습니다.실제 운영 언어가 얼마나 실제적인지 궁금해요. 제 모든 구현은 항상 몇 배 느리지만, 제 알고리즘은 상당히 영리하다고 생각됩니다(분명히 그렇지 않습니다).

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가 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. 분할이나 분기 없이 하는 방법은 꽤 영리합니다.근데 댓글이 하나도 없어요.

---

GCC/glibc의 이전 32비트 버전에서는fsin일부 입력에서는 놀라울 정도로 부정확한 명령입니다.2줄의 코드로 이것을 설명하는 매력적인 블로그 투고가 있습니다.

fdlibm의 실장sin순수한 C는 glibc보다 훨씬 단순하고 평이 좋다.소스 코드: fdlibm/s_sin.c 및 fdlibm/k_sin.c

사인이나 코사인 등의 기능은 마이크로프로세서 내부에서 마이크로코드로 구현됩니다.예를 들어 인텔 칩에는 조립 명령이 있습니다.C 컴파일러는 이러한 어셈블리 명령을 호출하는 코드를 생성합니다.(반대적으로 Java 컴파일러는 그렇지 않습니다.Java는 trig 함수를 하드웨어가 아닌 소프트웨어로 평가하므로 실행 속도가 훨씬 느립니다.)

칩은 Taylor 급수를 사용하여 트리거 함수를 계산하지 않습니다(적어도 전체가 아님).먼저 CORDIC을 사용하지만 CORDIC의 결과를 다듬기 위해 짧은 Taylor 시리즈를 사용하거나 매우 작은 각도에 대해 높은 상대적 정확도로 사인(sine)을 계산하는 등의 특수한 경우 사용할 수 있습니다.자세한 내용은 다음 StackOverflow 답변을 참조하십시오.

다음과 같은 삼각함수에 대하여sin(),cos(),tan()5년이 지난 지금까지 고품질 트리거 함수의 중요한 측면인 범위 축소에 대한 언급은 없었다.

이러한 기능의 초기 단계는 각도를 라디안 단위로 2*µ 간격의 범위로 줄이는 것입니다.하지만 is는 비합리적이기 때문에 다음과 같이 간단히 줄일 수 있습니다.x = remainder(x, 2*M_PI)과오를 초래하다M_PImachine 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

다양한 트리거 ID와remquo()한층 더 개선하다샘플: sind()

위해서sin특히 Taylor 확장을 사용하면 다음과 같은 이점을 얻을 수 있습니다.

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;

}

Blindy의 답변에서 개선된 코드 버전

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

체비셰프 다항식은 다른 답변에서 언급했듯이 함수와 다항식 사이의 가장 큰 차이가 가능한 한 작은 다항식이다.그것은 훌륭한 시작이다.

경우에 따라서는 최대 오차는 관심사가 아닌 최대 상대 오차입니다.예를 들어 사인 함수의 경우 x = 0에 가까운 오차는 큰 값보다 훨씬 작아야 합니다. 따라서 상대 오차는 작아야 합니다.따라서 sin x / x에 대한 체비셰프 다항식을 계산하고 해당 다항식에 x를 곱합니다.

다음으로 다항식을 평가하는 방법을 알아내야 합니다.중간 값이 작으므로 반올림 오차가 작도록 평가하려고 합니다.그렇지 않으면 반올림 오차가 다항식의 오차보다 훨씬 커질 수 있습니다.사인함수와 같은 함수의 경우, 부주의하면 x < y일 때에도 sin x에 대해 계산한 결과가 sin y에 대한 결과보다 클 수 있습니다.따라서 계산 순서를 신중하게 선택하고 반올림 오차에 대한 상한을 계산해야 합니다.

For example, sin x = x - x^3/6 + x^5 / 120 - x^7 / 5040... If you calculate naively sin x = x * (1 - x^2/6 + x^4/120 - x^6/5040...), then that function in parentheses is decreasing, and it will happen that if y is the next larger number to x, then sometimes sin y will be smaller than sin x. Instead, calculate sin x = x - x^3 * (1/6 - x^2 / 120 + x^4/5040...) where this cannot happen.

When calculating Chebyshev polynomials, you usually need to round the coefficients to double precision, for example. But while a Chebyshev polynomial is optimal, the Chebyshev polynomial with coefficients rounded to double precision is not the optimal polynomial with double precision coefficients!

For example for sin (x), where you need coefficients for x, x^3, x^5, x^7 etc. you do the following: Calculate the best approximation of sin x with a polynomial (ax + bx^3 + cx^5 + dx^7) with higher than double precision, then round a to double precision, giving A. The difference between a and A would be quite large. Now calculate the best approximation of (sin x - Ax) with a polynomial (b x^3 + cx^5 + dx^7). You get different coefficients, because they adapt to the difference between a and A. Round b to double precision B. Then approximate (sin x - Ax - Bx^3) with a polynomial cx^5 + dx^7 and so on. You will get a polynomial that is almost as good as the original Chebyshev polynomial, but much better than Chebyshev rounded to double precision.

Next you should take into account the rounding errors in the choice of polynomial. You found a polynomial with minimum error in the polynomial ignoring rounding error, but you want to optimise polynomial plus rounding error. Once you have the Chebyshev polynomial, you can calculate bounds for the rounding error. Say f (x) is your function, P (x) is the polynomial, and E (x) is the rounding error. You don't want to optimise | f (x) - P (x) |, you want to optimise | f (x) - P (x) +/- E (x) |. You will get a slightly different polynomial that tries to keep the polynomial errors down where the rounding error is large, and relaxes the polynomial errors a bit where the rounding error is small.

이렇게 하면 마지막 비트의 최대 0.55배의 오차를 쉽게 반올림할 수 있습니다.여기서 +,-*/는 마지막 비트의 최대 0.50배의 반올림 오차를 가집니다.

이러한 함수의 실제 GNU 실장을 C에서 확인하려면 glibc의 최신 트렁크를 확인하십시오.GNU C 라이브러리를 참조해 주세요.

Taylor 시리즈는 사용하지 마세요.위의 두 사람이 지적한 바와 같이 체비셰프 다항식은 더 빠르고 정확하다.실장(주로 ZX Spectrum ROM에서)을 다음에 나타냅니다.https://albertveli.wordpress.com/2015/01/10/zx-sine/

네가 원한다면sin그리고나서

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

네가 원한다면cos그리고나서

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

네가 원한다면sqrt그리고나서

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

기계 지시가 있는데 왜 부정확한 코드를 사용합니까?

출처를 찾아보고, 실제로 라이브러리에서 어떤 사람이 공통적으로 사용하고 있는지 보는 것만큼 좋은 것은 없습니다.특히 C 라이브러리의 구현에 대해 살펴보겠습니다.저는 uLibC를 선택했습니다.

죄 함수는 다음과 같습니다.

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

사인/코사인/탱겐트 계산은 실제로 Taylor 시리즈를 사용하여 코드를 통해 매우 쉽게 수행할 수 있습니다.직접 쓰는데 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));  
}

그것은 복잡한 질문이다.x86 패밀리의 인텔과 같은 CPU에는,sin()기능하지만 x87 FPU의 일부이며 64비트 모드에서는 더 이상 사용되지 않습니다(대신 SSE2 레지스터가 사용됩니다.이 모드에서는 소프트웨어 실장이 사용됩니다.

그러한 구현은 몇 가지 밖에 있다.하나는 fdlibm으로 Java에서 사용됩니다.제가 알기로는 glibc 구현에는 fdlibm의 일부와 IBM이 제공한 다른 부분이 포함되어 있습니다.

Software implementations of transcendental functions such as sin() typically use approximations by polynomials, often obtained from Taylor series.

Yes, there are software algorithms for calculating sin too. Basically, calculating these kind of stuff with a digital computer is usually done using numerical methods like approximating the Taylor series representing the function.

Numerical methods can approximate functions to an arbitrary amount of accuracy and since the amount of accuracy you have in a floating number is finite, they suit these tasks pretty well.

Use Taylor series and try to find relation between terms of the series so you don't calculate things again and again

Here is an example for 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;
}

using this we can get the new term of the sum using the already used one (we avoid the factorial and x^2p)

The actual implementation of library functions is up to the specific compiler and/or library provider. Whether it's done in hardware or software, whether it's a Taylor expansion or not, etc., will vary.

I realize that's absolutely no help.

They are typically implemented in software and will not use the corresponding hardware (that is, aseembly) calls in most cases. However, as Jason pointed out, these are implementation specific.

Note that these software routines are not part of the compiler sources, but will rather be found in the correspoding library such as the clib, or glibc for the GNU compiler. See http://www.gnu.org/software/libc/manual/html_mono/libc.html#Trig-Functions

If you want greater control, you should carefully evaluate what you need exactly. Some of the typical methods are interpolation of look-up tables, the assembly call (which is often slow), or other approximation schemes such as Newton-Raphson for square roots.

If you want an implementation in software, not hardware, the place to look for a definitive answer to this question is Chapter 5 of Numerical Recipes. My copy is in a box, so I can't give details, but the short version (if I remember this right) is that you take tan(theta/2) as your primitive operation and compute the others from there. The computation is done with a series approximation, but it's something that converges much more quickly than a Taylor series.

Sorry I can't rembember more without getting my hand on the book.

Whenever such a function is evaluated, then at some level there is most likely either:

• A table of values which is interpolated (for fast, inaccurate applications - e.g. computer graphics)
• The evaluation of a series that converges to the desired value --- probably not a taylor series, more likely something based on a fancy quadrature like Clenshaw-Curtis.

If there is no hardware support then the compiler probably uses the latter method, emitting only assembler code (with no debug symbols), rather than using a c library --- making it tricky for you to track the actual code down in your debugger.

As many people pointed out, it is implementation dependent. But as far as I understand your question, you were interested in a real software implemetnation of math functions, but just didn't manage to find one. If this is the case then here you are:

• Download glibc source code from http://ftp.gnu.org/gnu/glibc/
• Look at file dosincos.c located in unpacked glibc root\sysdeps\ieee754\dbl-64 folder
• Similarly you can find implementations of the rest of the math library, just look for the file with appropriate name

You may also have a look at the files with the .tbl extension, their contents is nothing more than huge tables of precomputed values of different functions in a binary form. That is why the implementation is so fast: instead of computing all the coefficients of whatever series they use they just do a quick lookup, which is much faster. BTW, they do use Tailor series to calculate sine and cosine.

I hope this helps.

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.

Note that most hardware implementations of transcendental operations in an FPU usually involve some microcode and operations like this (depends on the hardware). Chebyshev polynomials are used for most transcendentals but not all. e.g. Square root is faster to use a double iteration of Newton raphson method using a lookup table first. Again, that book "Computer Approximations" will tell you that.

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.

Enjoy.

I'll try to answer for the case of sin() in a C program, compiled with GCC's C compiler on a current x86 processor (let's say a Intel Core 2 Duo).

In the C language the Standard C Library includes common math functions, not included in the language itself (e.g. pow, sin and cos for power, sine, and cosine respectively). The headers of which are included in math.h.

Now on a GNU/Linux system, these libraries functions are provided by glibc (GNU libc or GNU C Library). But the GCC compiler wants you to link to the math library (libm.so) using the -lm compiler flag to enable usage of these math functions. I'm not sure why it isn't part of the standard C library. These would be a software version of the floating point functions, or "soft-float".

Aside: The reason for having the math functions separate is historic, and was merely intended to reduce the size of executable programs in very old Unix systems, possibly before shared libraries were available, as far as I know.

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.

Now if the program executed the software version of the sin() function, it would do so based on a CORDIC (COordinate Rotation DIgital Computer) or BKM algorithm, or more likely a table or power-series calculation which is commonly used now to calculate such transcendental functions. [Src: http://en.wikipedia.org/wiki/Cordic#Application]

Any recent (since 2.9x approx.) version of gcc also offers a built-in version of sin, __builtin_sin() that it will used to replace the standard call to the C library version, as an optimization.

I'm sure that is as clear as mud, but hopefully gives you more information than you were expecting, and lots of jumping off points to learn more yourself.

The essence of how it does this lies in this excerpt from Applied Numerical Analysis by Gerald Wheatley:

When your software program asks the computer to get a value of or , have you wondered how it can get the values if the most powerful functions it can compute are polynomials? It doesnt look these up in tables and interpolate! Rather, the computer approximates every function other than polynomials from some polynomial that is tailored to give the values very accurately.

A few points to mention on the above is that some algorithms do infact interpolate from a table, albeit only for the first few iterations. Also note how it mentions that computers utilise approximating polynomials without specifying which type of approximating polynomial. As others in the thread have pointed out, Chebyshev polynomials are more efficient than Taylor polynomials in this case.

ReferenceURL : https://stackoverflow.com/questions/2284860/how-does-c-compute-sin-and-other-math-functions