% Enumerate magic squares of order 4.
% Count number with rank = 3 and rank = 4.
% See Bill Press, http://www.nr.com/whp/notes/DurerSquare.pdf.
% Results: rank3: 5120 = 8*640, rank4: 1920 = 8*240.
% CBM, 11/9/2012. 

S = 34;  % Magic sum for order 4.
tic
R = zeros(1,4);
k = 0;
parfor a = 1:16
   for b = 1:16
      if b==a, continue, end
      for c = 1:16
         if c==b || c==a, continue, end
         for e = 1:16
            if e==c || e==b || e==a, continue, end
            for f = 1:16
               if f==e || f==c || f==b || f==a, continue, end
               for g = 1:16
                  if g==f || g==e ||g==c || g==b || g==a, continue, end
                  for i = 1:16
                     if i==g || i==f || i==e || i==c || i==b || i==a, continue, end
                     A = [a        b                      c                  S-a-b-c
                          e        f                      g                  S-e-f-g
                          i        2*a+b+c+e-g+i-S        2*S-2*a-b-c-e-f-i  f+g-i
                          S-a-e-i  2*S-2*a-2*b-c-e-f+g-i  2*a+b+e+f-g+i-S    a+b+c+e+i-S];
                     if isequal(sort(A(:)),(1:16)')
                        % A is a magic square
                        clc
                        k = k+1  % Should count to 7040
                        A
                        r = rank(A);
                        R(r) = R(r)+1
                        pause(eps)
                     end
                  end
               end
            end
         end
      end
   end
end
toc  % 415 seconds on an Intel 2.67 GHz i7 laptopg/
