Download or view IntervalSolve.frink in plain text format
/** 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]
Download or view IntervalSolve.frink in plain text format
This is a program written in the programming language Frink.
For more information, view the Frink
Documentation or see More Sample Frink Programs.
Alan Eliasen, eliasen@mindspring.com