Search code examples
prologprimesswi-prologclpfd

CLPFD constraint: is a prime number


I'm not even sure if this is possible, but I'm trying to write a predicate prime/1 which constrains its argument to be a prime number.

The problem I have is that I haven't found any way of expressing “apply that constraint to all integers less than the variable integer”.

Here is an attempt which doesn't work:

prime(N) :-
    N #> 1 #/\                       % Has to be strictly greater than 1
    (
        N #= 2                       % Can be 2
        #\/                          % Or
        (
            N #> 2 #/\               % A number strictly greater than 2
            N mod 2 #= 1 #/\         % which is odd
            K #< N #/\
            K #> 1 #/\
            (#\ (
                N mod K #= 0         % A non working attempt at expressing:
                                         “there is no 1 < K < N such that K divides N”
            ))
        )
    ).

I hoped that #\ would act like \+ and check that it is false for all possible cases but this doesn't seem to be the case, since this implementation does this:

?- X #< 100, prime(X), indomain(X).
X = 2 ;    % Correct
X = 3 ;    % Correct
X = 5 ;    % Correct
X = 7 ;    % Correct
X = 9 ;    % Incorrect ; multiple of 3
X = 11 ;   % Correct
X = 13 ;   % Correct
X = 15     % Incorrect ; multiple of 5
…

Basically this unifies with 2\/{Odd integers greater than 2}.

EDIT

Expressing that a number is not prime is very easy:

composite(N) :-
    I #>= J,
    J #> 1,
    N #= I*J.

Basically: “N is composite if it can be written as I*J with I >= J > 1”.

I am still unable to “negate” those constraints. I have tried using things like #==> (implies) but this doesn't seem to be implification at all! N #= I*J #==> J #= 1 will work for composite numbers, even though 12 = I*J doesn't imply that necessarily J = 1!


Solution

  • prime/1

    This took me quite a while and I'm sure it's far from being very efficient but this seems to work, so here goes nothing:

    We create a custom constraint propagator (following this example) for the constraint prime/1, as such:

    :- use_module(library(clpfd)).
    :- multifile clpfd:run_propagator/2.
    
    prime(N) :-
        clpfd:make_propagator(prime(N), Prop),
        clpfd:init_propagator(N, Prop),
        clpfd:trigger_once(Prop).
    
    clpfd:run_propagator(prime(N), MState) :-
        (
            nonvar(N) -> clpfd:kill(MState), prime_decomposition(N, [_])
            ;
            clpfd:fd_get(N, ND, NL, NU, NPs),
            clpfd:cis_max(NL, n(2), NNL),
            clpfd:update_bounds(N, ND, NPs, NL, NU, NNL, NU)
        ).
    

    If N is a variable, we constrain its lower bound to be 2, or keep its original lower bound if it is bigger than 2.

    If N is ground, then we check that N is prime, using this prime_decomposition/2 predicate:

    prime_decomposition(2, [2]).
    prime_decomposition(N, Z) :-
        N #> 0,
        indomain(N),
        SN is ceiling(sqrt(N)),
        prime_decomposition_1(N, SN, 2, [], Z).
    
    prime_decomposition_1(1, _, _, L, L) :- !.
    prime_decomposition_1(N, SN, D, L, LF) :-
        (   
            0 #= N mod D -> !, false
            ;
            D1 #= D+1,
            (    
                D1 #> SN ->
                LF = [N |L]
                ;
                prime_decomposition_2(N, SN, D1, L, LF)
            )
        ).
    
    prime_decomposition_2(1, _, _, L, L) :- !.
    prime_decomposition_2(N, SN, D, L, LF) :-
        (   
            0 #= N mod D -> !, false
            ;
            D1 #= D+2,
            (    
                D1 #> SN ->
                LF = [N |L]
                ;
                prime_decomposition_2(N, SN, D1, L, LF)
            )
        ).
    

    You could obviously replace this predicate with any deterministic prime checking algorithm. This one is a modification of a prime factorization algorithm which has been modified to fail as soon as one factor is found.

    Some queries

    ?- prime(X).
    X in 2..sup,
    prime(X).
    
    
    ?- X in -100..100, prime(X).
    X in 2..100,
    prime(X).
    
    
    ?- X in -100..0, prime(X). 
    false.
    
    
    ?- X in 100..200, prime(X).
    X in 100..200,
    prime(X).
    
    
    ?- X #< 20, prime(X), indomain(X).
    X = 2 ;
    X = 3 ;
    X = 5 ;
    X = 7 ;
    X = 11 ;
    X = 13 ;
    X = 17 ;
    X = 19.
    
    ?- prime(X), prime(Y), [X, Y] ins 123456789..1234567890, Y-X #= 2, indomain(Y).
    X = 123457127,
    Y = 123457129 ;
    X = 123457289,
    Y = 123457291 ;
    X = 123457967,
    Y = 123457969
    …
    
    
    ?- time((X in 123456787654321..1234567876543210, prime(X), indomain(X))).
    % 113,041,584 inferences, 5.070 CPU in 5.063 seconds (100% CPU, 22296027 Lips)
    X = 123456787654391 .
    

    Some problems

    This constraint does not propagate as strongly as it should. For example:

    ?- prime(X), X in {2,3,8,16}.
    X in 2..3\/8\/16,
    prime(X).
    

    when we should know that 8 and 16 are not possible since they are even numbers.

    I have tried to add other constraints in the propagator but they seem to slow it down more than anything else, so I'm not sure if I was doing something wrong or if it is slower to update constaints than check for primeness when labeling.