gistfile1.c

#include <stdio.h>
#include <omp.h>

#include <sys/time.h>

#ifdef RDRAND
typedef char seed_t;
static inline unsigned int
myrand(seed_t *seed)
{
  int rdrand;
  __asm__ __volatile__("\n"
                       "0:\n"
                       "\trdrand %0\n\t"
                       "\tjnc 0b\n\t" :
                       "=r"(rdrand));
  return rdrand;
}
const unsigned long long max2 = 0xffffffffull * 0xffffffffull;
#else
#ifdef SFMT_MEXP
#include "SFMT.h"

typedef sfmt_t seed_t;
#define myrand sfmt_genrand_uint32
const unsigned long long max2 = 0xffffffffull * 0xffffffffull;
#else
#include <stdlib.h>

typedef unsigned int seed_t;
#define myrand rand_r
const unsigned long long max2 = (unsigned long long) RAND_MAX * RAND_MAX;
#endif
#endif

int
main(int argc, char **argv)
{
  long long li, lpi;
  const long long ln = 0x1000000000ll;
  long double pi;
  struct timeval start, end;

  gettimeofday(&start, NULL);
  lpi = 0;
#pragma omp parallel
  {
    seed_t seed;
#ifndef RDRAND
#ifdef SFMT_MEXP
    sfmt_init_gen_rand(&seed, omp_get_thread_num());
#else
    seed = omp_get_thread_num();
#endif
#endif
#pragma omp for reduction(+:lpi)
    for (li = 0; li < ln; li++) {
      unsigned long long ux, uy;
      ux = myrand(&seed);
      uy = myrand(&seed);
      if (ux * ux < (unsigned long long) (max2 - (uy * uy))) {
        lpi++;
      }
    }
  }
  gettimeofday(&end, NULL);
  pi = (long double) lpi / ln * 4.0l;
  printf("Calc PI: %.60Lf\n", pi);
  printf("Real PI: "
         "3.141592653589793238462643383279502884197169399375105820974944\n");
  printf("Time:    %06f sec\n",
         end.tv_sec - start.tv_sec + (end.tv_usec - start.tv_usec) * 1.0e-6);
  return 0;
}