Handling sums of squares can be mostly automated using some of my mathematical research. For instance consider:
52 + 202 = 17*52
82 + 192 = 17*52
132 + 162 = 17*52
These were found using what I decided to call a Binary Quadratic Diophantine iterator, or BQD Iterator for short, which lets you find in general solutions for x and y with:
x2 + (m-1)y2 = F*mn
When F = x02 + (m-1)y02, my research proves there are non-zero integer solutions for x and y, where n is a count of iterations.
To get my examples used x0 = 1, and y0 = 2, and picked m = 5 so that I had 4 which I could pull into the square, to get a sum of two squares.
At each iteration you get a split point, where you can go positive or negative, which means as you iterate you may generate extra solutions for the same sum.
For my example I had a duplicate in the second iteration which is why there are only 3 distinct solutions instead of 4.
And the 5 being squared in this example is not an accident, as that's the second iteration. The first iteration has 17*5.
So yeah I could have kept going, and would have had a maximum of 6 cases of two squares summing to equal 17*53, but less if there were duplicates.
I have an earlier post which shows how the solutions were found.
And I have additional research using the BQD Iterator showing how you get a desired number of sums of squares. Which could mean for math hobbyists, maybe could be useful in generating magic squares of squares?