Montecarlo Pi

Python Numpy

from __future__ import print_function
from benchpress import util
import numpy as np

def montecarlo_pi_se(samples, B):

    return np.sum(np.sqrt(
        B.random_array((samples,), dtype=B.dtype)**2 + \
        B.random_array((samples,), dtype=B.dtype)**2
    ) <= 1.0) * 4.0 / samples

def montecarlo_pi(samples, B):

    x = B.random_array((samples,), dtype=B.dtype)
    y = B.random_array((samples,), dtype=B.dtype)
    m = np.sqrt(x*x + y*y) <= 1.0

    return np.sum(m) * 4.0 / samples

def solve(samples, iterations, B):
    acc=0.0
    for _ in range(iterations):
        acc += montecarlo_pi(samples, B)
	util.Benchmark().flush()
    acc /= iterations
    return acc

def main():
    B = util.Benchmark()
    samples, iterations = B.size
    B.start()
    R = solve(samples, iterations, B)
    B.stop()
    B.pprint()
    if B.verbose:
        print(R)

if __name__ == "__main__":
    main()

C99 Omp

#include <inttypes.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <Random123/philox.h>
#include <bp_util.h>

typedef union philox2x32_as_1x64 {
    philox2x32_ctr_t orig;
    uint64_t combined;
} philox2x32_as_1x64_t;

double montecarlo_pi(uint32_t samples, uint32_t x_index, uint32_t y_index, uint32_t key)
{
    uint64_t darts = 0;
    
    #pragma omp parallel for reduction(+:darts)
    for (uint32_t eidx=0; eidx<samples; ++eidx) {

        uint64_t x_raw = ((philox2x32_as_1x64_t)philox2x32(
          ((philox2x32_as_1x64_t)( (uint64_t)(x_index + eidx) )).orig,
          (philox2x32_key_t){ { key } }
        )).combined;
        double x = x_raw / 18446744073709551616.000000;

        uint64_t y_raw = ((philox2x32_as_1x64_t)philox2x32(
          ((philox2x32_as_1x64_t)( (uint64_t)(y_index + eidx) )).orig,
          (philox2x32_key_t){ { key } }
        )).combined;
        double y = y_raw / 18446744073709551616.000000;

        darts += sqrt(x*x + y*y) <= 1;
    }
    return (double)darts * (double)4.0 / (double)samples;
}

double run_montecarlo_pi(int64_t samples, int64_t iterations)
{
    uint64_t index = 0;                 // Philox index
    const uint64_t key = 0;             // Philox key
    uint64_t x_index = 0;
    uint64_t y_index = 0;

    double pi_accu = 0.0;               // Accumulation over PI-approximations.
    for(int64_t i=0; i<iterations; ++i) {
        x_index = index;
        index += samples;
        y_index = index;
        index += samples;
        pi_accu += montecarlo_pi(samples, x_index, y_index, key);
    }
    pi_accu /= iterations;              // Approximated value of PI

    return pi_accu;
}

int main(int argc, char** argv)
{
    bp_util_type bp = bp_util_create(argc, argv, 2);// Grab arguments
    if (bp.args.has_error) {
        return 1;
    }
    const int samples = bp.args.sizes[0];
    const int iterations = bp.args.sizes[1];

    bp.timer_start();
    double pi = run_montecarlo_pi(samples, iterations);
    bp.timer_stop();

    bp.print("montecarlo_pi(c99_omp)");
    if (bp.args.verbose) {
        printf("PI-approximation = %f\n", pi);
    }

    return 0;
}

C99 Seq

#include <inttypes.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <Random123/philox.h>
#include <bp_util.h>

typedef union philox2x32_as_1x64 {
    philox2x32_ctr_t orig;
    uint64_t combined;
} philox2x32_as_1x64_t;

double montecarlo_pi(uint32_t samples, uint32_t x_index, uint32_t y_index, uint32_t key)
{
    uint64_t darts = 0;
    
    for (uint32_t eidx=0; eidx<samples; ++eidx) {

        uint64_t x_raw = ((philox2x32_as_1x64_t)philox2x32(
          ((philox2x32_as_1x64_t)( (uint64_t)(x_index + eidx) )).orig,
          (philox2x32_key_t){ { key } }
        )).combined;
        double x = x_raw / 18446744073709551616.000000;

        uint64_t y_raw = ((philox2x32_as_1x64_t)philox2x32(
          ((philox2x32_as_1x64_t)( (uint64_t)(y_index + eidx) )).orig,
          (philox2x32_key_t){ { key } }
        )).combined;
        double y = y_raw / 18446744073709551616.000000;

        darts += sqrt(x*x + y*y) <= 1;
    }
    return (double)darts * (double)4.0 / (double)samples;
}

double run_montecarlo_pi(int64_t samples, int64_t iterations)
{
    uint64_t index = 0;                 // Philox index
    const uint64_t key = 0;             // Philox key
    uint64_t x_index = 0;
    uint64_t y_index = 0;

    double pi_accu = 0.0;               // Accumulation over PI-approximations.
    for(int64_t i=0; i<iterations; ++i) {
        x_index = index;
        index += samples;
        y_index = index;
        index += samples;
        pi_accu += montecarlo_pi(samples, x_index, y_index, key);
    }
    pi_accu /= iterations;              // Approximated value of PI

    return pi_accu;
}

int main(int argc, char** argv)
{
    bp_util_type bp = bp_util_create(argc, argv, 2);// Grab arguments
    if (bp.args.has_error) {
        return 1;
    }
    const int samples = bp.args.sizes[0];
    const int iterations = bp.args.sizes[1];

    bp.timer_start();
    double pi = run_montecarlo_pi(samples, iterations);
    bp.timer_stop();

    bp.print("montecarlo_pi(c99_seq)");
    if (bp.args.verbose) {
        printf("PI-approximation = %f\n", pi);
    }

    return 0;
}

Cpp11 Bxx

#include <iostream>
#include "bxx/bohrium.hpp"
#include <bp_util.h>

using namespace std;
using namespace bxx;

template <typename T>
multi_array<T>& monte_carlo_pi(int samples)
{
    multi_array<T> x, y, m;             // Operands

    x = randu<T>(samples);              // Sample random numbers
    y = randu<T>(samples);
    m = as<T>(sqrt(x*x + y*y) <= 1.0);  // Model

    return sum(m) * 4.0 / (T)samples;   // Count
}

template <typename T>
T solve(int samples, int iterations)
{
    multi_array<T> acc;                 // Acculumate across iterations
    acc = zeros<T>(1);
    for(int i=0; i<iterations; ++i) {
        acc += monte_carlo_pi<T>(samples);        
    }
    acc /= (T)iterations;
    return scalar(acc);
}

int main(int argc, char* argv[])
{
    bp_util_type bp = bp_util_create(argc, argv, 2);// Grab arguments
    if (bp.args.has_error) {
        return 1;
    }
    const int samples = bp.args.sizes[0];
    const int iterations = bp.args.sizes[1];

    Runtime::instance().flush();
    bp.timer_start();
    double pi = solve<double>(samples, iterations);
    Runtime::instance().flush();
    bp.timer_stop();
                                                    // Output timing
    bp.print("montecarlo_pi(cpp11_bxx)");
    if (bp.args.verbose) {                          // and values.
        cout << "PI-approximation = " << pi << endl;
    }

    return 0;
}

Cpp11 Omp

#include <inttypes.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <Random123/philox.h>
#include <bp_util.h>

typedef union philox2x32_as_1x64 {
    philox2x32_ctr_t orig;
    uint64_t combined;
} philox2x32_as_1x64_t;

double montecarlo_pi(uint32_t samples, uint32_t x_index, uint32_t y_index, uint32_t key)
{
    uint64_t darts = 0;
    
    #pragma omp parallel for reduction(+:darts)
    for (uint32_t eidx=0; eidx<samples; ++eidx) {

        philox2x32_as_1x64_t x_counter;
        x_counter.combined = x_index + eidx;

        philox2x32_as_1x64_t x_philox;
        x_philox.orig = philox2x32(
          x_counter.orig,
          (philox2x32_key_t){ { key } }
        );
        double x = x_philox.combined / 18446744073709551616.000000;

        philox2x32_as_1x64_t y_counter;
        y_counter.combined = y_index + eidx;
        
        philox2x32_as_1x64_t y_philox;
        y_philox.orig = philox2x32(
          y_counter.orig,
          (philox2x32_key_t){ { key } }
        );
        double y = y_philox.combined / 18446744073709551616.000000;

        darts += sqrt(x*x + y*y) <= 1;
    }
    return (double)darts * (double)4.0 / (double)samples;
}

double run_montecarlo_pi(int64_t samples, int64_t iterations)
{
    uint64_t index = 0;                 // Philox index
    const uint64_t key = 0;             // Philox key
    uint64_t x_index = 0;
    uint64_t y_index = 0;

    double pi_accu = 0.0;               // Accumulation over PI-approximations.
    for(int64_t i=0; i<iterations; ++i) {
        x_index = index;
        index += samples;
        y_index = index;
        index += samples;
        pi_accu += montecarlo_pi(samples, x_index, y_index, key);
    }
    pi_accu /= iterations;              // Approximated value of PI

    return pi_accu;
}

int main(int argc, char** argv)
{
    bp_util_type bp = bp_util_create(argc, argv, 2);// Grab arguments
    if (bp.args.has_error) {
        return 1;
    }
    const int samples = bp.args.sizes[0];
    const int iterations = bp.args.sizes[1];

    bp.timer_start();
    double pi = run_montecarlo_pi(samples, iterations);
    bp.timer_stop();

    bp.print("montecarlo_pi(cpp11_omp)");
    if (bp.args.verbose) {
        printf("PI-approximation = %f\n", pi);
    }

    return 0;
}

Cpp11 Seq

#include <inttypes.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <Random123/philox.h>
#include <bp_util.h>

typedef union philox2x32_as_1x64 {
    philox2x32_ctr_t orig;
    uint64_t combined;
} philox2x32_as_1x64_t;

double montecarlo_pi(uint32_t samples, uint32_t x_index, uint32_t y_index, uint32_t key)
{
    uint64_t darts = 0;
    
    for (uint32_t eidx=0; eidx<samples; ++eidx) {

        philox2x32_as_1x64_t x_counter;
        x_counter.combined = x_index + eidx;

        philox2x32_as_1x64_t x_philox;
        x_philox.orig = philox2x32(
          x_counter.orig,
          (philox2x32_key_t){ { key } }
        );
        double x = x_philox.combined / 18446744073709551616.000000;

        philox2x32_as_1x64_t y_counter;
        y_counter.combined = y_index + eidx;
        
        philox2x32_as_1x64_t y_philox;
        y_philox.orig = philox2x32(
          y_counter.orig,
          (philox2x32_key_t){ { key } }
        );
        double y = y_philox.combined / 18446744073709551616.000000;

        darts += sqrt(x*x + y*y) <= 1;
    }
    return (double)darts * (double)4.0 / (double)samples;
}

double run_montecarlo_pi(int64_t samples, int64_t iterations)
{
    uint64_t index = 0;                 // Philox index
    const uint64_t key = 0;             // Philox key
    uint64_t x_index = 0;
    uint64_t y_index = 0;

    double pi_accu = 0.0;               // Accumulation over PI-approximations.
    for(int64_t i=0; i<iterations; ++i) {
        x_index = index;
        index += samples;
        y_index = index;
        index += samples;
        pi_accu += montecarlo_pi(samples, x_index, y_index, key);
    }
    pi_accu /= iterations;              // Approximated value of PI

    return pi_accu;
}

int main(int argc, char** argv)
{
    bp_util_type bp = bp_util_create(argc, argv, 2);// Grab arguments
    if (bp.args.has_error) {
        return 1;
    }
    const int samples = bp.args.sizes[0];
    const int iterations = bp.args.sizes[1];

    bp.timer_start();
    double pi = run_montecarlo_pi(samples, iterations);
    bp.timer_stop();

    bp.print("montecarlo_pi(cpp11_seq)");
    if (bp.args.verbose) {
        printf("PI-approximation = %f\n", pi);
    }

    return 0;
}