login
Squares of powerful numbers.
2

%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