在 Fortran 中使用 Monte Carlo 方法估算 pi

Question

这是使用Monte Carlo 方法估算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

Answer 1

提示 - 2 的哪个次方是 16777216？这与单精度实数中的小数位数相比如何？

答案是它们是相同的 - 24。一旦你得到 16777216.0 加 1 就不会改变值，因为这是可以表示的最大整数加 1 IEEE 单精度。

解决方案是对累加使用双精度。您可能还想对随机数使用双精度，因为出于同样的原因，可以以单精度返回多少个不同的值有上限。

编辑：我觉得有必要扩展我的答案。一个 IEEE 单数确实有 24 位的小数，但选择的表示形式使得最高有效小数位始终为 1，因此暗示或 "hidden" 并且 23 位在二进制小数字段中。

16777216 确实是 IEEE 单个（十六进制 4B800000）中 epsilon（可表示值之间的距离）为 1 的最大整数 - 下一个更大的可表示整数是 16777218（十六进制 4B800001）。当您将 1 添加到 16777216 得到 16777217 时，这是不可表示的并且规则要求 "round to even"，因此再次得到 16777216。

在 Fortran 中使用 Monte Carlo 方法估算 pi

Estimation of pi using Monte Carlo method in fortran

algorithm

fortran

pi

montecarlo