      PROGRAM INTERPOL
      implicit none
      real*8 X,F,P,XX,A,B,h,O,E
      integer j,k,N,M,pyt
      dimension X(0:100),XX(0:100),F(0:100),P(0:100)
      
      O=real(0)
      E=real(1)
      
   10 open(unit=10,file='X')
      open(unit=11,file='F')
      open(unit=12,file='P')
      open(unit=13,file='FO')
      open(unit=14,file='XX')
      write(*,*)'PROGRAM INTERPOL'
      write(*,*)'PODAJ N: LICZBA WEZLOW N+1'
      read (*,*) N
      write(*,*)'PODAWAJ X'
      do 1,k=0,N,1
      read (*,*)  X(k)
    1 continue
      write(*,*)'X - KONIEC'
      write(*,*)'PODAWAJ F'
      do 2,k=0,N,1
      read (*,*)  F(k)
      write(13,*) F(k)
    2 continue
      write(*,*)'F - KONIEC'
      write(*,*)'PODAJ PRZEDZIAL TABLICOWANIA [A,B]'
      read (*,*) A,B
      write(*,*)'PODAJ M: LICZBA PUNKTOW TABLICOWANIA M+1>'
      read (*,*) M
c     START
      h=(B-A)/M
      XX(0)=A
      do 3,k=1,M,1
      XX(k)=k*h
    3 continue
      do 4,k=0,M,1
      write(14,*) XX(k)
    4 continue
      write(*,*)'<ENTER> TO START'
      read (*,*)
      
      

      call DIVDIF(N,X,F)

      do 5,k=0,N,1
      write(10,*)'k=',k,' X=',X(k)
    5 continue
      do 6,k=0,N,1
      write(11,*)'k=',k,' F=',F(k)
    6 continue
c     HORNER NEWTON
      do 7,j=0,M,1
      P(j)=F(N)
      do 8,k=N,1,-1
      P(j)=F(k-1)+(XX(j)-X(k-1))*P(j)
    8 continue
    7 continue
c     ZAPIS NA DYSKU
      do 9,j=0,M,1
      write(12,*)'j=',j,' P=',P(j)
    9 continue
      
      write(*,*)'GOTOWE!'
      write(*,*)' 0 - KONIEC'
      write(*,*)' 1 - DA CAPO'
      read (*,*) pyt
      if(pyt.eq.0)goto 100
      goto 10
  100 stop 100
      end
      
      
      
      subroutine DIVDIF(N,X,F)
      implicit none
      
      real*8 X,F,ROB
      integer N,k,p
      dimension X(0:N),F(0:N),ROB(0:N)
      
      do 1,p=1,N,1
      do 2,k=p,N,1
      ROB(k)=(F(k)-F(k-1))/(X(k)-X(k-p))
    2 continue
      do 3,k=p,N,1
      F(k)=ROB(k)
    3 continue
    1 continue
      return
      end
      
