{"id":378,"date":"2012-11-05T12:17:43","date_gmt":"2012-11-05T17:17:43","guid":{"rendered":"https:\/\/blogs.mathworks.com\/cleve\/?p=378"},"modified":"2016-12-05T13:55:02","modified_gmt":"2016-12-05T18:55:02","slug":"magic-squares-part-2-algorithms","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/cleve\/2012\/11\/05\/magic-squares-part-2-algorithms\/","title":{"rendered":"Magic Squares, Part 2, Algorithms"},"content":{"rendered":"&nbsp;\r\n<div class=\"content\"><!--introduction--><img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/images\/cleve\/magicshow.gif\" alt=\"\" hspace=\"5\" vspace=\"5\" \/>\r\n\r\nThe magic squares of odd order generated by MATLAB show a pattern with increasing elements generally moving diagonally up and to the right.\r\n\r\n<!--\/introduction-->\r\n<h3>Contents<\/h3>\r\n<div>\r\n<ul>\r\n \t<li><a href=\"#fa048e2a-6584-4194-abf5-158e8dc1143f\">Three Cases<\/a><\/li>\r\n \t<li><a href=\"#b872e807-a87a-44e4-9b84-a9885b78a77a\">Odd Order<\/a><\/li>\r\n \t<li><a href=\"#adc5d2f0-fe61-4877-8fd4-f1c8177476ad\">A New Algorithm<\/a><\/li>\r\n \t<li><a href=\"#beb9d329-7820-4467-99d4-4a3b426027e4\">Doubly Even Order<\/a><\/li>\r\n \t<li><a href=\"#6cb69d92-50a2-4fde-922d-424a27e9f96f\">Singly Even Order<\/a><\/li>\r\n \t<li><a href=\"#1316ab63-3184-4e3c-8342-39233c336717\">Further Reading<\/a><\/li>\r\n<\/ul>\r\n<\/div>\r\n<h4>Three Cases<a name=\"fa048e2a-6584-4194-abf5-158e8dc1143f\"><\/a><\/h4>\r\nThe algorithms used by MATLAB for generating magic squares of order <i>n<\/i> fall into three cases:\r\n<div>\r\n<ul>\r\n \t<li>odd, <i>n<\/i> is odd.<\/li>\r\n \t<li>doubly-even, <i>n<\/i> is divisible by 4.<\/li>\r\n \t<li>singly-even, <i>n<\/i> is divisible by 2, but not by 4.<\/li>\r\n<\/ul>\r\n<\/div>\r\n<h4>Odd Order<a name=\"b872e807-a87a-44e4-9b84-a9885b78a77a\"><\/a><\/h4>\r\nThe best known algorithm for generating magic squares of odd order is de La Loubere's method. Simon de La Loubere was the French ambassador to Siam in the late 17th century. I sometimes refer to his method as the \"nor'easter algorithm\", after the winter storms that move northeasterly up the coast of New England. You can see why if you follow the integers sequentially through <tt>magic(9)<\/tt>.\r\n<pre>   47    58    69    80     1    12    23    34    45\r\n   57    68    79     9    11    22    33    44    46\r\n   67    78     8    10    21    32    43    54    56\r\n   77     7    18    20    31    42    53    55    66\r\n    6    17    19    30    41    52    63    65    76\r\n   16    27    29    40    51    62    64    75     5\r\n   26    28    39    50    61    72    74     4    15\r\n   36    38    49    60    71    73     3    14    25\r\n   37    48    59    70    81     2    13    24    35<\/pre>\r\nThe integers from $1$ to $n^2$ are inserted along diagonals, starting in the middle of first row and heading in a northeasterly direction. When you go off an edge of the array, which you do at the very first step, continue from the opposite edge. When you bump into a cell that is already occupied, drop down one row and continue.\r\n\r\nWe used this algorithm in MATLAB for many years. Here is the code.\r\n<pre>   A = zeros(n,n);\r\n   i = 1;\r\n   j = (n+1)\/2;\r\n   for k = 1:n^2\r\n      is = i;\r\n      js = j;\r\n      A(i,j) = k;\r\n      i = n - rem(n+1-i,n);\r\n      j = rem(j,n) + 1;\r\n      if A(i,j) ~= 0\r\n         i = rem(is,n) + 1;\r\n         j = js;\r\n      end\r\n   end<\/pre>\r\nA big difficulty with this algorithm and resulting program is that it inserts the elements one at a time---it cannot be <i>vectorized<\/i>.\r\n<h4>A New Algorithm<a name=\"adc5d2f0-fe61-4877-8fd4-f1c8177476ad\"><\/a><\/h4>\r\nA few years ago we discovered an algorithm for generating the same magic squares of odd order as de La Loubere's method, but with just four MATLAB matrix operations.\r\n<pre>   [I,J] = ndgrid(1:n);\r\n   A = mod(I+J+(n-3)\/2,n);\r\n   B = mod(I+2*J-2,n);\r\n   M = n*A + B + 1;<\/pre>\r\nLet's see how this works with <tt>n = 5<\/tt>. The statement\r\n<pre>   [I,J] = ndgrid(1:n)<\/pre>\r\nproduces a pair of matrices whose elements are just the row and column indices, <tt>i<\/tt> and <tt>j<\/tt>.\r\n<pre>   I =\r\n        1     1     1     1     1\r\n        2     2     2     2     2\r\n        3     3     3     3     3\r\n        4     4     4     4     4\r\n        5     5     5     5     5<\/pre>\r\n<pre>   J =\r\n        1     2     3     4     5\r\n        1     2     3     4     5\r\n        1     2     3     4     5\r\n        1     2     3     4     5\r\n        1     2     3     4     5<\/pre>\r\nUsing these indices, we generate two more matrices. The statements\r\n<pre>   A = mod(I+J+1,n)\r\n   B = mod(I+2*J-2,n)<\/pre>\r\nproduce\r\n<pre>   A =\r\n        3     4     0     1     2\r\n        4     0     1     2     3\r\n        0     1     2     3     4\r\n        1     2     3     4     0\r\n        2     3     4     0     1<\/pre>\r\n<pre>   B =\r\n        1     3     0     2     4\r\n        2     4     1     3     0\r\n        3     0     2     4     1\r\n        4     1     3     0     2\r\n        0     2     4     1     3<\/pre>\r\nBoth <tt>A<\/tt> and <tt>B<\/tt> are fledgling magic squares. They have equal row, column and diagonal sums. But their elements are not the integers from $1$ to $n^2$. Each has duplicated elements between $0$ and $n-1$. The final statement\r\n<pre>   M = n*A+B+1<\/pre>\r\nproduces a matrix whose elements are integers between $1$ and $n^2$ and which has equal row, column and diagonal sums. What is not obvious, but is true, is that there are no duplicates. So <tt>M<\/tt> must contain <i>all<\/i> of the integers between $1$ and $n^2$ and consequently is a magic square.\r\n<pre>   M =\r\n       17    24     1     8    15\r\n       23     5     7    14    16\r\n        4     6    13    20    22\r\n       10    12    19    21     3\r\n       11    18    25     2     9<\/pre>\r\n<h4>Doubly Even Order<a name=\"beb9d329-7820-4467-99d4-4a3b426027e4\"><\/a><\/h4>\r\nThe doubly-even algorithm is also short and sweet, and tricky.\r\n<pre>   M = reshape(1:n^2,n,n)';\r\n   [I,J] = ndgrid(1:n);\r\n   K = fix(mod(I,4)\/2) == fix(mod(J,4)\/2);\r\n   M(K) = n^2+1 - M(K);<\/pre>\r\nLet's look at our friend <tt>magic(4)<\/tt>. The matrix <tt>M<\/tt> is initially just the integers from 1 to 16 stored sequentially in 4 rows.\r\n<pre>   M =\r\n        1     2     3     4\r\n        5     6     7     8\r\n        9    10    11    12\r\n       13    14    15    16<\/pre>\r\nThe logical array <tt>K<\/tt> is true for half of the indices and false for the other half in a pattern like this.\r\n<pre>   K =\r\n        1     0     0     1\r\n        0     1     1     0\r\n        0     1     1     0\r\n        1     0     0     1<\/pre>\r\nThe elements where <tt>K<\/tt> is false, that is 0, are left alone.\r\n<pre>        .     2     3     .\r\n        5     .     .     8\r\n        9     .     .    12\r\n        .    14    15     .<\/pre>\r\nThe elements where <tt>K<\/tt> is true, that is 1, are reversed.\r\n<pre>       16     .     .    13\r\n        .    11    10     .\r\n        .     7     6     .\r\n        4     .     .     1<\/pre>\r\nThe final result merges these two matrices to produce the magic square.\r\n<h4>Singly Even Order<a name=\"6cb69d92-50a2-4fde-922d-424a27e9f96f\"><\/a><\/h4>\r\nThe algorithm for singly even order is the most complicated and so we will give just a glimpse of how it works. If <i>n<\/i> is singly even, then $n\/2$ is odd and <tt>magic(n)<\/tt> can be constructed from four copies of <tt>magic(n\/2)<\/tt>. For example, <tt>magic(10)<\/tt> is obtained from <tt>A = magic(5)<\/tt> by forming a block matrix.\r\n<pre>   [  A      A+50\r\n     A+75    A+25 ]<\/pre>\r\nThe column sums are all equal because <tt>sum(A) + sum(A+75)<\/tt> equals <tt>sum(A+50) + sum(A+25)<\/tt>. But the rows sums are not quite right. The algorithm must finish by doing a few swaps of pieces of rows to clean up the row sums. For the details, issue the command.\r\n<pre>   type magic<\/pre>\r\n<h4>Further Reading<a name=\"1316ab63-3184-4e3c-8342-39233c336717\"><\/a><\/h4>\r\nThis post of a revision of a portion of the <a href=\"https:\/\/www.mathworks.com\/content\/dam\/mathworks\/mathworks-dot-com\/moler\/exm\/chapters\/magic.pdf\">Magic Squares<\/a> chapter of <a href=\"https:\/\/www.mathworks.com\/moler\/exm\"><i>Experiments with MATLAB<\/i><\/a>.\r\n\r\n<script>\/\/ <![CDATA[\r\nfunction grabCode_6a5d1e3258594e68a14447491e460481() {\r\n        \/\/ Remember the title so we can use it in the new page\r\n        title = document.title;\r\n\r\n        \/\/ Break up these strings so that their presence\r\n        \/\/ in the Javascript doesn't mess up the search for\r\n        \/\/ the MATLAB code.\r\n        t1='6a5d1e3258594e68a14447491e460481 ' + '##### ' + 'SOURCE BEGIN' + ' #####';\r\n        t2='##### ' + 'SOURCE END' + ' #####' + ' 6a5d1e3258594e68a14447491e460481';\r\n    \r\n        b=document.getElementsByTagName('body')[0];\r\n        i1=b.innerHTML.indexOf(t1)+t1.length;\r\n        i2=b.innerHTML.indexOf(t2);\r\n \r\n        code_string = b.innerHTML.substring(i1, i2);\r\n        code_string = code_string.replace(\/REPLACE_WITH_DASH_DASH\/g,'--');\r\n\r\n        \/\/ Use \/x3C\/g instead of the less-than character to avoid errors \r\n        \/\/ in the XML parser.\r\n        \/\/ Use '\\x26#60;' instead of '<' so that the XML parser\r\n        \/\/ doesn't go ahead and substitute the less-than character. \r\n        code_string = code_string.replace(\/\\x3C\/g, '\\x26#60;');\r\n\r\n        copyright = 'Copyright 2012 The MathWorks, Inc.';\r\n\r\n        w = window.open();\r\n        d = w.document;\r\n        d.write('\r\n\r\n<pre>\\n');\r\n        d.write(code_string);\r\n\r\n        \/\/ Add copyright line at the bottom if specified.\r\n        if (copyright.length > 0) {\r\n            d.writeln('');\r\n            d.writeln('%%');\r\n            if (copyright.length > 0) {\r\n                d.writeln('% _' + copyright + '_');\r\n            }\r\n        }\r\n\r\n        d.write('<\/pre>\r\n\r\n\r\n\\n');\r\n\r\n        d.title = title + ' (MATLAB code)';\r\n        d.close();\r\n    }\r\n\/\/ ]]><\/script>\r\n<p style=\"text-align: right; font-size: xx-small; font-weight: lighter; font-style: italic; color: gray;\">\r\n<a><span style=\"font-size: x-small; font-style: italic;\">Get\r\nthe MATLAB code<noscript>(requires JavaScript)<\/noscript><\/span><\/a>\r\n\r\nPublished with MATLAB\u00ae R2012b<\/p>\r\n<p class=\"footer\">\r\nPublished with MATLAB\u00ae R2012b<\/p>\r\n\r\n<\/div>\r\n<!--\r\n6a5d1e3258594e68a14447491e460481 ##### SOURCE BEGIN #####\r\n%% Magic Squares, Part 2, Algorithms\r\n%\r\n% <<magicshow.gif>>\r\n%\r\n% The magic squares of odd order generated by MATLAB show a pattern\r\n% with increasing elements generally moving diagonally up and to the right.\r\n\r\n%% Three Cases\r\n% The algorithms used by MATLAB for generating magic squares of order _n_\r\n% fall into three cases:\r\n%\r\n% * odd, _n_ is odd.\r\n% * doubly-even, _n_ is divisible by 4.\r\n% * singly-even, _n_ is divisible by 2, but not by 4.\r\n\r\n%% Odd Order\r\n% The best known algorithm for generating magic squares of odd order is\r\n% de La Loubere's method.  Simon de La Loubere was the French ambassador\r\n% to Siam in the late 17th century.  I sometimes refer to his method as the\r\n% \"nor'easter algorithm\", after the winter storms that move northeasterly\r\n% up the coast of New England.  You can see why if you follow\r\n% the integers sequentially through |magic(9)|.\r\n%\r\n%     47    58    69    80     1    12    23    34    45\r\n%     57    68    79     9    11    22    33    44    46\r\n%     67    78     8    10    21    32    43    54    56\r\n%     77     7    18    20    31    42    53    55    66\r\n%      6    17    19    30    41    52    63    65    76\r\n%     16    27    29    40    51    62    64    75     5\r\n%     26    28    39    50    61    72    74     4    15\r\n%     36    38    49    60    71    73     3    14    25\r\n%     37    48    59    70    81     2    13    24    35\r\n%\r\n%\r\n% The integers from $1$ to $n^2$ are inserted along diagonals, starting\r\n% in the middle of first row and heading in a northeasterly direction.\r\n% When you go off an edge of the array, which you do at the very first step,\r\n% continue from the opposite edge.  When you bump into a cell that is\r\n% already occupied, drop down one row and continue.\r\n%\r\n% We used this algorithm in MATLAB for many years.\r\n% Here is the code.\r\n%\r\n%     A = zeros(n,n);\r\n%     i = 1;\r\n%     j = (n+1)\/2;\r\n%     for k = 1:n^2\r\n%        is = i;\r\n%        js = j;\r\n%        A(i,j) = k;\r\n%        i = n - rem(n+1-i,n);\r\n%        j = rem(j,n) + 1;\r\n%        if A(i,j) ~= 0\r\n%           i = rem(is,n) + 1;\r\n%           j = js;\r\n%        end\r\n%     end\r\n%\r\n%\r\n% A big difficulty with this algorithm and resulting program is that\r\n% it inserts the elements one at a timeREPLACE_WITH_DASH_DASH-it cannot be _vectorized_.\r\n\r\n%% A New Algorithm\r\n% A few years ago we discovered an algorithm for generating\r\n% the same magic squares of odd order as de La Loubere's method,\r\n% but with just four MATLAB matrix operations.\r\n%\r\n%     [I,J] = ndgrid(1:n);\r\n%     A = mod(I+J+(n-3)\/2,n);\r\n%     B = mod(I+2*J-2,n);\r\n%     M = n*A + B + 1;\r\n%\r\n%\r\n% Let's see how this works with |n = 5|.  The statement\r\n%\r\n%     [I,J] = ndgrid(1:n)\r\n%\r\n% produces a pair of matrices whose elements are just the row\r\n% and column indices, |i| and |j|.\r\n%\r\n%     I =\r\n%          1     1     1     1     1\r\n%          2     2     2     2     2\r\n%          3     3     3     3     3\r\n%          4     4     4     4     4\r\n%          5     5     5     5     5\r\n%\r\n%     J =\r\n%          1     2     3     4     5\r\n%          1     2     3     4     5\r\n%          1     2     3     4     5\r\n%          1     2     3     4     5\r\n%          1     2     3     4     5\r\n%\r\n% Using these indices, we generate two more matrices.\r\n% The statements\r\n%\r\n%     A = mod(I+J+1,n)\r\n%     B = mod(I+2*J-2,n)\r\n%\r\n% produce\r\n%\r\n%     A =\r\n%          3     4     0     1     2\r\n%          4     0     1     2     3\r\n%          0     1     2     3     4\r\n%          1     2     3     4     0\r\n%          2     3     4     0     1\r\n%\r\n%     B =\r\n%          1     3     0     2     4\r\n%          2     4     1     3     0\r\n%          3     0     2     4     1\r\n%          4     1     3     0     2\r\n%          0     2     4     1     3\r\n%\r\n% Both |A| and |B| are fledgling magic squares.\r\n% They have equal row, column and diagonal sums.\r\n% But their elements are not the integers from $1$ to $n^2$.\r\n% Each has duplicated elements between $0$ and $n-1$.\r\n% The final statement\r\n%\r\n%     M = n*A+B+1\r\n%\r\n% produces a matrix whose elements are integers between $1$ and $n^2$\r\n% and which has equal row, column and diagonal sums.  What is not obvious,\r\n% but is true, is that there are no duplicates.  So |M| must contain\r\n% _all_ of the integers between $1$ and $n^2$ and consequently\r\n% is a magic square.\r\n%\r\n%     M =\r\n%         17    24     1     8    15\r\n%         23     5     7    14    16\r\n%          4     6    13    20    22\r\n%         10    12    19    21     3\r\n%         11    18    25     2     9\r\n\r\n%% Doubly Even Order\r\n% The doubly-even algorithm is also short and sweet, and tricky.\r\n%\r\n%     M = reshape(1:n^2,n,n)';\r\n%     [I,J] = ndgrid(1:n);\r\n%     K = fix(mod(I,4)\/2) == fix(mod(J,4)\/2);\r\n%     M(K) = n^2+1 - M(K);\r\n%\r\n% Let's look at our friend |magic(4)|.\r\n% The matrix |M| is initially just the integers from 1 to 16\r\n% stored sequentially in 4 rows.\r\n%\r\n%     M =\r\n%          1     2     3     4\r\n%          5     6     7     8\r\n%          9    10    11    12\r\n%         13    14    15    16\r\n%\r\n% The logical array |K| is true for half of the indices and\r\n% false for the other half in a pattern like this.\r\n%\r\n%     K =\r\n%          1     0     0     1\r\n%          0     1     1     0\r\n%          0     1     1     0\r\n%          1     0     0     1\r\n%\r\n% The elements where |K| is false, that is 0, are left alone.\r\n%\r\n%          .     2     3     .\r\n%          5     .     .     8\r\n%          9     .     .    12\r\n%          .    14    15     .\r\n%\r\n% The elements where |K| is true, that is 1, are reversed.\r\n%\r\n%         16     .     .    13\r\n%          .    11    10     .\r\n%          .     7     6     .\r\n%          4     .     .     1\r\n%\r\n% The final result merges these two matrices to produce\r\n% the magic square.\r\n\r\n%% Singly Even Order\r\n% The algorithm for singly even order is the most complicated\r\n% and so we will give just a glimpse of how it works.\r\n% If _n_ is singly even, then $n\/2$ is odd and |magic(n)| can be\r\n% constructed from four copies of |magic(n\/2)|.  For example,\r\n% |magic(10)| is obtained from |A = magic(5)| by forming\r\n% a block matrix.\r\n%\r\n%     [  A      A+50\r\n%       A+75    A+25 ]\r\n%\r\n% The column sums are all equal because\r\n% |sum(A) + sum(A+75)| equals |sum(A+50) + sum(A+25)|.  But the\r\n% rows sums are not quite right.  The algorithm must finish\r\n% by doing a few swaps of pieces of rows to clean up the row\r\n% sums.  For the details, issue the command.\r\n%\r\n%     type magic\r\n%\r\n\r\n%% Further Reading\r\n% This post of a revision of a portion of the\r\n% <https:\/\/www.mathworks.com\/moler\/exm\/chapters\/magic.pdf Magic Squares>\r\n% chapter of\r\n% <https:\/\/www.mathworks.com\/moler\/exm _Experiments with MATLAB_>.\r\n\r\n##### SOURCE END ##### 6a5d1e3258594e68a14447491e460481\r\n-->","protected":false},"excerpt":{"rendered":"... <a class=\"read-more\" href=\"https:\/\/blogs.mathworks.com\/cleve\/2012\/11\/05\/magic-squares-part-2-algorithms\/\">read more >><\/a><\/p>","protected":false},"author":78,"featured_media":0,"comment_status":"open","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":[],"categories":[9],"tags":[],"_links":{"self":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/378"}],"collection":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/users\/78"}],"replies":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/comments?post=378"}],"version-history":[{"count":8,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/378\/revisions"}],"predecessor-version":[{"id":2178,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/378\/revisions\/2178"}],"wp:attachment":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/media?parent=378"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/categories?post=378"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/tags?post=378"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}