用随机数估算圆周率的浮点精度问题及优化探究

我尝试用随机数估算圆周率π的值，参考的示例代码存在问题，修复部分问题后仍达不到理想精度。初始代码如下：

# Expected Result 3.141592653589793238
#
from math import sqrt
from random import random
# number of random 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

预期结果为3.141592653589793238，但实际测试结果一致性很差：

• 100,000次迭代结果：

  Pi=3.14736
  Pi=3.14448
  Pi=3.14424
  

• 1,000,000次迭代结果：

  Pi=3.141496
  Pi=3.141356
  Pi=3.138
  

我希望能提升计算精度。

UPDATE
我按建议移除了sqrt和**运算，将迭代次数提升至4,000,000次，但精度仍仅约2位小数：

• 2,000,000次迭代结果：

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

• 4,000,000次迭代结果：

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

UPDATE 2
我引入decimal模块优化计算：

import decimal

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

4,000,000次迭代结果虽能达到3位小数精度，但出现大量连续9或0的情况；即使迭代到100,000,000次、1,000,000,000次，精度提升依旧有限。

UPDATE 3
我修正了计算语句的括号位置（修正后为print("Pi=" + str(4*I/N))），测试1,000,000次、10,000,000次、100,000,000次、1,000,000,000次甚至10,000,000,000次迭代后，结果仍仅能稳定到3位小数。我怀疑这可能受Python伪随机数的特性限制，请问如何进一步提升估算的精度？

内容的提问来源于stack exchange，提问作者David P