# deriv.py3   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)                          
#             point is the term where derivative is computed                 
#             f'(x0) = (1/bh^order)*( a(0)*f(x0) + a(1)*f(x1)                
#                      + ... + a(points-1)*f(x points-1)                     
#             uniformly spaced points x1=x0+h, x2=x1+h=x1+2*h, ...           
#                                                                         
#             rderiv returns c[0]=(1/bh^order)*a[0]  ...                     
#             nuderiv returns c[0]=(1/bh^order)*a[0]  ...
#
# 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                         

from inverse import *
from math import *

def gcd(a, b):
  if a==0 or b==0:
    return 1

  if abs(a)>abs(b):
    a1 = abs(a)
    b1 = abs(b)
  else:
    a1 = abs(b)
    b1 = abs(a)

  r=1
  while r!=0:
    q = a1 / b1
    r = a1 - q * b1
    a1 = b1
    b1 = r

  return a1   # end gcd 


def deriv(order, npoints, point):  # returns a,b
  # 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,                 
  # 'b' returned with denominator of h^order                          
  
  a = [0 for i in range(npoints)]  # coefficients of derivative 
  h = [0 for i in range(npoints+1)]  # coefficients of h 
  x = [0 for i in range(npoints+1)]  # coefficients of x, for differentiating 
  numer = [[0 for j in range(npoints)] for i in range(npoints)]
  denom = [0 for i in range(npoints)]      # denominator of coefficient 
  debug = 0  # set to 1 for testing
  
  if npoints<=order:
    print("ERROR in call to deriv, npoints=",npoints," <= order=",order)
    return 1

  for term in range (0,npoints):
    denom[term] = 1
    j = 0
    for i in range(0,npoints-1):
      if j==term:
        j += 1 # no index j in jth term 
      numer[term][i] = -j # from (x-x[j]) 
      denom[term] = denom[term]*(term-j)
      j += 1
  # end term setting up numerator and denominator 

  if debug>0:
    for i in range(0,npoints):
      print("denom[",i,"]=",denom[i])
      for j in range(0,npoints-1):
        print("numer[",i,"][",j,"]=",numer[i][j])
        
  # 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 range(0,npoints):
    x[0] = 1 # initialize 
    x[1] = 0
    h[1] = 0
    k = 1
    for i in range(0,npoints-1):
      for j in range(0,k):
        h[j] = x[j]*numer[term][i] # multiply be constant x[j]

      for j in range(0,k):         # multiply by x 
        h[j+1] = h[j+1]+x[j]       # multiply by x and add 

      k += 1
      for j in range(0,k):
        x[j] = h[j]

      x[k] = 0
      h[k] = 0
    # end for i

    for j in range(0,k):
      numer[term][j] = x[j]
      if debug>0:
        print("p(x)= ",numer[term][j]," x^",j," at term=",term)
    # have p(x) for this 'term' 

    # differentiate 'order' number of times 
    for jorder in range(0,order):
      for i in range(1,k):  # differentiate p'(x) 
        numer[term][i-1] = i*numer[term][i]

      k -= 1
    # have p^(order) for this term
 
    if debug>0:
      for i in range(0,k):
        print("p^(",order,")(x)= ",numer[term][i]," x^",i," at term=",term)
  # end '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 
  for jterm in range(0,npoints): # get each term 
    ipower = iat
    isum = numer[jterm][0]
    for j in range(1,k):
      isum = isum + ipower * numer[jterm][j]
      ipower = ipower * iat

    a[jterm] = isum
    if debug>0:
      print("f^(",order,")(x[",iat,"]) = (1/h^",order,)
      print(") (",a[jterm],"/",denom[jterm]," f(x[",jterm,"]) +")

  if debug>0:
    print(" ")

  b = 0
  for jterm in range(0,npoints): # clean up each term 
    if a[jterm]==0:
      continue

    r = gcd(a[jterm],denom[jterm])
    a[jterm] = a[jterm] / r
    denom[jterm] = denom[jterm] /r
    if denom[jterm]<0:
      denom[jterm] = -denom[jterm]
      a[jterm] = - a[jterm]
    
    if denom[jterm]>b:
      b=denom[jterm]   # largest denominator 

  for jterm in range(0,npoints):  # check for divisor 
    r = (a[jterm] * b) / denom[jterm]
    if r * denom[jterm] != a[jterm] * b:
      r = gcd(b, denom[jterm])
      b = b * (denom[jterm] / r)
       
  for jterm in range (0,npoints):  # adjust for divisor 
    a[jterm] = (a[jterm] * b) / denom[jterm]
    if debug>0:
      print("f^(",order,")(x[",iat,"])=(1/",b,"h^",order)
      print(")(",a[jterm]," f(x[",jterm,"]) +  ")

  # end computing terms of coefficients 
  return a,b  # end deriv 


def rderiv(order, npoints, point, h):  # c array returned
  # 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,                    
  a = [0 for i in range(0,npoints)]
  c = [0 for i in range(0,npoints)]

  a,bb = deriv(order, npoints, point)
  hpower = h
  for i in range (1,order):
    hpower *= h

  hpower = 1.0/(bb*hpower)
  for i in range (0,npoints): 
    c[i]=hpower*a[i]

  return c  # end rderiv 

def nuderiv(order, npoints, point, xg):  # c array returned
  # compute double precision coefficients to numerically compute 
  # a derivative  of order 'order' using 'npoints'
  # at ordinate 'point', using x grid values in 'xg',xg[0]..xg[npoints-1]
  # where order>=1, npoints>=order+1, 0 <= point < npoints,
  # 'c' array returned with the coefficients.
  A  = [[0 for j in range(0,npoints)] for i in range(0,npoints)] 
  AI = [[0 for j in range(0,npoints)] for i in range(0,npoints)] 
  c  = [0 for i in range(0,npoints)]
  fct = [0 for i in range(0,npoints)]
  pwrjo = [0 for i in range(0,npoints)]

  for i in range (0,npoints): # build matrix that will be inverted 
    pwr=1.0
    for j in range (0,npoints):
      A[i][j] = pwr
      pwr = pwr*xg[i]


  AI = inverse(npoints,A)  # invert the matrix 
  # form derivative for order and point 
  for i in range (0,npoints): # compute factors for order
    fct[i] = 1.0

  for j in range (0,order):
    for i in range (0,npoints):
      fct[i] = fct[i]*(i-j)


  for i in range(0,npoints):
    c[i] = 0.0
    pwr = 1.0
    for j in range(order,npoints):
      c[i] = c[i] + fct[j]*AI[j][i]*pwr
      pwr = pwr * xg[point]


  return c # end nuderiv

# end deriv.py3