Hey here is sudo command for you, comment code that I will send you, explaining and documenting the functions of each step / line / function etc. Most importantly, don't change this code, only copy it and add there a comments. 

Code:

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

#define THREADS 512

__device__ double generate(curandState* globalState, int ind) {
    curandState localState = globalState[ind];
    double RANDOM = curand_uniform( &localState );
    globalState[ind] = localState;
    return RANDOM;
}

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

__global__ void approximate(double *numVector, curandState* globalState, int LOOPS){
    int id = threadIdx.x + blockIdx.x * blockDim.x;
    __shared__ int temp[THREADS];

    double _counter = 0;
    for(int i = 0; i < LOOPS; i++) {
      double x = generate(globalState, id), y = generate(globalState, id);
      double z = x * x + y * y;
      if(z <= 1) _counter++;
    }
    temp[threadIdx.x] = _counter;
    
    __syncthreads();

    if(threadIdx.x == THREADS - 1) {
      double block_result = 0.0;
      for(int i = 0; i < THREADS; i++) block_result += temp[i];
      atomicAdd(&(numVector[blockIdx.x]), block_result);
    }

}

int main(int argc, char *argv[]) {
  srand(time(0));
  clock_t start_time, end_time;
  double elapsed = 0;

  long long int BLOCKS = 16384, SAMPLES = 0, LOOPS = 1;

  for(int i = 0; i < argc; i++) {
    if(strcmp(argv[i],"-g") == 0) {
      BLOCKS = atoi(argv[++i]);
    } else if(strcmp(argv[i],"-l") == 0) {
      LOOPS = atoi(argv[++i]);
    } else if(strcmp(argv[i],"-s") == 0) {
      SAMPLES = atoi(argv[++i]);
    }
  }

  if(SAMPLES > 0) {
    LOOPS = max(SAMPLES / BLOCKS / (long long int)THREADS, LOOPS);
  }

  //printf("%lld %lld %lld %lld",SAMPLES, LOOPS, BLOCKS, (long long int)THREADS);

  int N = BLOCKS * THREADS, seed = rand();

  start_time = clock();
  printf("Allocating device and host memory for vectors..\n");
  double *hostVector, *devVector;
  hostVector = (double*) malloc(BLOCKS * sizeof(double));
  memset(hostVector, 0, BLOCKS * sizeof(double));

  cudaMalloc(&devVector, BLOCKS * sizeof(double));
  cudaMemcpy(devVector, hostVector, BLOCKS * sizeof(double), cudaMemcpyHostToDevice);
  end_time = clock();
  printf("Done! in %.2f ms\n\n", ((double) (end_time - start_time)) / CLOCKS_PER_SEC * 1000.0);
  elapsed += ((double) (end_time - start_time)) / CLOCKS_PER_SEC * 1000.0;

  start_time = clock();
  printf("Allocating device states for RNG...\n");
  curandState* devStates;
  cudaMalloc(&devStates, N * sizeof(curandState));

  setup_kernel<<<BLOCKS, THREADS>>>(devStates, seed);
  end_time = clock();
  printf("Done! in %.2f ms\n\n", ((double) (end_time - start_time)) / CLOCKS_PER_SEC * 1000.0);
  elapsed += ((double) (end_time - start_time)) / CLOCKS_PER_SEC * 1000.0;

  start_time = clock();
  printf("Executing PI approximation with %ld samples...\n", (long int)BLOCKS * (long int)THREADS * (long int)LOOPS);
  approximate<<<BLOCKS, THREADS>>>(devVector, devStates, LOOPS);
  cudaDeviceSynchronize();

  cudaMemcpy(hostVector, devVector, BLOCKS * sizeof(double), cudaMemcpyDeviceToHost);
  end_time = clock();
  printf("Done! in %.2f ms\n\n", ((double) (end_time - start_time)) / CLOCKS_PER_SEC * 1000.0);

  long double approx = 0;

  for(int i = 0; i < BLOCKS; i++) approx += hostVector[i];
  printf("Approximated PI: %.20Lf\nPredefined   PI: 3.14159265358979323846\n\n", approx / ((long double)BLOCKS * (long double)THREADS * (long double)LOOPS) * 4);

  elapsed += ((double) (end_time - start_time)) / CLOCKS_PER_SEC * 1000.0;

  printf("Terminating execution, overall elapsed time: %.2f ms\n\n", elapsed);

}