|
%I #9 Mar 08 2024 08:07:02
%S 1,4,8,9,16,25,32,36,49,64,72,81,100,121,125,128,144,169,196,200,225,
%T 256,288,289,324,361,392,400,441,484,500,512,529,576,625,648,676,729,
%U 784,800,841,900,961,968,1000,1024,1089,1125,1152,1156,1225,1296,1352,1369
%N Powerful numbers that are the sum of 2 squares.
%C Each term can be decomposed in a unique way as 2^m * i * j^2 where m >= 2, i is a powerful number whose prime factors are all of the form 4*k + 1 (A369563), and j is a number whose prime factors are all of the form 4*k + 3 (A004614).
%H Amiram Eldar, <a href="/A371010/b371010.txt">Table of n, a(n) for n = 1..10000</a>
%H Rafael Jakimczuk, <a href="http://dx.doi.org/10.13140/RG.2.2.27745.48487">Generalizations of Mertens's Formula and k-Free and s-Full Numbers with Prime Divisors in Arithmetic Progression</a>, ResearchGate, 2024.
%F The number of terms that do not exceed x is ~ c * sqrt(x), where c = (6/Pi^2) * (1 + 1/(3*(sqrt(2)-1))) * Product_{primes p == 1 (mod 4)} (1 + 1/((sqrt(p)-1)*(p+1))) * Product_{primes p == 3 (mod 4)} (1 + 1/(p^2-1)) = 1.58769... (Jakimczuk, 2024, Theorem 4.7, p. 50).
%F Sum_{n>=1} 1/a(n) = (3/2) * Product_{primes p == 1 (mod 4)} (1 + 1/(p*(p-1))) * Product_{primes p == 3 (mod 4)} (1 + 1/(p^2-1)) = (3*Pi^2/16) * A334424 = 1.86676402705119927669... .
%t Select[Range[1500], SquaresR[2, #] > 0 && (# == 1 || Min[FactorInteger[#][[;; , 2]]] > 1) &]
%o (PARI) is(n) = {my(f=factor(n)); for(i=1, #f~, if(f[i, 2] == 1 || (f[i, 2]%2 && f[i, 1]%4 == 3), return(0))); 1;}
%Y Intersection of A001481 and A001694.
%Y A371011 is a subsequence.
%Y Cf. A004614, A020893, A369563.
%K nonn,easy
%O 1,2
%A _Amiram Eldar_, Mar 08 2024
|
