A small performance tip. a <.@% b will do integer division (<.@%: integer square root), but <. a % b does not.
----- Original Message ----- From: Mike Day <[email protected]> To: [email protected] Cc: Sent: Wednesday, February 25, 2015 11:40 AM Subject: [Jprogramming] Calculation of primepi (similar to _1&p:) for reasonably large scalar arguments Some recent Euler Project problems (to be found, for example, at projecteuler.net/) have required several values of pi(n), a.k.a. primepi(n), the "prime-counting function" for "large" values of n. It is the number of primes less than or equal to n. A convenient J formulation is pi =: 1&p: + _1&p: Unfortunately, J's _1 & p: only works for n up to _1 + 2^31, ie about 2.15e9, while arguments approaching 1e12 are encountered. I solved the problem in question by downloading a public-domain executable, producing a table of required pi(n)'s and reading it into J, but have now produced a reasonable J primepi function to produce the goods. There are several methods for finding primepi; see for example http://mathworld.wolfram.com/PrimeCountingFunction.html I've also coded the Lehmer method in J (NB, there's a mistake in http://mathworld.wolfram.com/LehmersFormula.html , the limits on the second sum should be 1 <= i < j <= a ) . It turns out (in J at least) that my non-Lehmer function outperforms the Lehmer one. I don't claim originality; I adapted the former one from an Euler correspondent's Python code. I've also tried versions in Pari GP and C++, though my skills in both these are best not described as skills. Here's a little table of performance: NB. Times 1000000 1e7 1e8 1e9 1e10 1e11 1e12 1e13 NB. c++~ 1ms 5ms 14ms 150ms 1.1s 9.4s 86s 13m05s NB. J primepi 1ms 6ms 20ms 62ms 3.4s 16.6s 92s 8m15s NB. J Lehmer 4ms 22ms 190ms 2.2s 8m23s ................. NB. Pari GP 100ms 616ms 4.4s 34.6s 4m27s ................. Any thoughts? I'm appending the code for primepi below. Pari GP and c++ source available, as well as my version of Lehmer. Mike Apologies for any line-spacing problems. pi =: (1&p: + _1&p:) @ <. NB. primepi : left argument is optional lower limit for application of NB. primepi rather than pi primepi =: 3 : 0 y primepi~ 2^31 : if. x > n =. y do. pi y return. end. NB. use built-in _1 p: r =. <.%:n if. (n >: *:r) * n < *:r+1 do. V =. <. n%i =. >:i.r V =. }: (,i.@-@{: ) V NB. V+=list(range(... iV =. i.tV =. #v =. V P =. <: V NB. number of primes (+1) NB. S =. <. <:(-:*>:) V NB. sum of primes Vx =. V NB. truncatable arrays Pv =. P NB. ... pi0=. P{~ip =. <:#V NB. plist=. '' NB. can build up partial list of primes here for_p. 2+i.r do. ip =. <: ip pi =. ip{P if. pi > pi0 do. NB. p is prime NB. plist =. plist, p NB. append to partial list tV =. +/v >: <.*:p vx =. i.tV v =. tV {. v vpx=. V I. <.v%p Pv =. tV {. Pv Pv =. pi0 + Pv - vpx{P P =. Pv vx } P NB. smoutput p; pi; pi0;ip; P NB. Sv =. Sv + <. p * si0 - vpx{S NB. similar ops to update sum of primes NB. S =. Sv vx } S NB. si0 is previous ip{S, cf pi0 wrt P end. pi0 =. pi end. V,: P else. 0 end. ) --- This email has been checked for viruses by Avast antivirus software. http://www.avast.com ---------------------------------------------------------------------- For information about J forums see http://www.jsoftware.com/forums.htm ---------------------------------------------------------------------- For information about J forums see http://www.jsoftware.com/forums.htm
