! deriv.f90  coefficients 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 rderiv(order, npoints, point, hx, cx)
  implicit none
  integer, intent(in) :: order, npoints, point
  double precision, intent(in) :: hx
  double precision, dimension(0:npoints-1), intent(out) :: cx

  ! compute double precision coefficients to numerically compute 
  ! a derivative  of order 'order' using 'npoints'               
  ! at ordinate point,                                           
  ! where order>=1, npoints>=order+1, 0 <= point < npoints,      
  ! 'c' array returned with the coefficients,                    
  integer i, bb
  double precision hpower
  integer, dimension(0:npoints-1) ::  a

  call deriv(order, npoints, point, a, bb)
  hpower = hx
  do i=2,order
    hpower = hpower * hx
  end do
  hpower = 1.0/(bb*hpower)
  do i=0,npoints-1
    cx(i)=hpower*a(i)
  end do
end subroutine rderiv 
! end deriv.f90