program f90_CALPI
    !Gauss積分点の座標と重み
    !こちらを参考にしました
    !http://www7.ocn.ne.jp/~kawa1/numeric.pdf
    implicit none
    integer::n,i
    real(8),allocatable::tt(:),ww(:)
    real(8)::a,b,x,s

    a=0.0D0
    b=1.0D0

    n=100
    allocate(tt(1:n),ww(1:n))
    call GIP(n,tt,ww)
    s=0.0D0
    do i=1,n
        x=0.5D0*(b-a)*tt(i)+0.5D0*(b+a)
        s=s+ww(i)*FX(x)
    end do
    s=0.5D0*(b-a)*s
    write(6,*) s,atan(1.0D0)*4.0D0

stop

contains

    real(8) function FX(x)
        real(8),intent(in)::x
        FX=4.0D0/(1.0D0+x*x)
    end function FX


    subroutine GIP(n,tt,ww)
        integer,intent(in)::n
        real(8),intent(out)::tt(1:n),ww(1:n)

        real(8),parameter::pi=atan(1.0D0)*4.0D0
        integer::m,mm,i,k,nn
        real(8),allocatable::t(:),w(:)
        real(8)::xx,x0,pn,dp
        real(8)::eps=1.0D-16

        m=(n+1)/2
        mm=n/2
        allocate(t(1:m),w(1:m))
        do i=1,m
            xx=cos(0.5D0*pi*dble(4*i-1)/dble(2*n+1))
            do
                x0=xx
                call PNDP(n,x0,pn,dp)
                xx=x0-pn/dp
                if(abs(xx-x0)<eps)exit
            end do
            call PNDP(n-1,xx,pn,dp)
            w(i)=2.0D0*(1.0D0-xx**2.0D0)/(dble(n)*pn)**2.0D0
            t(i)=xx
        end do

        nn=0
        do i=1,m-1
            nn=nn+1
            tt(nn)=-t(i)
            ww(nn)=w(i)
        end do
        i=m
        nn=nn+1
        if(m==mm)then
            tt(nn)=-t(i)
            ww(nn)=w(i)
        else
            tt(nn)=0.0D0
            ww(nn)=w(i)
        end if
        k=m
        if(m/=mm)k=m-1
        do i=k,1,-1
            nn=nn+1
            tt(nn)=t(i)
            ww(nn)=w(i)
        end do
    end subroutine GIP


    subroutine PNDP(n,x,pn,dp)
        integer,intent(in)::n
        real(8),intent(in)::x
        real(8),intent(out)::pn
        real(8),intent(out)::dp

        integer::i
        real(8)::p0,p1

        p1=1.0D0
        pn=x
        do i=2,n
            p0=p1
            p1=pn
            pn=((2.0D0*dble(i)-1.0D0)*x*p1-(dble(i)-1.0D0)*p0)/dble(i)
        end do
        dp=dble(n)*(x*pn-p1)/(x**2.0D0-1.0D0)
    end subroutine PNDP

end program f90_CALPI