C++ Pi Calculator – Formula Leibniz (Revizuirea codului, C++, Performanță, Multithreading, Windows, Metode Numerice)

Nikko77 a intrebat.

Ieri am văzut un CodeTrain video în care Daniel Shiffman a încercat să aproximeze Pi folosind formula seria Leibniz. A fost interesant, așa că am decis să încerc și eu.

Am scris această aplicație de consolă cu un scop simplu: Să aproximeze cât mai repede câte cifre din Pi poate programul până când ajunge la numărul maxim pe care îl poate atinge contorul de cicluri(fără să înceapă să comită erori de calcul cu formula seriei) SAU programul se oprește.

Acesta începe prin a rula mai întâi două fire de execuție în fundal care execută două funcții care fac toată treaba: piLoop(), , pentru calcularea seriilor; și displayLoop(), , pentru reîmprospătarea ecranului și afișarea Pi calculat efectiv. Pentru a preveni problemele de concurență, se utilizează mutexuri pentru a împiedica firele să acceseze Pi și alte variabile partajate ale programului în același timp. Programul se execută până când piCycles atinge valoarea maximă a long long int poate obține SAU se apasă tasta End.

Programul final rulează, aruncând o privire (ceea ce înseamnă o măsurătoare foarte puțin fiabilă), la o rată aproximativă de 30 de milioane de cicluri pe secundă în calculatorul meu. Aceasta înseamnă că poate calcula a zecea cifră a PI (aproximativ 87 de miliarde de cicluri) în aproximativ o oră.

Vreau să știu dacă există unele optimizări care pot fi făcute pentru ca programul să poată rula la toată puterea lui.

main.cpp

#include <iostream>
#include <thread>
#include <mutex>
#include <chrono>
#include <windows.h>
#include <limits>

void piLoop();
void displayLoop();
void getShouldStop();

namespace {
    typedef std::numeric_limits<long double> ldbl;
    typedef std::numeric_limits<long long int> ullint;

    const char displayDigits = ldbl::max_digits10;
    const short int updateDelay_ms = 50; // You can modify this value to change the speed of console refreshing.
    const long long int maxCicles = ullint::max();

    bool shouldStop = false;
    long double PI = 0;
    unsigned long long int piCicles = 0;
    std::mutex piLock;
}

int main()
{
    std::thread piCalculatorThread( piLoop );
    std::thread displayUpdaterThread( displayLoop );
    piCalculatorThread.join();
    displayUpdaterThread.join();

    std::cin.get();

    return 0;
}

void getShouldStop() {
    shouldStop = GetAsyncKeyState( VK_END ) || shouldStop;
}

// If PI mutex is not locked, proceed to calculate the next PI, else wait
void piLoop() {
    while ( true ) {
        piLock.lock();

        if ( shouldStop ) {
            piLock.unlock();
            return;
        }

        // Here the series is calculated
        PI += ((long double)((piCicles&1)>0 ? -1 : 1)) / (piCicles*2 + 1);

        if ( piCicles == maxCicles ) {
            shouldStop = true;
            piLock.unlock();
            return;
        }

        piCicles++;

        piLock.unlock();
    }
}

// Wait until PI mutex is unlocked, then refresh the screen.
void displayLoop() {
    std::cout.precision( displayDigits );
    while ( true ) {
        piLock.lock();

        system("cls");
        std::cout << "

  ----- PI Calculator ----- 

  - 
Seeing how PI is calculated is more enjoyable while eating some snacks.
";
        std::cout << "


    PI:      " << PI * 4;
        std::cout << "
    Cicles:  " << piCicles << std::endl;
        if ( maxCicles == piCicles ) //  Just in case someone has the time to run this program for long enough
            std::cout << "
    Max precise 64-bit computable value reached!" << std::endl;

        if ( !shouldStop ) getShouldStop();

        piLock.unlock();

        if ( shouldStop ) return;

        std::this_thread::sleep_for( std::chrono::milliseconds( updateDelay_ms ) );
    }
}

Cred că acest cod poate fi rulat și pe alte sisteme de operare, eliminând windows.h și antetul getShouldStop() funcție de asemenea, așa că și utilizatorii de Linux îl pot încerca.

Iată un link pentru a descărca executabilul pentru Windows, dacă doriți să îl testați.

Comentarii

  • Doar o notă secundară, acesta este un canal YouTube excelent. Daniel este un profesor excelent, în ciuda faptului că este foarte împrăștiat în creier. Îl recomand cu căldură.-  > Por Carcigenicate.
  • @Carcigenicate Într-adevăr, așa este 🙂 Doar pentru ca acest comentariu să nu fie neinformativ, am descoperit că seria Leibniz poate fi accelerată folosind niște tehnici cunoscute sub numele de „Accelerarea seriilor”, , poate într-o postare viitoare voi îmbunătăți acest program cu recenzia pe care am putut-o obține de aici și cu seria Leibniz accelerată, dacă reușesc să o rezolv, și poate făcând și o interfață de utilizare frumoasă :)-.  > Por Nikko77.
  • Obții o accelerare mai bună folosind formulele Euler-Machin sau formule asemănătoare lui Machin. Acestea se obțin tot din geometrie relativ simplă în planul complex și din aceeași serie arcus tangentă care dă formula Leibniz-Georgy.-  > Por Lutz Lehmann.
1 răspunsuri
Jerry Coffin

Structura de filetare

Îmi displace mai degrabă structura generală a codului. În loc de a avea două fire de execuție care se dispută o singură locație în care estimarea curentă a $pi$ este stocată, aș prefera ca firul calculatorului să calculeze aproximări succesive și să le scrie într-o coadă. Firul de afișare ar aștepta apoi ca valoarea să apară în coadă, să o afișeze și să repete. Blocajul și altele ar trebui să fie în coada de așteptare.

De asemenea, ați creat nu doar unul, ci două fire secundare. Acest lucru nu are niciun rost. Ați început cu firul principal al procesului, apoi tot ceea ce face este să creeze alte două fire, apoi să aștepte ca acestea să se termine. Ar putea la fel de bine să facă ceva util și să creeze doar un singur fir.

Formula

Mi se pare că, dacă vrei să încerci să-l faci mai rapid, ai putea măcar să faci o mică îmbunătățire a formulei. Una care este destul de rapidă și încă absolut trivială (de fapt, probabil mai simplă decât formula simplă a lui Leibniz) este:

$ {pi pe 4} = {sum_{n=0}^{infty} {2 pe (4n+1)(4n+3)}}}}$

Suma

Această formulă are o particularitate interesantă: dacă folosiți doar o adunare naivă peste această formulă, veți constata că se „blochează” – cu o dublă tipică, indiferent de câte iterații faceți, nu va obține decât aproximativ 8 cifre corecte.

Problema este că aveți un număr de mărime relativ mare și încercați să adăugați la acesta numere de mărime relativ mică. Fiecare dintre ele, în mod individual, este prea mic pentru a afecta numărul mai mare, așa că acesta rămâne constant.

În acest caz, există un remediu destul de simplu: începeți de la numărul cel mai mic numere mai mici, și mergeți spre numere mai mari. În acest fel, suma pe care o aveți și suma pe care o adăugați de fiecare dată au mai aproape aceeași mărime, astfel încât mai mulți termeni continuă să îmbunătățească rezultatul general.

Pentru acest program special, acest lucru prezintă totuși o mică problemă: întrucât încercăm să afișăm aproximația pe măsură ce se îmbunătățește în timp, dorim ca aproximațiile să se îmbunătățească în timp – și dacă începem de la cel mai mic și adunăm spre cel mai mare, aproximația noastră este cu adevărat groaznică pentru aproape tot timpul, apoi, la sfârșit, ultimele câteva iterații o fac brusc să devină un mult mai bine în grabă.

Ca un compromis, putem alege o cale de mijloc: să lucrăm de la început până la sfârșit, dar să lucrăm în „bucăți” de (să zicem) 10’000’000 de termeni, astfel încât să avem o estimare generală și o sumă temporară doar a celor mai recenți termeni. La începutul fiecărei iterații a buclei interioare, o luăm de la capăt de la o valoare de 0,0, astfel încât să nu avem un termen drastic mai mare care să domine atunci când efectuăm adunarea. Apoi, după ce am adunat acești termeni, adăugăm rezultatul la estimarea generală.

Acest lucru funcționează bine și cu actualizarea afișajului – de fiecare dată când adăugăm o valoare intermediară la estimarea globală, putem trimite rezultatul pentru a fi afișat.

Portabilitate

Am omis citirea asincronă a tastaturii, astfel încât codul să poată fi portabil. Aș prefera să folosesc o formulă și mai bună decât să las o formulă care converge extrem de lent să ruleze săptămâni în șir (și apoi să adaug cod neportabil pentru a le permite să renunțe mai ușor atunci când se plictisesc).

Cod

Așadar, procedând în acest fel, am putea ajunge la un cod în această ordine generală:

#include <iostream>
#include <iomanip>
#include <deque>
#include <mutex>
#include <thread>
#include <condition_variable>

namespace concurrent {
    template<class T>
    class queue {
        std::deque<T> storage;
        std::mutex m;
        std::condition_variable cv;
    public:
        void push(T const &t) {
            std::unique_lock<std::mutex> L(m);
            storage.push_back(t);
            cv.notify_one();
        }

        // This is not exception safe--if copying T may throw,
        // this can/will lose data.
        T pop() {
            std::unique_lock<std::mutex> L(m);
            cv.wait(L, [&] { return !storage.empty(); });
            auto t = storage.front();
            storage.pop_front();
            return t;
        }
    };
}

int main() { 
    concurrent::queue<double> values;

    auto t = std::thread([&] {
        double pi4 = 0.0;
        for (double n = 0.0; n < 8'000'000'000.0; n += 10'000'000) {
            double temp = 0.0;
            for (double i = 0; i < 10'000'000; i++)
                temp += 2.0 / ((4.0 * (n + i) + 1.0) * (4.0 * (n + i) + 3.0));
            pi4 += temp;
            values.push(4.0 * pi4);
        }
        values.push(-1.0);
    });

    double pi;
    while ((pi=values.pop())> 0.0)    
        std::cout << "r" << std::setw(11) << std::setprecision(11) << std::fixed << pi;
    t.join();
}

Pe mașina mea, acesta calculează Pi până la 10 locuri în aproximativ 45 de secunde. Cu excepția cazului în care computerul tău este destul de vechi/lent (cum este al meu), probabil că va rula mai repede decât atât pentru tine.

Pe de altă parte, pentru a urmări convergența valorii, acest lucru are un fel de problemă opusă: găsește în jur de șase sau șapte cifre aproape instantaneu, apoi macină mult timp până la alte câteva cifre. Din punct de vedere vizual, unii ar putea găsi formula Leibniz mai atrăgătoare pentru faptul că merge mai întâi deasupra, apoi sub, apoi din nou deasupra, din nou sub și așa mai departe pe măsură ce se apropie de adevărata valoare a lui Pi (deși acest lucru este mai evident dacă desenezi un grafic în loc să tipărești valorile).

Comentarii

  • Ideea de a lua mai multe sume parțiale este legată de însumarea Kahan (a se vedea Wikipedia), care poate fi utilă în acest caz. În acest caz, ar putea fi util să aveți mai mult de un termen de corecție. Pe de altă parte, de îndată ce introduceți o singură valoare inexactă în sumă, toate pariurile sunt anulate pentru orice cifră din dreapta.  > Por Roland Illig.
  • Vă mulțumesc foarte mult 🙂 Nu sunt foarte familiarizat cu C++, așa că voi încerca să învăț ceva din codul pe care l-ai scris, de asemenea, mulțumesc pentru informații @RolandIllig . Mâine voi face niște teste de reverse engineering(e târziu acum) la codul tău și voi accepta răspunsul tău dacă acesta satisface nevoile programului meu :)-.  > Por Nikko77.
  • @RolandIllig: Este adevărat că ar putea folosiți însumarea Kahan în acest caz, dar este de asemenea adevărat că însumarea Kahan impune o cantitate destul de mare de supraîncărcare. Probabil că ar trebui să fac un test rapid pentru a vedea câtă diferență reală (dacă există) face, dar nu am făcut-o încă.-.  > Por Jerry Coffin.
  • Ai vreun motiv pentru care folosești un dublu ca numărător de buclă în loc să faci o turnare după aceea?  > Por JVApen.
  • Cred că metoda ta pop este thread safe, totuși, dacă copia aruncă, cel mai probabil te afli într-o buclă infinită dacă ai încerca din nou-.  > Por JVApen.