Fix prime_pi errors and make primecount standard

[sagetrac-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.

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

---

Here we take the route of replacing Cython code in sage.functions.prime_pi
with calls to functions in sage.interface.primecount

Will be fixed by #32894

Depends on #25009

CC: @slel

Component: number theory

Keywords: primes prime_pi primecount

Branch/Commit: u/dimpase/functions/prime_pi-fix @ 4f77073

Issue created by migration from https://trac.sagemath.org/ticket/24960

[videlec]

comment:1

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 http://doc.sagemath.org/html/en/reference/misc/sage/misc/stopgap.html].

[videlec]

Dependencies: #24966

[videlec]

comment:2

See #24966 for actually packaging primecount.

[sagetrac-gianluca-amato-74]

comment:3

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.

[sagetrac-gianluca-amato-74]

comment:4

Added some context information to the ticket description.

[ohanar]

comment:5

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.

[sagetrac-gianluca-amato-74]

comment:6

Fixed typos.

[videlec]

Changed dependencies from #24966 to #24966, #25009

[videlec]

comment:8

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
19439675999018

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
19439675999019

[miguelmarco]

comment:9

Since primecount has been an optional package for several years, maybe it would be a good idea to upgrade it to standard and just call it?

[williamstein]

comment:10

Or if nobody has time to work on this, maybe somebody could replace the message ""computation of prime_pi for x >= 2^50 has not been as thoroughly tested as for smaller values", which is only triggered for large input, by a similar message triggered by any x >= 10^6 say? The message could be "computation of prime_pi has known bugs (see #24960)".

[miguelmarco]

comment:11

Replying to @miguelmarco:

Since primecount has been an optional package for several years, maybe it would be a good idea to upgrade it to standard and just call it?

Of course, in that case we should reconsider #32412

[dimpase]

comment:12

I think we should upgrade primecount spkg and use it in prime_pi. And revert #32412

[mkoeppe]

comment:13

See #32412 comment:7 for how to do the primecount Cython interfaces properly

[miguelmarco]

comment:14

Replying to @dimpase:

I think we should upgrade primecount spkg and use it in prime_pi. And revert #32412

+1

The difference in speed in huge. Not to mention the correctness.

Is the cmake dependency really a no-go for making it standard?

[mkoeppe]

comment:15

Replying to @miguelmarco:

Is the cmake dependency really a no-go for making it standard?

No, that's outdated. cmake is a standard package now.

[miguelmarco]

comment:16

Replying to @mkoeppe:

See #32412 comment:7 for how to do the primecount Cython interfaces properly

So we should move the current cython code to a pip package, and then instruct sage to pip-install it instead of using a sage SPKG?

[mkoeppe]

comment:17

Yes, either the current cython code should become a pip-installable package; or switch to using https://pypi.org/project/primecount/ (but see hearot/primecount-python#3).

The pip-installable package can then become an SPKG.

[dimpase]

comment:18

I think it's more urgent to provide a bug-free prime_pi than to package the implementation nicely.

[mkoeppe]

comment:19

It's 10 minutes of work to package it.

[dimpase]

comment:20

if you know what to do, yes.
if not, who knows.

[mkoeppe]

comment:21

In this case, one copies and adjusts something that works, such as https://github.com/sagemath/memory_allocator

[dimpase]

comment:22

I've posted some not yet complete thing on #25009

[dimpase]

comment:26

one has to decide what to do with the optional (usless) 2nd parameter to prime_pi()

---

New commits:

1e24a4a	bump primecount to ver. 7.1, and standard
a660ee4	removed deprecations
398748d	prime_pi and legendre_phi calling primecount

[slel]

Changed keywords from primes prime_pi to primes prime_pi primecount

[dimpase]

comment:28

should primecount to be added to dependencies of, hmm, sagelib ?
Because on a slow 32-bit machine I see how building docs is failing (if I start make ptest right after ./configure), as it needs prime_pi, which in turn depends on primecount, or, more precisely, on sage.interfaces.primecount.

[mkoeppe]

comment:30

comment:21.

[dimpase]

comment:31

I finally built and tested this on a 32-bit host, the only one test error there, it is due to a different error message, int is too big vs Python int is too big.

Needs some dots to be added.

[dimpase]

comment:32

Replying to @mkoeppe:

comment:21.

I don't understand, do you mean that primecount is the only package that needs this treatment?

[sagetrac-git]

Changed commit from 398748d to 0c3be43

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. New commits:

0c3be43

allow for slightly different Overflow messages in doctest

[dimpase]

comment:34

Do we also need

--- a/build/pkgs/sagelib/dependencies
+++ b/build/pkgs/sagelib/dependencies
@@ -1,4 +1,4 @@
-FORCE $(SCRIPTS) arb boost_cropped $(BLAS) brial cliquer cypari cysignals cython ecl eclib ecm flint libgd gap giac givaro glpk gmpy2 gsl iml jinja2 jupyter_core lcalc lrcalc libbraiding libhomfly libpng linbox m4ri m4rie memory_allocator mpc mpfi mpfr $(MP_LIBRARY) ntl numpy pari pip pkgconfig planarity ppl pplpy pycygwin $(PYTHON) ratpoints rw sage_conf singular symmetrica zn_poly $(PCFILES) | $(PYTHON_TOOLCHAIN) sage_setup
+FORCE $(SCRIPTS) arb boost_cropped $(BLAS) brial cliquer cypari cysignals cython ecl eclib ecm flint libgd gap giac givaro glpk gmpy2 gsl iml jinja2 jupyter_core lcalc lrcalc libbraiding libhomfly libpng linbox m4ri m4rie memory_allocator mpc mpfi mpfr $(MP_LIBRARY) ntl numpy pari pip pkgconfig planarity ppl pplpy primecount pycygwin $(PYTHON) ratpoints rw sage_conf singular symmetrica zn_poly $(PCFILES) | $(PYTHON_TOOLCHAIN) sage_setup
 
 ----------
 All lines of this file are ignored except the first.
diff --git a/src/pyproject.toml.m4 b/src/pyproject.toml.m4
index 125da7ac1c..6b6402884a 100644
--- a/src/pyproject.toml.m4
+++ b/src/pyproject.toml.m4
@@ -18,6 +18,7 @@ requires = [
         numpy          \
         pkgconfig      \
         pplpy          \
+        primecount     \
         memory_allocator \
                     ')]
 build-backend = "setuptools.build_meta"

to avoid possible races? Still, this does not seem to be relevant for the issue in comment:28.

[mkoeppe]

comment:35

Replying to @dimpase:

Replying to @mkoeppe:

comment:21.

I don't understand, do you mean that primecount is the only package that needs this treatment?

You were asking how to do dependencies right. The answer is: The bindings need to go into a separate Python package -- that will look just like memory_allocator. In this way, sagelib will not acquire another compile-time dependency. Modularization = do not pile on additional compile-time dependencies.

[dimpase]

comment:36

I think we're talking about different issues. I saw an error during docbuild stage which came from
src/sage/interfaces/primecount.pyx not yet built at that moment. The external library spkg primecount was already built at that moment.

[mkoeppe]

comment:37

As I have explained, src/sage/interfaces/primecount.pyx should not be part of sagelib. This is why it was deprecated!

[dimpase]

comment:38

I asked you whether it's the only Cython file out there left that should not be a part of sagelib.
You didn't reply - I assume this means it is not.

I don't have mental capacity now to learn the ideal modern design, I just want this ticket done in the most direct way.
Modularisation can be done on another ticket.

[dimpase]

comment:39

I was relieved today from teaching computer graphics next term, so I have some spare mental capacity now :-)

Note that src/sage/functions/prime_pi.pyx relies heavily on sagelib (e.g. on SR).

What does not rely on sagelib (or on Sage at all) is src/sage/interfaces/primecount.pyx,
which is not even touched by this ticket (it's touched by #25009).

The latter indeed can be spun off, modelling on memory_allocator
(I suppose after it's done, one would only need to change one import statement in src/sage/functions/prime_pi.pyx.)

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. This was a forced push. Last 10 new commits:

2944718	check version of primesieve in the config of primecount
e8dbc40	distros info for primecount/sieve, correct m4s
c655136	remove sagemath-primecount declaration
9ca1f67	misc config changes for primesieve
d8e2cda	adjust deps,remove stuff from cfg and toml
3f1e441	no need for CMAKE_INSTALL_PREFIX (it's in sdh_cmake)
abce309	flyby patch: remove CMAKE_INSTALL_PREFIX, as it's set in sdh_cmake
a6de7f4	check external primesieve with cmake
733cebc	prime_pi and legendre_phi calling primecount
4f77073	allow for slightly different Overflow messages in doctest

[sagetrac-git]

Changed commit from 0c3be43 to 4f77073

[dimpase]

comment:41

we will proceed with this on #32894, which includes commits on this ticket's branch.

[mkoeppe]

Changed author from Dima Pasechnik to none