subroutine gauss_legendre(pi)
  implicit none
  real*16,intent(out) :: pi
  integer,parameter :: Nmax=10
  real*16,dimension(0:Nmax) :: a,b,t,p
  integer :: i
  
  a(0)=1q0;b(0)=1q0/(2q0**0.5q0);t(0)=0.25q0;p(0)=1q0
  do i=1,Nmax
     a(i)=(a(i-1)+b(i-1))*0.5q0
     b(i)=sqrt(a(i-1)*b(i-1))
     t(i)=t(i-1)-p(i-1)*(a(i-1)-a(i))**2
     p(i)=2q0*p(i-1)
     pi=0.25q0*(a(i)+b(i))**2/t(i)
!     write(6,'(I4,10(1x,es32.24))') i,pi,2q0*acos(0q0),pi-2q0*acos(0q0)
  end do

  return
end subroutine gauss_legendre
!=================================================
subroutine piarctan(nmax,pi)
  implicit none
  integer,intent(in) :: nmax
  real*16,intent(out) :: pi
  integer :: k,n
  real*16,dimension(0:1000) :: sumv
  integer,dimension(0:1000) :: sumv2

  pi=2q0
  do k=1,nmax
     pi=pi+cf(k)*(0.5q0)**(k-1)
  end do

  return
contains
  recursive real*16 function cf(k) result(res)
    implicit none
    integer,intent(in) :: k

    if(k>=2) then
       res = cf(k-1)*(2q0*(k)/(2q0*(k)+1q0))
    elseif(k==1) then
       res=2q0/3q0
    end if
    
  end function cf
end subroutine piarctan
!=================================================
program main
  implicit none
  integer :: n
  real*16 :: pi,arctan
  real*16,external :: fact

  call gauss_legendre(pi)
  write(6,'(1x,f34.32,1x,a)') pi,"by Gauss-Legendre formula."
  call piarctan(110,pi)
  write(6,'(1x,f34.32,1x,a)') pi,"by Euler's expansion formula of arctan."
  
end program main