Sieve of Atkin in Erlang with a mutable HiPE bitarray. Forgive me, I was bored.

atkin.erl

-module(atkin).

-export([primes_smaller_than/1]).

primes_smaller_than(Limit) ->
    primes_smaller_than(Limit, sieve(Limit), []).

primes_smaller_than(0, _, Acc) ->
    Acc;
primes_smaller_than(N, Bin, Acc) ->
    case hipe_bifs:bitarray_sub(Bin, N) of
        true -> primes_smaller_than(N - 1, Bin, [N|Acc]);
        false -> primes_smaller_than(N - 1, Bin, Acc)
    end.

sieve(Limit) when Limit > 3 ->
    Bin = init(Limit + 1),
    Bound = trunc(math:sqrt(Limit)) + 1,
    put_candidates(Limit, Bin, Bound, 1, 1),
    elim_composites(Limit, Bin, Bound, 5).

init(Limit) ->
    Bin = hipe_bifs:bitarray(Limit, false),
    hipe_bifs:bitarray_update(Bin, 2, true),
    hipe_bifs:bitarray_update(Bin, 3, true).

put_candidates(_, Bin, Bound, Bound, Bound) ->
    Bin;
put_candidates(Limit, Bin, Bound, X, Bound) ->
    put_candidates(Limit, Bin, Bound, X + 1, 1);
put_candidates(Limit, Bin, Bound, X, Y) ->
    X2 = X * X,
    Y2 = Y * Y,
    Candidates =
        [ N || N <- [4 * X2 + Y2],
               N =< Limit,
               N rem 12 =:= 1 orelse N rem 12 =:= 5 ] ++
        [ N || N <- [3 * X2 + Y2],
               N =< Limit,
               N rem 12 =:= 7 ] ++
        [ N || X > Y,
               N <- [3 * X2 - Y2],
               N =< Limit,
               N rem 12 =:= 11 ],
    lists:foldl(fun (C, B) ->
                        Prev = hipe_bifs:bitarray_sub(B, C),
                        hipe_bifs:bitarray_update(B, C, not Prev)
                end, Bin, Candidates),
    put_candidates(Limit, Bin, Bound, X, Y + 1).

elim_composites(_, Bin, Bound, N) when N >= Bound ->
    Bin;
elim_composites(Limit, Bin, Bound, N) ->
    case hipe_bifs:bitarray_sub(Bin, N) of
        true ->
            N2 = N * N,
            elim_composites(Limit, Bin, Bound, N, N2, N2);
        false ->
            elim_composites(Limit, Bin, Bound, N + 1)
    end.

elim_composites(Limit, Bin, Bound, N, _, K) when K > Limit ->
    elim_composites(Limit, Bin, Bound, N + 1);
elim_composites(Limit, Bin, Bound, N, N2, K) ->
    hipe_bifs:bitarray_update(Bin, K, false),
    elim_composites(Limit, Bin, Bound, N, N2, K + N2).