%I #45 Sep 10 2024 19:40:34

%S 36,92,88,104,40,68,148,27,125,64,104,4,153,27,171,29,20,196,232,144,

%T 56,312,280,108,188,199,113,67,189,72,344,16,112,232,268,63,45,392,

%U 292,32,76,8,80,587,50,147,456,184,288,488,115,772,137,36,40,212,248

%N Differences of consecutive Achilles numbers.

%C 29 is the first prime in this sequence, and it equals 1352 - 1323. Clearly, if the difference is prime, then these two Achilles numbers must be relatively prime, so primes appear in this sequence rarely. However, are there infinitely many n such that a(n) is prime?

%C The number 1 can also appear in this sequence, because it equals 5425069448 - 5425069447 = (2^3 * 26041^2) - (7^3 * 41^2 * 97^2). Does every natural number appear in this sequence? If so, do they appear infinitely often?

%H Amiram Eldar, <a href="/A247246/b247246.txt">Table of n, a(n) for n = 1..10000</a>

%H Carlos Rivera, <a href="http://www.primepuzzles.net/problems/prob_053.htm">Problem 53</a>, The Prime Puzzles and Problems Connection.

%H Eric Weisstein's World of Mathematics, <a href="http://mathworld.wolfram.com/AchillesNumber.html">Achilles number</a>.

%H Wikipedia, <a href="http://en.wikipedia.org/wiki/Achilles_number">Achilles number</a>.

%F a(n) = A052486(n+1) - A052486(n).

%p f:= proc(n) local E; E:= map(t->t[2], ifactors(n)[2]); min(E)>1 and igcd(op(E))=1 end proc:

%p Achilles:= select(f, [$1..10^5]):

%p seq(Achilles[i+1]-Achilles[i],i=1..nops(Achilles)-1); # _Robert Israel_, Dec 13 2014

%t achillesQ[n_] := With[{ee = FactorInteger[n][[All, 2]]}, Min[ee] > 1 && GCD @@ ee == 1];

%t Select[Range[10^4], achillesQ] // Differences (* _Jean-François Alcover_, Sep 26 2020 *)

%o (PARI) isA052486(n) = { n>1 & vecmin(factor(n)[, 2])>1 & !ispower(n); }

%o lista(nn) = {v = select(n->isA052486(n), vector(nn, i, i)); vector(#v-1, n, v[n+1] - v[n]);} \\ _Michel Marcus_, Nov 29 2014

%o (Python)

%o from math import isqrt

%o from sympy import mobius, integer_nthroot

%o def A247246(n):

%o def squarefreepi(n):

%o 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+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 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 return c

%o return -(a:=bisection(f,n,n))+bisection(lambda x:f(x)+1,a,a) # _Chai Wah Wu_, Sep 10 2024

%Y Cf. A052486, A076446, A053289.

%K nonn,easy

%O 1,1

%A _Eric Chen_, Nov 28 2014