User group for PFGW & PrimeForm programs Yahoo Group number of prime-index-primes < n =============================================== cino hilliard Message 1 of 15 Apr 2, 2006 ----------------------------------------------- Hi, I have found an interesting function which I call pipips(n) . It is the count of prime-index-primes (pips) less than n. We define pip(n) = prime(prime(n)) where prime(n) is the n-th prime. The interesting thing about this is that pipip(n) ~ R(R(n))= pipipR(n) where R(n) is the Riemann approx to Pi(n). Below is the know values and error terms of pipip(n) and pipipR(n) for n in powers of 10. n pipip(n) pipipR(n) error relative error 10^1 2 3 -1 - 0.5 10^2 9 9 0 0.0 10^3 39 39 0 0.0 10^4 201 200 1 0.004975124 10^5 1184 1182 2 0.001689189 10^6 7702 7707 -5 - 0.000649182 10^7 53911 53929 -18 - 0.000333884 10^8 397557 397578 -21 - 0.000052823 10^9 3048956 3048846 +110 0.000036078 10^10 24106416 24106273 +142 0.000005932 10^11 195296943 195297175 -232 - 0.000001188 10^12 1613846646 1613844824 +1822 0.000001129 10^13 13556756261 13556753359 +2902 0.000000214 I used Pari and a 105 gig file of primes < 300 billion to compute through 10^11 I then used the nth prime page to do 10^12 and 10^13. If we can compute prime(n) we can easily compute pipip(n) exactly by computing R_0 = pipipR(n) = R(R(n)) and comparing prime(prime(R_0)) with n and do binary search adjusting until it crosses n. This relationship continues with higher orders such as pipips3(n) ~ pipipsR3(n) ~ R(R(R(n))) With lesser and lesser error. It is incredible to see such regularity in primes while knowing that there exist arbitrarily large gaps of composite numbers between two primes. I am looking for a fast program to compute prime(n) so we can exted this table to see (heuristically) if the diminishing error continues or surprizes us. Also comments and elaborations on this are welcome. Have fun, Cino *********************Code ******************************* Pari definitions \\Brute force pipips(n) = \\The number of prime-index-primes (pips) less than or equal to n { local(x,y); x=0; y=0; while(y<=n, x++; y=prime(prime(x)); ); return(x-1); } pipipsR(n) = \\The number of prime-index-primes (pips) less than or equal to n \\Uses the Riemann Integral approximation to pi(n) { local(x,y,yR,ppyR); x=0; y=0; yR = R(R(n)); ppyR =prime2(prime2(yR)); if(ppyR > n, while(prime2(prime2(yR)) > n, yR -=1); ); if(ppyR < n, while(prime2(prime2(yR)) < n, yR +=1); if(n <= 10^5,yR--); ); return(yR); } R(x) = \\ Riemann's approx of Pi(x) { local(j); round(sum(j=1,200,moebius(j)*Li(x^(1/j))/j)) } =============================================== David Broadhurst Message 2 of 15 Apr 2, 2006 ----------------------------------------------- --- In primeform@yahoogroups.com, "cino hilliard" wrote: > 10^13 > 13556756261 > 13556753359 to which one may quickly add pi(pi(10^14)=115465507935 [R(R(10^14)]=115465491695 pi(pi(10^15)=995112599484 [R(R(10^15)]=995112567449 But why confine the comparison to the integer part of R(R(n))? Finer detail is revealed by comparing pi(n) directly with [R(n)], as in Ribenboim Table 27, with 28 entries from n=10^15 to n=10^18 David =============================================== David Broadhurst Message 3 of 15 Apr 2, 2006 ----------------------------------------------- --- In primeform@yahoogroups.com, "cino hilliard" wrote: > R(x) = \\ Riemann's approx of Pi(x) > { > local(j); > round(sum(j=1,200,moebius(j)*Li(x^(1/j))/j)) > } I believe that it is more efficient to use Gram's formula, which gives R(R(10^15)) to 28-digit accuracy in less than 0.01 seconds. \p28 {R(x)=local(n=1,L,s=1,r);L=log(x);r=L; while(s<10^30*r,s=s+r/zeta(n+1)/n;n=n+1;r=r*L/n);s} gettime;print(R(R(10^15)));print(gettime" ms"); 995112567449.7996774369207682 9 ms Note that there is no "round" or "floor" in R(x); it is an analytic function. David =============================================== Andrey Kulsha Message 4 of 15 Apr 2, 2006 ----------------------------------------------- > I believe that it is more efficient to use Gram's formula, > which gives R(R(10^15)) to 28-digit accuracy in > less than 0.01 seconds. > > \p28 > > {R(x)=local(n=1,L,s=1,r);L=log(x);r=L; > while(s<10^30*r,s=s+r/zeta(n+1)/n;n=n+1;r=r*L/n);s} > > gettime;print(R(R(10^15)));print(gettime" ms"); > > 995112567449.7996774369207682 > 9 ms > > Note that there is no "round" or "floor" in R(x); > it is an analytic function. But there's a bit better estimator than R(x): R(x) + ArcTan[Pi/L]/Pi - 1/L, where L = log[x], Pi = 3.14159265358979... Andrey :) [Non-text portions of this message have been removed] =============================================== David Broadhurst Message 5 of 15 Apr 3, 2006 ----------------------------------------------- --- In primeform@yahoogroups.com, "Andrey Kulsha" wrote: > But there's a bit better estimator than R(x): > > R(x) + ArcTan[Pi/L]/Pi - 1/L, > > where L = log[x], Pi = 3.14159265358979... That's a tiny O(1/log(x)^3) = o(1) addendum. The actual discrepancies between R(x) and pi(x) with x in [10^15,10^8] are far larger than your itsy-bitsy footnote :-) David =============================================== Andrey Kulsha Message 6 of 15 Apr 3, 2006 ----------------------------------------------- > That's a tiny O(1/log(x)^3) = o(1) addendum. > > The actual discrepancies between R(x) and pi(x) > with x in [10^15,10^8] are far larger > than your itsy-bitsy footnote :-) But for small x's it is quite significant. :) http://www.primefan.ru/stuff/primes/plot_r.gif There are normalized (i.e. multiplied by sqrt(x)/log(x)) differencies between pi(x) and our two estimators on log(x) scale. The lower one corresponds to R(x) ;) Andrey [Non-text portions of this message have been removed] =============================================== cino hilliard Message 7 of 15 Apr 3, 2006 ----------------------------------------------- >From: "David Broadhurst" >Reply-To: primeform@yahoogroups.com >To: primeform@yahoogroups.com >Subject: [primeform] Re: number of prime-index-primes < n >Date: Sun, 02 Apr 2006 13:59:05 -0000 > >--- In primeform@yahoogroups.com, "cino hilliard" >wrote: > > > R(x) = \\ Riemann's approx of Pi(x) > > { > > local(j); > > round(sum(j=1,200,moebius(j)*Li(x^(1/j))/j)) > > } > >I believe that it is more efficient to use Gram's formula, >which gives R(R(10^15)) to 28-digit accuracy in >less than 0.01 seconds. I now have the following functions in my util.gp which is called in the .gprc file. Rg(x) = \\ Gram's Riemann's Approx of Pi(x) { local(n=1,L,s=1,r); L=r=log(x); while(s<10^30*r, s=s+r/zeta(n+1)/n; n=n+1; r=r*L/n); round(s) } Rgb(x) = \\ Gram's Riemann's Approx of Pi(x) Broadhurst style { local(n=1,L,s=1,r); L=r=log(x); while(s<10^30*r, s=s+r/zeta(n+1)/n; n=n+1; r=r*L/n); s } >gettime;print(R(R(10^15)));print(gettime" ms"); > >995112567449.7996774369207682 >9 ms Indeed, P4 2.53 ghz 2 gig ram windows xp pro GP/PARI CALCULATOR Version 2.2.12 (beta) i686 running cygwin (ix86 kernel) 32-bit version compiled: Jan 3 2006, gcc-3.4.1 (cygming special) (readline v5.0 enabled, extended help available) (13:13) gp > Rg(Rg(10^30)) %47 = 230140690020794659839995596 (13:13) gp > ## *** last result computed in 16 ms. (13:13) gp > R(R(10^30)) %48 = 230140690020794659839995596 (13:14) gp > ## *** last result computed in 31 ms. > >Note that there is no "round" or "floor" in R(x); >it is an analytic function. True. Also true 3.13159... is Real and so is 3 and 3.000000... etc. I have defined this approx as such. > > R(x) = \\ Riemann's approx of Pi(x) When we say approximatioin, in my view, it does not matter that R(x) is analytic. Origionally I let the thing rip and ran into problems with real numbers in some Pari functions. (13:26) gp > B=Rgb(Rgb(10^15)) %53 = 995112567449.7996774369207682 (13:27) gp > isprime(B) *** isprime: not an integer argument in an arithmetic function (13:28) gp > ispseudoprime(B) *** ispseudoprime: not an integer argument in an arithmetic function (13:28) gp > B! *** this should be an integer: B! ^-- (13:28) gp > ispower(B) *** ispower: missing exponent. (13:29) gp > Etc My choice to round or truncate is arbitrary. I do not believe it plays a part in the grand scheme of things. I believe Real Rgb(x) does and is beyond my scope. Can there be x's such that round(Rgb(Rgb(x))) != Rg(Rg(x)) ? Indeed, there are many. The frequency of approximations to the exact sum of pips < x is about the same for Rg(Rg(x)) and round(Rgb(Rgb(x)). Cheers and Roebuck, Cino =============================================== David Broadhurst Message 8 of 15 Apr 3, 2006 ----------------------------------------------- Cino: > When we say approximation, in my view, > it does not matter that R(x) is analytic. Please do not forget that there is an /exact/ formula for pi(x) in terms of Riemann's analytic function. When you mutilated it, by rounding, I felt, perhaps irrationally, that you were despoiling one of the most impressive monuments of mathematics :-) David per proxy SPMM Society for the Preservation of Mathematical Monuments =============================================== cino hilliard Message 9 of 15 Apr 4, 2006 ----------------------------------------------- Hi David, Yo have got my cognitive juces are flowing. >From: "David Broadhurst" >Date: Tue, 04 Apr 2006 01:09:43 -0000 > > When we say approximation, in my view, > > it does not matter that R(x) is analytic. > >Please do not forget that there is an /exact/ formula >for pi(x) in terms of Riemann's analytic function. Ok, use it to compute pi(10^22) pi(10^23) ... ... pi(10^50) Show me the real numbers or provide the pari script to do it. I am from Missouri at this point. >one of the most impressive monuments of mathematics :-) This will be the case when the Riemann hypothesis is proved. I will have to say though, that guy of 40 years had one heck of a mind. He was almost as brilliant as your buddy, Newton.:-) I did a lot of estimating when I worked in Retai.l We had a complicated exponential smoothing forecast formula which I simplified to R = [(A/B - C)*alpha + C]*D + E Where R = Seasonal Revised Estimate through current period. A = Current period sales (week, 2-weeks, days, hours etc) B = Expected current period % contribution to total sales (from history or theoretical trends) C = Previous period Seasonal Revised estimate. Orgional est in the first period. D = Cumulative balance to go percent from current period to end (1 - percent done). E = Cumulative sales through the current period. alpha = exponential smoothing constant (usually 0.1 to start). For example lets say we have a sales event of 10 weeks and are 1 week into the event where A = 103 B = 0.05 C = 3000 D = 0.95 E = 103 alpha = 0.1 Then R = [(103/.05 -3000)*0.1 + 3000]*0.95 +103 = 2863.7. The system retained the decimals to 4 places. But we humans referred to the extimate as 2864 to upper management. It is folly to express estimate to a precision greater than the margin of errors in the data. In the case of estimates, the compound relative error in B,D and C could be > 30%. So even the 2864 is pure folly and we often rounded into hundreds or thousands. In words we project the current period sales to get a new seasonal projection. We compare this projection with the last estimate by taking the difference to see how far off we were from our last estimate. We are not sure how much we want to change the last estimate so we only change the last estimate by a percentage of this difference or alpha. Since we want to rely more heavily on actual sales through time, we project the seasons estimate on a balance to go basis by multiplying our projection by the balance to go %. We then add the actual sales to get the revised estimate. In this way we more and more smooth our estimate with the nitty grity - the actual sales. So as we get to the end of the event periodic sales cannot create a revised estimate greater than what we sold. It seems to me we should be able to treat P(x) as cumulative sales through periods x and use R(x) to build a trend for comparative projection. I have tried to apply this algorithm to estimate out pi(x) to little avail. Maybe someone can set it up better. If you recall, the origional table I presented I had for n = 10^13 pipips = 3555756261 actual sum of prime-index-primes and R(R(n)) = 3556753359 for a difference of +2902. This error completely swamps any mistake of +/- 1 wel will make by rounding. So it is folly to estimate the count of pips < 10^n for n > 6. This reminds me of one of my Re-Buyers in the 1970's who estimated his sales quantity by location using a slide rule and then went down the hall to use, at that time, one of the rare calculators to compute the percent contribution of each location to 4 decimal places. I tried to explain to him the folly but he did not understand. Hey, the calculator gave 4 places why not use it? >When you mutilated it, by rounding, I felt, >perhaps irrationally, that you were despoiling I did not intend to mutilate anything. I was estimating integer results with integers. That is consistent in any field of inquiry. However, it is ALSO correct to say there are approximately 201467286689188773625.1590116 primes < 10^22. But it is folly. We know the absolute error from the actual 201467286689315906290 is 127132664.84098835289. Rounding R(10^22) at this point is just as good as the real calculation. In fact, 201467286689000000000 is as good an estimate as any not knowing the exact number. " A chain is no stronger than its weakest link." I have uploaded a table of 30,000 p(x),Li(x),R(x) up to 30 trillion in increments of 1 billion. Notice I have used decimal values. http://groups.yahoo.com/group/primeform/files/Pi%28x%29%20data/ The table was built in Excel using Macro Pro, the Nth Prime Page for Pi(x) and Pari to build Li(x) and R(x). http://primes.utm.edu/nthprime/index.php Cheers and Roebuck, Cino =============================================== David Broadhurst Message 10 of 15 Apr 4, 2006 ----------------------------------------------- Cino: > > one of the most impressive monuments of mathematics :-) > This will be the case when the Riemann hypothesis is proved. The formula is true even if RH is not. Oh, a little joke, regarding your favourite state: http://www.sos.mo.gov/archives/history/slogan.asp I prefer the example of the more monumental Vandiver: MR0342346 (49 #7092) Lehmer, D. H. Harry Schultz Vandiver, 1822--1973. Bull. Amer. Math. Soc. 80 (1974), 817--818. 01A70 (10-03) A brief scientific biography of Vandiver who, though best known for his work on Fermat's last theorem, was a self-taught mathematician who enjoyed a long academic career (ill health forced his retirement from the University of Texas at the age of 84) during which he produced 174 publications covering a wide range of subjects in number theory and algebra. Hmm, forced to retire at 84. Poor fellow. David =============================================== cino hilliard Message 11 of 15 Apr 4, 2006 ----------------------------------------------- >From: "David Broadhurst" >Reply-To: primeform@yahoogroups.com >To: primeform@yahoogroups.com >Subject: [primeform] Re: number of prime-index-primes < n >Date: Tue, 04 Apr 2006 15:26:15 -0000 > >Cino: > > > > one of the most impressive monuments of mathematics :-) > > This will be the case when the Riemann hypothesis is proved. > >The formula is true even if RH is not. > >Oh, a little joke, regarding your favourite state: > >http://www.sos.mo.gov/archives/history/slogan.asp Interesting. I am glad "you showed me." Still waiting for the "exact" expansion for pi(10^50) > >I prefer the example of the more monumental Vandiver: > >MR0342346 (49 #7092) >Lehmer, D. H. >Harry Schultz Vandiver, 1822--1973. >Bull. Amer. Math. Soc. 80 (1974), 817--818. >01A70 (10-03) > >A brief scientific biography of Vandiver who, though best known for >his work on Fermat's last theorem, was a self-taught mathematician who >enjoyed a long academic career (ill health forced his retirement from >the University of Texas at the age of 84) during which he produced 174 >publications covering a wide range of subjects in number theory and >algebra. > >Hmm, forced to retire at 84. Poor fellow. David, you are a cornucopia of information. Most people will go through their entire life without knowing this. Cheers and Roebuck, Cino =============================================== David Broadhurst Message 12 of 15 Apr 4, 2006 ----------------------------------------------- --- In primeform@yahoogroups.com, "cino hilliard" wrote: > Still waiting for the "exact" expansion for pi(10^50) Ah, finally you asked precisely the right question :-) The exact expansion is pi(10^50) = R(10^50) - sum_i R(10^(50*z_i)) where the sum extends over all complex zeros, z_i, of the Riemann zeta function, each counted with its multiplicity. This /expansion/ is indeed /exact/ even if RH is false. RH concerns the issue of where those complex zeros actually lie :-) David =============================================== David Broadhurst Message 13 of 15 Apr 4, 2006 ----------------------------------------------- --- In primeform@yahoogroups.com, "cino hilliard" wrote: > > Hmm, forced to retire at 84. Poor fellow. > David, you are a cornucopia of information. Most people will > go through their entire life without knowing this. Cino: Here's another monumental curio. In 1896, proofs of the prime Number Theorem were given by Jacques Salomon Hadamard (1865-1963) Charles Jean Gustave Nicolas Baron de la Vallée Poussin (1866-1962) Observation [Ribenboim, p.225] If we take the average of their ages at death we obtain the largest prime with precisely 2 decimal digits. David =============================================== cino hilliard Message 14 of 15 Apr 5, 2006 ----------------------------------------------- >Cino: Here's another monumental curio. > >In 1896, proofs of the prime Number Theorem were given by > >Jacques Salomon Hadamard >(1865-1963) > >Charles Jean Gustave Nicolas Baron de la Vall�e Poussin >(1866-1962) > >Observation [Ribenboim, p.225] >If we take the average of their ages at death we obtain >the largest prime with precisely 2 decimal digits. Hi David, I like Fermat's birth and death date curios - born prime and died decomposable. Sorry about that.:-( Observation [El Cino, mathforum] Birth = 1601 = 40^2 +1^2 a prime = 801^2 - 800^2 Death = 1665 = 12^2 + 39^2 = 24^2 + 33^2 = 41^2 - 4^2 = 97^2 - 88^2 = 63^2 - 48^2 = 169^2 - 164^2 = 279^2 - 276^2 = 833^2 - 832^2 = 40^2 + 8^2 + 1^2 = 40^2 + 7^2 + 4^2 etc. Birth - Death = 8^2 = 4^3 = 2^6. And... 2^(2^6)+1 = 274177*67280421310721 contrary to Fermat's thinking. Whew! I will stop. Fermat was indeed blessed with numbers! Unfortunately, the poor chap was not able to deduce the last few cases.:-) El Cino =============================================== David Broadhurst Message 15 of 15 Apr 5, 2006 ----------------------------------------------- I wrote > Jacques Salomon Hadamard > (1865-1963) Note however that Paulo had made a mistake: Hadamard died at the prime age of 97, not 98: > Born: 8 Dec 1865 in Versailles, France > Died: 17 Oct 1963 in Paris, France print(nextprime(98)) 101 Cartan, fils, is I believe still with us, at the prime age of 101: http://www-gap.dcs.st-and.ac.uk/~history/Biographies/Cartan_Henri.html print(nextprime(102)) 103 An historian of mathematics, Dirk Jan Struik, reached 106: http://www-gap.dcs.st-and.ac.uk/~history/Biographies/Struik.html print(nextprime(106)) 107 Any candidate? David =============================================== Cached by Georg Fischer at Nov 14 2019 12:47 with clean_yahoo.pl V1.4