Add Pi Monte Carlo Parallel Simulation

README.md

@@ -35,6 +35,7 @@ Examples include (and will expand to):
     * Encoding
         * [rot47](./rot47/)
     * [prefix-sum](./prefix-sum/)
+    * [pi-monte-carlo](./pi-monte-carlo/)
 * Data Structures
     * [map-with-unknown-key](./map-with-unknown-key/)
 * OOP

pi-monte-carlo/Makefile

@@ -0,0 +1,30 @@
+# pull in shared compiler settings
+include ../common.mk
+
+# per-example flags
+# CXXFLAGS += -pthread
+
+## get it from the folder name
+TARGET := $(notdir $(CURDIR))
+SRCS   := $(wildcard *.cpp)
+OBJS   := $(SRCS:.cpp=.o)
+
+all: $(TARGET)
+
+$(TARGET): $(OBJS)
+	$(CXX) $(CXXFLAGS) -o $@ $^
+
+%.o: %.cpp
+	$(CXX) $(CXXFLAGS) -c $< -o $@
+
+run: $(TARGET)
+	./$(TARGET) $(ARGS)
+
+clean:
+	rm -f $(OBJS) $(TARGET)
+
+# Delegates to top-level Makefile
+check-format:
+	$(MAKE) -f ../Makefile check-format DIR=$(CURDIR)
+
+.PHONY: all clean run check-format

pi-monte-carlo/main.cpp

@@ -0,0 +1,61 @@
+#include <algorithm>
+#include <cstdint>
+#include <iostream>
+#include <numeric>
+#include <random>
+#include <thread>
+#include <vector>
+
+static std::uint64_t
+count_hits(std::uint64_t n, std::uint64_t seed)
+{
+    std::mt19937_64 rng(seed);
+    std::uniform_real_distribution<double> dist(0.0, 1.0);
+
+    std::uint64_t hits = 0;
+    for (std::uint64_t i = 0; i < n; ++i) {
+        const double x = dist(rng);
+        const double y = dist(rng);
+        if (x * x + y * y <= 1.0)
+            ++hits;
+    }
+    return hits;
+}
+
+int
+main(int argc, char* argv[])
+{
+    int threads = 2;
+    std::uint64_t n = 1'000'000;
+
+    if (argc > 1)
+        threads = std::max(1, std::atoi(argv[1]));
+    if (argc > 2)
+        n = std::max<std::uint64_t>(1, std::strtoull(argv[2], nullptr, 10));
+
+    const std::uint64_t base = n / threads;
+    const std::uint64_t rem = n % threads;
+
+    std::vector<std::thread> pool;
+    std::vector<std::uint64_t> hits(threads, 0);
+
+    // a single nondeterministic seed source; each thread gets a different seed
+    std::random_device rd;
+    const std::uint64_t master_seed = (static_cast<std::uint64_t>(rd()) << 32) ^ static_cast<std::uint64_t>(rd());
+
+    for (int i = 0; i < threads; ++i) {
+        const std::uint64_t points = base + (static_cast<std::uint64_t>(i) < rem ? 1 : 0);
+        const std::uint64_t seed = master_seed + 0x9e3779b97f4a7c15ULL * static_cast<std::uint64_t>(i + 1);
+
+        pool.emplace_back([&, i, points, seed] { hits[i] = count_hits(points, seed); });
+    }
+
+    for (auto& t : pool)
+        t.join();
+
+    const std::uint64_t total_hits = std::accumulate(hits.begin(), hits.end(), std::uint64_t{0});
+    const long double pi = 4.0 * static_cast<long double>(total_hits) / static_cast<long double>(n);
+
+    std::cout << pi << "\n";
+    return 0;
+}

pi-monte-carlo/tests.sh

@@ -0,0 +1,15 @@
+#!/bin/bash
+
+set -ex
+
+# Run the tests
+pi=$(./pi-monte-carlo)
+
+# Check that the result is within a reasonable range
+if (( $(echo "$pi < 3.0" | bc -l) )) ||
+   (( $(echo "$pi > 3.2" | bc -l) )); then
+    echo "Test failed: pi is out of range: $pi"
+    exit 1
+else
+    echo "Test passed: pi is within range: $pi"
+fi