! lagrange.f90  simple definition and build 
!             L_n,j(x) = product i=0,n i!=j (x-x_i)/(x_j-x_i)  
!             let 0 < x < 1 for points x_0, x_1, ... , x_n
!             let h = x_i - x_i-1,  x_i = i * h,  h = 1/n
!             L_n(x) = sum j=0,n f(x_j) L_n,j(x)
!             This L_n(x) fits f(x) with a n th order polynomial
!
!  generate L_n,j(x) for n=1,6  j=0,n

 
program lagrange
  implicit none
  integer :: NN = 6
  integer :: i, j, k, n, m
  double precision :: h, termx, termc
  double precision, dimension(0:20) :: p ! coefficients of polynomial 
  double precision, dimension(0:20) :: q ! temp 
  double precision, dimension(0:20) :: L ! accumulate fit of e^x 
  double precision :: x, fit, err, maxerr, rmserr

  print *, "lagrange.f90 running on 0<x<1 for N=", NN
  print *, " L_order,point  order is highest power of x, point 0..n "
  print *, "      L_n,j(x) = product i=0,n i!=j (x-x_i)/(x_j-x_i)  "
  print *, "      let 0 < x < 1 for points x_0, x_1, ... , x_n "
  print *, "      let h = x_i - x_i-1,  x_i = i * h,  h = 1/n "
  print *, "      L_n(x) = sum j=0,n f(x_j) L_n,j(x) "
  print *, "      This L_n(x) fits f(x) with a n th order polynomial "

  do n=1,NN
    h = 1.0d0/n
    do j=0,n
      ! print *, "computing terms of L_",n,",",j 
      p(0) = 1.0d0 ! identity polynomial to start product 
      do i=1,n
        p(i) = 0.0d0
      end do ! on i
      do i=0,n
        if(i .ne. j) then
          termx = 1.0/((j-i)*h) ! (x    )/(x_j-x_i) 
          termc = -h*i*termx    ! ( -x_i)/(x_j-x_i) 
          ! do product p() * (termx * x + termc) 
          do k=0,n
            q(k) = p(k) * termc
          end do ! on k
          do k=1,n
            q(k) = q(k) + p(k-1)*termx
          end do ! on k
          do k=0,n
            p(k) = q(k)
          end do ! on k
        end if
      end do ! on i
      ! print L_n,j(x) = p(...) 
      print *, "L_",n,",",j,"(x) = ",p(0)," x^0 +"
      do i=1,n
         print *, "           ",p(i)," x^",i," + "
      end do ! on i
      print *, " "
    end do ! on j 
    print *, " "
  end do ! on n 

  ! repeat generation of L_5,* and fit e^x  0<x<1 
  print *, " "
  n = 5
  h = 1.0d0/n
  do i=0,n
    L(i) = 0.0 ! accumulate sum
  end do ! on i 
  do j=0,n
    p(0) = 1.0 ! identity polynomial to start product 
    do i=1,n
      p(i) = 0.0
    end do ! on i
    do i=0,n
      if(i .ne. j) then
        termx = 1.0d0/(dble(j-i)*h) ! (x    )/(x_j-x_i) 
        termc = -h*dble(i)*termx    ! ( -x_i)/(x_j-x_i) 
        ! do product p() * (termx * x + termc) 
        do k=0,n
          q(k) = p(k) * termc
        end do ! on k
        do k=1,n
          q(k) = q(k) + p(k-1)*termx
        end do ! on k
        do k=0,n
          p(k) = q(k)
        end do ! on k
      end if
    end do ! on i
    do k=0,n
      L(k) = L(k)+exp(j*h)*p(k)
    end do ! on k
  end do ! on j 
  ! check results at computed points and intermediate points 
  print *, "check Lagrange 5 point fit of exp(x) in 0 <= x <=1 "
  print *, "L_",n,"(x) = ",L(0)," x^0 +"
  do i=1,n
      print *, "        ",L(i)," x^",i," + "
  end do ! on i
  print *, " "
  m = 2*n
  h = 1.0d0/m
  maxerr = 0.0d0
  rmserr = 0.0d0
  do k=0,m
    x = dble(k)*h
    fit = L(n)*x ! evaluate L at x 
    do i=n-1,1,-1
      fit = (L(i)+fit)*x
    end do ! on i
    fit = fit + L(0)
    err = fit - exp(x)
    maxerr = max(maxerr, abs(err))
    rmserr = rmserr + err*err
    print *, "check exp(",x,")=",exp(x),", fit=",fit,", err=", err
  end do ! on k
  rmserr = sqrt(rmserr/dble(m+1))
  print *, "maxerr=",maxerr,", rmserr", rmserr

  ! repeat generation of L_6,* and fit e^x  0<x<1 
  print *, " "
  n = 6
  h = 1.0d0/n
  do i=0,n
    L(i) = 0.0 ! accumulate sum
  end do ! on i 
  do j=0,n
    p(0) = 1.0 ! identity polynomial to start product 
    do i=1,n
      p(i) = 0.0
    end do ! on i
    do i=0,n
      if(i .ne. j) then
        termx = 1.0/((j-i)*h) ! (x    )/(x_j-x_i) 
        termc = -h*i*termx    ! ( -x_i)/(x_j-x_i) 
        ! do product p() * (termx * x + termc) 
        do k=0,n
          q(k) = p(k) * termc
        end do ! on k
        do k=1,n
          q(k) = q(k) + p(k-1)*termx
        end do ! on k
        do k=0,n
          p(k) = q(k)
        end do ! on k
      end if
    end do ! on i
    do k=0,n
      L(k) = L(k)+exp(j*h)*p(k)
    end do ! on k
  end do ! on j 
  ! check results at computed points and intermediate points 
  print *, "check Lagrange 6 point fit of exp(x) in 0 <= x <=1 "
  print *, "L_",n,"(x) = ",L(0)," x^0 +"
  do i=1,n
      print *, "        ",L(i)," x^",i," + "
  end do ! on i
  print *, " "
  m = 2*n
  h = 1.0d0/m
  maxerr = 0.0
  rmserr = 0.0
  do k=0,m
    x = k*h
    fit = L(n)*x ! evaluate L at x 
    do i=n-1,1,-1
      fit = (L(i)+fit)*x
    end do ! on i
    fit = fit + L(0)
    err = fit - exp(x)
    maxerr = max(maxerr, abs(err))
    rmserr = rmserr + err*err
    print *, "check exp(%g)=%g, fit=%g, err=%g ", x, exp(x), fit, err
  end do ! on k
  rmserr = sqrt(rmserr/dble(m+1))
  print *, "maxerr=",maxerr," rmserr=",rmserr

  print *, " end lagrange.f90 "
end program lagrange