function [R1, R2, R3, R4, k] = enumerate4_p
% [R1, R2, R3, R4, k] = enumerate4_p
% Three difference ways to parallelize enumerate4.m.
% Enumerate magic squares of order 4.
% Count number with rank = 3 and rank = 4.
% See Cleve's Corner Blog, 11/26/2012.

S = 34;  % Magic sum for order 4.
R1 = zeros(4, 1);
R2 = [];
R3 = 0;
R4 = 0;
k = 0;

tic
parfor ab = 0:255
   a = floor(ab/16) + 1;
   b = rem(ab, 16) + 1;
      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
                        r = rank(A);
                        k = k + 1;
                        x = zeros(4, 1);
                        x(r) = 1;
                        R1 = R1 + x;
                        R2 = [R2 r];
                        if r==3
                           R3 = R3+1
                        else
                           R4 = R4+1
                        end
                     end
                  end
               end
            end
         end
      end
 end
 toc
