Mostly to Eric Angelini:
Hello SeqFan and math-fun, a somehow self-describing sequence here, with a question: http://www.cetteadressecomportecinquantesignes.com/PrimePos.htm Best, �ric A.
In the limit, exactly half of the terms are prime. Here's a formula, found empirically, for a(n) for n>=5: Let pi(n) be the number of primes <= n and p(n) be the n'th prime. Then: if n is prime or (n is composite and n+pi(n) is even) then a(n) = p(floor((n+pi(n))/2)); if n is composite and n+pi(n) is odd and n+1 is composite then a(n) = n+1; if n is composite and n+pi(n) is odd and n+1 is prime then a(n) = n+2. Also, for n>=5, n is in the sequence iff either n is prime or n+pi(n) is even. I'm sure this can all be proved by induction on n; I'll leave it to you to check all of the cases if you want to. It follows from this that, for n>=4, the number of primes among a(1), ..., a(n) is exactly floor((n+pi(n))/2. Since pi(n)/n -> 0 as n -> infinity, this is asymptotic to n/2. Here's some Mathematica code that can be used to check this. a[n] is computed from your definition; b[n] is computed from the formula above: $RecursionLimit=Infinity; a[1]=2; inseq[n_]:=False; (* inseq[n] will be True if n is in the sequence *) inseq[2]=True; a[n_]:=a[n]= Module[{m}, a[n-1]; For[m=1, inseq[m] || ((inseq[n] || m==n) && !PrimeQ[m]) || (m<n && !PrimeQ[a[m]]), m++, Null]; inseq[m] = True; Return[m] ]; b[n_]:= If[n<=4, {2,3,5,1}[[n]], If[PrimeQ[n] || EvenQ[n+PrimePi[n]], Prime[Floor[(n+PrimePi[n])/2]], If[PrimeQ[n+1], n+2, n+1] ] ]; a/@Range[1000] == b/@Range[1000] (* Should return True *) countprimes[n_]:=Length[Select[a/@Range[n],PrimeQ[#]&]]; cp[n_]:=Floor[(n+PrimePi[n])/2]; countprimes/@Range[4,1000] == cp/@Range[4,1000] (* Should return True *) Eric, feel free to use any or all of this in your sequence entry. Dean Hickerson dean@math.ucdavis.edu