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

def gcd(a, b):  # needed by deriv  
  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): # return 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                          
  
  a=[0 for i in range(npoints)]     # numerator of terms returned 
  h=[0 for i in range(npoints+1)]   # coefficients of h 
  x=[0 for i in range(npoints+1)]   # coefficients of x,
                                    # variable for differentiating 
  numer=[[0 for j in range(npoints)] for i in range(npoints+1)]
                                    # numerator of a term numer[term][pos] 
  denom=[0 for i in range(npoints)] # denominator of coefficient 
  debug = 0
  
  if npoints<=order:
    print "ERROR in call to deriv, npoints=",
    print npoints,
    print " < order=",
    print order,
    print "+1=%d \n"
    return a,0
  for term in range(npoints):
    denom[term] = 1
    j = 0
    for i in range(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(npoints):
      print "denom[",
      print i,
      print "]=",
      print denom[i]
      for j in range(npoints-1):
        print "numer[",
        print i,
        print "][",
        print j,
        print "]=",
        print 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(npoints):
    x[0] = 1 # initialize 
    x[1] = 0
    h[1] = 0
    k = 1
    for i in range(npoints-1):
      for j in range(k):            # multiply by (x + j h) 
         h[j] = x[j]*numer[term][i] # multiply be constant x[j] 
     
      for j in range(k):         # multiply by x 
        h[j+1] = h[j+1]+x[j] # multiply by x and add 
    
      k+=1
      for j in range(k):
        x[j] = h[j]
      if debug>0:
        print "k=",
        print k,
        print "x,h len=",
        print len(x)
      x[k] = 0
      h[k] = 0
  
    for j in range(k):
      numer[term][j] = x[j]
      if debug>0:
        print "p(x)= ",
        print numer[term][j],
        print " x^",
        print j,
        print " at term=",
        print term
    # have p(x) for this 'term' 

    # differentiate 'order' number of times 
    for jorder in range(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(k):
        print "p^(%d)(x)= %d x^%d at term=%d \n" \
              %( order, numer[term][i], i, 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(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^(",
      print order,
      print ")(x[",
      print iat,
      print "]) = (1/h^",
      print order,
      print ") (",
      print a[jterm],
      print "/",
      print denom[jterm],
      print " f(x[",
      print  jterm,
      print "]) + "

  if debug>0:
     print " "
  b = 0
  for jterm in range(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(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(npoints): # adjust for divisor 
    a[jterm] = (a[jterm] * b) / denom[jterm]
    if debug>0:
      print "f^(",
      print order,
      print ")(x[",
      print iat,
      print "])=(1/",
      print b,
      print "h^",
      print order,
      print ")(",
      print  a[jterm],
      print " f(x[",
      print jterm,
      print "]) +"

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