6
\$\begingroup\$

Yesterday I saw a CodeTrain video in where Daniel Shiffman tried to approximate Pi using the Leibniz series. It was interesting so I decided to try it too.

I wrote this console application with a simple goal: Approximate as fast as possible as many digits of Pi the program can until it reaches the maximum number the cycle counter can achieve(without starting to commit calculation errors with the series formula) OR the program is stopped.

It starts by first running two background threads which execute two functions that do all the job: piLoop(), for calculating the series; and displayLoop(), for refreshing the screen and displaying the actual calculated Pi. To prevent concurrency problems, mutexes are used to prevent the threads from accessing Pi and other program shared variables at the same time. The program runs until piCycles reaches the maximum value a long long int can get OR the End key is pressed.

The final program runs, by glancing at it (that means a very unreliable measurement), at an approximated rate of 30 million cycles per second in my computer. That means it can calculate the 10th digit of PI (around 87 billion cycles) in about an hour.

I want to know if there are some optimizations that can be made so the program can run at all his power.

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 << "\n\n  ----- PI Calculator ----- \n\n  - \"Seeing how PI is calculated is more enjoyable while eating some snacks.\"";
        std::cout << "\n\n\n    PI:      " << PI * 4;
        std::cout << "\n    Cicles:  " << piCicles << std::endl;
        if ( maxCicles == piCicles ) //  Just in case someone has the time to run this program for long enough
            std::cout << "\n    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 ) );
    }
}

I think this code can also be run on other operating systems by removing the windows.h header and the getShouldStop() function as well, so Linux users may also give it a try.

Here is a link to download the Windows executable if you want to test it.

\$\endgroup\$
  • \$\begingroup\$ Just a side note, that's a excellent YouTube channel. Daniel is a great teacher despite being very scatter brained. I highly recommend it. \$\endgroup\$ – Carcigenicate Mar 20 at 0:17
  • \$\begingroup\$ @Carcigenicate Indeed it is :) Just for the sake of this comment not being uninformative, I discovered that the Leibniz series can be accelerated using some techniques known as "Series acceleration", maybe in a future post I will improve this program with the review I could get from here and the accelerated Leibniz series if I manage to work it out, and maybe making a nice UI too :) \$\endgroup\$ – Nikko77 Mar 20 at 4:18
  • \$\begingroup\$ You get a better acceleration using the Euler-Machin or Machin-like formulas. These are still obtained from relative simple geometry in the complex plane and the same arcus tangent series that gives the Leibniz-Georgy formula. \$\endgroup\$ – LutzL Apr 1 at 9:05
4
\$\begingroup\$

Threading structure

I rather dislike the overall structure of the code. Rather than having two threads contending over a single location where the current estimate of \$\pi\$ is stored, I'd rather have the calculator thread compute successive approximations, and write them to a queue. The display thread would then wait on value to show up in the queue, display it, and repeat. The locking and such should be in the queue itself.

You've also created not just one, but two secondary threads. There's no real point in this. You started with the process' primary thread, then all it does is create two other threads, then wait for them to finish. It might as well do something useful, and only create one other thread.

Formula

It seems to me that if you're going to try to make it faster, you could at least do a minor improvement to the formula. One that's quite a bit faster, and still utterly trivial (in fact, arguably simpler than the plain Leibniz formula) is:

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

Summation

This does have one interesting twist: if you just use a naive summation over this formula, you'll find that it "gets stuck"--with a typical double, no matter how many iterations you do, it will only ever get about 8 digits correct.

The problem is that you have a number of relatively large magnitude, and you're trying to add numbers of relatively small magnitude to it. Each of them individually is too small to affect the larger number, so it remains constant.

In this case, there's a pretty simple cure: start from the smallest numbers, and work your way toward larger numbers. That way, the sum you have and the amount you're adding to it each time have closer to the same magnitude, so more terms continue to improve the overall result.

For this particular program, that has a bit of a problem though: since we're trying to display the approximation as it gets better over time, we want approximations that actually do get better over time--and if we start from the smallest, and sum toward the largest, our approximation is truly terrible for almost the entire time, then at the very end, the last few iterations suddenly make it a lot better in a hurry.

As a compromise, we can take a middle ground: work from beginning to end, but work in "chunks" of (say) 10'000'000 terms, so we have an overall estimate, and we have a temporary sum of only the most recent terms. At the beginning of each iteration of the inner loop, we start over from a value of 0.0, so we don't have a drastically larger term dominating when we do the addition. Then when we've added those terms together, we add the result to the overall estimate.

This also works nicely with updating the display--each time we add an intermediate value to our overall estimate, we can send the result out to be displayed.

Portability

I've left out the asynchronous keyboard reading, so the code can be portable. I'd rather use an even better formula than just let one that converges extremely slowly run for weeks on end (and then add non-portable code to let them quit more easily when they get bored).

Code

So, doing things this way, we could end up with code on this general order:

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

On my machine, this calculates Pi to 10 places in about 45 seconds. Unless your computer is quite old/slow (like mine is) it'll probably run faster than that for you.

On the other hand, for watching the value converge, this does have kind of the opposite problem: it finds around six or seven digits nearly instantly, then grinds for a long time to a few more. Visually, some might find the Leibniz formula more appealing for the fact that it goes first above, then below, then back above, below again, and so on as it approaches the true value of Pi (though that's more apparent if you draw a graph rather than just printing out the values).

\$\endgroup\$
  • 1
    \$\begingroup\$ The idea of taking multiple partial sums is related to Kahan summation (see Wikipedia), which may be useful in this case. In this case it might be worth having more than one correction term. On the other hand, as soon as you feed a single inexact value to the sum, all bets are off for any digits to the right. \$\endgroup\$ – Roland Illig Mar 21 at 6:15
  • \$\begingroup\$ Thank you very much :) I'm not very familiar with C++ so I'll try to learn something from the code you wrote, also thanks for the info @RolandIllig . I'll make some reverse engineering tests tomorrow(its late now) to your code and accept your answer if it satisfies my program needs :) \$\endgroup\$ – Nikko77 Mar 21 at 6:32
  • 1
    \$\begingroup\$ @RolandIllig: It's true that you could use Kahan Summation in this case, but it's also true that Kahan summation imposes a fair amount of overhead. I should probably do a quick test to see how much (if any) real difference it makes, but I haven't yet. \$\endgroup\$ – Jerry Coffin Mar 21 at 6:35
  • \$\begingroup\$ Any reason you use a double as loop counter instead of casting afterwards? \$\endgroup\$ – JVApen Mar 21 at 7:22
  • \$\begingroup\$ I think your pop method is thread safe, however, if copying throws, you most likely are in an infinite loop if you would try again \$\endgroup\$ – JVApen Mar 21 at 7:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.