/** 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
   {
      [x1, x2] = narrowInteger[x1, x2]
      if x1 == undef
         return
   }

   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))
      {
         [x1, cx1, cx2, x2] = splitInt[x1, x2, intsols]
         if x1 != undef
            test1D[x1, cx1, solutions, eq, nextLevel, intsols]

         if splitx and cx2 != undef
            test1D[cx2, 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
   {
      [x1, x2] = narrowInteger[x1, x2]
      if x1 == undef
         return

      [y1, y2] = narrowInteger[y1, y2]
      if y1 == undef
         return
   }

   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)))
      {
         [x1, cx1, cx2, x2] = splitInt[x1, x2, intsols]
         [y1, cy1, cy2, y2] = splitInt[y1, y2, intsols]

         if x1 != undef and y1 != undef
            test2D[x1, cx1, y1, cy1, solutions, eq, nextLevel, intsols]

         if splitx and cx2 != undef and y1 != undef
            test2D[cx2, x2, y1, cy1, solutions, eq, nextLevel, intsols]

         if splity and x1 != undef and cy2 != undef
            test2D[x1, cx1, cy2, y2, solutions, eq, nextLevel, intsols]

         if splitx and splity and cx2 != undef and cy2 != undef
            test2D[cx2, x2, cy2, 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
   {
      [x1, x2] = narrowInteger[x1, x2]
      if x1 == undef
         return

      [y1, y2] = narrowInteger[y1, y2]
      if y1 == undef
         return

      [z1, z2] = narrowInteger[z1, z2]
      if z1 == undef
         return
   }

   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 splitz)))
      {
         [x1, cx1, cx2, x2] = splitInt[x1, x2, intsols]
         [y1, cy1, cy2, y2] = splitInt[y1, y2, intsols]
         [z1, cz1, cz2, z2] = splitInt[z1, z2, intsols]

         if x1 != undef and y1 != undef and z1 != undef
            test3D[x1, cx1, y1, cy1, z1, cz1, solutions, eq, nextLevel, intsols]

         if splitx and cx2 != undef and y1 != undef and z1 != undef
            test3D[cx2, x2, y1, cy1, z1, cz1, solutions, eq, nextLevel, intsols]

         if splity and x1 != undef and cy2 != undef and z1 != undef
            test3D[x1, cx1, cy2, y2, z1, cz1, solutions, eq, nextLevel, intsols]

         if splitx and splity and cx2 != undef and cy2 != undef and z1 != undef
            test3D[cx2, x2, cy2, y2, z1, cz1, solutions, eq, nextLevel, intsols]

         if splitz and x1 != undef and y1 != undef and cz2 != undef
            test3D[x1, cx1, y1, cy1, cz2, z2, solutions, eq, nextLevel, intsols]

         if splitx and splitz and cx2 != undef and y1 != undef and cz2 != undef
            test3D[cx2, x2, y1, cy1, cz2, z2, solutions, eq, nextLevel, intsols]

         if splity and splitz and x1 != undef and cy2 != undef and cz2 != undef
            test3D[x1, cx1, cy2, y2, cz2, z2, solutions, eq, nextLevel, intsols]

         if splitx and splity and splitz and cx2 != undef and cy2 != undef and cz2 != undef
            test3D[cx2, x2, cy2, y2, cz2, z2, solutions, eq, nextLevel, intsols]
      } else
           if (res)             // Valid point
           {
              solutions.push[[new interval[x1,x2], new interval[y1, y2], new interval[z1, z2]]]
           // println["putting $x1, $x2   $y1, $y2   $z1, $z2  splitx=$splitx splity=$splity splitz=$splitz"]
           } 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 undef.  If
    it does contain an integer, this narrows the interval to its integer
    bounds and returns an array [imin, imax] with the integer bounds. */
narrowInteger[x1, x2] :=
{
   imin = ceil[x1]
   imax = floor[x2]
   if imin > imax
      return undef
   return [imin, imax]
}

/** Tests if this interval contains an integer.  If not, it return undef.  If
    it does contain an integer, this narrows the interval to its integer
    bounds and returns an array [imin, imax] with the integer bounds. */
narrowIntegerUpper[cx, x2] :=
{
   if isInteger[cx]
      imin = cx + 1
   else
      imin = ceil[cx]
   
   imax = floor[x2]
   if imin > imax
      return undef
   return [imin, imax]
}

splitInt[x1, x2, intsols] :=
{
   cx = (x2-x1) / 2 + x1  // Done this way so dates work
   
   if intsols
   {
      [x1, cx1] = narrowInteger[x1, cx]
      [cx2, x2] = narrowIntegerUpper[cx, x2]
      return [x1, cx1, cx2, x2]
   }
   
   return [x1, cx, cx, x2]
}

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 = 10000"
println[]
println[f]
// This is a search space of about 64 quintillion integer points!
a = integerSolve[f, -2000000.1, 2000000, -2000000, 2000003, -2000000, 2000009]
for c = lexicalSort[a]
   println[c]
println[length[a] + " solutions"]

/*
// Integer factoring (inefficient but a good test.  Could
// possibly be more efficient with more constraints/"witnesses"?)
println[]
f = "x y = 1009 * 2003 and x <= y"
a = integerSolve[f, 1, billion + 10.7, 1, billion + 10.5]
println[formatTable[sort[a, byColumn[0]]]]
*/