Different ways to get pi in FORTRAN (apart from 4*atan(1.0))

leipniz.f90

use,intrinsic :: iso_fortran_env, only : output_unit
implicit none
intrinsic :: kind,command_argument_count,get_command,index,sum,abs,epsilon
integer,parameter :: sp = kind(1.0e0)
integer,parameter :: dp = kind(1.0d0)
integer,parameter :: qp = kind(1.0q0)
integer(kind=8)   :: i
integer(kind=8)   :: i_max
character(len=:),allocatable :: line
character(len=2)  :: ch
logical  :: inv
integer  :: argc,err
real(sp) :: spi,stmp
real(dp) :: dpi,dtmp
real(qp) :: qpi,qtmp

argc = command_argument_count()
if (argc.lt.2) error stop 'Usage: leipniz <sp|dp|qp> <max_iter|eps>'

call get_command(length=i,status=err)
if (err.ne.0) error stop 'could not read command line'
allocate(character(len=i) :: line)
call get_command(line)
i = index(line,' ')
line = line(i:)

if (index(line,'eps')) then
   inv = .false.
   read(line,*,iostat=err) ch
   if (err.ne.0) error stop 'could not read cmd arguments'
else
read(line,*,iostat=err) ch,i_max
   inv = .true.
   if (err.ne.0) error stop 'could not read cmd arguments'
   if (i_max.le.0) error stop 'max_iter must be greater then zero'
endif

select case(ch)
case('sp')
   if (inv) then
      spi =  sum( (/ (l_sp(i), i=i_max, 0, -2) /) )
   else
      i = 0
      spi = 0.0_sp
      do
         stmp = l_sp(i)
         i = i+2
         if (abs(stmp).lt.epsilon(spi)) exit
         spi = spi + stmp
      enddo
   endif
   write(output_unit,'(''Single:'',x)',advance='no')
   write(output_unit,*) 4*spi
case('dp')
   if (inv) then
      dpi = sum( (/ (l_dp(i), i=i_max, 0, -2) /) )
   else
      i = 0
      dpi = 0.0_dp
      do
         dtmp = l_dp(i)
         i = i+2
         if (abs(dtmp).lt.epsilon(dpi)) exit
         dpi = dpi + dtmp
      enddo
   endif
   write(output_unit,'(''Double:'',x)',advance='no')
   write(output_unit,*) 4*dpi
case('qp')
   if (inv) then
      qpi = sum( (/ (l_qp(i), i=i_max, 0, -2) /) )
   else
      i = 0
      qpi = 0.0_qp
      do
         qtmp = l_qp(i)
         i = i+2
         if (abs(qtmp).lt.epsilon(qpi)) exit
         qpi = qpi + qtmp
      enddo
   endif
   write(output_unit,'(''Quadr.:'',x)',advance='no')
   write(output_unit,*) 4*qpi
case default
   error stop 'please specify precision as sp, dp or qp'
end select
write(output_unit,'(3x''Ref:    3.141592653589793238462643383279502884197'')')

contains
elemental function l_sp(i) result(res)
   implicit none
   integer,parameter :: sp = kind(1.0e0)
   integer(kind=8),intent(in) :: i
   real(sp) :: res
   res = 1._sp/(2._sp*i+1._sp)-1._sp/(2._sp*(i+1)+1._sp)
end function l_sp
elemental function l_dp(i) result(res)
   implicit none
   integer,parameter :: dp = kind(1.0d0)
   integer(kind=8),intent(in) :: i
   real(dp) :: res
   res = 1._dp/(2._dp*i+1._dp)-1._dp/(2._dp*(i+1)+1._dp)
end function l_dp
elemental function l_qp(i) result(res)
   implicit none
   integer,parameter :: qp = kind(1.0q0)
   integer(kind=8),intent(in) :: i
   real(qp) :: res
   res = 1._qp/(2._qp*i+1._qp)-1._qp/(2._qp*(i+1)+1._qp)
end function l_qp

end

mcint.f90

use,intrinsic :: iso_fortran_env, only : output_unit
implicit none
intrinsic :: kind,command_argument_count,get_command
intrinsic :: index,random_seed,random_number,sum,real
integer,parameter :: sp = kind(1.0e0)
integer,parameter :: dp = kind(1.0d0)
integer,parameter :: qp = kind(1.0q0)
integer(kind=8)   :: i,j
integer(kind=8)   :: i_max
character(len=:),allocatable :: line
character(len=2)  :: ch
integer  :: argc,err
real(sp) :: sdum(2),spi
real(dp) :: ddum(2),dpi
real(qp) :: qdum(2),qpi

argc = command_argument_count()
if (argc.lt.2) error stop 'Usage: mcint <sp|dp|qp> <samples>'

call get_command(length=i,status=err)
if (err.ne.0) error stop 'could not read command line'
allocate(character(len=i) :: line)
call get_command(line)
i = index(line,' ')
line = line(i:)

read(line,*,iostat=err) ch,i_max
if (err.ne.0) error stop 'could not read cmd arguments'
if (i_max.le.0) error stop 'max_iter must be greater then zero'

call random_seed

j=0
select case(ch)
! 1 = x**2 + y**2 <=> y = sqrt( 1 - x**2 )
case('sp')
   do i = 1, i_max
      call random_number(sdum)
      if (sum(sdum**2).lt.1.0_sp) j = j+1
   enddo
   write(output_unit,'(''Single:'',x)',advance='no')
   write(output_unit,*) 4*(real(j,sp)/real(i_max,sp))
case('dp')
   do i = 1, i_max
      call random_number(ddum)
      if (sum(ddum**2).lt.1.0_dp) j = j+1
   enddo
   write(output_unit,'(''Double:'',x)',advance='no')
   write(output_unit,*) 4*(real(j,dp)/real(i_max,dp))
case('qp')
   do i = 1, i_max
      call random_number(qdum)
      if (sum(qdum**2).lt.1.0_qp) j = j+1
   enddo
   write(output_unit,'(''Quadr.:'',x)',advance='no')
   write(output_unit,*) 4*(real(j,qp)/real(i_max,qp))
case default
   error stop 'please specify precision as sp, dp or qp'
end select

write(output_unit,'(3x''Ref:    3.141592653589793238462643383279502884197'')')

end

trapzd.f90

use,intrinsic :: iso_fortran_env, only : output_unit
implicit none
intrinsic :: kind,command_argument_count,get_command,index,sum,sqrt,real
integer,parameter :: sp = kind(1.0e0)
integer,parameter :: dp = kind(1.0d0)
integer,parameter :: qp = kind(1.0q0)
integer(kind=8)   :: i
integer(kind=8)   :: i_max
character(len=:),allocatable :: line
character(len=2)  :: ch
integer  :: argc,err
real(sp) :: spi,stab
real(dp) :: dpi,dtab
real(qp) :: qpi,qtab

argc = command_argument_count()
if (argc.lt.2) error stop 'Usage: trapzd <sp|dp|qp> <tabl>'

call get_command(length=i,status=err)
if (err.ne.0) error stop 'could not read command line'
allocate(character(len=i) :: line)
call get_command(line)
i = index(line,' ')
line = line(i:)

read(line,*,iostat=err) ch,i_max
if (err.ne.0) error stop 'could not read cmd arguments'
if (i_max.le.0) error stop 'max_iter must be greater then zero'

select case(ch)
! 1 = x**2 + y**2 <=> y = sqrt( 1 - x**2 )
case('sp')
   stab = 1.0_sp/real(i_max,sp)
   spi=(sum((/(sqrt(1.0_sp-(stab*i)**2),i=i_max,1,-1)/))+0.5_sp)*stab
   write(output_unit,'(''Single:'',x)',advance='no')
   write(output_unit,*) 4*spi
case('dp')
   dtab = 1.0_dp/real(i_max,dp)
   dpi=(sum((/(sqrt(1.0_dp-(dtab*i)**2),i=i_max,1,-1)/))+0.5_dp)*dtab
   write(output_unit,'(''Double:'',x)',advance='no')
   write(output_unit,*) 4*dpi
case('qp')
   qtab = 1.0_qp/real(i_max,qp)
   qpi=(sum((/(sqrt(1.0_qp-(qtab*i)**2),i=i_max,1,-1)/))+0.5_qp)*qtab
   write(output_unit,'(''Quadr.:'',x)',advance='no')
   write(output_unit,*) 4*qpi
case default
   error stop 'please specify precision as sp, dp or qp'
end select

write(output_unit,'(3x''Ref:    3.141592653589793238462643383279502884197'')')

end