!###########################################################
! page.24
! ラグランジュ補間のプログラム
!###########################################################
PROGRAM prog03
IMPLICIT REAL(8)(a-h,o-z)
IMPLICIT INTEGER(i-n)

INTEGER,PARAMETER :: ii=4
REAL(8) :: x,y
REAL(8),DIMENSION(ii) :: fx0,x0

x=20.d0

fx0=0.d0;S=0.d0

fx0(1)=2.712
fx0(2)=2.714
fx0(3)=2.716
fx0(4)=2.718

x0=fx0**3

y=lagrange(x)

WRITE(*,'("f(",F8.3,")=",F15.12)') x,y

STOP

!###########################################################
CONTAINS

FUNCTION lagrange(x)
REAL(8) :: lagrange
REAL(8) :: bunbo,bunshi
REAL(8),INTENT(IN) :: x

lagrange=0.d0
DO i=1,ii
bunshi=1.d0
bunbo=1.d0
DO j=1,ii
IF(i==j)CYCLE
bunshi=bunshi*(x-x0(j))
bunbo=bunbo*(x0(i)-x0(j))
END DO
lagrange=lagrange+(bunshi/bunbo)*fx0(i)
END DO
RETURN
END FUNCTION lagrange

END PROGRAM prog03
!###########################################################