prime-sieve.g

# suggested by Reimer Behrends
AllPrimesUpTo := function(n)
    local i, j, sieve, result;
    sieve := BlistList([1..n], [1..n]);
    sieve[1] := false;
    for i in [2..Int(n/2)] do
        sieve[i*2] := false;
    od;
    i := 3;
    while i * i <= n do
        if sieve[i] then
            j := 3*i;
            while j <= n do
                sieve[j] := false;
                j := j + 2*i;
            od;
        fi;
        i := i + 2;
    od;
    result := [];
    for i in [2..n] do
        if sieve[i] then
            Add(result, i);
        fi;
    od;
    return result;
end;