How to approximate pi out of random numbers
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
And it’s pretty easy to see, that
So we need to take two steps:
- Check if two random numbers are coprime
- 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.
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
I present you both algorithms in three ways:
- The basic code [ pi_calc_X.X.rb ]
- Does the calculation 400 times and returns the most accurate one [ pi_calc_gen_X.X.rb ]
- Gives out a remaining time estimate and calculates performance [ pi_long_X.X.rb ]