#include <cuda.h>
#include <curand_kernel.h>
#include <stdio.h>
#include <time.h>
#include<stdio.h>
#include<math.h>
#include<time.h>
#include<stdlib.h>

#define N 1
const int n = 3;

__global__ void setup_kernel(curandState * state, unsigned long seed)
{
    int id = threadIdx.x;
    curand_init(seed, id, 0, &state[id]);
}

__global__ void generate(curandState* globalState, float* randomArray)
{
    int ind = threadIdx.x;
    curandState localState = globalState[ind];
    float RANDOM = curand_uniform(&localState);
    randomArray[ind] = RANDOM;
    globalState[ind] = localState;
}

double f(double x[])
{
    return x[0]*x[1]*x[2];
}

double r()
{
    dim3 tpb(N, 1, 1);
    curandState* devStates;
    float* randomValues = new float[N];
    float* devRandomValues;

    cudaMalloc(&devStates, N * sizeof(curandState));
    cudaMalloc(&devRandomValues, N * sizeof(*randomValues));// setup seeds
    setup_kernel <<< 1, tpb >>> (devStates, time(NULL));

    //printf("%s\n", cudaGetErrorString(cudaGetLastError()));// generate random numbers
    generate <<< 1, tpb >> > (devStates, devRandomValues);

    //printf("%s\n", cudaGetErrorString(cudaGetLastError()));

    cudaMemcpy(randomValues, devRandomValues, N * sizeof(*randomValues), cudaMemcpyDeviceToHost);

    cudaFree(devRandomValues);
    cudaFree(devStates);

    return randomValues[0];
}

void integral(int number_of_points, double m[]) {
    //srand(time(NULL));

    int in = 0;
    double k[n + 1];
    for (int i = 0; i < number_of_points; i++)
    {
        for (int i = 0; i < n + 1; i++)
        {
            k[i] = r()*(m[2 * i + 1] - m[2 * i]) + m[2 * i];
        }

        if (f(k) < k[n] && k[n] < 0)
        {
            in--;
        }
        if (f(k) > k[n] && k[n] > 0) {
            in++;
        }
    }

    double v = 1;
    for (int i = 0; i < n + 1; i++)
    {
        v = v * (m[2 * i + 1] - m[2 * i]);
    }
    double real = 9176.7;
    double result = v*in/number_of_points;
    printf(" result= %lf, error %lf\n ", result, fabs(result - real));
}

int main(int argc, char** argv)
{
    //const int n = 3;
    //printf("enter n: \n");
    //scanf("%d", &n);

    double m[2 * (n + 1)];

    for (int i = 0; i < n; i++)
    {
        printf("enter min x[%d]: \n", i);
        scanf("%lf", &m[2 * i]);

        printf("enter max x[%d]: \n", i);
        scanf("%lf", &m[2 * i + 1]);
    }

    for (int i = 1; i < 2; i++)
    {
        printf("enter min y: \n");
        scanf("%lf", &m[2 * n]);

        printf("enter max y: \n");
        scanf("%lf", &m[2 * n + i]);
    }


    for (int number_of_points = 100; number_of_points < 1000000000; number_of_points *= 10) {
        printf("For %d points\n", number_of_points);
        integral(number_of_points, m);
    }

}