cp's OEIS Frontend

A320966 Powerful numbers A001694 divisible by a cube > 1.

%I A320966 #35 Sep 15 2024 22:01:30
%S A320966 8,16,27,32,64,72,81,108,125,128,144,200,216,243,256,288,324,343,392,
%T A320966 400,432,500,512,576,625,648,675,729,784,800,864,968,972,1000,1024,
%U A320966 1125,1152,1296,1323,1331,1352,1372,1568,1600,1728,1800,1936,1944,2000,2025,2048,2187,2197,2304,2312,2401,2500
%N A320966 Powerful numbers A001694 divisible by a cube > 1.
%C A320966 Powerful numbers that are not squares of squarefree numbers. - _Amiram Eldar_, Jun 25 2022
%H A320966 Hugo Pfoertner, <a href="/A320966/b320966.txt">Table of n, a(n) for n = 1..10000</a>
%F A320966 Sum_{n>=1} 1/a(n) = zeta(2)*zeta(3)/zeta(6) - 15/Pi^2 = 0.4237786821... . - _Amiram Eldar_, Jun 25 2022
%t A320966 Select[Range[2500], (m = MinMax[FactorInteger[#][[;; , 2]]])[[1]] > 1 && m[[2]] > 2 &] (* _Amiram Eldar_, Jun 25 2022 *)
%o A320966 (PARI) isA001694(n)=n=factor(n)[, 2]; for(i=1, #n, if(n[i]==1, return(0))); 1 \\ from _Charles R Greathouse IV_
%o A320966 isA046099(n)=n=factor(n)[, 2]; for(i=1, #n, if(n[i]>2, return(1)));0
%o A320966 for (k=1,2500,if(isA001694(k)&&isA046099(k),print1(k,", ")))
%o A320966 (Python)
%o A320966 from math import isqrt
%o A320966 from sympy import mobius, integer_nthroot
%o A320966 def A320966(n):
%o A320966     def squarefreepi(n): return int(sum(mobius(k)*(n//k**2) for k in range(1, isqrt(n)+1)))
%o A320966     def bisection(f,kmin=0,kmax=1):
%o A320966         while f(kmax) > kmax: kmax <<= 1
%o A320966         while kmax-kmin > 1:
%o A320966             kmid = kmax+kmin>>1
%o A320966             if f(kmid) <= kmid:
%o A320966                 kmax = kmid
%o A320966             else:
%o A320966                 kmin = kmid
%o A320966         return kmax
%o A320966     def f(x):
%o A320966         c, l = n+x+squarefreepi(isqrt(x))-squarefreepi(integer_nthroot(x,3)[0]), 0
%o A320966         j = isqrt(x)
%o A320966         while j>1:
%o A320966             k2 = integer_nthroot(x//j**2,3)[0]+1
%o A320966             w = squarefreepi(k2-1)
%o A320966             c -= j*(w-l)
%o A320966             l, j = w, isqrt(x//k2**3)
%o A320966         return c+l
%o A320966     return bisection(f,n,n) # _Chai Wah Wu_, Sep 15 2024
%Y A320966 Cf. A000578, A320965.
%Y A320966 Intersection of A001694 and A046099.
%Y A320966 A001694 \ A062503.
%K A320966 nonn
%O A320966 1,1
%A A320966 _Hugo Pfoertner_, Oct 25 2018