Floating Point Accuracy Problems While Calculating Pi

I'm trying to use random numbers to estimate the value of Pi.
The example I started with obviously has some problems which is a good thing as it means I have to understand it so I can fix it.

# Expected Result 3.141592653589793238
#
from math import sqrt
from random import random
# number of fandom points
N=100000
# number of points inside
I=0
for i in range(N):
  #print("I="+str(I))
  #Generate random point in 1x1 square
  x=random()
  y=random()
  # Is the point inside the circle?
  # r=sqrt(x**2 + y**2)
  #print(str(r))
  #if r<1:
  If (x*x + y*y) <1:    # Update 2
    I+=1
# print("Pi=" + str(4*I/N))
print("Pi=" + str(4*I)/N)) # Update 3

Expected Result

3.141592653589793238

Actual Results

• 100,000 Iterations

    Pi=3.14736
    Pi=3.14448
    Pi=3.14424
  

• 1,000,000 Iterations

    Pi=3.141496
    Pi=3.141356
    Pi=3.138
  

I would have expected the results to be more consistent.
How can I add more precision to the calculations?

UPDATE <br>
I modified the code to remove sqrt and ** as recommended in the comments. I cranked up the iterations to 4,000,000 but it didn’t make a lot of difference. It is still only accurate to about 2 decimal places. How can I add more precision to the calculation. (That question was in the title but it was edited out)

2,000,000 iterations

Pi=3.141342
Pi=3.141328
Pi=3.143074
Pi=3.139084

4,000,000 iterations

Pi=3.142605
Pi=3.141509
Pi=3.140663
Pi=3.143194

UPDATE 2 <br>
I changed the final calculation to the following.

import decimal

result=decimal.Decimal(4*I/N)
print("Pi=" + str(result))

4,000,000 iterations with decimal.Decimal()

Pi=3.14080899999999996197175278211943805217742919921875
Pi=3.143496999999999985675458447076380252838134765625
Pi=3.14080899999999996197175278211943805217742919921875
Pi=3.141859999999999875086587053374387323856353759765625

40,000,000 iterations with

decimal.Decimal()
Pi=3.141386900000000093058361017028801143169403076171875
Pi=3.1414208999999999605279299430549144744873046875
Pi=3.1414591999999998961357050575315952301025390625
Pi=3.14168000000000002813749233609996736049652099609375

This latest test looks like I have achieved accuracy of 3 decimal places but the results still look odd with big rows of 99999 or 00000. Any other ideas what is going on here?

Last a couple of larger itterations but I think the increase in accuracy is diminishing

100,000,000 Iterations (yes 100 million)

Pi=3.141542439999999825062104719108901917934417724609375
Pi=3.141720879999999826992507223621942102909088134765625
Pi=3.141750519999999990972128216526471078395843505859375
Pi=3.14174343999999994281324688927270472049713134765625
Pi=3.14132400000000000517275111633352935314178466796875

1,000,000,000 Iterations (1 billion too a while on my old laptop :))

Pi=3.141603959999999862162667341181077063083648681640625

Update 3 <BR>
I moved the brackets in final calculation as described by @Michael Butscher however even with 1 billion iterations the I only consistently obtained 3 decimal places accurately. Maybe I have hit some other limitation perhaps the sudo randomness of Python random numbers? I'll call this exercise done and incorporate my findings into my actual project.
OK I couldn't resist I did 10 billion as well
Surprisingly it didn't make much difference

1,000,000 Iterations

Pi=3.142056
Pi=3.136428
Pi=3.1407
Pi=3.141612

10,000,000 Iterations

Pi=3.142806
Pi=3.142266
Pi=3.141996
Pi=3.1422232

100,000,000 Iterations

Pi=3.14151576
Pi=3.1417604
Pi=3.1415038
Pi=3.1413738

1,000,000,000

Pi=3.141553108
Pi=3.1415895
Pi=3.141629112

10,000,000,000 Iterations (10 Billion)
Took about 12 hours to execute

Pi=3.1416011832

The final division is currently done as

result=decimal.Decimal(4*I/N)

but this means that the division is using float arithmetic which produces imprecise results like 3.141603959999999862162667341181077063083648681640625 instead of 3.14160396.

Moreover the division can produce nonsense due to an overflow for very large integers to divide.

Instead of that, the int can be converted exactly to a Decimal first and then be divided by another int:

result=decimal.Decimal(4*I)/N

This time the division uses the implementation for Decimal which is exact for the defined precision (28 digits by default).

Here is my final solution all in one place.
I was hoping for a better result. Perhaps this needs more iterations.

from math import sqrt
from random import random
from decimal import Decimal

# number of random points
N=10000000
# number of points inside
I=0
for i in range(N):
  #Generate random point in 1x1 square
  x=Decimal(str(random()))
  y=Decimal(str(random()))
  # Is the point inside the circle?
  if (x*x + y*y) < Decimal("1"):    
    I+=1
result=Decimal(4*I)/N
print( "Expected Result= 3.141592653589793238 ")
print("          Result="+str(result))

# Exp    3.141592653589793238
# Result=3.141676  Iterations=10,000,000
# Result=3.1418548 Iterations=10,000,000