Message Boards Message Boards

A Diophantine Benchmark

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, 440

and column sums

2000, 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 >= 0

is 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})

1322852284046903937703609378206302241

cases 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
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 BY: Vitaliy Kaurov
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
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 BY: Vitaliy Kaurov
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract