OFFSET
1,1
COMMENTS
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
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
LINKS
Amiram Eldar, Table of n, a(n) for n = 1..10000 (terms 1..1000 from T. D. Noe)
Project Euler, Problem 302: Strong Achilles Numbers.
Robert Israel, log-log plot of a(n).
Eric Weisstein's World of Mathematics, Achilles Number.
OEIS Wiki, Achilles numbers.
FORMULA
a(n) = O(n^2). - Daniel Forgues, Aug 11 2015
a(n) = O(n^2 / log log n). - Daniel Forgues, Aug 12 2015
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
EXAMPLE
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.
Factorizations of a(1) to a(20):
72 = 2^3 3^2, 108 = 2^2 3^3, 200 = 2^3 5^2, 288 = 2^5 3^2,
392 = 2^3 7^2, 432 = 2^4 3^3, 500 = 2^2 5^3, 648 = 2^3 3^4,
675 = 3^3 5^2, 800 = 2^5 5^2, 864 = 2^5 3^3, 968 = 2^3 11^2,
972 = 2^2 3^5, 1125 = 3^2 5^3, 1152 = 2^7 3^2, 1323 = 3^3 7^2,
1352 = 2^3 13^2, 1372 = 2^2 7^3, 1568 = 2^5 7^2, 1800 = 2^3 3^2 5^2.
Examples for a(n) = (s(n))^2 * f(n): (see above comment)
s(n) = 6, 6, 10, 12, 14, 12, 10, 18, 15, 20, 12, 22, 18, 15, 24, 21,
f(n) = 2, 3, 2, 2, 2, 3, 5, 2, 3, 2, 6, 2, 3, 5, 2, 3,
MAPLE
filter:= proc(n) local E; E:= map(t->t[2], ifactors(n)[2]); min(E)>1 and igcd(op(E))=1 end proc:
select(filter, [$1..10000]); # Robert Israel, Aug 11 2014
MATHEMATICA
achillesQ[n_] := Block[{ls = Last /@ FactorInteger@n}, Min@ ls > 1 == GCD @@ ls]; Select[ Range@ 5500, achillesQ@# &] (* Robert G. Wilson v, Jun 10 2010 *)
PROG
(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
(Python)
from math import gcd
from itertools import count, islice
from sympy import factorint
def A052486_gen(startvalue=1): # generator of terms >= startvalue
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()))
(Python)
from math import isqrt
from sympy import mobius, integer_nthroot
def A052486(n):
def squarefreepi(n): return int(sum(mobius(k)*(n//k**2) for k in range(1, isqrt(n)+1)))
def bisection(f, kmin=0, kmax=1):
while f(kmax) > kmax: kmax <<= 1
while kmax-kmin > 1:
kmid = kmax+kmin>>1
if f(kmid) <= kmid:
kmax = kmid
else:
kmin = kmid
return kmax
def f(x):
c, l = n+x+1, 0
j = isqrt(x)
while j>1:
k2 = integer_nthroot(x//j**2, 3)[0]+1
w = squarefreepi(k2-1)
c -= j*(w-l)
l, j = w, isqrt(x//k2**3)
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()))
return c
return bisection(f, n, n) # Chai Wah Wu, Sep 10 2024
CROSSREFS
KEYWORD
nonn
AUTHOR
Henry Bottomley, Mar 16 2000
EXTENSIONS
Example edited by Mac Coombe (mac.coombe(AT)gmail.com), Sep 18 2010
Name edited by M. F. Hasler, Jul 17 2019
STATUS
approved