The problem with π (iii) – Pascal to Fortran

By: spqr

13 January 2023 at 15:31

A while back I posted on the Spigot algorithm in Pascal (and C) for calculating the first 1000 decimal places of π. Now converting it to Fortran is quite easy.

program piSpigot

   integer :: n, len
   integer :: i, j, k, q, x, nines, predigit
   integer, dimension(3500) :: a

   n = 1000
   len = (10 * n) / 3

   a = 2
   nines = 0
   predigit = 0
   999 format(I1)

   do j = 1,n
      q = 0
      do i = len,1,-1
         x = 10 * a(i) + q * i
         a(i) = mod(x,2*i-1)
         q = x / (2*i-1)
      end do
      a(1) = mod(q,10)
      q = q / 10
      if (q == 9) then
         nines = nines + 1
      elseif (q == 10) then
         write(*,999,advance='no') predigit+1
         do k = 1, nines
            write(*,999,advance='no') 0
         end do
         predigit = 0
         nines = 0
      else
         write(*,999,advance='no') predigit
         predigit = q
         if (nines /= 0) then
            do k = 1,nines
               write(*,999,advance='no') 9
               nines = 0
            end do
         end if
      end if
   end do
   write(*,999,advance='no') predigit

end program piSpigot

This algorithm is good, but it could be improved upon by making the array dynamic. All that really involves is reworking some code at the top of the program. In the snippet of code below, the array a is made allocatable (Line 5), I/O is added to prompt for the number of decimal points required to be calculated, and then on Line 11 the space for a is allocated. At the end of the program, the memory is deallocated.

program piSpigot

   integer :: n, len
   integer :: i, j, k, q, x, nines, predigit
   integer, dimension(:), allocatable :: a

   write(*,*) 'Number of decimal points? '
   read(*,*) n
   len = (10 * n) / 3

   allocate(a(len))
   ...
   deallocate(a)
end program piSpigot

Normal view