C ==========================================================
C
C  page.53 表3.2
C
C   ガウス・ルジャンドル積分公式のプログラム
C
C   f(x)=1.d0/sqrt(1.d0-x**2)
C
C ==========================================================
      PROGRAM prog08
      IMPLICIT REAL*8 (a-h,o-z)
      PARAMETER(NMAX=200)
      REAL*8 xx2(NMAX),ww2(NMAX)
      REAL*8 P(0:NMAX)
      REAL*8 I3
      pi=4.d0*datan(1.d0)
      DO 10 N=50,200,50
        xx2=0.d0
        ww2=0.d0
C
C ======================================================
C     ガウス・ルジャンドル公式
C ======================================================
        CALL ZEROTEN(N,NMAX,xx2)
C
        DO 20 i=1,N
          a=0.d0
          CALL legendre(N,xx2(i),P,NMAX)
          DO 30 k=0,N-1
            a=a+dble(2*k+1)*(P(k))**2
   30   CONTINUE
        ww2(i)=2.d0/a
   20   CONTINUE
C
        I3=0.d0
        DO 40 i=1,N
          I3=I3+ww2(i)*func(xx2(i))
   40   CONTINUE
        E3=dabs(I3-pi)
        WRITE(*,"('N=',I3,1p3E16.8)") N,E3
   10 CONTINUE
      STOP
      END
C
C ======================================================
C
C        ルジャンドル多項式
C       P_j(x) (j=0,1,2,...,N)
C
C ======================================================
      SUBROUTINE legendre(N,x,P,NMAX)
      IMPLICIT REAL*8 (a-h,o-z)
      REAL*8 P(0:NMAX)
      DO 10 I=0,NMAX
        P(I)=0.d0
   10 CONTINUE
      P(0)=1.d0
      P(1)=x
      DO 20 j=1,N-1
        P(j+1)=(dble(2*j+1)*x*P(j)-dble(j)*P(j-1))/dble(j+1)
   20 CONTINUE
      RETURN
      END
C
C ======================================================
C
      SUBROUTINE ZEROTEN(N,NMAX,xx)
      IMPLICIT REAL*8 (a-h,o-z)
      REAL*8 xx(NMAX)
      REAL*8 P(0:NMAX)
C
      pi=4.d0*datan(1.d0)
      DO 40 I=1,NMAX
        xx(I)=0.d0
   40 CONTINUE
      delta=1.d-10
      eps=1.d-10
      M=N/2
C
      DO 10 i=1,M
        a=dsin(pi*dble(N-2*i)/dble(2*N))
        b=dsin(pi*dble(N-2*i+2)/dble(2*N))
        CALL legendre(N,a,P,NMAX)
        fa=P(N)
        CALL legendre(N,b,P,NMAX)
        fb=P(N)
        IF(fa*fb>0) THEN
          WRITE(*,*) 'Error a,b,  i==',i
        ENDIF
C
        fc=1.d0
        DO 20 J=1,100000
          IF(dabs(fc)>delta.and.dabs(a-b)>eps) THEN
            diff=dabs(a-b)
            c=(a+b)/2.d0
            CALL legendre(N,c,P,NMAX)
            fc=P(N)
C
            IF(fc*fa>0.d0)THEN
              a=c
              fa=fc
            ELSE
              b=c
              fb=fc
            END IF
          ELSE
            GO TO 100
          ENDIF
   20   CONTINUE
C
  100 CONTINUE
      xx(N-i+1)=c
      xx(i)=-c
C
   10 CONTINUE
      RETURN
      END
C
C ======================================================
C
      FUNCTION func(x)
      REAL*8 x,func
      func=1.d0/dsqrt(1.d0-x**2)
      RETURN
      END
C ======================================================