20151009, 01:50  #1  
∂^{2}ω=0
Sep 2002
República de California
3·13^{2}·23 Posts 
Expected number of primes in OEIS A007908
Recently Neil Sloane (curator of the OEIS) sent a message to the NMBRTHRY mailing list re. the sequence A007908:
Quote:
Observations: 1. Based on the trivial observation that only terms ending in 1,3,7,9 have chance of being prime, of the first (say) 100 sequence terms, only 40 can possibly be prime, but in fact less than half of the 40 can be prime because 2. ...at least 2 of every 1/3/7/9ending quartet are divisible by 3, and for some quartets every member is divisible by 3. Specifically the divisibility pattern for such 1/3/7/9ending quartets is in the form of repeating triplets (where 0 indicates != 0 (mod 3), 1 indicates divisible by 3) ...0101 1010 1111..., thus precisely 4 of every 30 sequence terms starting with the first 30 can possibly be prime. (This is not terribly difficult to prove, but I'll let readers confirm it for themselves, as it's a fun little bit of maths.) 3. The factorizations of the remaining notdivby3 terms appear to be 'random', i.e. modelable by the statistics of randomly chosen odd integers of similar size. 4. Using [13] plus a few more simple observations and some basic number theory we can generate an expected number (or density, if one prefers) of primes for the sequence. However when I do this I get a result which is somewhat at odds with Charles' comment in the notes: Quote:
Here is what I get: The odds of a randomly selected odd integer x being prime is ~2/ln(x) ... summing this for the oddsnotdivisibleby3 for terms 130, 3160 and 6190, &c, each of which intervals contains 4 possiblyprime terms, we get the expected #primes for the first few of said intervals to be: 1 30: 0.22749 31 60: 0.05150 61 90: 0.02700 91120: 0.01809 121150: 0.01242 151180: 0.00939 181210: 0.00755 I did the above by hand (assisted by Pari) ... at this point in order to investigate the convergence (or not) of the estimated #primes I wrote some simple Pari code for playing with this sequence  uncomment the if(isprime...) code snip just below the update of nexpect if you want to check for prime terms, at the cost of drastically increased runtime: Code:
ilog10 = 1/log(10); n = 1; i = 2; nexpect = 0.; while(i < 1000000,\ ndd = ceil(log(i+0.5)*ilog10); /* Need + 0.5 so e.g. ndd(100) comes out = 3 rather than 2 */\ pow10 = 10^ndd;\ n = pow10^2*n + pow10*i + (i+1);\ if(i%1000 == 0,\ print("i = ",i,"; nexpect = ",nexpect);\ );\ if(i%10 != 4, /* Skip terms divisible by 5 */\ if(n%3 != 0, /* Skip terms divisible by 3 */\ nexpect += 2/log(n);\ /* if(isprime(n),\ print("i = ",i+1,": n prime!");\ ); */\ );\ );\ i += 2;\ );\ Here are the results for successive powers of 10 from 10^3 to 10^6  I use logarithmically constant increments here because if the resulting increments in the expectation value decrease from one power of 10 to the next we at least have hope that there may be a limit at infinity: 10^3: 0.4206922620678406265572242819 10^4: 0.4959359595134930290178514034 10^5: 0.5545675055579183966439241436 10^6: 0.6026039035873964125108005995 One might expect the summation to diverge as n > oo based on divergence of the harmonic series  note that even knocking out fixed patterns of terms from the harmonic as we do here using divisibility patterns  does not alter the divergence property. The reason I think the present summation may in fact converge is that due to increasing digit length of the appended numbers, the terms grow faster than log(T_n) ~ n. Thus, rather than the expected #primes being given by a harmonic sum(1/n), which diverges, it is rather given by something which is perhaps like sum(1/n^a) (a.k.a. the pseries or hyperharmonic series) or sum(1/(n * log(n)^a)) (a.k.a. the lnseries, where the summation starts at n = 2 rather than n = 1) with a > 1, both of which converge for all a > 0. Actually, on second thought, considering the logarithmic growth rate of appended numbers, perhaps log(T_n) ~ n log(n) (i.e. the expectation value governed by the lnseries with a = 1, which does in fact diverge, albeit slowly) is the correct asymptotic estimate. However, even if the sum does diverge, it does so sufficiently slowly that the absence of primes in the ranges tested to date should not be surprising. It's certainly an interesting problem, in any event  comments, corrections, further insights appreciated! Last fiddled with by ewmayer on 20151009 at 03:58 Reason: Added n = 10^6 #primes estimate 

20151009, 02:21  #2  
"Forget I exist"
Jul 2009
Dumbassville
20300_{8} Posts 
Quote:
and : Quote:


20151009, 04:00  #3 
∂^{2}ω=0
Sep 2002
República de California
3·13^{2}·23 Posts 
Thanks, sm88  did not see that thread previously, but see nothing useful about #primes estimation in it. The MacPhee comment restates my note about only 4 n's per 30 possibly yielding a prime.

20151009, 12:32  #4 
"Forget I exist"
Jul 2009
Dumbassville
2^{6}×131 Posts 
you can make your Pari code even simpler other than to knock out the remaining multiples of five you can do:
Code:
forstep(x=1,10000,6,...) Code:
parfor(x=1,10000,...;if(x%5!=4,...)) because all the indexes that aren't even and don't divide by 3 are of form 6y+1 and then you can show they are 1,7,13,19,25 mod 30 and 25 mod 30 can always be taken away for the index because it will always end in 5. 
20151017, 21:24  #5  
∂^{2}ω=0
Sep 2002
República de California
3·13^{2}·23 Posts 
The primeodds have just been raised! (Actually the odds remain the same, it is my estimate thereof which has increased.)
Quote:
Of course those are the naive odds  the conditional odds given no primes found below 10^5 are still quite low for the 10^5  10^6 range, around 10%. Last fiddled with by ewmayer on 20151017 at 21:26 

20151017, 21:32  #6  
If I May
"Chris Halsall"
Sep 2002
Barbados
9963_{10} Posts 
Quote:


20151110, 16:33  #7 
Aug 2006
13533_{8} Posts 
My reasoning for the heuristic given above:
A random number would be coprime to 30 for 8/30 of the time. But if you look at the concatenations mod 30 there are only 4 which lead to numbers coprime to 30, so these numbers would be expected to be prime half as often as random numbers. Then all you need is to find the size of the nth term, which is about (n times the average length of a number from 1 to n) decimal digits, or roughly 10^(n*log_{10} n) = exp(n log n), so putting the two together the expected 'chance' that the nth concatenation is prime is 0.5 n log n. Now integrate from n = 1 to x and you get 0.5 log log x. 
Thread Tools  
Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
Status page with expected number of Mersenne primes in each interval?  CRGreathouse  PrimeNet  2  20180110 06:13 
Primes in A048788 OEIS Sequence  carpetpool  Miscellaneous Math  9  20170317 22:57 
The expected number of primes  Batalov  Computer Science & Computational Number Theory  5  20160811 01:17 
I get 13% less primes than I expected:(  mart_r  Math  2  20101029 17:31 
Expected Number of Factors for numbers within a ra  grandpascorpion  Math  2  20071217 13:48 