!###########################################################
! page.31 表2.1
! ラグランジュ補間のε1からε4までを計算するプログラム!
! 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

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

! Nx 次元ラグランジュ多項式f_Nx(x2)を計算する
    DO Nx=1,4
	fflag=0.d0
    DO k=0,N-Nx,Nx
    DO ii=k*10,(k+Nx)*10
    fflag(ii)=lagrange(x2(ii))
    END DO
    END DO

eps=SUM((fftrue-fflag)**2)*h
WRITE(*,"('E',I1,'=',1PE15.8)") Nx,eps
    END DO

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