!###########################################################
! page.31 表2.1
! 不等間隔ラグランジュ補間のプログラム(εNを計算する)
!! f(x)=1/(1+25*x**2)
!###########################################################
PROGRAM prog04
IMPLICIT REAL(8)(a-h,o-z)
IMPLICIT INTEGER(i-n)

INTEGER,PARAMETER :: NMAX=48
REAL(8),DIMENSION(0:NMAX) :: x,y
REAL(8),DIMENSION(0:NMAX*10) :: fflag,fftrue,x2

pi=4.d0*datan(1.d0)

WRITE(*,*) '補間点の数N+1に対応するNとして(表2.1のNとして)Nを入力してください'
READ(*,*) N
WRITE(*,"('N=',I2)") N

fflag=0.d0;fftrue=0.d0
x=0.d0;y=0.d0;x2=0.d0

  dx=2.d0/dble(N)
  h=dx/10.d0

  DO i=0,N
  x(i)=-1.d0+dble(i)*dx
  y(i)=func(x(i))
  END DO

  DO i=0,N*10
  x2(i)=-1.d0+dble(i)*h
  fftrue(i)=func(x2(i))
  END DO
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
! N 次元ラグランジュ多項式f_N(x)を計算する
! 補間点は　x_i=cos(pi*(N-i)/N)
k=0
Nx=N
DO i=0,N
x(i)=dcos(pi*dble(N-i)/dble(N))
y(i)=func(x(i))
END DO
  
    fflag=0.d0
    DO ii=0,N*10
    fflag(ii)=lagrange(x2(ii))
    END DO

eps=SUM((fftrue-fflag)**2)*h
WRITE(*,"('EN=',1PE15.8)") eps

STOP
!###########################################################
CONTAINS
! Nx次元ラグランジュ多項式f_Nx(xx)を返す
! k番目の点からNx+1個の点を用いて補間を行う
FUNCTION lagrange(xx)
   REAL(8),INTENT(IN) :: xx
   REAL(8) :: lagrange
   REAL(8) :: bunbo,bunshi
   lagrange=0.d0
   DO i=k+0,k+Nx
     bunshi=1.d0
     bunbo=1.d0
        DO j=k+0,k+Nx
        IF(i==j)CYCLE
        bunshi=bunshi*(xx-x(j))
        bunbo=bunbo*(x(i)-x(j))
        END DO
     lagrange=lagrange+(bunshi/bunbo)*y(i)
   END DO
RETURN
END FUNCTION lagrange
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
FUNCTION func(xx)
   REAL(8),INTENT(IN) :: xx
   REAL(8) :: func
   func=1.d0/(1.d0+25.d0*xx**2)
RETURN
END FUNCTION func

END PROGRAM prog04