module para
   integer,parameter::n=12
end module para

program gauss_legendre
   use para
   implicit none
   integer i,m,k,kk
   double precision a,b,lambdak,wi,yi,ii
   double precision,dimension(0:n)::xi
   double precision pi,hpi
   double precision f,x,p,c

   f(x)=dsqrt(1.0d0-x*x)

   pi=dacos(-1.0d0)
   hpi=pi/2
   if(mod(n,2)==0)then
      m=n/2
   else
      m=(n-1)/2
   end if
   do i=1,m
      a=dsin((pi*dble(n-2*i))/dble(2*n))
      b=dsin((pi*dble(n+2-2*i))/dble(2*n))
      call nibun(a,b,c)
      xi(i)=c
      xi(i+m)=-xi(i)
      write(*,100) i,xi(i)
   end do
   100 format(1x,'x(',i2,')=',1x,f10.7)
   if(mod(n,2)==1) xi(n)=0.0d0
   ii=0.0d0
   do i=1,n
      wi=0.0d0
      do k=0,n-1
         lambdak=2.0d0/(2.0d0*dble(k)+1.0d0)
         call p_k(k,xi(i),p)
         wi=wi+p**2/lambdak
      end do
      wi=1.0d0/wi
      yi=f(xi(i))
      ii=ii+wi*yi
   end do
   write(*,200) n,dabs(ii-hpi)
   200 format(1x,'N=',1x,i4,1x,'|I-Pi/2|=',1x,f16.10)
   stop
   end

subroutine p_k(k,x,p)
   implicit none
   integer k
   double precision x,p
   integer n
   double precision p0,p1,p2
   p0=1.0d0
   p1=x
   if(k==0)then
      p=p0
   else
      do n=1,k-1
         p2=((2.0d0*dble(n)+1.0d0)*x*p1-dble(n)*p0)/dble(n+1)
         p0=p1
         p1=p2
      end do
      p=p1
   end if
   return
   end

subroutine nibun(a,b,c)
   use para
   implicit none
   double precision a,b,c,p,p1,p2
   integer k
   integer,parameter::kmax=1000
   double precision,parameter::delta=10.0d-10,epsilon=10.0d-10
   do k=1,kmax
      c=(a+b)/2.0d0
      call p_k(n,c,p)
      if(dabs(p)<delta.or.dabs(a-b)<epsilon)then
         exit
      else
         call p_k(n,c,p1)
         call p_k(n,a,p2)
         if(p1*p2<0.0d0)then
            b=c
         else
            a=c
         end if
      end if
   end do
   return
   end