A389551 - OEIS

%I #33 Dec 03 2025 17:54:57

%S 432,648,864,1728,1944,2000,2592,3456,3888,4000,5000,5184,5488,5832,

%T 6912,8000,10125,10368,10976,11664,13824,15552,16000,16875,17496,

%U 19208,20000,20736,21296,21952,23328,25000,27648,27783,30375,31104,32000,34992,35152,40000

%N Cubefull numbers that are not powers of squarefree numbers.

%C Excludes powers k^m of squarefree numbers k, i.e., A072777.

%C Intersection of A359280 and A036966, where A359280 = A001694 \ A072777.

%C Disjoint union of A383394, A388549, and A388293 = A359280 \ A389558.

%H Michael De Vlieger, <a href="/A389551/b389551.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 From _Amiram Eldar_, Nov 23 2025: (Start)

%F Equals (A036966 \ A390540) \ {1} .

%F Sum_{n>=1} 1/a(n) = A065483 - A368250 + A082020 - 2 = 0.01096804024453018557... . (End)

%e Let r = A388293, s = A383394, and t = A388549.

%e Table of n, a(n) for select n:

%e  n     a(n)

%e ------------------------------------

%e  1     432 = 2^4 * 3^3        = r(1)

%e  2     648 = 2^3 * 3^4        = r(2)

%e  3     864 = 2^5 * 3^3        = r(3)

%e  4    1728 = 2^6 * 3^3 = 12^3 = t(1)

%e  5    1944 = 2^3 * 3^5        = r(4)

%e  6    2000 = 2^4 * 5^3        = r(5)

%e  7    2592 = 2^5 * 3^4        = r(6)

%e  8    3456 = 2^7 * 3^3        = r(7)

%e  9    3888 = 2^4 * 3^5        = r(8)

%e 12    5184 = 2^6 * 5^4 = 72^2 = s(1)

%e 17   10125 = 3^4 * 5^3        = r(13)

%e 46   54000 = 2^4 * 3^3 * 5^3  = r(37)

%t With[{nn = 40000}, Union@ Flatten@ Table[If[! PrimePowerQ[#], If[CountDistinct[FactorInteger[#][[;;,-1]] ]>1, #, Nothing], Nothing] &[a^5 * b^4 * c^3], {c, Surd[nn, 3]}, {b,Surd[nn/(c^3), 4]}, {a, Surd[nn/(b^4 * c^3), 5] } ] ]

%o (Python)

%o from math import isqrt, gcd

%o from sympy import integer_nthroot, mobius, factorint

%o from oeis_sequences.OEISsequences import bisection, squarefreepi

%o def A389551(n):

%o     def g(x):

%o         c, l = squarefreepi(integer_nthroot(x,3)[0])+sum(mobius(k)*(integer_nthroot(x, k)[0]-1) for k in range(2, x.bit_length()))-1, 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     def f(x):

%o         c = n+1+x+sum(squarefreepi(integer_nthroot(x,k)[0])-1 for k in range(3, x.bit_length()))

%o         for w in range(1,integer_nthroot(x,5)[0]+1):

%o             if all(d<=1 for d in factorint(w).values()):

%o                 for y in range(1,integer_nthroot(z:=x//w**5,4)[0]+1):

%o                     if gcd(w,y)==1 and all(d<=1 for d in factorint(y).values()):

%o                         c -= integer_nthroot(z//y**4,3)[0]

%o         return c

%o     return bisection(f,n,n) # _Chai Wah Wu_, Dec 02 2025

%Y Cf. A001597, A001694, A036966, A072777, A126706, A131605, A286708, A359280, A383394, A388293, A388549, A389558, A390540.

%Y Cf. A065483, A082020, A368250.

%K nonn,easy

%O 1,1

%A _Michael De Vlieger_, Nov 17 2025