search

oscar-system / Oscar.jl

A comprehensive open source computer algebra system for computations in algebra, geometry, and number theory.
https://www.oscar-system.org
Other
366 stars 129 forks source link

Use primesieve for prime enumeration #1455

Open AndrewVSutherland opened 2 years ago

AndrewVSutherland commented 2 years ago

The Julia Primes package which provides a Julia wrapper for Kim Walisch's primesieve library is faster than the prime enumeration functions that Oscar currently provides.

I would suggest either incorporating Primes.jl into Oscar or modifying functions like PrimesSet to use primesieve when appropriate (whenever the bounds will fit in an Int64).

Some performance comparisons on my system (single threaded):

julia> @time length(PrimesSet(1,10^9))
 53.839917 seconds (22 allocations: 489.726 MiB, 0.68% gc time)
50847534

julia> @time length(Primes.primes(1,10^9))
  3.293949 seconds (7 allocations: 643.503 MiB, 1.26% gc time)
50847534

julia> @time length(PrimesSet(10^15,10^15+10^9))
 84.298921 seconds (20 allocations: 229.086 MiB, 0.06% gc time)
28946421

julia> @time length(Primes.primes(10^15,10^15+10^9))
  4.292900 seconds (9 allocations: 483.821 MiB, 0.49% gc time)
28946421
fieker commented 2 years ago

We will look into it, but note that the "normal" use of PrimesSet is an unbounded set, ie. no upper bound. Its used in some modular methods to iterate over primes congruent to ?? mod ?? until enough have been found. It should be improved.

On Mon, Jul 11, 2022 at 07:31:54AM -0700, Andrew Sutherland wrote:

The Julia Primes package which provides a Julia wrapper for Kim Walisch's primesieve library is faster than the prime enumeration functions that Oscar currently provides.

I would suggest either incorporating Primes.jl into Oscar or modifying functions like PrimesSet to use primesieve when appropriate (whenever the bounds will fit in an Int64).

Some performance comparisons on my system (single threaded):

julia> @time length(PrimesSet(1,10^9))
 53.839917 seconds (22 allocations: 489.726 MiB, 0.68% gc time)
50847534

julia> @time length(Primes.primes(1,10^9))
  3.293949 seconds (7 allocations: 643.503 MiB, 1.26% gc time)
50847534

julia> @time length(PrimesSet(10^15,10^15+10^9))
 84.298921 seconds (20 allocations: 229.086 MiB, 0.06% gc time)
28946421

julia> @time length(Primes.primes(10^15,10^15+10^9))
  4.292900 seconds (9 allocations: 483.821 MiB, 0.49% gc time)
28946421

-- Reply to this email directly or view it on GitHub: https://github.com/oscar-system/Oscar.jl/issues/1455 You are receiving this because you are subscribed to this thread.

Message ID: @.***>

AndrewVSutherland commented 2 years ago

Yes, I'm sure the primesieve interface doesn't support everything you might want PrimesSet to do (in particular, handling integers larger than 2^63), but having a 15-20x faster implementation for the most common use case (primes up to X or within a known interval below 2^63) without needing to write any new prime enumeration code yourself seems like an easy win. Of course people can just use Primes.jl directly, but it would be nice not to add another dependency (and to be able to call PrimesSet with confidence in every situation).

But I predict that whatever you do will not beat using primesieve for enumerating all primes up to X < 2^63. I recently ripped all the prime enumeration code out of smalljac and replaced it with calls to primesive because despite having put a fair amount of effort into it (including using multiple wheels) it wasn't even close to primesieve's speed (even though it was much faster than Magma, Sage, and Pari/GP).