Tricky Triads

An integer $n$ is expressible as the sum of two squares if and only if it is of the form $n = A B^2,$ where $A$ is divisible only by primes not congruent to $3 \bmod 4$ and $B$ is divisible only by primes congruent to $3 \bmod 4.$

Writing $A = 2^a p_1^{a_1} \cdots p_k^{a_k},$ the number of representations of $n$ as a sum of two squares is $r_2(n) = 4(a_1+1)(a_2+1)(\cdots)(a_k+1).$ This formula (and the previous paragraph) comes from facts about Gaussian integers .

If $n$ is a square and we want to write it as the sum of two positive squares in two essentially different ways, we'll want $r_2(n) \ge 20.$ This is because there are four ways to write $n = k^2 + 0^2 = (-k)^2 + 0^2 = 0^2+k^2 = 0^2 +(-k)^2,$ none of which help us, and then an actual representation of $n=a^2+b^2$ as a sum of two positive squares will give rise to eight different representations counted in $r_2(n)$ (by switching $a$ and $b$ and/or applying negative signs to one or both of them.)

So we want $r_2(n) \ge 20$ and all the $a_i$ even. Both $a_1=4$ and $a_1=a_2=2$ will work, and any other choices will lead to larger numbers. The smallest example of the former is $5^4 = 625,$ and the smallest example of the latter is $(5^2 \cdot 13^2 = 4225.$ So the answer is $\fbox{625}.$ The two representations are $625 = 15^2+20^2 = 7^2 + 24^2.$

Here is a little python code i came up with, probably it's not the most optimal solution but it works just fine ;)

def perfectSquare():

nums=[]

for i in range (1,10000):
    for j in range (1,100):
        for k in range (1,100):
            if i==j*j+k*k and i!=1 and \
            not j in nums and not k in nums and \
            (math.sqrt(i) % 1 == 0): 
                    nums.append(j)
                    nums.append(k)
    if len(nums) == 4 :
        print nums,
        print nums[0]*nums[0]+nums[1]*nums[1]
    del nums[:]