Question: how to find the pairs of solutions [x,y] for x^2+y^2=n where n is defined as a product of multiples of 5,13,17,29

I currently have a function quadsum(n) that determines the [x,y] solutions of the above equation for an integer n. :

quadsum:= proc(n::nonnegint)
local
k:= 0, mylist:= table(),
x:= isqrt(iquo(n,2)), y:= x, x2:= x^2, y2:= y^2;
if 2*x2 <> n then x:= x+1; x2:= x2+2*x-1; y:= x; y2:= x2; end if;
while x2 <= n do
y:= isqrt(n-x2); y2:= y^2;
if x2+y2 = n then k:= k+1; mylist[k]:= [x,y] end if;
x:= x+1; x2:= x2+2*x-1;
end do;
convert(mylist, list)
end proc:

How would I alter this so that I get [x,y] for n= (5^a).(13^b).(17^c)(29^d) for non-negative integers a,b,c,d?

Please Wait...