MC calculation of Pi (OpenMP)

pi_omp.cc

#include <iostream>
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <random>
//----------------------------------------------------------------------
const int N = 10000;
const int SAMPLE = 65536;
//----------------------------------------------------------------------
double
myrand(void){
  return static_cast<double>(rand())/static_cast<double>(RAND_MAX);
}
//----------------------------------------------------------------------
double
one_trial(int index){
  std::random_device rd;
  std::mt19937 g(rd());
  std::uniform_real_distribution<double> dist(0,1);
  g.seed(index);
  int sum = 0;
  for(int i=0;i<N;i++){
    const double x = dist(g);
    const double y = dist(g);
    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;
  double r[SAMPLE];
#pragma omp parallel for
  for(int i=0;i<SAMPLE;i++){
    r[i] =  one_trial(i);
  }
  for(int i=0;i<SAMPLE;i++){
    a += r[i];
    v += r[i]*r[i];
  }
  const double s = static_cast<double>(SAMPLE);
  a /= s;
  v /= s;
  v = sqrt((v -a*a)/(s-1.0));
  std::cout << "Total Samples " << N*SAMPLE << std::endl;
  std::cout << a << std::endl;
  std::cout << v << std::endl;
}