Estimation of pi using Monte Carlo method in fortran

Question

This is a a typical code to estimate pi (3.14...) value using Monte Carlo method. So I have trouble with the number of iterations of the do-while loop. Until the number of iterations less than or equal to 10000000, the approximate value of pi is correct, but when the number of iterations are more than that, the value of pi is wrong. These are the outputs for two different number of iterations.

Output 1. (for 10000000 iterations)
Approximate value of pi is 3.14104080
Number of points in circle 7852602.00
Number of points in square 10000000.0

Output 2. (for 100000000 iterations)
Approximate value of pi is 4.00000000
Number of points in circle 16777216.0
Number of points in square 16777216.0

Fortran Code:

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

A hint - which power of 2 is 16777216? How does that compare to the number of fraction bits in a single-precision real?

The answer is that they're the same - 24. Once you get to 16777216.0 adding 1 to it does not change the value as that's the largest integer that can be represented adding 1 in IEEE single precision.

The solution is to use double precision for your accumulations. You may want to also use double precision for your random number, since, for the same reason, there's an upper limit to how many different values can be returned in single precision.

Edit: I feel the need to expand on my answer. An IEEE single indeed has 24 bits of fraction, but the representation is chosen such that the most significant fraction bit is always 1, so it is implied or "hidden" and 23 bits are in the binary fraction field.

16777216 is indeed the largest integer with an epsilon (distance between representable values) of 1 in IEEE single (hex 4B800000) - the next larger representable integer is 16777218 (hex 4B800001). When you add 1 to 16777216 to get 16777217, that isn't representable and the rules call for "round to even", hence one gets 16777216 again.