A320966 - OEIS

A320966

Powerful numbers A001694 divisible by a cube > 1.

%I #35 Sep 15 2024 22:01:30

%S 8,16,27,32,64,72,81,108,125,128,144,200,216,243,256,288,324,343,392,

%T 400,432,500,512,576,625,648,675,729,784,800,864,968,972,1000,1024,

%U 1125,1152,1296,1323,1331,1352,1372,1568,1600,1728,1800,1936,1944,2000,2025,2048,2187,2197,2304,2312,2401,2500

%N Powerful numbers A001694 divisible by a cube > 1.

%C Powerful numbers that are not squares of squarefree numbers. - _Amiram Eldar_, Jun 25 2022

%H Hugo Pfoertner, <a href="/A320966/b320966.txt">Table of n, a(n) for n = 1..10000</a>

%F Sum_{n>=1} 1/a(n) = zeta(2)*zeta(3)/zeta(6) - 15/Pi^2 = 0.4237786821... . - _Amiram Eldar_, Jun 25 2022

%t Select[Range[2500], (m = MinMax[FactorInteger[#][[;; , 2]]])[[1]] > 1 && m[[2]] > 2 &] (* _Amiram Eldar_, Jun 25 2022 *)

%o (PARI) isA001694(n)=n=factor(n)[, 2]; for(i=1, #n, if(n[i]==1, return(0))); 1 \\ from _Charles R Greathouse IV_

%o isA046099(n)=n=factor(n)[, 2]; for(i=1, #n, if(n[i]>2, return(1)));0

%o for (k=1,2500,if(isA001694(k)&&isA046099(k),print1(k,", ")))

%o (Python)

%o from math import isqrt

%o from sympy import mobius, integer_nthroot

%o def A320966(n):

%o def squarefreepi(n): 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+squarefreepi(isqrt(x))-squarefreepi(integer_nthroot(x,3)[0]), 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 return c+l

%o return bisection(f,n,n) # _Chai Wah Wu_, Sep 15 2024

%Y Cf. A000578, A320965.

%Y Intersection of A001694 and A046099.

%Y A001694 \ A062503.

%K nonn

%O 1,1

%A _Hugo Pfoertner_, Oct 25 2018