/*
 * C program to compute the value of pi using Monte Carlo method.
 *
 * command-line arguments:  number_of_samples, seed
 *
 * Preliminary MPI parallel version.  Does not work well -- results
 * vary depending on the number of processes, and are worse for more
 * processes (for reasons as discussed in class).
 */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>   /* has fabs() */
/* copied from not-strictly-standard part of math.h */
#define M_PI		3.14159265358979323846
#include <mpi.h>

#include "timer.h"  /* has get_time() */
#include "utility.h"    /* has error_exit() */

/* main program */

int main(int argc, char* argv[]) {

    char* usage_fmt = "usage:  %s number_of_samples seed\n";
    char* end_ptr_for_strtol;

    int myID;
    int nprocs;

    /* initialize for MPI */
    if (MPI_Init(&argc, &argv) != MPI_SUCCESS) {
        fprintf(stderr, "MPI initialization error\n");
        error_exit(NULL);
    }
    MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
    MPI_Comm_rank(MPI_COMM_WORLD, &myID);

    /* process command-line arguments */
    if (argc != 3)  {
        fprintf(stderr, usage_fmt, argv[0]);
        error_exit(NULL);
    }
    long num_samples = strtol(argv[1], &end_ptr_for_strtol, 10);
    if (*end_ptr_for_strtol != '\0') {
        fprintf(stderr, usage_fmt, argv[0]);
        error_exit(NULL);
    }
    long seed = strtol(argv[2], &end_ptr_for_strtol, 10);
    if (*end_ptr_for_strtol != '\0') {
        fprintf(stderr, usage_fmt, argv[0]);
        error_exit(NULL);
    }

    /* record start time */
    double start_time = get_time();

    /* do calculation */ 
    srand((unsigned int) seed);
    long count;
    long local_count = 0;
    double x, y;
    for (long i = myID; i < num_samples; i += nprocs) {
        x = (double) rand() / (double) (RAND_MAX);
        y = (double) rand() / (double) (RAND_MAX);
        if ((x*x + y*y) <= 1.0)
            ++local_count;
    }
    MPI_Reduce(&local_count, &count, 1, MPI_LONG, MPI_SUM, 0, MPI_COMM_WORLD);
    if (myID == 0)
        pi = 4.0 * (double) count / (double) num_samples;

    /* record end time and print results */
    end_time = get_time();
    if (myID == 0) {
        printf("MPI program results with %d processes:\n", nprocs);
        printf("number of samples = %ld, seed = %ld\n", num_samples, seed);
        printf("estimated pi = %12.10f\n", pi);
        printf("difference between estimated pi and math.h M_PI = %12.10f\n", 
                fabs(pi - M_PI));
        printf("time to compute = %g seconds\n", end_time - start_time);
    }

    /* clean up and return */
    MPI_Finalize();
    return EXIT_SUCCESS;
}