# A Diophantine Benchmark

Posted 10 years ago
6704 Views
|
3 Replies
|
4 Total Likes
|
 In 1997, the following problem was posed to test a programme. How many different 4 x4 matrices having nonnegative integer entries have row sums 1000, 9000, 3000, 440and column sums2000, 4000, 200, 7240 ?The answer is exactly 160256589785681535184841, and it was found in 2 seconds real time on a typical workstation at MSRI, Berkeley, in 1997!. How did they came up with that number so rapidly?The following is one of those matrices Timing[MatrixForm[{{x1, x2, x3, x4}, {x5, x6, x7, x8}, {x9, x10, x11,       x12}, {x13, x14, x15, x16}} /.     First@FindInstance[(x1 + x2 + x3 + x4 ==          1000) && (x5 + x6 + x7 + x8 ==          9000) && (x9 + x10 + x11 + x12 ==          3000) && (x13 + x14 + x15 + x16 ==          440) && (x1 + x5 + x9 + x13 ==          2000) && (x2 + x6 + x10 + x14 ==          4000) && (x3 + x7 + x11 + x15 ==         200) && (x4 + x8 + x12 + x16 == 7240) && x1 >= 0 && x2 >= 0 &&       x3 >= 0 && x4 >= 0 && x5 >= 0 && x6 >= 0 && x7 >= 0 &&       x8 >= 0 && x9 >= 0 && x10 >= 0 && x11 >= 0 && x12 >= 0 &&       x13 >= 0 && x14 >= 0 && x15 >= 0 && x16 >= 0, {x1, x2, x3, x4,       x5, x6, x7, x8, x9, x10, x11, x12, x13, x14, x15, x16},      Integers]]]if we divide all the sums by 40 (which is the gcd of the given sums) we get the following example  Timing[MatrixForm[{{x1, x2, x3, x4}, {x5, x6, x7, x8}, {x9, x10, x11,       x12}, {x13, x14, x15, x16}} /.     First@FindInstance[(x1 + x2 + x3 + x4 ==          25) && (x5 + x6 + x7 + x8 == 225) && (x9 + x10 + x11 + x12 ==          75) && (x13 + x14 + x15 + x16 == 11) && (x1 + x5 + x9 + x13 ==          50) && (x2 + x6 + x10 + x14 ==          100) && (x3 + x7 + x11 + x15 == 5) && (x4 + x8 + x12 + x16 ==          181) && x1 >= 0 && x2 >= 0 && x3 >= 0 && x4 >= 0 && x5 >= 0 &&        x6 >= 0 && x7 >= 0 && x8 >= 0 && x9 >= 0 && x10 >= 0 &&       x11 >= 0 && x12 >= 0 && x13 >= 0 && x14 >= 0 && x15 >= 0 &&       x16 >= 0, {x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12,       x13, x14, x15, x16}, Integers]]]Note this solution is not the previous one divided by 40.In looking for a matrix having the first row adding to 25,  the above row {0, 22, 0, 3} must be among the following 3276 FrobeniusSolve[{1, 1, 1, 1}, 25]and there are 1949476 for the second row adding to 225.  There is no need to generate all solutions to know how many are there, we know that, in general, the number of solutions to the Diophantine equation x1 + x2 + x3 + x4 = n ,   x1, x2, x3, x4 integer values >= 0is Binomial[n + 3, n] so, in orther to test  all rows complying with their sums to see whether they comply with the sum of the columns as well we need to examine Times @@ (Binomial[# + 3, #] & /@ {1000, 9000, 3000, 440})1322852284046903937703609378206302241cases to end up with the amount quoted at the beginning. A massive task. So, again, how did they do it ?.As I do not know the answer, I thought I could explain something to give an idea of what is involved by examining the complexities of an easier problem. Let H (n, r) denote the number of square matrices n x n with non - negative elements wherein each row and column sum is equal to r. It is easy to see that H (2, r) = r + 1.  The problem to compute H (n, r) goes back to McMahon (1916) who showed that H (3, r) is equal to Binomial(r+2, 2)+3 Binomial(r+3, 4). This quantity grows as:Table[Binomial[r + 2, 2] + 3 Binomial[r + 3, 4], {r, 20}]{6, 21, 55, 120, 231, 406, 666, 1035, 1540, 2211, 3081, 4186, 5565, 7260, 9316, 11781, 14706, 18145, 22155, 26796}The 21 above corresponds to the number of elements in the following list MatrixForm /@ ({{x1, x2, x3}, {x4, x5, x6}, {x7, x8, x9}} /.    Solve[(x1 + x2 + x3 == x4 + x5 + x6 == x7 + x8 + x9 ==        x1 + x4 + x7 == x2 + x5 + x8 == x3 + x6 + x9 == 2) && x1 >= 0 &&      x2 >= 0 && x3 >= 0 && x4 >= 0 && x5 >= 0 && x6 >= 0 && x7 >= 0 &&      x8 >= 0 && x9 >= 0, {x1, x2, x3, x4, x5, x6, x7, x8, x9},     Integers])Length[%]by the way, the 9316 includes magic squares . if one wants to find out how many magic squares 3 x3 are there (disregarding equivalent ones under rotation/reflections), a better method is the following MatrixForm[Partition[#, 3]] & /@  Select[Permutations[   Range[9]], ({{x1, x2, x3}, {x4, x5, x6}, {x7, x8, x9}} =      Partition[#,       3]; (x1 + x2 + x3 == x4 + x5 + x6 == x7 + x8 + x9 ==        x1 + x4 + x7 == x2 + x5 + x8 == x3 + x6 + x9 == 15) && (x1 ==        Min[x1, x3, x7, x9]) && (x3 < x7)) &]What could be the formula for H (4, r)?.  Stanley and Ehrhart showed that H (n, r) is a polynomial in r of degree (n - 1)^2. For example, H (3, r) is  Simplify[Binomial[r + 2, 2] + 3 Binomial[r + 3, 4]]To obtain the expression for H(4,r) we proceed by interpolating the following values.  Timing[Table[   Length[({{x1, x2, x3, x4}, {x5, x6, x7, x8}, {x9, x10, x11,         x12}, {x13, x14, x15, x16}} /.       Solve[(x1 + x2 + x3 + x4 == x5 + x6 + x7 + x8 ==           x9 + x10 + x11 + x12 == x13 + x14 + x15 + x16 ==           x1 + x5 + x9 + x13 == x2 + x6 + x10 + x14 ==           x3 + x7 + x11 + x15 == x4 + x8 + x12 + x16 == i) && x1 >= 0 &&         x2 >= 0 && x3 >= 0 && x4 >= 0 && x5 >= 0 && x6 >= 0 &&         x7 >= 0 && x8 >= 0 && x9 >= 0 && x10 >= 0 && x11 >= 0 &&        x12 >= 0 && x13 >= 0 && x14 >= 0 && x15 >= 0 && x16 >= 0, {x1,        x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12, x13, x14, x15,        x16}, Integers])], {i, 0, 8}]]and the answer for H (4, r) is Expand[InterpolatingPolynomial[{{0, 1}, {1, 24}, {2, 282}, {3,     2008}, {4, 10147}, {5, 40176}, {6, 132724}, {7, 381424}, {8,     981541}}, r]]To date not even H (7, r) is known.
3 Replies
Sort By:
Posted 10 years ago
 Interesting. This reference seems to be quite popular online. I found Google Books link to the exact page where one can read up a bit more details and see the MSRI computation you mentioned in the footnote on the page 287.
Posted 10 years ago
 Great post, Jaime! Could you please post some links or references to where we can learn about that original computation done in MSRI, Berkeley 1997?
Posted 10 years ago
 Here is the reference. Unfortunately it only mentions the problem casually and moves on to other rather theoretical things ... apart from that I have no other reference.Approximate counting. D. Welsh . Surveys in Combinatorics, 1997. R. A. Bailey ed. Cambridge University Press 1997, pp. 287 - 323