// rderiv.swift
//           compute formulas for numerical derivatives on equal spaced points
//           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)  equal spaced    
//           term is the point where derivative is computed
//           f'(x0) = cx[0]*f(x0) + cx[1]*f(x1)
//                    + ... + cx[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   

func rderiv(_ order:Int,_ npoints:Int,_ point:Int,_ hx:Double) -> [Double]{
  // var cx = [Double](repeating:0.0,count:npoints) // just test call
  // for i in 0..<npoints {
  // cx[i] = Double(point) + hx*(Double(i)+1.0)
  // } // end i
  //return cx
  
  var a = [Int](repeating:0,count:npoints)    // integer coefficients of solution
  var h = [Int](repeating:0,count:npoints+1)  // coefficients of h 
  var x = [Int](repeating:0,count:npoints+1)  // coefficients of x, variable for differentiating 
  var numer = [[Int]](repeating:[Int](repeating:0,count:npoints),count:npoints) // numerator of a term numer[term][pos] 
  var denom = [Int](repeating:0,count:npoints)  // denominator of coefficient 
  var c = [Double](repeating:0.0,count:npoints)  // returned
  var jj = 0
  var b  = 0
  var k  = 0
  var ipower = 0
  var isum = 0
  var iat = 0
  var r = 0
  var bh = 0.0

  if npoints <= order {
      print("ERROR in call to rderiv, npoints=\(npoints) < order=\(order)");
      return c;
  }
  for term in 0..<npoints {
    denom[term] = 1
    jj = 0
    for i in 0..<npoints-1 {
      if jj==term {
        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 += 1
    } // end i
  } // end term setting up numerator and denominator 
  // 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 += 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]
    } // end j  have p(x) for this 'term' 

    // differentiate 'order' number of times 
    for _ in 0..<order { // jorder
      for i in 1..<k { // differentiate p'(x) 
        numer[term][i-1] = i*numer[term][i]
      } // end i
      k -= 1
    } // end _ jorder   have p^(order) for this 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 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
  } // end jterm
  b = 0
  for jterm in 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]
    } // 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]
  } // end jterm   end computing terms of coefficients 
  bh = 1.0/Double(b)
  for _ in 0..<order { // i
    bh = bh/hx
  } // end _ i
  for i in 0..<npoints {
    c[i] = bh*Double(a[i])
  } // end i
  return c
} // end rderiv 

  
func gcd(_ a:Int,_ b:Int) -> Int {
  var a1 = 0
  var b1 = 0
  var r  = 0
  var q  = 0
    
  if a == 0 || b == 0 {
    a1 = 1
    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 

// end  rderiv.swift