在 fortran 中使用蒙特卡罗方法估计 pi

Question

这是使用蒙特卡洛方法估计 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 不会改变值，因为这是 IEEE 单精度加 1 可以表示的最大整数。

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

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

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