Monte Carlo Pi Approximations: C++, Julia and Python

My coding background began with C++ in middle school competitive programming before shifting to Python for astrophysics and cosmology projects due to its rich ecosystem (numpy, scipy, matplotlib, astropy). However, with the growing adoption of Julia in scientific computing and a desire to refresh my C++ knowledge, I set out to test how these three languages compare on a classic numerical benchmark: estimating $\pi$ using a Monte Carlo algorithm.

Monte Carlo algorithms are ideal for performance comparisons because they are straightforward to implement, CPU-bound, and scale effortlessly from millions to billions of iterations.

Setup

Monte Carlo methods estimate deterministic quantities using random sampling. Consider a unit square with an inscribed quadrant (or full unit circle centered at the origin). By generating uniform random coordinates $(x,y)$ within the unit square $[0,1] \times [0,1]$, the probability of a point falling inside the circular arc $x^2 + y^2 \leq 1$ is equal to the ratio of their areas:

$$A_\circ = \frac{1}{4} \pi r^2 = \frac{\pi}{4}, \quad A_\square = r^2 = 1$$ $$\implies \pi = 4 \frac{A_\circ}{A_\square} \approx 4 \frac{N_\text{inside}}{N_\text{total}}$$

By counting how many points satisfy $x^2 + y^2 \leq 1$, we can approximate $\pi$ with increasing precision as $N_\text{total}$ increases.

Implementation

To ensure a rigorous comparison, each script uses high-quality random number generator and language optimisations:

Language	Python 3.12	C++ (g++ 13.3)	Julia 1.12
Random Generator	numpy.random.random	std::mt19937 (<random>)	rand()
Execution Paradigm	Vectorised Batching	Native Scalar loop with -O3	JIT Compilation Loop
Time Measurement	time	high_resolution_clock	@btime (BenchmarkTools)

Optimisation Strategies

1. C++: Standard compiled scalar loops execute directly in CPU registers with small memory usage. Compiling with g++ -O3 enables maximum compiler optimisation and instruction-level parallelism.
2. Julia: Julia relies on Just-In-Time (JIT) compilation. A warm-up run executes before timing to ensure LLVM compilation overhead is excluded from the benchmark.
3. Python (NumPy): Pure Python loops suffer from interpreter overhead. Switching to NumPy vectorisation pushes the loops into C-level routines. To avoid allocating gigabytes of RAM when scaling up to $N = 10^{10}$ points, computation is executed in memory-efficient batches of $10^7$ iterations.

All implementations and benchmarking automation scripts are available in the repository.

Results

Performance was evaluated across scale sizes ranging from $N=10^6$ to $N=10^{10}$ total points.

Execution Speed

The first observation is that the plot is approximately linear, which means $\log t_\text{exec} = \log N$. Though it is noteworthy, that for increasing $N_\text{total}$ the absolute difference between languages becomes a 2x, 4x factor which should be considered for large runs.

Unexpectedly, Julia achieved the fastest execution time, even against C++. Julia’s JIT compiler optimised scalar loops into native assembly appears to be slightly more efficient for this use case than the closely matched C++ compiled code.

Using NumPy vectorisation reduced Python’s execution times from hundreds of seconds down to competitive numbers compared with previous script drafts. This demonstrates how array batching overcomes interpreter overhead in Python’s case.

Accuracy and Convergence

According to the Central Limit Theorem, the error in Monte Carlo integration scales as $O(1/\sqrt{N})$. All three languages demonstrated identical convergence rates, with the absolute error $|\hat{\pi} - \pi|$ decreasing reliably as sample size $N$ approached $10^{10}$ points. Small variance between runs at identical $N$ values is purely a result of differing random number engine seeds.

Conclusions

While traditional Python interpreter loops are poorly suited for raw numerical loops, modern NumPy vectorisation makes Python competitive for array-based operations on relatively small $N_\text{total} < 10^7$.

For maximum performance without writing C++ boilerplate, Julia provides an impressive balance: its JIT compiler generates machine code on par with optimised C++ while retaining the readable syntax of a high-level dynamically typed language.

This is my first attempt to a blog-style post on my website’s new ideas directory. If you have feedback, suggestions or any constructive comments, I encourage you to contact me. Thank you in advance!