O szybszym zbliżeniu log (x)


10

Napisałem kodu jednocześnie temu, która polegała obliczyć bez korzystania z funkcji biblioteki. Wczoraj sprawdzałem stary kod i starałem się, aby był jak najszybciej (i poprawiony). Oto moja dotychczasowa próba:log(x)

const double ee = exp(1);

double series_ln_taylor(double n){ /* n = e^a * b, where a is an non-negative integer */
    double lgVal = 0, term, now;
    int i, flag = 1;

    if ( n <= 0 ) return 1e-300;
    if ( n * ee < 1 )
        n = 1.0 / n, flag = -1; /* for extremely small n, use e^-x = 1/n */

    for ( term = 1; term < n ; term *= ee, lgVal++ );
    n /= term;

    /* log(1 - x) = -x - x**2/2 - x**3/3... */
    n = 1 - n;
    now = term = n;
    for ( i = 1 ; ; ){
        lgVal -= now;
        term *= n;
        now = term / ++i;
        if ( now < 1e-17 ) break;
    }

    if ( flag == -1 ) lgVal = -lgVal;

    return lgVal;
}

Tutaj staram się znaleźć tak, e właśnie nad n, a następnie dodać wartość logarytmu z naea , który jest mniejszy niż 1. W tym momencie, dla ekspansji TayloralOg(1-x)może być stosowany bez obawy.nealog(1  x)

Ostatnio zainteresowałem się analizą numeryczną i dlatego nie mogę przestać zadawać pytania, o ile szybciej ten segment kodu można uruchomić w praktyce, a jednocześnie jest wystarczająco poprawny? Czy muszę przejść na inne metody, na przykład używając ciągłego ułamka, takiego jak ten ?

funkcję dostarczaną przez C standardowej biblioteki jest prawie 5.1 razy szybciej niż w tym wykonaniu.log(x)

AKTUALIZACJA 1 : Korzystając z hiperbolicznej serii arctan wspomnianej w Wikipedii , obliczenia wydają się prawie 2,2 razy wolniejsze niż standardowa funkcja dziennika biblioteki C. Chociaż nie sprawdziłem dokładnie wydajności, a w przypadku większych liczb moja obecna implementacja wydaje się NAPRAWDĘ powolna. Chcę sprawdzić zarówno moją implementację pod kątem związanego z błędem, jak i średniego czasu dla szerokiego zakresu liczb, jeśli mogę. Oto mój drugi wysiłek.

double series_ln_arctanh(double n){ /* n = e^a * b, where a is an non-negative integer */
    double lgVal = 0, term, now, sm;
    int i, flag = 1;

    if ( n <= 0 ) return 1e-300;
    if ( n * ee < 1 ) n = 1.0 / n, flag = -1; /* for extremely small n, use e^-x = 1/n */

    for ( term = 1; term < n ; term *= ee, lgVal++ );
    n /= term;

    /* log(x) = 2 arctanh((x-1)/(x+1)) */
    n = (1 - n)/(n + 1);

    now = term = n;
    n *= n;
    sm = 0;
    for ( i = 3 ; ; i += 2 ){
        sm += now;
        term *= n;
        now = term / i;
       if ( now < 1e-17 ) break;
    }

    lgVal -= 2*sm;

    if ( flag == -1 ) lgVal = -lgVal;
    return lgVal;
}

Każda sugestia lub krytyka jest mile widziana.

1e81e3084e15

double series_ln_better(double n){ /* n = e^a * b, where a is an non-negative integer */
    double lgVal = 0, term, now, sm;
    int i, flag = 1;

    if ( n == 0 ) return -1./0.; /* -inf */
    if ( n < 0 ) return 0./0.;   /* NaN*/
    if ( n < 1 ) n = 1.0 / n, flag = -1; /* for extremely small n, use e^-x = 1/n */

    /* the cutoff iteration is 650, as over e**650, term multiplication would
       overflow. For larger numbers, the loop dominates the arctanh approximation
       loop (with having 13-15 iterations on average for tested numbers so far */

    for ( term = 1; term < n && lgVal < 650 ; term *= ee, lgVal++ );
    if ( lgVal == 650 ){
        n /= term;
        for ( term = 1 ; term < n ; term *= ee, lgVal++ );
    }
    n /= term;

    /* log(x) = 2 arctanh((x-1)/(x+1)) */
    n = (1 - n)/(n + 1);

    now = term = n;
    n *= n;
    sm = 0;

    /* limiting the iteration for worst case scenario, maximum 24 iteration */
    for ( i = 3 ; i < 50 ; i += 2 ){
        sm += now;
        term *= n;
        now = term / i;
        if ( now < 1e-17 ) break;
    }

    lgVal -= 2*sm;

    if ( flag == -1 ) lgVal = -lgVal;

    return lgVal;
}

Odpowiedzi:


17

To nie jest tak naprawdę autorytatywna odpowiedź, a raczej lista problemów, które moim zdaniem powinieneś rozważyć, a ja nie przetestowałem twojego kodu.

log2.15.1

f(x)doublen12

n1.7976e+308term=infn=11017nterm *= e709.78266108405500745

1030000

Podejrzewam, że przy odrobinie wysiłku możesz poświęcić część tej solidności na wydajność, np. Ograniczając zakres argumentów lub zwracając nieco mniej dokładne wyniki.

3. Wydajność tego rodzaju kodu może w dużym stopniu zależeć od architektury procesora, na której działa. To głęboki i zaangażowany temat, ale producenci procesorów, tacy jak Intel, publikują przewodniki dotyczące optymalizacji, które wyjaśniają różne interakcje między twoim kodem a procesorem, na którym działa. Buforowanie może być stosunkowo proste, ale rzeczy takie jak przewidywanie rozgałęzień, równoległość na poziomie instrukcji i utknięcie potoku z powodu zależności danych są trudne do zobaczenia w kodzie wysokiego poziomu, ale mają duże znaczenie dla wydajności.

x~y~=f~(x~)y=f(x~) jakdokładny?). To nie jest to samo, co pokazanie, że szereg Taylora jest zbieżny z powodu obecności błędów zaokrąglania zmiennoprzecinkowego.

4.5 Dobrym sposobem na przetestowanie niesprawdzonej funkcji pod kątem dokładności byłoby oszacowanie jej dla każdego z czterech miliardów (mniej, jeśli poprawnie redukujesz argumenty, jak tutaj), liczby zmiennoprzecinkowe pojedynczej precyzji i porównanie błędów ze standardowym logem z libm. To zajmuje trochę czasu, ale przynajmniej jest dokładne.

5. Ponieważ od samego początku znasz dokładność podwójnych, nie musisz mieć nieograniczonej pętli: liczbę iteracji można ustalić z góry (prawdopodobnie około 50). Użyj tego, aby usunąć gałęzie z kodu lub przynajmniej z góry ustawić liczbę iteracji.

Obowiązują również wszystkie zwykłe pomysły dotyczące rozwijania pętli.

6. Możliwe jest zastosowanie technik aproksymacyjnych innych niż seria Taylora. Istnieją również serie Czebiszewa (z powtarzaniem Clenshawa), przybliżenia Pade, a czasem metody znajdowania korzeni, takie jak metoda Newtona, ilekroć twoja funkcja może zostać przekształcona jako pierwiastek prostszej funkcji (np. Słynna sztuczka sqrt ).

Ciągłe ułamki prawdopodobnie nie będą zbyt duże, ponieważ obejmują podział, który jest znacznie droższy niż pomnażanie / dodawanie. Jeśli spojrzeć _mm_div_ssna https://software.intel.com/sites/landingpage/IntrinsicsGuide/ , podział ma opóźnienie 13-14 cykli i przepustowość 5-14, w zależności od architektury, w porównaniu z 3-5 / 0,5-1 dla pomnóż / dodaj / madd. Tak więc ogólnie (nie zawsze) sensowne jest dążenie do wyeliminowania podziałów w jak największym stopniu.

Niestety matematyka nie jest tutaj tak świetnym przewodnikiem, ponieważ wyrażenia o krótkich formułach niekoniecznie są najszybszymi. Na przykład matematyka nie karze podziałów.

x=m×2em12<m1exfrexp

8. Porównaj swoje logz logw libmlub openlibm(np .: https://github.com/JuliaLang/openlibm/blob/master/src/e_log.c ). Jest to zdecydowanie najłatwiejszy sposób, aby dowiedzieć się, co inni ludzie już odkryli. Istnieją również specjalnie zoptymalizowane wersje libm specyficzne dla producentów procesorów, ale te zwykle nie mają opublikowanego kodu źródłowego.

Boost :: sf ma pewne funkcje specjalne, ale nie podstawowe. Warto jednak przyjrzeć się źródłu log1p: http://www.boost.org/doc/libs/1_58_0/libs/math/doc/html/math_toolkit/powers/log1p.html

Istnieją również biblioteki arytmetyczne o dowolnej precyzji typu open source, takie jak mpfr, które mogą wykorzystywać inne algorytmy niż libm ze względu na wymaganą wyższą precyzję.

9. Dokładność i stabilność algorytmów numerycznych Highama jest dobrym wstępem do analizy błędów algorytmów numerycznych na wyższym poziomie. W przypadku samych algorytmów aproksymacyjnych dobrym przykładem jest teoria aproksymacji. Praktyka aproksymacji autorstwa Trefethena.

10. Wiem, że mówi się to zbyt często, ale dość duże projekty oprogramowania rzadko zależą od czasu wykonywania jednej małej funkcji wywoływanej w kółko. Nie tak często trzeba się martwić o wydajność dziennika, chyba że profilowałeś swój program i upewniłeś się, że jest to ważne.


26414e15

1.13e13term

 1e8

1
k=11071lnk

2
frexp x=m×2elnx=eln2+lnm

5

Odpowiedź Kirilla dotyczyła już wielu istotnych kwestii. Chciałbym rozwinąć niektóre z nich w oparciu o praktyczne doświadczenia w projektowaniu bibliotek matematycznych. Uwaga z góry: projektanci bibliotek matematycznych zwykle używają każdej opublikowanej optymalizacji algorytmicznej, a także wielu optymalizacji specyficznych dla maszyny, z których nie wszystkie zostaną opublikowane. Kod często będzie pisany w języku asemblera, zamiast używać skompilowanego kodu. Jest zatem mało prawdopodobne, aby prosta i skompilowana implementacja osiągnęła ponad 75% wydajności istniejącej wysokiej jakości implementacji biblioteki matematycznej, przy założeniu identycznych zestawów funkcji (dokładność, obsługa przypadków specjalnych, raportowanie błędów, obsługa trybu zaokrąglania).

explogerfcΓ

Dokładność jest zazwyczaj oceniana przez porównanie z referencją o wyższej precyzji (innej firmy). Funkcje pojedynczej argumentacji o pojedynczej precyzji można łatwo przetestować wyczerpująco, inne funkcje wymagają testowania za pomocą (ukierunkowanych) losowych wektorów testowych. Oczywiście nie można obliczyć nieskończenie dokładnych wyników odniesienia, ale badania dylematu twórcy tabeli sugerują, że dla wielu prostych funkcji wystarczy obliczyć odniesienie z dokładnością około trzykrotnie większą niż precyzja docelowa. Zobacz na przykład:

Vincent Lefèvre, Jean-Michel Muller, „Najgorsze przypadki prawidłowego zaokrąglania funkcji elementarnych w podwójnej precyzji”. W postępowaniu 15 Sympozjum IEEE na temat arytmetyki komputerowej , 2001,111-118). (preprint online)

Pod względem wydajności należy odróżnić optymalizację pod kątem opóźnień (ważne, gdy patrzy się na czas wykonywania operacji zależnych) od optymalizacji pod kątem przepustowości (istotnej przy rozważaniu czasu wykonania niezależnych operacji). W ciągu ostatnich dwudziestu lat rozpowszechnianie sprzętowych technik paralelizacji, takich jak równoległość na poziomie instrukcji (np. Procesory superskalarne, procesory poza kolejnością), równoległość na poziomie danych (np. Instrukcje SIMD) i równoległość na poziomie wątku (np. Hiperwątkowość, procesory wielordzeniowe) doprowadziły do ​​położenia nacisku na przepustowość obliczeniową jako bardziej odpowiednią metrykę.

log(1+x)=p(x)log(x)=2atanh((x1)/(x+1))=p(((x1)/(x+1))2)p

Połączona operacja wielokrotnego dodawania ( FMA ), po raz pierwszy wprowadzona przez IBM 25 lat temu, a teraz dostępna we wszystkich głównych architekturach procesorów, jest kluczowym elementem składowym nowoczesnych implementacji bibliotek matematycznych. Zapewnia redukcję błędów zaokrąglania, daje ograniczoną ochronę przed odejmowaniem anulowania i znacznie upraszcza arytmetykę podwójnego podwójnego .

C99log()C99fma()233

#include <math.h>

/* compute natural logarithm

   USE_ATANH == 1: maximum error found: 0.83482 ulp @ 0.7012829191167614
   USE_ATANH == 0: maximum error found: 0.83839 ulp @ 1.2788954397331760
*/
double my_log (double a)
{
    const double LOG2_HI = 0x1.62e42fefa39efp-01; // 6.9314718055994529e-01
    const double LOG2_LO = 0x1.abc9e3b39803fp-56; // 2.3190468138462996e-17
    double m, r, i, s, t, p, f, q;
    int e;

    m = frexp (a, &e);
    if (m < 0.70703125) { // 181/256
        m = m + m;
        e = e - 1;
    }
    i = (double)e;

    /* m in [181/256, 362/256] */

#if USE_ATANH
    /* Compute q = (m-1) / (m+1) */
    p = m + 1.0;
    m = m - 1.0;
    q = m / p;

    /* Compute (2*atanh(q)/q-2*q) as p(q**2), q in [-75/437, 53/309] */
    s = q * q;
    r =             0x1.2f1da230fb057p-3;  // 1.4800574027992994e-1
    r = fma (r, s,  0x1.399f73f934c01p-3); // 1.5313616375223663e-1
    r = fma (r, s,  0x1.7466542530accp-3); // 1.8183580149169243e-1
    r = fma (r, s,  0x1.c71c51a8bf129p-3); // 2.2222198291991305e-1
    r = fma (r, s,  0x1.249249425f140p-2); // 2.8571428744887228e-1
    r = fma (r, s,  0x1.999999997f6abp-2); // 3.9999999999404662e-1
    r = fma (r, s,  0x1.5555555555593p-1); // 6.6666666666667351e-1
    r = r * s;

    /* log(a) = 2*atanh(q) + i*log(2) = LOG2_LO*i + p(q**2)*q + 2q + LOG2_HI*i.
       Use K.C. Ng's trick to improve the accuracy of the computation, like so:
       p(q**2)*q + 2q = p(q**2)*q + q*t - t + m, where t = m**2/2.
    */
    t = m * m * 0.5;
    r = fma (q, t, fma (q, r, LOG2_LO * i)) - t + m;
    r = fma (LOG2_HI, i, r);

#else // USE_ATANH

    /* Compute f = m -1 */
    f = m - 1.0;
    s = f * f;

    /* Approximate log1p (f), f in [-75/256, 106/256] */
    r = fma (-0x1.961d64ddd82b6p-6, f, 0x1.d35fd598b1362p-5); // -2.4787281515616676e-2, 5.7052533321928292e-2
    t = fma (-0x1.fcf5138885121p-5, f, 0x1.b97114751d726p-5); // -6.2128580237329929e-2, 5.3886928516403906e-2
    r = fma (r, s, t);
    r = fma (r, f, -0x1.b5b505410388dp-5); // -5.3431043874398211e-2
    r = fma (r, f,  0x1.dd660c0bd22dap-5); //  5.8276198890387668e-2
    r = fma (r, f, -0x1.00bda5ecdad6fp-4); // -6.2680862565391612e-2
    r = fma (r, f,  0x1.1159b2e3bd0dap-4); //  6.6735934054864471e-2
    r = fma (r, f, -0x1.2489f14dd8883p-4); // -7.1420614809115476e-2
    r = fma (r, f,  0x1.3b0ee248a0ccfp-4); //  7.6918491287915489e-2
    r = fma (r, f, -0x1.55557d3b497c3p-4); // -8.3333481965921982e-2
    r = fma (r, f,  0x1.745d4666f7f48p-4); //  9.0909266480136641e-2
    r = fma (r, f, -0x1.999999d959743p-4); // -1.0000000092767629e-1
    r = fma (r, f,  0x1.c71c70bbce7c2p-4); //  1.1111110722131826e-1
    r = fma (r, f, -0x1.fffffffa61619p-4); // -1.2499999991822398e-1
    r = fma (r, f,  0x1.249249262c6cdp-3); //  1.4285714290377030e-1
    r = fma (r, f, -0x1.555555555f03cp-3); // -1.6666666666776730e-1
    r = fma (r, f,  0x1.999999999759ep-3); //  1.9999999999974433e-1
    r = fma (r, f, -0x1.fffffffffff53p-3); // -2.4999999999999520e-1
    r = fma (r, f,  0x1.555555555555dp-2); //  3.3333333333333376e-1
    r = fma (r, f, -0x1.0000000000000p-1); // -5.0000000000000000e-1

    /* log(a) = log1p (f) + i * log(2) */
    p = fma ( LOG2_HI, i, f);
    t = fma (-LOG2_HI, i, p);
    f = fma ( LOG2_LO, i, f - t);
    r = fma (r, s, f);
    r = r + p;
#endif // USE_ATANH

    /* Handle special cases */
    if (!((a > 0.0) && (a <= 0x1.fffffffffffffp1023))) {
        r = a + a;  // handle inputs of NaN, +Inf
        if (a  < 0.0) r =  0.0 / 0.0; //  NaN
        if (a == 0.0) r = -1.0 / 0.0; // -Inf
    }
    return r;
}

(+1) Czy wiesz, czy popularne implementacje typu open source (takie jak openlibm) są tak dobre, jak mogą być, czy też można ulepszyć ich funkcje specjalne?
Kirill

1
@Kirill Ostatnio przyglądałem się implementacjom open source (wiele lat temu), nie wykorzystywały one zalet FMA. W tamtym czasie IBM Power i Intel Itanium były jedynymi architekturami, które obejmowały tę operację, teraz wsparcie sprzętowe jest wszechobecne. Również wtedy aproksymacje między tabelą a wielomianem były wówczas bardzo aktualne, teraz tabele są niekorzystne: dostęp do pamięci powoduje większe zużycie energii, mogą (i robią) zakłócać wektoryzację, a przepustowość obliczeniowa wzrosła bardziej niż przepustowość pamięci powodując potencjalny negatywny wpływ na wydajność tabel.
njuffa,
Korzystając z naszej strony potwierdzasz, że przeczytałeś(-aś) i rozumiesz nasze zasady używania plików cookie i zasady ochrony prywatności.
Licensed under cc by-sa 3.0 with attribution required.