Skip to content

Instantly share code, notes, and snippets.

@nox nox/atkin.erl
Last active Dec 20, 2015

Embed
What would you like to do?
Sieve of Atkin in Erlang with a mutable HiPE bitarray. Forgive me, I was bored.
-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).
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.