MC calculation of Pi (OpenACC)

pi_acc.cc

#include <iostream>
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
//----------------------------------------------------------------------
const int N = 10000;
const int SAMPLE = 65535;
//----------------------------------------------------------------------
double
myrand(void){
  return static_cast<double>(rand())/static_cast<double>(RAND_MAX);
}
//----------------------------------------------------------------------
double
one_trial(void){
  int sum = 0;
  for(int i=0;i<N;i++){
    const double x = myrand();
    const double y = myrand();
    if (x*x+y*y < 1.0) sum ++;
  }
  return static_cast<double>(sum)/static_cast<double>(N)*4.0;
}
//----------------------------------------------------------------------
double r[2*N];
//----------------------------------------------------------------------
double
one_trial_acc(){
  int sum = 0;
#pragma acc kernels pcopyin(r[1:2*N])
#pragma acc loop reduction(+:sum)
  for(int i=0;i<N;i++){
    const double x = r[i*2];
    const double y = r[i*2+1];
    if (x*x+y*y < 1.0) sum ++;
  }
  return static_cast<double>(sum)/static_cast<double>(N)*4.0;
}

//----------------------------------------------------------------------
int
main(void){
  double a = 0.0;
  double v = 0.0;
  for(int i=0;i<SAMPLE;i++){
    for(int i=0;i<2*N;i++){
      r[i] = myrand();
    }
    const double p = one_trial_acc();
    a += p;
    v += p*p;
  }
  const double s = static_cast<double>(SAMPLE);
  a /= s;
  v /= s;
  v = sqrt((v -a*a)/(s-1.0));
  std::cout << a << std::endl;
  std::cout << v << std::endl;
}