Method for finding 3x3 nearly-magic squares of squares with close sums
by Lee Morgenstern, January 2010.
(sent to Frank Rubin in 2009)

I looked up a suitable formula for three Pythagorean Triangles having equal areas. I found one in Dickson's History of the Theory of Numbers, Chapter IV, Rational Right Triangles, near the end of the section, "Right Triangles of Equal Area". It is Hillyer's formula.

The generators of the three triangles are
p = K^2 + KL + L^2, q = K^2 - L^2
r = K^2 + KL + L^2, s = 2KL + L^2
v = K^2 + KL + L^2, u = 2KL + K^2

Here is how to use Hillyer's formula to form a 3x3 nearly-magic square of squares.

A^2 B^2 C^2
D^2 E^2 F^2
G^2 H^2 I^2

A,E,I will be the magic diagonal.
C,E,G
will be the non-magic diagonal.

A = u^2 + v^2
B = |p^2 - q^2 - 2pq|
C = r^2 - s^2 + 2rs
D = p^2 - q^2 + 2pq
E = r^2 + s^2
F = |u^2 - v^2 - 2uv|
G = |r^2 - s^2 - 2rs|
H = u^2 - v^2 + 2uv
I = p^2 + q^2

For example, suppose L = 7, K = 10.

p = r = v = 100 + 70 + 49 = 219
q = K^2 - L^2 = 100 - 49  =  51
s = 2KL + L^2 = 140 + 49  = 189
u = 2KL + K^2 = 140 + 100 = 240

A = 240^2 + 219^2               = 105561
B = |219^2 - 51^2 - 2x219x51|   =  23022
C = 219^2 - 189^2 + 2x219x189   =  95022
D = 219^2 - 51^2 + 2x219x51     =  67698
E = 219^2 + 189^2               =  83682
F = |240^2 - 219^2 - 2x240x219| =  95481
G = |219^2 - 189^2 - 2x219x189| =  70542
H = 240^2 - 219^2 + 2x240x219   = 114759
I = 219^2 + 51^2                =  50562

The magic sum is about 2.07023 x 10^10 and the non-magic diagonal sum = 3E^2 is about 2.10080 x 10^10 for a ratio of about 1.015.

In order for the ratio to be exactly 1, we need to have

A^2 + I^2 = 2E^2.

Plugging in the expressions for K and L yields the polynomial

L^8 + 12L^7K + 34L^6K^2 + 44L^5K^3 + 20L^4K^4 - 8L^3K^5 - 12L^2K^6 - 8LK^7 - 2K^8 = 0

If we divide by K^8 and set x = L/K, we get

f(x) = x^8 + 12x^7 + 34x^6 + 44x^5 + 20x^4 - 8x^3 - 12x^2 - 8x - 2 = 0

The first derivative of this is

f'(x) = 8x^7 + 84x^6 + 204x^5 + 220x^4 + 80x^3 - 24x^2 - 24x - 8

Using Newton's iteration for root finding, if x[n] is a good guess at a root, then a better guess will be

x[n+1] = x[n] - f(x) / f'(x).

If you graph f(x) between 0 and 1, there is one root at about 0.7. Using a few Newton iterations, you can improve that to 0.68793. If you convert x back to L/K, which can be done simply by setting L = 68793, K = 100000, you get a nearly-magic square of integer squares with a diagonal sum ratio of 1.000004.

More Newton interations will get a better approximation and the diagonal sum ratio will get closer to 1 to any desired accuracy.

Note that Hillyer's formula can never produce a fully-magic square of squares.