Select Git revision
config.php
-
Kathrin Osswald authoredKathrin Osswald authored
intervals.mac 33.11 KiB
/* Author Chris Sangwin
University of Edinburgh
Copyright (C) 2020 Chris Sangwin
This program is free software: you can redistribute it or modify
it under the terms of the GNU General Public License version two.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for details.
You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>. */
/********************************************************************/
/* A package for manipulating intervals in Maxima. */
/* Based on code by Matthew James Read, 2012. */
/* Re-written, May 2020. Chris Sangwin, <C.J.Sangwin@ed.ac.uk> */
/* */
/* V1.0 May 2020 */
/* */
/********************************************************************/
/* Deal with unions. */
unionp(ex) := if safe_op(ex)="%union" or safe_op(ex)="union" then true else false;
intersectionp(ex) := if safe_op(ex)="%intersection" then true else false;
/* Define simple intervals. */
/* Defines the check functions for when intervals are entered: */
cc_num(x,y) := block([Ans],
Ans: 'cc(x,y), /* Makes Ans equal to the original interval. Note the ' to stop evaluation or else it would create an infinite loop. */
if not ev(real_numberp(x), simp) then /* Checks x is a real number. */
error("intervals: ",x," should be a real number"),
if not ev(real_numberp(y), simp) then /* Checks y is a real number. */
error("intervals: ",y," should be a real number"),
if y<x then Ans:{}, /* Our interval is the empty set if y<x. */
if x=y then Ans:{x}, /* Simply the set {x} is x=y. */
Ans
)$
oo_num(x,y) := block([Ans],
Ans: 'oo(x,y),
if ev(not real_numberp(x) and not(x=inf or x=-inf ), simp) then
error("intervals: ",x," should be a real number"),
if ev(not real_numberp(y) and not(y=inf or y=-inf ), simp) then
error("intervals: ",y," should be a real number"),
if y<x then Ans:{},
if x=y then Ans:{},
Ans
)$
co_num(x,y) := block([Ans],
Ans: 'co(x,y),
if ev(not real_numberp(x), simp) then
error("intervals: ",x," should be a real number"),
if ev((not real_numberp(y) and not(y=inf or y=-inf)), simp) then
error("intervals: ",y," should be a real number"),
if y<x then Ans:{},
if x=y then Ans:{},
Ans
)$
oc_num(x,y) := block([Ans],
Ans: 'oc(x,y),
if ev(not real_numberp(x) and not(x=inf or x=-inf), simp) then error("intervals: ",x," should be a real number"),
if ev(not real_numberp(y), simp) then
error("intervals: ",y," should be a real number"),
if y<x then Ans:{},
if x=y then Ans:{},
Ans
)$
/* Validate student's input. */
/* Return a list of errors for a single connected component. */
interval_validate_single_interval(ex) := block([ret, iop, il, ir],
ret:"",
if trivialintervalp(ex) then return(""),
if not(intervalp(ex)) then
return(StackAddFeedback("", "Interval_notinterval", stack_disp(ex, "i"))),
if not(is(length(args(ex))=2)) then
/* The tex functions only cope with two arguments, so we have to use a string here! */
return(StackAddFeedback("", "Interval_wrongnumargs", stack_disp(string(ex), "i"))),
iop:op(ex),
il:first(args(ex)),
ir:second(args(ex)),
if real_numberp(il) and real_numberp(ir) and is(ir<il) then
ret:StackAddFeedback(ret, "Interval_backwards", stack_disp(ex, "i"), stack_disp(apply(iop,[ir, il]), "i")),
return(ret)
)$
/* Validate a realset, mostly for student feedback, so no errors thrown. */
interval_validate_realset(ex) := block(
if trivialintervalp(ex) then return(""),
if setp(ex) then return(""),
if intervalp(ex) then return(interval_validate_single_interval(ex)),
if safe_op(ex)="%union" then return(apply(sconcat, maplist(interval_validate_realset, args(ex)))),
if safe_op(ex)="%intersection" then return(apply(sconcat, maplist(interval_validate_realset, args(ex)))),
return(StackAddFeedback("", "Interval_illegal_entries", stack_disp(ex, "i")))
)$
cc_interval_tex(ex) := block([a, b],
a:first(args(ex)),
b:second(args(ex)),
concat("\\left[ ", tex1(a),",\\, ",tex1(b), "\\right]")
)$
texput(cc, cc_interval_tex)$
/* Note, the mismatching square brackets play havoc with the PHP interface. */
co_interval_tex(ex) := block([a, b],
a:first(args(ex)),
b:second(args(ex)),
/*concat("\\left[ ", tex1(a),",\\, ",tex1(b), "\\right)")*/
concat("!LEFTSQ! ", tex1(a),",\\, ",tex1(b), "!RIGHTR!")
)$
texput(co, co_interval_tex)$
oc_interval_tex(ex) := block([a, b],
a:first(args(ex)),
b:second(args(ex)),
/*concat("\\left( ", tex1(a),",\\, ",tex1(b), "\\right]")*/
concat("!LEFTR! ", tex1(a),",\\, ",tex1(b), "!RIGHTSQ!")
)$
texput(oc, oc_interval_tex)$
oo_interval_tex(ex) := block([a, b],
a:first(args(ex)),
b:second(args(ex)),
concat("\\left( ", tex1(a),",\\, ",tex1(b), "\\right)")
)$
texput(oo, oo_interval_tex)$
realset_tex(ex) := block([a, b, c],
if not(length(args(ex))=2) then error("realset: this function must have two arguments."), a:first(args(ex)),
b:second(args(ex)),
c:ev(interval_complement(b), simp),
if safe_setp(c) then
concat("{", tex1(a), " \\not\\in {",tex1(c), "}}")
else
concat("{", tex1(a), " \\in {",tex1(b), "}}")
)$
texput(realset, realset_tex)$
/* Returns True if p is an element of A. False, otherwise: */
inintervalp(p, A) := block ([Ans, Args, x, y, Atemp, cc:cc(0,1), oo:oo(0,1), co:co(0,1), oc:oc(0,1), i:1, j:1, n],
cc:op(cc), oo:op(oo), co:op(co), oc:op(oc),
Ans:false,
if not ev(real_numberp(p), simp) then
error("intervals: ",p," should be a real number"),
if atom(A) then Ans:false
elseif op(A)=set then
(
Atemp:listify(A),
n:length(Atemp),
while i<(n+1) do
(
if p=Atemp[i] then Ans:true,
i:i+1 )
)
elseif not( op(A)="[" ) then
(
Args:args(A),
x:first(Args),
y:last(Args),
if op(A)=cc then
(
if (p>=x and p<=y) then Ans:true
),
if op(A)=oo then
(
if (p>x and p<y) then Ans:true
),
if op(A)=co then
(
if (p>=x and p<y) then Ans:true
),
if op(A)=oc then
(
if (p>x and p<=y) then Ans:true
)
)
elseif op(A)="[" then
(
n:length(A),
while j<n+1 do
(
Atemp:A[j],
Ans:inintervalp(p,Atemp),
if Ans=false then j:j+1 else j:n+1
)
)
else error("intervals: the interval, ",A,", is not of a recognised form"),
Ans
)$
intervalp(X) := if (safe_op(X)="cc" or safe_op(X)="oo" or safe_op(X)="oc" or safe_op(X)="co") then true else false$
realsetp(ex) := block(
if is(ex=all) then return(true),
if is(ex=none) then return(true), if atom(ex) then return(false),
if safe_setp(ex) then return(all_listp(real_numberp, args(ex))),
if intervalp(ex) then return(all_listp(extended_real_numberp, args(ex))),
if op(ex)=%union then return(all_listp(realsetp, args(ex))),
if op(ex)=%intersection then return(all_listp(realsetp, args(ex))),
return(false)
)$
/* Does not require all numbers to be actual real numbers. */
realset_soft_p(ex) := block(
if is(ex=all) then return(true),
if is(ex=none) then return(true),
if atom(ex) then return(false),
if safe_setp(ex) then return(true),
if intervalp(ex) then return(true),
if op(ex)=realset then return(true),
if op(ex)=%union then return(all_listp(realset_soft_p, args(ex))),
if op(ex)=%intersection then return(all_listp(realset_soft_p, args(ex))),
return(false)
)$
/* Only looks at the very top level, used for validation */
realset_surface_p(ex) := block(
if is(ex=all) then return(true),
if is(ex=none) then return(true),
if atom(ex) then return(false),
if safe_setp(ex) then return(true),
if intervalp(ex) then return(true),
if op(ex)=%union then return(true),
if op(ex)=%intersection then return(true),
return(false)
)$
/* Make a real set, taking edge cases into account. This is also a top level function to convert true/false into all/none. */
realsetmake(v, ex) := block(
if is(ex=false) then return(none),
if is(ex={}) then return(none),
if is(ex=%union()) then return(none),
if is(ex=%intersection()) then return(none),
if is(ex=true) then return(all),
if is(ex=all) or is(ex=none) or is(ex=unknown) then return(ex),
if atom(ex) then return(ex),
if is(safe_op(ex)="realset") then return(ex),
if not(realset_soft_p(ex)) then error("realsetmake: second argument must appear to be a real set."),
return(realset(v, ex))
)$
/* Predicate to remove trivial cases like oo(a,a) and co(-inf, -inf). */
trivialintervalp(ex) := block(
if is(ex=all) or is(ex=none) then return(true),
if safe_setp(ex) and ex={} then return(true),
if not(intervalp(ex)) then return(false),
if safe_op(ex)="oo" and first(ex)=second(ex) then return(true),
if safe_op(ex)="oc" and first(ex)=second(ex) then return(true),
if safe_op(ex)="co" and first(ex)=second(ex) then return(true),
if first(ex)=inf then return(true),
if second(ex)=-inf then return(true),
return(false)
)$
/* Return the number of separate connected components. */
interval_count_components(ex) := block(
if not(realsetp(ex)) then error("interval_count_components"),
if ex=all then return(1),
if trivialintervalp(ex) then return(0),
if intervalp(ex) then return(1),
if setp(ex) then return(cardinality(ex)),
ev(apply("+", map(interval_count_components, args(ex))), simp)
)$
interval_simple_union(X,Y) := block([A:X, B:Y, Ans, x1, x2, y1, y2, Args1, Args2, Aset, swap:false, setAns:[], cc:cc(0,1), oo:oo(0,1), co:co(0,1), oc:oc(0,1), i:1, j:1, n],
cc:op(cc), oo:op(oo), co:op(co), oc:op(oc),
if A=all then return(all),
if B=all then return(all),
if A=none or A={} then return(B),
if B=none or B={} then return(A),
if atom(A) then error("interval_simple_union: invalid first argument"),
if atom(B) then error("interval_simple_union: invalid second argument"),
if safe_setp(A) then (
if safe_setp(B) then
Ans:union(A,B)
else (
Args1:args(B),
x1:first(Args1),
y1:last(Args1),
Aset:listify(A),
n:length(Aset),
while i<(n+1) do (
if (Aset[i]<x1 or Aset[i]>y1) then
setAns:cons(Aset[i],setAns)
elseif Aset[i]=x1 then (
if op(B)=oc then B:cc(x1,y1),
if op(B)=oo then B:co(x1,y1)
)
elseif Aset[i]=y1 then (
if op(B)=co then B:cc(x1,y1),
if op(B)=oo then B:oc(x1,y1)
),
i:i+1
),
if length(setAns)>0 then (setAns:setify(setAns), Ans: [B,setAns] ) else Ans:B
)
)
elseif safe_setp(B) then (
Args1:args(A),
x1:first(Args1), y1:last(Args1),
Aset:listify(B),
n:length(Aset),
while i<(n+1) do (
if (Aset[i]<x1 or Aset[i]>y1) then
setAns:cons(Aset[i],setAns)
elseif Aset[i]=x1 then (
if op(A)=oc then A:cc(x1,y1),
if op(A)=oo then A:co(x1,y1)
)
elseif Aset[i]=y1 then (
if op(A)=co then A:cc(x1,y1),
if op(A)=oo then A:oc(x1,y1)
),
i:i+1
),
if length(setAns)>0 then (setAns:setify(setAns), Ans: [A,setAns] ) else Ans:A
),
if ( not atom(A) and not atom(B) ) then (
Args1:args(A),
Args2:args(B),
if not(atom(A) or safe_setp(A) or atom(B) or safe_setp(B)) then (
if first(Args1)<first(Args2) then
swap:false,
if first(Args1)=first(Args2) then (
if ( op(A)=co or op(A)=cc ) then
swap:false
elseif ( op(B)=co or op(B)=cc ) then
swap:true
else swap:false ),
if first(Args1)>first(Args2) then swap:true,
if swap=false then (
x1:first(Args1),
y1:last(Args1),
x2:first(Args2),
y2:last(Args2)
) else (
Atemp:A,
A:B,
B:Atemp,
x2:first(Args1),
y2:last(Args1),
x1:first(Args2),
y1:last(Args2)
),
if x2>y1 then
Ans:[A,B],
if (x2<y1 and y2>y1) then (
if (op(A)=cc or op(A)=co) then (
if (op(B)=oc or op(B)=cc) then
Ans:cc(x1,y2)
elseif (op(B)=oo or op(B)=co) then
Ans:co(x1,y2)
)
elseif (op(A)=oc or op(A)=oo) then (
if (op(B)=oc or op(B)=cc) then
Ans:oc(x1,y2)
elseif (op(B)=oo or op(B)=co) then
Ans:oo(x1,y2)
)
),
if (x2<y1 and y2=y1) then (
if (op(B)=oc or op(B)=cc) then
Ans:interval_simple_union( A , {y2} )
else
Ans:A
),
if (x2<y1 and y2<y1) then
Ans:A,
if x2=y1 then (
if ( (op(A)=co or op(A)=oo) and (op(B)=oo or op(B)=oc) ) then
Ans:[A,B]
else (
if (op(A)=cc or op(A)=co) then (
if (op(B)=oc or op(B)=cc) then
Ans:cc(x1,y2)
elseif (op(B)=oo or op(B)=co) then
Ans:co(x1,y2)
)
elseif (op(A)=oc or op(A)=oo) then (
if (op(B)=oc or op(B)=cc) then
Ans:oc(x1,y2)
elseif (op(B)=oo or op(B)=co) then
Ans:oo(x1,y2)
)
)
)
)
),
Ans
)$
/* Finds the intersection of two "simple" real sets. */
interval_simple_intersect(X,Y) := block([A:X, B:Y, Ans, x1, x2, y1, y2, Args1, Args2, Aset, swap:false, lopen:false, ropen:false, setAns:[], cc:cc(0,1), oo:oo(0,1), co:co(0,1), oc:oc(0,1), i:1, n],
cc:op(cc), oo:op(oo), co:op(co), oc:op(oc),
if not(realsetp(X)) then error("interval_simple_intersect expects its first argument to be a real set."),
if not(realsetp(Y)) then error("interval_simple_intersect expects its second argument to be a real set."),
if safe_setp(A) and safe_setp(B) then return(intersect(A,B)),
/* A & B are not both sets. */
if safe_setp(B) then (
A:Y,
B:X
),
if safe_setp(A) then (
Args1:args(B),
x1:first(Args1), y1:last(Args1),
Aset:listify(A),
n:length(Aset),
while i<(n+1) do (
if inintervalp(Aset[i],B) then setAns:cons(Aset[i],setAns),
i:i+1
),
if length(setAns)>0 then (
setAns:setify(setAns),
Ans:setAns
) else (
Ans:{}
),
return(Ans)
),
/* At this point we have both A & B not sets. */
if not(intervalp(A) and intervalp(B)) then error("interval_simple_intersect expects its arguments to be sets or simple intervals."),
Args1:args(A),
Args2:args(B),
if first(Args1)<first(Args2) then
swap:false,
if first(Args1)=first(Args2) then (
if (op(A)=co or op(A)=cc) then (
swap:false
) elseif (op(B)=co or op(B)=cc ) then (
swap:true
) else (
swap:false
)
),
if is(first(Args1)>first(Args2)) then (
swap:true
),
if swap=false then (
x1:first(Args1),
y1:last(Args1),
x2:first(Args2),
y2:last(Args2)
) else (
Atemp:A,
A:B,
B:Atemp,
x2:first(Args1),
y2:last(Args1),
x1:first(Args2),
y1:last(Args2)
),
if x2>y1 then (
Ans:{}
),
if (x2<y1 and y2>y1) then (
if (op(A)=cc or op(A)=oc) then (
if (op(B)=cc or op(B)=co) then
Ans:cc(x2,y1)
elseif (op(B)=oo or op(B)=oc) then
Ans:oc(x2,y1)
) elseif (op(A)=co or op(A)=oo) then (
if (op(B)=co or op(B)=cc) then Ans:co(x2,y1)
elseif (op(B)=oo or op(B)=oc) then (
Ans:oo(x2,y1)
)
)
),
if (x2<y1 and y2<y1) then
Ans:B,
if (x2<y1 and y2=y1) then (
if (op(B)=oc or op(B)=oo) then lopen:true,
if (op(B)=oo or op(B)=co or op(A)=oo or op(A)=co) then ropen:true,
if (lopen and ropen) then Ans:oo(x2,y1),
if (lopen and not ropen) then Ans:oc(x2,y1),
if (not lopen and ropen) then Ans:co(x2,y1),
if (not lopen and not ropen) then Ans:cc(x2,y1)
),
if x2=y1 then (
if ((op(A)=cc or op(A)=oc) and (op(B)=co or op(B)=cc)) then
Ans:{x2}
else
Ans:{}
),
Ans
)$
interval_disjointp(A, B) := if interval_simple_intersect(A, B)={} then true else false$
/* Is the ex1 contained within the real set ex2? */
interval_subsetp(ex1, ex2) := block(
if not(realsetp(ex1)) then error("interval_subsetp expects its first argument to be a real set."),
if not(realsetp(ex2)) then error("interval_subsetp expects its second argument to be a real set."),
if interval_intersect(ex1, ex2) = ex1 then true else false
)$
/* Is the simple interval ex a explicitly a subinterval of EX? */
interval_containsp(ex, EX) := block(
if not(intervalp(ex)) then error("interval_containsp expects its first argument to be a simple interval."),
if not(realsetp(EX)) then error("interval_containsp expects its second argument to be a real set."),
if is(ex=EX) then return(true),
if not(safe_op(EX)="%union" or safe_op(EX)="%intersection") then return(false),
if elementp(ex,setify(args(EX))) then return(true),
return(false)
)$
/* Top level intersection function which takes real sets, such as %unions. */
interval_intersect(X,Y) := block([A, B, Ans:[], temp, m, n, i:1, j:1],
A:X,
B:Y,
if safe_op(A)="%intersection" then A:interval_intersect_list(args(A)),
if safe_op(B)="%intersection" then B:interval_intersect_list(args(B)),
if is(A=all) then return(B),
if is(B=all) then return(A),
if atom(A) then return({}),
if atom(B) then return({}),
if is(A={}) then return({}),
if is(B={}) then return({}),
if op(A)=%union then A:args(A),
if op(B)=%union then B:args(B),
if not(listp(A)) and not(listp(B)) then return(interval_simple_intersect(A,B)),
/* Ensure we have lists to deal with, by making them lists of one element if needed. */
if not(listp(A)) then (temp:[], A:cons(A,temp) ),
if not(listp(B)) then (temp:[], B:cons(B,temp) ),
m:length(A),
n:length(B),
if (m=1 and n=1) then ( A:A[1],
B:B[1],
return(interval_simple_intersect(A,B))
) else (
while i<m+1 do (
while j<n+1 do (
temp:interval_simple_intersect(A[i], B[j]),
if not atom(temp) then (
Ans:append(Ans, [temp])
),
j:j+1
),
j:1,
i:i+1
)
),
if listp(Ans) then (
if length(Ans)=1 then Ans:Ans[1],
if length(Ans)=0 then Ans:{}
),
interval_tidy(Ans)
)$
/* Given a *list* of intervals, returns the intersection of all of them. */
interval_intersect_list(X) :=
block ( [A:X, Ans, n, i, simp],
simp:true,
if X=[] then return({}),
n:length(A),
if n=1 then return(first(A)),
Ans:A[1],
i:2,
while i<n+1 do
(
Ans:interval_intersect(Ans, A[i]),
i:i+1
),
Ans
);
interval_intersect_nary([X]) := interval_intersect_list(X)$
/* Given intervals, returns the same intervals but in ascending order of the first element in the interval. */
interval_sort(X) := block([A:X, Ans:[], x, n, i],
if safe_op(X) = "%union" then A:args(X),
n:length(A),
while n>0 do
(
x:A[1],
i:2,
while i<n+1 do block(
if is(first(A[i]) < first(x)) then x:A[i],
i:ev(i+1,simp)
),
Ans:append(Ans,[x]),
A:delete(x, A, 1),
n:ev(n-1, simp)
),
/* %union does things to its arguments like moving -inf to the right with simp:true. */
/* Return a list to avoid killing the order here. */
Ans
);
/* Given a union of disjoint intervals,
checks whether any intervals are connected, and if so, joins them up and returns the ammended union. */
interval_connect(X) := block([Ans, n, x, y, i:1],
if not(op(X)=%union or listp(X)) then error("interval_connect requires a %union or list of intervals"),
Ans:args(X),
n:length(Ans), while i<n do
(
i:ev(i,simp),
if last( Ans[i] ) >= first( Ans[ev(i+1, simp)] ) then
(
x:interval_simple_union( Ans[i], Ans[ev(i+1, simp)] ),
if ( not op(x) = "[" ) then
(
Ans:delete( Ans[ev(i+1, simp)], Ans, 1 ),
Ans:delete( Ans[i], Ans, 1 ),
Ans:append( Ans, [x] ),
i:ev(i-1, simp),
n:ev(n-1, simp)
)
),
i:i+1
),
if length(Ans) = 1 then return(Ans[1]),
Ans:apply(%union, Ans),
Ans
);
/* Given a union of disjoint sets, returns the "canonical form" of this union: */
interval_tidy(X) := block([A, Ans:[], n, setpart:{}, x, y, i:1],
if X=all then return(all),
if atom(X) then return(Ans:phi),
if listp(X) then X:apply(%union, X),
X:ev(X, %intersection=interval_intersect_nary),
if not(op(X)=%union or listp(X)) then (
Ans:X
) else (
A:args(X),
i:1,
n:length(A),
while i<ev(n+1, simp) do (
i:ev(i,simp),
if safe_setp(A[i]) then (
setpart:union(setpart, A[i] ),
A:delete( A[i], A, 1 ),
i:ev(i-1, simp),
n:ev(n-1, simp)
) else if trivialintervalp(A[i]) then (
A:delete( A[i], A, 1 ),
i:ev(i-1, simp),
n:ev(n-1, simp)
),
i:ev(i+1, simp)
),
A:interval_sort(A),
if is(length(A)>1) then
A:interval_connect(A),
if length(setpart)>0 then A:append( args(A), [setpart] ),
if is(A=[]) then
A:{}
elseif is(length(A)=1) then
A:first(A),
Ans:A
),
if Ans=oo(-inf,inf) then return(all),
Ans
)$
interval_complement_order_points(X):=
block( [A:X, Ans:[], setpart, n, i:1],
A:interval_tidy(A),
if safe_setp(last(A)) then (
setpart:listify(last(A)),
A:delete(last(A), A, 1),
n:length(A) + length(setpart),
while i<n+1 do
(
if length(setpart)>0 then
(
if length(A)=0 then
(
Ans:append( Ans, [ { setpart[1] } ] ),
setpart:delete( setpart[1], setpart, 1 )
)
else
(
if setpart[1] < first( A[1] ) then
(
Ans:append( Ans, [ { setpart[1] } ] ),
setpart:delete( setpart[1], setpart, 1 )
)
else
(
Ans:append( Ans, [ A[1] ] ),
A:delete( A[1], A, 1 )
)
)
),
i:i+1
)
)
else Ans:A,
Ans
)$
/* Return the set complement of a real set. */
interval_complement(X):= block([A:X, Ans:[], x, y, cc:cc(0,1), oo:oo(0,1), co:co(0,1), oc:oc(0,1), n, i:1],
cc:op(cc), oo:op(oo), co:op(co), oc:op(oc),
if atom(A) then return(oo(-inf,inf)),
if not (op(A) = "[" or op(A)=%union) then (
if safe_setp(A) then Ans:interval_set_complement(A)
elseif intervalp(A) then (
if op(A)=co then
(
Ans:append( Ans, [ oo(-inf, first(A) ) ] ),
Ans:append( Ans, [ co( last(A), inf) ] )
),
if op(A)=cc then
(
Ans:append( Ans, [ oo(-inf, first(A) ) ] ),
Ans:append( Ans, [ oo( last(A), inf) ] )
),
if op(A)=oc then
(
Ans:append( Ans, [ oc(-inf, first(A) ) ] ),
Ans:append( Ans, [ oo( last(A), inf) ] )
),
if op(A)=oo then
(
Ans:append( Ans, [ oc(-inf, first(A) ) ] ),
Ans:append( Ans, [ co( last(A), inf) ] )
)
)
) else (
A:interval_complement_order_points(A),
A:args(A),
/* Just use DeMorgan's laws. */
Ans:ev(interval_intersect_list(maplist(lambda([ex2], interval_tidy(interval_complement(ex2))), A)), simp),
if listp(Ans) and length(Ans)=1 then
Ans:Ans[1]
), if listp(Ans) then
Ans:apply(%union, Ans),
Ans
)$
/* Take a set of real numbers, and return the %union of intervals not containing these numbers. */
interval_set_complement(X):= block([A:X, Ans:[], temp, n, i:1],
if not(setp(X)) then error("interval_set_complement requires a set."),
A:listify(A),
n:length(A),
temp:oo(-inf, A[1]),
Ans:[temp],
while i<n do (
temp:oo( A[i], A[i+1] ),
temp:[temp],
Ans:append(Ans, temp),
i:i+1
),
temp:oo(A[n], inf),
temp:[temp],
Ans:append(Ans, temp),
apply(%union, Ans)
)$
/* Turns a single variable system over the reals into a set of real numbers, together with insoluable bits (if any). */
single_variable_solver_real(ex) := block([v, rs1, rs2],
if is(ex=false) then return(none),
if is(ex=true) then return(all),
if atom(ex) then return(ex),
v:listofvars(ex),
if is(length(v)=0) then block
(
if is(ratsimp(lhs(ex)-rhs(ex))=0) then
ex:all
else
ex:none
),
if not(length(v)=1) then return(ex),
v:first(v),
ex:abs_replace_eq(ex),
ex:stack_noteq_single_remove(ex),
ex:subst("%and", "nounand", ex),
ex:subst("%or", "nounor", ex),
/* %not is not an infix operator... */
ex:subst(%not, "not", ex),
ex:subst(%not, "nounnot", ex),
ex:subst("%and", "and", ex),
ex:subst("%or", "or", ex),
/* Notes,
(1) assume_pos automatically removes terms like v>=0 in the simplifier.
(2) we do need simplification here to reduce execution time.
*/
if assume_pos then
ex:block([assume_pos:false], ev(single_variable_solver_real_rec(ex %and (v>=0), v), simp))
else
ex:ev(single_variable_solver_real_rec(ex, v), simp),
if ((safe_op(ex)="[" or safe_op(ex)="%union") and is(length(args(ex))=1)) then ex:first(ex),
if is(ex={}) then return(none),
if is(ex={v}) then return(all),
if logic_edgep(ex) then return(ex),
if is(equal(ex,oo(-inf,inf))) then return(all),
rs1:ex,
rs2:false, if safe_op(ex)="%or" then block(
rs1:ev(sublist(args(ex), realset_soft_p), simp),
rs2:ev(sublist(args(ex), lambda([ex2], not realset_soft_p(ex2))), simp),
if is(length(rs1)=1) then rs1:first(rs1),
if is(rs1=none) then
ex:apply("%or", rs2)
else if is(rs1=all) then
ex:all
else
ex:(if realset_soft_p(rs1) then realsetmake(v, rs1) else rs1) %or apply("%or", rs2)
),
if safe_op(ex)="%union" or safe_setp(ex) then
ex:realsetmake(v, ex),
return(ex)
)$
single_variable_solver_real_rec(ex, v) := block([r0, r1, r2],
if atom(ex) then return(ex),
if intervalp(ex) then return(ex),
/* Equations should look real. */
if not(freeof(%i,ex)) then return(ex),
if equationp(ex) then ex:subst("%or", "nounor", pm_replace(ex)),
if equationp(ex) then return(ev(equation_to_intervals(ex, v), simp)),
if linear_inequalityp(ex) then return(ev(linear_inequality_to_interval(ex), simp)),
/* Possible recursion from here. */
if inequalityp(ex) then ex:ev(inequality_factor_solve(ex), simp),
if safe_op(ex)="%or" or safe_op(ex)="%and" then block(
r0:maplist(lambda([ex2], single_variable_solver_real_rec(ex2, v)), args(ex)),
r1:ev(sublist(r0, realset_soft_p), simp),
r2:ev(sublist(r0, lambda([ex2], not(realset_soft_p(ex2)))), simp)
),
if safe_op(ex)="%or" then return(ev(apply("%or", append([interval_tidy(r1)], r2)), simp)),
if safe_op(ex)="%and" then return(ev(apply("%and", append([interval_intersect_list(r1)], r2)), simp)),
return(ex)
)$
equation_to_intervals(ex, v) := block([sol0, sol1, sol2],
sol0:radcan(solve(ex, v)),
if sol0=[] then return({}),
if logic_edgep(sol0) then return(sol0),
/* We need the "freeof" clause to catch rearrangements of equations. */
sol1:sublist(sol0, lambda([ex2], is(lhs(ex2)=v) and freeof(v, rhs(ex2)))),
sol2:sublist(sol0, lambda([ex2], not(is(lhs(ex2)=v) and freeof(v, rhs(ex2))))),
sol1:maplist(rhs,sol1),
sol1:flatten(setify(sol1)),
if is(length(sol2)=1) then
sol2:first(sol2)
else
sol2:apply("%or", sol2),
if emptyp(sol1) then
return(sol2),
return(sol1 %or sol2)
)$
/* Calculate the natural domain of a single-variable term. */
natural_domain(ex) := block([v, ex2],
if atom(ex) then return(all),
v:listofvars(ex),
if is(v=[]) then return(all),
if ev(not(is(length(v)=1)), simp) then return(unknown),
/* We only work over real expressions. */
if not(is(freeof(%i, ex))) then return(unknown),
/* We only calculate domains of some things. */
if not(is(freeof(sum, ex))) then return(unknown),
if not(is(freeof(int, ex))) then return(unknown), v:first(v),
/* Recurse using true/false instead of all/none, then convert. */
ex2:natural_domain_rec(ex),
if realset_soft_p(ex2) then ex2:realsetmake(v, ex2),
ex2
)$
/* Calculate the natural domain of a single-variable term. */
natural_domain_rec(ex) := block([v, ex2],
if atom(ex) then return(all),
v:listofvars(ex),
if is(v=[]) then return(all),
if not(is(length(v)=1)) then return(unknown),
v:first(v),
if safe_op(ex)="sqrt" then
return(single_variable_solver_real(first(args(ex))>=0)),
if safe_op(ex)="ln" or safe_op(ex)="log" or safe_op(ex)="lg" then
return(single_variable_solver_real(first(args(ex))>0)),
if safe_op(ex)="/" then
ex2:[natural_domain_rec(first(args(ex))), single_variable_solver_real((second(args(ex))>0) %or (second(args(ex))<0))]
else
ex2:map(natural_domain_rec, args(ex)),
/* We have to strip of the realset bit before intersecting. */
ex2:map(lambda([ex3], if is(safe_op(ex3)="realset") then second(ex3) else ex3), ex2),
/* Only return a define value if we really have one. */
if any_listp(lambda([ex3], is(ex3=unknown) or not(realset_soft_p(ex3) or is(ex3=true) or is(ex3=false))), ex2) then
ex2:unknown
else
ex2:interval_intersect_list(ex2),
ex2
)$