#####################################################################
##PeaceableQueens.txt: Save this file as PeaceableQueens.txt         #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read `PeaceableQueens.txt`<Enter>                  #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Yukun Yao and Doron Zeilberger, Rutgers University      #
#yao at math dot rutgers dot edu; DoronZeil at gmail dot com         #
######################################################################
 
#Created: Nov. 2018 - Jan. 2019


with(plots):

print(`Created:  Nov. 2018-Jan. 2019`):
print(` This is PeaceableQueens.txt `):
print(`It is one of the packages that accompany the article `):
print(` On the Peaceable Queens Question`):
print(`by Yukun Yao and Doron Zeilberger`):

print(`and also available from Yao's and Zeilberger's websites`):
print(``):
print(`Please report bugs to DoronZeil at gmail dot com `):
print(``):
 print(`The most current version of this  package and paper`):
 print(` are  available from`):
 print(`http://sites.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/peaceable.html  .`):

print(`---------------------------------------`):
 print(`For a list of the Story procedures type ezraS();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

print(`---------------------------------------`):
 print(`For a list of the Shape procedures type ezraShape();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

print(`---------------------------------------`):
 print(`For a list of the Plotting procedures type ezraP();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

print(`---------------------------------------`):
 print(`For a list of the symbolic procedures type ezraSy();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

print(`---------------------------------------`):
 print(`For a list of the Continuous procedures type ezraC();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

print(`---------------------------------------`):
 print(`For a list of the Optimization procedures type ezraOp();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

print(`---------------------------------------`):
 print(`For a list of the procedures to find estimated solutions for specific configurations type ezraY();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);.`):
 print(``):
print(`---------------------------------------`):

print(`---------------------------------------`):
 print(`For a list of the Supporting procedures type ezra1();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

print(`---------------------------------------`):
 print(`For a list of the MAIN procedures type ezra();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

with(combinat):

ezraShape:=proc()

if args=NULL then
 print(` The Shape procedures are: GenLeftTri, GenRect, GenRightTri `):
 print(``):

else
ezra(args):
fi:

end:

ezraOp:=proc()

if args=NULL then
 print(` The Optimization procedures are: BestNei, HC , MaxC, BestHexagon `):
 print(``):

else
ezra(args):
fi:

end:

ezraSy:=proc()

if args=NULL then
 print(` The the symbolic procedures are: AvParaTrapS .`):
 print(``):

else
ezra(args):
fi:

end:

ezraC:=proc()

if args=NULL then
 print(` The Continuous procedures are: AreaC, AreaSef, BestGenConfC, BestLeftPentC`):
 print(` ConfCs, GenConfC, GenConfCS,GenConfCa, GenConfCSa, IsGoodParaC, LeftPentCabe, RightPentCdef  `):
 print(` TwoTriaConf`):

else
ezra(args):
fi:

end:

ezraP:=proc()

if args=NULL then
 print(` The Plotting procedures are: AvParaP,AvParaPw, AvParaTrapP, AvParaTrapPd, AvParaTrapPa, AvParaTrapPw,AvTrapP, AvTrapPw, `):
print(` GenConfCaP, GenConfCP, GenConfCSaP, GenConfCSP, PlotSet, PlotWB `):
 print(``):

else
ezra(args):
fi:

end:

ezraS:=proc()

if args=NULL then
 print(` The story procedures are: BestGenConfCv `):
 print(``):

else
ezra(args):
fi:

end:


ezraY:=proc()

if args=NULL then
 print(` The procedures to find estimatied solutions are: Plot2Square, NuA2Square, FindM2Square, Plot2Triangle, NuA2Triangle, FindM2Triangle, PlotSquareTriangle, NuASquareTriangle, FindMSquareTriangle `):
 print(``):

else
ezra(args):
fi:

end:



ezraDep:=proc()
if args=NULL then
 print(` The depercated procedures are: `):
 print(` JubinS,JubinSG,`):
else
ezra(args):
fi:

end:

ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are: `):
print(`AvPara,AvParaTrap,  , AvTrap, BestParTrap, Comp, Imp, IsGood, IsGoodPt, IsGoodPoly, LastB, LeftPentABE,`):
 print(` NuWB, ParaAB, RightPentABC, TrapAB, Hexagon, Rectangle, Triangle, MaxC1`):
 print(``):

else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are:  AreaWB1,AreaWBa, BestAB, BestC, GrS, ImpS , Jubin, NuA `):
 print(` `):


elif nops([args])=1 and op(1,[args])=AreaC then
print(`AreaC(a,b,e,v0,c,d,f): the area of GenConfC(a,b,e,v0,c,d,f) . Try:`):
print(`AreaC(1/4,1/3,1/12,1/2,1/4,1/6,1/12);`):

elif nops([args])=1 and op(1,[args])=AreaSef then
print(`AreaSef(e,f): the area of GenConfCS(e,f);. Try:`):
print(`AreaSef(1/12,1/12);`):


elif nops([args])=1 and op(1,[args])=AreaWB1 then
print(`AreaWB1(a,b,c,d,e): the area of the  the white region and the black region if the white area is the parralelogram`):
print(`[[0,0],[a,a],[a,a+b],[0,b]] where 0<a,b<1 union the trapezoid whose vertices are `):
print(` [[c,0],[c+d,0],[c+d,e+d],[c,e]] where 0<c,d,e<1 , and e+d<c under Regime 1. Try: `):
print(` AreaWB1(a,b,c,d,e); `):
print(`AreaWB1(1/5,1/4,1/2,1/6,1/8);`):

elif nops([args])=1 and op(1,[args])=AvPara then
print(`AvPara(a,b): the set of the two triangles where black queens can redise if the white queens are in the`):
print(`parallelogram whose  vertices are [[0,0],[a,a],[a,a+b],[0,b]] where 0<a,b<1 . Try:`):
print(`AvPara(1/4,1/3);`):



elif nops([args])=1 and op(1,[args])=AvParaP then
print(`AvParaP(a,b): plots AvPara(a,b) (q.v.). Try:`):
print(` AvParaP(1/5,1/4); `):



elif nops([args])=1 and op(1,[args])=AvParaPw then
print(`AvParaPw(a,b): plots AvPara(a,b) (q.v.) and the parallelogram  [[0,0],[a,a],[a,a+b],[0,b]]. Try:`):
print(` AvParaPw(1/5,1/4); `):

elif nops([args])=1 and op(1,[args])=AvParaTrap then
print(`AvParaTrap(a,b,c,d,e): the set of polygons  where black queens can redise if the white queens are in the`):
print(`AvPara(a,b) union AvTrap(c,d,e) . Try: `):
print(`AvParaTrap(1/4,1/3,1/2,1/4,1/6);`):

elif nops([args])=1 and op(1,[args])=AvParaTrapP then
print(` AvParaTrapP(a,b,c,d,e):   plots AvParaTrap(a,b,c,d,e). Try: `):
print(` AvParaTrapP(1/4,1/3,1/2,1/4,1/6);`):

elif nops([args])=1 and op(1,[args])=AvParaTrapPd then
print(`AvParaTrapPd(a,b,c,d,e): plots AvPara(a,b) (red) and AvTrap(c,d.e) (blue). DIRTY VERSION. Try:`):
print(`AvParaTrapPd(1/5,1/4,1/2,1/6,1/8);`):

elif nops([args])=1 and op(1,[args])=AvParaTrapPa then
print(`AvParaTrapPa(a,b,c,d,e,N): An approximation to the available region of AvParaTrapP(a,b,c,d,e) using N. Try: `):
print(`AvParaTrapPa(1/5,1/4,1/2,1/6,1/8,120);`):

elif nops([args])=1 and op(1,[args])=AvParaTrapPw then
print(` AvParaTrapPw(a,b,c,d,e):   plots AvParaTrap(a,b,c,d,e). Together with the white area. Try: `):
print(` AvParaTrapPw(1/4,1/3,1/2,1/4,1/6);`):

elif nops([args])=1 and op(1,[args])=AreaWBa then
print(`AreaWBa(a,b,c,d,e,N): An approximation to the  White area of AvParaTrapP(a,b,c,d,e) followed by the Black area, using N.  Try: `):
print(`AreaWBa(1/5,1/4,1/2,1/6,1/8,120);`):

elif nops([args])=1 and op(1,[args])=AvParaTrapS then
print(`AvParaTrapS(a,b,c,d,e): the set of polygons  where black queens can redise if the white queens are in the`):
print(`for SYMBOLIC a,b,c,d,e`):
print(`AvPara(a,b) union AvTrap(c,d,e) .`):
print(` Try: AvParaTrapS(a,b,c,d,e); `):

elif nops([args])=1 and op(1,[args])=AvTrap then
print(`AvTrap(c,d,e): the set of the four polygons (triangle,trapezoid, trapezoid,triangle), where black queens can redise if the white queens are in the`):
print(`trapezoid  vertices are [[c,0],[c+d,0],[c+d,e+d],[c,e]] where 0<c,d,e<1 , and e+d<c and possibly other conditions. Try:`):
print(`AvTrap(1/2,1/6,1/8);`):

elif nops([args])=1 and op(1,[args])=AvTrapP then
print(`AvTrapP(c,d,e): plots AvTrap(c,d,e) (q.v.). Try:`):
print(` AvTrapP(1/2,1/8,1/8);`):

elif nops([args])=1 and op(1,[args])=AvTrapPw then
print(`AvTrapPw(c,d,e): Like AvTrapP(c,d,e) but also plots the trapezoid [[c,0],[c+d,0],[c+d,e+d],[c,e]]. Try:`):
print(` AvTrapPw(1/2,1/8,1/8);`):



elif nops([args])=1 and op(1,[args])=BestAB then
print(`BestAB(N): the largest ParaAB(A,B) that is good for N. Try:`):
print(`BestAB(8); `):

elif nops([args])=1 and op(1,[args])=BestC then
print(`BestC(S,N): Given a subset S of [0,N-1]x[0,N-1] and a pos. integer N, finds the set of pt outside S with the largest value of`):
print(`NuA(S,N); It returns the record and the set of champions. Try: `):
print(`BestC({[0,0]},8); `):

elif nops([args])=1 and op(1,[args])=BestLeftPentC then
print(`BestLeftPentC(e1): The best choice in LeftPentC(a,b,e) woth e=e1, (the left pentagon  in the Jubin's construction).`):
print(`It returns the White pentagon and left pentagon`):
print(` Try: `):
print(`BestLeftPentC(.309915920); `):

elif nops([args])=1 and op(1,[args])=BestParTrap then
print(`BestParTrap(r1,r2,c1): finds the optimal placement of AreaWB1(a,r1*a,c1,d,d*r2). `):
print(`It returns the optimal configuration [a,r1*a,c1,d,d*r2], followed by the area of the white region. Try: `):
print(`BestParTrap(1,1,1/2); `):

elif nops([args])=1 and op(1,[args])=BestGenConfC then
print(`BestGenConfC(): The best choice in GenConfC(a,b,e,v0,c,d,f) (two pentagons as in Jubin's construction).`):
print(`It returns the two White pentagons and the two Black pentagons that are the best`):
print(`Try: `):
print(`BestGenConfC();`):

elif nops([args])=1 and op(1,[args])=BestGenConfCv then
print(`BestGenConfCv(): Verbose version of BestGenConfC(): `):
print(` The best choice in GenConfC(a,b,e,v0,c,d,f) (two pentagons as in Jubin's construction).`):
print(`It returns the two White pentagons and the two Black pentagons that are the best`):
print(`Try: `):
print(`BestGenConfCv();`):

elif nops([args])=1 and op(1,[args])=BestNei then
print(`BestNei(P,x,x0,eps) inputs a function P of the set of variables x, a current place x0, and a resolution eps, finds`):
print(`the best neighbor. Try:`):
print(` BestNei(x->x^3+3*x^2-9*x,[x],[-1],0.01);`):

elif nops([args])=1 and op(1,[args])=Comp then
print(`Comp(S,N): the complement [0,N-1]x[0,N-1] minus S `):

elif nops([args])=1 and op(1,[args])=ConfCs then
print(`ConfCs(a,b,e,v0,c,d,f): The symmetric configuration, followed by the free variables. Try:`):
print(`ConfCs(a,b,e,v0,c,d,f);`):


elif nops([args])=1 and op(1,[args])=GenConfC then
print(`GenConfC(a,b,e,v0,c,d,f): the generalized configuration with parameters a,b,e,v0,c,d,e where the`):
print(`left pentagon is LeftPentCabe(a,b,e) and the right pentagon is RightPentCdef(v0,c,d,f):`):
print(`It returns the pair [[LeftWhitePentagon,RightWhitePentagon], [LeftBlackPentagon,RightBlackPentagon]]`):
print(`Try: `):
print(` GenConfC(a,b,e,v0,c,d,f); `):
print(`GenConfC(1/4,1/3,1/12,1/2,1/4,1/6,1/12);`):


elif nops([args])=1 and op(1,[args])=GenConfCa then
print(`GenConfCa(a,b,e,v0,c,d,f,N): the white queen in the approximation to an N by N board for , GenConfC(a,b,e,v0,c,d,f) (q.v.) `):
print(`followed by the black queens. Try:`):
print(`GenConfCa(1/4,1/3,1/12,1/2,1/4,1/6,1/12,24);`):

elif nops([args])=1 and op(1,[args])=GenConfCaP then
print(`GenConfCaP(a,b,e,v0,c,d,f,N): Plotting the white queen in the approximation to an N by N board. Try:`):
print(`GenConfCaP(1/4,1/3,1/12,1/2,1/4,1/6,1/12,24);`):

elif nops([args])=1 and op(1,[args])=GenConfCP then
print(`GenConfCP(a,b,e,v0,c,d,f): Plots GenConfC(a,b,e,v0,c,d,f).`):
print(`Try:`):
print(`GenConfCP(1/4,1/3,1/12,1/2,1/4,1/6,1/12);`):

elif nops([args])=1 and op(1,[args])=GenConfCS then
print(`GenConfCS(e,f): The general SYMMETRIC Jubin-like configurations. Try:`):
print(`GenConfCS(1/12,1/12);`):

elif nops([args])=1 and op(1,[args])=GenConfCSa then
print(`GenConfCSa(e,f,N): the white queen in the approximation to an N by N board for , GenConfCS(e,f) (q.v.) `):
print(`followed by the black queens. Try:`):
print(`GenConfCSa(1/12,1/12,24);`):

elif nops([args])=1 and op(1,[args])=GenConfCSaP then
print(`GenConfCSaP(e,f): Plots GenConfCSa(e,f). Try:`):
print(`GenConfCSaP(1/12,1/12,24);`):

elif nops([args])=1 and op(1,[args])=GenConfCSP then
print(`GenConfCSP(e,f): Plots GenConfCS(e,f). Try:`):
print(`GenConfCSP(1/12,1/12);`):


elif nops([args])=1 and op(1,[args])=GenLeftTri then
print(`GenLeftTri(v0,A): The general equilateral triangle with vertex v0 and side of size A where the hypot. is pointing left.`):
print(`Try: `):
print(`GenLeftTri([0,0],30); `):

elif nops([args])=1 and op(1,[args])=GenRect then
print(`GenRect(v0,A,B): The general A by B rectangle with vertex v0. Try:`):
print(`GenRect([0,0],10,10);`):

elif nops([args])=1 and op(1,[args])=GrS then
print(`GrS(N): uses a greedy algorithm to find a good set of peacebale queens, followed by the set of black queens. Try:`):
print(`GrS(8);`):

elif nops([args])=1 and op(1,[args])=HC then
print(`HC(P,x,x0,eps,K) inputs a function P of the set of variables x, a current place x0, and a resolution eps, finds`):
print(`a local maximum with up to K iterations, followed by the value. Try:`):
print(`HC(x->-x^3-3*x^2+9*x,[x],[-0.9],0.01,300);,[x],[-1],0.01);`):
print(` HC((a,b,c)->evalf(BestParTrap(a,b,c))[2][1],[a,b,c],[1,1,0.5],0.01,100);`):
print(` HC((a,b,c)->evalf(BestParTrap(a,b,c))[2][1],[a,b,c],[1.16, 0.82, 0.51],0.001,100);`):

elif nops([args])=1 and op(1,[args])=Imp then
print(`Imp(pt,N): inputs a lattice point pt=[x,y] in [0,N-1]x[0,N-1] and outputs  the set of lattice points attacked by it. Try`):
print(`Inp([0,0],8);`):

elif nops([args])=1 and op(1,[args])=ImpS then
print(`ImpS(S,N): inputs a set of lattice point {pt=[x,y]} [0,N-1]x[0,N-1] and outputs the set of lattice points attacked by them. Try`):
print(`ImpS({[0,0]},8);`):

elif nops([args])=1 and op(1,[args])=IsGood then
print(`IsGood(S,N): Is nops(S)>=NuA(S,N)?`):

elif nops([args])=1 and op(1,[args])=IsGoodPoly then
print(`IsGoodPoly(Po): is the polygon Po good? , i.e. are all its vertices good? Try:`):
print(`IsGoodPoly([[1/2,1/3],[1/2,-1]]);`):

elif nops([args])=1 and op(1,[args])=IsGoodPt then
print(`IsGoodPt(P): is the point P good? Try:`):
print(`IsGoodPt([1/2,1/3]); IsGoodPt([1/2,3/2]);`):

elif nops([args])=1 and op(1,[args])=IsGoodParaC then
print(`IsGoodParaC(a,b): Is the continuous parallelogram whose vertices are [0,0],[a,a],[a,a+b],[0,b] a peaceable one?`):
print(`Try: `):
print(`IsGoodParaC(1/3,1/3);`):

elif nops([args])=1 and op(1,[args])=Jubin then
print(`Jubin(N): Jubin's shape for an N by N board. Try:`):
print(`Jubin(12);`):

elif nops([args])=1 and op(1,[args])=JubinS then
print(`JubinS(N): Jubin's shape before adjustment. Try: `):
print(` JubinS(12); `):

elif nops([args])=1 and op(1,[args])=JubinSG then
print(`JubinSG(N,A,B,C,D): Jubin's generalized shape before adjustment. Try`):
print(`JubinSG(12,12/4,12/3,12/4,12/6);`):

elif nops([args])=1 and op(1,[args])=LastB then
print(`LastB(A,N): the last B for which ParaAB(A,B) is good for N. Try:`):
print(`LastB(4,20);`):

elif nops([args])=1 and op(1,[args])=LeftPentABE then
print(`LeftPentABE(A,B,E): the left pentagon `):
print(`{seq(seq([i,j],j=i..i+B-1),i=0..A-E-1),seq(seq([i,j],j=i..A+B-E-1),i=A-E..A-1)}: Try:`):
print(`LeftPentABE(3,4,1);`):

elif nops([args])=1 and op(1,[args])=LeftPentCabe then
print(`LeftPentCabe(a,b,e): the left-pentagon  continuous version. Try:`):
print(`LeftPentCabe(1/4,1/3,1/12);`):


elif nops([args])=1 and op(1,[args])=MaxC then
print(`MaxC(L,v): given  a list of length 2, L, consisting of polynomials in the list of variables v`):
print(`finds all the extreme points of L[1], subject to the constraint L[1]=L[2], using`):
print(`Largrange multipliers. Try:`):
print(`MaxC([a*b,(1-a-b)^2],[a,b]);`):
print(`MaxC(AreaC(a,b,e,v0,c,d,f),[a,b,e,v0,c,d,f]);`):

elif nops([args])=1 and op(1,[args])=NuA then
print(`NuA(S,N): the number of lattice points in [0,N-1]x[0,N-1] not atttacked by S`):
print(`Try NuA(Jubin(12),12);`):

elif nops([args])=1 and op(1,[args])=NuWB then
print(`NuWB(S,N): [nops(S),NuS(S,N)]`):

elif nops([args])=1 and op(1,[args])=ParaAB then
print(`ParaAB(A,B): the parallelogram {seq(seq([i,i+j],j=0..B),i=0..A)}. Try:`):
print(`ParaAB(5,6);`):

elif nops([args])=1 and op(1,[args])=PlotSet then
print(`PlotSet(Q): plots the set Q.Try:`):
print(`PlotSet(ParaAB(5,6);`):

elif nops([args])=1 and op(1,[args])=PlotWB then
print(`PlotWB(Q,N): plots the set Q (in red) and Comp(ImpS(S,N),N) in blue. Try: `):
print(`PlotWB(ParaAB(5,6),20);`):

elif nops([args])=1 and op(1,[args])=RightPentABC then
print(`RightPentABC(V0,A,B,C): the right-pentagon with left vertex V0`):
print(`V0+{seq(seq([i,j],j=0..B+i),i=0..A-C), seq(seq([i,j],j=0..B+2*(A-C)-i,i=A-C+1..A-1)}`):
print(`Try: `):
print(`RightPentABC([6,0],3,2,1);`):

elif nops([args])=1 and op(1,[args])=RightPentCdef then
print(`RightPentCdef(v0,d,e,f): the right-pentagon with left vertex [v0,0]. The continuous analog of`):
print(`RightPentABC([V0,0],D,E,F); Try:`):
print(`RightPentCdef(1/2,1/4,1/6,1/12);`):

elif nops([args])=1 and op(1,[args])=TrapAB then
print(`TrapAB(V0,A,B): the trapezoid with bottom left vertex V`):
print(`V0+{seq(seq([i,j],j=0..B+i),i=0..A)}:`):
print(` Try: `):
print(` TrapAB([11,0],5,4);`):


elif nops([args])=1 and op(1,[args])=TwoTriaConf then
print(`TwoTriaConf(a,b): inputs numeric or symbolic a and b (such that 0<a<b<1)`):
print(`and output a list consistinbg of the two triangles with vertices`):
print(` (i)[[0,0],[a,0],[a,a]], [[b,0],[a+b,0],[a,+b,a]] (White area) `):
print(` (ii) The quadrilateral , triangle and parallelogram `):
print(` [[a,1],[a,a+b],[b,2*a],[b,1]], [[a+b,1],[a+b,a+b],[1,1]], [[a+b,b],]a+b,a],[1,1-b],[1,1-a]] `):
print(` (iii) The white area (a^2) `):
print(` (iv) The black area. try `):
print(` TwoTriaConf(a,b); `):
print(`TwoTriaConf(1/3,1/2);`):


elif nops([args])=1 and op(1,[args])=Plot2Square then
print(`Plot2Square(a,s,N): plot two squares of edge length a, the left one’s lower left corner is [0,0], the right one’s is [s,0] on N by N chess board. Try:`):
print(`Plot2Square(25, 50, 100);`):

elif nops([args])=1 and op(1,[args])=Plot2Triangle then
print(`Plot2Triangle(a,s,N): plot two triangles of edge length a (bottom and right), the left one’s lower left corner is [0,0], the right one’s is [s,0] on N by N chess board.Try:`):
print(`Plot2Triangle(25, 50, 100);`):

elif nops([args])=1 and op(1,[args])=PlotSquareTriangle then
print(`PlotSquareTriangle(a,s,N): plot one square and one triangle of edge length a, the left square’s lower left corner is [0,0], the right triangle’s is [s,0] on N by N chess board. Try:`):
print(`PlotSquareTriangle(25, 50, 100);`):

elif nops([args])=1 and op(1,[args])=NuA2Square then
print(`NuA2Square(a,s,N): the non-attacking area in the two identical squares configuration. Meaning of parameters are the same withPlot2Square. Try:`):
print(`NuA2Square(25, 50, 100);`):

elif nops([args])=1 and op(1,[args])=NuA2Triangle then
print(`NuA2Triangle(a,s,N): the non-attacking area in the two identical triangles configuration. Meaning of parameters are the same withPlot2Triangle. Try:`):
print(`NuA2Triangle(25, 50, 100);`):

elif nops([args])=1 and op(1,[args])=NuASquareTriangle then
print(`NuASquareTriangle(a,s,N): the non-attacking area in the square + triangle (the same side length) configuration. Meaning of parameters are the same with PlotSquareTriangle. Try:`):
print(`NuASquareTriangle(25, 50, 100);`):

elif nops([args])=1 and op(1,[args])=FindM2Square then
print(`FindM2Square(s,N): for different s which is the x-coordinate of the right square’s lower left corner, find out the interval [a, a+1] and [a/N, (a+1)/N] such that at a and a+1 of the two ratios white/black one is >1 and the other is <1. Try:`):
print(`FindM2Square(50, 100);`):

elif nops([args])=1 and op(1,[args])=FindM2Triangle then
print(`FindM2Triangle(s,N): for different s which is the x-coordinate of the right square’s lower left corner, find out the interval [a, a+1] and [a/N, (a+1)/N] such that at a and a+1 of the two ratios white/black one is >1 and the other is <1. Try:`):
print(`FindM2Triangle(50, 100);`):

elif nops([args])=1 and op(1,[args])=FindMSquareTriangle then
print(`FindMSquareTriangle(s,N): for different s which is the x-coordinate of the right triangle’s lower left corner, find out the interval [a, a+1] and [a/N, (a+1)/N] such that at a and a+1 of the two ratios white/black one is >1 and the other is <1. Try:`):
print(`FindMSquareTriangle(50, 100);`):

elif nops([args])=1 and op(1,[args])=BestHexagon then
print(`BestHexagon(): the hexagon configuration in the lower left corner which maximizes the number of queens. Try:`):
print(`BestHexagon();`):

elif nops([args])=1 and op(1,[args])=Hexagon then
print(`Hexagon(a,b,c,d): the hexagon in lower left corner. Try:`):
print(`Hexagon(0.2,0.2,0.2,0.2);`):

elif nops([args])=1 and op(1,[args])=Rectangle then
print(`Rectangle(s,a,b): a rectangle whose lower left point is [s,0] with width a and height b. Try:`):
print(`Rectangle(0.2, 0.3, 0.2);`):

elif nops([args])=1 and op(1,[args])=Triangle then
print(`Triangle(s,a): a triangle whose lower left point is [s,0] and upper right point is [s+a-1, a-1] with width a and height a. Try:`):
print(`Triangle(0.2,0.3);`):


else
print(`There is no ezra for`,args):
fi:
 
end:



#Imp(pt,N): inputs a lattice point pt=[x,y] in [0,N-1]x[0,N-1] and outputs  the set of lattice points attacked by it. Try
#Inp([0,0],8);
Imp:=proc(pt,N) local x,y,i:
x:=pt[1]: y:=pt[2]:
{
seq([x,i],i=0..N-1),seq([i,y],i=0..N-1),
seq([x+i,y+i],i=0..min(N-1-x,N-1-y)),seq([x-i,y-i],i=0..min(x,y)),
seq([x+i,y+i],i=0..min(N-1-x,N-1-y)),seq([x-i,y-i],i=0..min(x,y)),
seq([x+i,y-i],i=0..min(N-1-x,y)),seq([x-i,y+i],i=0..min(x,N-1-y))
}:
end:


#ImpS(S,N): inputs a set of lattice point {pt=[x,y]} [0,N-1]x[0,N-1] and outputs the set of lattice points attacked by them. Try
#ImpS({[0,0]},8);
ImpS:=proc(S,N) local pt:

{seq(op(Imp(pt,N)), pt in S)}:

end:

#Comp(S,N): the complement [0,N-1]x[0,N-1] minus S
Comp:=proc(S,N) local i,j:

{seq(seq([i,j],i=0..N-1),j=0..N-1)} minus S:

end:


#NuA(S,N): the number of lattice points in [0,N-1]x[0,N-1] not atttacked by S
NuA:=proc(S,N):
nops(Comp(ImpS(S,N),N)):
end:

#NuWB(S,N): [nops(S),NuS(S,N)]
NuWB:=proc(S,N):
[nops(S),NuA(S,N)]:
end:


#BestC(S,N): Given a subset S of [0,N-1]x[0,N-1] and a pos. integer N, finds the set of pt outside S with the largest value of
#NuA(S,N); It returns the record and the set of champions. Try
#BestC({[0,0]},8);
BestC:=proc(S,N) local i,j,gu,alufim,si,pt:
gu:={seq(seq([i,j],j=0..N-1),i=0..N-1)} minus S:
alufim:={}:

if nops(gu)=0 then
 RETURN({}):
fi:
gu:=convert(gu,list):
alufim:={gu[1]}:
si:=NuA(S union {gu[1]},N):

for i from 2 to nops(gu) do
 pt:=gu[i]:

 if NuA(S union {pt},N)> si then
  si:=NuA(S union {pt},N):
   alufim:={pt}:
 elif NuA(S union {pt},N)= si then
    alufim:= alufim union {pt}:
 fi:
od:
si,alufim:
end:

#GrS(N): uses a greedy algorithm to find a good set of peacebale queens, followed by the set of black queens. Try:
#GrS(8);
GrS:=proc(N) local S,S1:
S:={[0,0]}:

while nops(S)<=NuA(S,N) do
 S1:= S union {BestC(S,N)[2][1]}:
 if nops(S1)>NuA(S1,N) then
  RETURN([S,Comp(ImpS(S,N),N)]):
 fi:
S:=S1:
od:

end:

#ParaAB(A,B): the parallelogram {seq(seq([i,i+j],j=0..B),i=0..A)}. Try:
#ParaAB(5,6);
ParaAB:=proc(A,B) local i,j:
{seq(seq([i,j],j=i..i+B-1),i=0..A-1)}:
end:

#IsGood(S,N): Is nops(S)>=NuA(S,N)?
IsGood:=proc(S,N)
evalb(nops(S)<=NuA(S,N)):
end:

#LastB(A,N): the last B for which ParaAB(A,B) is good for N
LastB:=proc(A,N) local B:

for B from 1 while IsGood(ParaAB(A,B),N) do od:
B-1:
end:

#BestAB(N): the largest ParaAB(A,B) that is good for N
BestAB:=proc(N) local A,B,si,aluf,hal:
 B:=LastB(1,N):
 aluf:=0:
 si:=nops(ParaAB(1,B)):
 for A from 2 to N/2 do
   B:=LastB(A,N):
   hal:=nops(ParaAB(A,B)):
   if hal>si then
      si:=hal:
      aluf:=A:
   fi:
 od:
[aluf,LastB(aluf,N)]:
end:

#TrapAB(V0,A,B): the trapezoid with bottom left vertex V
#V0+{seq(seq([i,j],j=0..B+i),i=0..A)}:
#Try:
#TrapAB([11,0],5,4);
TrapAB:=proc(V0,A,B) local i,j:
{seq(seq([V0[1]+i,V0[2]+j],j=0..B+i-1),i=0..A-1)}:
end:

#JubinS(N): Jubin's shape before adjustment. Try
#JubinS(12);
JubinS:=proc(N)
if N mod 12<>0 then
 print(`N must be a multiple of 12`):
 RETURN(FAIL):
fi:
ParaAB(N/4,N/3)  union TrapAB([N/2,0],N/4,N/6):
end:



#JubinSG(N,A,B,C,D): Jubin's generalized shape before adjustment. Try
#JubinSG(12,12/4,12/3,12/4,12/6);
JubinSG:=proc(N,A,B,C,D)

ParaAB(A,B)  union TrapAB([trunc(N/2),0],C,D):
end:


#IsGoodParaC(a,b): Is the continuous parallelogram whose vertices are [0,0],[a,a],[a,a+b],[0,b] a peaceable one?
#Try
#IsGoodParaC(1/3,1/3);
IsGoodParaC:=proc(a,b)
evalb(0<=a and a<=1 and 0<=b and b<=1 and 2*a+b<=1 and a*b<=(1-a-b)^2):
end:

#PlotSet(Q): plots the set Q
PlotSet:=proc(Q)
#plot(Q,style=point,symbol=CIRCLE,color=black,axes=none):
plot(Q,style=point,symbol=CROSS,color=black,scaling=constrained):
end:

#PlotWB(Q,N): plots the set Q (in red) and Comp(ImpS(S,N),N) in blue. Try
#PlotWB(ParaAB(5,6),20);
PlotWB:=proc(Q,N) local Q1,pic:
Q1:=Comp(ImpS(Q,N),N):
pic:=plot(Q,style=point,symbol=CROSS,color=red,scaling=constrained),plot(Q1,style=point,symbol=CROSS,color=blue,scaling=constrained):
display(pic);
end:


#AvPara(a,b): the set of the two triangles where black queens can redise if the white queens are in the
#parallelogram whose  vertices are [[0,0],[a,a],[a,a+b],[0,b]] where 0<a,b<1 . Try:
#AvPara(1/4,1/3);
AvPara:=proc(a,b) local gu,mu,i:
mu:={[[a,a+b], [1-b,1],[a,1]],
[[a+b,a+b],[1,a+b],[1,1]]}:
gu:={}:
for i from 1 to nops(mu) do
 if IsGoodPoly(mu[i]) then
  gu:=gu union {mu[i]}:
 fi
od:
gu:

end:

#AvParaP(a,b): plots AvPara(a,b) (q.v.). Try:AvParaP(1/5,1/4);
AvParaP:=proc(a,b) local gu,pic,gu1:
gu:=AvPara(a,b):
if gu={} then
 RETURN(FAIL):
fi:

gu1:=gu[1]:
gu1:=[op(gu1),gu1[1]]:
pic:=plot(gu1,color=red,scaling=constrained):

if nops(gu)=2 then
gu1:=gu[2]:
gu1:=[op(gu1),gu1[1]]:
pic:=pic,plot(gu1,color=red,scaling=constrained):
pic:=pic,plot({[0,0]}	,style=point,symbol=CROSS):
fi:

display(pic);
end:

#AvParaPw(a,b): plots AvPara(a,b) and the parallelogram [[0,0],[a,a],[a,a+b],[0,b]] where 0<a,b<1 . Try:
#AvParaPw(1/5,1/4);
AvParaPw:=proc(a,b) local gu,pic,gu1:
gu:=AvPara(a,b):
 if gu={} then
   RETURN(FAIL):
 fi:

gu1:=gu[1]:
gu1:=[op(gu1),gu1[1]]:
pic:=plot(gu1,color=red,scaling=constrained):
if nops(gu)=2 then
gu1:=gu[2]:
gu1:=[op(gu1),gu1[1]]:
pic:=pic,plot(gu1,color=red,scaling=constrained):
pic:=pic,plot({[0,0]}	,style=point,symbol=CROSS):
fi:
pic:=pic,plot([[0,0],[a,a],[a,a+b],[0,b],[0,0]],color=blue):
display(pic);
end:





#AvTrap(c,d,e): the set of the four polygons (triangle,trapezoid, trapezoid,triangle), where black queens can redise if the white queens are in the
#trapezoid  vertices are [[c,0],[c+d,0],[c+d,e+d],[c,e]] where 0<c,d,e<1 , and e+d<c and possibly other conditions. Try:
#AvTrap(1/2,1/6,1/8);
AvTrap:=proc(c,d,e) local gu,mu,i:

mu:=[
[[0,e+d], [0,c],[c-e-d,e+d]],
[[0,c+e+2*d],[0,1],[c,1],[c,e+2*d]],
[[c+d,1], [1,1],[1,1-c+e],[c+d,e+d]],
[[c+e+2*d,e+d],[1,1-c-d],[1,e+d]]
]:

gu:=[]:
for i from 1 to nops(mu) do
  gu:=[op(gu),TrimPoly(mu[i])]:
od:
gu:



end:

#AvTrapP(c,d,e): plots AvTrap(c,d,e) (q.v.). Try:AvTrapP(1/2,1/8,1/8);
AvTrapP:=proc(c,d,e) local gu,pic,gu1,i:


gu:=AvTrap(c,d,e):
gu1:=gu[1]:
gu1:=[op(gu1),gu1[1]]:
pic:=plot(gu1,color=blue,scaling=constrained):

for i from 2 to nops(gu) do
gu1:=gu[i]:
gu1:=[op(gu1),gu1[1]]:
pic:=pic,plot(gu1,color=blue,scaling=constrained):
od:
pic:=pic,plot({[0,0]}	,style=point,symbol=CROSS):

display(pic);

end:

#AvTrapPw(c,d,e): plots AvTrap(c,d,e) (q.v.) but also the trapezoid. 
#[[c,0],[c+d,0],[c+d,e+d],[c,e]]
#Try:AvTrapPw(1/2,1/8,1/8);
AvTrapPw:=proc(c,d,e) local gu,pic,gu1,i:
#if c+e+2*d>1 then
# print(`Not yet implenented`):
# RETURN(FAIL):
#fi:

gu:=AvTrap(c,d,e):
gu1:=gu[1]:
gu1:=[op(gu1),gu1[1]]:
pic:=plot(gu1,color=red,scaling=constrained):

for i from 2 to nops(gu) do
gu1:=gu[i]:
gu1:=[op(gu1),gu1[1]]:
pic:=pic,plot(gu1,color=red,scaling=constrained):
od:
pic:=pic,plot({[0,0]}	,style=point,symbol=CROSS):
pic:=pic, plot([[c,0],[c+d,0],[c+d,e+d],[c,e],[c,0]],color=blue):
display(pic);
end:

#IsGoodPt(P): is the point P good? Try:
#IsGoodPt([1/2,1/3]); IsGoodPt([1/2,3/2]);

IsGoodPt:=proc(P)
if not (type(P[1],numeric) or type(P[2],numeric) ) then
 RETURN(true):
fi:

if type(P[1],numeric) then
 if P[1]<0 or P[1]>1 then
   RETURN(false):
 fi:
fi:

if type(P[2],numeric) then
 if P[2]<0 or P[2]>1 then
   RETURN(false):
 fi:
fi:
true:
end:




#IsGoodPoly(Po): is the polygon Po good? , i.e. are all its vertices good? Try:
#IsGoodPoly([[1/2,1/3],[1/2,-1]]);

IsGoodPoly:=proc(Po) local i:
evalb({seq(IsGoodPt(Po[i]),i=1..nops(Po))}={true}):
end:


#AvParaTrapPd(a,b,c,d,e): plots AvPara(a,b) (red) and AvTrap(c,d.e) (blue). Dirty Version Try
#AvParaTrapPd(1/5,1/4,1/2,1/6,1/8);
AvParaTrapPd:=proc(a,b,c,d,e) local gu,pic,gu1,i:
gu:=AvPara(a,b):
if gu={} then
 RETURN(FAIL):
fi:

gu1:=gu[1]:
gu1:=[op(gu1),gu1[1]]:
pic:=plot(gu1,color=red,scaling=constrained):

if nops(gu)=2 then
gu1:=gu[2]:
gu1:=[op(gu1),gu1[1]]:
pic:=pic,plot(gu1,color=red,scaling=constrained):
fi:

gu:=AvTrap(c,d,e):


gu1:=gu[1]:
gu1:=[op(gu1),gu1[1]]:
pic:=pic,plot(gu1,color=blue,scaling=constrained):

for i from 2 to nops(gu) do
gu1:=gu[i]:
gu1:=[op(gu1),gu1[1]]:
pic:=pic,plot(gu1,color=blue,scaling=constrained):
od:
pic:=pic,plot({[0,0]}	,style=point,symbol=CROSS):


display(pic);
end:

#AvParaTrapPa(a,b,c,d,e,N): An approximation to the available region of AvParaTrapP(a,b,c,d,e) using N. Try
#AvParaTrapPa(1/5,1/4,1/2,1/6,1/8,120);
AvParaTrapPa:=proc(a,b,c,d,e,N) local gu:
gu:=TrapAB([trunc(N*c),0],trunc(d*N),trunc(e*N)) union ParaAB(trunc(a*N),trunc(b*N)):
PlotWB(gu,N):
end:





#AvParaTrap(a,b,c,d,e): the set of polygons  where black queens can redise if the white queens are in the
#AvPara(a,b) union AvTrap(c,d,e) .
#Try: AvParaTrap(1/5,1/4,1/2,1/6,1/6);
AvParaTrap:=proc(a,b,c,d,e) local gu,mu:
gu:=[]:

mu:=TrimPoly([[a,e+c+2*d-a],[a,1],[c,1],[c,c+b],[(e+c-b)/2+d, (e+c+b)/2+d]]):

if nops(mu)>2 then
 gu:=[op(gu),mu]:
fi:

mu:=TrimPoly([[c+d,b+c+d],[c+d,1],[1-b,1]]):

if nops(mu)>2 then
 gu:=[op(gu),mu]:
fi:

 mu:=TrimPoly([[c+d,a+b],[c+d,c+d],[1,1],[1,1-(c-e)],[a+b+c-e,a+b],[c+d,a+b] ]):

if nops(mu)>2 then
 gu:=[op(gu),mu]:
fi:

gu:

end:


#AvParaTrapP(a,b,c,d,e):   plots AvParaTrap(a,b,c,d,e). Try:
# AvParaTrapP(1/5,1/4,1/2,1/6,1/6);
AvParaTrapP:=proc(a,b,c,d,e) local gu,pic,gu1,i:

gu:=AvParaTrap(a,b,c,d,e):


gu1:=gu[1]:
gu1:=[op(gu1),gu1[1]]:
pic:=plot(gu1,color=blue,scaling=constrained):

for i from 2 to nops(gu) do
gu1:=gu[i]:
gu1:=[op(gu1),gu1[1]]:
pic:=pic,plot(gu1,color=blue,scaling=constrained):
od:
pic:=pic,plot({[0,0]}	,style=point,symbol=CROSS):

display(pic);

end:


#TrimPoly(Po): kicks from Poly all the illegal points. Try:
#TrimPoly([[1/2,1/3],[1/2,-1]]);

TrimPoly:=proc(Po) local gu,i:
gu:=[]:

for i from 1 to nops(Po) do
if IsGoodPt(Po[i]) then
 gu:=[op(gu),Po[i]]:
fi:
od:
gu:
end:


#AvParaTrapS(a,b,c,d,e): the set of polygons  where black queens can redise if the white queens are in the
#for SYMBOLIC a,b,c,d,e
#AvPara(a,b) union AvTrap(c,d,e) .
#Try: AvParaTrapS(a,b,c,d,e);
AvParaTrapS:=proc(a,b,c,d,e) :

 [[[a,e+c+2*d-a],[a,1],[c,1],[c,c+b],[(e+c-b)/2+d, (e+c+b)/2+d]],[[c+d,b+c+d],[c+d,1],[1-b,1]],
[[c+d,a+b],[c+d,c+d],[1,1],[1,1-(c-e)],[a+b+c-e,a+b],[c+d,a+b] ]]:


end:


#AreaWB1(a,b,c,d,e): the area of the  the white region and the black region if the white area is the parralelogram
#[[0,0],[a,a],[a,a+b],[0,b]] where 0<a,b<1 union the trapezoid whose vertices are
#[[c,0],[c+d,0],[c+d,e+d],[c,e]] where 0<c,d,e<1 , and e+d<c under Regime 1. Try:
#AreaWB1(a,b,c,d,e);
AreaWB1:=proc(a,b,c,d,e) local guL,guR:

guL:=(1-a-b)^2/2-(e+c+2*d-2*a-b)*((e+c+2*d-b)/2-a)/2-(1-c-b)^2/2:

guR:=(1-a-b)^2/2-(c+d-a-b)^2/2-(1+e-c-a-b)^2/2:

[a*b+d*e+d^2/2,expand(guL+guR)]:

end:

#AreaWBa(a,b,c,d,e,N): An approximation to the  White area of AvParaTrapP(a,b,c,d,e) followed by the Black area, using N.  Try
#AreaWBa(1/5,1/4,1/2,1/6,1/8,120);
AreaWBa:=proc(a,b,c,d,e,N) local gu:

gu:=TrapAB([trunc(N*c),0],trunc(d*N),trunc(e*N)) union ParaAB(trunc(a*N),trunc(b*N)):

evalf([nops(gu)/N^2, NuA(gu,N)/N^2]):
end:



#BestParTrap(r1,r2,c1): finds the optimal placement of AreaWB1(a,r1*a,c1,d,d*r2). 
#It returns the pair: [optimal configuration [a,r1*a,c1,d,d*r2], the area of the white region]. Try
#BestParTrap(1,1,1/2); 
BestParTrap:=proc(r1,r2,c1) local a,d,hal,gu,lu,a1,mu,d1:
hal:=[a,r1*a,c1,d,d*r2]:
gu:=AreaWB1(op(hal)):
lu:=solve(gu[1]-gu[2],d):
if subs(a=0.2,lu[1])>0 then
 lu:=lu[1]:
else
lu:=lu[2]:
fi:
lu:

mu:=subs(d=lu,gu[1]):
a1:=solve(diff(mu,a),a):

if evalf(a1[1])<0 then
 a1:=a1[2]:
else
 a1:=a1[1]:
fi:

d1:=subs(a=a1,lu):

[subs({a=a1,d=d1},hal),subs({a=a1,d=d1},gu)]:

end:



ez:=proc() :print(`BestNei(x->x^3+3*x^2-9*x,[x],[-1],0.01);, HC(x->x^3+3*x^2-9*x,[x],[-1],0.01);`): end:

#BestNei(P,x,x0,eps) inputs a function P of the set of variables x, a current place x0, and a resolution eps, finds
#the best neighbor. Try:
# BestNei(x->x^3+3*x^2-9*x,[x],[-1],0.01);
BestNei:=proc(P,x,x0,eps) local gu,i,i1,aluf,si,gu1,hal:

if nops(x)<>nops(x0) then
print(`Bad input`):
RETURN(FAIL):
fi:


gu:={seq([seq(x0[i1],i1=1..i-1),x0[i]-eps,seq(x0[i1],i1=i+1..nops(x0))],i=1..nops(x0)),
seq([seq(x0[i1],i1=1..i-1),x0[i]+eps,seq(x0[i1],i1=i+1..nops(x0))],i=1..nops(x0))}:

si:=P(op(x0)):
aluf:=x0:

for gu1 in gu do
hal:=P(op(gu1)):
 
if hal>si then
 aluf:=gu1:
 si:=hal:
fi:
od:

aluf,si:

end:


#HC(P,x,x0,eps,K) applies BestNei(P,x,x0,eps) utp to K times
HC:=proc(P,x,x0,eps,K) local x1,x2,x3,i:
x1:=x0:
x2:=BestNei(P,x,x1,eps)[1]:

for i from 1 to K while x1<>x2 do
 x3:=BestNei(P,x,x2,eps)[1]:
 x1:=x2:
 x2:=x3:
od:

if i=K+1 then
 RETURN(FAIL,x3):
else
 RETURN(x2):
fi:

end:













#AvParaTrapPw(a,b,c,d,e):   plots AvParaTrap(a,b,c,d,e) with the parallelogram and trapezoid. Try:
# AvParaTrapPw(1/5,1/4,1/2,1/6,1/6);
AvParaTrapPw:=proc(a,b,c,d,e) local gu,pic,gu1,i:

gu:=AvParaTrap(a,b,c,d,e):


gu1:=gu[1]:
gu1:=[op(gu1),gu1[1]]:
pic:=plot(gu1,color=blue,scaling=constrained):

for i from 2 to nops(gu) do
gu1:=gu[i]:
gu1:=[op(gu1),gu1[1]]:
pic:=pic,plot(gu1,color=blue,scaling=constrained):
od:
pic:=pic,plot({[0,0]}	,style=point,symbol=CROSS):
pic:=pic,plot([[0,0],[a,a],[a,a+b],[0,b],[0,0]],color=red):
pic:=pic, plot([[c,0],[c+d,0],[c+d,e+d],[c,e],[c,0]],color=red):

display(pic);

end:


#LeftPentABE(A,B,E): the left-pentagon   
#{seq(seq([i,j],j=i..i+B-1),i=0..A-E-1),seq(seq([i,j],j=i..A+B-E-1),i=A-E..A)}: Try:
#LeftPentABE(3,4,1);
LeftPentABE:=proc(A,B,E) local i,j:
{seq(seq([i,j],j=i..i+B-1),i=0..A-E-1),seq(seq([i,j],j=i..A+B-E-1),i=A-E..A-1)}:
end:

#RightPentABC(V0,A,B,C): the right-pentagon with left vertex V0
#V0+{seq(seq([i,j],j=0..B+i),i=0..A-C), seq(seq([i,j],j=0..B+2*(A-C)-i,i=A-C+1..A)}
#Try:
#RightPentABC([6,0],3,2,1);
RightPentABC:=proc(V0,A,B,C) local i,j:
{seq(seq([V0[1]+i,V0[2]+j],j=0..B+i-1),i=0..A-C-1), seq(seq([V0[1]+i,V0[2]+j],j=0..B+2*(A-C-1)-i-1),i=A-C..A) }:
end:


#Jubin(N): Jubin's shape for an N by N board. Try:
#Jubin(12);
Jubin:=proc(N)
if N mod 12<>0 then
 print(`N must be a multiple of 12`):
 RETURN(FAIL):
fi:
LeftPentABE(N/4,N/3,N/12)  union RightPentABC([N/2,0],N/4,N/6,N/12):
end:


#LeftPentCabe(a,b,e): the left-pentagon  continuous version. Try:
#LeftPentCabe(1/4,1/3,1/12);
LeftPentCabe:=proc(a,b,e):
[[0,0],[a,a],[a,a+b-e],[a-e,a+b-e],[0,b],[0,0]]:
end:


#RightPentCdef(v0,c,d,f): the right-pentagon with left vertex [v0,0]. The continuous analog of
#RightPentABC([V0,0],C,D,F); Try:
#RightPentCdef(1/2,1/4,1/6,1/12);
RightPentCdef:=proc(v0,c,d,f):
[[v0,0],[v0+c,0],[v0+c,c-2*f+d],[v0+c-f,c-f+d],[v0,d],[v0,0]]:
end:


#GenConfC(a,b,e,v0,c,d,f): the generalized configuration with parameters a,b,e,v0,c,d,e where the
#left pentagon is LeftPentCabe(a,b,e) and the right pentagon is RightPentCdef(v0,c,d,f):
#It returns the pair [[LeftWhitePentagon,RightWhitePentagon], [LeftBlackPentagon,RightBlackPentagon]]
#Try:
#GenConfC(a,b,e,v0,c,d,f);
#GenConfC(1/4,1/3,1/12,1/2,1/4,1/6,1/12);
GenConfC:=proc(a,b,e,v0,c,d,f) local gu1,gu2,mu1,mu2,i1:
gu1:=LeftPentCabe(a,b,e):
gu2:=RightPentCdef(v0,c,d,f):

mu1:=[[v0,1],[v0,v0+b],[(v0+d-b)/2+c-f,(v0+d+b)/2+c-f],[a,v0+2*c-2*f+d-a],[a,1],[v0,1]]:
mu1:=[seq(mu1[7-i1],i1=1..6)]:
mu2:=[[1,1],[1,1+d-v0],[a+b-e+v0-d,a+b-e],[v0+c,a+b-e],[v0+c,v0+c],[1,1]]:
mu2:=[seq(mu2[7-i1],i1=1..6)]:
[[gu1,gu2],[mu1,mu2]]:
end:


#GenConfCP(a,b,e,v0,c,d,f): Plots GenConfC(a,b,e,v0,c,d,f).
#Try:
#GenConfCP(1/4,1/3,1/12,1/2,1/4,1/6,1/12);
GenConfCP:=proc(a,b,e,v0,c,d,f) local gu,pic:
gu:=GenConfC(a,b,e,v0,c,d,f):

pic:=plot({op(gu[1])},color=red,scaling=constrained):
pic:=pic,plot({op(gu[2])},color=blue,scaling=constrained):
display(pic);
end:

#GenConfCPlabeled(): Labelled version of GenConfCP for GenConfCP(1/4,1/3,1/12,1/2,1/4,1/6,1/12), and the length of each sides is labelled with parameters. 
GenConfCPlabeled:=proc() local P,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,t11,t12,t13,t14,t15,t16,t17,t18,t19:
P:=GenConfCP(1/4,1/3,1/12,1/2,1/4,1/6,1/12):
t1:= textplot([1/8,1/8, sqrt(2)*a], color = "DarkGray", align=right,thickness=0):
t2:= textplot([1/4,3/8, b-e], color = "DarkGray", align=right):
t3:= textplot([5/24,1/2, e], color = "DarkGray", align=above):
t4:= textplot([1/12,5/12, sqrt(2)*(a-e)], color = "DarkGray", align=above):
t5:= textplot([0,1/6, b], color = "DarkGray", align=right):
t6:= textplot([1/2,0, [g,0]], color = "DarkGray", align=below):
t7:= textplot([5/8,0,c], color = "DarkGray", align=above):
t8:= textplot([3/4,1/8,c-2*f+d], color = "DarkGray", align=right):
t9:= textplot([17/24,7/24,sqrt(2)*f], color = "DarkGray", align=right):
t10:= textplot([7/12,1/4, sqrt(2)*(c-f)], color = "DarkGray", align=left):
t11:= textplot([1/2,1/12, d], color = "DarkGray", align=left):
t12:=implicitplot(x=y, x=0.25..0.75, y=0.25..0.75, linestyle=dash, thickness=0, color=“Purple”, transparency=0.5):
t13:=implicitplot(x=1-y, x=1/3..2/3, y=1/3..2/3, linestyle=dash, thickness=0, color=“Purple”,transparency=0.5):
t14:=implicitplot(x=0.25, x=0.24..0.26, y=0.5..0.75, linestyle=dash, thickness=0, color=“Purple”,transparency=0.5):
t15:=implicitplot(x=0.75, x=0.74..0.76, y=0.25..0.5, linestyle=dash, thickness=0, color=“Purple”,transparency=0.5):
t16:=implicitplot(y=0.5, x=0.25..0.75, y=0.49..0.51, linestyle=dash, thickness=0, color=“Purple”,transparency=0.5):
t17:=implicitplot(x=0.5, x=0.49..0.51, y=1/6..5/6, linestyle=dash, thickness=0, color=“Purple”,transparency=0.5):
t18:=implicitplot(y-x=1/3, x=1/6..1/3, y=1/2..2/3, linestyle=dash, thickness=0, color=“Purple”,transparency=0.5):
t19:=implicitplot(y-x=-1/3, x=2/3..5/6, y=1/3..1/2, linestyle=dash, thickness=0, color=“Purple”,transparency=0.5):

display({P,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,t11,t12,t13,t14,t15,t16,t17,t18,t19}):
end:







#GenConfCa(a,b,e,v0,c,d,f,N): the white queen in the approximation to an N by N board of,
#followed by the black queens. Try:
#GenConfCa(1/4,1/3,1/12,1/2,1/4,1/6,1/12,24);
GenConfCa:=proc(a,b,e,v0,c,d,f,N) local gu,gu1,gu2,mu:
gu1:=LeftPentABE(trunc(N*a),trunc(N*b),trunc(N*e) ):
gu2:=RightPentABC([trunc(v0*N),0],trunc(c*N),trunc(d*N),trunc(f*N)):
gu:=gu1 union gu2:
mu:=Comp(ImpS(gu,N),N):
[gu,mu]:
end:

#GenConfCaP(a,b,e,v0,c,d,f,N): Plotting the white queen in the approximation to an N by N board. Try:
#GenConfCaP(1/4,1/3,1/12,1/2,1/4,1/6,1/12,24);
GenConfCaP:=proc(a,b,e,v0,c,d,f,N) local gu,gu1,gu2:
gu1:=LeftPentABE(trunc(N*a),trunc(N*b),trunc(N*e) ):
gu2:=RightPentABC([trunc(v0*N),0],trunc(c*N),trunc(d*N),trunc(f*N)):
gu:=gu1 union gu2:
PlotWB(gu,N):
end:



#ConfCs(a,b,e,v0,c,d,f): The symmetric configuration, followed by the free variables. Try:
#ConfCs(a,b,e,v0,c,d,f);
ConfCs:=proc(a,b,e,v0,c,d,f) local gu,eq,var,var1,var11,i1,j1:
gu:=GenConfC(a,b,e,v0,c,d,f);
var:={a,b,e,v0,c,d,f}:
eq:={seq(seq(gu[1][2][i1][j1]+gu[2][1][i1][j1]-1,i1=1..nops(gu[1][2])),j1=1..2),
seq(seq(gu[1][1][i1][j1]+gu[2][2][i1][j1]-1,i1=1..nops(gu[2][2])),j1=1..2)}:

var:=solve(eq,var):

if var=NULL then
 RETURN(FAIL):
fi:

var1:={}:

for var11 in var do
 if op(1,var11)=op(2,var11) then
   var1:=var1 union {op(1,var11)}:
 fi:
od:

subs(var,gu),var1:
end:





#GenConfCS(e,f): The general SYMMETRIC Jubin-like configurations. Try:
#GenConfCS(1/12,1/12);
GenConfCS:=proc(e,f):
[[[[0, 0], [1/2-e-2*f, 1/2-e-2*f], [1/2-e-2*f, 1/2], [-2*e-2*f+1/2, 1/2], [0, 2*e+2*f], [0, 0]], [[1/2, 0], [1/2+e+2*f, 0], [1/2+e+2*f, 1/2-e-2*f], [1/2+e+f, -e-f+1
/2], [1/2, -2*e-2*f+1/2], [1/2, 0]]], [[[1/2, 1], [1/2-e-2*f, 1], [1/2-e-2*f, 1/2+e+2*f], [-e-f+1/2, 1/2+e+f], [1/2, 1/2+2*e+2*f], [1/2, 1]], [[1, 1], [1/2+e+2*f, 1
/2+e+2*f], [1/2+e+2*f, 1/2], [1/2+2*e+2*f, 1/2], [1, 1-2*e-2*f], [1, 1]]]]:
end:



#GenConfCSP(e,f): Plots GenConfCS(e,f).
#Try:
#GenConfCSP(1/12,1/12);
GenConfCSP:=proc(e,f) local gu,pic:
gu:=GenConfCS(e,f):

pic:=plot({op(gu[1])},color=red,scaling=constrained):
pic:=pic,plot({op(gu[2])},color=black,scaling=constrained):
display(pic);
end:







#GenConfCSa(e,f,N): the white queen in the approximation to an N by N board of,
#followed by the black queens for GenConfC(e,f). Try:
#GenConfCSa(1/12,1/12,24);
GenConfCSa:=proc(e,f,N) local gu1,gu2,gu,mu,a,b,c,d,v0:
a:= 1/2-e-2*f:
b:= 2*e+2*f:
c:= e+2*f:
d:=-2*e-2*f+1/2:
v0:=1/2:
gu1:=LeftPentABE(trunc(N*a),trunc(N*b),trunc(N*e) ):
gu2:=RightPentABC([trunc(v0*N),0],trunc(c*N),trunc(d*N),trunc(f*N)):
gu:=gu1 union gu2:
mu:=Comp(ImpS(gu,N),N):
[gu,mu]:
end:







#GenConfCSaP(e,f,N): Plotting the white queen in the approximation to an N by N board. Try:
#GenConfCSaP(1/12,1/12,120);
GenConfCSaP:=proc(e,f,N) local gu:
gu:=GenConfCSa(e,f,N)[1]:
PlotWB(gu,N):
end:

#AreaSef(e,f): the area of GenConfCS(e,f);. Try:
#AreaSef(1/12,1/12);
AreaSef:=proc(e,f) local gu1,gu2,gu2b,gu2m,gu2t:
#gu1 is the area of the left Whitepentagon
gu1:=expand((e+f)*(1-2*e-4*f)-e^2/2):

#gu2b is the area of the bottom layer (a rectangle) of the right pentagon
gu2b:=(e+2*f)*(1/2-2*e-2*f):
#gu2m is the area of the middle layer of the right pentagon (a tapezoid)
gu2m:=(e+2*f)*e-e^2/2:
#gu2t is the area of the top layer of the right pentagon (a triangle)
gu2t:=f^2:
gu2:=gu2b+gu2m+gu2t:
expand(gu1+gu2):
end:



#AreaC(a,b,e,v0,c,d,f): the area of GenConfC(a,b,e,v0,c,d,f) . Try:
#AreaC(1/4,1/3,1/12,1/2,1/4,1/6,1/12);
AreaC:=proc(a,b,e,v0,c,d,f) local gu1,gu2,gu2b,gu2m, gu2t,guW,gu3,gu4,guB:

#GenConfC(a,b,e,v0,c,d,f) 

#gu1 is the area of the left Whitepentagon
gu1:=expand(a*b-e^2/2):

#gu2b is the area of the bottom layer (a rectangle) of the right pentagon
gu2b:=c*d:
#gu2m is the area of the middle layer of the right pentagon (a tapezoid)
gu2m:=(c+2*f)/2*(c-2*f):
#gu2t is the area of the top layer of the right pentagon (a triangle)
gu2t:=f^2:
gu2:=gu2b+gu2m+gu2t:
guW:=expand(gu1+gu2):



gu3:=(v0-a)*(1-v0-b)+(v0-a)^2/2-(-2*a-b+v0+2*c-2*f+d)*(1/2*v0+1/2*d-1/2*b+c-f-a)/2:


gu4:=(1-v0-c)*(v0-d)-(a+b-e-d-c)^2/2:
guB:=expand(gu3+gu4):

[guW,guB]:


end:





#MaxC(L,v): given  a list of length 2, L, consisting of polynomials in the list of variables v
#finds all the extreme points of L[1], subject to the constraint L[1]=L[2], using
#Largrange multipliers. Try
#MaxC([a^2+b^2,a^3+b^3],[a,b]);
#MaxC(AreaC(a,b,e,v0,c,d,f),[a,b,e,v0,c,d,f]);
MaxC:=proc(L,v) local g,gu,i1,eq,var,v1,var1:
gu:=L[1]+g*(L[2]-L[1]):
eq:={seq(diff(gu,v[i1]),i1=1..nops(v)), diff(gu,g)}:

var:=[solve(eq,{op(v),g})]:

var1:=[]:


for i1 from 1 to nops(var) do

  v1:=subs(var[i1],v):

  if  min(op(evalf(v1)))>0 and max(op(evalf(v1)))<1 then
    var1:=[op(var1),[v1,subs(var[i1],L[1])] ]:
  fi:
od:
var1:
end:



#BestGenConfC(): The best choice in GenConfC(a,b,e,v0,c,d,f) (two pentagons as in Jubin's construction).
#It returns the two White pentagons and the two Black pentagons that are the best
#Try:
#BestGenConfC();
BestGenConfC:=proc() local gu,mu,lu,a,b,e,v0,c,d,f,rec:

gu:=GenConfC(a,b,e,v0,c,d,f):

mu:=AreaC(a,b,e,v0,c,d,f):

lu:=MaxC(mu,[a,b,e,v0,c,d,f]):

if nops(lu)>1 then
  print(`There are more than one solution.`):
 RETURN(FAIL):
fi:

lu:=lu[1]:

rec:=lu[2]:

lu:=lu[1]:


gu:=subs({a=lu[1],b=lu[2],e=lu[3],v0=lu[4],c=lu[5],d=lu[6],f=lu[7]},gu):

gu,rec:

end:




#BestGenConfCv(): Verbose form of GenConfC(a,b,e,v0,c,d,f) (two pentagons as in Jubin's construction).
#It returns the two White pentagons and the two Black pentagons that are the best
#Try:
#BestGenConfCv();
BestGenConfCv:=proc() local gu,mu,lu,a,b,e,v0,c,d,f,rec:

gu:=GenConfC(a,b,e,v0,c,d,f):

print(`Consider all configurations where the White Queens are located in the following two pengagons`):
print(`inside the unit square [0,1]x[0,1] `):
print(``):
print(gu[1][1]):
print(``):
print(` and `):
print(``):
print(gu[1][2]):
print(``):
print(`then the Black Queens are located in the two pentagons`):
print(``):
print(gu[2][1]):
print(``):
print(` and `):
print(``):
print(gu[2][2]):

mu:=AreaC(a,b,e,v0,c,d,f):

print(`The area of the White Queens region is`):
print(``):
print(mu[1]):
print(``):
print(`and the area of the Black Queens region is`):
print(``):
print(mu[2]):
print(``):

lu:=MaxC(mu,[a,b,e,v0,c,d,f]):

if nops(lu)>1 then
  print(`There are more than one solution.`):
 RETURN(FAIL):
fi:

lu:=lu[1]:

rec:=lu[2]:

print(`Optimally, these should be the same. Maximizing them under this constaint shows that`):
print(`the maximal area (for each of them) is`):
print(``):
print(rec):
print(``):
lu:=lu[1]:


gu:=subs({a=lu[1],b=lu[2],e=lu[3],v0=lu[4],c=lu[5],d=lu[6],f=lu[7]},gu):

print(`and it is achieved with the following configuration for White followed by Black `):
print(``):
print(gu):

gu,rec:

end:


#TwoTriaConf(a,b): inputs numeric or symbolic a and b (such that 0<a<b<1)
#and output a list consistinbg of the two triangles with vertices
#(i)[[0,0],[a,0],[a,a]], [[b,0],[a+b,0],[a,+b,a]] (White area)
#(ii) The quadrilateral , triangle and parallelogram
#[[a,1],[a,a+b],[b,2*a],[b,1]], [[a+b,1],[a+b,a+b],[1,1]], [[a+b,b],]a+b,a],[1,1-b],[1,1-a]]
#(iii) The white area (a^2)
#(iv) The black area. try
#TwoTriaConf(a,b);
#TwoTriaConf(1/3,1/2);

TwoTriaConf:=proc(a,b) local gu:
gu:=factor((b-a)*(1-b-a)+(b-a)^2/2+(1-a-b)^2/2+(b-a)*(1-a-b)):
[{[[0,0],[a,0],[a,a]], [[b,0],[a+b,0],[a,+b,a]] },
{[[a,1],[a,a+b],[b,2*a],[b,1]], [[a+b,1],[a+b,a+b],[1,1]], [[a+b,b],[a+b,a],[1,1-b],[1,1-a]]},
a^2,gu]:
end:


#GenRect(v0,A,B): The general A by B rectangle with vertex v0. Try:
#GenRect([0,0],10,10);
GenRect:=proc(v0,A,B) local i,j:
{seq(seq(v0+[i,j],j=0..B-1),i=0..A-1)}:
end:

#GenLeftTri(v0,A): The general equilateral triangle with vertex v0 and side of size A where the hypot. is pointing left.
#Try
#GenLeftTri([0,0],30);
GenLeftTri:=proc(v0,A) local i,j:
{seq(seq(v0+[i,j],j=0..i),i=0..A-1)}:
end:



#GenRightTri(v0,A): The general equilateral triangle with vertex v0 and side of size A where the hypot. is pointing right.
#Try
#GenRightTri([0,0],30);
GenRightTri:=proc(v0,A) local i,j:
{seq(seq(v0+[i,j],j=i..A-1),i=0..A-1)}:
end:



#BestGenConfCvLP(): Verbose form of GenConfC(a,b,e,v0,c,d,f) (two pentagons as in Jubin's construction).
#It returns the two White pentagons and the two Black pentagons that are the best
#Try:
#BestGenConfCvLP();
BestGenConfCvLP:=proc() local gu,mu,lu,a,b,e,v0,c,d,f,rec:

gu:=GenConfC(a,b,e,v0,c,d,f):

print(`Consider all configurations where the White Queens are located in the following two pengagons`):
print(`inside the unit square [0,1]x[0,1] `):
print(``):
lprint(gu[1][1]):
print(``):
lprint(gu[1][1]):
print(` and `):
print(``):
lprint(gu[1][2]):
print(``):
print(`then the Black Queens are located in the two pentagons`):
print(``):
lprint(gu[2][1]):
print(``):
print(` and `):
print(``):
lprint(gu[2][2]):

mu:=AreaC(a,b,e,v0,c,d,f):

print(`The area of the White Queens region is`):
print(``):
lprint(mu[1]):
print(``):
print(`and the area of the Black Queens region is`):
print(``):
lprint(mu[2]):
print(``):

lu:=MaxC(mu,[a,b,e,v0,c,d,f]):

if nops(lu)>1 then
  print(`There are more than one solution.`):
 RETURN(FAIL):
fi:

lu:=lu[1]:

rec:=lu[2]:

print(`Optimally, these should be the same. Maximizing them under this constaint shows that`):
print(`the maximal area (for each of them) is`):
print(``):
lprint(rec):
print(``):
lu:=lu[1]:


gu:=subs({a=lu[1],b=lu[2],e=lu[3],v0=lu[4],c=lu[5],d=lu[6],f=lu[7]},gu):

print(`and it is achieved with the following configuration for White followed by Black `):
print(``):
lprint(gu):

gu,rec:

end:



#########Start BestLeftPent
#BestLeftPentC(e1): The best choice in LeftPentC(a,b,e) with e=e1. (the left pentagon  in the Jubin's construction).
#It returns the White pentagon and left pentagon
#Try:
#BestLeftPentC(0.292);
BestLeftPentC:=proc(e1) local a,b,e,gu,mu,lu:

mu:=[a*b-e^2/2, (1-a-b)^2/2+(1-a-b+e)^2/2]:

lu:=MaxC(subs(e=e1,mu),[a,b])[1]:
gu:=LeftPentCabe(a,b,e):

subs({a=lu[1][1],b=lu[1][2],e=e1},gu),lu:

end:
#########End BestLeftPent






#Hexagon(a,b,c,d): the hexagon in lower left corner.
Hexagon:=proc(a,b,c,d)
[[0,0], [a,0], [a+b,b], [a+b,b+c],[b+c-d,b+c],[0,d]]:
end:

#BestHexagon(): the hexagon configuration in the lower left corner which maximizes the number of queens as a local maximum. But actually it’s not a global optimum.
BestHexagon:=proc() local AreaB, AreaW,a,b,c,d:
AreaB:=1/2*((1-a-b-c)^2+(1-a-b-d)^2):
AreaW:=(a+b)*(b+c)-1/2*b^2-1/2*(b+c-d)^2:
MaxC1([a^2+2*a*b+a*c+b^2+b*c+(1/2)*c^2-2*a-2*b-c+1+a*d+b*d+(1/2)*d^2-d, a*b+a*c+b*d-(1/2)*c^2+c*d-(1/2)*d^2], [a,b,c,d]):
end:


#Rectangle(s,a,b): a rectangle whose lower left point is [s,0] with width a and height b.
Rectangle:=proc(s,a, b) local i,j:
{seq(seq([s+j,i],j=0..a-1),i=0..b-1)}:
end:


#Triangle(s,a): a triangle whose lower left point is [s,0] and upper right point is [s+a-1, a-1] with width a and height a.
Triangle:=proc(s,a) local i,j:
{seq(seq([s+j,i],j=i..a-1),i=0..a-1)}:
end:


#Triangle2(s,a): a triangle whose lower left point is [s,0], lower right point is [s+a-1,0] and upper left point is [s, a-1] with width a and height a.
Triangle2:=proc(s,a) local i,j:
{seq(seq([s+j,i],j=0..a-1-i),i=0..a-1)}:
end:




#Plot2Square(a,s,N): plot two squares of edge length a, the left one’s lower left corner is [0,0], the right one’s is [s,0] on N by N chess board.
Plot2Square:=proc(a,s,N)
PlotWB(Rectangle(0,a,a) union Rectangle(s,a,a), N):
end:


#NuA2Square(a,s,N): the non-attacking area in the two identical squares configuration. Meaning of parameters are the same withPlot2Square.
NuA2Square:=proc(a,s,N) local co:
co:=NuA(Rectangle(0,a,a) union Rectangle(s,a,a), N):
co, evalf(co/N^2), evalf(2*a^2/N^2):
end:


#FindM2Square(s,N): for different s which is the x-coordinate of the right square’s lower left corner, find out the interval [a, a+1] and [a/N, (a+1)/N] such that at a and a+1 of the two ratios white/black one is >1 and the other is <1.
FindM2Square:=proc(s, N) local i:
for i from floor(N/5) to min(s, N-s) do
  if NuA2Square(i, s, N)[3] / NuA2Square(i,s,N)[2] >=1 then
    return [i-1, i], [evalf((i-1)/N), evalf(i/N)]:
  fi:
od:
return fail:
end:


#Plot2Triangle(a,s,N): plot two triangles of edge length a (bottom and right), the left one’s lower left corner is [0,0], the right one’s is [s,0] on N by N chess board.
Plot2Triangle:=proc(a,s,N)
PlotWB(Triangle(0,a) union Triangle(s,a), N):
end:

#NuA2Triangle(a,s,N): the non-attacking area in the two identical triangles configuration. Meaning of parameters are the same withPlot2Triangle.
NuA2Triangle:=proc(a,s,N) local co:
co:=NuA(Triangle(0,a,a) union Triangle(s,a,a), N):
co, evalf(co/N^2), evalf(a^2/N^2):
end:

#FindM2Triangle(s,N): for different s which is the x-coordinate of the right square’s lower left corner, find out the interval [a, a+1] and [a/N, (a+1)/N] such that at a and a+1 of the two ratios white/black one is >1 and the other is <1.
FindM2Triangle:=proc(s, N) local i:
for i from floor(N/5) to min(s, N-s) do
  if NuA2Triangle(i, s, N)[3] / NuA2Triangle(i,s,N)[2] >=1 then
    return [i-1, i], [evalf((i-1)/N), evalf(i/N)]:
  fi:
od:
return fail:
end:

#PlotSquareTriangle(a,s,N): plot one square and one triangle of edge length a, the left square’s lower left corner is [0,0], the right triangle’s is [s,0] on N by N chess board.
PlotSquareTriangle:=proc(a,s,N)
PlotWB(Rectangle(0,a,a) union Triangle(s,a,a), N):
end:


#NuASquareTriangle(a,s,N): the non-attacking area in the square + triangle (the same side length) configuration. Meaning of parameters are the same with PlotSquareTriangle.
NuASquareTriangle:=proc(a,s,N) local co:
co:=NuA(Rectangle(0,a,a) union Triangle(s,a,a), N):
co, evalf(co/N^2), evalf(3/2*a^2/N^2):
end:


#FindMSquareTriangle(s,N): for different s which is the x-coordinate of the right triangle’s lower left corner, find out the interval [a, a+1] and [a/N, (a+1)/N] such that at a and a+1 of the two ratios white/black one is >1 and the other is <1.
FindMSquareTriangle:=proc(s, N) local i:
for i from floor(N/5) to min(s, N-s) do
  if NuASquareTriangle(i, s, N)[3] / NuASquareTriangle(i,s,N)[2] >=1 then
    return [i-1, i], [evalf((i-1)/N), evalf(i/N)]:
  fi:
od:
return fail:
end:

#MaxC1(L,v): given  a list of length 2, L, consisting of polynomials in the list of variables v
#finds all the extreme points of L[1], subject to the constraint L[1]=L[2], using
#Largrange multipliers. Compared to MaxC, this procedures doesn’t do pre-check. It will return all possible solutions including complex roots. Try
#MaxC1([a^2+b^2,a^3+b^3],[a,b]);
#MaxC1(AreaC(a,b,e,v0,c,d,f),[a,b,e,v0,c,d,f]);
MaxC1:=proc(L,v) local g,gu,i1,eq,var,v1,var1:
gu:=L[1]+g*(L[2]-L[1]):
eq:={seq(diff(gu,v[i1]),i1=1..nops(v)), diff(gu,g)}:

var:=[solve(eq,{op(v),g})]:

var1:=[]:


for i1 from 1 to nops(var) do

  v1:=subs(var[i1],v):

  
  var1:=[op(var1),[v1,subs(var[i1],L[1])] ]:
 
od:
var1:
end:


#Plot2Squares(a,b,s,N): plot two squares of edge length a and b, the left one’s lower left corner is [0,0], the right one’s is [s,0] on N by N chess board.
Plot2Squares:=proc(a,b,s,N)
PlotWB(Rectangle(0,a,a) union Rectangle(s,b,b), N):
end:


#NuA3Square(a,N): the non-attacking area in the three equidistance identical squares configuration. The side length is a. 
NuA3Square:=proc(a,N) local co:
co:=NuA(Rectangle(0,a,a) union Rectangle(floor(N/3),a,a) union Rectangle(floor(2*N/3),a,a), N):
co, evalf(co/N^2), evalf((3*a^2)/N^2):
end: