#include <stdio.h>
#include <stdlib.h>
#include <CL/cl.h>
#include <CL_2_0\cl.h>
#include <c:/altera/15.1/aclrte-windows64/host/include/CL/opencl.h>
#include <iostream>
#include <algorithm>    // Needed for the "max" function
#include <cmath>
#include "c:/altera/15.1/aclrte-windows64/board/de1soc/examples/common/inc/aocl_utils.h"



//Kernel
const char* programSource =
"__kernel								\n"
"void vecadd(__int num_elements,                        	\n"
"				__global double *A,			\n"
"              __global double *B,					\n"
"              __global double *C)					\n"
"{									\n"
"									\n"
"int idx = get_global_id(0);					\n"
"double euclidSq = 0.0;						\n"
"double x_select = 0.0;                                 	\n"
"for (int i = 0; i<num_elements; ++i){                  	\n"
"do {									\n"
"      euclidSq = A[i]*A[i] + B[i]*B[i];	            		\n"
"      x_select = A[i];                                 	\n"
"   } while (euclidSq >= 1.0);					\n"
"   C[i] = x_select*sqrt(-2 * log(euclidSq) / euclidSq);	\n"
"}                                                      	\n"
"}                                                      	\n"
;
using namespace aocl_utils;

int main() {

	//**********Problem Data***********************
	const int elements = 100000;
	double S = 100.0;  // Option price
	double K = 100.0;  // Strike price
	double r = 0.05;   // Risk-free rate (5%)
	double v = 0.2;    // Volatility of the underlying (20%)
	double T = 1.0;    // One year until expiry
	double call, put;


	size_t datasize = sizeof(double)*elements;

	double *A = (double*)malloc(datasize);
	double *B = (double*)malloc(datasize);
	double *C = (double*)malloc(datasize);

	//Loads the random x & y values into arrays A and B in host memory
	int i;
	for (i = 0; i < elements; i++) {
		A[i] = 2.0 * rand() / static_cast<double>(RAND_MAX) - 1;
		B[i] = 2.0 * rand() / static_cast<double>(RAND_MAX) - 1;
	}

	cl_int status;

	//Find the OpenCL Platform
	cl_platform_id platform;
	status = clGetPlatformIDs(1, &platform, NULL);


	//Get the Devices
	cl_device_id device;
	status = clGetDeviceIDs(platform, CL_DEVICE_TYPE_ALL, 1, &device, NULL);

	//Create the Context
cl_context context = clCreateContext(NULL, 1, &device, NULL, NULL, &status);

	//Create the command queue for the device chosen
cl_command_queue cmdQueue = clCreateCommandQueue(context, device, 0, &status);

	//Create the buffers for memory transfers between host and device
	cl_mem bufA = clCreateBuffer(context, CL_MEM_READ_ONLY, datasize,
		NULL, &status);
	cl_mem bufB = clCreateBuffer(context, CL_MEM_READ_ONLY, datasize,
		NULL, &status);
	cl_mem bufC = clCreateBuffer(context, CL_MEM_WRITE_ONLY, datasize,
		NULL, &status);

	//Write the data from host memory to the kernel buffers
	status = clEnqueueWriteBuffer(cmdQueue, bufA, CL_FALSE, 0,
		datasize, A, 0, NULL, NULL);
	status = clEnqueueWriteBuffer(cmdQueue, bufB, CL_FALSE, 0,
		datasize, B, 0, NULL, NULL);

	//Create the program by extracting the kernel
	cl_program program = clCreateProgramWithSource(context, 1,
		(const char**)&programSource, NULL, &status);

	//Build the program for the device
	status = clBuildProgram(program, 1, &device, NULL, NULL, NULL);

	//Create the kernel object
	cl_kernel kernel = clCreateKernel(program, "vecadd", &status);

	//Set the arguments for the kernel
	status = clSetKernelArg(kernel, 0, sizeof(cl_mem), &bufA);
	status = clSetKernelArg(kernel, 1, sizeof(cl_mem), &bufB);
	status = clSetKernelArg(kernel, 2, sizeof(cl_mem), &bufC);

	//Define an index space of wor-items for execution
	size_t indexSpaceSize[1], workGroupSize[1];

	//These are 'elements' wor-items
	indexSpaceSize[0] = elements;
	workGroupSize[0] = elements;

	//Execute the kernel
status = clEnqueueNDRangeKernel(cmdQueue, kernel, 1, NULL, indexSpaceSize,
		workGroupSize, 0, NULL, NULL);

	//Read the device output buffer to the host output array
	status = clEnqueueReadBuffer(cmdQueue, bufC, CL_TRUE, 0,
		datasize, C, 0, NULL, NULL);

	//******************call function**********************************
	double S_adjust = S * exp(T*(r - 0.5*v*v));
	double S_cur = 0.0;
	double payoff_sum = 0.0;

	for (int j = 0; j < elements; j++) {
		double gauss_bm = C[j]; //pulls the values returned by the kernel
		//printf("%d", gauss_bm);
		S_cur = S_adjust * exp(sqrt(v*v*T)*gauss_bm);
		payoff_sum += std::max(S_cur - K, 0.0);
	}

	call = (payoff_sum / static_cast<double>(elements)) * exp(-r*T);
	//****************************************************************

	//****************put function************************************
	S_adjust = S * exp(T*(r - 0.5*v*v));
	S_cur = 0.0;
	payoff_sum = 0.0;

	for (int j = 0; j < elements; j++) {
		double gauss_bm = C[j]; //pulls the values returned by the kernel
		//printf("%d", gauss_bm);
		S_cur = S_adjust * exp(sqrt(v*v*T)*gauss_bm);
		payoff_sum += std::max(K - S_cur, 0.0);
	}

	put = (payoff_sum / static_cast<double>(elements)) * exp(-r*T);
	//****************************************************************


	// Finally we output the parameters and prices**********
	std::cout << "Number of Paths: " << elements << std::endl;
	std::cout << "Underlying:      " << S << std::endl;
	std::cout << "Strike:          " << K << std::endl;
	std::cout << "Risk-Free Rate:  " << r << std::endl;
	std::cout << "Volatility:      " << v << std::endl;
	std::cout << "Maturity:        " << T << std::endl;

	std::cout << "Call Price:      " << call << std::endl;
	std::cout << "Put Price:       " << put << std::endl;
	//******************************************************


	//Free OpenCL resources
	clReleaseKernel(kernel);
	clReleaseProgram(program);
	clReleaseCommandQueue(cmdQueue);
	clReleaseMemObject(bufA);
	clReleaseMemObject(bufB);
	clReleaseMemObject(bufC);
	clReleaseContext(context);

	//Free host resources
	free(A);
	free(B);
	free(C);

	return 0;
}