Search method to find 3x3 magic squares with 7 distinct square entries
by Lee Morgenstern, October 2010.

This method is complete, non-redundant, and has the capability of finding solutions having entries of more than 30 digits in a reasonable amount of time. The method uses parametric formulas to directly generate all possible 3x3 magic squares with 5 square entries and then test each solution for two more squares.

Any 7-square solution can be formed from one of the following six configurations of 5-square solutions. Configurations 1,2,3 correspond to 7-square solutions where the central entry is a square. Configurations 4,5,6 correspond to 7-square solutions where the central entry is not necessarily a square. The formulas for the entries have subexpressions in common among the different configurations so that computation can be shared when testing all of them together.

Each of the configurations of 5 squares consists of two 3-square arithmetic progressions having one entry in common.

Any 3-square arithmetic progression is a scaling of a primitive arithmetic progression given by the squares of

Given any other different primitive 3-square arithmetic progression,

we can create two progressions having a common entry by scaling each progression by the common entry of the other progression.

For example, in configuration 1, the central entries are in common, thus we can produce all solutions using

In configuration 4, the smallest entries are in common, thus we can produce all solutions using

Note that, in general, the 5-square solution will be scaled. It is not necessary to reduce a solution to its lowest terms since no scaling of any solution will be a duplicate of any scaling of any other solution using this method.

But for the final output of a solution, it would be nice to have the smallest values. Divide all entries by the GCD of the two common entries. Thus for configuration 1, divide by GCD(Z1, Z2); for configuration 4, divide by GCD(Y1-X1, Y2-X2).

Here are the formulas for the other entries in terms of the entries of the 5-square solutions.

Note that the last three configurations need only test three other entries instead of four because any central square solutions will have been found by the first three configurations.

Configuration 1          Configuration 2          Configuration 3

B^2 = G^2 + I^2 - E^2;   A^2 = (F^2 + H^2) / 2;   C^2 = H^2 + E^2 - A^2;
D^2 = C^2 + I^2 - E^2;   C^2 = (D^2 + H^2) / 2;   D^2 = H^2 + I^2 - A^2;
F^2 = G^2 + A^2 - E^2;   G^2 = (F^2 + B^2) / 2;   F^2 = B^2 + A^2 - I^2;
H^2 = C^2 + A^2 - E^2;   I^2 = (D^2 + B^2) / 2;   G^2 = B^2 + E^2 - I^2;

Configuration 4          Configuration 5          Configuration 6

A^2 = F^2 + I^2 - B^2;   B^2 = D^2 + F^2 - H^2;   I^2 = B^2 + A^2 - F^2;
C^2 = D^2 + G^2 - B^2;   G^2 = C^2 + F^2 - H^2;   D^2 = B^2 + H^2 - F^2;
H^2 = D^2 + F^2 - B^2;   I^2 = D^2 + A^2 - H^2;   C^2 = H^2 + G^2 - F^2;

Here are the formulas in terms of subexpressions of the m,n values.

Note that configurations 3 and 6 need to be tested a second time with the m1,n1 values swapped with the m2,n2 values in order to be complete.

X1 = 2*m1*n1;  Y1 = m1^2-n1^2;  Z1 = m1^2+n1^2;
YX1dif = Y1-X1;  YX1sum = Y1+X1;
Z1sqr = Z1^2;  YX1difsqr = (Y1-X1)^2;  YX1sumsqr = (Y1+X1)^2;

X2 = 2*m2*n2;  Y2 = m2^2-n2^2;  Z2 = m2^2+n2^2;
YX2dif = Y2-X2;  YX2sum = Y2+X2;
Z2sqr = Z2^2;  YX2difsqr = (Y2-X2)^2;  YX2sumsqr = (Y2+X2)^2;

DZ = YX1difsqr * Z2sqr;   DD = YX1difsqr * YX2difsqr;
SZ = YX1sumsqr * Z2sqr;   DS = YX1difsqr * YX2sumsqr;
ZD = Z1sqr * YX2difsqr;   SD = YX1sumsqr * YX2difsqr;
ZS = Z1sqr * YX2sumsqr;   SS = YX1sumsqr * YX2sumsqr;
ZZ = Z1sqr * Z2sqr;

Configuration 1                 Configuration 2
B = Test_Sqrt(ZD + DZ - ZZ);    A = Test_Sqrt((ZS + SZ) / 2);
D = Test_Sqrt(ZS + DZ - ZZ);    C = Test_Sqrt((ZD + SZ) / 2);
F = Test_Sqrt(ZD + SZ - ZZ);    G = Test_Sqrt((ZS + DZ) / 2);
H = Test_Sqrt(ZS + SZ - ZZ);    I = Test_Sqrt((ZD + DZ) / 2);

Configuration 3                 Configuration 3 (swapped)
C = Test_Sqrt(ZS + ZZ - SZ);    C = Test_Sqrt(SZ + ZZ - ZS);
D = Test_Sqrt(ZS + DZ - SZ);    D = Test_Sqrt(SZ + ZD - ZS);
F = Test_Sqrt(ZD + SZ - DZ);    F = Test_Sqrt(DZ + ZS - ZD);
G = Test_Sqrt(ZD + ZZ - DZ);    G = Test_Sqrt(DZ + ZZ - ZD);

Configuration 4                 Configuration 5
A = Test_Sqrt(SD + DZ - DD);    G = Test_Sqrt(DS + SZ - SS);
C = Test_Sqrt(DS + ZD - DD);    I = Test_Sqrt(SD + ZS - SS);
B = Test_Sqrt(SD + DS - SS);    H = Test_Sqrt(DS + SD - DD);

Configuration 6                 Configuration 6 (swapped)
I = Test_Sqrt(DD + SZ - SD);    I = Test_Sqrt(DD + ZS - DS);
C = Test_Sqrt(SS + ZD - SD);    C = Test_Sqrt(SS + DZ - DS);
D = Test_Sqrt(SS + DD - SD);    D = Test_Sqrt(SS + DD - DS);

If two extra squares are found in a configuration, then the starting 5 entries can be computed as follows.

Configuration 1        Configuration 2
E = Z1 * Z2;           E = Z1 * Z2;
A = YX1sum * Z2;       H = YX1sum * Z2;
I = YX1dif * Z2;       B = YX1dif * Z2;
C = YX2sum * Z1;       F = YX2sum * Z1;
G = YX2dif * Z1;       D = YX2dif * Z1;

Configuration 3        Configuration 3 (swapped)
E = Z1 * Z2;           E = Z1 * Z2;
A = YX1sum * Z2;       A = YX2sum * Z1;
I = YX1dif * Z2;       I = YX2dif * Z1;
H = YX2sum * Z1;       H = YX1sum * Z2;
B = YX2dif * Z1;       B = YX1dif * Z2;

Configuration 4        Configuration 5
G = Z1 * YX2dif;       A = Z1 * YX2sum;
B = YX1dif * YX2dif;   F = YX1dif * YX2sum;
F = YX1sum * YX2dif;   H = YX1sum * YX2sum;
I = YX1dif * Z2;       C = YX1sum * Z2;
D = YX1dif * YX2sum;   D = YX1sum * YX2dif;

Configuration 6        Configuration 6 (swapped)
G = Z1 * YX2dif;       G = Z2 * YX1dif;
B = YX1dif * YX2dif;   B = YX2dif * YX1dif;
F = YX1sum * YX2dif;   F = YX2sum * YX1dif;
A = YX1sum * Z2;       A = YX2sum * Z1;
H = Yx1sum * YX2sum;   H = YX2sum * YX1sum;

Here is how to produce a list of appropriate m,n values from which m1,n1 and m2,n2 can be taken.

int putx, getx;
struct { int m,n; } MN_List[MANYMANY];
#define M_LIMIT 2400    // must be at least this big since
                        // all lower values have already
                        // been searched by me
---------------------
Make_MN_List()
{ int m,n,s;

  putx = 0;
  getx = 0;

  MN_List[putx].m = 2;
  MN_List[putx].n = 1;

  ++putx;

  while (getx < putx)
  { m = MN_List[getx].m;
    n = MN_List[getx].n;
    ++getx;

    s = m + n;
    if (s < M_LIMIT)
    { Put_MN(s,m);
      Put_MN(s,n);
    }
  }
}

----------------------

Put_MN(int m, int n)
{ int s;

  if ((n & 1) == 0)
  { MN_List[putx].m = m;
    MN_List[putx].n = n;
    ++putx;
  }
  else
  { s = m + n;
    if (s < M_LIMIT)
     {  MN_List[putx].m = s;
        MN_List[putx].n = m;
        ++putx;
        MN_List[putx].m = s;
        MN_List[putx].n = n;
        ++putx;
     }
  }
}

Test_Sqrt()

 Here is how to efficiently test a number for being a square.

Prepare a bit array containing 13x17x19x23 bits. Initialize the array once at the start of your program so that a bit set indicates that the corresponding number is a quadratic residue modulo 13x17x19x23.

 #define MODULUS 13*17*19*23
Init_QuadRes()
{
  for (b = 0; b < MODULUS; ++b) QuadRes[b] = 0;
  square = 0; square_inc = 1;
  for (b = 0; b < MODULUS; ++b)
  { QuadRes[square] = 1;
    square += square_inc;  square_inc += 2;
    while (square >= MODULUS) square -= MODULUS;
  }
}

Once QuadRes[] has been initialized, it can be used as follows.

// if the odd value sq is not a square, return -1;
// else return the positive square root of sq

bigint Test_Sqrt(bigint sq)
{
  if ((sq & 7) != 1) return -1;
  rem = sq % MODULUS;
  if (QuadRes[rem] == 0) return -1;
  root = sqrt(sq);
  if (root * root != sq) return -1;
  return root;
}

Return to the home page http://www.multimagie.com