Path: newsspool2.news.pas.earthlink.net!stamper.news.pas.earthlink.net!newsread3.news.pas.earthlink.net.POSTED!not-for-mail Message-ID: <3r4c21-42k.ln1@lexi2.athghost7038suus.net> From: The Ghost In The Machine Newsgroups: sci.physics,sci.math,sci.math.num-analysis,alt.writing Subject: Re: My prime research, focus on a feature REVISED Date: Mon, 01 Sep 2003 15:28:12 GMT References: <3c65f87.0308270833.66eedc8e@posting.google.com> <3c65f87.0308271550.325219f5@posting.google.com> <585ab5d8.0308281046.b3ebad5@posting.google.com> <3c65f87.0308300705.403324df@posting.google.com> <3c65f87.0308301524.4b58602@posting.google.com> <3F513626.771F2246@ix.netcom.com> <3c65f87.0308310458.3bb6b21d@posting.google.com> <2ss921-alh.ln1@lexi2.athghost7038suus.net> <3c65f87.0308311520.1ecf4660@posting.google.com> Lines: 171 Followup-To: sci.math X-face: "i;@/WO(?;[KC9sW;wG/4@H[_VFFH4?QHJ#O(?m}7fQMrJ,]0THA'\|e-EPG_>56Mi}_RRhBS'a2}u_7jm)0_+'=$V#E2r4#IQE/d)yMv3_4@hl<)mA&*tDN/ User-Agent: slrn/0.9.7.4 (Linux) NNTP-Posting-Host: 165.247.205.101 X-Complaints-To: abuse@earthlink.net X-Trace: newsread3.news.pas.earthlink.net 1062430092 165.247.205.101 (Mon, 01 Sep 2003 08:28:12 PDT) NNTP-Posting-Date: Mon, 01 Sep 2003 08:28:12 PDT Organization: EarthLink Inc. -- http://www.EarthLink.net X-Received-Date: Mon, 01 Sep 2003 08:28:12 PDT (newsspool2.news.pas.earthlink.net) Xref: lexi2.athghost7038suus.net sci.physics:76354 sci.math:35013 Followups reset. In sci.physics, Christian Bau wrote on Mon, 01 Sep 2003 09:39:11 +0100 : > In article , > The Ghost In The Machine wrote: > >> It probably is, since Legendre's Formula (according to >> >> http://mathworld.wolfram.com/LegendresFormula.html >> >> anyway) requires p_i to drive it (p_1 = 2, p_2 = >> 3, etc.). The recurrence relation mentioned: >> >> phi(x,n) = phi(x, n-1) - phi(x / p_n, n - 1) >> >> may be of some use to somebody, with the initial conditions >> phi(x, 0) = floor(x) for any x >= 0 and phi(x, n) = >> 0 for any x < 2 and n >= 0. I'm estimating at least >> x^(1/2) space complexity, mostly because of the prime >> table; however, initial runs suggest far more than that. >> Meissel's or Lehmer's formulas might yield better results >> but would still require the prime table. The author >> asserts phi(x,x) = pi(x); I don't know if phi(x,sqrt(x)) >> would work better or not. > > The recurrence relation is _exactly_ what you should use in an > implementation, it very nicely adds up all the terms in Legendre's > formula and conveniently leaves those out that evaluate to 0. > > Note that the definition of Legendre's formula on the webpage you quote > seems completely messed up: First, phi (x, a) is correctly defined as > the number of integers <= x that are not divisible by the first primes; > the summation formula then doesn't write how far the sums are supposed > to go. In the second formula, you don't choose a = x, that is nonsense; > you choose a = pi (floor (sqrt (x))). > > Once you fix the obvious bugs on the webpage, you get Legendres formula >= Harris' formula. > > As to space requirements: Who says that you should use ONE sieve to find > all the primes <= sqrt (x) ? Use a sieve to find all primes <= x ^ > (1/3), then observe that all primes > x^(1/3) are used only once if you > use the recursion formula, so there is no need to store them. Use > several sieves of size x^(1/3) to find those primes, and your space > requirement is O (N^(1/3)). Well, something's still messed up. I feel like a schoolkid again... :-) Here's *my* program, which is fairly fast but gives out the wrong answer. (The source should be fairly obvious even though I'm using STL.) ----------8< >8---------- #include #include #include std::vector primeArray; // Routine to compute the ix'th prime. static unsigned prime(unsigned ix) { while(ix >= primeArray.size()) { // glory be, we've run out of primes! unsigned last = primeArray[primeArray.size()-1]; bool foundPrime = false; while(!foundPrime) { if(last % 6 == 1) last += 4; else last += 2; // don't bother testing 2 and 3 for(unsigned j = 3; j < primeArray.size(); j++) { if(primeArray[j] * primeArray[j] > last) { // we've run the gauntlet and know 'last' is prime primeArray.push_back(last); foundPrime = true; break; } if(last % primeArray[j] == 0) { break; // not prime, sorry } } } } return primeArray[ix]; } // Routine to compute the number of primes less than or equal to lim // This routine has a number of latent bugs which are not elicited // by the rest of the implementation; a more general routine would // use a binary search if the last prime is greater than the limit. static unsigned stupidpi(unsigned lim) { while(primeArray[primeArray.size()-1] < lim) prime(primeArray.size()); return primeArray.size() - 2; } // Legendre's function. static double phi(double x, unsigned n) { if(x < 2) return 0; if(n == 0) return floor(x); // Since the primes are strictly monotonic, if we know // the next iteration is going to be phi(y,something) where // y < 2, don't bother explicitly computing prime(n). if(n > primeArray.size() && primeArray[primeArray.size() - 1] > x/2) return phi(x, n-1); return phi(x, n-1) - phi(x / prime(n), n-1); } // Main program. int main(int argc, char **argv) { primeArray.push_back(1); // p_0 is not defined so it doesn't matter anyway primeArray.push_back(2); primeArray.push_back(3); primeArray.push_back(5); int a = stupidpi(1000); std::cout << "Precalculation done; a = " << a << std::endl; std::cout << phi(1000000,a) << std::endl; std::cout << primeArray.size() << " " << primeArray[primeArray.size()-1] << std::endl; int b = stupidpi(1000000); std::cout << b << std::endl; return 0; } ----------8< >8---------- It gives me the wrong answers, so I've obviously missed something. Sigh. $ legendrephi Precalculation done; a = 168 81142 170 1009 78498 $ sieve 1000000 1000000 78498 $ -- #191, ewill3@earthlink.net -- but at least it's fast :-) It's still legal to go .sigless.