Semimagic Square

Question

I have to create function which generate random NxN magic square with all sums of columns and rows equal to S. Variables N and S have to be users input. I've already create function which draws squares and breaks the loop when meet given conditions. It's working fine for small squars like 3x3 or 4x4 but it`s extremly unefficient for larger ones. All numbers inside the matrix have to be integer, there is no need diagonals be equal rows sums. Any suggestion how to solve this problem?

Answer 1

As any magic square is a semimagic one too (which doesn't have to have diagonals to sum up to the magic constant), to answer the question it's enough to provide an efficient algorithm for generating truly magic squares. One of classical algorithms to solve magic squares is the following (citation by Robin K. S. Hankin):

In a square of order 4m, shade the long major diagonal. Then shade all major diagonals distant by a multiple of 4 cells from the long diagonal. Do the same with the minor diagonals. Then, starting with “1” at the top left corner and proceeding from left to right and top to bottom, count from 1 to n^2, filling in the shaded squares with the appropriate number and omitting the unshaded ones. Fill in the remaining (unshaded) squares in the same way, starting at the lower right corner, moving leftwards and upward.

You might find R package magic a useful resource for studying efficient magic square generation algorithms. It becomes clear from the package sources that the problem can be divided into three separate cases:

"magic" <-
function (n) 
{
    if(length(n)>1){
      return(sapply(n,match.fun(sys.call()[[1]])))
    }
    n <- round(n)
    if (n == 2) {
        stop("Normal magic squares of order 2 do not exist")
    }
    if (n%%2 == 1) {
        return(as.standard(magic.2np1(floor(n/2))))
    }
    if (n%%4 == 0) {
        return(as.standard(magic.4n(round(n/4))))
    }
    if (n%%4 == 2) {
        return(as.standard(magic.4np2(round((n - 2)/4))))
    }
    stop("This cannot happen")
}

The execution time of magic is surprisingly small even for large n:

> system.time(result <- magic(2017))
   user  system elapsed 
  1.256   0.316   1.573

Finally, if you really want just a semimagic square instead of a truly magic one, you can always exchange any two columns of resulting square matrix. This operation wouldn't affect sums in rows and columns, but will definitely break sums of diagonals.