cp's OEIS Frontend

A052486 Achilles numbers - powerful but imperfect: if n = Product(p_i^e_i) then all e_i > 1 (i.e., powerful), but the highest common factor of the e_i is 1, i.e., not a perfect power.

%I A052486 #97 Feb 16 2025 08:32:42
%S A052486 72,108,200,288,392,432,500,648,675,800,864,968,972,1125,1152,1323,
%T A052486 1352,1372,1568,1800,1944,2000,2312,2592,2700,2888,3087,3200,3267,
%U A052486 3456,3528,3872,3888,4000,4232,4500,4563,4608,5000,5292,5324,5400,5408,5488,6075
%N A052486 Achilles numbers - powerful but imperfect: if n = Product(p_i^e_i) then all e_i > 1 (i.e., powerful), but the highest common factor of the e_i is 1, i.e., not a perfect power.
%C A052486 Number of terms < 10^n: 0, 1, 13, 60, 252, 916, 3158, 10553, 34561, 111891, 359340, 1148195, 3656246, 11616582, 36851965, ..., A118896(n) - A070428(n). - _Robert G. Wilson v_, Aug 11 2014
%C A052486 a(n) = (s(n))^2 * f(n), s(n) > 1, f(n) > 1, where s(n) is not a power of f(n), and f(n) is squarefree and gcd(s(n), f(n)) = f(n). - _Daniel Forgues_, Aug 11 2015
%H A052486 Amiram Eldar, <a href="/A052486/b052486.txt">Table of n, a(n) for n = 1..10000</a> (terms 1..1000 from T. D. Noe)
%H A052486 Project Euler, <a href="https://projecteuler.net/problem=302">Problem 302: Strong Achilles Numbers</a>.
%H A052486 Robert Israel, <a href="/A052486/a052486.png">log-log plot of a(n)</a>.
%H A052486 Eric Weisstein's World of Mathematics, <a href="https://mathworld.wolfram.com/AchillesNumber.html">Achilles Number</a>.
%H A052486 OEIS Wiki, <a href="https://oeis.org/wiki/Achilles_numbers">Achilles numbers</a>.
%F A052486 a(n) = O(n^2). - _Daniel Forgues_, Aug 11 2015
%F A052486 a(n) = O(n^2 / log log n). - _Daniel Forgues_, Aug 12 2015
%F A052486 Sum_{n>=1} 1/a(n) = zeta(2)*zeta(3)/zeta(6) - Sum_{k>=2} mu(k)*(1-zeta(k)) - 1 = A082695 - A072102 - 1 = 0.06913206841581433836...  - _Amiram Eldar_, Oct 14 2020
%e A052486 a(3)=200 because 200=2^3*5^2, both 3 and 2 are greater than 1, and the highest common factor of 3 and 2 is 1.
%e A052486 Factorizations of a(1) to a(20):
%e A052486     72 = 2^3  3^2,  108 = 2^2 3^3,  200 = 2^3 5^2,  288 = 2^5  3^2,
%e A052486    392 = 2^3  7^2,  432 = 2^4 3^3,  500 = 2^2 5^3,  648 = 2^3  3^4,
%e A052486    675 = 3^3  5^2,  800 = 2^5 5^2,  864 = 2^5 3^3,  968 = 2^3 11^2,
%e A052486    972 = 2^2  3^5, 1125 = 3^2 5^3, 1152 = 2^7 3^2, 1323 = 3^3  7^2,
%e A052486   1352 = 2^3 13^2, 1372 = 2^2 7^3, 1568 = 2^5 7^2, 1800 = 2^3  3^2 5^2.
%e A052486 Examples for a(n) = (s(n))^2 * f(n): (see above comment)
%e A052486 s(n) = 6,  6, 10, 12, 14, 12, 10, 18, 15, 20, 12, 22, 18, 15, 24, 21,
%e A052486 f(n) = 2,  3,  2,  2,  2,  3,  5,  2,  3,  2,  6,  2,  3,  5,  2,  3,
%p A052486 filter:= proc(n) local E; E:= map(t->t[2], ifactors(n)[2]); min(E)>1 and igcd(op(E))=1 end proc:
%p A052486 select(filter,[$1..10000]); # _Robert Israel_, Aug 11 2014
%t A052486 achillesQ[n_] := Block[{ls = Last /@ FactorInteger@n}, Min@ ls > 1 == GCD @@ ls]; Select[ Range@ 5500, achillesQ@# &] (* _Robert G. Wilson v_, Jun 10 2010 *)
%o A052486 (PARI) is(n)=my(f=factor(n)[,2]); n>9 && vecmin(f)>1 && gcd(f)==1 \\ _Charles R Greathouse IV_, Sep 18 2015, replacing code by _M. F. Hasler_, Sep 23 2010
%o A052486 (Python)
%o A052486 from math import gcd
%o A052486 from itertools import count, islice
%o A052486 from sympy import factorint
%o A052486 def A052486_gen(startvalue=1): # generator of terms >= startvalue
%o A052486     return (n for n in count(max(startvalue,1)) if (lambda x: all(e > 1 for e in x) and gcd(*x) == 1)(factorint(n).values()))
%o A052486 A052486_list = list(islice(A052486_gen(),20)) # _Chai Wah Wu_, Feb 19 2022
%o A052486 (Python)
%o A052486 from math import isqrt
%o A052486 from sympy import mobius, integer_nthroot
%o A052486 def A052486(n):
%o A052486     def squarefreepi(n): return int(sum(mobius(k)*(n//k**2) for k in range(1, isqrt(n)+1)))
%o A052486     def bisection(f,kmin=0,kmax=1):
%o A052486         while f(kmax) > kmax: kmax <<= 1
%o A052486         while kmax-kmin > 1:
%o A052486             kmid = kmax+kmin>>1
%o A052486             if f(kmid) <= kmid:
%o A052486                 kmax = kmid
%o A052486             else:
%o A052486                 kmin = kmid
%o A052486         return kmax
%o A052486     def f(x):
%o A052486         c, l = n+x+1, 0
%o A052486         j = isqrt(x)
%o A052486         while j>1:
%o A052486             k2 = integer_nthroot(x//j**2,3)[0]+1
%o A052486             w = squarefreepi(k2-1)
%o A052486             c -= j*(w-l)
%o A052486             l, j = w, isqrt(x//k2**3)
%o A052486         c -= squarefreepi(integer_nthroot(x,3)[0])-l+sum(mobius(k)*(integer_nthroot(x, k)[0]-1) for k in range(2, x.bit_length()))
%o A052486         return c
%o A052486     return bisection(f,n,n) # _Chai Wah Wu_, Sep 10 2024
%Y A052486 Cf. A001597, A001694, A007916, A072102, A082695.
%K A052486 nonn
%O A052486 1,1
%A A052486 _Henry Bottomley_, Mar 16 2000
%E A052486 Example edited by Mac Coombe (mac.coombe(AT)gmail.com), Sep 18 2010
%E A052486 Name edited by _M. F. Hasler_, Jul 17 2019