icp/oer/courses/c-openmp/sections/04-problem-solving/01-montecarlo/solution.c

52 lines
1.3 KiB
C

#include <stdio.h>
#include <omp.h>
#include <math.h>
#include <stdlib.h>
#include <stdbool.h>
#include <string.h>
static inline bool in_circle(double x, double y)
{
return (x*x + y*y < 1.0);
}
int main(int argc, const char *argv[])
{
omp_set_num_threads(10);
// A constant seed is used to keep our results reproducible
srand48(42);
double x;
double y;
int in_circle_count = 0;
int total_count = 0;
unsigned short rand_state[3];
double start = omp_get_wtime();
// We use now use the Monte Carlo method of approximating pi
// TODO: parallelize this while avoiding data races:
#pragma omp parallel private(x, y, rand_state) reduction(+:total_count, in_circle_count)
{
rand_state[0] = rand_state[0] * (unsigned short)omp_get_num_threads();
#pragma omp for
for (int i = 0; i < 90000000; i++) {
x = erand48(rand_state);
y = erand48(rand_state);
if (in_circle(x, y)) {
in_circle_count++;
}
total_count++;
}
}
double quater_pi = (double)in_circle_count / (double)total_count;
double end = omp_get_wtime();
printf("Pi is: %f\n", quater_pi * 4);
printf("That took a total of %f seconds.\n", end - start);
return 0;
}