I know, I waste far to much of my time on youtube watching random videos. But out of luck I sometimes find inspiring stuff. Like this video by standupmaths:

And so I got on the journey to waist even more of my time.

How does it work?

The system of the approximation works pretty simple:

The probability, that two random numbers are coprime (Greatest common divisor is 1) to each other, is

The proof.

And it’s pretty easy to see, that

So we need to take two steps:

1. Check if two random numbers are coprime
2. Generate many large numbers and check them

The possibility x is based on random numbers between 1 and inf. If we limit the maximum random number we decrease the value of pi. So we need to maximize the max random number, what increases the calculation time. So we need a fast coprime algorithm.

The code

My first try is based on the euclidean algorithm.

def coprime?(x,y)  # Finds out if two numbers are coprime (Euclidean algorithm)
  if y == 0 then
    return x == 1
  else
    return coprime?(y, x % y)
  end
end

def pi(time,height)  # Approximate Pi out of random numbers
  truestate = 0.0
  time.times do
      if coprime?(rand(1..height),rand(1..height)) then truestate += 1 end
  end
  return Math.sqrt(6 /(truestate / time))
end

t = Time.new
p pi(100000000,1000000)
p Time.new - t

Real Value:         3. 141592653589793
Calculated Value:   3. 141484581244406
def. Deviation:    -0. 00010807234538701138
rel.  Deviation:    0. 0034400495959707733%
Runtime:           65. 7191823 Seconds

On my PC the calculation time of this is around 66 seconds.
The value for pi changes in each runthrough, but its still pretty accurate.

The optimized version

Thanks to reddit I found a few changes, that improve the runtime by cutting in half:

def pi(time,height)  # Calculates Pi out of random numbers
  truestate = 0
  range = 1..height
  i = 0
  while i < time
    truestate += 1 if rand(range).gcd(rand(range)) == 1
    i += 1
  end
  Math.sqrt(6 / (truestate.to_f / time))
end

t = Time.new
p pi(100000000,1000000)
p Time.new - t

Real Value:         3. 141592653589793
Calculated Value:   3. 141516024003613
def. Deviation:    -7. 662958618004367e-05
rel.  Deviation:    0. 002439195485526291%
Runtime:           28. 2471332 Seconds

Lets look at the changes:

• With using a while loop instead of a for loop we reduce the runtime by a bit
• Using the truestate iterator as an integer we reduce the internal effort
• With introducing a range variable we reduce the runtime by a lot
• With the internal gcd() function we reduce the runtime by a lot, even more than any of the above together

The programs

I present you both algorithms in three ways:

1. The basic code [ pi_calc_X.X.rb ]
2. Does the calculation 400 times and returns the most accurate one [ pi_calc_gen_X.X.rb ]
3. Gives out a remaining time estimate and calculates performance [ pi_long_X.X.rb ]