# Deriv.rb  methods:  nderiv, rderiv, nuderiv
require_relative 'Inverse'

class Deriv

  def pown(x, pwr)
    val = 1.0
    for i in 0..pwr-1
      val = val * x
    end
    return val
  end  # pown 

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

  def nderiv(order, npoints, point, 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                          
    # nderiv  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   

    h=Array.new(40)      # coefficients of h 
    x=Array.new(40)      # coefficients of x, variable for differentiating 
    numer=Array.new(40){Array.new(40)}  # numerator of a term numer[term][pos] 
    denom=Array.new(40) # denominator of coefficient 
    debug = 0
  
    if npoints<=order
      puts "ERROR in call to deriv, npoints=#{npoints} < order+1=#{order+1}"
      return
    end
    for term in 0..npoints-1
      denom[term] = 1
      j = 0
      for i in 0..npoints-2
        if j==term
          j=j+1   # no index j in jth term 
        end
        numer[term][i] = -j   # from (x-x[j]) 
        denom[term] = denom[term]*(term-j)
        j=j+1
      end
    end  # term setting up numerator and denominator 
    if debug>0
      for i in 0..npoints-1
        puts "denom[#{i}]=#{denom[i]}"
        for j in 0..npoints-2
          puts "numer[#{i}][#{j}]=#{numer[i][j]}"
        end
      end
    end
        
    # 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
      if debug>0
        puts "working on term=#{term}"
      end  
      x[0] = 1   # initialize 
      x[1] = 0
      h[1] = 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 in 0..k-1          # multiply by x
          h[j+1] = h[j+1]+x[j]   # multiply by x and add 
        end
        k=k+1
        for j in 0..k-1
          x[j] = h[j]
        end
        x[k] = 0
        h[k] = 0
      end # i
    
      if debug>0
        puts "working on k=#{k}"
      end  

      for j in 0..k-1
        numer[term][j] = x[j]
        if debug>0
          puts "p(x)=#{numer[term][j]} x^#{j} at term=#{term}"
        end
      end   # have p(x) for this 'term' 

      # differentiate 'order' number of times 
      for jorder in 0..order-1
        for i in 1..k-1   # differentiate p'(x) 
          numer[term][i-1] = i*numer[term][i]
        end
        k=k-1
      end   # have p^(order) for this term 
      if debug>0
        for i in 0..k-1
          puts "p^(#{order})(x)= #{numer[term][i]} x^#{i} at term=#{term}"
        end
      end
    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 0..npoints-1   # get each term 
      ipower = iat
      isum = numer[jterm][0]
      for j in 1..k-1
        isum = isum + ipower * numer[jterm][j]
        ipower = ipower * iat
      end
      a[jterm] = isum
      if debug>0
        puts "f^(#{order})(x[#{iat}]) = (1/h^#{order}) (#{a[jterm]}/#{denom[jterm]} f(x[#{jterm}]) + \n"
      end
    end
    if debug>0
      puts " "
    end
    b = 0
    for jterm in 0..npoints-1   # clean up each term 
      if a[jterm]==0
        continue
      end
      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]
      end
      if denom[jterm]>b
         b=denom[jterm]   # largest denominator
      end 
    end
    for jterm in 0..npoints-1   # 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)
      end
    end
    for jterm in 0..npoints-1  # adjust for divisor 
      a[jterm] = (a[jterm] * b) / denom[jterm]
      if debug>0 
        puts "f^(#{order})(x[#{iat}])=(1/#{b}h^#{order})(#{a[jterm]} f(x[#{jterm}]) +"
      end
    end   # computing terms of coefficients 
    bb = b
  end   #  nderiv 

  def nuderiv(order, npoint, point, x, c) # x array in, c array out
    #  non uniformly spaced derivative coefficients 
    #    given  x0, x1, x2, x3  non uniformly spaced
    #    find the coefficients of y0, y1, y2, y3
    #    in order to compute the derivatives:
    #      y'(x0) = c00*y0 + c01*y1 + c02*y2 + c03*y3
    #      y'(x1) = c10*y0 + c11*y1 + c12*y2 + c13*y3
    #      y'(x2) = c20*y0 + c21*y1 + c22*y2 + c23*y3
    #      y'(x3) = c30*y0 + c31*y1 + c32*y2 + c33*y3 
    #
    #  Method:  fit data to  y(x)    = a + b*x + c*x^2 + d*x^3
    #                        y'(x)   =     b   + 2*c*x + 3*d*x^2
    #                        y''(x)  =           2*c   + 6*d*x
    #                        y'''(x) =                   6*d
    #
    #           with y0, y1, y2 and y3 variables:
    #           (Remember the y values are unknown at this time.
    #            The y values are the unknowns in the PDE.
    #            We are finding the c coefficients in order to
    #            build a system of simultaneous equations.)
    #
    #                              y0  y1  y2  y3 
    #  form:  | 1  x0  x0^2  x0^3   1   0   0   0 | = a
    #         | 1  x1  x1^2  x1^3   0   1   0   0 | = b
    #         | 1  x2  x2^2  x2^3   0   0   1   0 | = c
    #         | 1  x3  x3^2  x3^3   0   0   0   1 | = d
    #
    # reduce  | 1   0     0     0  a0  a1  a2  a3 | = a
    #         | 0   1     0     0  b0  b1  b2  b3 | = b
    #         | 0   0     1     0  c0  c1  c2  c3 | = c
    #         | 0   0     0     1  d0  d1  d2  d3 | = d
    #
    #                              this is just the inverse
    #
    #   y'(x0) =  b0*y0        + b1*y1        + b2*y2        + b3*y3        +
    #             2*c0*y0*x0   + 2*c1*y1*x0   + 2*c2*y2*x0   + 2*c3*y3*x0   +
    #             3*d0*y0*x0^2 + 3*d1*y1*x0^2 + 3*d2*y2*x0^2 + 3*d3*y3*x0^2
    #
    #   or  c00 = b0 + 2*c0*x0 + 3*d0*x0^2
    #       c01 = b1 + 2*c1*x0 + 3*d1*x0^2
    #       c02 = b2 + 2*c2*x0 + 3*d2*x0^2
    #       c03 = b3 + 2*c3*x0 + 3*d3*x0^2
    #
    #   y'(x0) = c00*y0 + c01*y1 + c02*y2 +c03*y3
    #
    #   y'(x1) =  b0*y0        + b1*y1        + b2*y2        + b3*y3        +
    #             2*c0*y0*x1   + 2*c1*y1*x1   + 2*c2*y2*x1   + 2*c3*y3*x1   +
    #             3*d0*y0*x1^2 + 3*d1*y1*x1^2 + 3*d2*y2*x1^2 + 3*d3*y3*x1^2
    #
    #   or  c10 = b0 + 2*c0*x1 + 3*d0*x1^2
    #       c11 = b1 + 2*c1*x1 + 3*d1*x1^2
    #       c12 = b2 + 2*c2*x1 + 3*d2*x1^2
    #       c13 = b3 + 2*c3*x1 + 3*d3*x1^2
    #
    #       cij = bj + 2*cj*xi + 3*dj*xi^2
    #
    #
    #   y'(x1) = c10*y0 + c11*y1 + c12*y2 +c13*y3
    #
    #   y''(x0) = 2*c0*y0    + 2*c1*y1    + 2*c2*y2    + 2*c3*y3    +
    #             6*d0*y0*x0 + 6*d1*y1*x0 + 6*d2*y2*x0 + 6*d3*y3*x0
    #
    #   or  c00 = 2*c0 + 6*d0*x0
    #       c01 = 2*c1 + 6*d1*x0
    #       c02 = 2*c2 + 6*d2*x0
    #       c03 = 2*c3 + 6*d3*x0
    #
    #   y'''(x0) = 6*d0*y0 + 6*d1*y1 + 6*d2*y2 + 6*d3*y3
    #
    #   or  c00 = 6*d0    independent of x
    #       c01 = 6*d1
    #       c02 = 6*d2
    #       c03 = 6*d3
    #
    n = npoint
    a = Array.new(n){Array.new(n)}
    ai = Array.new(n){Array.new(n)} # inverse
    fct = Array.new(n)
    puts "Deriv.new.nuderin running, npoint=#{n}"
    for i in 0...n  # build matrix that will be inverted 
      pwr=1.0
      for j in 0...n
        a[i][j] = pwr
        pwr = pwr*x[i]
      end # j
    end # i
    Inverse.new.invert(a, ai) # invert the matrix 
    # form derivative for order and point 
    for i in 0...n # compute factors for order 
      fct[i] = 1.0
    end # i
    for j in 0...order
      for i in 0...n
        fct[i] = fct[i]*(i-j)
      end # i
    end # j
    for i in 0...n
      c[i] = 0.0
      pwr = 1.0
      for j in order...n   
        c[i] = c[i] + fct[j]*ai[j][i]*pwr
        pwr = pwr*x[point]
      end # j 
    end # i
  end # end nuderiv 

  def rderiv(order, npoints, point, hx, c) # c[] coefficents returned
    # 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.new(npoints)    # integer coefficients of solution 
    h = Array.new(npoints+1)  # coefficients of h 
    x = Array.new(npoints+1)  # coefficients of x, variable for differentiating
    numer = Array.new(npoints){Array.new(npoints)} # numerator of a
                                                  # term numer[term][pos] 
    denom = Array.new(npoints)        # denominator of coefficient 
    debug = 0
    if npoints<=order
      puts "ERROR in call to rderiv, npoints=#{npoints} <= #{order}"
      return
    end # if
    for term in 0...npoints
      denom[term] = 1
      jj = 0
      for i in 0...(npoints-1)
        if jj==term
          jj = jj+1  # no index j in jth term
        end # if
        numer[term][i] = -jj # from (x-x[j]) 
        denom[term] = denom[term]*(term-jj)
        jj = jj+1
      end # i
    end # term setting up numerator and denominator 
    if debug>0
      for i in 0...npoints 
        puts "denom[#{i}]=#{denom[i]}"
        for j in 0...(npoints-1)
          puts "numer[#{i}][#{j}]=#{numer[i][j]}"
        end # j
      end # 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
      x[0] = 1 # initialize 
      x[1] = 0
      h[1] = 0
      k = 1
      for i in 0...(npoints-1)
        for j in 0...k            # multiply by (x + j h) 
          h[j] = x[j]*numer[term][i] # multiply be constant x[j] 
        end # j
        for j in 0...k         # multiply by x 
          h[j+1] = h[j+1]+x[j] # multiply by x and add 
        end # j
        k = k+1
        for j in 0...k
          x[j] = h[j]
        end # j
        x[k] = 0
        h[k] = 0
      end # i
      for j in 0...k
        numer[term][j] = x[j]
        if debug>0
          puts "p(x)= #{numer[term][j]} x^#{j} at term=#{term}"
        end # if
      end # j  # have p(x) for this 'term' 

      # differentiate 'order' number of times 
      for jorder in 0...order
        for i in 1...k  # differentiate p'(x) 
          numer[term][i-1] = i*numer[term][i]
        end # i
        k = k-1
      end # jorder   # have p^(order) for this term 
      if debug>0
        for i in 0...k
          puts "p^(#{order})(x)= #{numer[term][i]} x^#{i} at term=#{term}"
        end # i
      end # if
    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 0...npoints  # get each term 
      ipower = iat
      isum = numer[jterm][0]
      for j in 1...k
        isum = isum + ipower * numer[jterm][j]
        ipower = ipower * iat
      end # j
      a[jterm] = isum
      if debug>0
        puts "f^(#{order})(x[#{iat}]) = (1/h^#{order}) (#{a[jterm]}/#{denom[jterm]} f(x[#{jterm}])"
      end # if
      if debug>0
        puts " "
      end # if
    end # jterm
    b = 0
    for jterm in 0...npoints  # clean up each term 
      if a[jterm]==0
        next
      end # if
      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]
      end # if
      if denom[jterm]>b
        b=denom[jterm] # largest denominator
      end # if
    end # jterm
    for jterm in 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)
      end # if
    end # jterm
    for jterm in 0...npoints  # adjust for divisor 
      a[jterm] = (a[jterm] * b) / denom[jterm]
      if debug>0
        puts "f^(#{order})(x[#{iat}])=(1/#{b}h^#{order})(#{a[jterm]} f(x[#{jterm}])"
      end # if
    end # jterm computing terms of coefficients 
    bh = 1.0/b
    for i in 0...order
      bh = bh/hx
    end # i
    for i in 0...npoints
      c[i] = bh*a[i]
    end # i
  end # rderiv 

end # Deriv