// Copyright David Abrahams 2008. Distributed under the Boost
// Software License, Version 1.0. (See accompanying
// file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
#include "boost/thread.hpp"
#include "boost/iterator/counting_iterator.hpp"
#include "boost/lambda/lambda.hpp"
#include "boost/lambda/bind.hpp"
#include "boost/timer.hpp"
#include "boost/type_traits/remove_const.hpp"
#include "boost/type_traits/remove_reference.hpp"
#include "boost/array.hpp"
#include "boost/date_time/local_time/local_time.hpp"
#include <vector>
#include <iterator>
#include <numeric>
#include <algorithm>
#include <cmath>
#include <iostream>

namespace my
{
double sqrt(double x)
{
    return std::sqrt(x);
}
}

// sum terms [i-j) of the power series for pi/4
double calc(std::size_t i, std::size_t j)
{
    double sum = 0.0;
    
    while (i < j)
    {
        sum += (i % 2 == 0 ? 1.0 : -1.0) / (2*i + 1);
        ++i;
    }
    return sum;
}

std::size_t const max_parallelism = 20;
typedef boost::array<boost::thread, max_parallelism> thread_pool;
typedef boost::array<double, max_parallelism> result_pool;

double calc_pi(std::size_t nterms, std::size_t parallelism, thread_pool& threads, result_pool& results)
{
    assert(parallelism > 0 && parallelism <= threads.size());  // make sure we've allocated enough space
    assert(threads.size() == results.size());  // make sure we've allocated enough space

    using namespace boost::lambda;
    double const step = nterms * 1.0 /  parallelism;
    for (int i = 1; i < parallelism; ++i)
    {
        threads[i] = boost::thread(
            var(results[i]) = bind(calc, i * step, (i+1) * step)
        );
    }

    // do a calculation in this thread
    results[0] = calc(0, step);

    // wait for all threads to finish
    for (int i = 1; i < parallelism; ++i)
        threads[i].join();

    return std::accumulate( results.begin(), results.begin() + parallelism, 0.0 );
}



int main()
{
    std::size_t const n = 200000000;


    thread_pool threads;
    result_pool results;

    std::cout << "warming up the cache..." << std::endl;
    double pi = 4*calc_pi(n, 4, threads, results);
    std::cout << "start"  << std::endl;
    
    for (std::size_t p = 1; p < max_parallelism; ++p)
    {
        using namespace boost::gregorian;
        using namespace boost::posix_time;
        using namespace boost::local_time;
        
        ptime start = microsec_clock::local_time();
        boost::timer time;

        pi = 4*calc_pi(n, p, threads, results);
        microsec_clock::time_duration_type t0 = microsec_clock::local_time() - start;
        double t1 = time.elapsed();
        
        std::cout << "calculating pi = " << pi << " using " << p << " threads took: " << t0 << " but timer reports: " << t1 << std::endl;
    }
}