The Gerrymandering Problem

OEIS A348456: Number of ways to dissect a 2*n X 2*n chessboard into two polyominoes each of area 2*n^2.

Code

SetDirectory[NotebookDirectory[]];
Quiet[<<Combinatorica`];

(* Determine chunks (runs) of same elements in a {0,1}-vector. *)
Chunks[vec_]:=Map[Last,SplitBy[Transpose[{vec,Range[Length[vec]]}],First],{2}];
(* Construct all non-crossing arc configurations on n nodes with alternating colors. *)
NoncrossingConnections[n_Integer]:=NoncrossingConnections[n]=
Select[Tuples[SetPartitions[Table[i,{i,#,n,2}]]&/@{1,2}],Cases[Subsets[Join@@#,{2}],{{___,a1_,___,a2_,___},{___,b1_,___,b2_,___}}/;(a1<b1<a2<b2||b1<a1<b2<a2)]==={}&];
NoncrossingConnections[chunks_List]:=Map[Flatten,NoncrossingConnections[Length[chunks]]/.a_Integer:>chunks[[a]],{3}];

(* Rows and columns of the transfer matrix are labeled by states of the form {col, co0, co1}. *)
(* Example: {col, co0, co1} = {{1,0,0,1,0,1,1,0}, {{2,3},{5,8}}, {{1},{4},{6,7}}} means 
that the 0's at positions 5 and 8 are connected, while the three chunks of 1's are not. *)
(* Step (1): Construct the set of all possible (reachable) states. *)
MatrixLabels1[n_]:=Join[Flatten[
Function[v,If[v[[1]]===0,{v,#1,#2},{v,#2,#1}]&@@@NoncrossingConnections[Chunks[v]]]/@Tuples[{0,1},2n],1],
(* The following two extra cases represent situations where a unicolor column has fields of the other color before it, hence it can only be followed by the same column. *)
{{Table[0,{2n}],{},{}},{Table[1,{2n}],{},{}}}];
(* Step (2): Discard all uninteresting states. *)
MatrixLabels[n_]:=MatrixLabels[n]=
DeleteCases[MatrixLabels1[n],
{_,{___,c0a_,___,c0b_,___},{___,c1a_,___,c1b_,___}}/;
(Max[c0a]+1==Min[c1a]&&Max[c1a]==Min[c0b]-1&&Max[c0b]<Max[c1b])||
(Max[c0a]+1==Min[c1b]&&Max[c1b]==Min[c0b]-1&&Min[c1a]<Min[c0a])||
(Max[c1a]+1==Min[c0a]&&Max[c0a]==Min[c1b]-1&&Max[c1b]<Max[c0b])||
(Max[c1a]+1==Min[c0b]&&Max[c0b]==Min[c1b]-1&&Min[c0a]<Min[c1a])];

(* Determine which chunks of the new column are connected, based on the connection information of the preceding column. *)
(* Examples:
ComposeConnections[{{1,3},{5,7},{9}},{{1},{3,4,5},{7,8},{10}}];
 ComposeConnections[{{1,3,9},{5},{7}},{{1,2},{5,6},{8,9}}];
 ComposeConnections[{{3,4}},{{1,2},{5,6}}]; *)
ComposeConnections[co_,ch1_]:=
Module[{ch=ch1,s,s1,conn={}},
While[ch=!={},
AppendTo[conn,
FixedPoint[Function[s,
s1=Flatten[Select[co,Intersection[#,s]=!={}&]];
Union[s,Flatten[Select[ch,Intersection[#,s1]=!={}&]]]],First[ch]]];
ch=Select[ch,Intersection[#,Last[conn]]==={}&];
];
Return[conn];
];

(* Returns a sparse square {0,1}-matrix of dimension Length[MatrixLabels]. *)
(* The true transfer matrix (containing powers of x) is XVector[n]*BuildMatrix[n].
Note that each row contains only entries 0 or x^k with the same k. *)
BuildMatrix[n_]:=
Module[{ml=MatrixLabels[n],dim,chs,mlspl,pos,mat,vecs=Tuples[{0,1},2n],row,nl,p1},
dim=Length[ml];
chs=Function[v,Take[Chunks[v],{#,-1,2}]&/@If[v[[1]]===0,{1,2},{2,1}]]/@vecs;
mlspl=Append[SplitBy[Drop[ml,-3],First],Take[ml,-3]];
pos=Most[FoldList[Plus,0,Length/@mlspl]];
mat=SparseArray[{},{1,dim}];
Function[{vv,c0,c1},
row=SparseArray[{},{dim}];
MapThread[Function[{v,ch,mlv,p},
If[(Or@@(MatchQ[v[[#]],{(1)..}]&/@c0))||(Or@@(MatchQ[v[[#]],{(0)..}]&/@c1)),
Which[
MatchQ[v,{(0)..}]&&Length[c1]≤1,row[[-2]]+=1,
MatchQ[v,{(1)..}]&&Length[c0]≤1,row[[-1]]+=1]
,
If[c0===c1==={},
Which[
vv===v===Table[0,{2n}],row[[-2]]+=1,
vv===v===Table[1,{2n}],row[[-1]]+=1]
,
nl=Prepend[MapThread[ComposeConnections,{{c0,c1},ch}],v];
If[(p1=FirstPosition[mlv,nl])=!=Missing["NotFound"],row[[p+p1[[1]]]]+=1];
];
];
],{vecs,chs,mlspl,pos}];
AppendTo[mat,row];
]@@@ml;
Return[Rest[mat]];
];
XVector[n_]:=(x^Count[#[[1]],0])&/@MatrixLabels[n];
(* Encodes possible start configurations (i.e. those with no nontrivial connections). *)
StartVector[n_]:=If[Length[Join[#2,#3]]===Length[Chunks[#1]],x^Count[#1,0],0]&@@@MatrixLabels[n];
(* {0,1}-vector that picks only those configurations that are fully connected. *)
EndVector[n_]:=If[MatchQ[#2,{}|{_}]&&MatchQ[#3,{}|{_}],1,0]&@@@MatrixLabels[n];

n=1

MatrixForm[BuildMatrix[1].DiagonalMatrix[XVector[1]]]

StartVector[1]

EndVector[1]

Expand[StartVector[1].BuildMatrix[1]*XVector[1]]

%.EndVector[1]

n=2

Print/@MatrixLabels[2];

MatrixForm[mat2=BuildMatrix[2]]

StartVector[2]

EndVector[2]

Nest[Expand[(#.mat2)*XVector[2]]&,StartVector[2],3].EndVector[2]

(* Verify this number by a naive, brute-force computation *)
MyConnectedQ[set_]:=FixedPoint[Function[s,Intersection[Join@@({#,#+{1,0},#+{-1,0},#+{0,1},#+{0,-1}}&/@s),set]],{First[set]}]===Sort[set];
With[{t=Tuples[Range[4],2]},
test=Select[Subsets[t,{8}],MyConnectedQ[#]&&MyConnectedQ[Complement[t,#]]&];
Length[test]
]

n=3

Length[MatrixLabels[3]]

Timing[mat3=BuildMatrix[3];]

svec3=StartVector[3];
evec3=EndVector[3];
xvec3=XVector[3];
mat3x=mat3.SparseArray[DiagonalMatrix[xvec3]];

(* Most stupid way: powering the matrix, then multiply with vector *)
Timing[test1=Expand[svec3.MatrixPower[mat3x,5].evec3];]

(* Better: do only matrix-vector multiplications *)
Timing[test2=Nest[Expand[#.mat3x]&,svec3,5].evec3;]

(* Even better: use {0,1} matrix and x-vector *)
Timing[test3=Nest[Expand[(#.mat3)*xvec3]&,svec3,5].evec3;]

(* Using PolynomialMod *)
Timing[Nest[PolynomialMod[(#.mat3)*xvec3,x^19]&,svec3,5].evec3]

(* Using evaluation-interpolation *)
Timing[test4=Expand[InterpolatingPolynomial[Table[xvec0=xvec3/.x->x0;Nest[((#.mat3)*xvec0)&,svec3/.x->x0,5].evec3,{x0,40}],x]]]

test1===test2===test3===test4

(* How much memory is saved by using SparseArray and 0/1's. *)
ByteCount/@{mat3,mat3x,Normal[mat3],Normal[mat3x]}

n=4

(* Size of the transfer matrix if uninteresting states are not discarded. *)
Timing[Length[MatrixLabels1[4]]]

(* Improved version: discarding uninteresting states yields a smaller matrix. *)
Timing[Length[MatrixLabels[4]]]

Timing[mat4=BuildMatrix[4];]

Timing[xvec4=XVector[4];svec4=StartVector[4];evec4=EndVector[4];]

Timing[test1=Nest[Expand[(#.mat4)*xvec4]&,svec4,7].evec4;]

(* Comparison to using a normal matrix (not SparseArray). *)
mat4a=Normal[mat4];
Print[ByteCount/@{mat4,mat4a}];
Timing[test2=Nest[Expand[(#.mat4a)*xvec4]&,svec4,7].evec4;]

(* Putting the x-powers into the matrix (not using SparseArray). *)
mat4a=Normal[mat4].DiagonalMatrix[xvec4];
Print[ByteCount/@{mat4,mat4a}];
Timing[test3=Nest[Expand[#.mat4a]&,svec4,7].evec4;]

(* Using sparse matrix and (palindromic) evaluation-interpolation,
modulo a single, sufficiently large prime. *)
Timing[
pr=NextPrime[2^41,-1];
vals=Table[xvmod=Mod[xvec4/.x->xx,pr];{xx,Nest[Mod[(#.mat4)*xvmod,pr]&,svec4/.x->xx,7].evec4},{xx,2,34}];
vals=PolynomialMod[Join[vals,{1/#1,#2/#1^64}&@@@vals],pr];
test4=PolynomialMod[InterpolatingPolynomial[vals,x,Modulus->pr],pr];
]

(* Sparsity: 96% of the matrix entries are zeros. *)
N[Total[Count[#,0]&/@mat4a]/(Times@@Dimensions[mat4a])]

test1===test2===test3===test4

test1

n=5

(* Size of the transfer matrix if uninteresting states are not discarded. *)
Timing[Length[MatrixLabels1[5]]]

(* Size of the transfer matrix after discarding uninteresting states. *)
Timing[Length[MatrixLabels[5]]]

Timing[mat5=BuildMatrix[5];]

(* Sparsity: 98.4% of the matrix entries are zeros. *)
N[Count[Normal[mat5],0,{2}]/(Times@@Dimensions[mat5])]

Timing[xvec5=XVector[5];svec5=StartVector[5];evec5=EndVector[5];]

Timing[res5=Nest[Expand[(#.mat5)*xvec5]&,svec5,9].evec5;]

(* Interesting: this computation takes half the time, although it is basically the same! *)
mat5a=Transpose[mat5];
Timing[res5a=Nest[Expand[(mat5a.#)*xvec5]&,svec5,9].evec5;]

res5

n=6

(* Size of the transfer matrix if uninteresting states are not discarded. *)
(* The corresponding BuildMatrix took 42h on radon1 (single processor). *)
Timing[Length[MatrixLabels1[6]]]

(* Size of the transfer matrix after discarding uninteresting states. *)
Timing[Length[MatrixLabels[6]]]

BuildMatrix[6] took about 20min on compute-grantley using 15 parallel processes.
For parallelization, we divided the matrix into chunks of 1000 rows.
Matrix (SparseArray) has ByteCount 264,509,768 and 195M file size.

(* Assemble the whole matrix from the parallelized chunks. *)
dim=Length[MatrixLabels[6]];
mat6=Transpose[Join@@Table[Get["mat6_"<>ToString[p]<>".m"],{p,1,dim,1000}]];

xvec=XVector[6];
(* took 2956s on compute-grantley *)
Timing[res6=Nest[Expand[xvec*(mat6.#)]&,StartVector[6],11].EndVector[6];]

(* took 2651s on compute-grantley *)
Timing[res6=Nest[PolynomialMod[xvec*(mat6.#),x^73]&,StartVector[6],11].EndVector[6];]

With evaluation-interpolation and chinese remaindering, it took 20 minutes (without exploiting palindromicity) on compute-grantley, using 3 primes.

res6=Get["res6.m"]

n=7

All computations have been done on radon1 (compute-grantley, Intel Xeon E5-2630v3 processors at 2.4 GHz).

Old version (without discarding uninteresting states):
MatrixLabels[7] took 36s, matrix dimension is 1,848,612.
BuildMatrix[7] took about 46h using 15 parallel processes.
Matrix (SparseArray) has ByteCount 34,944,460,920 and 29G file size.
Matrix-vector multiplications take about 17h for a single prime (one processor).

New version (after discarding uninteresting states):
MatrixLabels[7] took 229s, matrix dimension is 393,878.
BuildMatrix[7] took about 8h using 15 parallel processes.
Matrix (SparseArray) has ByteCount 5,400,473,096 and 4.2G file size.
The matrix contains 337,332,580 nonzero entries, which is 0.2% of all entries.

Timings for the matrix-vector multiplications:
- computing over Z[x] without any tricks took 28.7h
- over Z[x] using PolynomialMod[..., x^99]: took 25h
- using evaluation-interpolation and Chinese remaindering: we used 5 primes , each took 2.5h (without exploiting palindromicity), but actually 4 primes would suffice.

(* Comparing the timings for matrix-vector mults without (blue) and with (orange) PolynomialMod. *)
time1={12,293,1579,3172,4526,5925,7309,8981,10923,12753,14403,16083,17541};
time2={13,302,1817,3724,5218,6648,8035,9426,10424,10967,11148,11152,11157};
ListPlot[{time1,time2},Joined->True]

res7=Get["res7.m"]

Sequence terms

{4,140,161036,14314228378,99704315229167288,56578717186086829451888706,2670113158846160742372913777087464324}/2