cp's OEIS Frontend

A378768 Squares of powerful numbers that are not prime powers.

%I A378768 #16 Dec 10 2024 09:58:47
%S A378768 1296,5184,10000,11664,20736,38416,40000,46656,50625,82944,104976,
%T A378768 153664,160000,186624,194481,234256,250000,331776,419904,455625,
%U A378768 456976,614656,640000,746496,810000,937024,944784,1000000,1185921,1265625,1327104,1336336,1500625,1679616
%N A378768 Squares of powerful numbers that are not prime powers.
%C A378768 Contained in A286708, which is a proper subset of A126706.
%H A378768 Michael De Vlieger, <a href="/A378768/b378768.txt">Table of n, a(n) for n = 1..10000</a>
%F A378768 a(n) = A286708(n)^2.
%F A378768 Intersection of A000290 and A286708.
%F A378768 Intersection of A000290 and A372695.
%F A378768 Sum_{n>=1} 1/a(n) = zeta(4)*zeta(6)/zeta(12) - Sum_{p prime} (1/(p^4-p^2)) - 1 = 0.0013772572536044025109... . - _Amiram Eldar_, Dec 10 2024
%t A378768 With[{nn = 2000}, Select[Rest@ Union[Flatten@ Table[a^2*b^3, {b, Surd[nn, 3]}, {a, Sqrt[nn/b^3]}] ], Not@*PrimePowerQ]^2]
%o A378768 (Python)
%o A378768 from math import isqrt
%o A378768 from sympy import integer_nthroot, primepi, mobius
%o A378768 def A378768(n):
%o A378768     def squarefreepi(n): return int(sum(mobius(k)*(n//k**2) for k in range(1, isqrt(n)+1)))
%o A378768     def bisection(f, kmin=0, kmax=1):
%o A378768         while f(kmax) > kmax: kmax <<= 1
%o A378768         while kmax-kmin > 1:
%o A378768             kmid = kmax+kmin>>1
%o A378768             if f(kmid) <= kmid:
%o A378768                 kmax = kmid
%o A378768             else:
%o A378768                 kmin = kmid
%o A378768         return kmax
%o A378768     def f(x):
%o A378768         c, l = n+x, 0
%o A378768         j = isqrt(x)
%o A378768         while j>1:
%o A378768             k2 = integer_nthroot(x//j**2, 3)[0]+1
%o A378768             w = squarefreepi(k2-1)
%o A378768             c -= j*(w-l)
%o A378768             l, j = w, isqrt(x//k2**3)
%o A378768         c -= squarefreepi(integer_nthroot(x, 3)[0])-l
%o A378768         return c+1+sum(primepi(integer_nthroot(x, k)[0]) for k in range(2, x.bit_length()))
%o A378768     return bisection(f, n, n)**2 # _Chai Wah Wu_, Dec 08 2024
%Y A378768 Cf. A000290, A001694, A126706, A286708, A372695.
%Y A378768 Cf. A013662, A013664, A013670.
%K A378768 nonn,easy
%O A378768 1,1
%A A378768 _Michael De Vlieger_, Dec 06 2024