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

Question

这是使用蒙特卡罗方法估计 pi ​​(3.14...) 值的典型代码。所以我在 do-while 循环的迭代次数上遇到了麻烦。直到迭代次数小于或等于 10000000 时，pi 的近似值是正确的，但是当迭代次数大于此值时，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。