#include <stdio.h>
#include <time.h>
#include <cassert>
#include <cstdio>
#include <cstring>
#include <math.h>
#include <cstdlib>
#include <cctype>
#include "tbb/parallel_reduce.h"
#include "tbb/task_scheduler_init.h"
#include "tbb/blocked_range.h"

using namespace std;
using namespace tbb;

#define Nthreads 2

#define GrainSize      10000

long long num_steps =  1000000000;

class MyPi {
                double *const my_step;
public:
                double sum;
                void operator()( const blocked_range<size_t>& r ) {
                                double step = *my_step;
                                double x;
                    for (size_t i=r.begin(); i!=r.end(); ++i)
                                {
                                                x = (i + .5)*step;
                                                sum = sum + 4.0/(1.+ x*x);
                                }
                }

                MyPi( MyPi& x, split ) : my_step(x.my_step), sum(0) {}
                void join( const MyPi& y ) {sum += y.sum;}
                MyPi(double *const step) : my_step(step), sum(0) {}
};

int main(int argc, char* argv[])
{
                clock_t start, stop;
                double pi;
                double width = 1./(double)num_steps;
                printf("Approximating pi with %ld rectangles\n\n", num_steps);

                MyPi step((double *const)&width);
                task_scheduler_init init(task_scheduler_init::deferred);

                start = clock();
                init.initialize(Nthreads);
                parallel_reduce(blocked_range<size_t>(0,num_steps,GrainSize), step);
                pi = step.sum*width;
                stop = clock();

                printf("The value of PI is %15.12f\n",pi);

                printf("The time to calculate PI was %f seconds\n",((double)(stop - start)/1000.0));

                return 0;

}