# Speed up Julia code for simple Monte Carlo Pi estimation (compared to Numba)

**URL:** <https://discourse.julialang.org/t/speed-up-julia-code-for-simple-monte-carlo-pi-estimation-compared-to-numba/59808>\
**Category:** Performance\
**Tags:** performance\
**Created:** [April 22, 2021, 2:25pm UTC](https://discourse.julialang.org/t/speed-up-julia-code-for-simple-monte-carlo-pi-estimation-compared-to-numba/59808 "2021-04-22T14:25:16Z")\
**Posts on this page:** 1\
**Showing post:** 19

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [April 22, 2021, 4:23pm UTC](https://discourse.julialang.org/t/speed-up-julia-code-for-simple-monte-carlo-pi-estimation-compared-to-numba/59808/19 "2021-04-22T16:23:32Z")

</div>

Indeed, that speeds up by a factor of 2 both alternatives. Copying the solution from here: [Random number and parallel execution - #16 by greg\_plowman](https://discourse.julialang.org/t/random-number-and-parallel-execution/57024/16)

We have:

```julia
julia> Threads.nthreads()
4

julia> @btime estimate_pi_thread(10000000)
  35.675 ms (22 allocations: 1.98 KiB)
3.1413444

julia> @btime estimate_pi_floop(10000000)
  12.244 ms (57 allocations: 2.77 KiB)
3.1405184

```

But the code is quite boring for that:

> **Code**
>
> ```julia
> using Random, Future, BenchmarkTools, FLoops
> 
> function parallel_rngs(rng::MersenneTwister, n::Integer)
> step = big(10)^20
> rngs = Vector{MersenneTwister}(undef, n)
> rngs[1] = copy(rng)
> for i = 2:n
> rngs[i] = Future.randjump(rngs[i-1], step)
> end
> return rngs
> end
> 
> const N = Threads.nthreads() * 10^8
> const rng = MersenneTwister();
> const rngs = parallel_rngs(MersenneTwister(), Threads.nthreads());
> 
> function estimate_pi_thread(nMC)
> radius = 1.
> diameter = 2. * radius 
> n_circle = zeros(Int,Threads.nthreads())
> Threads.@threads for i in 1:nMC
> rng = rngs[Threads.threadid()]
> x = (rand(rng) - 0.5) * diameter
> y = (rand(rng) - 0.5) * diameter
> r = sqrt(x^2 + y^2) 
> if r <= radius 
> n_circle[Threads.threadid()] += 1
> end    
> end    
> return (sum(n_circle) / nMC) * 4.
> end 
> 
> function estimate_pi_floop(nMC)
> radius = 1.
> diameter = 2. * radius
> @floop for i in 1:nMC
> rng = rngs[Threads.threadid()]               
> x = (rand(rng) - 0.5) * diameter
> y = (rand(rng) - 0.5) * diameter
> r = sqrt(x^2 + y^2)
> if r <= radius
> @reduce(n_circle += 1)
> end
> end
> return (n_circle / nMC) * 4.
> end
> 
> ```

---

_[View the full topic](https://discourse.julialang.org/t/speed-up-julia-code-for-simple-monte-carlo-pi-estimation-compared-to-numba/59808)._