Section 14: Pi Estimation

Section 14: Pi Estimation#

Adapted from: “Introduction to Programming Using Fortran 95/2003/2008” by Ed Jorgensen (March 2018 / Version 3.0.51).

Program to Estimate the Value of Pi#

program piestimation

! declare variabes

implicit none

integer                             :: count, alstat, i, incount
real                                :: x, y, pi_est, pt
real, allocatable, dimension (:,:)  :: points

! display inital header
    write (*, '(/a/)') "Program Example - Pi estimation"

! prompt for and obtain count value
    ! This section will be skipped for the Jupyter Notebook
    !do
        ! prompt for count value
        ! write (*, '(a)', advance="no") "Enter Count (100 - 1,000,000): "

        ! read count value
        ! read (*,*) count

        ! if count is correct, exit loop
        !if ( count >= 100 .and. count <= 1000000 ) exit

        ! Otherwise, display error message
        !write (*, '(a,a,/a)') "Error, count must be ", &
        !    "between 100 and 1,000,000.",              &
        !    "Please re-enter."

    !end do

! Set number of estimation points (i.e. count) = 1 million
count = 1000000

! allocate two dimensional array
    allocate (points(count,2), stat=alstat)

    ! Check for allocation errors
    if ( alstat /= 0 ) then
        write (*, '(a,a/a)') "Error, unable to",        &
            " allocate memory.", "Program terminated."
        stop
    end if


! generate points
    call random_seed()

    ! loop count times
    do i = 1, count

        ! generate x and y values
        call random_number(x)
        call random_number(y)

        ! place (x,y) values in array
        points(i,1) = x
        
        points(i,2) = y
    end do

! perform monte carlo estimation

    ! set count of samples inside circle = 0
    incount = 0

    ! loop count times
    do i = 1, count

    ! if [ sqrt (x(i)^2 + y(i)^2) < 1.0 ]
    ! increment count of samples inside circle

        pt = points(i,1)**2 + points(i,2)**2
        if (sqrt(pt) < 1.0) incount = incount + 1

    end do

    pi_est = 4.0 * real (incount) / real (count)

! display results
    write (*, '(a, i10)') "The number of points to process is -> ", count
    print *
    write (*, '(a, i10)') "The number of points inside the circle were -> ", incount
    print *
    write (*, '(a, f16.8)') "Estimated Pi Value: ", pi_est
    print *
end program piestimation

The following explanation is given by “Code Copilot” within ChatGPT 4.

The Fortran program piestimation is designed to estimate the value of \(\pi\) using a Monte Carlo method. This method is a statistical technique that uses random sampling to approximate mathematical constants or solve mathematical problems. Here’s a breakdown of the program’s structure and functionality:

Variable Declaration#

  • count: The number of random points to be generated for the estimation.

  • alstat: Status variable for array allocation.

  • i: Loop control variable.

  • incount: Counter for the number of points that fall inside the unit circle.

  • x, y: Coordinates of a randomly generated point.

  • pi_est: The estimated value of \(\pi\).

  • pt: Temporary variable to store the sum of squares of x and y.

  • points: A 2-dimensional allocatable array to store the x and y coordinates of each point.

User Interaction#

  • The program begins by printing a header message to indicate its function.

  • It initially contains commented-out code that prompts the user to enter a value for count with validations to ensure the input is within a specified range (100 to 1,000,000). This interactive portion is omitted in this version, likely to accommodate automated execution environments like Jupyter Notebooks.

Setting Parameters#

  • Directly sets count to 1,000,000 for simplicity, meaning one million points will be used for the estimation.

Memory Allocation#

  • Allocates space for the points array, which will hold the x and y coordinates of the randomly generated points. It checks for successful allocation and terminates the program if memory allocation fails.

Random Point Generation#

  • Initializes the random number generator.

  • Loops count times to generate random x and y coordinates for each point. These coordinates are stored in the points array.

Monte Carlo Estimation#

  • Resets the incount (inside count) to zero.

  • Loops through each point, calculating the distance from the origin using the Euclidean distance formula (i.e., checking if \(x^2 + y^2 < 1\)). If the condition is true, the point is inside the unit circle, and incount is incremented.

  • Calculates the estimated value of \(\pi\) using the formula:

\[ \Large \pi_{\text{est}} = 4 \times \frac{\text{incount}}{\text{count}} \]

This formula derives from the ratio of the area of the unit circle to the area of the bounding square (since the radius of the circle is 1, the area of the circle is \(\pi\), and the area of the square is \(2^2 = 4\)).

Display Results#

  • Prints the total number of points processed.

  • Displays the number of points that fell inside the unit circle.

  • Outputs the estimated value of \(\pi\).

This program is an excellent example of using random sampling methods for numerical approximation and demonstrates fundamental concepts of probabilistic simulation in Fortran.

Program Compilation and Execution#

The above program is compiled and run using Fortran Package Manager (fpm):

import os
root_dir = os.getcwd()

code_dir = root_dir + "/" + "Fortran_Code/Section_14_Pi_Estimation"

os.chdir(code_dir)

build_status = os.system("fpm build 2>/dev/null")

exec_status = os.system("fpm run 2>/dev/null")

Program Example - Pi estimation

The number of points to process is ->    1000000

The number of points inside the circle were ->     784823

Estimated Pi Value:       3.13929200