%I #16 Sep 11 2024 00:35:02
%S 1,16,64,81,256,625,729,1024,1296,2401,4096,5184,6561,10000,11664,
%T 14641,15625,16384,20736,28561,38416,40000,46656,50625,59049,65536,
%U 82944,83521,104976,117649,130321,153664,160000,186624,194481,234256,250000,262144,279841,331776
%N Squares of powerful numbers.
%C First differs from A340588 at n = 12.
%C 4-full (or 3-full) squares.
%C Numbers whose exponents in their prime factorization are all even numbers >= 4.
%C This sequence is closed under multiplication.
%C The sequence {A000290(n)*A078615(A000290(n)), n>=1} is a permutation of this sequence, and the sequence {a(n)/A078615(a(n)), n>=1} is a permutation of {A000290(n), n>=1}.
%C The sequence {A335988(n)*A007947(A335988(n)), n>=1} is a permutation of this sequence, and the sequence {a(n)/A007947(a(n)), n>=1} is a permutation of A335988.
%H Amiram Eldar, <a href="/A374291/b374291.txt">Table of n, a(n) for n = 1..10000</a>
%H <a href="/index/Pow#powerful">Index entries for sequences related to powerful numbers</a>.
%F a(n) = A000290(A001694(n)) = A001694(n)^2.
%F Sum_{n>=1} 1/a(n) = Product_{p prime} (1 + 1/(p^2*(p^2-1))) = zeta(4)*zeta(6)/zeta(12) = 15015/(1382*Pi^2) = 1.10082313486953808844... .
%F Sum_{n>=1} 1/a(n)^s = Product_{p prime} (1 + 1/(p^(2*s)*(p^(2*s)-1))) = zeta(4*s)*zeta(6*s)/zeta(12*s), for s > 1/4.
%t powQ[n_] := n==1 || AllTrue[FactorInteger[n][[;; , 2]], # > 1 &]; Select[Range[600], powQ]^2
%o (PARI) is(k) = issquare(k) && ispowerful(sqrtint(k));
%o (Python)
%o from math import isqrt
%o from sympy import mobius, integer_nthroot
%o def A374291(n):
%o def squarefreepi(n):
%o return int(sum(mobius(k)*(n//k**2) for k in range(1, isqrt(n)+1)))
%o def bisection(f, kmin=0, kmax=1):
%o while f(kmax) > kmax: kmax <<= 1
%o while kmax-kmin > 1:
%o kmid = kmax+kmin>>1
%o if f(kmid) <= kmid:
%o kmax = kmid
%o else:
%o kmin = kmid
%o return kmax
%o def f(x):
%o c, l = n+x, 0
%o j = isqrt(x)
%o while j>1:
%o k2 = integer_nthroot(x//j**2, 3)[0]+1
%o w = squarefreepi(k2-1)
%o c -= j*(w-l)
%o l, j = w, isqrt(x//k2**3)
%o c -= squarefreepi(integer_nthroot(x, 3)[0])-l
%o return c
%o return bisection(f,n,n)**2 # _Chai Wah Wu_, Sep 10 2024
%Y Intersection of A000290 and A036967 (or A036966).
%Y Intersection of A000290 and A337050.
%Y Subsequence of A322449.
%Y Cf. A000290, A001694, A007947, A078615, A335988, A340588.
%K nonn,easy
%O 1,2
%A _Amiram Eldar_, Jul 02 2024