/** This library contains functions to solve functions using interval arithmetic
    techniques.

    THINK ABOUT:  Mix initial interval solve with a Newton's or secant method
    solver to find precise values after a certain depth?

    THINK ABOUT:  Find integer solutions only?

    THINK ABOUT:  When finding integer solutions, the same point can be checked
    up to 8 times (from different sides) in the 3-D case if the interval bound
    lies on an integer.  Figure out a way to only check it from one side.

    TODO: Make 2-D splitx look like 3D code

    TODO:  Remove duplicate results from 2D and 3D integer results.
*/

/** Solve a (1-argument) function, returning a sequence of intervals that
    contain solutions for the equation. */
intervalSolve[func, xmin, xmax, levels=53] :=
{
   eq = strToIntervalFunction[func]
   solutions = new OrderedList[{|x,y| infimum[x] <=> infimum[y]}]
   test1D[xmin, xmax, solutions, eq, levels, false]
   return coalesce1D[solutions]
}

/** Solve a (1-argument) function, returning an array of integers that 
    contain integer solutions for the equation. */
integerSolve[func, xmin, xmax, levels=53] :=
{
   eq = strToIntervalFunction[func]
   solutions = new OrderedList[{|x,y| infimum[x] <=> infimum[y]}]
   test1D[xmin, xmax, solutions, eq, levels, true]
   return solutions
}

/** Solve a (2-argument) function, returning a sequence of intervals that
    contain solutions for the equation. */
intervalSolve[func, xmin, xmax, ymin, ymax, levels=53] :=
{
   eq = strToIntervalFunction[func]
   solutions = new array
   test2D[xmin, xmax, ymin, ymax, solutions, eq, levels, false]
   return solutions
}

/** Solve a (2-argument) function, returning an array of integers that 
    contain integer solutions for the equation. */
integerSolve[func, xmin, xmax, ymin, ymax, levels=53] :=
{
   eq = strToIntervalFunction[func]
   solutions = new array
   test2D[xmin, xmax, ymin, ymax, solutions, eq, levels, true]
   return solutions
}

/** Solve a (3-argument) function, returning a sequence of intervals that
    contain solutions for the equation. */
intervalSolve[func, xmin, xmax, ymin, ymax, zmin, zmax, levels=53] :=
{
   eq = strToIntervalFunction[func]
   solutions = new array
   test3D[xmin, xmax, ymin, ymax, zmin, zmax, solutions, eq, levels, false]
   return solutions
}

/** Solve a (3-argument) function, returning an array of integers that 
    contain integer solutions for the equation. */
integerSolve[func, xmin, xmax, ymin, ymax, zmin, zmax, levels=53] :=
{
   eq = strToIntervalFunction[func]
   solutions = new array
   test3D[xmin, xmax, ymin, ymax, zmin, zmax, solutions, eq, levels, true]
   return solutions
}


/** Recursive function to test 1-dimensional interval.  This should not be
    called directly, but is called by intervalSolve above.
*/
test1D[x1, x2, solutions, eq, level, intsols] :=
{
   nextLevel = level - 1
   
   if intsols
   {
      x = new interval[x1, x2]
      xr = containsInteger[x]
      if xr == false
         return

      [x1, x2] = xr
   }

   xwidth = x2 - x1
   x = new interval[x1, x2]
   splitx = (xwidth != 0 xwidth)
   
   // Test the interval.  If it possibly contains solutions, recursively
   // subdivide.
   res = eval[eq]
   
   if res or res==undef
   {
      if (nextLevel >= 0) and ((!intsols) or (intsols and splitx))
      {
         cx = xwidth/2 + x1   // Done this way to allow date intervals
         test1D[x1, cx, solutions, eq, nextLevel, intsols]

         if splitx
            test1D[cx, x2, solutions, eq, nextLevel, intsols]
      } else
           if (res)             // Valid point
              solutions.insertUnique[new interval[x1,x2]]
           else
           {
              // Error in evaluating point but at bottom.  Possible solution?
              //solutions.insert[new interval[x1,x2]]
           }
   }
}

/** Recursive function to test 2-dimensional interval.  This should not be
    called directly, but is called by intervalSolve above.
*/
test2D[x1, x2, y1, y2, solutions, eq, level, intsols] :=
{
   nextLevel = level - 1
   
   if intsols
   {
      x = new interval[x1, x2]
      xr = containsInteger[x]
      if xr == false
         return

      [x1, x2] = xr

      y = new interval[y1, y2]
      yr = containsInteger[y]
      if yr == false
         return

      [y1, y2] = yr
   }

   xwidth = x2 - x1
   x = new interval[x1, x2]
   splitx = (xwidth != 0 xwidth)
   
   ywidth = y2 - y1
   y = new interval[y1, y2]
   splity = (ywidth != 0 ywidth)
   
   // Test the interval.  If it possibly contains solutions, recursively
   // subdivide.
   res = eval[eq]

   // if res == true
   //   println["x=$x, y=$y, res=$res"]
   
   if res or res==undef
   {
      if (nextLevel >= 0) and ((!intsols) or (intsols and (splitx or splity)))
      {
         cx = xwidth/2 + x1   // Done this way to allow date intervals
         cy = ywidth/2 + y1   // Done this way to allow date intervals

         test2D[x1, cx, y1, cy, solutions, eq, nextLevel, intsols]

         if splitx
            test2D[cx, x2, y1, cy, solutions, eq, nextLevel, intsols]

         if splity
            test2D[x1, cx, cy, y2, solutions, eq, nextLevel, intsols]

         if splitx and splity
            test2D[cx, x2, cy, y2, solutions, eq, nextLevel, intsols]
      } else
           if (res)             // Valid point
              solutions.push[[new interval[x1,x2], new interval[y1, y2]]]
           else
           {
              // Error in evaluating point but at bottom.  Possible solution?
              //solutions.insert[new interval[x1,x2]]
           }
   }
}

/** Recursive function to test 2-dimensional interval.  This should not be
    called directly, but is called by intervalSolve above.
*/
test3D[x1, x2, y1, y2, z1, z2, solutions, eq, level, intsols] :=
{
   nextLevel = level - 1
   
   if intsols
   {
      x = new interval[x1, x2]
      xr = containsInteger[x]
      if xr == false
         return

      [x1, x2] = xr

      y = new interval[y1, y2]
      yr = containsInteger[y]
      if yr == false
         return

      [y1, y2] = yr

      z = new interval[z1, z2]
      zr = containsInteger[z]
      if zr == false
         return

      [z1, z2] = zr
   }

   xwidth = x2 - x1
   x = new interval[x1, x2]
   splitx = (xwidth != 0 xwidth)
   
   ywidth = y2 - y1
   y = new interval[y1, y2]
   splity = (ywidth != 0 ywidth)
   
   zwidth = z2 - z1
   z = new interval[z1, z2]
   splitz = (zwidth != 0 zwidth)
   
   // Test the interval.  If it possibly contains solutions, recursively
   // subdivide.
   res = eval[eq]

   // if res == true
   //   println["x=$x, y=$y, z=$z, res=$res"]
   
   if res or res==undef
   {
      if (nextLevel >= 0) and ((!intsols) or (intsols and (splitx or splity or zplitz)))
      {
         cx = xwidth/2 + x1   // Done this way to allow date intervals
         cy = ywidth/2 + y1   // Done this way to allow date intervals
         cz = zwidth/2 + z1   // Done this way to allow date intervals

         test3D[x1, cx, y1, cy, z1, cz, solutions, eq, nextLevel, intsols]

         if splitx
            test3D[cx, x2, y1, cy, z1, cz, solutions, eq, nextLevel, intsols]

         if splity
            test3D[x1, cx, cy, y2, z1, cz, solutions, eq, nextLevel, intsols]

         if splitx and splity
            test3D[cx, x2, cy, y2, z1, cz, solutions, eq, nextLevel, intsols]

         if splitz
            test3D[x1, cx, y1, cy, cz, z2, solutions, eq, nextLevel, intsols]

         if splitx and splitz
            test3D[cx, x2, y1, cy, cz, z2, solutions, eq, nextLevel, intsols]

         if splity and splitz
            test3D[x1, cx, cy, y2, cz, z2, solutions, eq, nextLevel, intsols]

         if splitx and splity and splitz
            test3D[cx, x2, cy, y2, cz, z2, solutions, eq, nextLevel, intsols]
      } else
           if (res)             // Valid point
              solutions.push[[new interval[x1,x2], new interval[y1, y2], new interval[z1, z2]]]
           else
           {
              // Error in evaluating point but at bottom.  Possible solution?
              //solutions.insert[[new interval[x1,x2], new interval[y1, y2], new interval[z1, z2]]]
           }
   }
}

/** Coalesce an array of intervals, coalescing overlapping or touching
    intervals when possible.  It first sorts the intervals by their infimum
    and then performs a single sweep to coalesce intervals.
*/
coalesce1D[unions] :=
{
   len = length[unions]
   if len <= 1
      return unions

   ret = new array
   u = unions.shallowCopy[]
   sort[u, {|a,b| infimum[a] <=> infimum[b]}]
   curr = u@0
   
   for i = 1 to len-1
   {
      item = u@i
      if intersection[curr, item] == undef
      {
         // No intersection, push current interval
         ret.push[curr]
         curr = item
      }  else  // We have an intersection, update end to greater
      {
         si = supremum[item]
         sc = supremum[curr]
         if  si > sc
            curr = new interval[infimum[curr], si]
      }
   }

   ret.push[curr]

   return ret
}

/** Parses a string like "x = y" into an expression with interval-aware
    "possibly-equals" operators.  This also handles inequalities like "x <= y"
*/
strToIntervalFunction[func] :=
{
   func =~ %s/>=/ PGE /g    // Replace >= with possibly greater than or equals
   func =~ %s/<=/ PLE /g    // Replace <= with possibly less than or equals
   func =~ %s/!=/ PNE /g    // Replace != with possibly not equals
   func =~ %s/>/ PGT /g     // Replace > with possibly greater than
   func =~ %s/</ PLT /g     // Replace < with possibly less than
   func =~ %s/===/ PEQ /g   // Replace === with possibly equals
   func =~ %s/=/ PEQ /g     // Replace = with possibly equals
   eq = parseToExpression[func]
   //println["eq is \"" + inputForm[eq] + "\""]
   return eq
}

/** Tests if this interval contains an integer.  If not, it return false.  If
    it does contain an integer, this narrows the interval to its integer
    bounds and returns an array [imin, imax] with the integer bounds. */
containsInteger[ival] :=
{
   imin = ceil[infimum[ival]]
   imax = floor[supremum[ival]]
   //println["imin = $imin"]
   //println["imax = $imax"]
   if imin > imax
      return false
   return [imin, imax]
}


f = "cos[x] = x"
f = "cos[x] === 0"
f = "abs[x^2 - 100 x] = 0"
//f = "sin[40/x] = 0"  // This has an infinite number of solutions around 0
println[f]
//a = intervalSolve[f, -9., 10]
a = intervalSolve[f, -100000000.5, 100000000.5]
for c = a
   println[c]

f = "abs[x^2 - 100 x] = 0"
println[]
println[f]
a = integerSolve[f, -100000000.5, 100000000]
for c = a
   println[c]

// Two-variable solutions.
f = "x = y"
// Famous math puzzle, https://x.com/aap03102/status/849182047390363648
// Note that this uses an "and" conjunction to solve simultaneous equations!
f = "x^3 + x y^2 = 4640 y and x^2 y - y^3 = 537.6 x"
//f = "x (y + z)^-1 + (x + z)^-1 y + (x + y)^-1 z = 4"
//f = "(x-million)^2 + (y + 2 million)^2 = 10 billion"
println[]
println[f]
//a = integerSolve[f, -100.5, 100,5, -100.5, 100.5, 53]
a = integerSolve[f, -quadrillion - 13.5, quadrillion + 10.7, -quadrillion - 101.2, quadrillion + 10.5]
for c = a
   println[c]

// 3-dimensional integer solution
f = "x^2 + y^2 + z^2 = 8000"
println[]
println[f]
a = integerSolve[f, -20000.1, 20000, -20000, 20003, -20000, 20009]
for c = a
   println[c]