#include <omp.h>
#include <stdio.h>
#include <stdlib.h>
#include <limits.h>
#include <getopt.h>
#include <string.h>
#include <stdint.h>

#include "cycletimer.h"

static void usage(char *name) {
    printf("Usage: %s [-h] [-t THD] [-n NVAL]\n", name);
    printf("    -h       Print this message\n");
    printf("    -t THD   Set number of threads\n");
    printf("    -N NVAL  Set number of values\n");
    printf("    -n LOGNVAL  Set log of number of values\n");
    exit(0);
}


static inline double dmin(double x, double y) {
    return x < y ? x : y;
}

static inline double dmax(double x, double y) {
    return x > y ? x : y;
}

/* Random number generation */
typedef uint32_t random_t;

/* Standard parameters */
#define GROUPSIZE 2147483647
#define MVAL  48271
#define VVAL  16807
#define INITSEED  418

static inline random_t rnext(random_t *seedp, random_t x) {
    uint64_t s = (uint64_t) *seedp;
    random_t val = ((x+1) * VVAL + s * MVAL) % GROUPSIZE;
    *seedp = (random_t) val;
    return val;
}

/* Reinitialize seed based on list of seeds, where list has length len */
void reseed(random_t *seedp, random_t seed_list[], size_t len) {
    *seedp = INITSEED;
    size_t i;
    for (i = 0; i < len; i++)
	rnext(seedp, seed_list[i]);
}

/* Generate double in range [0.0, upperlimit) */
double next_random_float(random_t *seedp, double upperlimit) {
    random_t val = rnext(seedp, 0);
    return ((double) val / (double) GROUPSIZE) * upperlimit;
}


typedef struct {
    unsigned long count;
    double min;
    double max;
    double sum;
} stat_t;

/*
  Partition computation into blocks of nearly equal size
  Require 0 <= id < nblock
*/
void block_partition(unsigned long total_count, unsigned long nblock, unsigned long id, unsigned long *local_count, unsigned long *local_start) {
    unsigned long lcount = total_count / nblock;
    unsigned long xcount = total_count % nblock;
    unsigned long lstart = id * lcount;
    if (id < xcount) {
	lcount += 1;
        lstart += id;
    } else {
        lstart += xcount;
    }
    *local_count = lcount;
    *local_start = lstart;
}


void gen_stat(unsigned long start, unsigned long count, stat_t *stat) {
    double min = 1e50;
    double max = -1e50;
    double sum = 0.0;
    unsigned long i;
    random_t seed = 0;
    for (i = 0; i < count; i++) {
	random_t s[1] = { i + start };
	reseed(&seed, s, 1);
	double val = next_random_float(&seed, 1.0);
	sum += val;
	min = dmin(min, val);
	max = dmax(max, val);
    }
    stat->count = count;
    stat->min = min;
    stat->max = max;
    stat->sum = sum;
}

void run(unsigned long count, int nthread) {
    stat_t s[nthread];
    #pragma omp parallel
    {
	unsigned long local_count, local_start;
	unsigned long id = omp_get_thread_num();
	block_partition(count, nthread, id, &local_count, &local_start);
	gen_stat(local_start, local_count, &s[id]);
    }
    int t;
    double global_min = s[0].min;
    double global_max = s[0].max;
    double global_sum = s[0].sum;

    for (t = 1; t < nthread; t++) {
	global_min = dmin(global_min, s[t].min);
	global_max = dmax(global_max, s[t].max);
	global_sum += s[t].sum;
    }
    double global_average = global_sum / count;
    printf("Global result %d threads over %lu values.  Avg = %f, Max = %f, Min = %f\n",
	   nthread, count, global_average, global_max, global_min);

}

int main(int argc, char **argv) {
    int c;
    unsigned long n = 1L << 35;
    int nthread = 0;
    while ((c = getopt(argc, argv, "ht:N:n:")) != -1) {
	switch (c) {
	case 'h':
	    usage(argv[0]);
	    break;
	case 't':
	    nthread = atoi(optarg);
	    break;
	case 'N':
	    n = strtoul(optarg, NULL, 10);
	    break;
	case 'n':
	    n = 1LU << atoi(optarg);
	    break;
	default:
	    printf("Unknown flag '%c'\n", c);
	    usage(argv[0]);
	}
    }
    double startTime = currentSeconds();
#ifdef RUN_MIC
    #pragma offload target(mic)
#endif
    {
	if (nthread < 1)
	    nthread = omp_get_max_threads();
	omp_set_num_threads(nthread);
	run(n, nthread);
    }
    double delta = currentSeconds() - startTime;
    double t_op = 1e9 * delta / n;
    double tpt_op = t_op * nthread;
    printf("Ran with n = %lu in %.1f ms. nthread = %d Total time per op = %.3f ns.  Time per op per thread = %.3f ns\n",
	   n, 1000.0 * delta, nthread, t_op, tpt_op);
    return 0;
}