algorithm - 在 fortran 中使用蒙特卡罗方法估计 pi
问题描述
这是使用蒙特卡罗方法估计 pi (3.14...) 值的典型代码。所以我对 do-while 循环的迭代次数有疑问。直到迭代次数小于等于10000000,pi的近似值是正确的,但是当迭代次数多于10000000时,pi的值是错误的。这些是两个不同迭代次数的输出。
输出 1.(10000000 次迭代)
pi 的近似值为 3.14104080
圆圈中的点数 7852602.00
正方形中的点数 10000000.0
输出 2.(对于 100000000 次迭代)
pi 的近似值为 4.00000000
圆圈
中的点数 16777216.0 正方形中的点数 16777216.0
Fortran 代码:
program estimate_pi
implicit none
real :: length, r, rand_x, rand_y, radius, area_square, square_points, circle_points, origin_dist, approx_pi
integer :: i, iterations
! length of side of a square
length = 1
! radius of the circle
radius = length/2
! area of the square considered
area_square = length * length
square_points = 0
circle_points = 0
! number of iterations
iterations = 10000000
i = 0
do while (i < iterations)
! generate random x and y values
call random_number(r)
rand_x = - length/2 + length * r
call random_number(r)
rand_y = - length/2 + length * r
! calculate the distance of the point from the origin
origin_dist = rand_x * rand_x + rand_y * rand_y
! check whether the point is within the circle
if (origin_dist <= radius * radius) then
circle_points = circle_points + 1
end if
! total number of points generated
square_points = square_points + 1
! approximate value of pi
approx_pi = 4 * (circle_points/square_points)
! increase the counter by +1
i = i + 1
end do
print*, 'Approximate value of pi is', approx_pi
print*, 'Number of points in circle', circle_points
print*, 'Number of points in square', square_points
end program estimate_pi
解决方案
提示 - 2 的哪个幂是 16777216?这与单精度实数中的小数位数相比如何?
答案是它们是相同的 - 24。一旦你到达 16777216.0,加 1 不会改变值,因为这是可以表示的最大整数,在 IEEE 单精度中加 1。
解决方案是对累积使用双精度。您可能还希望对随机数使用双精度,因为出于同样的原因,单精度可以返回多少个不同的值也有上限。
编辑:我觉得有必要扩展我的答案。一个 IEEE 单曲确实有 24 位小数,但选择的表示形式是最高有效小数位始终为 1,因此它是隐含的或“隐藏的”,并且 23 位在二进制小数字段中。
16777216 确实是 IEEE 单一(十六进制 4B800000)中 epsilon(可表示值之间的距离)为 1 的最大整数 - 下一个更大的可表示整数是 16777218(十六进制 4B800001)。当您将 1 加到 16777216 以得到 16777217 时,这是不可表示的,并且规则要求“四舍五入”,因此再次得到 16777216。
推荐阅读
- python - 将 dcc.Slider 连接到 Dash 中的 clientside_callback 以获取动画图
- sql - 如何让 Rails 5 与带有句点的主键很好地配合?
- selenium - 如何在节点中使用 selenium 设置 scraperapi 代理
- javascript - Javascript - 如何使用变量作为键检查数组的值
- google-tag-manager - 有没有办法跟踪在谷歌标签管理器上填写的任何表格?
- javascript - 如何在映射后映射另一个数组对象?
- php - 如何根据选择选项中的主键从数据库中获取某些数据?
- java - 如何将 LinkedList 写入 .txt 文件
- .net-core - Hangfire (.net core) - 有没有办法计算一个作业或重复作业执行了多少次?
- python - 如何在 Python 中遍历 Outlook 文件夹