%I #13 Apr 30 2016 18:01:16
%S 10,17,40,52,58,73,89,90,97,106,145,153,193,202,232,233,241,250,298,
%T 313,337,338,346,360,409,416,424,449,457,468,505,521,522,538,577,586,
%U 634,640,657,673,680,724,730,745,778,801,808,809,810,845,865,873,881,890,953,954,976,986,1000
%N Values of a^2 + b^2 such that sigma(a^2 + b^2) is of the form x^2 + y^2 where a, b, x, y are nonzero integers.
%H Charles R Greathouse IV, <a href="/A268183/b268183.txt">Table of n, a(n) for n = 1..10000</a>
%e 10 is a term because 10 = 1^2 + 3^2 and 10 is divisible by 1, 2, 5, 10 and 1 + 2 + 5 + 10 = 3^2 + 3^2.
%o (PARI) isA000404(n)=for( i=1, #n=factor(n)~%4, n[1, i]==3 && n[2, i]%2 && return); n && ( vecmin(n[1, ])==1 || (n[1, 1]==2 && n[2, 1]%2))
%o lista(nn) = for(n=1, nn, if(isA000404(n) && isA000404(sigma(n)), print1(n, ", ")));
%o (PARI) isA000404(n)= for( i=1, #n=factor(n)~%4, n[1, i]==3 && n[2, i]%2 && return); n && ( vecmin(n[1, ])==1 || (n[1, 1]==2 && n[2, 1]%2))
%o list(lim)=my(v=List(),x2,t); lim\=1; for(x=1,sqrtint(lim-1), x2=x^2; for(y=1,sqrtint(lim-x2), if(isA000404(sigma(t=x2+y^2)), listput(v,t)))); Set(v) \\ _Charles R Greathouse IV_, Apr 30 2016
%Y Cf. A000203, A000404.
%K nonn
%O 1,1
%A _Altug Alkan_, Apr 30 2016