I tried your primepi, and I'm not certain I understand it. It seems to
me that its behavior above the threshold is very different from that
of pi:
pi i. 10
0 0 1 2 2 3 3 4 4 4
pi 2 3 4
1 2 2
primepi 2 3 4
1 2 2
1 primepi 2 3 4
|length error: primepi
| V=.<.n %i=.>:i.r
1 primepi("0) 2 3 4
2 1 0
1 0 0
3 2 1
2 1 0
4 2 1
2 1 0
Was this the code that you meant to post?
Thanks,
--
Raul
On Wed, Feb 25, 2015 at 11:40 AM, Mike Day <[email protected]> wrote:
> 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
