# nderiv.py print descrete derivatives of various orders at various terms
#           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)    

from deriv import *  # a[],b = deriv(order, npoints, point)
                     # c[]  = rderiv(order, npoints, point)

def pown(x, pwr):
  val = 1.0

  for i in range(0<pwr):
    val *= x
  return val
# end pow 

def check(order, npoints, point, a, b):
  hh = [0.125, 0.0625, 0.03125, 0.015625]

  
  m = 2
  for degree in range(npoints-2,npoints+2):
    if order>3 and degree==npoints+1:
      m=3 # extra check on error 
    for k in range(0,m):
      if k>0 and degree<npoints:
        break
      h = hh[k]
      if order>3:
        h=h/2.0 # cumulative, for reasonable error 
      if order>4:
        h=h/2.0
      if order>5:
        h=h/2.0
      pwr = degree
      sum = 0.0
      x = 1.0 # let x = 1.0, compute derivative using this x 
      for i in range(0,npoints):
        term = a[i]*pown(x, pwr) # x^pwr + x^pwr-1 + ... 1 
        sum = sum + term
        x = x + h

      sum = sum / b
      sum = sum / pown(h, order) # approximate derivative at 'point' 
      xpoint = 1.0+point*h # xpoint is x offset by 'point' 
      deriv = 1.0
      for i in range(0,order):
        deriv = deriv * pwr # analytic derivative 
        pwr = pwr -1

      deriv = deriv * pown(xpoint, pwr)
      err = deriv - sum
      print "  err=",err,", h=",h,", deriv=",deriv,
      print ", x=",xpoint,", pwr=",degree

# end check 


# nderiv main
norder = 6
npoints = 12
a = [0 for i in range(npoints)]  # coefficients of f(x0), f(x1), ... 
  
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' "
  
for order in range(1,norder+1):
  npoints = npoints -2 # fewer for larger derivatives 
  if npoints<=order+1:
    npoints = order + 2
  for points in range(order+1,npoints+1):
    for term in range(0,points):
      print "computing order=",order,", npoints=",points,", at term=",term
      a,b = deriv(order, points, term)
        
      for jterm in range(0,points):
        if jterm==0:
          print "f^(",order,")(x[",term,"])=(1/",b,"h^",order,
          print ")( ",a[jterm]," * f(x[",jterm,"]) +"

        elif jterm==points-1:
          print "                               ",a[jterm]," * f(x[",
          print jterm,"] ) + O(h^",points-order,")f^(",points,")(z)"

        else:
          print "                               ",a[jterm],
          print " * f(x[",jterm,"] +"
            
      # end jterm loop 

      check(order, points, term, a, b)
      print " "
    # end term loop 

    print " "
  # end points loop 

  print " "
# end order loop 

print "many formulas checked against Abramowitz and Stegun,"
print "Handbook of Mathematical Functions Table 25.2"

# end  nderiv.py