  We present a randomized algorithm, which, given positive integers n and t and
a real number 0< epsilon <1, computes the number Sigma(n, t) of n x n
non-negative integer matrices (magic squares) with the row and column sums
equal to t within relative error epsilon. The computational complexity of the
algorithm is polynomial in 1/epsilon and quasi-polynomial in N=nt, that is, of
the order N^{log N}. A simplified version of the algorithm works in time
polynomial in 1/epsilon and N and estimates Sigma(n,t) within a factor of
N^{log N}. This simplified version has been implemented. We present results of
the implementation, state some conjectures, and discuss possible
generalizations.
