Prolog : CHR library for Project Euler

primality_test.swi

:- module(primality_test,[is_prime/1]).
:- use_module(library(chr)).
:- chr_option(optimize, full).
:- chr_type list(T) ---> []; [T|list(T)].
:- chr_constraint
	is_prime(+int),
	is_prime(+list(int),+int,+int),
	is_sub(+int,+int,+int),
	min_odd(+int,-int).

%%	Miller Rabin primality test

is_prime(2) <=> true.
is_prime(N) <=> N < 2 | false.
is_prime(N) <=> N rem 2 =:= 0 | false.
is_prime(N) <=> N < 4759123141 |
                N1 is N-1, min_odd(N1,T), is_prime([2,7,61],N,T).
is_prime(N) <=> N < 341550071728321 |
                N1 is N-1, min_odd(N1,T), is_prime([2,3,5,7,11,13,17],N,T).

is_prime([H|R],N,T) <=> H < N | Y is powm(H,T,N), is_sub(N,T,Y), is_prime(R,N,T).
is_prime(_,_,_)     <=> true.

is_sub(N,T,Y) <=> T =\= N-1, Y \= 1, Y =\= N-1 |
		  T1 is T*2, Y1 is (Y*Y) rem N, is_sub(N,T1,Y1).
is_sub(N,T,Y) <=> Y =\= N-1, T rem 2 =:= 0 | false.
is_sub(_,_,_) <=> true.

min_odd(N,X) <=> N rem 2 =:= 0 | N1 is N//2, min_odd(N1,X).
min_odd(N,X) <=> X = N.

prime_assoc.swi

:- module(prime_assoc,[prime_assoc/2]).
:- dynamic(prime/1).

prime_assoc(N,Assoc) :-
	set_prime(N), sieve(N),
	findall(P-true, prime(P), L),
	ord_list_to_assoc(L, Assoc),
	retractall(prime(_)).

set_prime(Max)   :- assertz(prime(2)), set_prime(3,Max).
set_prime(N,Max) :- N > Max, !.
set_prime(N,Max) :- assertz(prime(N)), N1 is N+2, set_prime(N1,Max).

sieve(Max) :- sieve(3,3,Max).
sieve(M,_,Max) :- M > Max, !.
sieve(M,N,Max) :- M*N > Max, !, M1 is M+2, sieve(M1,M1,Max).
sieve(M,N,Max) :- P is M*N, retractall(prime(P)), N1 is N+2, sieve(M,N1,Max).

prime_factors.swi

:- module(prime_factors,[prime_factors/2,count_factors/2,sum_factors/2]).
:- use_module(library(chr)).
:- chr_option(optimize, full).
:- chr_type list(T) ---> []; [T|list(T)].
:- chr_type pair    ---> int-int.
:- chr_constraint
        prime_factors(+int), num(+int), fact(+int), prime(+int,+int),
        list(+list(pair)), count, sum, mul(+int),
        return/1.

prime_factors(N,X) :- prime_factors(N), list([]), return(X).
count_factors(N,X) :- prime_factors(N), count,    return(X).
sum_factors(N,X)   :- prime_factors(N), sum,      return(X).

return(X), list(L) <=> X = L.
return(X), mul(N)  <=> X = N.

prime_factors(N) <=> num(N), fact(2).

fact(M), num(N)  <=> M*M > N | prime(N,1).
fact(M) \ num(N) <=> N rem M =:= 0 | prime(M,1), N1 is N//M, num(N1).
fact(2) <=> fact(3).
fact(N) <=> N1 is N+2, fact(N1).

prime(M,N), prime(M,N1) <=> X is N+N1, prime(M,X).

list(L), prime(M,N) # passive <=> list([M-N|L]).

count \ prime(_,N) # passive <=> N1 is N+1, mul(N1).
count <=> true.

sum \ prime(M,N) # passive
<=> numlist(0,N,L), maplist(pow(M),L,L1), sumlist(L1,X), mul(X).
sum <=> true.

mul(M), mul(N) <=> X is M*N, mul(X).