1########################################################################
2##
3## Copyright (C) 1999-2021 The Octave Project Developers
4##
5## See the file COPYRIGHT.md in the top-level directory of this
6## distribution or <https://octave.org/copyright/>.
7##
8## This file is part of Octave.
9##
10## Octave is free software: you can redistribute it and/or modify it
11## under the terms of the GNU General Public License as published by
12## the Free Software Foundation, either version 3 of the License, or
13## (at your option) any later version.
14##
15## Octave is distributed in the hope that it will be useful, but
16## WITHOUT ANY WARRANTY; without even the implied warranty of
17## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
18## GNU General Public License for more details.
19##
20## You should have received a copy of the GNU General Public License
21## along with Octave; see the file COPYING.  If not, see
22## <https://www.gnu.org/licenses/>.
23##
24########################################################################
25
26## -*- texinfo -*-
27## @deftypefn {} {} magic (@var{n})
28##
29## Create an @var{n}-by-@var{n} magic square.
30##
31## A magic square is an arrangement of the integers @code{1:n^2} such that the
32## row sums, column sums, and diagonal sums are all equal to the same value.
33##
34## Note: @var{n} must be a scalar greater than or equal to 3.  If you supply
35## @var{n} less than 3, magic returns either a nonmagic square, or else the
36## degenerate magic squares 1 and [].
37## @end deftypefn
38
39function A = magic (n)
40
41  if (nargin != 1)
42    print_usage ();
43  endif
44
45  n = fix (n);
46  if (n < 0)
47    error ("magic: N must be non-negative");
48  elseif (n < 1)
49    A = [];
50  elseif (mod (n, 2) == 1)
51
52    shift = floor ((0:n*n-1)/n);
53    c = mod ([1:n*n] - shift + (n-3)/2, n);
54    r = mod ([n*n:-1:1] + 2*shift, n);
55    A(c*n+r+1) = 1:n*n;
56    A = reshape (A, n, n);
57
58  elseif (mod (n, 4) == 0)
59
60    A = reshape (1:n*n, n, n)';
61    I = [1:4:n, 4:4:n];
62    J = fliplr (I);
63    A(I,I) = A(J,J);
64    I = [2:4:n, 3:4:n];
65    J = fliplr (I);
66    A(I,I) = A(J,J);
67
68  elseif (mod (n, 4) == 2)
69
70    m = n/2;
71    A = magic (m);
72    A = [A, A+2*m*m; A+3*m*m, A+m*m];
73    k = (m-1)/2;
74    if (k > 1)
75      I = 1:m;
76      J = [2:k, n-k+2:n];
77      A([I,I+m],J) = A([I+m,I],J);
78    endif
79    I = [1:k, k+2:m];
80    A([I,I+m],1) = A([I+m,I],1);
81    I = k + 1;
82    A([I,I+m],I) = A([I+m,I],I);
83
84  endif
85
86endfunction
87
88
89%!test
90%! for i = 3:30
91%!   A = magic (i);
92%!   assert (norm(diff([sum(diag(A)),sum(diag(flipud(A))),sum(A),sum(A')])),0);
93%! endfor
94
95## Not a magic square but we must return something (bug #46672).
96## While one day we may change the actual return of magic (2),
97## this properties still must be true.
98%!test <*46672>
99%! m = magic (2);
100%! assert (size (m), [2 2]);
101%! assert (m, [4 3; 1 2]);
102
103%!assert (isempty (magic (0)))
104%!assert (magic (1), 1)
105%!assert (magic (1.5), 1)
106
107## Test input validation
108%!error magic ()
109%!error magic (1, 2)
110%!error <N must be non-negative> magic (-5)
111