# deriv.jl  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   

using printf

function gcd(a::Int, b::Int)::Int
  if a==0 || b==0
    return 0
  end
  if abs(a)>abs(b)
    a1 = abs(a)
    b1 = abs(b)
  else
    a1 = abs(b)
    b1 = abs(a)
  end # if
  r = 1
  while r!=0
    q = floor(a1 / b1)
    r = a1 - q * b1
    a1 = b1
    b1 = r
  end 
  return a1
end # gcd 

function deriv(order::Int, npoints::Int, point::Int, a, 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                          

  h = zeros(Float64, npoints)  # coefficients of h 
  x = zeros(Float64, npoints)  # coefficients of x,
                               # variable for differentiating 
  numer = zeros(Float64, npoints, npoints)
                           # numerator of a term numer[term][pos] 
  denom = zeros(Float64, npoints)  # denominator of coefficient 
  debug = 1

  
  if npoints <= order
    print("ERROR in call to deriv, npoints=")
    print(npoints)
    print(" <= order=")
    println(order)
    return nothing
  end # if 
  for term in 0:npoints-1
    denom[term] = 1.0
    j = 0
    for i in 0:npoints-2
      if j == term
        j = j+1 # no index j in jth term
      end # if
      numer[term][i] = -j # from (x-x(j)) 
      denom[term] = denom[term]*(term-j)
      j = j+1
    end # for i
  end  # end term setting up numerator and denominator 
  if debug>0 
    for i in 0:npoints-1
      print("denom(", i, ")= ", denom[i])
      for j in 0:npoints-2
        print("numer[", i, "][", j, "]=", numer[i][j])
      end # for j
    end # for i
  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 
  for term in 0:npoints-1
    x[0] = 1.0 # initialize 
    x[1] = 0.0
    h[1] = 0.0
    k = 1
    for i in 0:npoints-2
      for j in 0:k-1            # multiply by (x + j h) 
        h[j] = x[j]*numer[term][i] # multiply be constant x(j) 
      end # for j
      for j in 0:k-1         # multiply by x 
        h[j+1] = h[j+1]+x[j] # multiply by x and add 
      end # for j
      k = k+1
      for j in 0:k-1
        x[j)] = h[j]
      end # end j
      x[k] = 0.0
      h[k] = 0.0
    end # for 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
  return nothing
end # deriv 

function rderiv(order::Int, npoints::Int, point::Int, hx, cx)
  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

  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 

function rderiv1(order::Int, npoints::Int, point::Int, hx, cx)
  double precision, intent(in) :: hx
  double precision, dimension(npoints), 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, 1 <= 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-1, a, bb)
  hpower = hx
  for i in 2:order
    hpower = hpower * hx
  end # for
  hpower = 1.0/(bb*hpower)
  for i in 0:npoints-1
    cx(i+1)=hpower*a(i)
  end # for
  return nothing
end # rderiv1 
# end deriv.jl