! nderiv.f90 compute formulas for numerical derivatives                     
!            order is order of derivative, 1 = first derivative, 2 = second 
!            points is number of points where value of function is known    
!            f(x0), f(x1), f(x2) ... f(x points-1)                          
!            term is the point where derivative is computer                 
!            f'(x0) = (1/bh)*( a(0)*f(x0) + a(1)*f(x1)                      
!                     + ... + a(points-1)*f(x points-1)                     

! algorithm: use divided differences to get polynomial p(x) that           
!            approximates  f(x).  f(x)=p(x)+error term                     
!            f'(x) = p'(x) + error term'                                   
!            substitute xj = x0 + j*h                                      
!            substitute x = x0 to get  p'(x0) etc   

function gcd(a, b) result(a1)
  implicit none
  integer, intent(in) :: a, b
  integer :: a1
  integer :: b1, q, r

  a1 = 0
  if(a .eq. 0 .or. b .eq. 0) return
  if(abs(a)>abs(b)) then
    a1 = abs(a)
    b1 = abs(b)
  else
    a1 = abs(b)
    b1 = abs(a)
  end if
  r=1
  do while(r .ne. 0)
    q = a1 / b1
    r = a1 - q * b1
    a1 = b1
    b1 = r
  end do
end function gcd 


subroutine deriv(order, npoints, point, a, bb)
  implicit none
  integer, intent(in) :: order, npoints, point
  integer, dimension(0:*), intent(out) :: a
  integer, intent(out) :: bb
  ! compute the exact coefficients to numerically compute a derivative 
  ! of order 'order' using 'npoints' at ordinate point,                
  ! where order>=1, npoints>=order+1, 0 <= point < npoints,            
  ! 'a' array returned with numerator of coefficients,                 
  ! 'bb' returned with denominator of h^order                          
  
  integer, dimension(0:99) :: h   ! coefficients of h 
  integer, dimension(0:99) :: x   ! coefficients of x,
                                  ! variable for differentiating 
  integer, dimension(0:99,0:99) :: numer ! numerator of a term numer(term,pos) 
  integer, dimension(0:99) :: denom      ! denominator of coefficient 
  integer :: i, j, k, term, b
  integer :: jorder, ipower, isum, iat, jterm, r
  integer :: debug = 0

  interface
    function gcd(a, b) result(a1)
      integer, intent(in) :: a, b
      integer :: a1
    end function gcd
  end interface
  
  if(npoints .le. order) then
    print *, 'ERROR in call to deriv, npoints=', npoints, &
             ' <= order=', order
    return
  end if
  do term=0,npoints-1
    denom(term) = 1
    j = 0
    do i=0,npoints-2
      if(j .eq. term) j=j+1 ! no index j in jth term 
      numer(term,i) = -j ! from (x-x(j)) 
      denom(term) = denom(term)*(term-j)
      j=j+1
    end do
  end do ! term setting up numerator and denominator 
  if(debug>0) then
    do i=0,npoints-1
      print *, 'denom(', i, ')= ', denom(i)
      do j=0,npoints-2
        print *, 'numer(', i, ',', j, ')=', numer(i,j)
      end do
    end do
  end if

  ! substitute xj = x0 + j*h, just keep j value in numer(term,i) 
  ! don't need to count x0 because same as x's 
  ! done by above setup 
      
  ! multiply out polynomial in x with coefficients of h 
  ! starting from (x+j1*h)*(x+j2*h) * ... 
  ! then differentiate d/dx 
  do term=0,npoints-1
    x(0) = 1 ! initialize 
    x(1) = 0
    h(1) = 0
    k = 1
    do i=0,npoints-2
      do j=0,k-1            ! multiply by (x + j h) 
        h(j) = x(j)*numer(term,i) ! multiply be constant x(j) 
      end do
      do j=0,k-1         ! multiply by x 
        h(j+1) = h(j+1)+x(j) ! multiply by x and add 
      end do
      k=k+1
      do j=0,k-1
        x(j) = h(j)
      end do
      x(k) = 0
      h(k) = 0
    end do ! on i
    do j=0,k-1
      numer(term,j) = x(j)
      if(debug>0)print *, 'p(x)= ', numer(term,j),' x^', j, ' at term=', term
    end do ! have p(x) for this 'term' 

    ! differentiate 'order' number of times 
    do jorder=0,order-1
      do i=1,k-1 ! differentiate p'(x) 
        numer(term,i-1) = i*numer(term,i)
      end do
      k=k-1
    end do ! have p^(order) for this term 
    if(debug>0) then
      do i=0,k-1
        print *, 'p^(', order, ')(x)= ', numer(term,i) , &
                 ' x^', i,' at term=', term
      end do
    end if
  end do ! 'term' loop  for one 'order', for one 'npoints' 

  ! substitute each point for x, evaluate, get coefficients of f(x(j)) 
  iat = point ! evaluate at x(iat), substitute 
  do jterm=0,npoints-1 ! get each term 
    ipower = iat
    isum = numer(jterm,0)
    do j=1,k-1
      isum = isum + ipower * numer(jterm,j)
      ipower = ipower * iat
    end do
    a(jterm) = isum
    if(debug>0) then
      print *, 'f^(', order, ')(x(', iat, ')) = (1/h^', order, ') (', &
               a(jterm), '/', denom(jterm), ' f(x(', jterm, ')) + '
    end if
  end do
  if(debug>0) print *, ' '
  b = 0
  do jterm=0,npoints-1 ! clean up each term 
    if(a(jterm) .eq. 0) cycle
    r = gcd(a(jterm),denom(jterm))
    a(jterm) = a(jterm) / r
    denom(jterm) = denom(jterm) /r
    if(denom(jterm)<0) then
      denom(jterm) = -denom(jterm)
      a(jterm) = - a(jterm)
    end if
    if(denom(jterm)>b) b=denom(jterm) ! largest denominator 
  end do
  do jterm=0,npoints-1 ! check for divisor 
    r = (a(jterm) * b) / denom(jterm)
    if(r * denom(jterm) .ne. a(jterm) * b) then
      r = gcd(b, denom(jterm))
      b = b * (denom(jterm) / r)
    end if
  end do
  do jterm=0,npoints-1 ! adjust for divisor 
    a(jterm) = (a(jterm) * b) / denom(jterm)
    if(debug>0) then
      print *, 'f^(', order, ')(x(', iat, '))=(1/', b, 'h^', order, &
               ',', a(jterm), ' f(x(', jterm, ')) + '
    end if
  end do ! computing terms of coefficients 
  bb = b
end subroutine deriv 


subroutine check(order, npoints, point, a, b)
  implicit none
  integer, intent(in) :: order, npoints, point, b
  integer, dimension(0:*), intent(in) :: a
  integer :: pwr ! simple test using x^pwr 
  double precision :: sum, x, term, err, xpoint, deriv
  double precision :: h
  integer :: i, degree
  integer :: debug = 0

  if(debug>0) print *, 'check order=', order, ', points=', npoints, &
                       ', point=', point, ', a(0)=', a(0), ', b=', b
  print *, ' '
  h = 0.125
  if(order>3) h=h/2.0d0 ! cumulative
  if(order>4) h=h/2.0d0
  if(order>5) h=h/2.0d0
  do degree=npoints-2,npoints+1
    pwr = degree
    sum = 0.0d0
    x = 1.0d0
    do i=0,npoints-1
      term = dble(a(i))*(x**pwr)
      sum = sum + term
      x = x + h
    end do
    sum = sum / dble(b)
    sum = sum / (h**order)
    xpoint = 1.0d0+dble(point)*h
    deriv = 1.0d0
    do i=0,order-1
      deriv = deriv * pwr
      pwr = pwr -1
    end do
    deriv = deriv * (xpoint**pwr)
    err = deriv - sum
    print '("               err=",g9.3,", h=",g9.3," deriv=",g9.3," x=",g9.3," pwr=",i2)', &
            err, h, deriv, xpoint, degree
  end do
end subroutine check 


program nderiv
  implicit none
  integer :: norder = 6
  integer :: npoints = 12
  integer :: order, points, term, jterm
  integer, dimension(0:50) :: a   ! coefficients of f(x0), f(x1), ... 
  integer :: b

  interface
    subroutine deriv(order, npoints, point, a, bb)
      integer, intent(in) :: order, npoints, point
      integer, dimension(0:*), intent(out) :: a
      integer, intent(out) :: bb
    end subroutine deriv
    subroutine check(order, npoints, point, a, b)
      integer, intent(in) :: order, npoints, point, b
      integer, dimension(0:*), intent(in) :: a
    end subroutine check
  end interface

  print *, 'Formulas for numerical computation of derivatives '
  print *, '  f^(3) means third derivative '
  print *, '  f^(2)(x(1)) means second derivative computed at point x(1) '
  print *, '  f(x(2)) means function evaluated at point x(2) '
  print *, '  Points are x(0), x(1)=x(0)+h, x(2)=x0+2h, ... '
  print *, '  The error term gives the power of h and '
  print *, '  derivative order with z chosen for maximum value.'
  print *, '  The error terms, err=, are for test polynomial sum(x^pwr) '
  print *, '  with h = 0.125 or smaller, and order = "pwr" '
  
  do order=1,norder
    npoints = npoints -2 ! fewer for larger derivatives 
    if(npoints .le. order+1) npoints = order + 2
    do points=order+1,npoints
      do term=0,points-1
        print '("computing order=",I2,", npoints=",I3,", at term=",I3)', order, points, term 
        call deriv(order, points, term, a, b)
        
        do jterm=0,points-1
          if(jterm .eq. 0) then
            print *, 'f^(', order, ')(x(', term, '))=(1/', b, 'h^', order, &
                   ',', a(jterm), '* f(x(', jterm, ')) + '
          else if(jterm .eq. points-1) then
            print *, '                             ', a(jterm), &
                     '* f(x(', jterm, ') ) + O(h^', points-order, &
                     ')f^(', points, ')(z) '
          else
            print *, '                             ', a(jterm), &
                     '* f(x(', jterm, ') + '
          end if
        end do ! jterm loop
        call check(order, points, term, a, b)
        print *, ' '
      end do ! term loop
      print *, ' '
    end do ! points loop 
    print *, ' '
  end do ! order loop 
  print *, 'many formulas checked against Abramowitz and Stegun,'
  print *, 'Handbook of Mathematical Functions Table 25.2'
end program nderiv