Different speed when estimating pi

General Usage
1

I used two sets of codes to estimate the value of pi as a simple example of Monte Carlo methods.

Code 1:


function compute_pi(N::Int)
  
  n_landed_in_circle = 0  

  for i = 1:N
      x = rand()   
      y = rand()   

      r2 = x*x + y*y  
      if r2 < 1.0
          n_landed_in_circle += 1
      end
  end

  return n_landed_in_circle / N * 4.0    
end

Code2:

function compute_pi1(N::Int) 
  return 4 * sum([1 for _ in filter(_ -> rand()^2 + rand()^2 < 1, 1:N)]) / N 
end

I tested the speed with N = 1_000_000_000, and found that code1 is 2 times faster than code 2. I am not sure what makes the difference.

I used Julia 6.2 on Fedora 25 OS with I7 core and 16gb ram.

1 Like
2

The second one allocates an array of ones. You may want to try a generator.

4 Likes
3 
julia> function sum1() s=0; for i in 1:3 s+=1 end;s;end;
julia> function sum2() sum(1 for i in 1:3); end;

julia> @code_llvm sum1()
julia> @code_llvm sum2()

Is optimization in sum1 Julia compiler’s work or llvm’s?

4

ifelse should speed this up around 40% (on my pc)

  function compute_pi_ifelse(N::Int)
  
    n_landed_in_circle = 0  
  
    for i = 1:N
        x = rand()   
        y = rand()   
  
        r2 = x*x + y*y  
        n_landed_in_circle += ifelse(r2<1.0,1,0)        
    end
  
    return n_landed_in_circle / N * 4.0    
  end
1 Like
5

Also, if you can (and want to) use memory, you may be able to speed things up further by generating multiple random numbers at once.
The last function below is a little faster on my system (you will need to tune the second argument; here 10’000 seems to work reasonably well).


using Compat, StatsBase,StaticArrays
using Revise
using BenchmarkTools

function compute_pi(N::Int)
  
    n_landed_in_circle = 0  
  
    for i = 1:N
        x = rand()   
        y = rand()   
  
        r2 = x*x + y*y  
        if r2 < 1.0
            n_landed_in_circle += 1
        end
    end
  
    return n_landed_in_circle / N * 4.0    
  end


  function compute_pi_ifelse(N::Int)
  
    n_landed_in_circle = 0  
  
    for i = 1:N
        x = rand()   
        y = rand()   
  
        r2 = x*x + y*y  
        n_landed_in_circle += ifelse(r2<1.0,1,0)        
    end
  
    return n_landed_in_circle / N * 4.0    
  end


  function compute_pi_vector(N::Int,vec_size::Int)
    n_landed_in_circle = 0 
    randvec = zeros(vec_size) 
     
    counter = vec_size
    for i = 1:N
        if counter >= vec_size-1
            rand!(randvec)
            counter = 1
        end
        @inbounds x = randvec[counter]           
        @inbounds y = randvec[counter+1]   
        counter += 2   
  
        r2 = x*x + y*y  
        n_landed_in_circle += ifelse(r2<1.0,1,0)                
    end
  
    return n_landed_in_circle / N * 4.0   
  end

  
  @benchmark compute_pi_ifelse(1_000_000)

  @benchmark compute_pi_vector(1_000_000,100)
  @benchmark compute_pi_vector(1_000_000,1000)
  @benchmark compute_pi_vector(1_000_000,10000)
2 Likes