Jak C oblicza sin () i inne funkcje matematyczne?

Przeglądałem dezasemblacje .NET i kod źródłowy GCC, ale wydaje się, że nigdzie nie mogę znaleźć faktycznej implementacji sin()i innych funkcji matematycznych ... zawsze wydają się odnosić do czegoś innego.

Czy ktoś może mi pomóc je znaleźć? Wydaje mi się, że jest mało prawdopodobne, aby WSZYSTKIE urządzenia, na których będzie działał C, wspierały funkcje wyzwalające w sprzęcie, więc musi być gdzieś algorytm oprogramowania , prawda?

Zdaję sobie sprawę z kilku sposobów, które funkcje mogą być obliczone i pisało własne procedury do obliczania funkcji wykorzystujących szereg Taylora dla zabawy. Ciekawi mnie, jak robią to prawdziwe języki produkcyjne, ponieważ wszystkie moje implementacje są zawsze o kilka rzędów wielkości wolniejsze, chociaż uważam, że moje algorytmy są dość sprytne (oczywiście nie są).

W GNU libm implementacja sinzależy od systemu. Dlatego implementację można znaleźć dla każdej platformy, gdzieś w odpowiednim podkatalogu sysdeps .

Jeden katalog zawiera implementację w C, wniesioną przez IBM. Od października 2011 r. Jest to kod, który faktycznie działa, gdy wywołujesz sin()typowy system Linux x86-64. Jest najwyraźniej szybszy niż fsininstrukcja montażu. Kod źródłowy: sysdeps / ieee754 / dbl-64 / s_sin.c , poszukaj __sin (double x).

Ten kod jest bardzo złożony. Żaden algorytm programowy nie jest tak szybki, jak to możliwe, a także dokładny w całym zakresie wartości x , więc biblioteka implementuje kilka różnych algorytmów, a jego pierwszym zadaniem jest przyjrzenie się x i zdecydowanie, którego algorytmu użyć.

• Gdy x jest bardzo bardzo bliskie zeru, sin(x) == xjest właściwą odpowiedzią.

• Nieco dalej sin(x)wykorzystuje znaną serię Taylor. Jest to jednak dokładne tylko w pobliżu 0, więc ...

• Gdy kąt jest większy niż około 7 °, stosuje się inny algorytm, obliczając przybliżenia szeregów Taylora zarówno dla sin (x), jak i cos (x), a następnie używając wartości ze wstępnie obliczonej tabeli w celu uściślenia przybliżenia.

• Kiedy | x | > 2, żaden z powyższych algorytmów nie działałby, więc kod zaczyna się od obliczenia wartości bliższej 0, którą można podać sinlub coszamiast niej.

• Jest jeszcze jedna gałąź, w której x jest NaN lub nieskończonością.

W tym kodzie wykorzystano kilka numerycznych hacków, których nigdy wcześniej nie widziałem, choć dla wszystkich wiem, że mogą być dobrze znane wśród ekspertów zmiennoprzecinkowych. Czasami wyjaśnienie wymaga kilku wierszy kodu. Na przykład te dwie linie

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

są używane (czasami) w zmniejszaniu x do wartości bliskiej 0, która różni się od x wielokrotnością π / 2, a konkretnie xn× π / 2. Sposób, w jaki odbywa się to bez podziału lub rozgałęzienia, jest dość sprytny. Ale w ogóle nie ma komentarza!

Starsze 32-bitowe wersje GCC / glibc korzystały z fsininstrukcji, która jest zaskakująco niedokładna w przypadku niektórych danych wejściowych. Jest fascynujący post na blogu ilustrujący to za pomocą zaledwie 2 linii kodu .

Implementacja fdlibm sinw czystym C jest znacznie prostsza niż glibc i jest ładnie komentowana. Kod źródłowy: fdlibm / s_sin.c i fdlibm / k_sin.c

Funkcje takie jak sinus i cosinus są zaimplementowane w mikrokodzie wewnątrz mikroprocesorów. Na przykład układy Intel mają instrukcje montażu. Kompilator prądu przemiennego wygeneruje kod wywołujący instrukcję montażu. (W przeciwieństwie do tego kompilator Java nie. Java ocenia funkcje triggera w oprogramowaniu, a nie w sprzęcie, dlatego działa znacznie wolniej.)

Chipy nie używają szeregów Taylora do obliczania funkcji triggera, przynajmniej nie do końca. Przede wszystkim używają CORDIC , ale mogą również użyć krótkiej serii Taylora, aby dopracować wynik CORDIC lub w szczególnych przypadkach, takich jak obliczenie sinusoidy z wysoką względną dokładnością dla bardzo małych kątów. Aby uzyskać więcej wyjaśnień, zobacz odpowiedź StackOverflow .

OK dzieciaki, czas na profesjonalistów ... To jedna z moich największych skarg z niedoświadczonymi inżynierami oprogramowania. Przychodzą one do obliczania funkcji transcendentalnych od zera (przy użyciu serii Taylora), tak jakby nikt nigdy wcześniej nie dokonywał tych obliczeń w swoim życiu. Nie prawda. Jest to dobrze zdefiniowany problem, do którego tysiące razy podchodzili bardzo sprytni inżynierowie oprogramowania i sprzętu i ma dobrze zdefiniowane rozwiązanie. Zasadniczo większość funkcji transcendentalnych korzysta z wielomianów Czebyszewa do ich obliczania. To, które wielomiany są stosowane, zależy od okoliczności. Po pierwsze, Biblia na ten temat jest książką Hart and Cheney zatytułowaną „Computer Approximations”. W tej książce możesz zdecydować, czy masz sumator sprzętowy, mnożnik, dzielnik itp. I zdecydować, które operacje są najszybsze. np. jeśli miałeś naprawdę szybki dzielnik, najszybszym sposobem obliczenia sinusa może być P1 (x) / P2 (x), gdzie P1, P2 to wielomiany Czebyszewa. Bez szybkiego dzielnika może to być po prostu P (x), gdzie P ma znacznie więcej terminów niż P1 lub P2 .... więc byłoby wolniej. Pierwszym krokiem jest więc określenie sprzętu i jego możliwości. Następnie wybierasz odpowiednią kombinację wielomianów Czebiszewa (zwykle ma postać cos (ax) = aP (x) na przykład dla cosinusa, ponownie, gdzie P jest wielomianem Czebyszewa). Następnie decydujesz, jakiej precyzji dziesiętnej chcesz. np. jeśli chcesz precyzji 7 cyfr, spójrz na to w odpowiedniej tabeli we wspomnianej książce i da ci (dla precyzji = 7,33) liczbę N = 4 i wielomian 3502. N jest rzędu wielomian (więc jest to p4.x ^ 4 + p3.x ^ 3 + p2.x ^ 2 + p1.x + p0), ponieważ N = 4. Następnie sprawdzasz rzeczywistą wartość p4, p3, p2, p1, wartości p0 z tyłu książki poniżej 3502 (będą w liczbach zmiennoprzecinkowych). Następnie implementujesz swój algorytm w oprogramowaniu w postaci: (((p4.x + p3) .x + p2) .x + p1) .x + p0 .... i tak obliczysz cosinus do 7 miejsc po przecinku miejsca na tym sprzęcie.

Zauważ, że większość sprzętowych implementacji operacji transcendentalnych w FPU zwykle wymaga jakiegoś mikrokodu i podobnych operacji (zależy od sprzętu). Wielomiany Czebyszewa są używane w przypadku większości transcendentałów, ale nie wszystkich. Np. pierwiastek kwadratowy jest szybszy w użyciu podwójnej iteracji metody Newtona raphsona przy użyciu najpierw tabeli odnośników. Znów powie ci to książka „Computer Approximations”.

Jeśli planujesz wdrożyć te funkcje, polecam każdemu, aby otrzymał kopię tej książki. To naprawdę Biblia dla tego rodzaju algorytmów. Zauważ, że istnieją wiązki alternatywnych sposobów obliczania tych wartości, takich jak cordics itp., Ale zwykle są one najlepsze dla określonych algorytmów, w których potrzebujesz tylko niskiej precyzji. Aby zagwarantować precyzję za każdym razem, wielomian Czebyszewa jest właściwą drogą. Tak jak powiedziałem, dobrze zdefiniowany problem. Został rozwiązany od 50 lat ... i tak to się robi.

To powiedziawszy, istnieją techniki, dzięki którym wielomiany Czebyszewa mogą być użyte do uzyskania pojedynczej precyzji wyniku z wielomianem niskiego stopnia (jak przykład dla cosinusa powyżej). Następnie istnieją inne techniki interpolacji między wartościami w celu zwiększenia dokładności bez konieczności przechodzenia do znacznie większego wielomianu, takie jak „metoda dokładnych tabel Gal”. Ta ostatnia technika odnosi się do postu odnoszącego się do literatury ACM. Ale ostatecznie wielomian Czebeszewa służy do uzyskania 90% drogi.

Cieszyć się.

W sinszczególności użycie rozszerzenia Taylor dałoby ci:

sin (x): = x - x ^ 3/3! + x ^ 5/5! - x ^ 7/7! + ... (1)

dodawałeś terminy, dopóki różnica między nimi nie będzie niższa niż akceptowany poziom tolerancji lub tylko na skończoną liczbę kroków (szybciej, ale mniej precyzyjnie). Przykładem może być coś takiego:

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

Uwaga: (1) działa z powodu aproksymacji sin (x) = x dla małych kątów. W przypadku większych kątów musisz obliczyć coraz więcej terminów, aby uzyskać akceptowalne wyniki. Możesz użyć argumentu while i kontynuować dla pewnej dokładności:

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;

}

Tak, istnieją algorytmy programowe do obliczania sin. Zasadniczo obliczanie tego rodzaju rzeczy za pomocą komputera cyfrowego jest zwykle wykonywane przy użyciu metod numerycznych, takich jak aproksymacja szeregu Taylora reprezentującego funkcję.

Metody numeryczne mogą aproksymować funkcje z dowolną dokładnością, a ponieważ dokładność, którą masz w liczbach zmiennoprzecinkowych, jest skończona, całkiem dobrze nadają się do tych zadań.

Użyj serii Taylora i spróbuj znaleźć relację między terminami serii, aby nie obliczać rzeczy raz po raz

Oto przykład dla cosinusa:

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

używając tego możemy uzyskać nowy warunek sumy używając już używanego (unikamy silni i x ^2p )

To złożone pytanie. Procesor podobny do Intela z rodziny x86 ma sprzętową implementację tej sin()funkcji, ale jest częścią FPU x87 i nie jest już używany w trybie 64-bitowym (zamiast tego używane są rejestry SSE2). W tym trybie używana jest implementacja oprogramowania.

Istnieje kilka takich wdrożeń. Jeden jest w fdlibm i jest używany w Javie. O ile mi wiadomo, implementacja glibc zawiera części fdlibm i inne części dostarczone przez IBM.

Implementacje programowe funkcji transcendentalnych, takie jak sin()zwykle wykorzystują aproksymacje wielomianami, często uzyskiwane z serii Taylora.

Wielomiany Czebyszewa, jak wspomniano w innej odpowiedzi, to wielomiany, w których największa różnica między funkcją a wielomianem jest jak najmniejsza. To doskonały początek.

W niektórych przypadkach maksymalny błąd nie jest tym, co Cię interesuje, ale maksymalny błąd względny. Na przykład dla funkcji sinus, błąd w pobliżu x = 0 powinien być znacznie mniejszy niż dla większych wartości; chcesz mały błąd względny . Obliczymy zatem wielomian Czebiszewa dla sin x / x i pomnożymy ten wielomian przez x.

Następnie musisz dowiedzieć się, jak ocenić wielomian. Chcesz to ocenić w taki sposób, że wartości pośrednie są małe, a zatem błędy zaokrąglania są małe. W przeciwnym razie błędy zaokrąglania mogą stać się znacznie większe niż błędy wielomianu. A w przypadku funkcji takich jak funkcja sinusoidalna, jeśli jesteś nieostrożny, może się zdarzyć, że wynik obliczony dla sin x będzie większy niż wynik dla sin y, nawet gdy x <y. Potrzebny jest więc ostrożny wybór kolejności obliczeń i obliczanie górnych granic błędu zaokrąglania.

Na przykład sin x = x - x ^ 3/6 + x ^ 5/120 - x ^ 7/5040 ... Jeśli naiwnie obliczysz sin x = x * (1 - x ^ 2/6 + x ^ 4 / 120 - x ^ 6/5040 ...), wówczas funkcja w nawiasach maleje i stanie się tak, że jeśli y będzie następną większą liczbą od x, to czasami sin y będzie mniejszy niż sin x. Zamiast tego obliczyć sin x = x - x ^ 3 * (1/6 - x ^ 2/120 + x ^ 4/5040 ...) tam, gdzie to nie może się zdarzyć.

Obliczając wielomiany Czebyszewa, zwykle trzeba na przykład zaokrąglić współczynniki do podwójnej precyzji. Ale chociaż wielomian Czebyszewa jest optymalny, wielomian Czebyszewa o współczynnikach zaokrąglonych do podwójnej precyzji nie jest optymalnym wielomianem o współczynnikach podwójnej precyzji!

Na przykład dla sin (x), gdzie potrzebujesz współczynników dla x, x ^ 3, x ^ 5, x ^ 7 itd., Wykonaj następujące czynności: Oblicz najlepsze przybliżenie sin x wielomianem (ax + bx ^ 3 + cx ^ 5 + dx ^ 7) z wyższą niż podwójna precyzja, a następnie zaokrąglamy a do podwójnej precyzji, dając A. Różnica między a i A byłaby dość duża. Teraz obliczyć najlepsze przybliżenie (sin x - Ax) za pomocą wielomianu (bx ^ 3 + cx ^ 5 + dx ^ 7). Otrzymujesz różne współczynniki, ponieważ dostosowują się one do różnicy między a i A. Okrąg b do podwójnej precyzji B. Następnie przybliż (sin x - Ax - Bx ^ 3) za pomocą wielomianu cx ^ 5 + dx ^ 7 i tak dalej. Otrzymasz wielomian, który jest prawie tak dobry, jak oryginalny wielomian Czebeszewa, ale znacznie lepszy niż Czebyshev zaokrąglony do podwójnej precyzji.

Następnie należy wziąć pod uwagę błędy zaokrąglania przy wyborze wielomianu. Znalazłeś wielomian z minimalnym błędem w wielomianu ignorującym błąd zaokrąglania, ale chcesz zoptymalizować wielomian plus błąd zaokrąglania. Po uzyskaniu wielomianu Czebyszewa można obliczyć granice błędu zaokrąglania. Powiedz, że f (x) to twoja funkcja, P (x) to wielomian, a E (x) to błąd zaokrąglania. Nie chcesz optymalizować | f (x) - P (x) |, chcesz zoptymalizować | f (x) - P (x) +/- E (x) | Otrzymasz nieco inny wielomian, który próbuje utrzymać błędy wielomianu na dole tam, gdzie błąd zaokrąglania jest duży, i nieco rozluźnia błędy wielomianu, gdy błąd zaokrąglania jest mały.

Wszystko to pozwoli ci łatwo zaokrąglić błędy co najwyżej 0,55 razy ostatni bit, gdzie +, -, *, / mają błędy zaokrąglenia co najwyżej 0,50 razy ostatni bit.

Jeśli chodzi o funkcje trygonometryczne jak sin(), cos(), tan()nie było wzmianki, po 5 latach, ważnego aspektu wysokiej jakości funkcji trygonometrycznych: redukcją zasięgu .

Wczesnym krokiem każdej z tych funkcji jest zmniejszenie kąta w radianach do zakresu 2 * π. Ale π jest nieracjonalne, więc proste redukcje, takie jak x = remainder(x, 2*M_PI)wprowadzenie błędu as M_PI, lub maszyna pi, są przybliżeniem π. Jak to zrobić x = remainder(x, 2*π)?

Wczesne biblioteki korzystały z rozszerzonej precyzji lub spreparowanego programowania, aby dawać wysokiej jakości wyniki, ale nadal w ograniczonym zakresie double. Gdy zażądano dużej wartości sin(pow(2,30)), wyniki były bez znaczenia lub 0.0być może z flagą błędu ustawioną na coś takiego jak TLOSScałkowita utrata precyzji lub PLOSSczęściowa utrata precyzji.

Dobra redukcja dużych wartości do przedziału, takiego jak -π do π, stanowi trudny problem, który rywalizuje z wyzwaniami podstawowej funkcji trig, jak sin()sama w sobie.

Dobrym raportem jest redukcja argumentów w przypadku ogromnych argumentów: dobry do ostatniego kawałka (1992). Dobrze obejmuje ten problem: omawia potrzebę i sposób działania na różnych platformach (SPARC, PC, HP, 30+ i inne) i zapewnia algorytm rozwiązania, który daje wyniki jakościowe dla wszystkich double od -DBL_MAXdo DBL_MAX.

Jeśli oryginalne argumenty są wyrażone w stopniach, ale mogą mieć dużą wartość, użyj fmod()najpierw dla zwiększenia precyzji. Dobra fmod()nie wprowadzi żadnego błędu, a zatem zapewni doskonałe zmniejszenie zasięgu.

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

Różne tożsamości wyzwalaczy i remquo()oferują jeszcze więcej ulepszeń. Przykład: sind ()

Rzeczywista implementacja funkcji bibliotecznych zależy od konkretnego kompilatora i / lub dostawcy biblioteki. Niezależnie od tego, czy odbywa się to sprzętowo, czy programowo, niezależnie od tego, czy jest to rozszerzenie Taylora, czy nie, itd.

Zdaję sobie sprawę, że to absolutnie nie pomoże.

Zazwyczaj są one zaimplementowane w oprogramowaniu i w większości przypadków nie będą używać odpowiednich wywołań sprzętowych (czyli tak zwanych). Jednak, jak zauważył Jason, są one specyficzne dla implementacji.

Zauważ, że te procedury oprogramowania nie są częścią źródeł kompilatora, ale raczej można je znaleźć w bibliotece odpowiadającej kodowania, takiej jak clib lub glibc dla kompilatora GNU. Zobacz http://www.gnu.org/software/libc/manual/html_mono/libc.html#Trig-Functions

Jeśli chcesz mieć większą kontrolę, powinieneś dokładnie ocenić, czego dokładnie potrzebujesz. Niektóre z typowych metod to interpolacja tabel przeglądowych, wywołanie asemblera (które często jest wolne) lub inne schematy aproksymacji, takie jak Newton-Raphson dla pierwiastków kwadratowych.

Jeśli potrzebujesz implementacji w oprogramowaniu, a nie w sprzęcie, miejscem, w którym można znaleźć ostateczną odpowiedź na to pytanie, jest rozdział 5 przepisów numerycznych . Moja kopia jest w pudełku, więc nie mogę podać szczegółów, ale krótka wersja (jeśli dobrze to pamiętam) polega na tym, że bierzesz tan(theta/2)za swoją prymitywną operację i obliczasz pozostałe stamtąd. Obliczenia są wykonywane z aproksymacją serii, ale jest to coś, co zbiega się znacznie szybciej niż seria Taylora.

Przepraszam, że nie mogę zapamiętać więcej, nie wspierając książki.

Nie ma nic lepszego niż trafienie do źródła i zobaczenie, jak ktoś to zrobił w powszechnie używanej bibliotece; spójrzmy w szczególności na jedną implementację biblioteki C. Wybrałem uLibC.

Oto funkcja grzechu:

http://git.uclibc.org/uClibc/tree/libm/s_sin.c

który wygląda na to, że obsługuje kilka specjalnych przypadków, a następnie dokonuje pewnej redukcji argumentów, aby odwzorować dane wejściowe na zakres [-pi / 4, pi / 4], (dzielenie argumentu na dwie części, dużą część i ogon) przed dzwonieniem

http://git.uclibc.org/uClibc/tree/libm/k_sin.c

który następnie działa na tych dwóch częściach. Jeśli nie ma ogona, przybliżona odpowiedź jest generowana przy użyciu wielomianu stopnia 13. Jeśli istnieje ogon, otrzymujesz niewielki dodatek korygujący oparty na zasadzie, żesin(x+y) = sin(x) + sin'(x')y

Ilekroć taka funkcja jest oceniana, wówczas na pewnym poziomie najprawdopodobniej istnieje:

• Tabela wartości, która jest interpolowana (dla szybkich, niedokładnych aplikacji - np. Grafika komputerowa)
• Ocena szeregu, który zbiega się do pożądanej wartości --- prawdopodobnie nie jest to szereg Taylora, bardziej prawdopodobne coś opartego na fantazyjnej kwadraturze, takiej jak Clenshaw-Curtis.

Jeśli nie ma wsparcia sprzętowego, kompilator prawdopodobnie używa tej drugiej metody, emitując tylko kod asemblera (bez symboli debugowania), zamiast używać biblioteki ac --- co utrudnia śledzenie rzeczywistego kodu w debuggerze.

Jak zauważyło wiele osób, zależy to od wdrożenia. Ale o ile rozumiem twoje pytanie, byłeś zainteresowany rzeczywistą implementacją oprogramowania funkcji matematycznych, ale po prostu nie udało ci się go znaleźć. Jeśli tak jest, to tutaj jesteś:

• Pobierz kod źródłowy glibc ze strony http://ftp.gnu.org/gnu/glibc/
• Spójrz na plik dosincos.cznajdujący się w rozpakowanym folderze głównym \ sysdeps \ ieee754 \ dbl-64 glibc
• Podobnie możesz znaleźć implementacje reszty biblioteki matematycznej, po prostu poszukaj pliku o odpowiedniej nazwie

Możesz także zajrzeć do plików z .tblrozszerzeniem, ich zawartość to nic innego jak ogromne tabele wstępnie obliczonych wartości różnych funkcji w formie binarnej. Właśnie dlatego implementacja jest tak szybka: zamiast obliczać wszystkie współczynniki z jakiejkolwiek serii, której używają, po prostu dokonują szybkiego wyszukiwania, co jest znacznie szybsze. BTW, używają serii Tailor do obliczania sinusa i cosinusa.

Mam nadzieję, że to pomoże.

Spróbuję odpowiedzieć na przypadek sin()programu C skompilowanego z kompilatorem C GCC na bieżącym procesorze x86 (powiedzmy Intel Core 2 Duo).

W języku C Standard C Library zawiera typowych funkcji matematycznych, nie ujęte w samym języku (np pow, sini cosdla władzy, sinus i cosinus, odpowiednio). Nagłówki są zawarte w matematyce . H.

Obecnie w systemie GNU / Linux funkcje bibliotek są dostarczane przez glibc (GNU libc lub GNU C Library). Ale kompilator GCC chce, abyś łączył się z biblioteką matematyczną ( libm.so) przy użyciu -lmflagi kompilatora, aby umożliwić korzystanie z tych funkcji matematycznych. Nie jestem pewien, dlaczego nie jest częścią standardowej biblioteki C. Byłaby to wersja oprogramowania funkcji zmiennoprzecinkowych lub „soft-float”.

Poza tym: Powód oddzielenia funkcji matematycznych jest historyczny i miał na celu jedynie zmniejszenie wielkości programów wykonywalnych w bardzo starych systemach uniksowych, być może zanim udostępnione biblioteki będą dostępne, o ile mi wiadomo.

Teraz kompilator może zoptymalizować standardową funkcję biblioteki C sin()(dostarczoną przez libm.so), aby zastąpić ją wywołaniem instrukcji natywnej do wbudowanej funkcji sin () procesora / FPU, która istnieje jako instrukcja FPU ( FSINdla x86 / x87) na nowsze procesory, takie jak seria Core 2 (jest to poprawne już od czasów i486DX). Zależy to od flag optymalizacji przekazywanych do kompilatora gcc. Gdyby kompilatorowi nakazano napisać kod, który byłby wykonywany na dowolnym procesorze i386 lub nowszym, nie dokonałby takiej optymalizacji. -mcpu=486Flag kompilatora by poinformować, że to bezpieczne, aby dokonać takiej optymalizacji.

Teraz, jeśli program wykonałby wersję oprogramowania funkcji sin (), zrobiłby to w oparciu o algorytm CORDIC (komputer cyfrowy z rotacją współrzędnych) lub BKM , lub bardziej prawdopodobne jest obliczenie tabeli lub szeregu mocy, które jest obecnie powszechnie używane do obliczania takie funkcje transcendentalne. [Src: http://en.wikipedia.org/wiki/Cordic#Application]

Każda najnowsza (od około 2,9x) wersja gcc oferuje również wbudowaną wersję sin, __builtin_sin()która zostanie użyta do zastąpienia standardowego wywołania do wersji biblioteki C, jako optymalizacja.

Jestem pewien, że jest to jasne jak błoto, ale mam nadzieję, że dostarczy ci więcej informacji, niż się spodziewałeś, i wiele skoków do punktów, aby dowiedzieć się więcej.

Jeśli chcesz spojrzeć na faktyczną implementację GNU tych funkcji w C, sprawdź najnowszy pakiet glibc. Zobacz biblioteka GNU C .

Nie używaj serii Taylor. Wielomiany Czebyszewa są zarówno szybsze, jak i dokładniejsze, jak wskazało kilka osób powyżej. Oto implementacja (oryginalnie z ZX Spectrum ROM): https://albertveli.wordpress.com/2015/01/10/zx-sine/

Obliczanie sinus / cosinus / tangens jest w rzeczywistości bardzo łatwe do zrobienia za pomocą kodu za pomocą serii Taylora. Samo napisanie zajmuje około 5 sekund.

Cały proces można podsumować za pomocą tego równania:

Oto kilka procedur, które napisałem dla 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));  
}

Poprawiona wersja kodu z odpowiedzi 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);
}

Istota tego, jak to robi, leży w tym fragmencie Applied Numerical Analysis Geralda Wheatleya:

Kiedy twój program prosi komputer o uzyskanie wartości lub , czy zastanawiałeś się, w jaki sposób można uzyskać te wartości, jeśli najpotężniejszymi funkcjami, jakie może obliczyć, są wielomiany? Nie wyszukuje ich w tabelach i interpoluje! Zamiast tego komputer aproksymuje każdą funkcję inną niż wielomiany z jakiegoś wielomianu, który jest dostosowany do bardzo dokładnych wartości.

Kilka punktów, o których należy wspomnieć powyżej, to fakt, że niektóre algorytmy faktycznie interpolują dane z tabeli, chociaż tylko dla kilku pierwszych iteracji. Zwróć też uwagę na to, jak wspomniano, że komputery używają wielomianów aproksymacyjnych bez określania, jaki typ wielomianu aproksymacyjnego. Jak zauważyli inni w wątku, wielomiany Czebyszewa są w tym przypadku bardziej wydajne niż wielomiany Taylora.

jeśli chcesz, sinto

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

jeśli chcesz, costo

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

jeśli chcesz, sqrtto

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

więc po co używać niedokładnego kodu, gdy zrobią to instrukcje maszyny?