How should we make this compatible wth current functionality, so that users don't see a bunch of method overwritten warnings?

Some options:

• Make the package 0.5 only: this will be a problem if we want to extend the package beyond the functionality available in Base (which was part of the reason of making it a package...)
• Don't export methods on 0.4, e.g. so that a user on 0.4 would need to explicitly import Primes: isprime.

Combinatorics.jl has trod this path already (though it was already a package beforehand); @jiahao any thoughts here?

Today primes([lo,] hi) returns a collection, an iterator looks more suitable.

To me, it would make sense if a range could be an argument to primes by
primes(r::UnitRange) = primes(r.start, r.stop)
to have things such as

> primes(1:10)
4-element Array{Int64,1}:
 2
 3
 5
 7

Problem Statement
primes(0) returns an error.

Desired Behavior
primes(0) returns a 0-element array like primes(0,0) and primes(0,1) do.

Example

julia> primes(0)
ERROR: ArgumentError: The condition lo ≤ hi must be met.
Stacktrace:
 [1] primes(::Int64, ::Int64) at C:\Users\amgough\.juliapro\JuliaPro_v1.4.0-1\packages\Primes\uaYlp\src\Primes.jl:111
 [2] primes(::Int64) at C:\Users\amgough\.juliapro\JuliaPro_v1.4.0-1\packages\Primes\uaYlp\src\Primes.jl:126
 [3] top-level scope at none:0

julia> primes(0,0)
0-element Array{Int64,1}

julia> primes(0,1)
0-element Array{Int64,1}

Julia version and Primes.jl version:

julia> versioninfo()
Julia Version 1.4.0
Commit b8e9a9ecc6 (2020-03-21 16:36 UTC)
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: Intel(R) Core(TM) i7-4960X CPU @ 3.60GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-8.0.1 (ORCJIT, ivybridge)

julia> Pkg.status("LinearAlgebra")
 [37e2e46d] LinearAlgebra

The docstrings for these aren't terribly long, so maybe the easiest thing to do is just copy them into the readme. There's also Documenter.jl but that may be overkill here.

idea comes from https://en.algorithmica.org/hpc/algorithms/factorization/. Lenstra will still be better for numbers larger than Int64, but getting extra speed is always nice.

The n-th prime function of julia is actually really slow.
For example in julia prime(10^7) runs in 14 seconds in my machine.

There is a well known C++ library with a 2-clause BSD license called primecount which can calculate
the prime(10^14) in 2 seconds (of course prime(10^7) is instantaneous.

Mathematica's Prime function is comparable to primecount (it is actually 2 times slower).

So given that, why not use primecount as a backend for the n-th prime function?

I could help with a patch if something like this sounds interesting.

using Primes;
factor(-52150815751994411270420247094986245419003171173880)

takes forever to compute. Any idea?

https://twitter.com/fermatslibrary/status/1364580859202965512?s=20

julia> using Primes
julia> isprime(82818079787776757473727170696867666564636261605958575655545352515049484746454443424140393837363534333231302928272625242322212019181716151413121110987654321)
true

or

julia> ss=""
""
julia> for x = 82:-1:1; ss = ss * string(x); end
julia> ii = parse(BigInt, ss)
82818079787776757473727170696867666564636261605958575655545352515049484746454443424140393837363534333231302928272625242322212019181716151413121110987654321
julia> isprime(ii)
true

Enhancement Request

Add a function for the next nsmooth number as suggested here. Suggested code is given. Note I have never made a pull request, so am intimidated by the thought.

I have several packages (CyclotomicNumbers, Combinat, Gapjm) where heavy use of factor is made -- mostly for other number-theoretic functions, like totient, divisors, prime_residues, primitive roots. The main application is various algorithms in finite group theory.
I need to call many times factor for small integers (for my applications small is <300 but this can be discussed of course).
I thus implemented caching:

const dict_factor=Dict(2=>Primes.factor(2))
factor(n::Integer)=if n<300 get!(()->Primes.factor(n),dict_factor,n) else Primes.factor(n) end

This way factor.(1:300) goes from 40μs to 3μs which is already better for my applications. But this redefinition causes problems when reexporting factor (conflict with Primes) so I wish this caching was in Primes.

The advantage of caching with a Dict is that only factorizations I need are stored (usually numbers with many factors which appears as orders of groups like 24, 120, etc...).

If just caching up to 300, it seems to me caching with a precomputed array is possible: this would just add 40μs to the precompilation time of Primes, and would bring the timing of factor.(1:300) to 500ns.

Anyway, I feel that Primes would be really more useful with such a caching, so this is my feature request.

Right now the docstring just says:

"""
    isprime(n::Integer) -> Bool

Returns `true` if `n` is prime, and `false` otherwise.

```julia
julia> isprime(3)
true

"""

But that's misleading because for large inputs, the result is only probabilistic.

Conversely, it *seems* to me that for $n < 2^{64}$ the test is proven to be accurate, i.e., deterministic. That would also be good to know and state explicitly.

I would provide a PR if that's desirable and if someone would confirm my guess above.

https://juliamath.github.io/Primes.jl/stable/ still is for version 0.5.0

This is because there is no working CI workflow set up for updating the website

The following PRs were never merged into Base, and so weren't copied over.

• JuliaLang/julia#16314
• JuliaLang/julia#16333
• JuliaLang/julia#16349

cc: @mschauer, @pabloferz

Hey!

primes(lo::Int, hi::Int) fails when lo <= 0 and hi >= 7

Probably a check in wheel_index(n) for n <= 1 to return 0 would work, assuming no negative primes.
Or, disallowing lo <= 0 as done in primesmask(lo::Int, hi::Int).

Here are the demos:

julia> Primes.primes(0, 6)
3-element Array{Int64,1}:
 2
 3
 5

julia> Primes.primes(0, 7)
ERROR: BoundsError: attempt to access 31-element Array{Int64,1}:
 0
 0
 0
 0
 0
 0
 0
 1
 1
 1
 1
 2
 2
 3
 3
 3
 3
 4
 4
 5
 5
 5
 5
 6
 6
 6
 6
 6
 6
 7
 7
  at index [0]
 in primes at /Users/myrmidon/.julia/v0.4/Primes/src/Primes.jl:124

A very small issue: in the doc of radical it is rightly stated radical(n::Integer). But in the source the function is radical(n) without any type restriction. It would be good to put one, to enable other packages to overload radical for other types with less clashes possible.

Implement a prime iterator. This was discussed in #8 . There was also a related question on stackoverflow indicating interest in such a thing, https://stackoverflow.com/q/37142821/2556061

At times, it is desirable to incorporate the default REPL output of factor(n) into strings for use, for example, within a println command. It would be nice to have a type variant - factor(String,n) - which outputs in string form (eg/ factor(String,90) = "2 * 3^2 * 5"), as well as a method for converting a Factorization form into a nice string (eg/ prodform(factor(90)) = "2 * 3^2 * 5").

I noticed something:

julia> h = Dict{Int,Int}(1=>10, 4=>6)
Dict{Int64,Int64} with 2 entries:
  4 => 6
  1 => 10

julia> factor(100, h)  # Undocumented usage
Dict{Int64,Int64} with 4 entries:
  4 => 6
  2 => 2
  5 => 2
  1 => 10

Despite the fact that such usage is undocumented, a user can still easily see that the method exists with methods(factor). Here's what I propose we do: follow the convention of using ! at the end of the function name, e.g. factor!, then define factor to call factor! on an empty dictionary, export factor, and not export factor!. That approach seems to be common in other packages. Thoughts?

cc @rfourquet

Ref: https://discourse.julialang.org/t/package-compatibility-caps/15301

Here is a example of two number theory functions made fast by eachfactor

"""
`prime_residues(n)` the numbers less than `n` and prime to `n`

julia> [prime_residues(24)]
1-element Vector{Vector{Int64}}:
 [1, 5, 7, 11, 13, 17, 19, 23]
"""
function prime_residues(n)
  if n==1 return [0] end
  pp=trues(n) # use a sieve to go fast
  for (p,np) in eachfactor(n)
    pp[p:p:n].=false
  end
  (1:n)[pp]
end

"""
`divisors(n)` the increasing list of divisors of `n`.

julia> [divisors(24)]
1-element Vector{Vector{Int64}}:
 [1, 2, 3, 4, 6, 8, 12, 24]
"""
function divisors(n::Integer)::Vector{Int}
  sort(vec(map(prod,Iterators.product((p.^(0:m) for (p,m) in eachfactor(n))...))))
end

I am considering adding some basic prime-factorization-related functions to 'Primes.jl', in particular some multiplicative functions (in the sense of number theory, i.e, such that f(a*b) = f(a)f(b) if a and b are relative prime).

Primes.jl already has totient (but it does not appear to be documented; I might fix this). I would add at least the Möbius function, the divisor function, and the Liouville function. Some of these functions I have used to help solve Project Euler problems. In Python, such functions are provided by SmyPy .

Similar to the totient function I would add two versions each, to have efficiency in case the
factorization is already known, e.g.:

moebiusmu(f::Factorization{T}) where T <: Integer
moebiusmu(n::Integer)

Also, I would add divisors(n::Integer) and divisors(f::Factorization{T}) which return an array (or an iterator) giving all the positive divisors of n.

What do you think about this?

Using Primes v0.5.0 on Julia Version 1.7.0-DEV.1108 (2021-05-16) we get:

$ julia
               _
   _       _ _(_)_     |  Documentation: https://docs.julialang.org
  (_)     | (_) (_)    |
   _ _   _| |_  __ _   |  Type "?" for help, "]?" for Pkg help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 1.7.0-DEV.1108 (2021-05-16)
 _/ |\__'_|_|_|\__'_|  |  Commit b1a4129338 (0 days old master)
|__/                   |

julia> using Primes

julia> nextprime(5)
5

julia> nextprime(97)
97

For non-primes it works fine:

julia> nextprime(8)
11

In certain situations, it is nice to have a totient function. I have a naive implementation below, but I'm not very familiar with the theory and so perhaps someone else can suggest a better one.

function totient(n::Integer)
    if iszero(n)
        throw(ArgumentError("ϕ(0) is not defined"))
    end
    result = one(n)
    for (p, k) in factor(abs(n))
        result *= p^(k-1) * (p - 1)
    end
    result
end

playing around with Factorizations I noticed a bug

julia> print(factor(24))
Primes.Factorization(2 => 3, 3 => 1)
julia> Primes.Factorization(2 => 3, 3 => 1)
ERROR: MethodError: no method matching Primes.Factorization(::Pair{Int64, Int64}, ::Pair{Int64, Int64})
Stacktrace:
 [1] top-level scope
   @ REPL[3]:1

you should either print factor(24) as Primes.Factorization(Dict(2 => 3, 3 => 1)) or add a method

Primes.Factorization(l::Pair{T,Int}...) where T=Primes.Factorization(Dict(l))

so that one can round-trip a Factorization

In Julia 1.9.2 using Primes v0.5.4 and SaferIntegers v3.4.0, the code

using Primes, SaferIntegers
factor(SafeInt128(8))

returns

ERROR: MethodError: no method matching increment!(::Primes.Factorization{SafeInt128}, ::SafeInt64, ::SafeInt128)

Closest candidates are:
  increment!(::Primes.Factorization{T}, ::Int64, ::Any) where T
   @ Primes ~/.julia/packages/Primes/sH6Tw/src/factorization.jl:48
  increment!(::AbstractDict, ::Int64, ::Any)
   @ Primes ~/.julia/packages/Primes/sH6Tw/src/factorization.jl:59

Stacktrace:
 [1] factor!(n::SafeInt128, h::Primes.Factorization{SafeInt128})
   @ Primes ~/.julia/packages/Primes/sH6Tw/src/Primes.jl:394
 [2] factor(n::SafeInt128)
   @ Primes ~/.julia/packages/Primes/sH6Tw/src/Primes.jl:429
 [3] top-level scope
   @ REPL[3]:1

Consider the following code.

julia> eltype(eachfactor(Int32(901800900)))
Tuple{Int32, Int64}

julia> for (p, e) in eachfactor(Int32(901800900))
       println()
       println((p, typeof(p)))
       println((e, typeof(e)))
       end

(2, Int32)
(2, Int64)

(3, Int64)
(2, Int64)

(5, Int64)
(2, Int64)

(7, Int64)
(2, Int64)

(11, Int32)
(2, Int64)

(13, Int32)
(2, Int64)

The eachfactor iterator should return (at least based on it's eltype) Tuple{Int32,Int64} but it returns a mixture of Tuple{Int32,Int64} and Tuple{Int64,Int64}.

I have written code for factoring a univariate polynomial over a finite field, and over the integers (modules Fact and FFfac in my package Gapjm). I thought it would be convenient to use Primes.Factorization to store these factorizations. However, I cannot since it is written in several places T<:Integer.

Are these annotations really useful? It seems to me they are not...

It would be cool if Factorization could be a representation of an abstract integer so that one could do arithmetics with them, etc.

f1e2 = factor(100)
f1e4 = f100 * f100

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

julia> nextprime(Int256(2)^186)
ERROR: InexactError: trunc(Int256, 776096500322445173294029667436835618525483108491085984157372499045210168189608446442688888866620606620268602202040084040440848044843)

We're special casing Int128 to use the BigInt method. Perhaps we should do that more generally for larger integer types?

There hasn't been one since #93 so factor is really slow still.

The second call doesn't look fine. The arguments are close to typemax(Int).

julia> primes(2^62-100, 2^62-1)
2-element Array{Int64,1}:
 4611686018427387817
 4611686018427387847

julia> primes(2^63-100, 2^63-1)
ERROR: BoundsError: attempt to access 31-element Array{Int64,1} at index [-26]
Stacktrace:
 [1] getindex at ./array.jl:758 [inlined]
 [2] wheel_index at /home/crusius/.julia/packages/Primes/uaYlp/src/Primes.jl:26 [inlined]
 [3] _primesmask(::Int64, ::Int64) at /home/crusius/.julia/packages/Primes/uaYlp/src/Primes.jl:68
 [4] primes(::Int64, ::Int64) at /home/crusius/.julia/packages/Primes/uaYlp/src/Primes.jl:119
 [5] top-level scope at REPL[91]:1

@simonbyrne mentioned in #13 that if we were doing this from scratch, he would have suggested primefactor in lieu of factor as it's more immediately descriptive. While we may not be doing things from scratch, we are in a situation in which we could make changes as needed. As such, I figured we should discuss this a bit.

I see three potential courses of action for factor in particular:

1. Leave it as-is and continue to export it
2. Deprecate it in favor of another name, e.g. primefactor
3. Leave it unexported, so one always calls Primes.factor (which itself is descriptive)

Surely there are other options as well. What does everyone think?

Moved from JuliaLang/julia#16780

The function isprimes of #16349, or more precisely of #9 with type hint for the witnesses observed a significant drop in speed recently, see the discussion at #9. When I was writing the comment JuliaLang/julia#16349 (comment) (using a slightly outdated master) no allocation took place in the inner loops at all. This seems to have changed recently as well and allocation takes place.

See this discourse comment.

It seems like this could be fixed with a little refactoring. Instead of having witnesses(n) return a variably sized tuple (which is not type stable), it would be better to inline this into isprime and do something like:

return n < 4294967296  ? checkwitnesses(n, _witnesses(UInt64(n))) :
        n < 2152302898747   ? checkwitnesses(n, (2, 3, 5, 7, 11)) :
        n < 3474749660383   ? checkwitnesses(n, (2, 3, 5, 7, 11, 13)):
                              checkwitnesses(n, (2, 325, 9375, 28178, 450775, 9780504, 1795265022))

where checkwitnesses(n, witnesses) is simply this loop refactored into a function.

Hello.
I believe that after finding every consecutive prime factor, performing isprime on resultant is harmful for performance of such function. Consider such snippet:

number = reduce(*, 1, Primes.PRIMES)
tic(); factor(number); toc()
# already running for 15 minutes, still didn't halt

Compare it with equivalent Python:

from operator import mul
from sympy import primerange, factorint
number = reduce(mul, primerange(1, 2**16), 1)
%timeit factorint(number)
# 1 loops, best of 3: 16.8 s per loop

Considering the fact that sympy is written in pure Python, this is pathetic.

Profiling Julia version brings such result (this time taking the product of first 3000 primes):

WARNING: The profile data buffer is full; profiling probably terminated
before your program finished. To profile for longer runs, call Profile.init
with a larger buffer and/or larger delay.
18816 ./event.jl:68; (::Base.REPL.##3#4{Base.REPL.REPLBackend})()
 18816 ./REPL.jl:95; macro expansion
  18816 ./REPL.jl:64; eval_user_input(::Any, ::Base.REPL.REPLBackend)
   18816 ./<missing>:?; anonymous
    18816 ./profile.jl:16; macro expansion;
     18816 /home/enedil/.julia/v0.5/Primes/src/Primes.jl:312; factor(::BigInt)
      1     ./dict.jl:701; factor!(::BigInt, ::Dict{BigInt,Int64})
      11    /home/enedil/.julia/v0.5/Primes/src/Primes.jl:267; factor!(::BigInt, ::Dict{BigInt,Int64})
       11 ./gmp.jl:110; convert
      3     /home/enedil/.julia/v0.5/Primes/src/Primes.jl:268; factor!(::BigInt, ::Dict{BigInt,Int64})
       1 ./gmp.jl:255; rem
       2 ./gmp.jl:256; rem
      2     /home/enedil/.julia/v0.5/Primes/src/Primes.jl:269; factor!(::BigInt, ::Dict{BigInt,Int64})
       1 ./dict.jl:644; setindex!(::Dict{BigInt,Int64}, ::Int64, ::BigInt)
      3     /home/enedil/.julia/v0.5/Primes/src/Primes.jl:270; factor!(::BigInt, ::Dict{BigInt,Int64})
       3 ./gmp.jl:256; div
      6     /home/enedil/.julia/v0.5/Primes/src/Primes.jl:271; factor!(::BigInt, ::Dict{BigInt,Int64})
       6 ./gmp.jl:255; rem
      18789 /home/enedil/.julia/v0.5/Primes/src/Primes.jl:276; factor!(::BigInt, ::Dict{BigInt,Int64})
       18789 /home/enedil/.julia/v0.5/Primes/src/Primes.jl:190; isprime
        18789 /home/enedil/.julia/v0.5/Primes/src/Primes.jl:190; isprime

Ha! I suspected that checking is number is prime after excluding one particular factor is superfluous, and I know it's basically one of the worst cases for such design, but after all, how such difference in time (with sympy) could be explained?

I propose the call to isprime to be supressed.

I recently set up documentation for this package using Documenter.jl. That worked well, though there are a couple of odd issues:

• All functions are listed with the prefix Base. This is likely because, even as of Julia 0.6-dev, isdefined(Base, :isprime) is true. Not sure what we can do about this.
• The "latest" page works, but "stable" does not; it's a 404. Is this because there hasn't been a tagged version since Documenter has been added?

Currently, memory usage of the prime sieve increases linearly with the size of the interval [lo, hi], while the number of primes only grows asymptotically as (hi) / log (hi). Furthermore, L1 memory caching is not exploited.

By implementing a segmented sieve both problems can be addressed at once: lowering memory consumption and increasing speed (in C++ I noticed up to 10x 5x better performance).

The current algorithm is basically:

1. Discover prime numbers p less than √hi
2. Cross off multiples of p starting at p * p.
3. Collect all prime numbers in [√hi, hi]

The suggested procedure is:

1. Discover all prime numbers p less than √hi
2. Divide [√hi, hi] into fixed size segments / smaller intervals (fitting L1 cache)
3. Outer loop: iterate over each segment
   • Inner loop: cross off all multiples of all primes in the current segment
   • Collect all primes in the current segment

You only need to store one segment (constant memory allocation) in step 2/3 and if the segment fits in L1 cache it is extremely fast, since you're iterating over it multiple times when sieving.

For reference see http://primesieve.org/

E.g.

julia> factor(20)
2^2 ⋅ 5

Before, the ⋅ operator (\cdot) was defined in Base (equivalent to dot), so, the expression 2^2 ⋅ 5 could be pasted back in the terminal and produce the input number, ie. 20. Now it works only if LinearAlgebra is imported. I see few solutions:

1. document that LinearAlgebra must be imported for that to work
2. change ⋅ to *
3. change ⋅ to * in the show method when we can detect that LinearAlgebra is not loaded
4. have the show(factor(20)) use * and keep ⋅ for the display at the REPL (suggested by @HarrisonGrodin, if I understood correctly)
5. have Primes depend on LinearAlgebra (I hope no-one will support this one!)

I favor 1. and then 2.. I don't like 3 because the output depends on global state, and I find 4. to be slighly overkill. But after all, the star * is not so ugly, so 2. may be the most pragmatic solution. Other ideas?

So this would create a factorization of 10*20 ...

myfact = factor(10)
Primes.factor!(20, myfact)

Could this line ...

be changed to:

isprime(n) && (h[n] = get(h, n, 0) + 1; return h) 

?

for p a prime, nextprime(p) == p

this should be changed so that for i>0 an integer nextprime(i) > i

thanks to leiteiro, here for noticing