login
Achilles numbers that are not cubefull.
8

%I #22 Dec 05 2025 23:38:26

%S 72,108,200,288,392,500,675,800,968,972,1125,1152,1323,1352,1372,1568,

%T 1800,2312,2700,2888,3087,3200,3267,3528,3872,4232,4500,4563,4608,

%U 5292,5324,5400,5408,6075,6125,6272,6728,7200,7688,7803,8575,8712,8748,8788,9000,9248

%N Achilles numbers that are not cubefull.

%C Intersection of A052486 and A362147 = A052486 \ A036966.

%H Michael De Vlieger, <a href="/A390539/b390539.txt">Table of n, a(n) for n = 1..10000</a>

%F From _Amiram Eldar_, Nov 18 2025: (Start)

%F Equals A325240 \ (A000290 \ A374291).

%F Sum_{n>=1} 1/a(n) = A082695 - A065483 - A013661 + zeta(4)*zeta(6)/zeta(12) = 0.059701351267723610429... . (End)

%e Let s = A052486.

%e a(1) = s(1) = 72 = 2^3 * 3*2 is the first term, since among prime power factor exponents, there exists an exponent m < 3.

%e a(2) = s(2) = 108 = 2^2 * 3*3, since among prime power factor exponents, there exists an exponent m < 3.

%e s(6) = 432 = 2^4 * 3^3 is not a term, since 432 is the smallest Achilles number with all prime power factor exponents greater than 2.

%p filter:= proc(n) local F;

%p F:= ifactors(n)[2][..,2];

%p min(F) = 2 and igcd(F) = 1

%p end proc:

%p select(filter, [$1..10000]); # _Robert Israel_, Nov 16 2025

%t With[{nn = 10^4}, Select[Rest@ Union@ Flatten@ Table[a^2*b^3, {b, Surd[nn, 3]}, {a, Sqrt[nn/b^3]}], And[Min[#] < 3, GCD @@ # == 1] &[FactorInteger[#][[;; , -1]] ] &]]

%o (Python)

%o from math import isqrt, gcd

%o from sympy import integer_nthroot, factorint

%o from oeis_sequences.OEISsequences import bisection, squarefreepi

%o def A390539(n):

%o def f(x):

%o c, l = n+x, 0

%o j = z = 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 += l+z-squarefreepi(integer_nthroot(x,3)[0])-sum(isqrt(z//k**3) for k in range(1, integer_nthroot(z, 3)[0]+1) if all(d<=1 for d in factorint(k).values()))

%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 04 2025

%Y Cf. A001694, A013929, A024619, A036966, A052486, A126706, A286708, A362147, A388293.

%Y Cf. A000290, A013661, A065483, A082695, A325240, A374291.

%K nonn,easy

%O 1,1

%A _Michael De Vlieger_, Nov 16 2025