# Multithreaded Monte Carlo pi approximation with own pseudorandom number generator

I made a Monte Carlo pi approximation program, that makes use of multithreading and a pseudorandom number generator I wrote (the one from big_wheel.hpp, which I tested using TestU01).

### Code:

monte_carlo_pi.cpp:

#include <iostream>
#include <iomanip>
#include <cstdlib>
#include <cstdint>
#include <cmath>
#include <mutex>
#include <random>
#include <vector>
#include <chrono>
#include "big_wheel.hpp"

typedef struct
{
double x, y, z;
} point_t;

static const double PI = 3.1415926535897932384626433832795;
static const std::uint32_t prng_max_value = big_wheel_engine<>::max();

double approximate_pi(std::uint32_t pc,std::uint32_t tp)
{
return 4.0*(double(pc)/double(tp));
}

double pi_error(double x)
{
return 100.0*((x-PI)/PI);
}

std::mutex mtx;

// Points within circle
std::uint32_t pcnt = 0;

template <class RandEng>
void calculate_pi(std::uint32_t nppt,RandEng rgen)
{
std::lock_guard<std::mutex> guard(mtx);

point_t rpt;
for (std::uint32_t i=0;i<nppt;i++)
{
rpt.x = double(rgen())/static_cast<double>(prng_max_value);
rpt.y = double(rgen())/static_cast<double>(prng_max_value);
rpt.z = (rpt.x*rpt.x)+(rpt.y*rpt.y);
if (rpt.z <= 1.0) pcnt++;
}
}

int main(int argc,char **argv)
{
// Check number of cores

// Generate seeds
std::random_device rd;
std::uint32_t seed[16] = {0x0};

for (int i=0;i<16;i++)
seed[i] = rd();

// Initialize PRNG
big_wheel_engine<> rnd(
seed[ 0],seed[ 1],seed[ 2],seed[ 3],
seed[ 4],seed[ 5],seed[ 6],seed[ 7],
seed[ 8],seed[ 9],seed[10],seed[11],
seed[12],seed[13],seed[14],seed[15]
);

std::uint32_t num_points = 3221225472;

for (std::size_t i=0;i<num_cores;i++)
calculate_pi<big_wheel_engine<>>,
rnd
)
);

// Calculated Pi value
double gen_pi = 0;

// Calculation error
double pi_err = 0;

// Duration
double dur = 0;

// Print info
std::cout << "Seeds:\n";
for (int iy=0;iy<4;iy++)
{
std::cout << " ";
for (int ix=0;ix<4;ix++)
std::cout << std::setw(11) << std::setfill(' ') << seed[4*iy+ix] << " ";
std::cout << "\n";
}
std::cout << "Number of points: " << num_points << "\n";
std::cout << "Number of cores: " << num_cores << "\n";
std::cout << " ==> Number of points per thread: " << num_points_per_thread << "\n\n";

std::cout << "Generating random points...\n";

std::chrono::high_resolution_clock::time_point tpa = std::chrono::high_resolution_clock::now();

for (std::size_t i=0;i<num_cores;i++)

std::chrono::high_resolution_clock::time_point tpb = std::chrono::high_resolution_clock::now();

std::chrono::duration<double> tm = std::chrono::duration_cast<std::chrono::duration<double>>(tpb-tpa);
dur = tm.count();

// Approximate pi and calculate error
gen_pi = approximate_pi(pcnt,num_points);
pi_err = pi_error(gen_pi);

std::cout << "\nResults:\n";
std::cout << "  Number of points within circle: " << pcnt << "\n";
std::cout << "  Pi approximation:               " << std::setprecision(32) << gen_pi << "\n";
std::cout << "  Error:                          " << std::setprecision(10) << pi_err << "%\n";

std::cout << "Time elapsed: " << std::setprecision(10) << dur << " seconds\n";
}

big_wheel.hpp:

#ifndef BIG_WHEEL_HPP
#define BIG_WHEEL_HPP

#include <iostream>
#include <cstdlib>
#include <cstdint>
#include <climits>
#include <limits>

template <
std::uint32_t m = 0x3655accc,
std::size_t rs1 = 11,std::size_t rs2 = 20,std::size_t rs3 = 14,std::size_t rs4 = 18,
std::size_t rm1 = 16,std::size_t rm2 = 16,
std::size_t rf = 7
>
class big_wheel_engine
{
public:
typedef std::uint32_t result_type;

static constexpr std::size_t   state_size = 4096;
static constexpr std::uint32_t mix_mask  = m;
static constexpr std::size_t   rot_sel1  = rs1;
static constexpr std::size_t   rot_sel2  = rs2;
static constexpr std::size_t   rot_sel3  = rs3;
static constexpr std::size_t   rot_sel4  = rs4;
static constexpr std::size_t   rot_mix1  = rm1;
static constexpr std::size_t   rot_mix2  = rm2;
static constexpr std::size_t   rot_final = rf;

big_wheel_engine(
const result_type& s01 = 0x0,const result_type& s02 = 0x0,const result_type& s03 = 0x0,const result_type& s04 = 0x0,
const result_type& s05 = 0x0,const result_type& s06 = 0x0,const result_type& s07 = 0x0,const result_type& s08 = 0x0,
const result_type& s09 = 0x0,const result_type& s10 = 0x0,const result_type& s11 = 0x0,const result_type& s12 = 0x0,
const result_type& s13 = 0x0,const result_type& s14 = 0x0,const result_type& s15 = 0x0,const result_type& s16 = 0x0)
{
m_i = ~m;

seed(s01,s02,s03,s04,
s05,s06,s07,s08,
s09,s10,s11,s12,
s13,s14,s15,s16);
}

big_wheel_engine(const big_wheel_engine<m,rs1,rs2,rs3,rs4,rm1,rm2,rf>& x)
{
for (std::size_t i=0;i<state_size;i++)
state[i] = x.state[i];
counter = x.counter;
}

big_wheel_engine<m,rs1,rs2,rs3,rs4,rm1,rm2,rf>& operator=(const big_wheel_engine<m,rs1,rs2,rs3,rs4,rm1,rm2,rf>& x)
{
for (std::size_t i=0;i<state_size;i++)
state[i] = x.state[i];
counter = x.counter;
return *this;
}

void seed(
const result_type& s01 = 0x0,const result_type& s02 = 0x0,const result_type& s03 = 0x0,const result_type& s04 = 0x0,
const result_type& s05 = 0x0,const result_type& s06 = 0x0,const result_type& s07 = 0x0,const result_type& s08 = 0x0,
const result_type& s09 = 0x0,const result_type& s10 = 0x0,const result_type& s11 = 0x0,const result_type& s12 = 0x0,
const result_type& s13 = 0x0,const result_type& s14 = 0x0,const result_type& s15 = 0x0,const result_type& s16 = 0x0)
{
counter = 0;
p1 = (counter+  0)%4096; p2 = (counter+  8)%4096;
p3 = (counter+ 18)%4096; p4 = (counter+ 30)%4096;
p5 = (counter+ 44)%4096; p6 = (counter+ 60)%4096;
p7 = (counter+ 78)%4096; p8 = (counter+ 98)%4096;
p9 = (counter+120)%4096;

// Based on decimal part of sqrt(61/19)
result_type off1[16] = {
0xcab305ab,0x22132737,0xbda768dd,0x2438d3ab,
0x00f4ca5a,0xd4681d61,0xc6420bc7,0x2828bdba,
0x64806a32,0x2fa2391f,0xbc61cbca,0xcc6eb7a0,
0xfa01a2cd,0xb766ce25,0x5499507a,0x19c839a0
};
// Based on decimal part of sqrt(73/29)
result_type off2[16] = {
0x962a374b,0x0fc75efe,0x86dfdeff,0xddb14a03,
0xf61361c6,0xc908fde9,0x27df3b7d,0x04ae7492,
0x5ee38142,0xfb21ee51,0x56fae74d,0x4f688f90,
0x7889f34c,0x44e4aa33,0x547e3778,0xb6201d7b
};
// Based on decimal part of sqrt(79/31)
result_type off3[16] = {
0x98ab7f56,0x38ca9f10,0xa2f3d406,0x322c4113,
0xb8d42dd6,0xe6490633,0x4c5bdcf7,0x454c3122,
0x7b22e315,0x1a205e44,0xc6d1b31b,0x06eb7dea,
0x875f1cf1,0xbe9bec2a,0x0350aa26,0x8b3d5aed
};

result_type sset[16] = {
s01,s02,s03,s04,
s05,s06,s07,s08,
s09,s10,s11,s12,
s13,s14,s15,s16
};

result_type st_a[16] = {0x0};
result_type st_b[256] = {0x0};
result_type st_c[4096] = {0x0};
result_type rn_1[4096] = {0x0};
result_type rn_2[4096] = {0x0};
result_type rn_3[4096] = {0x0};
result_type rn_4[4096] = {0x0};

// Expansion
for (std::size_t i=0;i<16;i++)
{
st_a[i] = sset[i]+off1[i];
st_a[i] = smix(st_a[i],sbox_i1);
}

for (std::size_t i=0;i<16;i++)
{
for (std::size_t j=0;j<16;j++)
{
st_b[16*i+j] = st_a[i]+off2[j];
st_b[16*i+j] = smix(st_b[16*i+j],sbox_i2);
}
}

for (std::size_t i=0;i<256;i++)
{
for (std::size_t j=0;j<16;j++)
{
st_c[16*i+j] = st_b[i]+off3[j];
st_c[16*i+j] = smix(st_c[16*i+j],sbox_i3);
}
}

// Mixing
for (std::size_t i=0;i<4096;i++)
{
rn_1[i] = (st_c[(i+ 0)%4096] << 13) | (st_c[(i+ 4)%4096] >> 19);
rn_1[i] = smix(rn_1[i],sbox_i4);
}
for (std::size_t i=0;i<4096;i++)
{
rn_2[i] = (rn_1[(i+10)%4096] << 11) | (rn_1[(i+18)%4096] >> 21);
rn_2[i] = smix(rn_2[i],sbox_i4);
}
for (std::size_t i=0;i<4096;i++)
{
rn_3[i] = (rn_2[(i+28)%4096] << 17) | (rn_2[(i+40)%4096] >> 15);
rn_3[i] = smix(rn_3[i],sbox_i4);
}
for (std::size_t i=0;i<4096;i++)
{
rn_4[i] = (rn_3[(i+54)%4096] << 10) | (rn_3[(i+70)%4096] >> 22);
rn_4[i] = smix(rn_4[i],sbox_i4);
}

for (std::size_t i=0;i<4096;i++)
{
state[i] = (rn_4[(i+88)%4096] << 16) | (rn_4[(i+108)%4096] >> 16);
state[i] = smix(state[i],sbox_i5);
}
}

// Get next random number
result_type operator()()
{
// Elements from state to mix
p1 = (counter+  0)%4096; p2 = (counter+  8)%4096;
p3 = (counter+ 18)%4096; p4 = (counter+ 30)%4096;
p5 = (counter+ 44)%4096; p6 = (counter+ 60)%4096;
p7 = (counter+ 78)%4096; p8 = (counter+ 98)%4096;
p9 = (counter+120)%4096;

// Start mixing
tmp_11 = (state[p1] << rs1) | (state[p2] >> (32-rs1));
tmp_12 = (state[p3] << rs2) | (state[p4] >> (32-rs2));
tmp_13 = (state[p5] << rs3) | (state[p6] >> (32-rs3));
tmp_14 = (state[p7] << rs4) | (state[p8] >> (32-rs4));

tmp_21 = ((tmp_11 << rm1) | (tmp_12 >> (32-rm1)));
tmp_22 = ((tmp_13 << rm2) | (tmp_14 >> (32-rm2)));

tmp_3 = (tmp_21 & m) | (tmp_22 & m_i);

tmp_4 = rot_r(state[p9] ^ tmp_3,rf);

state[p9] = smix(tmp_4,sbox_r);

// For jumping to next elements
counter++;

return state[p9];
}

{
for (unsigned long long i=0;i<z;i++)
(*this)();
}

// Minimum and maximum
static constexpr result_type min()
{
return std::numeric_limits<result_type>::min();
}

static constexpr result_type max()
{
return std::numeric_limits<result_type>::max();
}

// Comparision
friend bool operator==(
const big_wheel_engine<m,rs1,rs2,rs3,rs4,rm1,rm2,rf>& lhs,
const big_wheel_engine<m,rs1,rs2,rs3,rs4,rm1,rm2,rf>& rhs)
{
if (lhs.counter != rhs.counter)
return false;

for (std::size_t i=0;i<state_size;i++)
{
if (lhs.state[i] != rhs.state[i])
return false;
else
continue;
}

return true;
}

friend bool operator!=(
const big_wheel_engine<m,rs1,rs2,rs3,rs4,rm1,rm2,rf>& lhs,
const big_wheel_engine<m,rs1,rs2,rs3,rs4,rm1,rm2,rf>& rhs)
{
return !(lhs == rhs);
}

// Stream ops
template <class CharT,class Traits>
friend std::basic_ostream<CharT,Traits>& operator<<(
std::basic_ostream<CharT,Traits>& os,
const big_wheel_engine<m,rs1,rs2,rs3,rs4,rm1,rm2,rf>& bw)
{
os << bw.counter;
for (std::size_t i=0;i<bw.state_size;i++)
os << ' ' << bw.state[i];

return os;
}

template <class CharT,class Traits>
friend std::basic_istream<CharT,Traits>& operator>>(
std::basic_istream<CharT,Traits>& is,
const big_wheel_engine<m,rs1,rs2,rs3,rs4,rm1,rm2,rf>& bw)
{
is >> bw.counter >> std::ws;
for (std::size_t i=0;i<bw.state_size;i++)
is >> bw.state[i] >> std::ws;

return is;
}
private:
// State
result_type state[4096];

result_type m_i;

// Temporal variables
result_type tmp_11, tmp_12, tmp_13, tmp_14;
result_type tmp_21, tmp_22;
result_type tmp_3;
result_type tmp_4;

// For picking elements from state
std::size_t counter; // Counter
std::size_t p1, p2, p3, p4, p5, p6, p7, p8, p9; // Positions

// S-boxes:
// - sbox_i{1,2,3,4,5}: For initialization
// -- i{1,2,3}: Seed expansion
// -- i{4,5}:   Mixing
// - sbox_r: For generating the next random number
static std::uint32_t sbox_i1[256];
static std::uint32_t sbox_i2[256];
static std::uint32_t sbox_i3[256];
static std::uint32_t sbox_i4[256];
static std::uint32_t sbox_i5[256];
static std::uint32_t sbox_r[256];

// Left/right bit rotation
template <typename D> D rot_l(D x,unsigned r) { return ((x << r) | (x >> (CHAR_BIT*sizeof(D)-r))); }
template <typename D> D rot_r(D x,unsigned r) { return ((x >> r) | (x << (CHAR_BIT*sizeof(D)-r))); }

// Mix x with an S-box sb
result_type smix(result_type x,std::uint32_t *sb)
{
return result_type(
(sb[(x >> 24)&0xff] << 24)|
(sb[(x >> 16)&0xff] << 16)|
(sb[(x >>  8)&0xff] <<  8)|
(sb[(x >>  0)&0xff] <<  0)
);
}
};

// The S-boxes (randomly generated)
template <
std::uint32_t m,
std::size_t rs1,std::size_t rs2,std::size_t rs3,std::size_t rs4,
std::size_t rm1,std::size_t rm2,
std::size_t rf
>
std::uint32_t big_wheel_engine<m,rs1,rs2,rs3,rs4,rm1,rm2,rf>::sbox_i1[256] = {
0x36,0x09,0x2d,0x7f,0x28,0x8a,0x51,0x81,0x2b,0x69,0x90,0x58,0xaf,0x21,0xda,0x71,
0x49,0x6a,0x94,0x0a,0xaa,0x30,0xf3,0x73,0x37,0xfa,0xcf,0x7e,0x76,0xdf,0xb1,0x56,
0xd5,0x78,0x12,0x08,0xe3,0xf1,0x07,0x9d,0x98,0xc5,0xe6,0xcd,0x34,0xf8,0xfd,0xd0,
0x2f,0x84,0x66,0x4b,0x3c,0x4e,0xc6,0x63,0x99,0x57,0x60,0x54,0x97,0x9f,0x39,0xf0,
0x59,0xa1,0xcb,0x18,0xed,0x0b,0xd1,0xdb,0x44,0x1e,0x0f,0x72,0x2c,0x5b,0x27,0xc9,
0x55,0x1c,0x41,0xfc,0x11,0x50,0x4d,0x9b,0x88,0xea,0xd6,0xc8,0x45,0xa2,0xcc,0x8c,
0xa6,0x61,0xb5,0x01,0xb4,0x6d,0x93,0xc3,0xab,0xdd,0x9e,0x6c,0x38,0xe7,0x47,0xe9,
0x70,0x31,0xe8,0x03,0x52,0x75,0xe0,0xb6,0x53,0x40,0x46,0x0e,0x3a,0x5e,0x79,0x5c,
0xc2,0x7c,0x7a,0x6f,0xe1,0xd4,0xb8,0x85,0xee,0x22,0xb0,0x80,0x42,0xf6,0xec,0x24,
0xe2,0x9a,0x14,0x8f,0xac,0x15,0x10,0x82,0xa3,0x20,0x3b,0xf2,0x68,0x4c,0xd3,0x9c,
0xca,0x48,0xfe,0x00,0x04,0xc0,0x0c,0xba,0xb7,0xc7,0x3e,0xff,0xd8,0xfb,0x6b,0x65,
0x26,0x17,0x8b,0xf9,0x02,0xb2,0xe4,0xbf,0x8e,0xc1,0x62,0xd7,0x1d,0x77,0x2e,0x96,
0x25,0x33,0x5f,0xef,0x35,0x43,0xf5,0xf7,0xa0,0x64,0x05,0xc4,0x1b,0xd2,0x0d,0x16,
0xdc,0x86,0x7d,0xeb,0x87,0xa5,0x83,0xf4,0x23,0xd9,0x4a,0x95,0x1f,0x7b,0x74,0xa9,
0x2a,0x13,0x06,0x29,0xae,0xce,0x89,0xa7,0x3d,0x19,0x67,0xa8,0x32,0xb3,0xbe,0xe5
};

template <
std::uint32_t m,
std::size_t rs1,std::size_t rs2,std::size_t rs3,std::size_t rs4,
std::size_t rm1,std::size_t rm2,
std::size_t rf
>
std::uint32_t big_wheel_engine<m,rs1,rs2,rs3,rs4,rm1,rm2,rf>::sbox_i2[256] = {
0x2d,0xa9,0x38,0x42,0x24,0xd1,0x55,0x8a,0x3c,0x22,0x84,0x7a,0xd6,0xbc,0x2a,0xb7,
0x99,0xd7,0xe6,0x2f,0x9e,0x1f,0x1e,0x0a,0xd0,0x97,0x45,0xe9,0x90,0x8f,0xe7,0x35,
0x40,0x6e,0x9a,0xeb,0x28,0x75,0x03,0x5a,0x29,0x34,0x6b,0x17,0xaa,0x4a,0x6d,0x00,
0x36,0xd8,0x71,0x9b,0xf1,0x25,0x5b,0x46,0x12,0x69,0x9d,0xb9,0x66,0xf7,0x74,0x95,
0x61,0x48,0x1a,0x23,0xcb,0x5c,0xcd,0x70,0x18,0xb3,0x76,0x0e,0x37,0xc8,0x3a,0x92,
0x59,0xdd,0x26,0xaf,0xb0,0x8e,0x44,0xce,0xe3,0xcf,0x4f,0x0d,0x4c,0x19,0x49,0xca,
0x07,0xc7,0xd2,0x72,0x7b,0x3e,0xc5,0xc1,0xc2,0x91,0x86,0x4b,0x57,0xe0,0x02,0xb6,
0xf4,0x7e,0x39,0xd9,0x20,0x65,0xac,0x87,0xb8,0x33,0xe5,0xb2,0xa1,0x96,0x64,0xa6,
0x2b,0x01,0x85,0x50,0xed,0x27,0xff,0x62,0x30,0x80,0xc9,0xa4,0x31,0xde,0x21,0x7c,
0xdc,0xf6,0xae,0x52,0x8c,0x51,0x32,0xbf,0x0c,0xa8,0xb4,0xe2,0xe8,0x94,0x14,0x43,
0xdb,0x54,0xa3,0x47,0x2e,0xfd,0x67,0x53,0x60,0xee,0x6a,0xa2,0xd4,0xfe,0x8b,0x82,
0x4e,0xc6,0xc4,0x98,0xb5,0xbe,0x06,0x3d,0x13,0xf8,0xe1,0x3b,0x81,0xa7,0x79,0x8d,
0xf2,0xb1,0xd5,0x2c,0xc3,0x1d,0x1b,0x73,0x9f,0x6f,0xba,0xcc,0xfc,0xea,0x09,0xa0,
0xec,0x10,0xc0,0x4d,0x11,0x0f,0x88,0x08,0x5f,0x56,0xf0,0x5e,0x63,0xfa,0x05,0x78,
0x3f,0x7d,0x41,0xfb,0x68,0xa5,0x89,0xf3,0xef,0x83,0x6c,0x16,0xe4,0xda,0x9c,0x93
};

template <
std::uint32_t m,
std::size_t rs1,std::size_t rs2,std::size_t rs3,std::size_t rs4,
std::size_t rm1,std::size_t rm2,
std::size_t rf
>
std::uint32_t big_wheel_engine<m,rs1,rs2,rs3,rs4,rm1,rm2,rf>::sbox_i3[256] = {
0x06,0x38,0x28,0x6a,0x51,0xb1,0x0a,0x55,0xf9,0xa3,0x73,0x11,0xed,0x3a,0x43,0x83,
0x54,0x74,0x89,0xe0,0x85,0x10,0xcc,0x6f,0x64,0x2b,0x71,0xc1,0xa4,0x5b,0x42,0x92,
0xc8,0x99,0x4a,0x8f,0xda,0xbc,0xb7,0x5c,0x81,0x16,0x52,0x50,0x82,0x61,0x3f,0xf5,
0x1b,0x12,0x90,0x98,0x02,0x05,0x8d,0xa6,0x5a,0xbf,0xe9,0xa7,0x24,0x56,0xfd,0x40,
0x33,0xe6,0x17,0x72,0xf8,0x18,0xd1,0xac,0xa2,0xfe,0x0f,0x88,0x39,0x09,0x6c,0x0e,
0x67,0x7c,0x84,0xeb,0x8a,0xba,0x63,0x03,0x75,0xcd,0x93,0x97,0xab,0x6d,0x76,0xdc,
0xcb,0xc3,0xb4,0x9c,0xa9,0xd8,0x9d,0xb5,0xca,0x7a,0x1a,0xc4,0x0d,0x19,0xc6,0x30,
0xf0,0x5d,0x65,0x79,0x7e,0x2c,0x31,0x13,0x9e,0xbd,0x44,0xbe,0x57,0xb2,0xc5,0x80,
0xe3,0x27,0xd0,0xec,0x48,0x4e,0xae,0x86,0x8c,0xf4,0x3b,0xfb,0xb3,0xa8,0x69,0xb8,
0x1f,0x0b,0xcf,0x4c,0x45,0xaa,0xdf,0xf7,0xe8,0x68,0xef,0x58,0x34,0x7d,0xbb,0xf6,
0x2e,0x22,0xdb,0x35,0xd6,0x6b,0x1d,0x5e,0xc0,0xa1,0x3d,0xfa,0xb0,0xff,0x60,0x2a,
0x41,0x5f,0xc7,0xd5,0x4d,0x1e,0x2f,0x01,0xea,0x49,0xdd,0x26,0xb6,0x78,0xde,0x25,
0x36,0xf1,0x0c,0xe1,0x04,0x32,0x14,0x3e,0x87,0x6e,0x96,0xe5,0xce,0x95,0xe2,0xaf,
0xf3,0xa0,0x46,0x37,0x15,0x1c,0x00,0x2d,0x7f,0xe4,0x62,0xee,0x3c,0x08,0x47,0x9b,
0x66,0xfc,0x21,0xd3,0x94,0x59,0x8b,0x07,0xd4,0x29,0xc9,0x23,0xc2,0x9a,0xe7,0xd9
};

template <
std::uint32_t m,
std::size_t rs1,std::size_t rs2,std::size_t rs3,std::size_t rs4,
std::size_t rm1,std::size_t rm2,
std::size_t rf
>
std::uint32_t big_wheel_engine<m,rs1,rs2,rs3,rs4,rm1,rm2,rf>::sbox_i4[256] = {
0x1a,0x6b,0xfd,0xd0,0x0d,0x11,0x50,0x9e,0x1f,0xd3,0x95,0xca,0xbc,0x0c,0x80,0x29,
0x47,0xa4,0x3f,0x3c,0x2c,0x76,0x3a,0x91,0x60,0x8d,0x39,0x9b,0x7e,0x6d,0x1b,0x12,
0x7f,0xf3,0xab,0xb6,0x63,0x7b,0x53,0x31,0xd6,0x4d,0xbb,0xcb,0x56,0xda,0x15,0xbe,
0x67,0xfb,0xa9,0xdb,0x69,0x97,0x42,0xf1,0x5c,0xc5,0xb3,0x34,0xa3,0x6f,0x51,0x4f,
0xc2,0x70,0xf8,0x52,0x92,0x81,0xd1,0x0f,0xba,0x3d,0x4c,0xb0,0x48,0xff,0xd2,0x7d,
0xaf,0x1c,0x65,0xe8,0x4b,0x5f,0x94,0xcc,0x84,0x09,0x64,0x2e,0xb1,0x77,0x96,0xea,
0x79,0x24,0xb8,0x0e,0x90,0x33,0xf9,0xb2,0x57,0x28,0xe9,0x04,0x03,0xa1,0x17,0xe7,
0x82,0xaa,0x19,0x83,0x43,0x9a,0xfe,0xfc,0xc6,0xef,0x1e,0xa6,0x5b,0xbd,0xdd,0x08,
0xb9,0x16,0x9d,0xf2,0x20,0x78,0x06,0xae,0xf0,0xc3,0x59,0x45,0xb7,0xcd,0xa8,0x05,
0xeb,0x26,0x02,0xf6,0x6e,0x8a,0x14,0x98,0x0b,0x25,0xdc,0x7c,0xc1,0x6a,0x86,0x71,
0xed,0x68,0xec,0x85,0x35,0x87,0x18,0x93,0x46,0x23,0xd8,0x6c,0xbf,0x61,0x4a,0xac,
0xb4,0x8f,0xc7,0x36,0xe5,0x41,0xde,0xee,0xa7,0x9f,0xce,0xd5,0x21,0x8e,0xa5,0x32,
0x2b,0x54,0x99,0x7a,0xe0,0x58,0xfa,0x30,0x5d,0x44,0x8b,0x8c,0x4e,0x72,0x66,0xc4,
0x40,0x2f,0xe6,0x5a,0x2d,0x37,0x01,0x27,0xa2,0xe1,0x22,0x10,0xd7,0x3e,0x89,0x74,
0xc0,0x38,0x88,0xd9,0x62,0x49,0xdf,0x3b,0xd4,0x0a,0x13,0xf4,0xe4,0x07,0xe2,0xc9
};

template <
std::uint32_t m,
std::size_t rs1,std::size_t rs2,std::size_t rs3,std::size_t rs4,
std::size_t rm1,std::size_t rm2,
std::size_t rf
>
std::uint32_t big_wheel_engine<m,rs1,rs2,rs3,rs4,rm1,rm2,rf>::sbox_i5[256] = {
0x28,0x21,0xde,0xd7,0x8f,0x86,0xba,0xb3,0xeb,0xe2,0x7c,0x75,0x2d,0x24,0x18,0x11,
0x49,0x40,0x46,0x4f,0x17,0x1e,0x22,0x2b,0x73,0x7a,0xe4,0xed,0xb5,0xbc,0x80,0x89,
0xd1,0xd8,0x4b,0x42,0x1a,0x13,0x2f,0x26,0x7e,0x77,0xe9,0xe0,0xb8,0xb1,0x8d,0x84,
0xdc,0xd5,0xd3,0xda,0x82,0x8b,0xb7,0xbe,0xe6,0xef,0x71,0x78,0x20,0x29,0x15,0x1c,
0x0f,0x06,0x00,0x09,0x51,0x58,0x64,0x6d,0x35,0x3c,0xa2,0xab,0xf3,0xfa,0xc6,0xcf,
0x97,0x9e,0x0d,0x04,0x5c,0x55,0x69,0x60,0x38,0x31,0xaf,0xa6,0xfe,0xf7,0xcb,0xc2,
0x9a,0x93,0x95,0x9c,0xc4,0xcd,0xf1,0xf8,0xa0,0xa9,0x37,0x3e,0x66,0x6f,0x53,0x5a,
0x02,0x0b,0xf4,0xfd,0xa5,0xac,0x90,0x99,0xc1,0xc8,0x56,0x5f,0x07,0x0e,0x32,0x3b,
0x63,0x6a,0x6c,0x65,0x3d,0x34,0x08,0x01,0x59,0x50,0xce,0xc7,0x9f,0x96,0xaa,0xa3,
0xfb,0xf2,0x61,0x68,0x30,0x39,0x05,0x0c,0x54,0x5d,0xc3,0xca,0x92,0x9b,0xa7,0xae,
0xf6,0xff,0xf9,0xf0,0xa8,0xa1,0x9d,0x94,0xcc,0xc5,0x5b,0x52,0x0a,0x03,0x3f,0x36,
0x6e,0x67,0xb2,0xbb,0xe3,0xea,0xd6,0xdf,0x87,0x8e,0x10,0x19,0x41,0x48,0x74,0x7d,
0x25,0x2c,0x2a,0x23,0x7b,0x72,0x4e,0x47,0x1f,0x16,0x88,0x81,0xd9,0xd0,0xec,0xe5,
0xbd,0xb4,0x27,0x2e,0x76,0x7f,0x43,0x4a,0x12,0x1b,0x85,0x8c,0xd4,0xdd,0xe1,0xe8,
0xb0,0xb9,0xbf,0xb6,0xee,0xe7,0xdb,0xd2,0x8a,0x83,0x1d,0x14,0x4c,0x45,0x79,0x70
};

template <
std::uint32_t m,
std::size_t rs1,std::size_t rs2,std::size_t rs3,std::size_t rs4,
std::size_t rm1,std::size_t rm2,
std::size_t rf
>
std::uint32_t big_wheel_engine<m,rs1,rs2,rs3,rs4,rm1,rm2,rf>::sbox_r[256] = {
0xe9,0xe0,0xdc,0xd5,0x8d,0x84,0x1a,0x13,0x4b,0x42,0x7e,0x77,0x2f,0x26,0x20,0x29,
0x71,0x78,0x44,0x4d,0x15,0x1c,0x82,0x8b,0xd3,0xda,0xe6,0xef,0xb7,0xbe,0x6b,0x62,
0xa2,0xab,0x97,0x9e,0xc6,0xcf,0x51,0x58,0x00,0x09,0x35,0x3c,0x64,0x6d,0xfe,0xf7,
0xaf,0xa6,0x9a,0x93,0xcb,0xc2,0x5c,0x55,0x0d,0x04,0x38,0x31,0x69,0x60,0x66,0x6f,
0x37,0x3e,0x02,0x0b,0x53,0x5a,0xc4,0xcd,0x95,0x9c,0xa0,0xa9,0xf1,0xf8,0x07,0x0e,
0x56,0x5f,0x63,0x6a,0x32,0x3b,0xa5,0xac,0xf4,0xfd,0xc1,0xc8,0x90,0x99,0x9f,0x96,
0xce,0xc7,0xfb,0xf2,0xaa,0xa3,0x3d,0x34,0x6c,0x65,0x59,0x50,0x08,0x01,0x92,0x9b,
0xc3,0xca,0xf6,0xff,0xa7,0xae,0x30,0x39,0x61,0x68,0x54,0x5d,0x05,0x0c,0x0a,0x03,
0x5b,0x52,0x6e,0x67,0x3f,0x36,0xa8,0xa1,0xf9,0xf0,0xcc,0xc5,0x9d,0x94,0x41,0x48,
0x10,0x19,0x25,0x2c,0x74,0x7d,0xe3,0xea,0xb2,0xbb,0x87,0x8e,0xd6,0xdf,0xd9,0xd0,
0x88,0x81,0xbd,0xb4,0xec,0xe5,0x7b,0x72,0x2a,0x23,0x1f,0x16,0x4e,0x47,0xd4,0xdd,
0x85,0x8c,0xb0,0xb9,0xe1,0xe8,0x76,0x7f,0x27,0x2e,0x12,0x1b,0x43,0x4a,0x4c,0x45,
0x1d,0x14,0x28,0x21,0x79,0x70,0xee,0xe7,0xbf,0xb6,0x8a,0x83,0xdb,0xd2,0x2d,0x24,
0x7c,0x75,0x49,0x40,0x18,0x11,0x8f,0x86,0xde,0xd7,0xeb,0xe2,0xba,0xb3,0xb5,0xbc,
0xe4,0xed,0xd1,0xd8,0x80,0x89,0x17,0x1e,0x46,0x4f,0x73,0x7a,0x22,0x2b,0xb8,0xb1
};

#endif

So, I have some specific questions about it:

• How can I improve the performance of the program?
• How can the use of multithreading be improved?

Parallelism:
There are two major problems with your code

1. Your code doesn't run in parallel at all! Locking the mutex at the beginning of calculate_pi means you are always only run one instance of calculate_pi at the same time.
2. You are passing a copy of the same random number engine to each thread. As this is only a pseudo rando number generator, all threads will operate on the same sequence of numbers (you are doing the same work multiple times).

Solution:

1. The only thing you share among the codes is pcnt just define it as an atomic (std::atomic<std::uint32_t>). Even Better: Don't share anything: Let each task use it's own hit counter and sum them up at the end (easiest way to do this would be to use std::async).
2. Create and seed a separate random number engine for each thread

On my 4 Core machine with VS2015U3, this reduced the execution time from 120 to 20 Seconds and reduced the error roughly by a factor of two.

Small style tips:

1. Don't define variables before you use them. e.g.:
• gen_pi should be defined (as a constant) just when you need to store the result of approximate_pi.
• p1..9 and tmp_x should be (constant) function local variables
2. Const correctness: A lot of variables are only initialized and never changed afterwards - make them const
3. range based for: You could use more range based for loops (e.g. when joining the threads, creating the seed, or some loops that go over the state array of the random number engine.
• The idea of it was using the same PRNG for all the threads, so, is it a good idea to pass the generator by reference and put the lock guard before calling it? I mean, std::thread(calculate_pi<big_wheel_engine<>>,num_points_per_thread,std::ref(rnd)) for initialization and std::uint32_t rndx; { std::lock_guard<std::mutex> guard(mtx); rndx = rgen(); } ... for calling. Dec 26, 2016 at 19:24
• @alvaroceballos1: The idea of it was using the same PRNG for all the threads. That is most likely a very, very bad Idea. If you want to parallelize your work efficiently YOU MUST NOT SHARE MUTABLE STATE. I haven't put your code through a profiler, but most likely, generating the random number is a costly operation. Serializing that work and adding the overhead of locking/unlocking a contested mutex will likely leave you with no better performance than running your program on a single core. But feel free to try (and report the results). Dec 26, 2016 at 20:18
• Now I used std::async with separated random number engines as you suggested, and it works fine. Thanks! Dec 29, 2016 at 6:46