!###########################################################
! page.31 表2.1
! ラグランジュ補間とスプライン補間の誤差のプログラム
!
! 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
REAL(8),DIMENSION(4,NMAX+1) :: C

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

DO N=12,48,12
WRITE(*,*)
WRITE(*,"('N=',I2)") N

fflag=0.d0;fftrue=0.d0
x=0.d0;y=0.d0;x2=0.d0
C=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
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
! 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
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
! スプライン補間多項式P_i(x)を計算する
DO i=0,N
x(i)=-1.d0+dble(i)*dx
y(i)=func(x(i))
END DO
fflag=0.d0

CALL make_spline(N,x,y,C,NMAX)

fflag(0)=spline(x2(0),1)
DO i=1,N
DO i2=1,10
ii=(i-1)*10+i2
fflag(ii)=spline(x2(ii),i)
END DO
END DO

eps=SUM((fftrue-fflag)**2)*h
WRITE(*,"('ES=',1PE15.8)") 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
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
FUNCTION spline(xx,i)
   INTEGER,INTENT(IN) :: i
   REAL(8),INTENT(IN) :: xx
   REAL(8) :: spline,xd
   xd=xx-x(i-1)
   spline=C(1,i)+C(2,i)*xd+C(3,i)*xd**2+C(4,i)*xd**3
RETURN
END FUNCTION spline
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  

END PROGRAM prog04
!###########################################################

SUBROUTINE make_spline(N,x,y,C,NMAX)
IMPLICIT REAL(8)(a-h,o-z)
IMPLICIT INTEGER(i-n)

INTEGER,INTENT(IN) :: N,NMAX
REAL(8),INTENT(IN),DIMENSION(0:NMAX) :: x,y
REAL(8),INTENT(OUT),DIMENSION(4,NMAX+1) :: C

REAL(8),DIMENSION(NMAX+1) :: al,be,ga,yy
REAL(8),DIMENSION(NMAX) :: delx,dely

al=0.d0;be=0.d0;ga=0.d0;yy=0.d0
delx=0.d0;dely=0.d0;C=0.d0

DO i=1,N
delx(i)=x(i)-x(i-1)
dely(i)=y(i)-y(i-1)
END DO

al(1)=2.d0*delx(1)
al(N+1)=2.d0*delx(N)
ga(1)=delx(1)
be(N+1)=delx(N)
yy(1)=3.d0*dely(1)
yy(N+1)=3.d0*dely(N)

DO i=2,N
al(i)=2.d0*(delx(i)+delx(i-1))
be(i)=delx(i)
ga(i)=delx(i-1)
yy(i)=3.d0*(delx(i-1)*dely(i)/delx(i)+delx(i)*dely(i-1)/delx(i-1))
END DO

DO i=1,N
al(i+1)=al(i+1)-ga(i)*be(i+1)/al(i)
yy(i+1)=yy(i+1)-yy(i)*be(i+1)/al(i)
END DO

DO i=N,1,-1
yy(i)=yy(i)-yy(i+1)*ga(i)/al(i+1)
END DO
DO i=1,N+1
C(2,i)=yy(i)/al(i)
END DO

DO i=1,N
C(1,i)=y(i-1)
C(3,i)=(-2.d0*C(2,i)-C(2,i+1)+3.d0*dely(i)/delx(i))/delx(i)
C(4,i)=(C(2,i)+C(2,i+1)-2.d0*dely(i)/delx(i))/(delx(i)**2)
END DO

RETURN
END SUBROUTINE make_spline

!###########################################################