Jak C oblicza sin () i inne funkcje matematyczne?

248

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ą).

Motek
źródło
2
Należy pamiętać, że ta implementacja zależy. Powinieneś określić, która implementacja najbardziej Cię interesuje.
jason
3
Oznaczyłem .NET i C, ponieważ szukałem w obu miejscach i nie mogłem się zorientować. Chociaż patrząc na demontaż .NET wygląda na to, że może wywoływać niezarządzane C, o ile wiem, że mają tę samą implementację.
Hank

Odpowiedzi:

213

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

Jason Orendorff
źródło
35
Aby zobaczyć, że tak naprawdę jest to kod działający na x86: skompiluj program, który wywołuje sin(); wpisz gdb a.out, a break sinnastępnie run, a następnie disassemble.
Jason Orendorff
5
@Henry: nie popełnij jednak błędu, myśląc, że to dobry kod. To naprawdę okropne , nie ucz się kodować w ten sposób!
Thomas Bonini
2
@Andreas Hmm, masz rację, kod IBM wygląda dość okropnie w porównaniu z fdlibm. Zredagowałem odpowiedź, aby dodać linki do procedury sinusoidalnej fdlibm.
Jason Orendorff
3
@Henry: __kernel_sinjest jednak zdefiniowany w k_sin.c i jest czysty C. Kliknij go jeszcze raz - za pierwszym razem spartaczyłem adres URL.
Jason Orendorff,
3
Powiązany kod sysdeps jest szczególnie interesujący, ponieważ jest poprawnie zaokrąglony. To znaczy, najwyraźniej daje najlepszą możliwą odpowiedź dla wszystkich wartości wejściowych, co stało się możliwe całkiem niedawno. W niektórych przypadkach może to być powolne, ponieważ może być konieczne obliczenie wielu dodatkowych cyfr w celu zapewnienia prawidłowego zaokrąglenia. W innych przypadkach jest niezwykle szybki - w przypadku wystarczająco małych liczb odpowiedź jest tylko kątem.
Bruce Dawson
66

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 .

John D. Cook
źródło
10
Transcendentalne funkcje matematyczne, takie jak sinus i cosinus, mogą być zaimplementowane w mikrokodzie lub jako instrukcje sprzętowe w obecnych 32-bitowych procesorach stacjonarnych i serwerowych. Nie zawsze tak było, dopóki w i486 (DX) wszystkie obliczenia zmiennoprzecinkowe nie były wykonywane w oprogramowaniu („soft-float”) dla serii x86 bez osobnego koprocesora. Nie wszystkie z nich (FPU) obejmowały funkcje transcendentalne (np. Weitek 3167).
mctylr
1
Czy mógłbyś to sprecyzować? Jak „dopracować” przybliżenie za pomocą serii Taylora?
Hank
4
Jeśli chodzi o „dopracowanie” odpowiedzi, załóżmy, że obliczasz zarówno sinus, jak i cosinus. Załóżmy, że znasz dokładną wartość obu w jednym punkcie (np. Z CORDIC), ale chcesz tę wartość w pobliskim punkcie. Następnie, dla małej różnicy h, możesz zastosować przybliżenia Taylora f (x + h) = f (x) + h f '(x) lub f (x + h) = f (x) + h f' (x) + h ^ 2 f '' (x) / 2.
John D. Cook
6
Układy x86 / x64 mają instrukcję asemblera do obliczania sinusa (fsin), ale ta instrukcja bywa czasami niedokładna i dlatego jest rzadko używana. Szczegółowe informacje można znaleźć na stronie randomascii.wordpress.com/2014/10/09/ ... Większość innych procesorów nie ma instrukcji dla sinusa i cosinusa, ponieważ ich obliczanie w oprogramowaniu daje większą elastyczność, a może nawet jest szybsze.
Bruce Dawson
3
Kordowe elementy w chipach Intel na ogół NIE są używane. Po pierwsze, dokładność i rozdzielczość operacji jest niezwykle ważna dla wielu aplikacji. Cordic jest notorycznie niedokładny, gdy dojdziesz do siódmej cyfry, i jest nieprzewidywalny. Po drugie, słyszałem, że w ich implementacji występuje błąd, który powoduje jeszcze więcej problemów. Przyjrzałem się funkcji sin dla Linuksa gcc i oczywiście używa Chebyshev. wbudowane rzeczy nie są używane. Och, także algorytm kordowy w układzie jest wolniejszy niż rozwiązanie programowe.
Donald Murray,
63

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ę.

Donald Murray
źródło
6
Nie mogłem się bardziej zgodzić z kilkoma pierwszymi zdaniami. Warto również przypomnieć, że obliczanie specjalnych funkcji z gwarantowaną precyzją jest trudnym problemem . Mądrzy ludzie, o których wspominasz, spędzają na tym większość życia. Ponadto, bardziej technicznie, wielomian min-max są poszukiwanymi graalami, a wielomiany Czebyszewa są dla nich prostszymi proxy.
Alexandre C.
161
-1 za nieprofesjonalny i chaotyczny (i nieco niegrzeczny) ton oraz za fakt, że faktyczna, niepotrzebna treść tej odpowiedzi, pozbawiona bełkotu i protekcjonalności, sprowadza się w zasadzie do: „Często używają wielomianów Czebyszewa; zobacz tę książkę po więcej szczegółów, to naprawdę dobrze! ” Co, jak wiecie, może być absolutnie poprawne, ale tak naprawdę nie jest to taka samodzielna odpowiedź, jakiej chcemy tutaj na SO. Skondensowane w ten sposób byłoby jednak przyzwoitym komentarzem do tego pytania.
Ilmari Karonen,
2
We wczesnych latach tworzenia gier zwykle odbywało się to przy użyciu tabel przeglądowych, które wymagały szybkości). Zazwyczaj nie używaliśmy standardowych funkcji lib do tych rzeczy.
topspin
4
Często używam tabel odnośników w systemach osadzonych i Bittians (zamiast radianów), ale jest to do specjalistycznej aplikacji (takiej jak twoje gry). Myślę, że facet jest zainteresowany tym, jak kompilator c oblicza grzech dla liczb zmiennoprzecinkowych ....
Donald Murray,
1
Ach, 50 lat temu. Zacząłem grać z takimi na Burroughs B220 z serią McLaren. Później sprzęt CDC, a potem Motorola 68000. Arcsin był nieuporządkowany - wybrałem iloraz dwóch wielomianów i opracowałem kod, aby znaleźć optymalne współczynniki.
Rick James
15

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;

}
Blindy
źródło
1
Jeśli nieco poprawisz współczynniki (i kodujesz je na wielomian), możesz zatrzymać około 2 iteracji wcześniej.
Rick James
14

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ń.

Mehrdad Afshari
źródło
12
Prawdziwa implementacja prawdopodobnie nie użyje serii Taylora, ponieważ istnieją bardziej wydajne sposoby. Musisz tylko poprawnie aproksymować w domenie [0 ... pi / 2], a istnieją funkcje, które zapewnią dobre przybliżenie bardziej efektywnie niż seria Taylora.
David Thornley,
2
@David: Zgadzam się. Byłem wystarczająco ostrożny, by w mojej odpowiedzi wspomnieć słowo „lubię”. Ale ekspansja Taylora jest prosta, aby wyjaśnić ideę metod przybliżających funkcje. To powiedziawszy, widziałem implementacje oprogramowania (nie jestem pewien, czy zostały zoptymalizowane), które korzystały z serii Taylor.
Mehrdad Afshari
1
W rzeczywistości aproksymacje wielomianowe są jednym z najbardziej wydajnych sposobów obliczania funkcji trygonometrycznych.
Jeremy Salwen,
13

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 )

wyjaśnienie

Hannoun Yassir
źródło
2
Czy wiesz, że możesz używać interfejsu API Google Chart do tworzenia takich formuł za pomocą TeXa? code.google.com/apis/chart/docs/gallery/formulas.html
Gab Royer
11

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.

Thomas Pornin
źródło
3
Rejestry SSE2 nie są używane do obliczania sin (), ani w trybie x86, ani w trybie x64, i oczywiście sin jest obliczany sprzętowo niezależnie od trybu. Hej, w 2010 roku mieszkamy :)
Igor Korkhov
7
@Igor: to zależy od biblioteki matematyki, na którą patrzysz. Okazuje się, że najbardziej zoptymalizowane biblioteki matematyczne na x86 używają implementacji oprogramowania SSE sini cossą szybsze niż instrukcje sprzętowe na FPU. Prostsze, bardziej naiwne biblioteki zwykle używają instrukcji fsini fcos.
Stephen Canon
@Stephen Canon: Czy te szybkie biblioteki mają 80-bitową precyzję, podobnie jak rejestry FPU? Mam bardzo podstępne podejrzenie, że preferują szybkość zamiast precyzji, co oczywiście jest uzasadnione w wielu scenariuszach, na przykład w grach. I wierzę, że obliczanie sinusa z 32-bitową precyzją przy użyciu SSE i wstępnie obliczonych tabel pośrednich może być szybsze niż przy użyciu FSINpełnej precyzji. Byłbym bardzo wdzięczny, gdybyś podał mi nazwy tych szybkich bibliotek, ciekawe jest zajrzeć.
Igor Korkow
@Igor: na x86 w trybie 64-bitowym, przynajmniej na wszystkich znanych mi systemach uniksowych, precyzja jest ograniczona do 64 bitów, a nie 79 bitów FPU x87. Implementacja oprogramowania sin()jest około dwa razy szybsza niż to, co fsinoblicza (właśnie dlatego, że odbywa się to z mniejszą precyzją). Zauważ, że x87 ma nieco mniej dokładną dokładność niż zapowiadane 79 bitów.
Thomas Pornin
1
Rzeczywiście, zarówno 32-bitowe, jak i 64-bitowe implementacje sin () w bibliotekach wykonawczych msvc nie używają instrukcji FSIN. W rzeczywistości dają różne wyniki, na przykład grzech (0.70444454416678126). Spowoduje to 0,64761068800896837 (z tolerancją 0,5 * (eps / 2)) w programie 32-bitowym i spowoduje 0,64761068800896848 (źle) w programie 64-bitowym.
e.tadeu
9

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.

gnasher729
źródło
1
Jest to ładne wytłumaczenie jak jeden może obliczyć sin (x) sprawnie, ale to naprawdę nie wydaje się, aby odpowiedzieć na pytanie PO, która jest specjalnie o tym, jak wspólna bibliotek C / kompilatory należy go obliczyć.
Ilmari Karonen,
Wielomiany Czebyszewa minimalizują maksymalną wartość bezwzględną w danym przedziale, ale nie minimalizują największej różnicy między funkcją docelową a wielomianem. Robią to wielomiany Minimax.
Eric Postpischil,
9

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 ()

chux - Przywróć Monikę
źródło
6

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.

John Bode
źródło
5

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.

mnemosyn
źródło
5

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.

Norman Ramsey
źródło
5

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

Moschops
źródło
4

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.

James
źródło
4

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.

Igor Korkow
źródło
4

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.

Mctylr
źródło
3

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

Chris Tonkinson
źródło
3

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/

Albert Veli
źródło
2
To naprawdę nie odpowiada na zadane pytanie. OP pyta, w jaki sposób funkcje wyzwalania obliczane przez popularne kompilatory / biblioteki C (i jestem prawie pewien, że ZX Spectrum się nie kwalifikuje), a nie jak powinny być obliczane. Może to być jednak przydatny komentarz do niektórych wcześniejszych odpowiedzi.
Ilmari Karonen,
1
Ach, masz rację. Powinien to być komentarz, a nie odpowiedź. Od dawna nie używałem SO i zapomniałem, jak działa system. W każdym razie myślę, że implementacja Spectrum jest istotna, ponieważ miała naprawdę wolny procesor, a szybkość była najważniejsza. Najlepszy algorytm jest więc z pewnością całkiem niezły, dlatego dobrym pomysłem byłoby, aby biblioteki C implementowały funkcje triggera przy użyciu wielomianów Czebyszewa.
Albert Veli
2

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:

ekspansja grzechu i kosztów

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));  
}
użytkownik1432532
źródło
4
Jest to raczej zła implementacja, ponieważ nie korzysta z tego, że kolejne terminy z serii sinus i cosinus mają bardzo proste ilorazy. Co oznacza, że ​​można zmniejszyć liczbę multiplikacji i podziałów z O (n ^ 2) tutaj na O (n). Dalsze redukcje są osiągane przez zmniejszenie o połowę i podniesienie do kwadratu, jak na przykład w bibliotece matematycznej bc (kalkulator wieloprecyzyjny POSIX).
Lutz Lehmann
2
Wydaje się również, że nie odpowiada na zadane pytanie; OP pyta, w jaki sposób funkcje wyzwalania są obliczane przez popularne kompilatory / biblioteki C, a nie dla niestandardowych ponownych wdrożeń.
Ilmari Karonen,
2
Myślę, że to dobra odpowiedź, ponieważ odpowiada duchowi pytania, które (i mogę się tylko domyślać) ciekawość o „czarnej skrzynce”, która w przeciwnym razie działałaby jak sin (). Jest to jedyna odpowiedź, która daje szansę na szybkie zrozumienie, co się dzieje, przez przerzucenie go w ciągu kilku sekund zamiast odczytywania zoptymalizowanego kodu źródłowego C.
Mike M.
w rzeczywistości biblioteki używają znacznie bardziej zoptymalizowanej wersji, uświadamiając sobie, że kiedy już masz pewien termin, możesz uzyskać następny, mnożąc niektóre wartości. Zobacz przykład w odpowiedzi Blindy .
Obliczasz
0

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);
}
Rugnar
źródło
0

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 wprowadź opis zdjęcia tutajlub wprowadź opis zdjęcia tutaj, 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.

Dean P.
źródło
-1

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?

użytkownik80998
źródło
4
Być może dlatego, że instrukcje maszyny są również bardzo niedokładne.
Ilmari Karonen,