Quiet[Remove["Global`*"],{Remove::rmnsm}];Print["Mathematica $Version = “",$Version,"”"];Print["Execution time = ",DateString[DateList[],{"Hour",":","Minute"," on ","DayNameShort"," ","Day"," ","MonthNameShort"," ","Year"}]];

(*
	http://www.jdawiseman.com/papers/placemat/placemat.html#PermittedPackingStyles 
	http://www.jdawiseman.com/papers/placemat/PackingStyles_formulae.nb
	http://www.jdawiseman.com/papers/placemat/PackingStyles_formulae.pdf
*)

(*
	Two stage process. 

	Stage one, calculate maximum possible r. Generally this is done by having adjacent circles touch, by solving polynomial(s) of the form something == 4 (r^2). 
	But other constraints about non-adjacent circles not overlapping are also needed. E.g., in PostsAndLintel, also r ≤ w/4, otherwise in all narrow pages the first and last circles might overlap even though ‘adjacent’ circles only touch.
	Typically this stage entails creating and solving a polynomial in r, with coefficients involving the likes of h, w, nR, nC, etc. 

	Stage 2, r is further lowered by the likes of MaxRadius (a hard maximum) or ShrinkRadii (pages with different layouts to have same radius). 
	Then, using this smaller r, circles are positioned. 
	The broad principle is that the gap all pairs of ‘adjacent’ circles are the same distance apart.
*)

(* r = radius; h = usable Height; w = usable Width; nC= numCols; nSC = numShortColumns; nR= numRows; nG = numGlasses (or numGlassesAtEdge); gC = gap between Cols; gR = gap between Rows; gRD = gap between rows at dislocation; gCD = gap between columns in Diamond pattern *)
$Assumptions=Apply[And,Map[#∈Reals&&#>0&,{w,h,r,nR,nC,nSC,nG,gC,gR,gRD,gCD,x,y,x0,y0,x1,y2,x2,y2,x3,y3,x4,y4,x5,y5,x6,y6}]];

With[{mgn=24.,header=0.},  (* For analytic solutions set mgn and header to exact value *)
	(({{PortraitA4}, {LandscapeA4}, {LandscapeA3}, {LandscapeLegal}, {LandscapeLedger}, {LandscapePoint6}}))=(({{w->210 360/127-2mgn, h->297 360/127-2mgn-header}, {w->297 360/127-2mgn, h->210 360/127-2mgn-header}, {w->420 360/127-2mgn, h->297 360/127-2mgn-header}, {w->14 72-2mgn, h->(8+1/2) 72-2mgn-header}, {w->17 72-2mgn, h->11 72-2mgn-header}, {w->1000., s->600.}}));
];

(* PostScriptForm[] *)
(*
	http://mathematica.stackexchange.com/questions/101954/postscriptform-or-forthform
	http://mathematica.stackexchange.com/questions/102894/multi-case-function-many-single-case-delayed-assignments-or-one-which
*)
Remove[PostScriptForm];
PostScriptForm[thing_Rational]:=ToString[N[thing,20],InputForm,NumberMarks->False];
PostScriptForm[thing_?AtomQ]:=ToString[thing];
PostScriptForm[thing_List]:=StringJoin@@Riffle[Map[PostScriptForm,thing],"\r\n"];
PostScriptForm[MatrixForm[thing_]]:=PostScriptForm[thing];
PostScriptForm[Times[-1,thing_]]:=StringJoin[PostScriptForm[thing]," neg"];

PostScriptForm[thing_Power]:=(
	psExponent:=Which[
		#>5&&Divisible[#,3],psExponent[#/3]<>" dup dup mul mul",
		#≥5&&OddQ[#],"dup "<>psExponent[(#-1)/2]<>" dup mul mul" ,
		#≥4&&EvenQ[#],psExponent[#/2]<>" dup mul",
		#==3,"dup dup mul mul",
		#==2,"dup mul",
		#==1/2,"sqrt",
		(Rational===Head[#])&&(Log[2,#//Denominator]//IntegerQ),psExponent[Simplify[2#]]<>" sqrt",
		Not[IntegerQ[#]],PostScriptForm[#]<>" exp",
		#==1,"",
		True," !•Error•! "
	]&;
	Which[
		thing[[2]]>0||Not[IntegerQ[thing[[2]]]],PostScriptForm[thing[[1]]]<>" "<>psExponent[thing[[2]]],
		thing[[2]]==-1,"1 "<>PostScriptForm[thing[[1]]]<>" div",
		thing[[2]]==0,"1",
		True,"1 "<>PostScriptForm[thing[[1]]]<>" "<>psExponent[-thing[[2]]]<>" div"
	]);

PostScriptForm[thing_Times]:=StringJoin@Riffle[Reap[If[MatchQ[thing[[1]],Power[_,n_Integer/;n<0]],(Sow["1 "<>PostScriptForm[thing[[1,1]]]<>" div"];),(Sow[PostScriptForm[thing[[1]]]];)];Map[(If[MatchQ[#,Power[_,n_Integer/;n<0]],(Sow[PostScriptForm[#[[1]]^(-#[[2]])]<>" div"];),(Sow[PostScriptForm[#]<>" mul"]; )])&,Drop[List@@thing,1]]][[2,1]]," "];

PostScriptForm[thing_Plus]:=
StringJoin@@If[FreeQ[thing,_^n_],
	(* Simple expression, no powers, to be summed one item at a time *)
	Module[{i},
		i=Position[thing,Except[Times[-1,_]|(_?Negative)],1,Heads->False];
		If[Length[i]>0,i=i[[1,1]],(i=Position[thing,Not[MatchQ[#,Times[-1,_]]]&,1,Heads->False];i=If[Length[i]>0,i[[1,1]],1])];	Prepend[Map[(" "<>Replace[#,{(n_Integer/;n<0:>ToString[-n]<>" sub"),(Times[-1,_]:>PostScriptForm[Times@@Drop[#,1]]<>" sub"),(Times[n_/;n<0,_]:>PostScriptForm[Times@@Drop[#,1]]<>" "<>ToString[-#[[1]]]<>" mul sub"),(Times[n_/;n>0,_]:>PostScriptForm[Times@@Drop[#,1]]<>" "<>ToString[#[[1]]]<>" mul add"),(_:>PostScriptForm[#]<>" add")}])&,Drop[List@@thing,{i}]],Replace[thing[[i]],{Times[-1,_]:>PostScriptForm[-thing[[i]]]<>" neg",_:>PostScriptForm[thing[[i]]]}]]  ],
	(* Polynomial *)
	Module[{vars,exps,v,rcl,i,firstMul},
		vars=Variables[thing];
		exps=Exponent[thing,vars];
		v=Select[Transpose[{vars,exps}],(#[[2]]==Max@@exps)&][[1,1]];
		rcl=Reverse[Map[Factor,CoefficientList[thing,v]]];
		Reap[
			i=1;firstMul=True;If[rcl[[1]]=!=1,Sow[PostScriptForm[rcl[[1]]]]];
Map[If[#===0,i++,(Sow[If[firstMul&&rcl[[1]]===1,PostScriptForm[v^i]<>" "," "<>PostScriptForm[v^i]<>" mul "]<>If[MatchQ[#,(Times[_?Negative,_]|(_?Negative))],PostScriptForm[-#]<>" sub",PostScriptForm[#]<>" add"]];i=1;firstMul=False)]&,Drop[rcl,1]];
			If[i>1,Sow[" "<>PostScriptForm[v^(i-1)]<>" mul "]];
		][[2,1]]
]];

Diamonds

Print[];

(* http://www.jdawiseman.com/papers/placemat/placemat.html#PermittedPackingStyles_Diamonds *)

(* Also r bounded above by h/(nR+1) and by w/(nC+1) *)

Module[{poly},
	poly=((((((w-2r)/(nC-1))))^2)+(((((h-2r)/(nR-1))))^2)-4 (r^2))//Together//Numerator//Factor;
	Print[poly];
	Print[Map[Factor,CoefficientList[If[SameQ[Head[poly],Times],poly[[-1]],poly],r]]//MatrixForm];
	Print[Min[Select[r/.Solve[0==poly/.LandscapeLegal/.{nR->3,nC->6},{r}],#∈Reals&&#>0&]]];
];

With[{params={gC,gR}},Do[(
	Module[{poly},
	poly=Factor[GroebnerBasis[{
		h==2r+(nR-1)gR,
		w==2r+(nC-1)gC},params[[i]],Drop[params,{i}]][[1]]];
	Print[params[[i]],":\t",Map[Factor,CoefficientList[If[SameQ[Head[poly],Times],poly[[-1]],poly],params[[i]]]]//MatrixForm];
	Print[Solve[0==poly/.LandscapeLegal/.{nR->3,nC->6,r->110},params[[i]]]]
];),{i,Length[params]}]]

RectangularDislocation

Print[];

(* http://www.jdawiseman.com/papers/placemat/placemat.html#PermittedPackingStyles_RectangularDislocation *)

Module[{poly},
	poly=Factor[GroebnerBasis[{
		h==(2nR-2)r+gR,
		w==2r+2(nC-1)gC,
		(gR^2)+(gC^2)==4 (r^2)},{r},{gR,gC}][[1]]];
	Print[poly];
	Print[Map[Factor,CoefficientList[If[SameQ[Head[poly],Times],poly[[-1]],poly],r]]//MatrixForm];
	Print[Min[Select[r/.Solve[0==poly/.LandscapeA3/.{nR->3,nC->4},{r}],#∈Reals&&#>0&]]];
];

(* But r might be shrunk by ShrinkRadii or by MaxRadius. So need to solve separate equations for the positions (gC,gR,gRD) which will depend on r. *)

With[{params={gC,gR,gRD}},Do[(
	Module[{poly},
	poly=Factor[GroebnerBasis[{
		h==2r+(nR-2)gR+gRD,
		w==2r+2(nC-1)gC,
		(gRD^2)+(gC^2)==(gR^2)},params[[i]],Drop[params,{i}]][[1]]];
	Print[params[[i]],":\t",Map[Factor,CoefficientList[If[SameQ[Head[poly],Times],poly[[-1]],poly],params[[i]]]]//MatrixForm];
	Print[Solve[0==poly/.LandscapeA3/.{nR->3,nC->4,r->135},params[[i]]]]
];),{i,Length[params]}]]

DiamondsAndRectangular

Print[];

(* http://www.jdawiseman.com/papers/placemat/placemat.html#PermittedPackingStyles_DiamondsAndRectangular *)

(* Also r bounded above by h/(nR+1) and by w/(nC-nSC)/2 *)

(* It is obvious that this cannot be quite the maximum r. Imagine moving C2 and H7 vertically towards each other, and nudging M12 up by the same distance. Then there would be gaps either side of that column, which could be filled by increasing r. But only by a mite. And that equation is not solved here. *)

Module[{poly},
	poly=Factor[GroebnerBasis[{
		(* Equations need doing *)
		h==2r+gR(nR-1),
		w==2gC nSC+2r(nC-2nSC),
		(gR^2)+(gC^2)==4 (r^2)},{r},{gR,gC}][[1]]];
	Print[poly];
	Print[Map[Factor,CoefficientList[If[SameQ[Head[poly],Times],poly[[-1]],poly],r]]//MatrixForm];
	Print[Min[Select[r/.Solve[0==poly/.LandscapeA3/.{nR->5,nC->5,nSC->1},{r}],#∈Reals&&#>0&]]];
];

With[{params={gR,gC,gCD}},Do[(
	Module[{poly},
	poly=Factor[GroebnerBasis[{
		h==2r+(nR-1)gR,
		w==2r+2nSC gCD+(nC-1-2nSC)gC,
		(gCD^2)+(gR^2)==(gC^2)},params[[i]],Drop[params,{i}]][[1]]];
	Print[params[[i]],":\t",Map[Factor,CoefficientList[If[SameQ[Head[poly],Times],poly[[-1]],poly],params[[i]]]]//MatrixForm];
	Print[Solve[0==poly/.LandscapeA3/.{nR->5,nC->5,nSC->1,r->120},params[[i]]]]
];),{i,Length[params]}]]

RectangularAlternateNudge

Print[];

(* http://www.jdawiseman.com/papers/placemat/placemat.html#PermittedPackingStyles_RectangularAlternateNudge *)

(* Also r bounded above by h/(2+(Sqrt[3])(nR-1)) and by w/(2 nC) *)

Module[{poly},
	poly=Factor[GroebnerBasis[{
		h==2r+gR(nR-1),
		w==2nC r+gC,
		(gR^2)+(gC^2)==4 (r^2)},{r},{gR,gC}][[1]]];
	Print[poly];
	Print[Map[Factor,CoefficientList[If[SameQ[Head[poly],Times],poly[[-1]],poly],r]]//MatrixForm];
	Print[Min[Select[r/.Solve[0==poly/.LandscapeA3/.{nR->3,nC->4},{r}],#∈Reals&&#>0&]]];
];

With[{params={gR,gC}},Do[(
	Module[{poly},
	poly=Factor[GroebnerBasis[{
		h==2r+gR(nR-1),
		w==2nC r+gC},params[[i]],Drop[params,{i}]][[1]]];
	Print[params[[i]],":\t",Map[Factor,CoefficientList[If[SameQ[Head[poly],Times],poly[[-1]],poly],params[[i]]]]//MatrixForm];
	Print[Solve[0==poly/.LandscapeA3/.{nR->3,nC->4,r->130},params[[i]]]]
];),{i,Length[params]}]]

RectangularAlternateSplitNudge

Print[];

(* http://www.jdawiseman.com/papers/placemat/placemat.html#PermittedPackingStyles_RectangularAlternateSplitNudge *)

(* Also r bounded above by h/(2+(Sqrt[3])(nR-1)) and by w/(2 nC) *)

Module[{poly},
	poly=Factor[GroebnerBasis[{
		h==2r+gR(nR-1),
		w==2nC r+2gC,
		(gR^2)+(gC^2)==4 (r^2)},{r},{gR,gC}][[1]]];
	Print[poly];
	Print[Map[Factor,CoefficientList[If[SameQ[Head[poly],Times],poly[[-1]],poly],r]]//MatrixForm];
	Print[Min[Select[r/.Solve[0==poly/.LandscapeLedger/.{nR->3,nC->4},{r}],#∈Reals&&#>0&]]];
];

With[{params={gR,gC}},Do[(
	Module[{poly},
	poly=Factor[GroebnerBasis[{
		h==2r+gR(nR-1),
		w==2nC r+2gC},params[[i]],Drop[params,{i}]][[1]]];
	Print[params[[i]],":\t",Map[Factor,CoefficientList[If[SameQ[Head[poly],Times],poly[[-1]],poly],params[[i]]]]//MatrixForm];
	Print[Solve[0==poly/.LandscapeA3/.{nR->3,nC->4,r->120},params[[i]]]]
];),{i,Length[params]}]]

DiamondsPlus

Print[];

(* www.jdawiseman.com/papers/placemat/placemat.html#PermittedPackingStyles_DiamondsPlus *)

(* Also r bounded above by h/(nR+(Sqrt[3])) and by w/4 *)

Module[{poly},
	poly=Factor[GroebnerBasis[{
		w==2r+2gC,
		h==2r+(nR-2)gR+gRD,
		(gR^2)+(gC^2)==4 (r^2),
		(gRD^2)+(((w/2-2r))^2)==4 (r^2)},{r},{gR,gC,gRD}][[1]]];
	Print[poly];
	Print[Map[Factor,CoefficientList[If[SameQ[Head[poly],Times],poly[[-1]],poly],r]]//MatrixForm];
	Print[Min[Select[r/.Solve[0==poly/.PortraitA4/.{nR->4},{r}],#∈Reals&&#>0&]]];
];

With[{params={gR,gC,gRD,gCD}},Do[(
	Module[{poly},
	poly=Factor[GroebnerBasis[{
		w==2r+2gC,
		h==2r+(nR-2)gR+gRD,
		(gR^2)+(gC^2)==(gRD^2)+(gCD^2),
		(gR^2)+(gC^2)==(0^2)+(((w-2r-2gCD))^2)},params[[i]],Drop[params,{i}]][[1]]];
	Print[params[[i]],":\t",Map[Factor,CoefficientList[If[SameQ[Head[poly],Times],poly[[-1]],poly],params[[i]]]]//MatrixForm];
	Print[Solve[0==poly/.PortraitA4/.{nR->4,r->110},params[[i]],Reals]];
];),{i,Length[params]}]]

PostsAndLintel

(* Also r bounded above by h/(nR-1)/2 and by w/4, and if one central glass by w/6, and if two or three by w/8. Central glass(es) also cause bounds in h that vary according to whether odd or even on edge. *)

Print[];

(* http:///www.jdawiseman.com/papers/placemat/placemat.html#PermittedPackingStyles_Arch_PostsAndLintel *)

With[{
	x2=r,
	y2=(2 nR-3)r ,
	x3=w/2-(nG+1-2nR)r,
	y3=h-r},
	Module[{poly},
		poly=Factor[GroebnerBasis[{
			(((x2-x3))^2)+(((y2-y3))^2)==4 (r^2)},{r},{}][[1]]];
		Print[poly];
		Print[Map[Factor,CoefficientList[If[SameQ[Head[poly],Times],poly[[-1]],poly],r]]//MatrixForm];
		Print[Min[Select[r/.Solve[0==poly/.LandscapeA3/.{nR->3,nG->7},{r}],#∈Reals&&#>0&]]];
]];

With[{params={gR,gC,gRD,gCD}},Do[(
	Module[{poly},
	poly=Factor[GroebnerBasis[{
		h==2r+(nR-2)gR+gRD,
		w==2r+2gCD+(nG+1-2nR)gC,
		gC==gR,
		(gR^2)==(gRD^2)+(gCD^2)},params[[i]],Drop[params,{i}]][[1]]];
	Print[params[[i]],":\t",Map[Factor,CoefficientList[If[SameQ[Head[poly],Times],poly[[-1]],poly],params[[i]]]]//MatrixForm];
	Print[Min[Select[params[[i]]/.Solve[0==poly/.LandscapeA4/.{nR->3,nG->7,r->130},params[[i]]],#∈Reals&&#>0&]]];
];),{i,Length[params]}]]

Bespoke5

Print[];

(* http://www.jdawiseman.com/papers/placemat/placemat.html#PermittedPackingStyles_Bespoke *)

(* Also r bounded above by h/(2+(Sqrt[3])) and by w/(2+2(Sqrt[3])) *)

With[{
	x0=r,
	y0=h-r,
	x2=w/2,
	y3=r},
	Module[{poly},
		poly=Factor[GroebnerBasis[{
			(((x0-x2))^2)+(((y0-y2))^2)==4 (r^2),
			(((x2-x3))^2)+(((y2-y3))^2)==4 (r^2),
			(((x3-x0))^2)+(((y3-y0))^2)==4 (r^2)},{r},{y2,x3}][[1]]];
		Print[poly];
		Print[Map[Factor,CoefficientList[If[SameQ[Head[poly],Times],poly[[-1]],poly],r]]//MatrixForm];
		Print[Min[Select[r/.Solve[0==poly/.LandscapeA4,{r}],#∈Reals&&#>0&]]];
]];

Bespoke7

Print[];

(* http://www.jdawiseman.com/papers/placemat/placemat.html#PermittedPackingStyles_Bespoke *)

(* Also r bounded above by h/(2+(Sqrt[3])) and by w/6 *)

With[{
	y0=h-r,
	x2=w/2,
	x3=r,
	x5=w/2-r,
	y5=r},
	Module[{poly},
		poly=Factor[GroebnerBasis[{
			(((x0-x2))^2)+(((y0-y2))^2)==4 (r^2),
			(((x2-x5))^2)+(((y2-y5))^2)==4 (r^2),
			(((x5-x3))^2)+(((y5-y3))^2)==4 (r^2),
			(((x3-x0))^2)+(((y3-y0))^2)==4 (r^2)},{r},{x0,y2,y3}][[1]]];
		Print[poly];
		Print[Map[Factor,CoefficientList[If[SameQ[Head[poly],Times],poly[[-1]],poly],r]]//MatrixForm];
		Print[Min[Select[r/.Solve[0==poly/.LandscapeLegal,{r}],#∈Reals&&#>0&]]];
]]

Temple

Print[];

(* http://www.jdawiseman.com/papers/placemat/placemat.html#PermittedPackingStyles_Temple *)

(* Octic has a solution in a narrow range of aspect ratios. Other aspect ratios need spaced-out variant, coping with which is fiddly. *)

Module[{boundariesDist,polyRadiusEqns,polyR,specialR},
	With[{dist=2r },With[{
		row0X=(l-dist(n-1))/2,
		row0Y=s-r,
		row1X=r,
		(* row1Y *)
		row2X=l/2,
		(* row2Y *)
		(* row3X *)
		row3Y=r,
		centreOffsetX=dist(n-2)/2
	},
		boundariesDist={r>0,row1Y>0,row1Y<s,row2Y>0,row2Y<s,row3X>0};
		polyRadiusEqns={
			(((row1Y-row0Y))^2)+(((row1X-row0X))^2)==(dist^2),
			(((row3Y-row1Y))^2)+(((row3X-row1X))^2)==(dist^2),
			(((row1Y-row2Y))^2)+(((row1X+centreOffsetX-row2X))^2)==(dist^2),
			(((row3Y-row2Y))^2)+(((row3X+centreOffsetX-row2X))^2)==(dist^2)
		};
	]];
	Print[polyRadiusEqns//MatrixForm];
	polyR=Factor[GroebnerBasis[polyRadiusEqns,{r},{row1Y, row2Y,row3X},MonomialOrder->EliminationOrder][[1]]];
	Print[polyR];
	Print[Map[Factor,polyRcfl=CoefficientList[If[Head[polyR]===Times,polyR[[-1]],polyR],r]]//MatrixForm];
	Print[N[polyRcfl/.LandscapePoint6/.{n->4},12]//MatrixForm];
	Print["0.6×1: ",N[specialR=Min[Select[r/.Solve[0==polyR/.(LandscapePoint6/.{w->l,h->s})/.{n->4},{r}],#∈Reals&&#>0&]],12]];
Print[specialLedgerR=Min[Select[r/.Solve[0==polyR/.{n->4}/.LandscapeLedger/.{w->l,h->s},{r}],#∈Reals&&#>0&]]];
StringJoin@@Riffle[StringJoin/@Transpose[{Map[StringTrim[#//Simplify//Factor//PostScriptForm]&,polyRcfl/.{s->S,l->L,n->N}],Table["  % Coeff "<>ToString[i],{i,0,Length[polyRcfl]-1}]}],"\r\n"]
];

Module[{boundariesDist,polyDistEqns,specialTestValues,polyDist,cflDist},
	With[{
		row0X=(l-dist(n-1))/2,
		row0Y=s-r,
		row1X=r,
		(* row1Y *)
		row2X=l/2,
		(* row2Y *)
		(* row3X *)
		row3Y=r,
		centreOffsetX=dist(n-2)/2
	},
		boundariesDist={dist>0,row1Y>0,row1Y<s,row2Y>0,row2Y<s,row3X>0};
		polyDistEqns={
			(((row1Y-row0Y))^2)+(((row1X-row0X))^2)==(dist^2),
			(((row3Y-row1Y))^2)+(((row3X-row1X))^2)==(dist^2),
			(((row1Y-row2Y))^2)+(((row1X+centreOffsetX-row2X))^2)==(dist^2),
			(((row3Y-row2Y))^2)+(((row3X+centreOffsetX-row2X))^2)==(dist^2)
		};
	];
	specialTestValues=Flatten[{(LandscapePoint6/.{w->l,h->s}),n->4,r->103.38797850731751591628004827224490933}];
	Print[Solve[Join[boundariesDist,polyDistEqns]/.specialTestValues,{dist,row1Y, row2Y,row3X},Reals]];
	polyDist=Factor[GroebnerBasis[polyDistEqns,{dist},{row1Y, row2Y,row3X},MonomialOrder->EliminationOrder][[1]]];
	Print[Solve[0==polyDist/.specialTestValues,dist,Reals]];
	Print[polyDist];
	cflDist=Factor/@Simplify/@CoefficientList[If[SameQ[Head[polyDist],Times],polyDist[[-1]],polyDist],dist];
	Print[cflDist//MatrixForm];
	Print[cflDist/.specialTestValues//MatrixForm];
	Print[Min[Select[(dist/r/.specialTestValues)/.Solve[0==polyDist/.specialTestValues,{dist},Reals],#∈Reals&&#>0&]]];
	StringJoin@@Riffle[StringJoin/@Transpose[{Map[StringTrim[#//PostScriptForm]&,(cflDist/.{s->S,l->L,r->R,n->N})],Table["  % Coeff "<>ToString[i],{i,0,Length[cflDist]-1}]}],"\r\n"]
]

Arch

Print[];

(* http://www.jdawiseman.com/papers/placemat/placemat.html#PermittedPackingStyles_Arch_PostsAndLintel *)

(* No analytic solution. ☹ *)

(* But r bounded above by h/2 and by w/4, and if one central glass by w/6, and if two or three by w/8. Central glass(es) also cause bounds in h that vary according to whether odd or even on edge. *)

(* http://www.jdawiseman.com/2015/20151128_GroebnerBasis_slow.nb *)