π in C++

Recently I’ve helped out with a couple of queries about calculating PI.  They’re programming assingments for students rather than real mathematical problems. I thought I’d publish them here.

The Monte Carlo Method

The idea is to use a game of chance to esitmate a value for PI.

You start by drawing a square of length 1 with a quarter of a circle of radius 1 inside the square.

If you were to throw a dart randomly at this shape, what are the chances of it landing inside the quarter circle?

The answer is (area of quarter circle) / (area of square). As the area of a square of length 1 is 1, the chance becomes (area of quarter circle) which is (PI * radius squared / 4). And as the radius of 1 squared is still 1, the answer is PI/4.

So the chance of your dart landing in the quarter circle is PI/4.

If we do this once, we’ll get a value of 0 or 4 for PI. But as we do it more and more, we get better approximations for PI.

You could ask 3 friends to help you do this test be throwing darts and checking how many were in the quarter circle for a fixed number of throws. At the end, you’d find the ratio of how many were in, then multiply that by 4 and get a value for PI.

Back to the computer program. We throw the dart by getting random (x,y) coordinates. We know it’s in the circle if the distance from the origin (0,0) is less than or equal to 1. And we calculate the distance using Pythagoras’ theorem, c2 = a2 + b2. But as 12 is still 1, we can just work with c2 and not attempt the square root.

Each thread does the same test and you just check the results at the end.

#include <pthread.h>
#include <iostream>
#include <vector>
#include <stdlib.h>
#include <stdint.h>
namespace MonteCarloPI
{
    struct Info
    {
        typedef uint32_t value_type;

        Info() : inPoints(0), totalPoints(0) {}

        value_type inPoints;
        value_type totalPoints;
    };

    void* func(void* param)
    {
        if (Info* info = reinterpret_cast<Info*>(param))
        {
            // do some work
            for (size_t i = 0; i != info->totalPoints; ++i)
            {
                std::pair<double, double> pt(rand()/(double)RAND_MAX, rand()/(double)RAND_MAX);
                double d = pt.first*pt.first + pt.second*pt.second;
                if (d <= 1.0) 
                    ++(info->inPoints); 
            } 
        }

        return param;
    }

    double pi(size_t nthreads, size_t niters)
    {
        std::vector<pthread_t> handles(nthreads);
        std::vector<Info> info(nthreads);

        // start the threads
        for (size_t i = 0; i != nthreads; ++i)
        {
            info[i].totalPoints = niters;
            pthread_create(&handles[i], NULL, func, &info[i]);
        }

        // wait for the threads to complete
        for (size_t i = 0; i != nthreads; ++i)
            pthread_join(handles[i], NULL);

        // process the reply
        double inPoints = 0.0;
        double totalPoints = 0.0;
        for (size_t i = 0; i != nthreads; ++i)
        {
            inPoints += info[i].inPoints;
            totalPoints += info[i].totalPoints;
        }
        const double pi = (4.0 * inPoints / totalPoints);

        return pi;
    }
}

int main()
{
    srand(2012);

    const size_t NTHREADS = 4;
    const size_t NITERS = 100*1000*1000;

    double pi = MonteCarloPI::pi(NTHREADS, NITERS);
    std::cout << "pi = " << pi << std::endl;

    return 0;
}

Leibniz Method

It’s probably best to start with a description of the algorithm. http://en.wikipedia.org/wiki/Leibniz_formula_for_%CF%80

#include <iostream>
#include <stdlib.h>
#include <stdint.h>
namespace LiebnizPI
{
    typedef uint32_t value_type;

    double pi(value_type niters)
    {
        double pi = 0.0;
        for (size_t i = 0; i != niters; ++i)
        {
            if (i % 2)
                pi -= 1.0 / (2*i + 1);
            else
                pi += 1.0 / (2*i + 1);
        }

        return 4.0 * pi;
    }
}

int main()
{
    srand(2012);

    const LiebnizPI::value_type NITERS = 1*1000*1000;

    double pi = LiebnizPI::pi(NITERS);
    std::cout << "pi = " << pi << std::endl;

    return 0;
}
Advertisements

~ by kbw333 on July 19, 2012.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

 
%d bloggers like this: