Paul Dreik pi estimation

piestimate_pauldreik.cpp

/*
 * Entry for the pi estimate blog entry, see
 * http://www.fluentcpp.com/2017/03/02/the-pi-day-challenge-for-expressive-code-in-c/
 * By Paul Dreik 20170303
 * License: BOOST 1.0
 */
#include <random>
#include <iostream>
#include <iomanip>
#include <cmath>

//the underlying type for floating point calculations
using Real=double;

//do we have constexpr support for math functions?
#if defined(__GNUC__) && !defined(__clang__)
#  define CONSTEXPR constexpr
#else
#  define CONSTEXPR
#endif

//official value of pi, to compare with
CONSTEXPR const auto pi=std::atan(Real{1})*4;

/*
 * Estimates pi by random sampling. We will only sample a small part
 * of the circle, to improve the estimate for the same amount of cpu
 * time. The area to be sampled has area (1/sqrt(2)-1/2)r^2 instead of
 * 4r^2, which is approximately 5%.
 */
Real estimatePi(const Real radius, const std::size_t Npoints) {
  
  CONSTEXPR const Real sqrt2=std::sqrt(Real{2});
  CONSTEXPR const Real sqrt2inv=1/sqrt2;
  
  std::random_device rd;
  std::mt19937 gen(rd());    
  std::uniform_real_distribution<Real> xdist(0, radius*sqrt2inv);
  std::uniform_real_distribution<Real> ydist(radius*sqrt2inv,radius);
  
  std::size_t insidecount=0;
  const Real radius2=radius*radius;    
  for (std::size_t i=0; i<Npoints; ++i) {
    const auto x=xdist(gen);
    const auto y=ydist(gen);
    if(x*x+y*y < radius2) {
      ++insidecount;
    }
  }
  const Real insidefraction=static_cast<Real>(insidecount)/Npoints;
  const Real piest= insidefraction*4*(sqrt2-1)+2;
  return piest;
}

int main()
{
  const bool verbose=false;

  for(int n=0; n<=9; ++n) {
    const Real radius=std::pow(Real{10},n);
    for(int m=0; m<=7; ++m) {
      const std::size_t Npoints=std::pow(10,m);
      const Real piest=estimatePi(radius,Npoints);
      const Real err=std::abs(piest-pi);
      if(verbose) {
        std::cout << "radius="<<radius<<" Npoints="<<Npoints<<" error=";
      }
      std::cout<<std::setw(16)<<std::setfill(' ') << err;
    }
    std::cout<<'\n';
  }
}