#include <iostream>
#include <cstdlib>
#include <cmath>
#include <xmmintrin.h>
using namespace std;
#include <getopt.h>

struct option long_options[] {
	{"method", required_argument, 0, 'm'},
	{"size", required_argument, 0, 's'},
	{"iterations", required_argument, 0, 'i'},
	{0,0,0,0}
};

int method = 1;
int vsize = 10*1024*1024;
int iterations = 1;
	
	
	
typedef float real;

real prod1(real * x, real * y, real * z, int size) {
	real total = 0.0;
	for (int i=0; i<size; ++i) {
		z[i] = (z[i] + x[i]) / y[i];
		total += z[i];
//		cout << "z[" << i << "]=" << z[i] << " " << total << endl;
	}
	
	return total;
}


real prod2(real * x, real * y, real * z, int size) {
	real total = 0.0;

	#pragma omp parallel for reduction(+:total) 
	for (int i=0; i<size; ++i) {
		z[i] = (z[i] + x[i]) / y[i];
		total += z[i];
	}
	return total;
}

real prod3(real * x, real * y, real * z, int size) {
	real total = 0.0;

	#pragma omp simd reduction(+:total) 
	for (int i=0; i<size; ++i) {
		z[i] = (z[i] + x[i]) / y[i];
		total += z[i];
	}
	return total;
}

#pragma omp declare simd aligned(x,y,z:16) simdlen(4)
real prod4(real * x, real * y, real * z, int size) {
	real total = 0.0;

	#pragma omp parallel for simd reduction(+:total)
	for (int i=0; i<size; ++i) {
		z[i] = (z[i] + x[i]) / y[i];
		total += z[i];
	}
	
	return total;
}

typedef real (*ProductType)(real *x, real *y, real *z, int size);

extern "C" {
real prod_asm(real * x, real * y, real * z, int size);
real prod_asm_unroll2(real * x, real * y, real * z, int size);
real prod_asm_unroll4(real * x, real * y, real * z, int size);
}

ProductType products[] = {
	nullptr,
	prod1,
	prod2,
	prod3,
	prod4,
	prod_asm,
	prod_asm_unroll2,
	prod_asm_unroll4
};

int main(int argc, char *argv[]) {
	int cc, index;
	
	while ((cc = getopt_long(argc, argv, "m:s:i:", long_options, &index)) != -1) {
		switch (cc) {
			case 'm':
				method = atoi(optarg);
				break;
			case 's':
				vsize = atoi(optarg);
				break;
			case 'i':
				iterations = atoi(optarg);
				break;	
			default: /* '?' */
				exit(EXIT_FAILURE);
		}
	}
	real *a, *b, *c;
	

	a = (real *) _mm_malloc(vsize * sizeof(real), 16);
	b = (real *) _mm_malloc(vsize * sizeof(real), 16);
	c = (real *) _mm_malloc(vsize * sizeof(real), 16);
	
	for (int i=0; i<vsize; ++i) {
		//a[i] = 1 + (i % 10);
		a[i] = 0.2;
	}
	for (int i=0; i<vsize; ++i) {
		//b[i] = 1 + (i % 133);
		b[i] = 5;
	}
	for (int i=0; i<vsize; ++i) {
		c[i] = 0;
	}
	
	real result = 0;

	cout << "method=" << method << endl;	
	for (int k=0; k<iterations; ++k) {	
		cout << ".";
		result = (*products[method])(a,b,c,vsize);
	}
	cout << endl;
	cout << "result=" << std::fixed << result << endl;
	
	_mm_free(a);
	_mm_free(b);
	_mm_free(c);
	
	return EXIT_SUCCESS;
}