Can't accurately calculate pi on Python

Question

I am new member here and I'm gonna drive straight into this as I've spent my whole Sunday trying to get my head around it.

I'm new to Python, having previously learned coding on C++ to a basic-intermediate level (it was a 10-week university module).

I'm trying a couple of iterative techniques to calculate Pi but both are coming up slightly inaccurate and I'm not sure why.

The first method I was taught at university - I'm sure some of you have seen it done before.

x=0.0
y=0.0
incircle = 0.0
outcircle = 0.0
pi = 0.0
i = 0
while (i<100000):
    x = random.uniform(-1,1)
    y = random.uniform(-1,1)
    if (x*x+y*y<=1):
        incircle=incircle+1
    else:
        outcircle=outcircle+1
    i=i+1
pi = (incircle/outcircle)
print pi

It's essentially a generator for random (x,y) co-ordinates on a plane from -1 to +1 on both axes. Then if x^2+y^2 <= 1, we know the point rests inside a circle of radius 1 within the box formed by the co-ordinate axes.

Depending on the position of the point, a counter increases for incircle or outcircle.

The value for pi is then the ratio of values inside and outside the circle. The co-ordinates are randomly generated so it should be an even spread.

However, even at very high iteration values, my result for Pi is always around the 3.65 mark.

The second method is another iteration which calculates the circumference of a polygon with increasing number of sides until the polygon is almost a circle, then, Pi=Circumference/diameter. (I sort of cheated because the coding has a math.cos(Pi) term so it looks like I'm using Pi to find Pi, but this is only because you can't easily use degrees to represent angles on Python). But even for high iterations the final result seems to end around 3.20, which again is wrong. The code is here:

S = 0.0
C = 0.0
L = 1.0

n = 2.0
k = 3.0
while (n<2000):
    S = 2.0**k
    L = L/(2.0*math.cos((math.pi)/(4.0*n)))
    C = S*L
    n=n+2.0
    k=k+1.0

pi = C/math.sqrt(2.0)
print pi

I remember, when doing my C++ course, being told that the problem is a common one and it isn't due to the maths but because of something within the coding, however I can't remember exactly. It may be to do with the random number generation, or the limitations of using floating point numbers, or... anything really. It could even just be my mathematics...

Can anyone think what the issue is?

TL;DR: Trying to calculate Pi, I can get close to it but never very accurately, no matter how many iterations I do.

(Oh and another point - in the second code there's a line saying S=2.0**k. If I set 'n' to anything higher than 2000, the value for S becomes too big to handle and the code crashes. How can I fix this?)

Thanks!

This is a math problem. The method of Monte-Carlo gives an approximation of pi, not pi itself. [This](http://rosettacode.org/wiki/Pi#Python) should be more accurate. — Rolbrok, Jan 17 '16 at 16:48
I've also noticed that python is sometimes a little off on its calculations. For example, when applying `tan(45)` degrees it returns 0.99999... instead of 1. — Ashwin Gupta, Jan 17 '16 at 16:53
@AshwinGupta This isn't a shortcoming of just Python, but any language that implements floating point arithmetic. Also, it's tan(45) that is equal to 1. — Reti43, Jan 17 '16 at 16:56
By the way, Python internally uses float64, which have a maximum exponent of +307. This means you can't represent 2.0**1024. — Reti43, Jan 17 '16 at 17:12

Mike Müller · Accepted Answer · 2016-01-17T17:40:01.470

7

The algorithm for your first version should look more like this:

from __future__ import division, print_function

import sys
if sys.version_info.major < 3:
    range = xrange

import random 


incircle = 0
n = 100000
for n in range(n):
    x = random.random()
    y = random.random()
    if (x*x + y*y <= 1):
        incircle += 1
pi = (incircle / n) * 4
print(pi)

Prints:

3.14699146991

This is closer. Increase n to get even closer to pi.

The algorithm takes into account only one quarter of the unit circle, i.e. with a radius of 1.

The formula for the area of a quarter circle is:

area_c = (pi * r **2) / 4

That for the area of the square containing this circle:

area_s = r **2

where r is the radius of the circle.

Now the ratio is:

area_c / area_s

substitute the equations above, re-arange and you get:

pi = 4 * (area_c / area_s)

Going Monte Carlo, just replace both areas by a very high number that represents them. Typically, the analogy of darts thrown randomly is used here.

edited Jan 17 '16 at 17:40

answered Jan 17 '16 at 16:50

Mike Müller

82,630
20
166
161

Damn, was just posting that :). Oh well, not really a python guy anyway. – Mike Wise Jan 17 '16 at 16:51
It is (four times) the ratio of the points in the middle to all of the points, just just the points on the outside... You should mention that. – Mike Wise Jan 17 '16 at 16:52
Ahhh right, thanks! I get it now. I admit I was just going over the coding from memory (more than a year ago) and, rather stupidly without stopping to think about the basic maths behind it, I seemed to think it was simply (incircle/outcircle) but now, having rummaged through some old files and found the original C++ code, it does indeed do use the method explained here. Thanks again. – Stuart Aitken Jan 17 '16 at 23:55

score 2 · Answer 2 · answered Jan 17 '16 at 17:24

For the first one, your calculation should be

pi = incircle/1000000*4  # 3.145376..

This is the number of points that landed inside of the circle over the number of total points (approximately 0.785671 on my run).

With a radius of 1 (random.uniform(-1,1)), the total area is 4, so if you multiple 4 by the ratio of points that landed inside of the circle, you get the correct answer.

Can't accurately calculate pi on Python

2 Answers2