Add a prime sieve algorithm

[ssiccha]

We want to find primes in a given range [min .. max] or simply in [1 .. max] somewhat efficiently. Such an algorithm is not implemented in GAP yet.

For a first draft of an implementation see this gist (provided by @rbehrends).

Todo:

• implementation
• documentation
• tests

[RussWoodroofe]

Comment that the best algorithms for this are segmented sieves: they divide [1..n] into buckets of size sqrt(n), and use the fact that a non-prime has a divisor smaller than its square root. If you're just going to print the results and not store them, this reduced memory usage to O(sqrt(n)). Maybe the right thing for memory usage would be to implement a collection of primes? (Where the collection would store the primes in [1..sqrt(n)], and use that to quickly calculate the primes in [1..n].)

[Stefan-Kohl]

Finding primes is one thing; sometimes it suffices to count the primes in a certain interval, which can be done notably faster -- actually in time less than linear in the number of primes to be counted. Several algorithms for this are implemented in the package 'primecount' by Kim Walisch -- see https://github.com/kimwalisch/primecount

[fingolfin]

While there are all kinds of advanced things one could do, I really think for starters a completely basic stock prime sieve would be sufficient, esp. for the numbers @ssiccha needs it for (which are probably all well under 1 million).

[frankluebeck]

Stumbling over this discussion I remembered to have two functions which I used in an exercise class in number theory:

# returns for positive integer n a Blist of length n/2-1 or (n+1)/2-1 where 
# the i-th entry is true if and only if 2*i+1 is a prime
PrimeBlist := function(n)
  local n2, l, falses, i, p, s, e;
  if n mod 2 = 1 then
    n := n+1;
  fi;
  n2 := n/2-1;
  l := BlistList([1..n2],[1..n2]);
  falses := BlistList([1..QuoInt(n2,3)+1],[]);
  i := 0;
  while (2*i+1)^2 <= n do
    i := Position(l, true, i);
    p := 2*i+1;
    s := (p*p-1)/2;
    e := n - (n mod p);
    if e mod 2 = 0 then
      e := e-p;
    fi;
    e := (e-1)/2;
    l{[s, s+p..e]} := falses{[1..(e-s)/p+1]};
  od;
  return l;
end;
# list of primes <= n, where odd primes are read off from result 
# of previous function
PrimesUpTo := function(n)
  local l;
  l := PrimeBlist(n);
  return Concatenation([2], Positions(l, true)*2+1);
end;

The result of the first function is more memory efficient. PrimesUpTo(10^9); takes about 16 seconds on my computer.
Feel free to put the code somewhere in the library.

Of course, one could write a similar function which computes the primes in an interval [k..n] using the result of PrimesUpTo(RootInt(n,2));.