## CONTINGENCY TABLE and MAGIC SQUARE ENUMERATION

This is a supplemental discussion complementing the paper

by Alexander Barvinok, Alexander Samorodnitsky and Alexander Yong.

Below, we record estimates obtained by the Maple software contingency.txt available here.

Some of this data appears in the paper. Our two main motivations for providing more data here is for further investigation of the following aspects:

(1) We conjecture magic squares are estimated within a factor of N^{O(1)} by the algorithm, or possibly n^{0(1)} or even O(1) (but can only prove N^{O(log(n))}).

(2) We conjecture contingency tables are also well estimated by the algorithm.

(3) The algorithm allows one to "bootstrap" estimates of the number of contingency tables of size n with margins (R,C) to estimates for contingency tables of size n with margins (R',C') where R' and C' are component wise larger than R and C. (So for example, to compute the number of 10x10 magic squares with line sum 100, one can use the estimate for 10x10 magic squares with line sum 50.)

(4) To provide a cross comparisons of various ways to compute/estimate the number of contingency tables

LattE is a computer package for lattice point enumeration, which in particular can be applied to enumeration of contingency tables. For n roughly no larger than 10, one can expect to obtain exact enumeration. This is our first source of comparison (we thank Jesus DeLoera for providing LattE computations):

Comparisons to LattE
5 5 1.51*10^7 2.21*10^7
10 7.96*10^10 7.93*10^10
125 1.55*10^27 1.10*10^27
6 6 4.70*10^11 6.02*10^11
12 7.12*10^16 2.28*10^17
36 2.21*10^27 5.62*10^27
216 3.76*10^46 3.07*10^46
7 7 1.70*10^17 2.16*10^17
14 1.75*10^25 2.46*10^25
49 3.00*10^42 4.00*10^42
343 0.39*10^72 1.28*10^72

In addition, we compared out results to cases where the Diaconis-Efron heuristic is available (when t is "large" compared to n). (This heuristic is implemented into contingency.txt as Diaconis_Efron(R,C); where R and C are the row and column margins.)

Comparisons to Diaconis-Efron
15 100 0.765*10^237 2.71*10^237
20 70 0.567*10^350 0.791*10^351

Finally, we compared to cases where the asymptotic formula of A.~B\'{e}k\'{e}ssy, P.~B\'{e}k\'{e}ssy and J.~Koml\'{o}s is available (roughly when t is small compared to n. (Their formula is ${N! e^{(t-1)^2/2} \over (t!)^{2n}}$):

Comparisons to asymptotic formula
25 5 2.89*10^108 6.17*10^108
30 5 0.1485*10^142 0.302256*10^142

Here are some cases where no known comparisons are possible

Cases without comparison
12 8 2.72*10^48 ?
12 20 4.5458*10^81 ?
15 20 0.83*10^120 ?

Here are some non magic squares cases, including the 4x4 "hard" example

General contingency tables
r c R C "contingency.txt" estimate Answer
4 4 [220,215,93,64] [108,286,71,127] 1.02*10^15 1.22*10^15 (LattE)
12 12 [80,60,80,60,80,60,80,60,80,60,80,60] [80,60,80,60,80,60,80,60,80,60,80,60] 3.00*10^137 3.29*10^137 (Diaconis-Efron)

## Computation time:

All of the above computations were done in the order of a few days (more specific details available upon request). The mixing time to get stable answers grows with n and t. However, the fact that the algorithm is highly parallelizable and the opportunity to "bootstrap" mentioned above can bring down the computation time significantly. (For example, to compute 7x7--343 above, we actually can compute all the other 7x7 cases presented above, at once.)

# Future work

We would like to sharpen the theoretical understanding of our algorithm in terms of mixing time results for the Markov chain Monte Carlo hit-and-run algorithm used, as well as to generalize such results to more general cases, e.g., general contingency tables and integer flows.