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

    call getarg(1,dummy)
    read(dummy,*) n
    allocate(tt(1:n),ww(1:n))
    call GIP(n,tt,ww)
    do i=1,n
        write(6,'("gu(",i2,")=",f20.16,"D0 ; ","gw(",i2,")=",f20.16,"D0")') i,tt(i),i,ww(i)
    end do

stop

contains

    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_GAUSSIP