print(`This is SCHAAL`); print(`The main program is dan(R,B,t,b,c,tt,tipe,slak)`); print(`Investigating {tx+by=z,tx+cy=z}, b > c > tt-1.`); print(`Input tt is the lower bound for b and c.`); print(`So, for example, since t=c=1 is done via human means,`); print(`we use tt=2 when t=1.`); print(`Use tipe=0 to keep both b and c as parameters.`); print(`Use tipe=1 to have b a parameter and c numeric.`); print(`Use tipe=2 to have b numeric and c a parameter.`); print(`Use tipe=3 to have both b and c numeric.`); print(``); print(`slak is the amount to add to the proposed upper bound.`); print(``); print(`The diagonal case, b=c, is handled`); print(`by the program diagdan(R,B,t,b).`); print(`Here, b is a parameter.`); print(``); print(`In each, R is the initial set of red elements,`); print(`B is the initial set of blue elements,`); print(``); # Rstart starting SET of red elements with(PolynomialTools): dan:=proc(Rstart,Bstart,t,b,c,tt,tipe,slak) local i,j,R,B,C,v,n1,n2,TB,TR,Rpost,Bpost: R:=Rstart: B:=Bstart: if nops(Rstart) < nops(Bstart) then if tipe=0 then B:=Bstart: R:=Rstart union Build(B,t,c,b,c,tt,slak)[1]: B:=Bstart union Build(R,t,b,b,c,tt,slak)[1]: elif tipe=1 then B:=Bstart: R:=Rstart union BuildCNumeric(B,t,c,b,c,tt,slak)[1]: B:=Bstart union BuildCNumeric(R,t,b,b,c,tt,slak)[1]: elif tipe=2 then B:=Bstart: R:=Rstart union BuildBNumeric(B,t,c,b,c,tt,slak)[1]: B:=Bstart union BuildBNumeric(R,t,b,b,c,tt,slak)[1]: else B:=Bstart: R:=Rstart union BuildBothNumeric(B,t,c,b,c,tt,slak)[1]: B:=Bstart union BuildBothNumeric(R,t,b,b,c,tt,slak)[1]: fi: else if tipe=0 then R:=Rstart: B:=Bstart union Build(R,t,b,b,c,tt,slak)[1]: R:=Rstart union Build(B,t,c,b,c,tt,slak)[1]: elif tipe=1 then R:=Rstart: B:=Bstart union BuildCNumeric(R,t,b,b,c,tt,slak)[1]: R:=Rstart union BuildCNumeric(B,t,c,b,c,tt,slak)[1]: elif tipe=2 then R:=Rstart: B:=Bstart union BuildBNumeric(R,t,b,b,c,tt,slak)[1]: R:=Rstart union BuildBNumeric(B,t,c,b,c,tt,slak)[1]: else R:=Rstart: B:=Bstart union BuildBothNumeric(R,t,b,b,c,tt,slak)[1]: R:=Rstart union BuildBothNumeric(B,t,c,b,c,tt,slak)[1]: fi: fi: print('R',R): print('B',B): n1:=-1: n2:=-1: while (n1+n2 <> 0) do if tipe=0 then TB:=Build(R,t,b,b,c,tt,slak): elif tipe=1 then TB:=BuildCNumeric(R,t,b,b,c,tt,slak): elif tipe=2 then TB:=BuildBNumeric(R,t,b,b,c,tt,slak): else TB:=BuildBothNumeric(R,t,b,b,c,tt,slak): fi: Bpost:=B union TB[1]: n1:=nops(Bpost minus B): if tipe=0 then TR:=Build(B,t,c,b,c,tt,slak): elif tipe=1 then TR:=BuildCNumeric(B,t,c,b,c,tt,slak): elif tipe=2 then TR:=BuildBNumeric(B,t,c,b,c,tt,slak): else TR:=BuildBothNumeric(B,t,c,b,c,tt,slak): fi: Rpost:=R union TR[1]: n2:=nops(Rpost minus R): R:=Rpost: B:=Bpost: C:=B intersect R: print('Ris',R): print('Bis',B): if nops(C)>0 then print(t*b*c+t^2*b+(t^2+1)*c+t^3+slak): RETURN('Intersection',C): fi: if n1=0 and n2=0 then RETURN('NoMore'): fi: od: 0: end: Build:=proc(X,t,v,b,c,tt,slak) local r,gg,eq,d,cnt,q,i,j,T,x,y,z,m,Q,w,f,g,h,s,cc: T:={}: # PROPOSED BOUND: tbc+(t^2)b+(t^2+1)c+t^3: for i from 1 to nops(X) do for j from 1 to nops(X) do x:=X[i]: y:=X[j]: T:=T union {simplify(t*x+v*y)}: # ASSUME b variable, c variable, and t integer if rem(y-t*x,v,v)=0 then T:=T union {simplify((y-t*x)/v)}: fi: eq:=expand(y-v*x): cc:=[coeffs(eq)]: # list of coefficients d:=indets(eq): # set of indeterminants if type(cc,list(integer)) then if (eq mod t = 0) then T:=T union {simplify(eq/t)}: elif type(eq mod t,quadratic(v)) and d={v} then q:=unapply(eq,v): if q(tt)>0 then cnt:=0: for z from 0 to t-1 do if type(q(z)/t,integer) then cnt:=cnt+1:fi: od: if cnt=t then T:=T union {simplify(eq/t)}:fi: fi: fi: fi: od: if rem(x,v+t,v)=0 then T:=T union {simplify((x)/(v+t))}:fi: od: Q:={}: for i from 1 to nops(T) do w:=T[i]: if (type(w,integer)=false and type(w,linear)=false and type(w,quadratic)=false) or type(w,quadratic(b))=true then Q:=Q union {w}: else f:=unapply(w,b,c): if type(w,integer) and (w<1 or w>t*(tt)^2+(2*t^2+1)*tt+t^3+slak) then Q:=Q union {w}:fi: # may assume b,c both at least 2 if type(w,linear) and (f(tt,tt)<1 or f(0,0)>t^3 or diff(w,b)<0 or (diff(w,c)<0 and (diff(w,c)+diff(w,b)<0)) or diff(w,c)+diff(w,b)>2*t^2+t*tt+1 or diff(w,b)>t^2+t*tt or diff(w,c)>t^2+t*tt+1) then Q:=Q union {w}: fi: if type(w,quadratic) then #either c^2 term or bc term r:=unapply(diff(w,b),c): s:=unapply(diff(w,c),c): if (f(0,0)>t^3 or f(tt,tt)<1 or coeff(w,c^2)<0 or coeff(s(0),b)<0 or (coeff(w,c^2)+coeff(s(0),b))>t or r(tt)<0 or (r(tt)=0 and f(0,tt)<1) ) then Q:=Q union {w}: else g:=f(b,0): gg:=f(0,c): if (coeff(g,b)>t^2 or coeff(gg,c)>t^2+1) then Q:=Q union {w}:fi: fi: fi: fi: od: T:=T minus Q: [T,nops(T)]: end: BuildBNumeric:=proc(X,t,v,b,c,tt,slak) local i,j,T,x,y,Q,w,f,bnd: T:={}: # PROPOSED BOUND: tbc+(t^2)b+(t^2+1)c+t^3: for i from 1 to nops(X) do for j from 1 to nops(X) do x:=X[i]: y:=X[j]: T:=T union {simplify(t*x+v*y)}: # ASSUME b integer, c parameter, and t integer if type(v,integer) then if (y-t*x mod v = 0) then T:=T union {simplify((y-t*x)/v)}: fi: else if rem(y-t*x,v,v)=0 then T:=T union {simplify((y-t*x)/v)}: fi: fi: if (y-v*x mod t = 0) then T:=T union {simplify((y-v*x)/t)}:fi: od: if type(v,integer) then if (x mod (v+t) = 0) then T:=T union {simplify((x)/(v+t))}:fi: else if rem(x,v+t,v)=0 then T:=T union {simplify((x)/(v+t))}:fi: fi: od: Q:={}: for i from 1 to nops(T) do w:=T[i]: if (type(w,integer)=false and type(w,linear)=false) then Q:=Q union {w}: else bnd:=t*b*tt+t^2*b+(t^2+1)*tt+t^3+slak: f:=unapply(w,c): if type(w,integer) and (w<1 or w>bnd) then Q:=Q union {w}:fi: if type(w,linear) and (f(tt)<1 or f(0)>t^2*b+t^3 or diff(w,c)<0 or diff(w,c)>t*b+t^2+1) then Q:=Q union {w}: fi: fi: od: T:=T minus Q: [T,nops(T)]: end: BuildCNumeric:=proc(X,t,v,b,c,tt,slak) local i,j,T,x,y,Q,w,f,bnd: T:={}: # PROPOSED BOUND: tbc+(t^2)b+(t^2+1)c+t^3: for i from 1 to nops(X) do for j from 1 to nops(X) do x:=X[i]: y:=X[j]: T:=T union {simplify(t*x+v*y)}: # ASSUME b variable, c integer, and t integer if type(v,integer) then if (y-t*x mod v = 0) then T:=T union {simplify((y-t*x)/v)}: fi: else if rem(y-t*x,v,v)=0 then T:=T union {simplify((y-t*x)/v)}: fi: fi: if (y-v*x mod t = 0) then T:=T union {simplify((y-v*x)/t)}:fi: od: if type(v,integer) then if (x mod (v+t) = 0) then T:=T union {simplify((x)/(v+t))}:fi: else if rem(x,v+t,v)=0 then T:=T union {simplify((x)/(v+t))}:fi: fi: od: Q:={}: for i from 1 to nops(T) do w:=T[i]: if (type(w,integer)=false and type(w,linear)=false) then Q:=Q union {w}: else bnd:=t*c^2+(2*t^2+1)*c+t^3+slak: f:=unapply(w,b): if type(w,integer) and (w<1 or w>bnd) then Q:=Q union {w}:fi: if type(w,linear) and (f(c)<1 or f(0)>(t^2+1)*c+t^3 or diff(w,b)<0 or diff(w,b)>t*c+t^2) then Q:=Q union {w}: fi: fi: od: T:=T minus Q: [T,nops(T)]: end: BuildBothNumeric:=proc(X,t,v,b,c,tt,slak) local i,j,T,x,y,Q,w,f,bnd: T:={}: bnd:=t*b*c+(t^2)*b+(t^2+1)*c+t^3+slak: for i from 1 to nops(X) do for j from 1 to nops(X) do x:=X[i]: y:=X[j]: T:=T union {simplify(t*x+v*y)}: # ASSUME b,c,t all integers if (y-t*x mod v = 0) then T:=T union {simplify((y-t*x)/v)}: fi: if (y-v*x mod t = 0) then T:=T union {simplify((y-v*x)/t)}: fi: od: if (x mod (v+t) = 0) then T:=T union {simplify((x)/(v+t))}:fi: od: Q:={}: for i from 1 to nops(T) do w:=T[i]: if type(w,integer)=false then Q:=Q union {w}: else if type(w,integer) and (w<1 or w>bnd) then Q:=Q union {w}:fi: fi: od: T:=T minus Q: [T,nops(T),bnd]: end: anom:=proc(Rstart,Bstart,t,tt) local f,X,i,j,R,B,C,v,n1,n2,TB,TR,Rpost,Bpost,Cnum,Rnum,Bnum: R:=Rstart: B:=Bstart: if nops(Rstart) < nops(Bstart) then B:=Bstart: R:=Rstart union BuildAnom(B,t,t,tt)[1]: #(X,t,v,tt) B:=Bstart union BuildAnom(R,t,2*t+1,tt)[1]: else R:=Rstart: B:=Bstart union BuildAnom(R,t,2*t+1,tt)[1]: R:=Rstart union BuildAnom(B,t,t,tt)[1]: fi: print('R',R): print('B',B): n1:=-1: n2:=-1: while (n1+n2 <> 0) do TB:=BuildAnom(R,t,2*t+1,tt): Bpost:=B union TB[1]: n1:=nops(Bpost minus B): TR:=BuildAnom(B,t,t,tt): Rpost:=R union TR[1]: n2:=nops(Rpost minus R): R:=Rpost: B:=Bpost: C:=B intersect R: print('Ris',R): print('Bis',B): if nops(C)>0 then RETURN('Intersection',C): fi: if (n1=0 and n2=0) then #X:=[]: #for j from 3 to c do #Rnum:={}: Bnum:={}: #for i from 1 to nops(R) do #f:=unapply(R[i],t): #Rnum:=Rnum union {f(j)}: #od: #for i from 1 to nops(B) do #f:=unapply(B[i],t): #Bnum:=Bnum union {f(j)}: #od: #Cnum:=Rnum intersect Bnum: #if nops(Cnum)>0 then #X:=[op(X),j]: #print('IntersectionWhenTis',j,Cnum):fi: #od: #RETURN('ValuesOfTwithIntersection',X): RETURN('NoMore'): fi: od: 0: end: BuildAnom:=proc(X,t,v,tt) local cc,i,j,T,x,y,Q,w,f,bnd,dif,dfn,L,m,g,N,n: # b=2t+1, c=t, tt is bound on t. Start with tt=3. # t is variable in this Build function T:={}: # PROPOSED BOUND: 6t^3+2t^2+4t: for i from 1 to nops(X) do for j from 1 to nops(X) do x:=X[i]: y:=X[j]: T:=T union {simplify(t*x+v*y)}: # ASSUME b,c linear in t, and t integer if (y-v*x mod t = 0) then T:=T union {simplify((y-v*x)/t)}:fi: if rem(y-t*x,v,t)=0 then T:=T union {simplify((y-t*x)/v)}: fi: od: if rem(x,v+t,t)=0 then T:=T union {simplify((x)/(v+t))}:fi: od: Q:={}: for i from 1 to nops(T) do w:=T[i]: if (type(w,integer)=false and type(w,linear)=false and type(w,quadratic)=false and type(w,cubic)=false) then Q:=Q union {w}: else bnd:=6*t^3+2*t^2+4*t: g:=unapply(bnd,t): f:=unapply(w,t): if type(w,integer) and (w<1 or w>g(tt)) then Q:=Q union {w}: else cc:=CoefficientList(w,t): dif:=simplify(bnd-w): dfn:=unapply(dif,t): N:=realroot(w,1/2): L:=realroot(dif,1/2): n:=max(seq(N[i][2],i=1..nops(N)) ): m:=max(seq(L[i][2],i=1..nops(L)) ): if (f(tt)<1 or dfn(tt)<0 or n>tt or m>tt or type(cc,list(integer))=false) then Q:=Q union {w}: fi: fi: fi: od: T:=T minus Q: [T,nops(T)]: end: anomNum:=proc(Rstart,Bstart,t) local i,j,R,B,C,v,n1,n2,TB,TR,Rpost,Bpost,bnd: if nops(Rstart) < nops(Bstart) then B:=Bstart: R:=Rstart union BuildAnomNum(B,t,t)[1]: B:=Bstart union BuildAnomNum(R,t,2*t+1)[1]: else R:=Rstart: B:=Bstart union BuildAnomNum(R,t,2*t+1)[1]: R:=Rstart union BuildAnomNum(B,t,t)[1]: fi: print('R',R): print('B',B): n1:=-1: n2:=-1: while (n1+n2 <> 0) do TB:=BuildAnomNum(R,t,2*t+1): Bpost:=B union TB[1]: n1:=nops(Bpost minus B): TR:=BuildAnomNum(B,t,t): Rpost:=R union TR[1]: n2:=nops(Rpost minus R): R:=Rpost: B:=Bpost: C:=B intersect R: print('Ris',R): print('Bis',B): if nops(C)>0 then RETURN('Intersection',C): fi: if (n1=0 and n2=0) then bnd:=6*t^3+2*t^2+4*t: print('Missing',{seq(i,i=1..bnd)} minus (R union B)): RETURN('NoMore'): fi: od: 0: end: BuildAnomNum:=proc(X,t,v) local i,j,T,x,y,Q,w,bnd: # b=2t+1, c=t # t is variable in this Build function T:={}: # PROPOSED BOUND: 6t^3+2t^2+4t: for i from 1 to nops(X) do for j from 1 to nops(X) do x:=X[i]: y:=X[j]: T:=T union {simplify(t*x+v*y)}: # ASSUME b,c,t all integers if (y-t*x mod v = 0) then T:=T union {simplify((y-t*x)/v)}: fi: if (y-v*x mod t = 0) then T:=T union {simplify((y-v*x)/t)}: fi: od: if (x mod (v+t) = 0) then T:=T union {simplify((x)/(v+t))}:fi: od: Q:={}: for i from 1 to nops(T) do w:=T[i]: if type(w,integer)=false then Q:=Q union {w}: else bnd:=6*t^3+2*t^2+4*t: if (w<1 or w>bnd) then Q:=Q union {w}: fi: fi: od: T:=T minus Q: [T,nops(T)]: end: diagdan:=proc(Rstart,Bstart,t,b) local i,j,R,B,C,n1,n2,TB,TR,Rpost,Bpost: R:=Rstart: B:=Bstart: if nops(Rstart) < nops(Bstart) then B:=Bstart: R:=Rstart union diagBuild(B,t,b)[1]: B:=Bstart union diagBuild(R,t,b)[1]: else R:=Rstart: B:=Bstart union diagBuild(R,t,b)[1]: R:=Rstart union diagBuild(B,t,b)[1]: fi: print('R',R): print('B',B): n1:=-1: n2:=-1: while (n1+n2 <> 0) do TB:=diagBuild(R,t,b): Bpost:=B union TB[1]: n1:=nops(Bpost minus B): TR:=diagBuild(B,t,b): Rpost:=R union TR[1]: n2:=nops(Rpost minus R): R:=Rpost: B:=Bpost: C:=B intersect R: print('Ris',R): print('Bis',B): if nops(C)>0 then RETURN('Intersection',C): fi: if n1=0 and n2=0 then RETURN('NoMore'): fi: od: 0: end: diagBuild:=proc(X,t,b) local c,cnt,q,i,j,T,x,y,z,ans,m,Q,w,f,g,h: T:={}: m:=t*b^2+(2*t^2+1)*b+t^3: for i from 1 to nops(X) do for j from 1 to nops(X) do x:=X[i]: y:=X[j]: T:=T union {simplify(t*x+b*y)}: if type(b,integer) then if (y-t*x mod b = 0) then T:=T union {(y-t*x)/b}: fi: else if rem(y-t*x,b,b,'ans')=0 then cnt:=0: q:=unapply(ans,b): for z from 0 to t-1 do if type(q(z),integer) then cnt:=cnt+1: fi: od: if cnt=t then T:=T union {simplify(ans)}: fi: fi: fi: if type(t,integer) then c:=[coeffs(expand(y-b*x))]: if type(c,list(integer)) then if (y-b*x mod t = 0) then T:=T union {simplify((y-b*x)/t)}: elif type(y-b*x mod t,quadratic(b)) then q:=unapply(y-b*x,b): if q(t)>0 then cnt:=0: for z from 0 to t-1 do if type(q(z)/t,integer) then cnt:=cnt+1:fi: od: if cnt=t then T:=T union {simplify((y-b*x)/t)}: fi: fi: fi: fi: else if rem(y-b*x,t,t)=0 then T:=T union {simplify((y-b*x)/t)}: fi: fi: od: if rem(x,b+t,b)=0 then T:=T union {simplify((x)/(b+t))}:fi: od: Q:={}: for i from 1 to nops(T) do w:=T[i]: #if w < 1 or w > m then Q:=Q union {w}: fi: if (type(w,integer)=false and type(w,linear(b))=false and type(w,quadratic(b))=false) then Q:=Q union {w}: else f:=unapply(w,b): if type(w,integer) and (w<1 or w>4*t^3+2*t) then Q:=Q union {w}:fi: if type(w,linear(b)) and (f(t)<1 or f(0)>t^3 or diff(w,b)<0 or diff(w,b)>2*t^2+2*t+1) then Q:=Q union {w}: fi: if type(w,quadratic(b)) then h:=diff(diff(w,b),b): if (f(0)>t^3 or f(t)<1 or h>2*t or h<0) then Q:=Q union {w}: else g:=unapply(diff(w,b),b): if g(0)>2*t^2+1 then Q:=Q union {w}:fi: fi: fi: fi: od: T:=T minus Q: [T,nops(T)]: end: