// Secant method solver.
//
// This function finds a root of the function f.
//  In other words, it returns the value of x for which f[x] = 0.
//
// The arguments are:
//   f: A function that takes a single argument
//   x1, x2:  initial guesses
//   maxDelta:  the maximum error in x
secant[f, x1, x2, maxDelta = 1e-14] :=
{
   f1 = f[x1]
   f2 = f[x2]
   x = undef
   while (true)
   {
      diff = f1 - f2
      if diff == 0
         return x1
      
      x = x1 - (f1 * (x1 - x2)) / diff
      // println[x]

      if abs[x - x1] < maxDelta
         return x

      x2 = x1
      x1 = x
      f2 = f1
      f1 = f[x]
   }
}

// This uses the secant method to invert the function y = f[x].
// This will essentially find an inverse function for f[x] and return a value
//  of x for which f[x] = y.
//  other parameters:
//   x1,x2:  initial guesses that hopefully bound the desired result.
//   maxDelta:  maximum error in y

//   TODO:  Use interval techniques to make this more rigorous and powerful?
//   TODO:  Automatically make guesses for x1 and x2?  Somehow?
secantInvert[f, y, xmin, xmax, maxDelta = 1e-14] :=
{
   x1 = xmin
   x2 = xmax
   y1 = f[x1]
   y2 = f[x2]
   xnew = (x2-x1)/2 + x1
   while true
   {
      ydiff = y2 - y1
//      println["ydiff is $ydiff, x1 is $x1, x2 is $x2"]
      if ydiff == 0 y   // Degenerate case to avoid dividing by zero.
         return xnew    // This may not be always a correct solution?
      
      invSlope = (x2-x1) / ydiff
      xnew = x1 + (y - y1) invSlope
      if xnew < xmin
         xnew = xmin
      if xnew > xmax
         xnew = xmax
      
      ynew = f[xnew]
//      println["xnew=$xnew\tynew=$ynew"]

      if ynew == 0 y    // Degenerate case to avoid dividing by zero.
         return xnew      // This may not be always a correct solution?
      
      if abs[(ynew - y) / ynew] < maxDelta
         return xnew

      y2 = y1
      y1 = ynew
      x2 = x1
      x1 = xnew
   }
}

// This uses the secant method to invert the function y = f[x], assuming that
// y is an angle.  This prevents some overcorrections when angles are
// negative.
// This will essentially find an inverse function for f[x] and return a value
//  of x for which f[x] = y.
//  other parameters:
//   x1,x2:  initial guesses that hopefully bound the desired result.
//   maxDelta:  maximum error in y

//   TODO:  Use interval techniques to make this more rigorous and powerful?
//   TODO:  Automatically make guesses for x1 and x2?  Somehow?
secantInvertAngle[f, y, xmin, xmax, maxDelta = 0.015 arcsec] :=
{
   x1 = xmin
   x2 = xmax
   y1 = f[x1]
   y2 = f[x2]
   xnew = (x2-x1)/2 + x1
   while true
   {
      y1e = y1 - y
      y2e = y2 - y
      if y1e > 180 deg
         y1 = y1 - circle
      if y2e > 180 deg
         y2 = y2 - circle
      if y1e < -180 deg
         y1 = y1 + circle
      if y2e < -180 deg
         y2 = y2 + circle
      
      ydiff = y2 - y1
      
      // did we wrap around the circle?
      if abs[y1 - y2] > 180 degrees
         if (y1 < y2)
            y2 = y2 - circle
         else
            y1 = y1 - circle

//      println["ydiff is $ydiff, x1 is $x1, x2 is $x2, y1 is $y1, y2 is $y2"]
      if ydiff == 0 y   // Degenerate case to avoid dividing by zero.
         return xnew    // This may not be always a correct solution?
      
      invSlope = (x2-x1) / ydiff
      xnew = x1 + (y - y1) invSlope
      if xnew < xmin
         xnew = xmin
      if xnew > xmax
         xnew = xmax
      
      ynew = f[xnew]
//      println["xnew=$xnew\tynew=$ynew"]

      if ynew == 0 y    // Degenerate case to avoid dividing by zero.
         return xnew      // This may not be always a correct solution?
      
      if abs[(ynew - y) / ynew] < maxDelta
         return xnew

      y2 = y1
      y1 = ynew
      x2 = x1
      x1 = xnew
   }
}

// Minimize a function using the secant method.  This doesn't really work yet.
secantMinimize[f, xmin, xmax, minStepX] :=
{
   x1 = xmin
   x2 = xmax
   y1 = f[x1]
   y2 = f[x2]
   while true
   {
      println["x1=$x1\t x2=$x2"]
      diff = x2-x1
      if diff == 0
         return f[x1]
      
      slope = (y2-y1)/diff
      xnew = x1 + slope (x1+x2)/2
      ynew = f[xnew]
      println["ynew=$ynew\txnew= $xnew"]

      if (abs[x2-x1] < minStepX)
         return ynew
      
      y2 = y1
      y1 = ynew
      x2 = x1
      x1 = xnew
      if x1 > x2
         [x1, x2] = [x2, x1]
      if x1 < xmin
         x1 = xmin
      if x2 > xmax
         x2 = xmax
   }
}

// Sample root-finding:
//   Define a procedure block that represents the equation
// (this is just a function without a name, or think of it
//  as a reference to a function.)
//f = { |x| ln[x] - 1}
//println["Solution: " + secant[f, 1, 3]]


// Sample inverse-finding:
//   Find an inverse for the following function.
//   The call below finds a value x such that log[x]=2
//     in other words, calculates 10^2
// f = { |x| log[x] }
// println[secantInvert[f, 2, 1, 200, 1e-20]]

// Example secant method using arbitrary precision to calculate sin[x]
/*
use ArbitraryPrecision.frink
use pi2.frink
digits = 75
setPrecision[digits]
f ={|x| arbitrarySin[x, 75]}
y = secant[f, 3, 4, 10^-digits]
println["Solution is $y"]
println["Difference from pi is " + (y - Pi.getPi[digits])]
*/