// L-11 MCS 572 Fri 3 Feb 2012 : powers_tbb.cpp
// Simple C++ program to raise every element of an array
// of n complex doubles to the power d,
// using the parallel_for of the Intel TBB library.
// For timing purposes n, d, and the verbose level
// can be provided at the command line.

#include <cstdlib>
#include <cmath>
#include <complex>
#include <ctime>
#include <iostream>
#include <iomanip>
#include "tbb/tbb.h"
#include "tbb/blocked_range.h"
#include "tbb/parallel_for.h"
#include "tbb/task_scheduler_init.h"

using namespace std;
using namespace tbb;

typedef complex<double> dcmplx;

dcmplx random_dcmplx ( void );
// generates a random complex number
// on the complex unit circle

void write_numbers ( int n, dcmplx *x );
// writes the array of n doubles in x

void compute_powers ( int n, dcmplx *x,
                      dcmplx *y, int d );
// for arrays x and y of length n,
// on return y[i] equals x[i]**d

class ComputePowers
{
   dcmplx *const c; // numbers on input
   int d;           // degree
   dcmplx *result;  // output
   public:
      ComputePowers(dcmplx x[], int deg, dcmplx y[])
         : c(x), d(deg), result(y) { }
      void operator()
         ( const blocked_range<size_t>& r ) const
      {
         for(size_t i=r.begin(); i!=r.end(); ++i)
         {
            dcmplx z(1.0,0.0);
            for(int j=0; j < d; j++) z = z*c[i];
            result[i] = z;
         }
      }
};

int main ( int argc, char *argv[] )
{
   int v = 1;    // verbose if > 0
   if(argc > 3) v = atoi(argv[3]);
   int dim;      // get the dimension
   if(argc > 1)
      dim = atoi(argv[1]);
   else
   {
      cout << "how many numbers ? ";
      cin >> dim;
   }
   // fix seed to compare with serial
   srand(20120203); // srand(time(0));
   dcmplx r[dim];
   for(int i=0; i<dim; i++)
      r[i] = random_dcmplx();
   if(v > 0) write_numbers(dim,r);

   int deg;      // get the degree
   if(argc > 1)
      deg = atoi(argv[2]);
   else
   {
      cout << "give the power : ";
      cin >> deg;
   }
   dcmplx s[dim];
   task_scheduler_init init(task_scheduler_init::automatic);
   parallel_for(blocked_range<size_t>(0,dim),
                ComputePowers(r,deg,s));
   if(v > 0) write_numbers(dim,s);
   return 0;
}

dcmplx random_dcmplx ( void )
{
   int r = rand();
   double d = ((double) r)/RAND_MAX;
   double e = 2*M_PI*d;
   dcmplx c(cos(e),sin(e));
   return c;
}

void write_numbers ( int n, dcmplx *x )
{
   for(int i=0; i<n; i++)
      cout << scientific << setprecision(4)
           << "x[" << i << "] = ( " << x[i].real()
           << " , " << x[i].imag() << ")\n";
}

void compute_powers ( int n, dcmplx *x,
                      dcmplx *y, int d )
{
   for(int i=0; i < n; i++) // y[i] = pow(x[i],d);
   {                        // pow is too efficient
      dcmplx r(1.0,0.0);
      for(int j=0; j < d; j++) r = r*x[i];
      y[i] = r;
   }
}