Opened 2 years ago

Last modified 22 months ago

#24960 new defect

prime_pi errors

Reported by: gianluca.amato.74 Owned by:
Priority: minor Milestone: sage-8.2
Component: number theory Keywords: primes prime_pi
Cc: Merged in:
Authors: Reviewers:
Report Upstream: N/A Work issues:
Branch: Commit:
Dependencies: #24966, #25009 Stopgaps:

Description (last modified by gianluca.amato.74)

The prime_pi function sometimes computes a wrong result. In particular, I found the following errors:

  • prime_pi(14783545363230719) returns 408351706116074, it should be 408351706116078
  • prime_pi(5685979153167559) returns 161319877461134, it should be 161319877461136
  • prime_pi(642763101936913) returns 19439675999018, it should be 19439675999019

The correct countings are computed using Kim Walisch's primecount library, with different algorithms: Deleglise-Rivat, Lagarias-Miller-Odlyzko and Lehmer's formula. All of them return the same result.

Note that while the arguments to prime_pi in the first two examples exceed the value 2^50, causing the appearance of a warning ("computation of prime_pi for x >= 2^50 has not been as thoroughly tested as for smaller values"), the argument of the last example is within the range of inputs considered to be safe.

Actually, instead of fixing the current code, it could be easier to rewrite the prime_pi function and let it use the primecount library. This would have the additional advantage of being much faster than the current code. The table below compares the execution times (in seconds) of primecount with Deleglise-Rivat, primecount with Lagarias-Miller-Odlyzko and SageMath prime_pi function.

All results have been obtained with an Intel i5-2500K. For primecount, a single thread has been used.

algorithm 642763101936913 5685979153167559 14783545363230719
primecount 2.99 11.168 20.052
primecount --lmo 5.287 21.883 41.286
prime_pi 331 2489 6083

Change History (8)

comment:1 Changed 2 years ago by vdelecroix

Could you add timings in the ticket description? Did you check with PARI/GP?

First step would be to make primecount a proper optional package and interface it. In a second time it can be discussed whether this needs to be a standard package (ie coming with default Sage installation). As a parallel task, the fact that prime_pi is wrong is really bad and could be temporarilly fixed with a [stopgap].

comment:2 Changed 2 years ago by vdelecroix

  • Dependencies set to #24966

See #24966 for actually packaging primecount.

comment:3 Changed 2 years ago by gianluca.amato.74

  • Description modified (diff)

I changed the description of the ticket adding execution times. I did not try using Pari/GP because Pari primepi function is very very slow (at least until the latest release, there is some minor improvement in the development code). It works by enumerating one by one all prime numbers, together with a table of pre-computed values. Just to give an example, although primepi(200000000507)=8007105083 is pre-computed. computing primepi(251000000000) requires 372 seconds, while SageMath and primecount only take less than a second.

comment:4 Changed 2 years ago by gianluca.amato.74

  • Description modified (diff)

Added some context information to the ticket description.

comment:5 Changed 2 years ago by ohanar

When I wrote the current prime_pi code, there was no open source implementations of fast prime counting. I agree that the best long term solution is to use the primecount library.

comment:6 Changed 2 years ago by gianluca.amato.74

  • Description modified (diff)

Fixed typos.

comment:7 Changed 2 years ago by vdelecroix

  • Dependencies changed from #24966 to #24966, #25009

comment:8 Changed 22 months ago by vdelecroix

I can confirm the behavior of the last example

sage: %time prime_pi(642763101936913)
CPU times: user 4min 29s, sys: 16.1 ms, total: 4min 29s
Wall time: 4min 29s

and thanks to #24966 this can be checked with

sage: from sage.interfaces import primecount
sage: %time primecount.prime_pi(642763101936913)
CPU times: user 4.88 s, sys: 61.8 ms, total: 4.94 s
Wall time: 1.24 s
Note: See TracTickets for help on using tickets.