Generate a random magic square on a large field If consecutive natural numbers 1 .. n^2 are placed in a n*n grid so that sum of numbers in each column and in each row is a constant then we call it a magic square.
I would like to generate magic squares of size 16x16 and 256x256 for not skewed byte mapping. Maps arranged as magic squares tend to have better distribution properties. Algo should produce a map that depends on the seed, rather than yield the same result each time.

expectedSum = ((1+n*n)*n)/2). I can shuffle an array using some form of SecureRandom that accepts a seed. Then in case of a small square I can go for a brute force for loop.I can just brute force 3x3 magic square. Considering symmetries I need to check less than 8*(7!) or about 40k possibilities. 

This approach wouldn't work for 16x16 board because 256! > 10^500 is quite large. My best guess is to find two row_column pairs such that sums along row and column are greater than the expected sum  in first pair and lower in second pair. If shared element in first pair is greater than in the second pair - swap elements.  This algo seems better than bruteforce (can complete in reasonable time), but there is no guarantee that it halts. Consider 2x2 square to see how it fails. 

What is a good Monte-Carlo approach to solve a large magic square?
Here's a decent algorithm that will quickly generate large magic squares.

First Step: Establish Rows
1. Put all numbers in a bucket.
2. Loop through the rows, 0 to n.
3. Take out n (where n is the size of row) numbers randomly from the bucket.
4. Sum those numbers. If they don't equal the solution sum, goto 3.
5. If they do equal the solution sum, stick them into your grid at the row you are on, remove the numbers used from the bucket, and goto 2.

You'll now -- very quickly -- have a square with each row having the correct sum.

Second Step: Establish Columns
1. Loop through the columns, 0 to n.
2. Shuffle the elements in each row from the beginning of the row + the column you are on, to the end of the row. (That is, if we are working on the first column, shuffle all the elements in each row. If we are working on the second column, shuffle all the elements but those in the first column, etc.)
3. Sum the column. If it doesn't equal solution sum, goto 2. 
4. Goto 1 and proceed to the next column.

Now, in this step, you'll need a bail counter... that is, after 1000 tries or so, give up and return to the First Step. The number of attempts needs to increase as n increases... For a 25x25 grid, for example, let it try 100,000 times or so before giving up and creating a new set of rows.

Third Step: Establishing Diagonals
1. Swap any two rows and swap any two columns.
2. Check diagonals sums.
3. If the sum equals solution sum, you're finished. If not, go to 1.

Once again, stick in a bail counter... after 10,000 tries or so, go back to First Step.