#11475: improve prime_pi (speedup + small fixes)
-----------------------------------------------------------+----------------
Reporter: rohana | Owner:
was
Type: enhancement | Status:
positive_review
Priority: major | Milestone:
sage-5.1
Component: number theory | Resolution:
Keywords: primes, prime counting, prime_pi sd40.5 | Work issues:
Report Upstream: N/A | Reviewers:
Yann Laigle-Chapuy, Leif Leonhardy, William Stein, Karl-Dieter Crisman
Authors: R. Andrew Ohana | Merged in:
Dependencies: | Stopgaps:
-----------------------------------------------------------+----------------
Old description:
> I have rewritten the `prime_pi` method from scratch, it is much cleaner
> and less hacky this time (as I have learned more about coding), however
> it is still based on the same algorithm as before (no LMO yet, although I
> intend to take a stab at it later this year). This was developed in
> parallel with a method to count primes in residue classes (ticket
> forthcoming).
>
> The overhead for small input `>= 2` is rather large in the current
> `prime_pi`, this has been dramatically improved:
>
> New:
> {{{
> sage: timeit('prime_pi(100)')
> 625 loops, best of 3: 475 ns per loop
> sage: timeit('prime_pi(10000)')
> 625 loops, best of 3: 483 ns per loop
> }}}
> Old:
> {{{
> sage: timeit('prime_pi(100)')
> 625 loops, best of 3: 1.71 µs per loop
> sage: timeit('prime_pi(10000)')
> 625 loops, best of 3: 5.56 µs per loop
> }}}
>
> There is also a ~5x speedup (conservatively) for larger input:
>
> New:
> {{{
> sage: timeit('prime_pi(10**8)')
> 625 loops, best of 3: 530 µs per loop
> sage: timeit('prime_pi(10**10)')
> 5 loops, best of 3: 27.1 ms per loop
> sage: time prime_pi(10**12)
> 37607912018
> Time: CPU 1.53 s, Wall: 1.53 s
> }}}
> Old:
> {{{
> sage: timeit('prime_pi(10**8)')
> 125 loops, best of 3: 3.66 ms per loop
> sage: timeit('prime_pi(10**10)')
> 5 loops, best of 3: 161 ms per loop
> sage: time prime_pi(10**12)
> 37607912018
> Time: CPU 8.65 s, Wall: 8.64 s
> }}}
>
> The user can also give a second argument to use additional memory for a
> potential speedup. All primes less than the second argument are used in
> computing pi(x). For example:
>
> {{{
> sage: time prime_pi(10**12, 10**8)
> 37607912018
> Time: CPU 0.92 s, Wall: 0.92 s
> sage: time prime_pi(10**12, 10**8) # pari's prime table is now cached
> 37607912018
> Time: CPU 0.58 s, Wall: 0.58 s
> }}}
>
> An earlier version of this implementation started to fail to give correct
> output on 64-bit systems somewhere between `9.5*10**15` and
> `9.75*10**15`, for issues that I had previously been unable to determine,
> so I had initially limited input to be `< 2**49`, which I am fairly
> confident is safe, as I have done a fair bit of testing in this range
> trying to debug this issue (despite it taking hours per call). I would
> appreciate testing on 32-bit platforms.
>
> A lot of the speed of this algorithm relies on the `__cached_count`
> method, which is essentially a binary search (with a simple adjustment in
> the beginning that takes advantage of the density of the primes). This is
> the fastest type of binary search I know of, but if anyone knows of a way
> to speed it up, it could have a significant effect.
>
> The other big contributor to the speed is the `__smallPi` table, which
> stores the values of pi(x) for `x < 2**16`. This is used in the same
> manner as `__cached_count`, but is of course much faster.
>
> Finally, I have introduced the `legendres_formula` method, which takes
> two arguments, `x` and `a`, and returns the number of positive integers
> not larger than `x` that are not divisible by any -- i.e. co-prime to
> each -- of the first `a` primes. This function is core to all
> combinatorial methods for computing the prime counting function, so I
> figured that I would make this publicly accessible now that I have clean
> code.
>
> ----
>
> Apply [attachment:trac_11475mk2.patch] to the '''Sage library'''.
New description:
I have rewritten the `prime_pi` method from scratch, it is much cleaner
and less hacky this time (as I have learned more about coding), however it
is still based on the same algorithm as before (no LMO yet, although I
intend to take a stab at it later this year). This was developed in
parallel with a method to count primes in residue classes (ticket
forthcoming).
The overhead for small input `>= 2` is rather large in the current
`prime_pi`, this has been dramatically improved:
New:
{{{
sage: timeit('prime_pi(100)')
625 loops, best of 3: 475 ns per loop
sage: timeit('prime_pi(10000)')
625 loops, best of 3: 483 ns per loop
}}}
Old:
{{{
sage: timeit('prime_pi(100)')
625 loops, best of 3: 1.71 µs per loop
sage: timeit('prime_pi(10000)')
625 loops, best of 3: 5.56 µs per loop
}}}
There is also a ~5x speedup (conservatively) for larger input:
New:
{{{
sage: timeit('prime_pi(10**8)')
625 loops, best of 3: 530 µs per loop
sage: timeit('prime_pi(10**10)')
5 loops, best of 3: 27.1 ms per loop
sage: time prime_pi(10**12)
37607912018
Time: CPU 1.53 s, Wall: 1.53 s
}}}
Old:
{{{
sage: timeit('prime_pi(10**8)')
125 loops, best of 3: 3.66 ms per loop
sage: timeit('prime_pi(10**10)')
5 loops, best of 3: 161 ms per loop
sage: time prime_pi(10**12)
37607912018
Time: CPU 8.65 s, Wall: 8.64 s
}}}
The user can also give a second argument to use additional memory for a
potential speedup. All primes less than the second argument are used in
computing pi(x). For example:
{{{
sage: time prime_pi(10**12, 10**8)
37607912018
Time: CPU 0.92 s, Wall: 0.92 s
sage: time prime_pi(10**12, 10**8) # pari's prime table is now cached
37607912018
Time: CPU 0.58 s, Wall: 0.58 s
}}}
An earlier version of this implementation started to fail to give correct
output on 64-bit systems somewhere between `9.5*10**15` and `9.75*10**15`,
for issues that I had previously been unable to determine, so I had
initially limited input to be `< 2**49`, which I am fairly confident is
safe, as I have done a fair bit of testing in this range trying to debug
this issue (despite it taking hours per call). I would appreciate testing
on 32-bit platforms.
A lot of the speed of this algorithm relies on the `__cached_count`
method, which is essentially a binary search (with a simple adjustment in
the beginning that takes advantage of the density of the primes). This is
the fastest type of binary search I know of, but if anyone knows of a way
to speed it up, it could have a significant effect.
The other big contributor to the speed is the `__smallPi` table, which
stores the values of pi(x) for `x < 2**16`. This is used in the same
manner as `__cached_count`, but is of course much faster.
Finally, I have introduced the `legendres_formula` method, which takes two
arguments, `x` and `a`, and returns the number of positive integers not
larger than `x` that are not divisible by any -- i.e. co-prime to each --
of the first `a` primes. This function is core to all combinatorial
methods for computing the prime counting function, so I figured that I
would make this publicly accessible now that I have clean code.
----
Apply [attachment:trac_11475mk2.patch] and
[attachment:trac_11475-reviewer.patch] to the '''Sage library'''.
--
Comment (by kcrisman):
Patchbot: Apply trac_11475mk2.patch and trac_11475-reviewer.patch.
--
Ticket URL: <http://trac.sagemath.org/sage_trac/ticket/11475#comment:62>
Sage <http://www.sagemath.org>
Sage: Creating a Viable Open Source Alternative to Magma, Maple, Mathematica,
and MATLAB
--
You received this message because you are subscribed to the Google Groups
"sage-trac" group.
To post to this group, send email to [email protected].
To unsubscribe from this group, send email to
[email protected].
For more options, visit this group at
http://groups.google.com/group/sage-trac?hl=en.
