nns シリーズ更新
・一変数半代数系の高速ソルバー cineqs2 の正式公開.
・solvex から呼ぶ不等式ソルバー sineq も一変数なら高速になりました(cineqs2 の一世代前のアルゴリズムなので少し遅いですが).
load(qepmax)$ showtime:on$ display2d:false$ prtoff():=kill(prt)$ prton():=(prt([x]):=apply(print,x))$ /*********/ /* tools */ /*********/ /* compare length function */ comparelength(f,g):=length(f)<=length(g)$ comparelength2(f,g):=length(g)<=length(f)$ /* cnf dnf */ infix(Or,60,40)$ infix(And,60,40)$ matchdeclare([aa,bb,cc,dd,ee,ff,gg,aa1,bb1,cc1,dd1,ee1,ff1,gg1],true)$ defrule(ru05r, (aa) Or ((bb) And (cc)), ((aa) Or (bb)) And ((aa) Or (cc)))$ defrule(ru08r, ((bb) And (cc)) Or (aa), ((aa) Or (bb)) And ((aa) Or (cc)))$ defrule(ru05a, (aa) And ((bb) Or (cc)), ((aa) And (bb)) Or ((aa) And (cc)))$ defrule(ru08a, ((bb) Or (cc)) And (aa), ((aa) And (bb)) Or ((aa) And (cc)))$ defrule(ru00r, (aa) Or (bb), (aa) %or (bb))$ defrule(ru00a, (aa) And (bb), (aa) %and (bb))$ orand2orand(f):=block([g],g:f, if atom(g) then return(g), if op(g)="%or" then g:xreduce("Or",args(g)), if op(g)="%and" then g:xreduce("And",args(g)),g)$ cnf(f):=applyb2(applyb2(scanmap(orand2orand,f),ru05r,ru08r),ru00r,ru00a)$ dnf(f):=applyb2(applyb2(scanmap(orand2orand,f),ru05a,ru08a),ru00r,ru00a)$ cnfb(f):=applyb2(applyb2(scanmap(orand2orand,f,bottomup),ru05r,ru08r),ru00r,ru00a)$ dnfb(f):=applyb2(applyb2(scanmap(orand2orand,f,bottomup),ru05a,ru08a),ru00r,ru00a)$ /* necessary condition */ defrule(ne2inq, (aa) # (bb), "%or"((aa)<(bb),(bb)<(aa)))$ defrule(ru00n, (aa) # (bb), "%not"((aa)=(bb)))$ deli2(f):=(applyb2(f,ru001,ru002,ru003,ru004))$ ncondeq(f):=block([exs,L,exk], if is(ncondtype=off) then return(f), econd:[],L:true, prt("necessary eq. condition"), if is(ncondtype=x) then ( exs:fullmap(lambda([x],[E,x]), sort(rest(full_listify(powerset(setify(ncondvars)))),comparelength2)), for exk in exs do (q:qe(exk,(f)),prt([L,q]), if is(not(qe([],(L)%implies(q))=true)) then L:(L)%and(q)) ) else ( exs:fullmap(lambda([x],[E,x]), sort(rest(full_listify(powerset(setify(listofvars(f))))),comparelength2)), for exk in exs do (if is(ncondtype=full) then q:qe(exk,(f)) else q:deli2(qe(exk,(f))), prt([L,q]), if is(not(qe([],(L)%implies(q))=true)) then L:(L)%and(q)) ), prt(L), if is(not(ncondsimp=off)) then L:nnscandc(L), /* if is(qe([],(L)%implies(f))=true) and (tdeg(L)=1) then (prt("equivalent condition"),econd:[dnf(L)],return(dnf(L))) else */ return((L)%and(f)) )$ ncond(f):=block([exs,L,exk],L:true, prt("necessary condition"), exs:fullmap(lambda([x],[E,x]), sort(rest(full_listify(powerset(setify(listofvars(f))))),comparelength2)), for exk in exs do (q:qe(exk,f),prt([L,q]), if not(qe([],(L)%implies(q))=true) then L:(L)%and(nnscan(q))), prt(L), return((L)%and(f)) )$ /***************/ /* nns section */ /***************/ /* non-nonsense */ nns(f):=block([F,L,M], if atom(f) then return(f) elseif op(f)="%and" then (F:sort(args(f),comparelength2),L:setify(F), for g in F do (M:disjoin(g,L),prt(M), if p2emt((substpart("%and",M,0))%implies(g))=true then L:M), return(substpart("%and",L,0))) elseif op(f)="%or" then (F:sort(args(f),comparelength2),L:setify(F), for g in F do (M:disjoin(g,L),prt(M), if p2emt((g)%implies(substpart("%or",M,0)))=true then L:M), return(substpart("%or",L,0))) else f)$ /* non-nonsense scan */ nnscan(f):=(prt("nnscanning..."),scanmap(nns,f))$ nnscanb(f):=(prt("bottomup nnscanning..."),scanmap(nns,f,bottomup))$ /* non-nonsense scan with cnf */ nnscanc(f):=(prt("cnf nnscanning..."),scanmap(nns,cnf(f)))$ nnscancb(f):=(prt("cnf bottomup nnscanning..."),scanmap(nns,cnf(f),bottomup))$ /* non-nonsense scan with dnf */ nnscand(f):=(prt("dnf nnscanning..."),scanmap(nns,dnf(f)))$ nnscandb(f):=(prt("dnf bottomup nnscanning..."),scanmap(nns,dnf(f),bottomup))$ nnscand2(f):=(prt("dnf nnscanning..."),nns(map(nns,dnf(f))))$ /* non-nonsense scan with cnf,dnf */ nnscancd(f):=nnscanc(nnscand(f))$ nnscancdb(f):=nnscancb(nnscandb(f))$ nnscandc(f):=nnscand(nnscanc(f))$ nnscandcb(f):=nnscandb(nnscancb(f))$ /* non-nonsense set */ nnss(f):=block([f1,f2,f3,q,nns4], if mode=d then nns4(g):=nnscand(g) else nns4(g):=nns(g), prt("pre-simplified formula"),f1: dnf (qe ([],f) ) , if atom(f1) then return(f1), prt(length(f1)),prt(f1),prt("necessary condition"), if op(f1)="%or" then (f2:[], for k:1 thru length(f1) do (q:ncond(part(f1,k)),prt(k),prt(q),f2:append(f2,[q])), f2:substpart("%or",f2,0)) else (f2:dnf(ncond(f1)),prt(f2)), prt("simplifing..."), if op(f2)="%or" then (f3:[], for k:1 thru length(f2) do (q:qe([],part(f2,k)),prt(k),prt(q),f3:append(f3,[q])), prt("nnscanning..."),f3:substpart("%or",f3,0),prt(f3),return(nns4(f3))) else q:qe([],f2),prt(q),prt("nnscanning..."),return(nns4(q)) )$ /* nnss with double negation */ nnssdn(f):=qe([],%not(nnss(%not(f))))$ /**************/ /* p2t section */ /**************/ defrule(powertk1, (aa)^(bb), makeprtk(aa,num(bb),denom(bb)))$ defrule(powertk2, abs(aa), makeprtk(aa^2,1,2))$ domain:complex$ defrule(powertk2, abs(aa), (aa^2)^(1/2))$ defrule(powertk3, max2(aa,bb), (aa+bb+abs(aa-bb))/2)$ defrule(powertk4, min2(aa,bb), (aa+bb-abs(aa-bb))/2)$ p2e0(f):=block([L,makeprtk,g,q,rt,gg],L:[], makeprtk(x,a,b):=block([rtk], if b=1 then return(x^a) else (rtk:concat(rt,length(L)+1), if oddp(b) then L:append(L,[(rtk^b=x^a)]) else L:append(L,[(rtk^b=x^a)%and(rtk>=0)]), return(rtk))), g:applyb1(f,powertk4,powertk3,powertk2,powertk1), q:makelist([E,concat(rt,k)],k,length(L)), if L=[] then return(f) else (gg:(substpart("%and",L,0))%and(g), ev(qe(q,gg),noeval)) )$ /* 1-1 type */ p2e1(f):=block([L,M,makeprtk,g,q,rt,gg],M:[],L:[], makeprtk(x,a,b):=block([rtk], if b=1 then return(x^a) elseif is(assoc(x^(a/b),M)#false) then return(concat(rt,assoc(x^(a/b),M))) else (M:append(M,[[x^(a/b),length(L)+1]]),rtk:concat(rt,length(L)+1), if oddp(b) then L:append(L,[(rtk^b=x^a)]) else L:append(L,[(rtk^b=x^a)%and(rtk>=0)]), return(rtk))), g:applyb1(f,powertk4,powertk3,powertk2,powertk1), q:makelist([E,concat(rt,k)],k,length(L)), if L=[] then return(f) else (gg:(substpart("%and",L,0))%and(g), ev(qe(q,gg),noeval)) )$ /* base64 ver */ p2e64(f):=block([sL,L,makeprtk,g,q,gg,rt],sL:[],L:[], makeprtk(x,a,b):=block([strg,rtk], if b=1 then return(x^a) else ( strg:ssubst("equal","=",base64(string(x^(a/b)))), strg:ssubst("plus","+",strg), strg:ssubst("slash","/",strg), rtk:eval_string(strg), if not(member(rtk,sL)) then sL:append(sL,[rtk]), if oddp(b) then L:append(L,[(rtk^b=x^a)]) else L:append(L,[(rtk^b=x^a)%and(rtk>=0)]), return(rtk))), g:applyb1(f,powertk4,powertk3,powertk2,powertk1), q:makelist([E,s],s,sL),/* print([g,sL,q]), if is(is(p2etime>0)=true) then p2etime:p2etime+1 else p2etime:1, */ if L=[] then return(f) else (gg:(substpart("%and",L,0))%and(g),prt(sL), subst(makelist(part(sL,k)=concat(rt,k),k,length(sL)), ev(qe(q,gg),noeval)) ) )$ p2e(f):=if is(p2etype=64) then p2e64(f) elseif is(p2etype=0) then p2e0(f) else p2e1(f)$ /* p2e map */ p2em(f):=block([p2etlocal], if atom(f) then return(f) else p2etlocal(g):=block([], if atom(g) then g else if atom(g) then g elseif op(g)="#" then "%not"(p2e(substpart("=",g,0))) elseif (op(g)="=" or op(g)="<" or op(g)=">" or op(g)="<=" or op(g)=">=") then p2e(g) else g), scanmap(p2etlocal,f) /*qe([],scanmap(p2etlocal,f))*/ )$ /* p2e map tarski */ p2emt(f):=qe([],ev(p2em(f),eval))$ p2t(f):=p2emt(f)$ /*******************/ /* nnsolve section */ /*******************/ /* nnsolvexx */ algexact:true$ /* is imaginary */ imgp(x):=not(is(x>=0) or is(x<0)) or (listofvars(x)=[] and imagpart(x)#0)$ /* del formulas containnig imginary unit */ delimg(f):=block([sel3], if atom(f) then return(f), sel3(g):=block([], if atom(g) then g elseif (op(g)="#" or op(g)="=" or op(g)="<" or op(g)=">" or op(g)="<=" or op(g)=">=") and (ratsimp(imagpart(lhs(g)))=0) and (ratsimp(imagpart(rhs(g)))=0) then ratsimp(realpart(g)) elseif (op(g)="#" or op(g)="=" or op(g)="<" or op(g)=">" or op(g)="<=" or op(g)=">=") and (imgp(lhs(g))=true or imgp(rhs(g))=true) then false else g), scanmap(sel3,f))$ Lintersection(l,m):=listify(intersection(setify(l),setify(m)))$ Lminus(l,m):=listify(setdifference(setify(l),setify(m)))$ /* nnsolve: solver for cnj. of eq. */ /* new ver */ delp(L):=block([LL,ek],LL:L, for ek in L do (if member(rhs(ek),%rnum_list) then LL:subst(rhs(ek)=lhs(ek),LL)), LL )$ /* new ver */ nnsolve(f,[vv]):=block([fs,v,g,gg], linsolvewarn:false, if atom(f) then return(f) elseif op(f)="%and" then fs:args(f) elseif op(f)="=" then fs:[f] else return(f), prt("nnsolving..."),prt(fs),%rnum_list:[], v:Lintersection(vv,listofvars(fs)), if v=[] then v:listofvars(fs), g:delimg(algsys(fs,v)),prt([v,g]), if g=[] then (v:listofvars(fs),g:delimg(algsys(fs,v)),prt([v,g])), if g=[] or setify(g)={[false]} then return(false), solvar:v, return(substpart("%or",map(lambda([l],substpart("%and",l,0)),map(delp,g)),0)) )$ defrule(ru000, (aa) # (bb), true)$ defrule(ru001, (aa) > (bb), true)$ defrule(ru002, (aa) >= (bb), true)$ defrule(ru003, (aa) < (bb), true)$ defrule(ru004, (aa) <= (bb), true)$ defrule(ru006, (aa) = (bb), true)$ /* delete ineq. if f has no eq. then return true */ deli1(f):=(applyb2(f,ru000,ru001,ru002,ru003,ru004))$ deli(f):=block([g], if atom(f) then return(f), g:dnf(f), if op(g)="%or" then substpart("%or",delete(true,map(deli1,args(g))),0) else deli1(g))$ /* delete eq. if f has no ineq. then return true */ dele1(f):=(applyb2(f,ru006))$ dele(f):=block([g], if atom(f) then return(f), g:dnf(f), if op(g)="%or" then substpart("%or",delete(true,map(dele1,args(g))),0) else dele1(g))$ mapply(s,t):=if atom(t) then s(t) elseif op(t)="%or" then map(s,t) else s(t)$ /* solution selector for eq.s and ineq.s */ defrule(ru007, (aa) = (bb), if Lintersection(v,listofvars(expand(aa-bb)))=[] then true else (aa)=(bb))$ defrule(ru008, (aa) = (bb), if Lintersection(v,listofvars(expand(aa-bb)))=[] then (aa)=(bb) else true)$ delec1(f):=(applyb1(f,ru007))$ delec(f):=block([g], if atom(f) then return(f), g:dnf(f), if op(g)="%or" then substpart("%or",delete(true,map(delec1,args(g))),0) else delec1(g))$ delev1(f):=(applyb1(f,ru008))$ delev(f):=block([g], if atom(f) then return(f), g:dnf(f), if op(g)="%or" then substpart("%or",delete(true,map(delev1,args(g))),0) else delev1(g))$ /* new ver for ineq. and conj. of eq., separate ng type */ solsel(f,[v]):=block([eqss,eqs,ins,sel2,ng], eqss:deli(f), if v=[] or length(listofvars(eqss))=1 then v:listofvars(eqss), eqs:delec(eqss), prt(eqs), if atom(eqs) then return(f), if op(eqs)="%or" then (prt("I can't solve disj. of eq."),return(f)), ins:(dele(f))%and(delev(eqss)), sel2(g):=block([gg], /*(g)%and(qe(map(lambda([x],[E,x]),listofvars(g)),p2et((g)%and(ins))))*/ if /* ins=true then g elseif */ solvar=[] or atom(g) then (g)%and(ins) else (prt("verifing..."),prt(g), solvar:map(lhs,(if op(g)="%and" then args(g) else [g])), gg:(g)%and(qe(map(lambda([x],[E,x]),solvar),p2t((g)%and(ins)))), ng:(ng)%and(qe([],"%not"(qe(map(lambda([x],[E,x]),solvar),p2t(g))))), prt("sel2 returned..."),prt(gg),gg ) ), qq:apply(nnsolve,flatten([eqs,v])),prt(qq), if ins=true then prt("no ins"), ng:f, qq:mapply(sel2,qq), if is((ng)#false) then (prt("other type..."),prt(ng), /* ncondtype:full,ncondsimp:off, */ inner:true, qq:(qq)%or(ora(apply(nnsolvexx,flatten([ng,solvar])))), /* ncondtype:eq,ncondsimp:nn, */ inner:false), qq /* for nnsolvexx */ /* nnscan(qq) */ /* for stand alone */ )$ /* nnsolvexx: solver for semi-algebraic system */ nnsolvexx0(f,[v]):=block([g,gg], prt("pre-simplified formula"), g:dnf(qe([],p2t(f))),prt(g), if atom(g) then return(g) else (prt("necessary condition"),g:dnf(ncond(g)),prt(length(g)),prt(g)), if op(g)="%or" then (prt("simplifing..."),gg:dnf(map(ncond,g)),prt(length(gg)),prt(gg),nnscan(map(lambda([t],apply(solsel,flatten([t,v]))),gg))) else nnscan(apply(solsel,flatten([g,v]))) )$ /* total deg */ tdeg0(f):=block([tdegx], hipow(expand(subst(map(lambda([s],s=tdegx),listofvars(f)),f)),tdegx) )$ tdeg(f):=block([g],g:expand( lhs(f)-rhs(f)), if atom(g) then tdeg0(g) else apply(max,map(tdeg0,args(g)) ) )$ /* new ver */ nnsolvexx(f,[vv]):=block([t0,g,v,gg,q,dele1q], if not(inner=true) then (thisin:f,t0:elapsed_real_time()), if is(presimp=off) then prt("pre-simplify:off") elseif is(presimp=l) then prt("pre-simplify:lineq with nnscan") else prt("pre-simplify:p2t->qe"), if is(ncondtype=off) then prt("ncond-type:off") elseif is(ncondtype=full) then prt("ncond-type:full formula") elseif is(ncondtype=x) then prt("ncond-type:x") else prt("ncond-type:eq. only"), if is(ncondsimp=off) then prt("ncond-post-simplify:off") else prt("ncond-post-simplify:nnscandc"), if is(sineqsimp=nn) then prt("sineq-simplify:nnscan") else prt("sineq-simplify:lineq"), if is(postsimp=off) then prt("post-simplify:off") else prt("post-simplify:nnscanc->nnscand->nnscan"), prt(""),prt("given formula"),prt(f), if is(presimp=off) then g:dnf(f) elseif is(presimp=l) then (prt("pre-simplified dnf"),g:dnf(applyb1(lineq(f),ne2inq))) else (prt("pre-simplified dnf"),g:dnf(applyb1(qe([],p2t(f)),ne2inq))),prt(g), if atom(g) then return(bra(g)),prt("var for solve"), v:vv,for xxx in Lminus(vv,listofvars(g)) do v:delete(xxx,v), if v=[] then v:listofvars(f),prt(v),ncondvars:v, if atom(g) then return(bra(g)) else (if is(ncondsimp=off) then gg:dnf(mapply(ncondeq,g)) else gg:nnscan(dnf(mapply(ncondeq,g))), prt(length(gg)),prt(gg),prt("solving..."), q:scanmap(lambda([t],apply(fineq,flatten([t,last(v)]))), mapply(lambda([t],apply(solsel,flatten([t,v]))),gg) ,bottomup), prt("solved formula"),prt(q),prt(elapsed_real_time()-t0), collect(ff):=expand(ff), /* collect(ff):=block([fff,vc],fff:expand(f),vc:listofvars(fff), sum(factor(coeff(fff,vc,d))*x^d,d,0,hipow(fff,vc))), */ deleq:dele(q), /* if length(Lminus(listofvars(dele1q),listofvars(facts())))=1 then */ if length(listofvars(deleq))=1 then (thisout:mapply(lambda([g],(cineqs(dele1(g))) %and (deli1(g))),dnf(q)), if is(deleq=dnf(q)) then thisout:(cineqs(thisout)), thisout:bra(l2r(thisout)),thisout) elseif is(postsimp=off) then (thisout:scanmap(collect,l2rx(bra(l2r(q)),last(v))),thisout) else (q:mapply(nnscanc,q),prt("reduced formula"),prt(q), q:mapply(nnscand,q),prt("reduced formula"),prt(q), q:nnscan(q), thisout:scanmap(collect,l2rx(bra(l2r(q)),last(v))),thisout) ) )$ thischeck():=fullall((qex(thisin))%eq(qex(ora(thisout))))$ /* for scanmap */ sineq(ineq,[vv]):=block([p,v,xxx,sineqx,q,r,s,lc,sol,sols,sk,rt,qq,excond], if atom(ineq) then return(ineq), if not(op(ineq)=">" or op(ineq)="<" or op(ineq)=">=" or op(ineq)="<=") then return(ineq), p:ratsimp(lhs(ineq)-rhs(ineq)), if imgp(p)=true then return(ineq), if tdeg(num(factor(p)))<=1 and listofvars(denom(factor(p)))=[] then return(ineq), if length(vv)=1 then (if hipow(p,last(vv))=1 then return(ineq)), v:vv,for xxx in Lminus(vv,listofvars(p)) do v:delete(xxx,v), if v=[] then v:listofvars(p),sineqx:last(v), q:1,r:1,sols:[], sol:delimg(ratsimp(solve(lhs(ineq)-rhs(ineq),sineqx))), for k:1 thru length(sol) do (sk:part(sol,k),if not(sk=false) then (q:q*(lhs(sk)-rhs(sk))^(part(multiplicities,k)), r:r*(lhs(sk)-concat(rt,k))^(part(multiplicities,k)), sols:append(sols,makelist(lhs(sk)=rhs(sk),j,part(multiplicities,k))))), s:factor(quotient(p,q,sineqx)), lc:ratcoef(expand(s),sineqx,hipow(expand(s),sineqx)),/* lc:lc/abs(lc), */ prt("factor type"),prt([p,sineqx,r,s,lc]), if listofvars(q)=[sineqx] then qq:cineq(sols,lc,op(ineq)) else ( if s=lc then qq:qe([],qe([],apply(op(ineq),[r*lc,0]))) elseif ((op(ineq)=">") or (op(ineq)="<")) then (prt("nnsolving..."), qq:(qe([],qe([],apply(op(ineq),[r*lc,0]))) %and nnsolvexx0(s#0))) elseif ((op(ineq)=">=") or (op(ineq)="<=")) then (prt("nnsolving..."), qq:(qe([],qe([],apply(op(ineq),[r*lc,0]))) %or nnsolvexx0(s=0))), prt(qq), /* prt(qe([],cnf(qq))),*/ for k:1 thru length(sol) do (sk:part(sol,k),if not(sk=false) then qq:subst(concat(rt,k)=rhs(sk),qq)), if is(sineqsimp=off) then qq:qq elseif is(sineqsimp=ex) then qq:(map(lambda([t],(t)%and(fullex(t))),qq)) elseif is(sineqsimp=l) then qq:lineq(qq) else qq:nnscan(qq), excond:qe([[E,sineqx]],p2t(qq)),prt(excond), if not(is(excond=true)) then (inner:true, qq:(qq) %or (ora(nnsolvexx((ineq)%and(%not(excond)),sineqx))), inner:false) ), prt("sineq returned..."),prt(qq), qq )$ /* with cineq */ sineq(ineq,[vv]):=block([p,v,xxx,sineqx,q,r,s,lc,sol,sols,sk,rt,qq,excond], if atom(ineq) then return(ineq), if not(op(ineq)=">" or op(ineq)="<" or op(ineq)=">=" or op(ineq)="<=") then return(ineq), p:ratsimp(lhs(ineq)-rhs(ineq)), if imgp(p)=true then return(ineq), if listofvars(p)=[] then (prt("const."),return("%and"(apply(op(ineq),[p,0])))), if (tdeg(num(factor(p)))=1 and listofvars(denom(factor(p)))=[] and length(listofvars(expand(p)))=1) then (prt("linear"),/*********************/ v:first(listofvars(p)),lc:ratcoef(expand(p),v,hipow(expand(p),v)), if lc>0 then return(apply(op(ineq),[v,expand((lc*v-p)/lc)])) else return(apply(op(ineq),[expand((lc*v-p)/lc),v]))), /* if length(vv)=1 then (if hipow(p,last(vv))=1 then return(ineq)), */ v:vv,for xxx in Lminus(vv,listofvars(p)) do v:delete(xxx,v), if v=[] then v:listofvars(p),sineqx:last(v), q:1,r:1,sols:[], sol:delimg(ratsimp(solve(lhs(ineq)-rhs(ineq),sineqx))), for k:1 thru length(sol) do (sk:part(sol,k),if not(sk=false) then (q:q*(lhs(sk)-rhs(sk))^(part(multiplicities,k)), r:r*(lhs(sk)-concat(rt,k))^(part(multiplicities,k)), sols:append(sols,makelist(lhs(sk)=rhs(sk),j,part(multiplicities,k))))), s:factor(quotient(p,q,sineqx)), lc:ratcoef(expand(s),sineqx,hipow(expand(s),sineqx)),/* lc:lc/abs(lc), */ prt("factor type"),prt([p,sineqx,r,s,lc]), if listofvars(q)=[sineqx] then qq:cineq(sols,lc,op(ineq)) /* elseif length(delete(sineqx,listofvars(q)))=1 then qq:tineq(sols,lc,op(ineq)) */ else ( if s=lc then qq:qe([],qe([],apply(op(ineq),[r*lc,0]))) elseif ((op(ineq)=">") or (op(ineq)="<")) then (prt("nnsolving..."), qq:(qe([],qe([],apply(op(ineq),[r*lc,0]))) %and nnsolvexx0(s#0))) elseif ((op(ineq)=">=") or (op(ineq)="<=")) then (prt("nnsolving..."), qq:(qe([],qe([],apply(op(ineq),[r*lc,0]))) %or nnsolvexx0(s=0))), prt(qq), /* prt(qe([],cnf(qq))),*/ for k:1 thru length(sol) do (sk:part(sol,k),if not(sk=false) then qq:subst(concat(rt,k)=rhs(sk),qq)), if is(sineqsimp=off) then qq:qq elseif is(sineqsimp=ex) then qq:(map(lambda([t],(t)%and(fullex(t))),qq)) elseif is(sineqsimp=l) then qq:lineq(qq) else qq:nnscan(qq), excond:qe([[E,sineqx]],p2t(qq)),prt(excond), if not(is(excond=true)) then (inner:true, qq:(qq) %or (ora(nnsolvexx((ineq)%and(%not(excond)),sineqx))), inner:false) ), prt("sineq returned..."),prt(qq), qq )$ defrule(delsame, ((aa) < (bb)) %and ((bb) < (aa)), false)$ defrule(eqsame, ((aa) <= (bb)) %and ((bb) <= (aa)), (bb)=(aa))$ sineq2(ineq,sineqx):=block([p,sol,sols,conds,cfs,cqs,lc,crs,css,cts], /* sineq2(ineq,sineqx):=block([], */ p:(lhs(ineq)-rhs(ineq)),prt(p), sol:delimg(ratsimp(solve(lhs(ineq)-rhs(ineq),sineqx))),prt(sol), sol:makelist([part(sol,k),part(multiplicities,k),qex([[E,sineqx]],part(sol,k))],k,length(sol)),prt(sol), sols:[],for c in sol do if not(last(c)=false) then sols:append(sols,[ [(lhs(first(c))-rhs(first(c)))^(second(c)), (lhs(first(c))-concat(rt,length(sols)+1))^(second(c)), concat(rt,length(sols)+1)=rhs(first(c)), last(c)]]),prt(sols), conds:listify(setify(makelist(last(s),s,sols))),prt(conds), conds:delete(false,map(qex, delete(false,flatten(subst("Or"="[",applyb2(applyb2(xreduce("And",map(lambda([f],(f) Or (qex(%not(f)))),append(conds,[true]))),ru05a,ru08a),ru00a)))))),prt("case conds (conds):"),prt(l2r(conds)), cfs:makelist([c,append([[[1,1],t0=t0]],delete(false,makelist((if (fullall((c)%implies(last(s)))) then [rest(s,-2),part(s,3)]),s,sols)))],c,conds),prt(cfs), cqs:map(lambda([s],([first(s),substpart("*",map(first,last(s)),0),delete(t0=t0,map(last,last(s)))])),cfs),prt(cqs), lc:ratcoef(expand(p),sineqx,hipow(expand(p),sineqx)), crs:map(lambda([s],[first(s), lc*second(second(s)) ,factor(quotient(p,first(second(s)),sineqx))/lc ,last(s)]),cqs),prt(crs), css:map(lambda([s],apply(if op(ineq)="<" or op(ineq)=">" then "%and" else "%or",[(first(s)) %and (subst(last(s), dele( subst(if op(ineq)="<" or op(ineq)=">" then [] else [">"=">=","<"="<="], qex(apply(op(ineq),[second(s),0])))))), if op(ineq)="<" or op(ineq)=">" then ora(nnsolvexx((first(s))%and(third(s)#0))) else ora(nnsolvexx((first(s))%and(third(s)=0))) ])),crs),prt("raw sols:"),prt(css), cts:map(nns,l2rx(map(nnscand,css),sineqx)),prt("reduced sols (cts):"),prt(cts), substpart("%or",cts,0) )$ /* new ver */ /* sineq3(ineq,sineqx):=block([], */ sineq3(ineq,sineqx):=block([p,k,c,sol,sols,conds,cfs,crs,cqs,lc,css,cts], p:(lhs(ineq)-rhs(ineq)),prt(p), sol:delimg(ratsimp(solve(lhs(ineq)-rhs(ineq),sineqx))),prt(sol), sol:makelist([part(sol,k),part(multiplicities,k),qex([[E,sineqx]],part(sol,k))],k,length(sol)),prt(sol), sols:[],for c in sol do if not(last(c)=false) then sols:append(sols,[ [[(lhs(first(c))-rhs(first(c)))^(second(c))], [if second(c)=1 then lhs(first(c))=rhs(first(c)) else makelist(lhs(first(c))=rhs(first(c)),second(c))], last(c)]]),prt(sols), conds:listify(setify(makelist(last(s),s,sols))),prt(conds), conds:delete(false,map(qex, delete(false,flatten(subst("Or"="[",applyb2(applyb2(xreduce("And",map(lambda([f],(f) Or (qex(%not(f)))),append(conds,[true]))),ru05a,ru08a),ru00a)))))),prt("case conds (conds):"),prt(l2r(conds)), cfs:makelist([c,append([append,[[1],[]]],delete(false,makelist((if (fullall((c)%implies(last(s)))) then rest(s,-1)),s,sols)))],c,conds), cqs:map(lambda([s],([first(s), substpart("*",first(apply(map,last(s))),0), last(apply(map,last(s)))])),cfs),prt(cqs), lc:ratcoef(expand(p),sineqx,hipow(expand(p),sineqx)), crs:map(lambda([s],[first(s), factor(quotient(p,second(s),sineqx))/lc, last(s)]),cqs),prt("crs:"),prt(crs), ncondtype_out:ncondtype,postsimp_out:postsimp,ncondtype:off,postsimp:off, css:(apply(append,l2r(map(lambda([s], if op(ineq)="<" or op(ineq)=">" then block([a1],a1:nnssdn(p2t((first(s))%and(second(s)>0))), brand( tineq2(third(s),lc,op(ineq)), sepaxa(if member(sineqx,listofvars(a1)) then bra(a1) else [[a1]] ,sineqx), sineqx ) ) else block([a1],a1:nnss(p2t((first(s))%and(second(s)=0))), append( brand(tineq2(third(s),lc,op(ineq)), [[ first(s) ]], sineqx), sepaxa(if member(sineqx,listofvars(a1)) then bra(a1) else [[a1]] ,sineqx) ) ) ),crs)))), ncondtype:ncondtype_out,postsimp:postsimp_out, prt("raw sols:"),prt(css),css /*cts:map(nns,l2r(map(nnscand,css))),prt("reduced sols (cts):"),prt(cts), substpart("%or",cts,0)*/ )$ /* sineq3(x^4-2*a*x^2<b,x); */ sepaxa(braform,xxx):= map(lambda([l], block([L,M,e],L:[],M:[],for e in l do if member(xxx,listofvars(e)) then L:append(L,[e]) else M:append(M,[e]), [L,M])),map(flatten,braform))$ brand(bra1,bra2,v):=block([s1,s2], sepaxa(apply(append,makelist(makelist(append(s1,s2),s2,bra2),s1,bra1)),v) )$ braor(bra1,bra2):=append(bra1,bra2)$ /* sineq5(ineq,[vv]):=block([], */ sineq5(ineq,[vv]):=block([sineqx,p,k,c,sol,sols,conds,cfs,crs,cqs,lc,css,cts], if length(listofvars(ineq))<=1 then return(map(lambda([s],[s,[]]),cineqs2("%and"(ineq)))), if vv#[] then sineqx:first(vv),p:(lhs(ineq)-rhs(ineq)),prt(p), sol:delimg(ratsimp(solve(lhs(ineq)-rhs(ineq),sineqx))),prt(sol), sol:makelist([part(sol,k),part(multiplicities,k),qex([[E,sineqx]],part(sol,k))],k,length(sol)),prt(sol), sols:[],for c in sol do if not(last(c)=false) then sols:append(sols,[ [[(lhs(first(c))-rhs(first(c)))^(second(c))], [if second(c)=1 then lhs(first(c))=rhs(first(c)) else makelist(lhs(first(c))=rhs(first(c)),second(c))], last(c)]]),prt(sols), conds:listify(setify(makelist(last(s),s,sols))),prt(conds), conds:delete(false,map(qex, delete(false,flatten(subst("Or"="[",applyb2(applyb2(xreduce("And",map(lambda([f],(f) Or (qex(%not(f)))),append(conds,[true]))),ru05a,ru08a),ru00a)))))),prt("case conds (conds):"),prt(l2r(conds)), crs:makelist([c,append([],delete(false,makelist((if (fullall((c)%implies(last(s)))) then second(s)),s,sols)))],c,conds), prt("crs"),prt(crs), lc:ratcoef(expand(p),sineqx,hipow(expand(p),sineqx)), ncondtype_out:ncondtype,postsimp_out:postsimp,ncondtype:off,postsimp:off, css:(apply(append,l2r( map(lambda([s],brand(tineq2(second(s),lc,op(ineq)),[[first(s)]],sineqx)) ,crs)))), ncondtype:ncondtype_out,postsimp:postsimp_out, prt("raw sols:"),prt(css),css )$ /* sineq4(ineq,[vv]):=block([], */ sineq4(ineq,[vv]):=block([p,v,sineqx,hd,k,c,sol,sols,conds,cfs,cqs,lc,crs,css,cts], /* pre-checks */ if atom(ineq) then return([[ineq]]), if not(op(ineq)=">" or op(ineq)="<" or op(ineq)=">=" or op(ineq)="<=") then return([[ineq]]), p:(lhs(ineq)-rhs(ineq)),prt(p), if imgp(p)=true then return([[ineq]]), if listofvars(p)=[] then (prt("const."),return([["%and"(apply(op(ineq),[p,0]))]])), /* def. of sineqx */ v:vv,for xxx in Lminus(vv,listofvars(p)) do v:delete(xxx,v), if v=[] then v:listofvars(p),sineqx:first(v), /* linear */ hd:hipow(expand(p),sineqx),lc:ratcoef(expand(p),sineqx,hd), if hd=1 and listofvars(lc)=[] then (prt("linear"), if lc>0 then return([[apply(op(ineq),[sineqx,expand((lc*sineqx-p)/lc)])]]) else return([[apply(op(ineq),[expand((lc*sineqx-p)/lc),sineqx])]])), /* main */ sol:delimg(ratsimp(solve(lhs(ineq)-rhs(ineq),sineqx))),prt(sol), if listofvars(p)=[sineqx] then return(bra(cineqs(cineq( delete(false,flatten( makelist(if not(part(sol,k)=false) then makelist(lhs(part(sol,k))=rhs(part(sol,k)),part(multiplicities,k)), k,length(sol)))), lc,op(ineq))))), sol:makelist([part(sol,k),part(multiplicities,k),qex([[E,sineqx]],part(sol,k))],k,length(sol)),prt(sol), sols:[],for c in sol do if not(last(c)=false) then sols:append(sols,[ [[(lhs(first(c))-rhs(first(c)))^(second(c))], [if second(c)=1 then lhs(first(c))=rhs(first(c)) else makelist(lhs(first(c))=rhs(first(c)),second(c))], last(c)]]),prt(sols), conds:listify(setify(makelist(last(s),s,sols))),prt(conds), conds:delete(false,map(qex, delete(false,flatten(subst("Or"="[",applyb2(applyb2(xreduce("And",map(lambda([f],(f) Or (qex(%not(f)))),append(conds,[true]))),ru05a,ru08a),ru00a)))))),prt("case conds (conds):"),prt(l2r(conds)), cfs:makelist([c,append([append,[[1],[]]],delete(false,makelist((if (fullall((c)%implies(last(s)))) then rest(s,-1)),s,sols)))],c,conds), cqs:map(lambda([s],([first(s), substpart("*",first(apply(map,last(s))),0), last(apply(map,last(s)))])),cfs),prt(cqs), crs:map(lambda([s],[first(s), factor(quotient(p,second(s),sineqx))/lc, last(s)]),cqs),prt(crs), ncondtype:off,postsimp:off, css:l2r(map(lambda([s], if op(ineq)="<" or op(ineq)=">" then (ora(tineq2(third(s),lc,op(ineq)))) %and ora(nnsolvexx((first(s))%and(second(s)#0))) else ( (first(s)) %and (ora(tineq2(third(s),lc,op(ineq)))) ) %or (ora(nnsolvexx((first(s))%and(second(s)=0))))),crs)), kill(ncondtype,postsimp),prt("raw sols:"),prt(css), map("[",css) /*cts:map(nns,l2r(map(nnscand,css))),prt("reduced sols (cts):"),prt(cts), substpart("%or",cts,0)*/ )$ /* one type */ cineq(sols,lc,gl):=block([L,M], L:sort(sols,lambda([s,t],rhs(s)>=rhs(t))), if (lc>0 and (gl="<" or gl="<=")) or (lc<0 and (gl=">" or gl=">=")) then L:append([first(L)],L), L:makelist( if k=1 then rhs(part(L,k))<lhs(part(L,k)) elseif evenp(k) and k<length(L) then (rhs(part(L,k+1))<lhs(part(L,k+1))) %and (lhs(part(L,k))<rhs(part(L,k))) elseif evenp(k) and k=length(L) then lhs(part(L,k))<rhs(part(L,k)) else false, k,length(L)), if (lc>0 and (gl="<" or gl="<=")) or (lc<0 and (gl=">" or gl=">=")) then L:rest(L), if gl=">=" or gl="<=" then L:subst("<"="<=",L), L:reverse(L),prt(L), apply1(substpart("%or",L,0),delsame,eqsame) )$ cineqsd(l):=block([v,ll,lr], if listofvars(l)=[] then return(first(l)), v:first(listofvars(l)),ll:{},lr:{}, for f in l2r(l) do (if rhs(f)=v then ll:append(ll,{f}) else lr:append(lr,{f})), ll:sort(listify(ll),lambda([s,t],lhs(s)<lhs(t) or (lhs(s)=lhs(t) and op(s)="<=" and op(t)="<"))), lr:sort(listify(lr),lambda([s,t],rhs(s)<rhs(t) or (rhs(s)=rhs(t) and op(s)="<" and op(t)="<="))), prt([ll,lr]), if ll=[] then return(first(lr)) elseif lr=[] then return(last(ll)) elseif lhs(last(ll))<rhs(first(lr)) then (last(ll)) %and (first(lr)) elseif (op(last(ll))="<=" and op(first(lr))="<=" and lhs(last(ll))=rhs(first(lr))) then lhs(first(lr))=rhs(first(lr)) /* (last(ll)) %and (first(lr)) */ else false)$ cineqsc(l):=block([v,ll,lr], if listofvars(l)=[] then return(first(l)), v:first(listofvars(l)),ll:{},lr:{}, for f in l2r(l) do (if rhs(f)=v then ll:append(ll,{f}) else lr:append(lr,{f})), ll:sort(listify(ll),lambda([s,t],lhs(s)<lhs(t) or (lhs(s)=lhs(t) and op(s)="<=" and op(t)="<"))), lr:sort(listify(lr),lambda([s,t],rhs(s)<rhs(t) or (rhs(s)=rhs(t) and op(s)="<" and op(t)="<="))), prt([ll,lr]), if ll=[] then return(last(lr)) elseif lr=[] then return(first(ll)) elseif rhs(last(lr))<lhs(first(ll)) then (last(lr)) %or (first(ll)) elseif (op(last(lr))="<" and op(first(ll))="<" and rhs(last(lr))=lhs(first(ll))) then (last(lr)) %or (first(ll)) /* lhs(last(lr))#rhs(last(lr)) */ else true)$ cineqsdm(L):=(prt("cineqs-dnf returned..."), substpart("%or",map(cineqsd,bra(L)),0))$ cineqscm(L):=(prt("cineqs-cnf returned..."), substpart("%and",map(cineqsc,brac(ora(L))),0))$ defrule(eq2inq, (aa) = (bb), (((aa)<=(bb)) %and ((bb)<=(aa))) )$ cineqs(L):=cineqsdm(cineqscm(apply1(L,ne2inq,eq2inq)))$ cineqs2(F):=block([lv,v,imagunit,L,M,N,k,K,J], lv:listofvars(F),if length(lv)#1 then return([[F]]) else v:first(lv), L:sort(listify(setify(flatten(scanmap( lambda([e], if atom(e) then e elseif member(op(e),["%and","%or"]) then substpart("[",e,0) elseif member(op(e),["<",">","<=",">=","=","#"]) then map(rhs,delete(false,delimg(solve(substpart("=",e,0),v)))) else e), F)))),"<"), if L=[] then return([[is(subst(v=0,F))]]), L:sort(append(L,(append([first(L)-1],L)+append(L,[last(L)+1]))/2),"<"), /* */ M:map(lambda([e],is(ratsimp((subst(v=e,F))))),L),prt([L,M]), /* fpprec:64, M:map(lambda([e],is(rationalize(bfloat(subst(v=e,F))))),L),prt([L,M]), fpprec:16, M:map(lambda([e],is(subst(v=float(e),F))),L),prt([L,M]), */ if length(setify(M))=1 then return([[is(subst(v=0,F))]]), L:append([0,0],L,[0,0]),M:append([true,true],M,[true,true]), K:[],for k:3 thru length(M)-2 do ( if oddp(k) then if part(M,k) then (if not(part(M,k-1)) then K:append(K,[part(L,k-1)<v]) elseif not(part(M,k-2)) then K:append(K,[part(L,k-1)<=v]) else K:K, if not(part(M,k+1)) then K:append(K,[v<part(L,k+1)]) elseif not(part(M,k+2)) then K:append(K,[v<=part(L,k+1)]) else K:K) else K:K else if not(part(M,k-1)) and part(M,k) and not(part(M,k+1)) then K:append(K,[v=part(L,k)]) else K:K), N:[],J:[],for k in K do if rhs(k)=v then J:[k] else (N:append(N,[append(J,[k])]),J:[]), delete([],append(N,[J])) )$ ratprint:false$ /* asksign */ cineqs2f(F):=block([lv,v,imagunit,L,M,N,k,K], lv:listofvars(F),if length(lv)#1 then return(F) else v:first(lv), L:sort(listify(setify(flatten(scanmap( lambda([e], if atom(e) then e elseif member(op(e),["%and","%or"]) then substpart("[",e,0) elseif member(op(e),["<",">","<=",">=","=","#"]) then map(lambda([t],if member(imagunit,listofvars(t)) then [] else rhs(t)),subst(%i=imagunit,solve(substpart("=",e,0),v))) else e), F)))),"<"), if L=[] then return([[is(subst(v=0,F))]]), L:sort(append(L,(append([first(L)-1],L)+append(L,[last(L)+1]))/2),"<"), M:map(lambda([e],is(subst(v=float(e),F))),L),prt([L,M]), if length(setify(M))=1 then return([[is(subst(v=0,F))]]), L:append([0,0],L,[0,0]),M:append([true,true],M,[true,true]), N:[],for k:3 thru length(M)-2 do (K:[], if oddp(k) then if part(M,k) then (if not(part(M,k-1)) then K:append(K,[part(L,k-1)<v]) elseif not(part(M,k-2)) then K:append(K,[part(L,k-1)<=v]) else K:K, if not(part(M,k+1)) then K:append(K,[v<part(L,k+1)]) elseif not(part(M,k+2)) then K:append(K,[v<=part(L,k+1)]) else K:K) else K:K else if not(part(M,k-1)) and part(M,k) and not(part(M,k+1)) then K:append(K,[v=part(L,k)]) else K:K, N:delete([],append(N,[K]))), N)$ /* factor for polynomial inequation */ /* new ver */ fineq(ineq,[vv]):=block([pp,v,xx,yy,dg,lc], if atom(ineq) or listofvars(ineq)=[] then return("%and"(ineq)), if not(op(ineq)=">" or op(ineq)="<" or op(ineq)=">=" or op(ineq)="<=") then return(ineq), pp:ratsimp(lhs(ineq)-rhs(ineq)), if atom(pp) then return(ineq), prt(vv), v:vv,for yy in Lminus(vv,listofvars(pp)) do v:delete(yy,v), if v=[] then v:listofvars(pp),xx:last(v),dg:hipow(pp,xx),prt(xx), if imgp(pp)=true or dg=0 then return(ineq),lc:ratcoef(pp,xx,dg),/*prt(lc),*/ if listofvars(lc)=[] then (prt("case1"),sineq(apply(op(ineq),[lhs(ineq)-rhs(ineq),0]),xx)) elseif is(fineqfull=on) then ((prt("case2-1"),(nnsolvexx(lc=0) %and fineq(apply(op(ineq),[expand(pp-lc*xx^dg),0]),xx))) %or (prt("case2-2"),(fineq(lc>0) %and (sineq(apply(op(ineq),[pp/lc,0]),xx)))) %or (prt("case2-3"),(fineq(lc<0) %and (sineq(apply(op(ineq),[(-1)*pp/lc,0]),xx))))) else ((prt("case2-1"),((lc=0) %and fineq(apply(op(ineq),[expand(pp-lc*xx^dg),0]),xx))) %or (prt("case2-2"),((lc>0) %and (sineq(apply(op(ineq),[pp/lc,0]),xx)))) %or (prt("case2-3"),((lc<0) %and (sineq(apply(op(ineq),[(-1)*pp/lc,0]),xx))))) )$ /* automatic quantify */ fullex(f):=qe(map(lambda([x],[E,x]),listofvars(f)),p2t(f))$ fullall(f):=qe(map(lambda([x],[A,x]),listofvars(f)),p2t(f))$ /* lineq */ defrule(tolineq1, (aa)^(bb), lineqdec(1,aa^bb,bb))$ defrule(tolineq2, (aa)/(bb), lineqdec(2,aa/bb,bb))$ lineq(f):=block([L,lineqdec,q], if is(lsimp=off) then (prt("lineq off"),return(f)), prt("linear simplifing..."), L:[], lineqdec(m,g,bb):=block([strg], if (m=1 and denom(bb)=1) or (m=2 and listofvars(bb)=[]) then g else ( strg:ssubst("equal","=",base64(string(g))), strg:ssubst("plus","+",strg), strg:ssubst("slash","/",strg), L:append(L,[eval_string(strg)=g]),eval_string(strg))), q:qe([],apply1(f,tolineq1,tolineq2)),L:listify(setify(L)), prt([q,L]), /* nnscan(nnscan(subst(L,q)))*/ /* nnscan(subst(L,q)) */ nnscan(mapply(lambda([t],(t)%and(fullex(t))),dnf(subst(L,q)))) /* subst(L,q) */ )$ /* l2r */ defrule(l2r1, (aa) > (bb),block([p,l,v],p:expand(bb-aa),l:listofvars(p), if length(l)=1 then (v:first(l), if hipow(p,v)=1 and coeff(p,v,1)>0 then v<-coeff(p,v,0)/coeff(p,v,1) elseif hipow(p,v)=1 and coeff(p,v,1)<0 then -coeff(p,v,0)/coeff(p,v,1)<v else bb<aa) else bb<aa ))$ defrule(l2r2, (aa) >= (bb),block([p,l,v],p:expand(bb-aa),l:listofvars(p), if length(l)=1 then (v:first(l), if hipow(p,v)=1 and coeff(p,v,1)>0 then v<=-coeff(p,v,0)/coeff(p,v,1) elseif hipow(p,v)=1 and coeff(p,v,1)<0 then -coeff(p,v,0)/coeff(p,v,1)<=v else bb<=aa) else bb<=aa ))$ defrule(l2r3, (bb) < (aa),block([p,l,v],p:expand(bb-aa),l:listofvars(p), if length(l)=1 then (v:first(l), if hipow(p,v)=1 and coeff(p,v,1)>0 then v<-coeff(p,v,0)/coeff(p,v,1) elseif hipow(p,v)=1 and coeff(p,v,1)<0 then -coeff(p,v,0)/coeff(p,v,1)<v else bb<aa) else bb<aa ))$ defrule(l2r4, (bb) <= (aa),block([p,l,v],p:expand(bb-aa),l:listofvars(p), if length(l)=1 then (v:first(l), if hipow(p,v)=1 and coeff(p,v,1)>0 then v<=-coeff(p,v,0)/coeff(p,v,1) elseif hipow(p,v)=1 and coeff(p,v,1)<0 then -coeff(p,v,0)/coeff(p,v,1)<=v else bb<=aa) else bb<=aa ))$ defrule(l2r5, (bb) = (aa),block([p,l,v],p:expand(bb-aa),l:listofvars(p), if length(l)=1 then (v:first(l), if hipow(p,v)=1 then v=-coeff(p,v,0)/coeff(p,v,1) else bb=aa) else bb=aa ))$ defrule(l2r1x, (aa) > (bb),block([p,v],p:expand(bb-aa), if length(vv)=1 then (v:first(vv), if hipow(p,v)=1 and coeff(p,v,1)>0 then v<-coeff(p,v,0)/coeff(p,v,1) elseif hipow(p,v)=1 and coeff(p,v,1)<0 then -coeff(p,v,0)/coeff(p,v,1)<v else bb<aa) else bb<aa ))$ defrule(l2r2x, (aa) >= (bb),block([p,v],p:expand(bb-aa), if length(vv)=1 then (v:first(vv), if hipow(p,v)=1 and coeff(p,v,1)>0 then v<=-coeff(p,v,0)/coeff(p,v,1) elseif hipow(p,v)=1 and coeff(p,v,1)<0 then -coeff(p,v,0)/coeff(p,v,1)<=v else bb<=aa) else bb<=aa ))$ defrule(l2r3x, (bb) < (aa),block([p,v],p:expand(bb-aa), if length(vv)=1 then (v:first(vv), if hipow(p,v)=1 and coeff(p,v,1)>0 then v<-coeff(p,v,0)/coeff(p,v,1) elseif hipow(p,v)=1 and coeff(p,v,1)<0 then -coeff(p,v,0)/coeff(p,v,1)<v else bb<aa) else bb<aa ))$ defrule(l2r4x, (bb) <= (aa),block([p,v],p:expand(bb-aa), if length(vv)=1 then (v:first(vv), if hipow(p,v)=1 and coeff(p,v,1)>0 then v<=-coeff(p,v,0)/coeff(p,v,1) elseif hipow(p,v)=1 and coeff(p,v,1)<0 then -coeff(p,v,0)/coeff(p,v,1)<=v else bb<=aa) else bb<=aa ))$ defrule(l2r5x, (bb) = (aa),block([p,v],p:expand(bb-aa), if length(vv)=1 then (v:first(vv), if hipow(p,v)=1 then v=-coeff(p,v,0)/coeff(p,v,1) else bb=aa) else bb=aa ))$ defrule(l2r6, (bb) < (aa) %and (cc) < (bb), (cc) < (bb) %and (bb) < (aa))$ l2r(f):=block([],apply1(f,l2r1,l2r2,l2r3,l2r4,l2r5))$ l2rx(f,[vv]):=block([],apply1(f,l2r1x,l2r2x,l2r3x,l2r4x,l2r5x,l2r6))$ l2rxsort(f,g):=member(x,listofvars(f)) or not(member(x,listofvars(g)))$ /* bra, ora */ bra(f):=block([g], if listp(f) then return(f), if atom(f) then return([[f]]), g:dnf(f), if op(g)="%or" then g:args(g) else g:[g], map(lambda([t],if op(t)="%and" then args(t) else [t]),g) )$ brac(f):=block([g], if listp(f) then return(f), if atom(f) then return([[f]]), g:cnf(f), if op(g)="%and" then g:args(g) else g:[g], map(lambda([t],if op(t)="%or" then args(t) else [t]),g) )$ ora(f):=block([], if not(listp(f)) then return(f), if atom(f) then return(f) else substpart("%or",map(lambda([g],substpart("%and",g,0)),f),0) )$ orac(f):=block([], if not(listp(f)) then return(f), if atom(f) then return(f) else substpart("%and",map(lambda([g],substpart("%or",g,0)),f),0) )$ /* lfactor */ lfactor0(ff):=block([f,L,M,cf],L:[],M:[], if not(listp(ff)) then f:bra(ff) else f:ff, cf:delete({},map(lambda([t],apply(intersection,args(t))),sort(rest(listify( setdifference(powerset(setify(map(setify,f))),map(lambda([s],{s}),fullsetify(f))) )),comparelength2))), if cf=[] then return(f), cf:first(cf), for fk in f do if subsetp(cf,setify(fk)) then L:append([substpart("%and",setdifference(setify(fk),cf),0)],L) else M:append([fk],M), append(M,[flatten([substpart("%or",L,0),args(cf)])]) )$ lfactor(ff):=block([f,afr,bfr], if not(listp(ff)) then f:bra(ff) else f:ff, afr:f,bfr:[], for cc:1 unless length(bfr)=length(afr) do (bfr:afr,afr:lfactor0(bfr)), factor(afr))$ /* plineq */ plineq(f,[p]):=sort(append(Lminus(f,part(f,p)),bra(lineq(ora(part(f,p))))))$ /* addex */ addex(L,[v]):=l2r(map( lambda([Lk],append(Lk,first(bra(qe(makelist([E,vk],vk,v),p2t(substpart("%and",Lk,0))))))), L))$ /* qex, solvex */ qex(L,[f]):=if listp(L) then qe(L,apply(p2t,f)) else p2t(L)$ solvex(f,[v]):=apply(nnsolvexx,flatten([ora(bra(f)),v]))$ /* tineq */ sorted_cineq(sols,lc,gl):=block([L,M], L:sols, if (lc>0 and (gl="<" or gl="<=")) or (lc<0 and (gl=">" or gl=">=")) then L:append([first(L)],L), L:makelist( if k=1 then rhs(part(L,k))<lhs(part(L,k)) elseif evenp(k) and k<length(L) then (rhs(part(L,k+1))<lhs(part(L,k+1))) %and (lhs(part(L,k))<rhs(part(L,k))) elseif evenp(k) and k=length(L) then lhs(part(L,k))<rhs(part(L,k)) else false, k,length(L)), if (lc>0 and (gl="<" or gl="<=")) or (lc<0 and (gl=">" or gl=">=")) then L:rest(L), if gl=">=" or gl="<=" then L:subst("<"="<=",L), L:reverse(L),prt(L), apply1(substpart("%or",L,0),delsame,eqsame) )$ sorted_cineq5(sols,lc,gl):=block([L,M], L:sols, if (lc>0 and (gl="<" or gl="<=")) or (lc<0 and (gl=">" or gl=">=")) then L:append([first(L)],L), L:makelist( if k=1 then rhs(part(L,k))<lhs(part(L,k)) elseif evenp(k) and k<length(L) then (rhs(part(L,k+1))<lhs(part(L,k+1))) %and (lhs(part(L,k))<rhs(part(L,k))) elseif evenp(k) and k=length(L) then lhs(part(L,k))<rhs(part(L,k)) else false, k,length(L)), if (lc>0 and (gl="<" or gl="<=")) or (lc<0 and (gl=">" or gl=">=")) then L:rest(L), if gl=">=" or gl="<=" then L:subst("<"="<=",L), L:reverse(L),prt(L), apply1(substpart("%or",L,0),delsame,eqsame) )$ tineq2(sols,lc,gl):=block([x,PLx,delx,PL,IL,CL,M,N], /*tineq2(sols,lc,gl):=block([],*/ if sols=[] then return([[["%and"(apply(gl,[lc,0]))],[ ]]]), x:lhs(first(sols)), PLx:map(flatten,listify(permutations(sols))),prt(PLx), delx(l):=map(lambda([xa],rhs(xa)),l),PL:map(delx,PLx),prt(PL), IL:flatten(map(lambda([s],substpart("%and",makelist(part(s,k)-part(s,k+1)>=0,k,length(s)-1),0)),PL)),prt(IL), ncondtype_out:ncondtype,postsimp_out:postsimp,ncondtype:off,postsimp:off, CL:map(lambda([s],ora(l2r(solvex(s)))),IL), ncondtype:ncondtype_out,postsimp:postsimp_out,prt("case condition:"),prt(CL), CL:map(lambda([s],if atom(s) then s elseif op(s)="=" then false else s),CL), if not(casecheck=off) then ( M:[],CL:makelist(if member(e,M) then false elseif is(substpart("%or",makelist(qex((e)%implies(m)),m,M),0)=true) then false else (M:append(M,[e]),e),e,CL), M:[],CL:reverse(makelist(if member(e,M) then false elseif is(substpart("%or",makelist(qex((e)%implies(m)),m,M),0)=true) then false else (M:append(M,[e]),e),e,reverse(CL))) ),prt(CL), N:delete(false,makelist(if part(CL,j)#false then [[sorted_cineq(part(PLx,j),lc,gl)],[part(CL,j)]],j,length(CL))), /* N:sepaxa(map(flatten,bra(ora(N))),x), */ thistineq2:[PLx,PL,IL,CL,N], prt("tineq2 returned..."),prt(N),N /* M:map(lambda([l],sorted_cineq(l,lc,gl)),PLx), prt(M),map("[",map("%and",M,CL))*/ )$