# This is a program for real root isolation of bivariate polynomial systems with two polynomials. # We provide two commands: # 1: [R,SR]=birootfinding_cotrans(sys,var,B0,err,k) # 2: [R,SR]=birealrootfindingpoly(sys,var,B0,err,k) # The first one is based on coordinate transformation and then we need only to find real roots in subboxes of [-1,1]*[-1,1] # The second one is that finding real roots in the given box B0 ####################### # Parameters: # Input: # sys: A bivariate polynomial system (f1,f2) # var: list of variables # B0: An initial box # err: A termination precision and err>0 # k: An integer for spliting the x-axis # Output: # R: The isolating box set of the given systen # SR: The suspected root box set of the given systen ###################### # f1:=x^2+y^2-1;f2:=y-x^2; #A:=birootfinding_cotrans([f1,f2],[x,y],[[-10,10],[-10,10]],10^(-4),11); # [[[[-0.818181818181818, -0.636363636363636], [0.574959516525269, 0.669421553611755]], [ [0.63636363636364, 0.81818181818182], # [0.574959516525269, 0.669421553611755]]], []] # R:=A[1];SR:=A[2]; with(RootFinding): with(ListTools): with(VectorCalculus): with(LinearAlgebra): Digits:=15: birootfinding_cotrans:=proc(sys,var,B0,err) local uB,Q,T,i,j,B,R,SR,t,L,temp,F,k; if nargs<5 then k:=11; else k:=args[5]; fi; if nops(sys)<>2 or nops(var)<>2 then ERROR("The input system must be a one with two functions and two variables!"); fi; R:=[];SR:=[]; L:=allbinary(2); uB:=unitbox(B0); for i from 1 to 4 do temp:=transformsys(sys,var,uB,L[i]); F:=temp[1];B:=temp[2]; for j from 1 to nops(B) do T:=birealrootfindingpoly(F,var,B[j],err,k); Q:=transformbox(T[1],L[i]); R:=[op(R),op(Q)]; Q:=transformbox(T[2],L[i]); SR:=[op(SR),op(Q)]; od; od; T:=SR;SR:=[]; for i from 1 to nops(T) do Q:=birealrootfindingpoly(sys,var,T[i],err,1); R:=[op(R),op(Q[1])]; SR:=[op(SR),op(Q[2])]; od; return [R,SR]; end: transformbox:=proc(B,L) local n,Q,T,i,j,t; n:=nops(L);Q:=[]; for i from 1 to nops(B) do; T:=B[i]; for j from 1 to nops(L) do if L[j]=1 then t:=[1/T[j][2],1/T[j][1]]; T[j]:=t; fi; od; Q:=[op(Q),T]; od; return Q; end: unitbox:=proc(B) local n,i,V,Q,T; n:=nops(B);V:=[]; for i from 1 to n do T:=B[i]; if T[1]<-1 and T[2]>1 then Q:=[[-1,1/T[1]],[-1,1],[1/T[2],1]]; elif T[1]<-1 then if T[2]<=-1 then Q:=[[1/T[2],1/T[1]],[],[]]; else Q:=[[-1,1/T[1]],[-1,T[2]],[]]; fi; elif T[2]>1 then if T[1]>=1 then Q:=[[],[],[1/T[2],1/T[1]]]; else Q:=[[],[T[1],1],[1/T[2],1]]; fi; else Q:=[[],T,[]]; fi; V:=[op(V),Q]; od; return V; end: transformsys:=proc(sys,var,B,L) local n,F,T,Q,t,i,j,k,a; n:=nops(sys);F:=sys;Q:=[];k:=1; if L[1]=0 then if B[1][2]=[] then k:=0; else Q:=[[B[1][2]]]; fi; fi; if L[1]=1 then F:=subs(var[1]=1/var[1],F); if B[1][1]=[] and B[1][3]=[] then k:=0 else if B[1][1]<>[] then Q:=[[B[1][1]]]; fi; if B[1][3]<>[] then Q:=[op(Q),[B[1][3]]]; fi; fi; fi; i:=2; while i<=n and k=1 do if L[i]=0 then if B[i][2]=[] then k:=0; else T:=[]; for j from 1 to nops(Q) do a:=[op(Q[j]),B[i][2]]; T:=[op(T),a]; od; Q:=T; fi; fi; if L[i]=1 then F:=subs(var[i]=1/var[i],F); if B[i][1]=[] and B[i][3]=[] then k:=0 else T:=[]; if B[i][1]<>[] then for j from 1 to nops(Q) do a:=[op(Q[j]),B[i][1]]; T:=[op(T),a]; od; fi; if B[i][3]<>[] then for j from 1 to nops(Q) do a:=[op(Q[j]),B[i][3]]; T:=[op(T),a]; od; fi; Q:=T; fi; fi; i:=i+1 od; F:=simplify(F); F:=numer(F); if k=0 then Q:=[]; fi; Q:=evalf(Q); return [F,Q]; end: allbinary:=proc(n) local i,j,T,Q; T:=[[0],[1]]; for i from 2 to n do Q:=[]; for j from 1 to nops(T) do Q:=[op(Q),[op(T[j]),0]]; Q:=[op(Q),[op(T[j]),1]]; od; T:=Q; od; return T; end: birealrootfindingpoly:=proc(sys,var,Box0,err,k) local w,R,SR,Q,wx,i,ax,bx,F1,F2,F,a1,mx,ind,Q0,B,B0,c_1,c_2; ################# B0:=evalf(Box0); R:=[];SR:=[];Q:=[]; wx:=(B0[1][2]-B0[1][1])/k; for i from 1 to k do ax:=B0[1][1]+(i-1)*wx; bx:=B0[1][1]+i*wx; F1:=Boundpoly(sys[1],var[1],var[2],[ax,bx]); F2:=Boundpoly(sys[2],var[1],var[2],[ax,bx]); c_1:=candidate_interval(F1,var[2],B0[2]); c_2:=candidate_interval(F2,var[2],B0[2]); Q0:=candidates(c_1[3],c_2[3],[ax,bx]); Q:=[op(Q),op(evalf(Q0))]; od; while nops(Q)>0 do B:=Q[1];ind:=0; Q:=subsop(1=NULL,Q); F:=preconditioner(sys,var,B); a1:=IsOMsys(F,var,B); if a1>0 then ind:=Existence(F,var,B); if ind[1]=1 then R:=[op(R),B]; elif ind[1]=1/2 then SR:=[op(SR),B]; fi; else w:=B[1][2]-B[1][1]; if w>err then mx:=(B[1][1]+B[1][2])/2; F1:=Boundpoly(sys[1],var[1],var[2],[B[1][1],mx]); F2:=Boundpoly(sys[2],var[1],var[2],[B[1][1],mx]); c_1:=candidate_interval(F1,var[2],B[2]); c_2:=candidate_interval(F2,var[2],B[2]); Q0:=candidates(c_1[3],c_2[3],[B[1][1],mx]); Q:=[op(Q),op(evalf(Q0))]; F1:=Boundpoly(sys[1],var[1],var[2],[mx,B[1][2]]); F2:=Boundpoly(sys[2],var[1],var[2],[mx,B[1][2]]); c_1:=candidate_interval(F1,var[2],B[2]); c_2:=candidate_interval(F2,var[2],B[2]); Q0:=candidates(c_1[3],c_2[3],[mx,B[1][2]]); Q:=[op(Q),op(evalf(Q0))]; else SR:=[op(SR),B]; fi; fi; od; if nops(SR)>1 then SR:=collectsuspectedbox(SR,err); fi; return [R,SR]; end: exclusion:=proc(sys,var,B) local fB,ind; ind:=0; fB:=Intervalfun(sys[1],var,B); if fB[1]>0 or fB[2]<0 then ind:=1; else fB:=Intervalfun(sys[2],var,B); if fB[1]>0 or fB[2]<0 then ind:=1; fi; fi; return ind; end: splitB:=proc(B,z) local n,T,m,Q,temp,i,j,a,k,w; n:=nops(B);T:=[B]; for i from 1 to n do m:=nops(T);w:=(B[i][2]-B[i][1])/z;Q:=[]; for j from 1 to m do for k from 0 to z-1 do; a:=[B[i][1]+k*w,B[i][1]+(k+1)*w]; temp:=subsop(i=a,T[j]); Q:=[op(Q),temp]; od; od; T:=Q; od; return T; end: shiftbox:=proc(F,var,t,B) local w1,w2,w,B0; w1:=width(B[1]);w2:=width(B[2]); w:=min(w1,w2)/5; if t=1 or t=12 then B0:=[[B[1][1]-w,B[1][2]],[B[2][1],B[2][2]+w]]; fi; if t=2 or t=23 then B0:=[[B[1][1],B[1][2]+w],[B[2][1],B[2][2]+w]]; fi; if t=3 or t=34 then B0:=[[B[1][1],B[1][2]+w],[B[2][1]-w,B[2][2]]]; fi; if t=4 or t=41 then B0:=[[B[1][1]-w,B[1][2]],[B[2][1]-w,B[2][2]]]; fi; return B0; end: Intervalfun:=proc(f,var,B) local n,A,i,a,g,fp,b,Lc,Lt,newg; n:=nops(B);A:=[]; if type(f,numeric)=true then a:=f;b:=f; else for i from 1 to n do a:=INTERVAL(B[i][1]..B[i][2]); A:=[op(A),var[i]=a]; od; g:=convert(f,'horner'); fp:=evalr(subs(op(A),g)); if nops(fp)=1 then fp:=op(fp); if nops(fp)=2 then a:=op(1,fp);b:=op(2,fp); else g:=expand(f); Lc:=lcoeff(g, order = plex(op(var)), 'm'); Lt:=Lc*m; newg:=g-Lc*m+newvar*Lc*m; fp:=op(evalr(subs(op(A),newvar=INTERVAL(1..1),newg))); if nops(fp)=2 then a:=op(1,fp);b:=op(2,fp); else a:=min(fp);b:=max(fp); fi; fi; else a:=min(fp);b:=max(fp); fi; fi; return [a,b]; end: Intervalmulti:=proc(A,B) local a1,a2,a3,a4,a,b; a1:=A[1]*B[1]; a2:=A[2]*B[1]; a3:=A[1]*B[2]; a4:=A[2]*B[2]; a:=min([a1,a2,a3,a4]); b:=max([a1,a2,a3,a4]); return [a,b]; end: Intervalfun_mean:=proc(f,var,B) local n,P,p,A,Dif,DifB,T,i,temp,fp; n:=nops(var);P:=[];A:=[0,0]; for i from 1 to n do p:=(B[i][1]+B[i][2])/2; P:=[op(P),var[i]=p]; Dif:=diff(f,var[i]); DifB:=Intervalfun(Dif,var,B); temp:=B[i]-[p,p]; T:=Intervalmulti(DifB,temp); A:=A+T; od; fp:=subs(op(P),f); A:=[fp,fp]+A; return A; end: preconditioner:=proc(sys0,var,B) local p,Jp,A,sys,E,a; p:=[(B[1][1]+B[1][2])/2,(B[2][1]+B[2][2])/2]; Jp:=Jacobian(sys0 ,var = p); if Rank(Jp)=2 then a:=Jp[1][1]*Jp[2][2]-Jp[1][2]*Jp[2][1]; E:=Matrix(2, 2, [[Jp[2][2]/a,-Jp[1][2]/a],[-Jp[2][1]/a,Jp[1][1]/a]]); A:=Matrix(2, 2, [[E[1][1]+E[2][1],E[1][2]+E[2][2]],[E[1][1]-E[2][1],E[1][2]-E[2][2]]]); else A:=Matrix([[1,0],[0,1]]); fi; sys:=A.; sys:=convert(sys,list); return sys; end: IsOMsys:=proc(sys,var,B) local fx,fy,gx,gy,i,ind,i1,i2,i3,i4; ind:=0; fx:=diff(sys[1],var[1]); fy:=diff(sys[1],var[2]); i1:=Intervalfun(fx,var,B); i2:=Intervalfun(fy,var,B); gx:=diff(sys[2],var[1]); gy:=diff(sys[2],var[2]); i3:=Intervalfun(gx,var,B); i4:=Intervalfun(gy,var,B); if i1[1]>0 then if i2[1]>0 then if i3[1]>0 then if i4[2]<0 then ind:=1; fi; fi; fi; fi; return ind; end: Findindex:=proc(A,a) local n,i,t; n:=nops(A);i:=1;t:=0; while i<=n and t=0 do if A[i]=a then t:=i; fi; i:=i+1; od; return t; end: whichcase:=proc(f,g) local sf,sg,t,ind; sf:=sort(f); sg:=sort(g); if sf[1]>0 or sf[4]<0 or sg[1]>0 or sg[4]<0 then ind:=0; elif (sf[1]=0 and sf[2]>0) or (sf[4]=0 and sf[3]<0) then t:=Findindex(f,0); if g[t]=0 then ind:=1; else ind:=0; fi; elif (sg[1]=0 and sg[2]>0) or (sg[4]=0 and sg[3]<0) then t:=Findindex(g,0); if f[t]=0 then ind:=1; else ind:=0; fi; else ind:=2; fi; return [ind,t]; end: bisection:=proc(f,g,var,L0,F0,G0) local F,G,s,L,w,fp,gp,p; F:=F0;G:=G0;s:=0;L:=L0; if G[1]*G[2]>0 then s:=G[1]; else w:=abs(L[2]-L[1]); while w>10^(-10) and s=0 do p:=(L[1]+L[2])/2; fp:=evalf(subs(var=p,f)); gp:=evalf(subs(var=p,g)); w:=abs(L[2]-L[1]); if F[1]*fp<=0 and G[1]*gp<=0 then L[2]:=p; elif F[2]*fp<=0 and G[2]*gp<=0 then L[1]:=p; else s:=gp; fi; od; fi; return s; end: Existence:=proc(sys,var,B) local f,g,uf,ug,p1,p2,p3,p4,f1,f2,f3,f4,g1,g2,g3,g4,s,c,t,ind,P,G,F,L,S; f:=sys[1];g:=sys[2]; p1:=[B[1][1],B[2][2]]; p2:=[B[1][2],B[2][2]]; p3:=[B[1][2],B[2][1]]; p4:=[B[1][1],B[2][1]]; P:=[p1,p2,p3,p4]; f1:=evalf(subs(var[1]=p1[1],var[2]=p1[2],f)); f2:=evalf(subs(var[1]=p2[1],var[2]=p2[2],f)); f3:=evalf(subs(var[1]=p3[1],var[2]=p3[2],f)); f4:=evalf(subs(var[1]=p4[1],var[2]=p4[2],f)); g1:=evalf(subs(var[1]=p1[1],var[2]=p1[2],g)); g2:=evalf(subs(var[1]=p2[1],var[2]=p2[2],g)); g3:=evalf(subs(var[1]=p3[1],var[2]=p3[2],g)); g4:=evalf(subs(var[1]=p4[1],var[2]=p4[2],g)); F:=[f1,f2,f3,f4];G:=[g1,g2,g3,g4]; c:=whichcase([f1,f2,f3,f4],[g1,g2,g3,g4]); ################################################ if c[1]=0 then ind:=0; elif c[1]=1 then ind:=1/2; t:=c[2]; else S:=[]; if f1=0 then f1:=1; fi; if f2=0 then f2:=1; fi; if f3=0 then f3:=1; fi; if f4=0 then f4:=1; fi; if f1*f2<0 then L:=B[1]; uf:=subs(var[2]=B[2][2],f); ug:=subs(var[2]=B[2][2],g); s:=bisection(uf,ug,var[1],L,F[1..2],G[1..2]); S:=[op(S),s]; if s=0 then t:=12; fi; fi; if f2*f3<0 then L:=[B[2][2],B[2][1]]; uf:=subs(var[1]=B[1][2],f); ug:=subs(var[1]=B[1][2],g); s:=bisection(uf,ug,var[2],L,F[2..3],G[2..3]); S:=[op(S),s]; if s=0 then t:=23; fi; fi; if f3*f4<0 then L:=[B[1][2],B[1][1]]; uf:=subs(var[2]=B[2][1],f); ug:=subs(var[2]=B[2][1],g); s:=bisection(uf,ug,var[1],L,F[3..4],G[3..4]); S:=[op(S),s]; if s=0 then t:=34; fi; fi; if f4*f1<0 then L:=B[2]; uf:=subs(var[1]=B[1][1],f); ug:=subs(var[1]=B[1][1],g); s:=bisection(uf,ug,var[2],L,[F[4],F[1]],[G[4],G[1]]); S:=[op(S),s]; if s=0 then t:=41; fi; fi; if nops(S)>1 then if S[1]*S[2]<0 then ind:=1; elif S[1]*S[2]>0 then ind:=0; else ind:=1/2; fi; else ind:=1/2; t:=0; fi; fi; return [ind,t]; end: collectsuspectedbox:=proc(B,err) local a,A,b,c,MR,i; A:=sort(B); a:=A[1];MR:=[]; for i from 2 to nops(A) do b:=A[i]; if abs(a[1][1]-b[1][1])=b[2][1] and a[2][1]<=b[2][2] or a[2][2]>=b[2][1] and a[2][2]<=b[2][2] then c:=[[min(op(a[1]),op(b[1])),max(op(a[1]),op(b[1]))],[min(op(a[2]),op(b[2])),max(op(a[2]),op(b[2]))]]; a:=c; elif b[2][1]>=a[2][1] and b[2][1]<=a[2][2] or b[2][2]>=a[2][1] and b[2][2]<=a[2][2] then c:=[[min(op(a[1]),op(b[1])),max(op(a[1]),op(b[1]))],[min(op(a[2]),op(b[2])),max(op(a[2]),op(b[2]))]]; a:=c; else MR:=[op(MR),a]; a:=b; fi; else MR:=[op(MR),a]; a:=b; fi; od; MR:=[op(MR),a]; return MR; end: candidates:=proc(c1,c2,Bx) local Q,T,temp,i; T:=[]; Q:=intersect_interval_set(c1,c2); for i from 1 to nops(Q) do temp:=[Bx,Q[i]]; T:=[op(T),temp]; od; return T; end: Boundpoly:=proc(f,x,y,Bx) local i,n,g,pfu,pfd,nfu,nfd,T,a,b; n:=degree(f,y); pfu:=0; pfd:=0; nfu:=0; nfd:=0; for i from 0 to n do g:=coeff(f,y,i); if type(g,numeric)=true then a:=g;b:=g; else T:=Intervalfun(g,[x],[Bx]); a:=T[1];b:=T[2]; fi; pfu:=pfu+b*y^i; pfd:=pfd+a*y^i; if (i mod 2 = 1) then nfu:=nfu+a*y^i; nfd:=nfd+b*y^i; else nfu:=nfu+b*y^i; nfd:=nfd+a*y^i; fi; od; return [pfu,pfd,nfu,nfd]; end: sort_real_root:=proc(A) local n,temp,T,H,minT,maxT,i,s; n:=nops(A);H:=[]; if n=1 then T:=op(A[1])[2]; if T[1]=T[2] then H:=[[T[1]-10^(-8),T[2]+10^(-8)]]; else H:=[T]; fi; else for i from 1 to n do T:=op(A[i])[2]; if T[1]=T[2] then if i=1 then minT:=T[1]-10^(-8); temp:=op(A[i+1])[2]; s:=temp[1]-T[2]; if s<10^(-8) then maxT:=T[2]+s/2; else maxT:=T[2]+10^(-8); fi; elif i=n then maxT:=T[2]+10^(-8); temp:=op(A[i-1])[2]; s:=T[1]-temp[2]; if s<10^(-8) then minT:=T[1]-s/2; else minT:=T[1]-10^(-8); fi; else temp:=op(A[i+1])[2]; s:=temp[1]-T[2]; if s<10^(-8) then maxT:=T[2]+s/2; else maxT:=T[2]+10^(-8); fi; temp:=op(A[i-1])[2]; s:=T[1]-temp[2]; if s<10^(-8) then minT:=T[1]-s/2; else minT:=T[1]-10^(-8); fi; fi; T:=[minT,maxT]; fi; H:=[op(H),T]; od; fi; return H; end: MyIsolate:=proc(f,x) local big_0,less_0,contain_0,a,p,T,A,v,n,i; Digits:=5; if f=0 then big_0:=[];less_0:=[];contain_0:=[-infinity,infinity]; else A:=Isolate(f,output='interval'); v:=infinity;n:=nops(A); big_0:=[];less_0:=[];contain_0:=[]; Digits:=15; if n=0 then a:=subs(x=0,f); if a>0 then big_0:=[[-v,v]]; else less_0:=[[-v,v]]; fi; else A:=sort_real_root(A);contain_0:=A; for i from 1 to n do if i=1 then a:=subs(x=A[1][1]-abs(A[1][1])-1,f); if a>0 then big_0:=[op(big_0),[-v,A[1][1]]]; else less_0:=[op(less_0),[-v,A[1][1]]]; fi; else p:=(A[i-1][2]+A[i][1])/2; a:=subs(x=p,f); if a>0 then big_0:=[op(big_0),[A[i-1][2],A[i][1]]]; else less_0:=[op(less_0),[A[i-1][2],A[i][1]]]; fi; fi; od; a:=subs(x=A[n][2]+abs(A[n][2])+1,f); if a>0 then big_0:=[op(big_0),[A[n][2],v]]; else less_0:=[op(less_0),[A[n][2],v]]; fi; fi; fi; return [big_0,less_0,contain_0]; end: mymin_set:=proc(B) local a,ans; if B<>[] then a:=sort(B); if a[-1][1]=-infinity then ans:=-infinity; else ans:=a[1][1]; fi; else ans:=infinity; fi; return ans; end: mymax_set:=proc(B) local a,ans; if B<>[] then a:=sort(B,key = (x->x[2])); ans:=a[-1][2]; else ans:=-infinity; fi; return ans; end: intersect_interval:=proc(A,B) local T,a,b; if A[2]B[2] then T:=[]; else a:=sort([A[1],B[1]]);b:=sort([A[2],B[2]]); T:=[a[2],b[1]]; fi; return T; end: mysort:=proc(A) local T,a; if A<>[] then T:=sort(A,key = (x->x[1])); if T[-1][1]=-infinity then a:=T[-1]; T:=subsop(-1=NULL,T); T:=[a,op(T)]; fi; else T:=[]; fi; return T; end: intersect_interval_set:=proc(A1,B1) local A,B,T,a,j,n,ind,temp,m,i,t,s,b; a:=mymin_set(A1);b:=mymin_set(B1); if a>=b then A:=mysort(A1);B:=mysort(B1); else A:=mysort(B1);B:=mysort(A1); fi; T:=[]; i:=1;n:=nops(B);m:=nops(A);j:=1; while i<=m and j<=n do a:=A[i];ind:=0;t:=1; while j<=n and t=1 do temp:=intersect_interval(a,B[j]); if nops(temp)>0 then T:=[op(T),temp]; ind:=1; s:=j; fi; j:=j+1; if ind=1 and nops(temp)=0 then t:=0; fi; od; i:=i+1; if ind=0 then j:=1; else j:=s; fi; od; return T; end: complementary_set:=proc(A,B) local n,i,T; n:=nops(B); if nops(B)=0 then T:=[A]; else if A[1]1 then for i from 1 to n-1 do T:=[op(T),[B[i][2],B[i+1][1]]]; od; fi; if B[n][2]=0 then fu:=F[1];fd:=F[2]; c_u:=MyIsolate(fu,x); c_d:=MyIsolate(fd,x); big_0:=intersect_interval_set([T],c_d[1]); less_0:=intersect_interval_set([T],c_u[2]); temp:=mysort([op(big_0),op(less_0)]); contain_0:=complementary_set(T,temp); elif T[2]<=0 then fu:=F[3];fd:=F[4]; c_u:=MyIsolate(fu,x); c_d:=MyIsolate(fd,x); big_0:=intersect_interval_set([T],c_d[1]); less_0:=intersect_interval_set([T],c_u[2]); temp:=mysort([op(big_0),op(less_0)]); contain_0:=complementary_set(T,temp); else T1:=[0,T[2]]; fu:=F[1];fd:=F[2]; c_u:=MyIsolate(fu,x); c_d:=MyIsolate(fd,x); tb:=intersect_interval_set([T1],c_d[1]); tl:=intersect_interval_set([T1],c_u[2]); temp:=mysort([op(tb),op(tl)]); tc:=complementary_set(T1,temp); T1:=[T[1],0]; fu:=F[3];fd:=F[4]; c_u:=MyIsolate(fu,x); c_d:=MyIsolate(fd,x); xb:=intersect_interval_set([T1],c_d[1]); xl:=intersect_interval_set([T1],c_u[2]); temp:=mysort([op(xb),op(xl)]); xc:=complementary_set(T1,temp); if nops(xb)>0 and nops(tb)>0 then if xb[-1][2]=tb[1][1] then temp:=[xb[-1][1],tb[1][2]]; xb:=subsop(-1=NULL,xb); tb:=subsop(1=NULL,tb); big_0:=[op(xb),temp,op(tb)]; else big_0:=[op(xb),op(tb)]; fi; else big_0:=[op(xb),op(tb)]; fi; if nops(xl)>0 and nops(tl)>0 then if xl[-1][2]=tl[1][1] then temp:=[xl[-1][1],tl[1][2]]; xl:=subsop(-1=NULL,xl); tl:=subsop(1=NULL,tl); less_0:=[op(xl),temp,op(tl)]; else less_0:=[op(xl),op(tl)]; fi; else less_0:=[op(xl),op(tl)]; fi; if nops(xc)>0 and nops(tc)>0 then if xc[-1][2]=tc[1][1] then temp:=[xc[-1][1],tc[1][2]]; xc:=subsop(-1=NULL,xc); tc:=subsop(1=NULL,tc); contain_0:=[op(xc),temp,op(tc)]; else contain_0:=[op(xc),op(tc)]; fi; else contain_0:=[op(xc),op(tc)]; fi; fi; return [big_0,less_0,contain_0]; end: