// Deriv.scala  compute weights for numerical derivatives
//          order is order of derivative, 1 = first derivative, 2 = second 
//          npoints is number of points where value of function is known    
//          f(x0), f(x1), f(x2) ... f(x npoints-1)                          
//          point is the term where derivative is computed                 
//          (b,a) = D.deriv(order, npoints, point)
//          f'(x point) = (1/bh^order)*( a(0)*f(x0) + a(1)*f(x1)        
//                       + ... + a(npoints-1)*f(x npoints-1) 
//          uniformly spaced points x1=x0+h, x2=x1+h=x1+2*h, ...           
//                                                                         
//          c = rderiv(order, npoints, point, h)
//          f'(x point) = a(0)*f(x0) + a(1)*f(x1)     
//                       + ... + a(npoints-1)*f(x npoints-1)
//
//          given non-uniformly spaced xi  in incresing order
//          c = nuderiv(order, npoints, point, x)
//          f'(x point) = a(0)*f(x0) + a(1)*f(x1)
//                       + ... + a(npoints-1)*f(x npoints-1)
//                                  
// algorithm: deriv uses 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                         
//            nuderiv solves for coefficients

object Deriv
{
   var S = Simeq


  def gcd(a:Int, b:Int):Int=
  {
    var a1:Int = 0
    var b1:Int = 0
    var q:Int = 0
    var r:Int = 0

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

  def deriv(order:Int, npoints:Int, point:Int):Array[Int]=
  {
    var A = Array.ofDim[Int](npoints+1)
    var b:Int = 0 // returned as A(npoints)

    // 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                          
  
    var h = Array.ofDim[Int](npoints+1) // coefficients of h 
    var x = Array.ofDim[Int](npoints+1) // coefficients of x, for differentiating 
    var numer = Array.ofDim[Int](npoints,npoints) // numerator of a term
                                                  // numer(term)(pos) 
    var denom = Array.ofDim[Int](100)  // denominator of coefficient 
    var jorder = 0
    var ipower = 0
    var isum = 0
    var iat = 0
    var r = 0
    var jj = 0
    var k = 0
    var debug = 0
  
    if(npoints <= order)
    {
      println("ERROR in call to deriv, npoints="+npoints+" < order+1="+
              order+1)
      return A
    }
    for(term <- 0 until npoints)
    {
      denom(term) = 1
      jj = 0
      for(i <- 0 until npoints-1)
      {
        if(jj == term)
        { 
          jj = jj + 1 // no index jj in jjth term
        } 
        numer(term)(i) = -jj // from (x-x(jj)) 
        denom(term) = denom(term)*(term-jj)
        jj = jj + 1
      }
    } // end term setting up numerator and denominator 
    if(debug>0)
    {
      for(i <- 0 until npoints)
      {
        println("denom("+i+")="+denom(i))
        for(j <- 0 until npoints-1)
        {
          println("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 <- 0 until npoints)
    {
      x(0) = 1 // initialize 
      x(1) = 0
      h(0) = 0
      h(1) = 0
      k = 1
      for(i <- 0 until npoints-1)
      {
        for(j <- 0 until k)          // multiply by (x + j h) 
        {
          h(j) = x(j)*numer(term)(i) // multiply be constant x(j) 
        }
        for(j <- 0 until k)          // multiply by x 
        {
          h(j+1) = h(j+1)+x(j)       // multiply by x and add 
        }
        k = k + 1
        for(j <- 0 until k)
        {
          x(j) = h(j)
        }
        x(k) = 0
        h(k) = 0
      }
      for(j <- 0 until k)
      {
        numer(term)(j) = x(j)
        if(debug>0)
        {
          println("p(x)="+numer(term)(j)+" x^"+j+" at term="+term)
        }
      } // have p(x) for this 'term' 

      // differentiate 'order' number of times 
      for(jorder <- 0 until order)
      {
        for(i <- 1 until k) // differentiate p'(x) 
        {
          numer(term)(i-1) = i*numer(term)(i)
        }
        k = k -1
      } // have p^(order) for this term 
      if(debug>0)
      {
        for(i <- 0 until k)
        {
          println("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 <- 0 until npoints) // get each term 
    {
      ipower = iat
      isum = numer(jterm)(0)
      for(j <- 1 until k)
      {
        isum = isum + ipower * numer(jterm)(j)
        ipower = ipower * iat
      }
      A(jterm) = isum
      if(debug>0)
      {
        println("f^("+order+")(x("+iat+")) = (1/h^"+order+") ("+
                A(jterm)+"/"+denom(jterm)+" f(x("+jterm+")) +")
      }
    }
    if(debug>0)
    {
      println("")
    }
    b = 0
    for(jterm <- 0 until npoints) // clean up each term 
    {
      if(A(jterm) != 0)
      {
        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 <- 0 until 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 <- 0 until npoints) // adjust for divisor 
    {
      A(jterm) = (A(jterm) * b) / denom(jterm)
      if(debug>0)
      {
        println("f^("+order+")(x("+iat+"))=(1/"+b+"h^"+
                A(jterm)+")("+jterm+" f(x("+jterm+")) +")
      }
    } // end computing terms of coefficients 
  
    // return (b, A) not working
    A(npoints) = b
    return A
  } // end deriv

  def rderiv(order:Int, npoints:Int, point:Int, hx:Double):Array[Double]=
  {
    var C = Array.ofDim[Double](npoints)
    var b:Int = 0
    var A = Array.ofDim[Int](npoints+1)
    var AA = Array.ofDim[Double](npoints) // for nuderiv
    var hpower:Double = 0.0

    if(order+npoints <= 11)
    {
      A = deriv(order, npoints, point)
      b = A(npoints)
      hpower = hx
      for(i <- 1  until order)
      {
        hpower = hpower * hx
      }
      hpower = 1.0/((b.toDouble)*hpower)
      for(i <- 0 until npoints)
      {
        C(i) = hpower*(A(i).toDouble)
      }
    }
    else
    {
      for(i <- 0 until npoints)
      {
        AA(i) = i.toDouble
      }
      C = nuderiv(order, npoints, point, AA)
      hpower = 1.0/hx
      for(i <- 1 until order)
      {
        hpower = hpower / hx
      }
      for(i <- 0 until npoints)
      {
        C(i)=hpower*C(i)
      }
    }
    return C
  } // end rderiv 


  // algorithm nuderiv  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
  //
  def nuderiv(order:Int, npoints:Int, point:Int, X:Array[Double]):Array[Double]=
  {
    var C = Array.ofDim[Double](npoints)
    var A = Array.ofDim[Double](npoints,npoints)
    var AI = Array.ofDim[Double](npoints,npoints)
    var fct = Array.ofDim[Double](npoints)
    var pwr:Double = 0.0
    var n:Int = npoints // X must be this or larger 
    var debug:Int = 0

    for(i <- 0 until n) // build matrix that will be inverted 
    {
      pwr = 1.0
      for(j <- 0 until n)
      {
        A(i)(j) = pwr
        pwr = pwr*X(i)
      }
    }
    AI = S.inverse(A) // invert the matrix 
    // form derivative for order and point 
    for(i <- 0 until n) // compute factors for order 
    {
      fct(i) = 1.0
    }
    for(j <- 0 until order)
    {
      for(i <- 0 until n)
      {
        fct(i) = fct(i)*((i-j).toDouble)
      }
    }
    for(i <- 0 until n)
    {
      C(i) = 0.0
      pwr = 1.0
      for(j <- order until n)
      { 
        C(i) = C(i) + fct(j)*AI(j)(i)*pwr
        pwr = pwr * X(point) 
      }
    }
    return C
  } // end nuderiv 

} // end Deriv.scala