A PI Estimation example using plain C++ imperative abstractions

pi_estimation_plain_imperative.cpp

//
//  Created by Kjell-Olov Högdahl on 2017-03-08.
//  This work is licensed under
//  Creative Commons Attribution 4.0 International License.
//  https://creativecommons.org/licenses/by/4.0/
//

#include <iostream>
#include <cmath> // std::pow
#include <random> // random_device, std::mt19937, uniform_int_distribution

const auto PI_REFERENCE = 3.14159265359;
struct Fixture {double radius; long point_count;};
struct PiEstimateRecord {
    Fixture fixture;
    double pi_estimate;
    double deviation(const double pi_reference) const {
        return std::abs(pi_estimate-pi_reference);
    }
    bool operator<(const PiEstimateRecord& other) const {
        return std::abs(pi_estimate - PI_REFERENCE)
        < std::abs(other.pi_estimate - PI_REFERENCE);
    }
};
using PiEstimateRecords = std::vector<PiEstimateRecord>;
void print(const PiEstimateRecord& pi_estimate_record) {
    std::cout
    << '\n'
    << "PI Estimate"
    << "(radius:" << pi_estimate_record.fixture.radius
    << ",point count:" <<  pi_estimate_record.fixture.point_count
    << ") = " << pi_estimate_record.pi_estimate
    << " diff from " << PI_REFERENCE
    << " = " << pi_estimate_record.deviation(PI_REFERENCE);
}
void print(const PiEstimateRecords& pi_estimate_records) {
    for (auto& pi_estimate_record : pi_estimate_records) {
        print(pi_estimate_record);
    }
}

namespace plain_imperative {
    template <int N,int M>
    PiEstimateRecords do_pi_estimations(std::mt19937& random_generator) {
        PiEstimateRecords result;
        for (int n = 0;n<=N;++n) {
            for (int m = 0;m<=M;++m) {
                // radius:10^0..10^N, point_count:10^0..10^M
                double radius = std::pow(10,n);
                long point_count = std::pow(10,m);
                // place origo at center of square/circle
                std::uniform_real_distribution<> distribution(-radius, radius);
                long points_in_circle_count = 0;
                // Count number of random points within circle radius
                for (long i=0;i<point_count;++i) {
                    // random x,y in square -radius...radious
                    double x = distribution(random_generator);
                    double y = distribution(random_generator);
                    // Use pythagoras to tell points within radius
                    if ((x*x + y*y) <= radius*radius) {
                        points_in_circle_count += 1;
                    }
                };
                // Estimate Based on Circle Area / Square Area = Pi/4
                double pi_estimate =
                4.0
                * static_cast<double>(points_in_circle_count)
                / point_count;
                PiEstimateRecord entry{{radius,point_count},pi_estimate};
                result.push_back(entry);
                print(entry);
            }
        }
        return result;
    }
}

// See http://www.fluentcpp.com/2017/03/02/the-pi-day-challenge-for-expressive-code-in-c/
int main(int argc, const char * argv[]) {
    std::cout << "\nHello from PI Estimator!";
    std::random_device random_device;
    std::mt19937 random_generator(random_device());
    const int N = 9;
    const int M = 7;
    print(plain_imperative::do_pi_estimations<N,M>(random_generator));
    std::cout << "\nBye from PI Estimator!";
    return 0;
}