View source with raw comments or as raw
    1:- encoding(utf8).
    2
    3/*  Part of SWI-Prolog
    4
    5    Author:        Markus Triska
    6    E-mail:        triska@metalevel.at
    7    WWW:           http://www.swi-prolog.org
    8    Copyright (C): 2007-2017 Markus Triska
    9    All rights reserved.
   10
   11    Redistribution and use in source and binary forms, with or without
   12    modification, are permitted provided that the following conditions
   13    are met:
   14
   15    1. Redistributions of source code must retain the above copyright
   16       notice, this list of conditions and the following disclaimer.
   17
   18    2. Redistributions in binary form must reproduce the above copyright
   19       notice, this list of conditions and the following disclaimer in
   20       the documentation and/or other materials provided with the
   21       distribution.
   22
   23    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
   24    "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
   25    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
   26    FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
   27    COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
   28    INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
   29    BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
   30    LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
   31    CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
   32    LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
   33    ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
   34    POSSIBILITY OF SUCH DAMAGE.
   35*/
   36
   37/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
   38   Thanks to Tom Schrijvers for his "bounds.pl", the first finite
   39   domain constraint solver included with SWI-Prolog. I've learned a
   40   lot from it and could even use some of the code for this solver.
   41   The propagation queue idea is taken from "prop.pl", a prototype
   42   solver also written by Tom. Highlights of the present solver:
   43
   44   Symbolic constants for infinities
   45   ---------------------------------
   46
   47   ?- X #>= 0, Y #=< 0.
   48   %@ X in 0..sup,
   49   %@ Y in inf..0.
   50
   51   No artificial limits (using GMP)
   52   ---------------------------------
   53
   54   ?- N #= 2^66, X #\= N.
   55   %@ N = 73786976294838206464,
   56   %@ X in inf..73786976294838206463\/73786976294838206465..sup.
   57
   58   Often stronger propagation
   59   ---------------------------------
   60
   61   ?- Y #= abs(X), Y #\= 3, Z * Z #= 4.
   62   %@ Y in 0..2\/4..sup,
   63   %@ Y#=abs(X),
   64   %@ X in inf.. -4\/ -2..2\/4..sup,
   65   %@ Z in -2\/2.
   66
   67   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
   68
   69   Development of this library has moved to SICStus Prolog. If you
   70   need any additional features or want to help, please file an issue at:
   71
   72                    https://github.com/triska/clpz
   73                    ==============================
   74
   75- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
   76
   77:- module(clpfd, [
   78                  op(760, yfx, #<==>),
   79                  op(750, xfy, #==>),
   80                  op(750, yfx, #<==),
   81                  op(740, yfx, #\/),
   82                  op(730, yfx, #\),
   83                  op(720, yfx, #/\),
   84                  op(710,  fy, #\),
   85                  op(700, xfx, #>),
   86                  op(700, xfx, #<),
   87                  op(700, xfx, #>=),
   88                  op(700, xfx, #=<),
   89                  op(700, xfx, #=),
   90                  op(700, xfx, #\=),
   91                  op(700, xfx, in),
   92                  op(700, xfx, ins),
   93                  op(450, xfx, ..), % should bind more tightly than \/
   94                  (#>)/2,
   95                  (#<)/2,
   96                  (#>=)/2,
   97                  (#=<)/2,
   98                  (#=)/2,
   99                  (#\=)/2,
  100                  (#\)/1,
  101                  (#<==>)/2,
  102                  (#==>)/2,
  103                  (#<==)/2,
  104                  (#\/)/2,
  105                  (#\)/2,
  106                  (#/\)/2,
  107                  (in)/2,
  108                  (ins)/2,
  109                  all_different/1,
  110                  all_distinct/1,
  111                  sum/3,
  112                  scalar_product/4,
  113                  tuples_in/2,
  114                  labeling/2,
  115                  label/1,
  116                  indomain/1,
  117                  lex_chain/1,
  118                  serialized/2,
  119                  global_cardinality/2,
  120                  global_cardinality/3,
  121                  circuit/1,
  122                  cumulative/1,
  123                  cumulative/2,
  124                  disjoint2/1,
  125                  element/3,
  126                  automaton/3,
  127                  automaton/8,
  128                  transpose/2,
  129                  zcompare/3,
  130                  chain/2,
  131                  fd_var/1,
  132                  fd_inf/2,
  133                  fd_sup/2,
  134                  fd_size/2,
  135                  fd_dom/2
  136                 ]).  137
  138:- public                               % called from goal_expansion
  139        clpfd_equal/2,
  140        clpfd_geq/2.  141
  142:- use_module(library(apply)).  143:- use_module(library(apply_macros)).  144:- use_module(library(assoc)).  145:- use_module(library(error)).  146:- use_module(library(lists)).  147:- use_module(library(pairs)).  148
  149:- op(700, xfx, cis).  150:- op(700, xfx, cis_geq).  151:- op(700, xfx, cis_gt).  152:- op(700, xfx, cis_leq).  153:- op(700, xfx, cis_lt).

CLP(FD): Constraint Logic Programming over Finite Domains

Development of this library has moved to SICStus Prolog.

Please see CLP(Z) for more information.

Introduction

This library provides CLP(FD): Constraint Logic Programming over Finite Domains. This is an instance of the general CLP(X) scheme, extending logic programming with reasoning over specialised domains. CLP(FD) lets us reason about integers in a way that honors the relational nature of Prolog.

Read The Power of Prolog to understand how this library is meant to be used in practice.

There are two major use cases of CLP(FD) constraints:

  1. declarative integer arithmetic
  2. solving combinatorial problems such as planning, scheduling and allocation tasks.

The predicates of this library can be classified as:

In most cases, arithmetic constraints are the only predicates you will ever need from this library. When reasoning over integers, simply replace low-level arithmetic predicates like (is)/2 and (>)/2 by the corresponding CLP(FD) constraints like #=/2 and #>/2 to honor and preserve declarative properties of your programs. For satisfactory performance, arithmetic constraints are implicitly rewritten at compilation time so that low-level fallback predicates are automatically used whenever possible.

Almost all Prolog programs also reason about integers. Therefore, it is highly advisable that you make CLP(FD) constraints available in all your programs. One way to do this is to put the following directive in your ~/.swiplrc initialisation file:

:- use_module(library(clpfd)).

All example programs that appear in the CLP(FD) documentation assume that you have done this.

Important concepts and principles of this library are illustrated by means of usage examples that are available in a public git repository: github.com/triska/clpfd

If you are used to the complicated operational considerations that low-level arithmetic primitives necessitate, then moving to CLP(FD) constraints may, due to their power and convenience, at first feel to you excessive and almost like cheating. It isn't. Constraints are an integral part of all popular Prolog systems, and they are designed to help you eliminate and avoid the use of low-level and less general primitives by providing declarative alternatives that are meant to be used instead.

When teaching Prolog, CLP(FD) constraints should be introduced before explaining low-level arithmetic predicates and their procedural idiosyncrasies. This is because constraints are easy to explain, understand and use due to their purely relational nature. In contrast, the modedness and directionality of low-level arithmetic primitives are impure limitations that are better deferred to more advanced lectures.

We recommend the following reference (PDF: metalevel.at/swiclpfd.pdf) for citing this library in scientific publications:

@inproceedings{Triska12,
  author    = {Markus Triska},
  title     = {The Finite Domain Constraint Solver of {SWI-Prolog}},
  booktitle = {FLOPS},
  series    = {LNCS},
  volume    = {7294},
  year      = {2012},
  pages     = {307-316}
}

More information about CLP(FD) constraints and their implementation is contained in: metalevel.at/drt.pdf

The best way to discuss applying, improving and extending CLP(FD) constraints is to use the dedicated clpfd tag on stackoverflow.com. Several of the world's foremost CLP(FD) experts regularly participate in these discussions and will help you for free on this platform.

Arithmetic constraints

In modern Prolog systems, arithmetic constraints subsume and supersede low-level predicates over integers. The main advantage of arithmetic constraints is that they are true relations and can be used in all directions. For most programs, arithmetic constraints are the only predicates you will ever need from this library.

The most important arithmetic constraint is #=/2, which subsumes both (is)/2 and (=:=)/2 over integers. Use #=/2 to make your programs more general. See declarative integer arithmetic.

In total, the arithmetic constraints are:

Expr1 #= Expr2Expr1 equals Expr2
Expr1 #\= Expr2Expr1 is not equal to Expr2
Expr1 #>= Expr2Expr1 is greater than or equal to Expr2
Expr1 #=< Expr2Expr1 is less than or equal to Expr2
Expr1 #> Expr2Expr1 is greater than Expr2
Expr1 #< Expr2Expr1 is less than Expr2

Expr1 and Expr2 denote arithmetic expressions, which are:

integerGiven value
variableUnknown integer
?(variable)Unknown integer
-ExprUnary minus
Expr + ExprAddition
Expr * ExprMultiplication
Expr - ExprSubtraction
Expr ^ ExprExponentiation
min(Expr,Expr)Minimum of two expressions
max(Expr,Expr)Maximum of two expressions
Expr mod ExprModulo induced by floored division
Expr rem ExprModulo induced by truncated division
abs(Expr)Absolute value
Expr // ExprTruncated integer division
Expr div ExprFloored integer division

where Expr again denotes an arithmetic expression.

The bitwise operations (\)/1, (/\)/2, (\/)/2, (>>)/2, (<<)/2, lsb/1, msb/1, popcount/1 and (xor)/2 are also supported.

Declarative integer arithmetic

The arithmetic constraints #=/2, #>/2 etc. are meant to be used instead of the primitives (is)/2, (=:=)/2, (>)/2 etc. over integers. Almost all Prolog programs also reason about integers. Therefore, it is recommended that you put the following directive in your ~/.swiplrc initialisation file to make CLP(FD) constraints available in all your programs:

:- use_module(library(clpfd)).

Throughout the following, it is assumed that you have done this.

The most basic use of CLP(FD) constraints is evaluation of arithmetic expressions involving integers. For example:

?- X #= 1+2.
X = 3.

This could in principle also be achieved with the lower-level predicate (is)/2. However, an important advantage of arithmetic constraints is their purely relational nature: Constraints can be used in all directions, also if one or more of their arguments are only partially instantiated. For example:

?- 3 #= Y+2.
Y = 1.

This relational nature makes CLP(FD) constraints easy to explain and use, and well suited for beginners and experienced Prolog programmers alike. In contrast, when using low-level integer arithmetic, we get:

?- 3 is Y+2.
ERROR: is/2: Arguments are not sufficiently instantiated

?- 3 =:= Y+2.
ERROR: =:=/2: Arguments are not sufficiently instantiated

Due to the necessary operational considerations, the use of these low-level arithmetic predicates is considerably harder to understand and should therefore be deferred to more advanced lectures.

For supported expressions, CLP(FD) constraints are drop-in replacements of these low-level arithmetic predicates, often yielding more general programs. See n_factorial/2 for an example.

This library uses goal_expansion/2 to automatically rewrite constraints at compilation time so that low-level arithmetic predicates are automatically used whenever possible. For example, the predicate:

positive_integer(N) :- N #>= 1.

is executed as if it were written as:

positive_integer(N) :-
        (   integer(N)
        ->  N >= 1
        ;   N #>= 1
        ).

This illustrates why the performance of CLP(FD) constraints is almost always completely satisfactory when they are used in modes that can be handled by low-level arithmetic. To disable the automatic rewriting, set the Prolog flag clpfd_goal_expansion to false.

If you are used to the complicated operational considerations that low-level arithmetic primitives necessitate, then moving to CLP(FD) constraints may, due to their power and convenience, at first feel to you excessive and almost like cheating. It isn't. Constraints are an integral part of all popular Prolog systems, and they are designed to help you eliminate and avoid the use of low-level and less general primitives by providing declarative alternatives that are meant to be used instead.

Example: Factorial relation

We illustrate the benefit of using #=/2 for more generality with a simple example.

Consider first a rather conventional definition of n_factorial/2, relating each natural number N to its factorial F:

n_factorial(0, 1).
n_factorial(N, F) :-
        N #> 0,
        N1 #= N - 1,
        n_factorial(N1, F1),
        F #= N * F1.

This program uses CLP(FD) constraints instead of low-level arithmetic throughout, and everything that would have worked with low-level arithmetic also works with CLP(FD) constraints, retaining roughly the same performance. For example:

?- n_factorial(47, F).
F = 258623241511168180642964355153611979969197632389120000000000 ;
false.

Now the point: Due to the increased flexibility and generality of CLP(FD) constraints, we are free to reorder the goals as follows:

n_factorial(0, 1).
n_factorial(N, F) :-
        N #> 0,
        N1 #= N - 1,
        F #= N * F1,
        n_factorial(N1, F1).

In this concrete case, termination properties of the predicate are improved. For example, the following queries now both terminate:

?- n_factorial(N, 1).
N = 0 ;
N = 1 ;
false.

?- n_factorial(N, 3).
false.

To make the predicate terminate if any argument is instantiated, add the (implied) constraint `F #\= 0` before the recursive call. Otherwise, the query n_factorial(N, 0) is the only non-terminating case of this kind.

The value of CLP(FD) constraints does not lie in completely freeing us from all procedural phenomena. For example, the two programs do not even have the same termination properties in all cases. Instead, the primary benefit of CLP(FD) constraints is that they allow you to try different execution orders and apply declarative debugging techniques at all! Reordering goals (and clauses) can significantly impact the performance of Prolog programs, and you are free to try different variants if you use declarative approaches. Moreover, since all CLP(FD) constraints always terminate, placing them earlier can at most improve, never worsen, the termination properties of your programs. An additional benefit of CLP(FD) constraints is that they eliminate the complexity of introducing (is)/2 and (=:=)/2 to beginners, since both predicates are subsumed by #=/2 when reasoning over integers.

In the case above, the clauses are mutually exclusive if the first argument is sufficiently instantiated. To make the predicate deterministic in such cases while retaining its generality, you can use zcompare/3 to reify a comparison, making the different cases distinguishable by pattern matching. For example, in this concrete case and others like it, you can use zcompare(Comp, 0, N) to obtain as Comp the symbolic outcome (<, =, >) of 0 compared to N.

Combinatorial constraints

In addition to subsuming and replacing low-level arithmetic predicates, CLP(FD) constraints are often used to solve combinatorial problems such as planning, scheduling and allocation tasks. Among the most frequently used combinatorial constraints are all_distinct/1, global_cardinality/2 and cumulative/2. This library also provides several other constraints like disjoint2/1 and automaton/8, which are useful in more specialized applications.

Domains

Each CLP(FD) variable has an associated set of admissible integers, which we call the variable's domain. Initially, the domain of each CLP(FD) variable is the set of all integers. CLP(FD) constraints like #=/2, #>/2 and #\=/2 can at most reduce, and never extend, the domains of their arguments. The constraints in/2 and ins/2 let us explicitly state domains of CLP(FD) variables. The process of determining and adjusting domains of variables is called constraint propagation, and it is performed automatically by this library. When the domain of a variable contains only one element, then the variable is automatically unified to that element.

Domains are taken into account when further constraints are stated, and by enumeration predicates like labeling/2.

Example: Sudoku

As another example, consider Sudoku: It is a popular puzzle over integers that can be easily solved with CLP(FD) constraints.

sudoku(Rows) :-
        length(Rows, 9), maplist(same_length(Rows), Rows),
        append(Rows, Vs), Vs ins 1..9,
        maplist(all_distinct, Rows),
        transpose(Rows, Columns),
        maplist(all_distinct, Columns),
        Rows = [As,Bs,Cs,Ds,Es,Fs,Gs,Hs,Is],
        blocks(As, Bs, Cs),
        blocks(Ds, Es, Fs),
        blocks(Gs, Hs, Is).

blocks([], [], []).
blocks([N1,N2,N3|Ns1], [N4,N5,N6|Ns2], [N7,N8,N9|Ns3]) :-
        all_distinct([N1,N2,N3,N4,N5,N6,N7,N8,N9]),
        blocks(Ns1, Ns2, Ns3).

problem(1, [[_,_,_,_,_,_,_,_,_],
            [_,_,_,_,_,3,_,8,5],
            [_,_,1,_,2,_,_,_,_],
            [_,_,_,5,_,7,_,_,_],
            [_,_,4,_,_,_,1,_,_],
            [_,9,_,_,_,_,_,_,_],
            [5,_,_,_,_,_,_,7,3],
            [_,_,2,_,1,_,_,_,_],
            [_,_,_,_,4,_,_,_,9]]).

Sample query:

?- problem(1, Rows), sudoku(Rows), maplist(writeln, Rows).
[9,8,7,6,5,4,3,2,1]
[2,4,6,1,7,3,9,8,5]
[3,5,1,9,2,8,7,4,6]
[1,2,8,5,3,7,6,9,4]
[6,3,4,8,9,2,1,5,7]
[7,9,5,4,6,1,8,3,2]
[5,1,9,2,8,6,4,7,3]
[4,7,2,3,1,9,5,6,8]
[8,6,3,7,4,5,2,1,9]
Rows = [[9, 8, 7, 6, 5, 4, 3, 2|...], ... , [...|...]].

In this concrete case, the constraint solver is strong enough to find the unique solution without any search. For the general case, see search.

Residual goals

Here is an example session with a few queries and their answers:

?- X #> 3.
X in 4..sup.

?- X #\= 20.
X in inf..19\/21..sup.

?- 2*X #= 10.
X = 5.

?- X*X #= 144.
X in -12\/12.

?- 4*X + 2*Y #= 24, X + Y #= 9, [X,Y] ins 0..sup.
X = 3,
Y = 6.

?- X #= Y #<==> B, X in 0..3, Y in 4..5.
B = 0,
X in 0..3,
Y in 4..5.

The answers emitted by the toplevel are called residual programs, and the goals that comprise each answer are called residual goals. In each case above, and as for all pure programs, the residual program is declaratively equivalent to the original query. From the residual goals, it is clear that the constraint solver has deduced additional domain restrictions in many cases.

To inspect residual goals, it is best to let the toplevel display them for us. Wrap the call of your predicate into call_residue_vars/2 to make sure that all constrained variables are displayed. To make the constraints a variable is involved in available as a Prolog term for further reasoning within your program, use copy_term/3. For example:

?- X #= Y + Z, X in 0..5, copy_term([X,Y,Z], [X,Y,Z], Gs).
Gs = [clpfd: (X in 0..5), clpfd: (Y+Z#=X)],
X in 0..5,
Y+Z#=X.

This library also provides reflection predicates (like fd_dom/2, fd_size/2 etc.) with which we can inspect a variable's current domain. These predicates can be useful if you want to implement your own labeling strategies.

Using CLP(FD) constraints to solve combinatorial tasks typically consists of two phases:

  1. Modeling. In this phase, all relevant constraints are stated.
  2. Search. In this phase, enumeration predicates are used to search for concrete solutions.

It is good practice to keep the modeling part, via a dedicated predicate called the core relation, separate from the actual search for solutions. This lets us observe termination and determinism properties of the core relation in isolation from the search, and more easily try different search strategies.

As an example of a constraint satisfaction problem, consider the cryptoarithmetic puzzle SEND + MORE = MONEY, where different letters denote distinct integers between 0 and 9. It can be modeled in CLP(FD) as follows:

puzzle([S,E,N,D] + [M,O,R,E] = [M,O,N,E,Y]) :-
        Vars = [S,E,N,D,M,O,R,Y],
        Vars ins 0..9,
        all_different(Vars),
                  S*1000 + E*100 + N*10 + D +
                  M*1000 + O*100 + R*10 + E #=
        M*10000 + O*1000 + N*100 + E*10 + Y,
        M #\= 0, S #\= 0.

Notice that we are not using labeling/2 in this predicate, so that we can first execute and observe the modeling part in isolation. Sample query and its result (actual variables replaced for readability):

?- puzzle(As+Bs=Cs).
As = [9, A2, A3, A4],
Bs = [1, 0, B3, A2],
Cs = [1, 0, A3, A2, C5],
A2 in 4..7,
all_different([9, A2, A3, A4, 1, 0, B3, C5]),
91*A2+A4+10*B3#=90*A3+C5,
A3 in 5..8,
A4 in 2..8,
B3 in 2..8,
C5 in 2..8.

From this answer, we see that this core relation terminates and is in fact deterministic. Moreover, we see from the residual goals that the constraint solver has deduced more stringent bounds for all variables. Such observations are only possible if modeling and search parts are cleanly separated.

Labeling can then be used to search for solutions in a separate predicate or goal:

?- puzzle(As+Bs=Cs), label(As).
As = [9, 5, 6, 7],
Bs = [1, 0, 8, 5],
Cs = [1, 0, 6, 5, 2] ;
false.

In this case, it suffices to label a subset of variables to find the puzzle's unique solution, since the constraint solver is strong enough to reduce the domains of remaining variables to singleton sets. In general though, it is necessary to label all variables to obtain ground solutions.

Example: Eight queens puzzle

We illustrate the concepts of the preceding sections by means of the so-called eight queens puzzle. The task is to place 8 queens on an 8x8 chessboard such that none of the queens is under attack. This means that no two queens share the same row, column or diagonal.

To express this puzzle via CLP(FD) constraints, we must first pick a suitable representation. Since CLP(FD) constraints reason over integers, we must find a way to map the positions of queens to integers. Several such mappings are conceivable, and it is not immediately obvious which we should use. On top of that, different constraints can be used to express the desired relations. For such reasons, modeling combinatorial problems via CLP(FD) constraints often necessitates some creativity and has been described as more of an art than a science.

In our concrete case, we observe that there must be exactly one queen per column. The following representation therefore suggests itself: We are looking for 8 integers, one for each column, where each integer denotes the row of the queen that is placed in the respective column, and which are subject to certain constraints.

In fact, let us now generalize the task to the so-called N queens puzzle, which is obtained by replacing 8 by N everywhere it occurs in the above description. We implement the above considerations in the core relation n_queens/2, where the first argument is the number of queens (which is identical to the number of rows and columns of the generalized chessboard), and the second argument is a list of N integers that represents a solution in the form described above.

n_queens(N, Qs) :-
        length(Qs, N),
        Qs ins 1..N,
        safe_queens(Qs).

safe_queens([]).
safe_queens([Q|Qs]) :- safe_queens(Qs, Q, 1), safe_queens(Qs).

safe_queens([], _, _).
safe_queens([Q|Qs], Q0, D0) :-
        Q0 #\= Q,
        abs(Q0 - Q) #\= D0,
        D1 #= D0 + 1,
        safe_queens(Qs, Q0, D1).

Note that all these predicates can be used in all directions: We can use them to find solutions, test solutions and complete partially instantiated solutions.

The original task can be readily solved with the following query:

?- n_queens(8, Qs), label(Qs).
Qs = [1, 5, 8, 6, 3, 7, 2, 4] .

Using suitable labeling strategies, we can easily find solutions with 80 queens and more:

?- n_queens(80, Qs), labeling([ff], Qs).
Qs = [1, 3, 5, 44, 42, 4, 50, 7, 68|...] .

?- time((n_queens(90, Qs), labeling([ff], Qs))).
% 5,904,401 inferences, 0.722 CPU in 0.737 seconds (98% CPU)
Qs = [1, 3, 5, 50, 42, 4, 49, 7, 59|...] .

Experimenting with different search strategies is easy because we have separated the core relation from the actual search.

Optimisation

We can use labeling/2 to minimize or maximize the value of a CLP(FD) expression, and generate solutions in increasing or decreasing order of the value. See the labeling options min(Expr) and max(Expr), respectively.

Again, to easily try different labeling options in connection with optimisation, we recommend to introduce a dedicated predicate for posting constraints, and to use labeling/2 in a separate goal. This way, we can observe properties of the core relation in isolation, and try different labeling options without recompiling our code.

If necessary, we can use once/1 to commit to the first optimal solution. However, it is often very valuable to see alternative solutions that are also optimal, so that we can choose among optimal solutions by other criteria. For the sake of purity and completeness, we recommend to avoid once/1 and other constructs that lead to impurities in CLP(FD) programs.

Related to optimisation with CLP(FD) constraints are library(simplex) and CLP(Q) which reason about linear constraints over rational numbers.

Reification

The constraints in/2, #=/2, #\=/2, #</2, #>/2, #=</2, and #>=/2 can be reified, which means reflecting their truth values into Boolean values represented by the integers 0 and 1. Let P and Q denote reifiable constraints or Boolean variables, then:

#\ QTrue iff Q is false
P #\/ QTrue iff either P or Q
P #/\ QTrue iff both P and Q
P #\ QTrue iff either P or Q, but not both
P #<==> QTrue iff P and Q are equivalent
P #==> QTrue iff P implies Q
P #<== QTrue iff Q implies P

The constraints of this table are reifiable as well.

When reasoning over Boolean variables, also consider using CLP(B) constraints as provided by library(clpb).

Enabling monotonic CLP(FD)

In the default execution mode, CLP(FD) constraints still exhibit some non-relational properties. For example, adding constraints can yield new solutions:

?-          X #= 2, X = 1+1.
false.

?- X = 1+1, X #= 2, X = 1+1.
X = 1+1.

This behaviour is highly problematic from a logical point of view, and it may render declarative debugging techniques inapplicable.

Set the Prolog flag clpfd_monotonic to true to make CLP(FD) monotonic: This means that adding new constraints cannot yield new solutions. When this flag is true, we must wrap variables that occur in arithmetic expressions with the functor (?)/1 or (#)/1. For example:

?- set_prolog_flag(clpfd_monotonic, true).
true.

?- #(X) #= #(Y) + #(Z).
#(Y)+ #(Z)#= #(X).

?-          X #= 2, X = 1+1.
ERROR: Arguments are not sufficiently instantiated

The wrapper can be omitted for variables that are already constrained to integers.

Custom constraints

We can define custom constraints. The mechanism to do this is not yet finalised, and we welcome suggestions and descriptions of use cases that are important to you.

As an example of how it can be done currently, let us define a new custom constraint oneground(X,Y,Z), where Z shall be 1 if at least one of X and Y is instantiated:

:- multifile clpfd:run_propagator/2.

oneground(X, Y, Z) :-
        clpfd:make_propagator(oneground(X, Y, Z), Prop),
        clpfd:init_propagator(X, Prop),
        clpfd:init_propagator(Y, Prop),
        clpfd:trigger_once(Prop).

clpfd:run_propagator(oneground(X, Y, Z), MState) :-
        (   integer(X) -> clpfd:kill(MState), Z = 1
        ;   integer(Y) -> clpfd:kill(MState), Z = 1
        ;   true
        ).

First, make_propagator/2 is used to transform a user-defined representation of the new constraint to an internal form. With init_propagator/2, this internal form is then attached to X and Y. From now on, the propagator will be invoked whenever the domains of X or Y are changed. Then, trigger_once/1 is used to give the propagator its first chance for propagation even though the variables' domains have not yet changed. Finally, run_propagator/2 is extended to define the actual propagator. As explained, this predicate is automatically called by the constraint solver. The first argument is the user-defined representation of the constraint as used in make_propagator/2, and the second argument is a mutable state that can be used to prevent further invocations of the propagator when the constraint has become entailed, by using kill/1. An example of using the new constraint:

?- oneground(X, Y, Z), Y = 5.
Y = 5,
Z = 1,
X in inf..sup.

Applications

CLP(FD) applications that we find particularly impressive and worth studying include:

Acknowledgments

This library gives you a glimpse of what SICStus Prolog can do. The API is intentionally mostly compatible with that of SICStus Prolog, so that you can easily switch to a much more feature-rich and much faster CLP(FD) system when you need it. I thank Mats Carlsson, the designer and main implementor of SICStus Prolog, for his elegant example. I first encountered his system as part of the excellent GUPU teaching environment by Ulrich Neumerkel. Ulrich was also the first and most determined tester of the present system, filing hundreds of comments and suggestions for improvement. Tom Schrijvers has contributed several constraint libraries to SWI-Prolog, and I learned a lot from his coding style and implementation examples. Bart Demoen was a driving force behind the implementation of attributed variables in SWI-Prolog, and this library could not even have started without his prior work and contributions. Thank you all!

CLP(FD) predicate index

In the following, each CLP(FD) predicate is described in more detail.

We recommend the following link to refer to this manual:

http://eu.swi-prolog.org/man/clpfd.html

author
- Markus Triska */
  929:- create_prolog_flag(clpfd_monotonic, false, []).  930
  931/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  932   A bound is either:
  933
  934   n(N):    integer N
  935   inf:     infimum of Z (= negative infinity)
  936   sup:     supremum of Z (= positive infinity)
  937- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  938
  939is_bound(n(N)) :- integer(N).
  940is_bound(inf).
  941is_bound(sup).
  942
  943defaulty_to_bound(D, P) :- ( integer(D) -> P = n(D) ; P = D ).
  944
  945/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  946   Compactified is/2 and predicates for several arithmetic expressions
  947   with infinities, tailored for the modes needed by this solver.
  948- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  949
  950goal_expansion(A cis B, Expansion) :-
  951        phrase(cis_goals(B, A), Goals),
  952        list_goal(Goals, Expansion).
  953goal_expansion(A cis_lt B, B cis_gt A).
  954goal_expansion(A cis_leq B, B cis_geq A).
  955goal_expansion(A cis_geq B, cis_leq_numeric(B, N)) :- nonvar(A), A = n(N).
  956goal_expansion(A cis_geq B, cis_geq_numeric(A, N)) :- nonvar(B), B = n(N).
  957goal_expansion(A cis_gt B, cis_lt_numeric(B, N))   :- nonvar(A), A = n(N).
  958goal_expansion(A cis_gt B, cis_gt_numeric(A, N))   :- nonvar(B), B = n(N).
  959
  960% cis_gt only works for terms of depth 0 on both sides
  961cis_gt(sup, B0) :- B0 \== sup.
  962cis_gt(n(N), B) :- cis_lt_numeric(B, N).
  963
  964cis_lt_numeric(inf, _).
  965cis_lt_numeric(n(B), A) :- B < A.
  966
  967cis_gt_numeric(sup, _).
  968cis_gt_numeric(n(B), A) :- B > A.
  969
  970cis_geq(inf, inf).
  971cis_geq(sup, _).
  972cis_geq(n(N), B) :- cis_leq_numeric(B, N).
  973
  974cis_leq_numeric(inf, _).
  975cis_leq_numeric(n(B), A) :- B =< A.
  976
  977cis_geq_numeric(sup, _).
  978cis_geq_numeric(n(B), A) :- B >= A.
  979
  980cis_min(inf, _, inf).
  981cis_min(sup, B, B).
  982cis_min(n(N), B, Min) :- cis_min_(B, N, Min).
  983
  984cis_min_(inf, _, inf).
  985cis_min_(sup, N, n(N)).
  986cis_min_(n(B), A, n(M)) :- M is min(A,B).
  987
  988cis_max(sup, _, sup).
  989cis_max(inf, B, B).
  990cis_max(n(N), B, Max) :- cis_max_(B, N, Max).
  991
  992cis_max_(inf, N, n(N)).
  993cis_max_(sup, _, sup).
  994cis_max_(n(B), A, n(M)) :- M is max(A,B).
  995
  996cis_plus(inf, _, inf).
  997cis_plus(sup, _, sup).
  998cis_plus(n(A), B, Plus) :- cis_plus_(B, A, Plus).
  999
 1000cis_plus_(sup, _, sup).
 1001cis_plus_(inf, _, inf).
 1002cis_plus_(n(B), A, n(S)) :- S is A + B.
 1003
 1004cis_minus(inf, _, inf).
 1005cis_minus(sup, _, sup).
 1006cis_minus(n(A), B, M) :- cis_minus_(B, A, M).
 1007
 1008cis_minus_(inf, _, sup).
 1009cis_minus_(sup, _, inf).
 1010cis_minus_(n(B), A, n(M)) :- M is A - B.
 1011
 1012cis_uminus(inf, sup).
 1013cis_uminus(sup, inf).
 1014cis_uminus(n(A), n(B)) :- B is -A.
 1015
 1016cis_abs(inf, sup).
 1017cis_abs(sup, sup).
 1018cis_abs(n(A), n(B)) :- B is abs(A).
 1019
 1020cis_times(inf, B, P) :-
 1021        (   B cis_lt n(0) -> P = sup
 1022        ;   B cis_gt n(0) -> P = inf
 1023        ;   P = n(0)
 1024        ).
 1025cis_times(sup, B, P) :-
 1026        (   B cis_gt n(0) -> P = sup
 1027        ;   B cis_lt n(0) -> P = inf
 1028        ;   P = n(0)
 1029        ).
 1030cis_times(n(N), B, P) :- cis_times_(B, N, P).
 1031
 1032cis_times_(inf, A, P)     :- cis_times(inf, n(A), P).
 1033cis_times_(sup, A, P)     :- cis_times(sup, n(A), P).
 1034cis_times_(n(B), A, n(P)) :- P is A * B.
 1035
 1036cis_exp(inf, n(Y), R) :-
 1037        (   even(Y) -> R = sup
 1038        ;   R = inf
 1039        ).
 1040cis_exp(sup, _, sup).
 1041cis_exp(n(N), Y, R) :- cis_exp_(Y, N, R).
 1042
 1043cis_exp_(n(Y), N, n(R)) :- R is N^Y.
 1044cis_exp_(sup, _, sup).
 1045cis_exp_(inf, _, inf).
 1046
 1047cis_goals(V, V)          --> { var(V) }, !.
 1048cis_goals(n(N), n(N))    --> [].
 1049cis_goals(inf, inf)      --> [].
 1050cis_goals(sup, sup)      --> [].
 1051cis_goals(sign(A0), R)   --> cis_goals(A0, A), [cis_sign(A, R)].
 1052cis_goals(abs(A0), R)    --> cis_goals(A0, A), [cis_abs(A, R)].
 1053cis_goals(-A0, R)        --> cis_goals(A0, A), [cis_uminus(A, R)].
 1054cis_goals(A0+B0, R)      -->
 1055        cis_goals(A0, A),
 1056        cis_goals(B0, B),
 1057        [cis_plus(A, B, R)].
 1058cis_goals(A0-B0, R)      -->
 1059        cis_goals(A0, A),
 1060        cis_goals(B0, B),
 1061        [cis_minus(A, B, R)].
 1062cis_goals(min(A0,B0), R) -->
 1063        cis_goals(A0, A),
 1064        cis_goals(B0, B),
 1065        [cis_min(A, B, R)].
 1066cis_goals(max(A0,B0), R) -->
 1067        cis_goals(A0, A),
 1068        cis_goals(B0, B),
 1069        [cis_max(A, B, R)].
 1070cis_goals(A0*B0, R)      -->
 1071        cis_goals(A0, A),
 1072        cis_goals(B0, B),
 1073        [cis_times(A, B, R)].
 1074cis_goals(div(A0,B0), R) -->
 1075        cis_goals(A0, A),
 1076        cis_goals(B0, B),
 1077        [cis_div(A, B, R)].
 1078cis_goals(A0//B0, R)     -->
 1079        cis_goals(A0, A),
 1080        cis_goals(B0, B),
 1081        [cis_slash(A, B, R)].
 1082cis_goals(A0^B0, R)      -->
 1083        cis_goals(A0, A),
 1084        cis_goals(B0, B),
 1085        [cis_exp(A, B, R)].
 1086
 1087list_goal([], true).
 1088list_goal([G|Gs], Goal) :- foldl(list_goal_, Gs, G, Goal).
 1089
 1090list_goal_(G, G0, (G0,G)).
 1091
 1092cis_sign(sup, n(1)).
 1093cis_sign(inf, n(-1)).
 1094cis_sign(n(N), n(S)) :- S is sign(N).
 1095
 1096cis_div(sup, Y, Z)  :- ( Y cis_geq n(0) -> Z = sup ; Z = inf ).
 1097cis_div(inf, Y, Z)  :- ( Y cis_geq n(0) -> Z = inf ; Z = sup ).
 1098cis_div(n(X), Y, Z) :- cis_div_(Y, X, Z).
 1099
 1100cis_div_(sup, _, n(0)).
 1101cis_div_(inf, _, n(0)).
 1102cis_div_(n(Y), X, Z) :-
 1103        (   Y =:= 0 -> (  X >= 0 -> Z = sup ; Z = inf )
 1104        ;   Z0 is X // Y, Z = n(Z0)
 1105        ).
 1106
 1107cis_slash(sup, _, sup).
 1108cis_slash(inf, _, inf).
 1109cis_slash(n(N), B, S) :- cis_slash_(B, N, S).
 1110
 1111cis_slash_(sup, _, n(0)).
 1112cis_slash_(inf, _, n(0)).
 1113cis_slash_(n(B), A, n(S)) :- S is A // B.
 1114
 1115
 1116/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1117   A domain is a finite set of disjoint intervals. Internally, domains
 1118   are represented as trees. Each node is one of:
 1119
 1120   empty: empty domain.
 1121
 1122   split(N, Left, Right)
 1123      - split on integer N, with Left and Right domains whose elements are
 1124        all less than and greater than N, respectively. The domain is the
 1125        union of Left and Right, i.e., N is a hole.
 1126
 1127   from_to(From, To)
 1128      - interval (From-1, To+1); From and To are bounds
 1129
 1130   Desiderata: rebalance domains; singleton intervals.
 1131- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1132
 1133/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1134   Type definition and inspection of domains.
 1135- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1136
 1137check_domain(D) :-
 1138        (   var(D) -> instantiation_error(D)
 1139        ;   is_domain(D) -> true
 1140        ;   domain_error(clpfd_domain, D)
 1141        ).
 1142
 1143is_domain(empty).
 1144is_domain(from_to(From,To)) :-
 1145        is_bound(From), is_bound(To),
 1146        From cis_leq To.
 1147is_domain(split(S, Left, Right)) :-
 1148        integer(S),
 1149        is_domain(Left), is_domain(Right),
 1150        all_less_than(Left, S),
 1151        all_greater_than(Right, S).
 1152
 1153all_less_than(empty, _).
 1154all_less_than(from_to(From,To), S) :-
 1155        From cis_lt n(S), To cis_lt n(S).
 1156all_less_than(split(S0,Left,Right), S) :-
 1157        S0 < S,
 1158        all_less_than(Left, S),
 1159        all_less_than(Right, S).
 1160
 1161all_greater_than(empty, _).
 1162all_greater_than(from_to(From,To), S) :-
 1163        From cis_gt n(S), To cis_gt n(S).
 1164all_greater_than(split(S0,Left,Right), S) :-
 1165        S0 > S,
 1166        all_greater_than(Left, S),
 1167        all_greater_than(Right, S).
 1168
 1169default_domain(from_to(inf,sup)).
 1170
 1171domain_infimum(from_to(I, _), I).
 1172domain_infimum(split(_, Left, _), I) :- domain_infimum(Left, I).
 1173
 1174domain_supremum(from_to(_, S), S).
 1175domain_supremum(split(_, _, Right), S) :- domain_supremum(Right, S).
 1176
 1177domain_num_elements(empty, n(0)).
 1178domain_num_elements(from_to(From,To), Num) :- Num cis To - From + n(1).
 1179domain_num_elements(split(_, Left, Right), Num) :-
 1180        domain_num_elements(Left, NL),
 1181        domain_num_elements(Right, NR),
 1182        Num cis NL + NR.
 1183
 1184domain_direction_element(from_to(n(From), n(To)), Dir, E) :-
 1185        (   Dir == up -> between(From, To, E)
 1186        ;   between(From, To, E0),
 1187            E is To - (E0 - From)
 1188        ).
 1189domain_direction_element(split(_, D1, D2), Dir, E) :-
 1190        (   Dir == up ->
 1191            (   domain_direction_element(D1, Dir, E)
 1192            ;   domain_direction_element(D2, Dir, E)
 1193            )
 1194        ;   (   domain_direction_element(D2, Dir, E)
 1195            ;   domain_direction_element(D1, Dir, E)
 1196            )
 1197        ).
 1198
 1199/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1200   Test whether domain contains a given integer.
 1201- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1202
 1203domain_contains(from_to(From,To), I) :- From cis_leq n(I), n(I) cis_leq To.
 1204domain_contains(split(S, Left, Right), I) :-
 1205        (   I < S -> domain_contains(Left, I)
 1206        ;   I > S -> domain_contains(Right, I)
 1207        ).
 1208
 1209/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1210   Test whether a domain contains another domain.
 1211- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1212
 1213domain_subdomain(Dom, Sub) :- domain_subdomain(Dom, Dom, Sub).
 1214
 1215domain_subdomain(from_to(_,_), Dom, Sub) :-
 1216        domain_subdomain_fromto(Sub, Dom).
 1217domain_subdomain(split(_, _, _), Dom, Sub) :-
 1218        domain_subdomain_split(Sub, Dom, Sub).
 1219
 1220domain_subdomain_split(empty, _, _).
 1221domain_subdomain_split(from_to(From,To), split(S,Left0,Right0), Sub) :-
 1222        (   To cis_lt n(S) -> domain_subdomain(Left0, Left0, Sub)
 1223        ;   From cis_gt n(S) -> domain_subdomain(Right0, Right0, Sub)
 1224        ).
 1225domain_subdomain_split(split(_,Left,Right), Dom, _) :-
 1226        domain_subdomain(Dom, Dom, Left),
 1227        domain_subdomain(Dom, Dom, Right).
 1228
 1229domain_subdomain_fromto(empty, _).
 1230domain_subdomain_fromto(from_to(From,To), from_to(From0,To0)) :-
 1231        From0 cis_leq From, To0 cis_geq To.
 1232domain_subdomain_fromto(split(_,Left,Right), Dom) :-
 1233        domain_subdomain_fromto(Left, Dom),
 1234        domain_subdomain_fromto(Right, Dom).
 1235
 1236/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1237   Remove an integer from a domain. The domain is traversed until an
 1238   interval is reached from which the element can be removed, or until
 1239   it is clear that no such interval exists.
 1240- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1241
 1242domain_remove(empty, _, empty).
 1243domain_remove(from_to(L0, U0), X, D) :- domain_remove_(L0, U0, X, D).
 1244domain_remove(split(S, Left0, Right0), X, D) :-
 1245        (   X =:= S -> D = split(S, Left0, Right0)
 1246        ;   X < S ->
 1247            domain_remove(Left0, X, Left1),
 1248            (   Left1 == empty -> D = Right0
 1249            ;   D = split(S, Left1, Right0)
 1250            )
 1251        ;   domain_remove(Right0, X, Right1),
 1252            (   Right1 == empty -> D = Left0
 1253            ;   D = split(S, Left0, Right1)
 1254            )
 1255        ).
 1256
 1257%?- domain_remove(from_to(n(0),n(5)), 3, D).
 1258
 1259domain_remove_(inf, U0, X, D) :-
 1260        (   U0 == n(X) -> U1 is X - 1, D = from_to(inf, n(U1))
 1261        ;   U0 cis_lt n(X) -> D = from_to(inf,U0)
 1262        ;   L1 is X + 1, U1 is X - 1,
 1263            D = split(X, from_to(inf, n(U1)), from_to(n(L1),U0))
 1264        ).
 1265domain_remove_(n(N), U0, X, D) :- domain_remove_upper(U0, N, X, D).
 1266
 1267domain_remove_upper(sup, L0, X, D) :-
 1268        (   L0 =:= X -> L1 is X + 1, D = from_to(n(L1),sup)
 1269        ;   L0 > X -> D = from_to(n(L0),sup)
 1270        ;   L1 is X + 1, U1 is X - 1,
 1271            D = split(X, from_to(n(L0),n(U1)), from_to(n(L1),sup))
 1272        ).
 1273domain_remove_upper(n(U0), L0, X, D) :-
 1274        (   L0 =:= U0, X =:= L0 -> D = empty
 1275        ;   L0 =:= X -> L1 is X + 1, D = from_to(n(L1), n(U0))
 1276        ;   U0 =:= X -> U1 is X - 1, D = from_to(n(L0), n(U1))
 1277        ;   between(L0, U0, X) ->
 1278            U1 is X - 1, L1 is X + 1,
 1279            D = split(X, from_to(n(L0), n(U1)), from_to(n(L1), n(U0)))
 1280        ;   D = from_to(n(L0),n(U0))
 1281        ).
 1282
 1283/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1284   Remove all elements greater than / less than a constant.
 1285- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1286
 1287domain_remove_greater_than(empty, _, empty).
 1288domain_remove_greater_than(from_to(From0,To0), G, D) :-
 1289        (   From0 cis_gt n(G) -> D = empty
 1290        ;   To cis min(To0,n(G)), D = from_to(From0,To)
 1291        ).
 1292domain_remove_greater_than(split(S,Left0,Right0), G, D) :-
 1293        (   S =< G ->
 1294            domain_remove_greater_than(Right0, G, Right),
 1295            (   Right == empty -> D = Left0
 1296            ;   D = split(S, Left0, Right)
 1297            )
 1298        ;   domain_remove_greater_than(Left0, G, D)
 1299        ).
 1300
 1301domain_remove_smaller_than(empty, _, empty).
 1302domain_remove_smaller_than(from_to(From0,To0), V, D) :-
 1303        (   To0 cis_lt n(V) -> D = empty
 1304        ;   From cis max(From0,n(V)), D = from_to(From,To0)
 1305        ).
 1306domain_remove_smaller_than(split(S,Left0,Right0), V, D) :-
 1307        (   S >= V ->
 1308            domain_remove_smaller_than(Left0, V, Left),
 1309            (   Left == empty -> D = Right0
 1310            ;   D = split(S, Left, Right0)
 1311            )
 1312        ;   domain_remove_smaller_than(Right0, V, D)
 1313        ).
 1314
 1315
 1316/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1317   Remove a whole domain from another domain. (Set difference.)
 1318- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1319
 1320domain_subtract(Dom0, Sub, Dom) :- domain_subtract(Dom0, Dom0, Sub, Dom).
 1321
 1322domain_subtract(empty, _, _, empty).
 1323domain_subtract(from_to(From0,To0), Dom, Sub, D) :-
 1324        (   Sub == empty -> D = Dom
 1325        ;   Sub = from_to(From,To) ->
 1326            (   From == To -> From = n(X), domain_remove(Dom, X, D)
 1327            ;   From cis_gt To0 -> D = Dom
 1328            ;   To cis_lt From0 -> D = Dom
 1329            ;   From cis_leq From0 ->
 1330                (   To cis_geq To0 -> D = empty
 1331                ;   From1 cis To + n(1),
 1332                    D = from_to(From1, To0)
 1333                )
 1334            ;   To1 cis From - n(1),
 1335                (   To cis_lt To0 ->
 1336                    From = n(S),
 1337                    From2 cis To + n(1),
 1338                    D = split(S,from_to(From0,To1),from_to(From2,To0))
 1339                ;   D = from_to(From0,To1)
 1340                )
 1341            )
 1342        ;   Sub = split(S, Left, Right) ->
 1343            (   n(S) cis_gt To0 -> domain_subtract(Dom, Dom, Left, D)
 1344            ;   n(S) cis_lt From0 -> domain_subtract(Dom, Dom, Right, D)
 1345            ;   domain_subtract(Dom, Dom, Left, D1),
 1346                domain_subtract(D1, D1, Right, D)
 1347            )
 1348        ).
 1349domain_subtract(split(S, Left0, Right0), _, Sub, D) :-
 1350        domain_subtract(Left0, Left0, Sub, Left),
 1351        domain_subtract(Right0, Right0, Sub, Right),
 1352        (   Left == empty -> D = Right
 1353        ;   Right == empty -> D = Left
 1354        ;   D = split(S, Left, Right)
 1355        ).
 1356
 1357/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1358   Complement of a domain
 1359- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1360
 1361domain_complement(D, C) :-
 1362        default_domain(Default),
 1363        domain_subtract(Default, D, C).
 1364
 1365/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1366   Convert domain to a list of disjoint intervals From-To.
 1367- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1368
 1369domain_intervals(D, Is) :- phrase(domain_intervals(D), Is).
 1370
 1371domain_intervals(split(_, Left, Right)) -->
 1372        domain_intervals(Left), domain_intervals(Right).
 1373domain_intervals(empty)                 --> [].
 1374domain_intervals(from_to(From,To))      --> [From-To].
 1375
 1376/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1377   To compute the intersection of two domains D1 and D2, we choose D1
 1378   as the reference domain. For each interval of D1, we compute how
 1379   far and to which values D2 lets us extend it.
 1380- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1381
 1382domains_intersection(D1, D2, Intersection) :-
 1383        domains_intersection_(D1, D2, Intersection),
 1384        Intersection \== empty.
 1385
 1386domains_intersection_(empty, _, empty).
 1387domains_intersection_(from_to(L0,U0), D2, Dom) :-
 1388        narrow(D2, L0, U0, Dom).
 1389domains_intersection_(split(S,Left0,Right0), D2, Dom) :-
 1390        domains_intersection_(Left0, D2, Left1),
 1391        domains_intersection_(Right0, D2, Right1),
 1392        (   Left1 == empty -> Dom = Right1
 1393        ;   Right1 == empty -> Dom = Left1
 1394        ;   Dom = split(S, Left1, Right1)
 1395        ).
 1396
 1397narrow(empty, _, _, empty).
 1398narrow(from_to(L0,U0), From0, To0, Dom) :-
 1399        From1 cis max(From0,L0), To1 cis min(To0,U0),
 1400        (   From1 cis_gt To1 -> Dom = empty
 1401        ;   Dom = from_to(From1,To1)
 1402        ).
 1403narrow(split(S, Left0, Right0), From0, To0, Dom) :-
 1404        (   To0 cis_lt n(S) -> narrow(Left0, From0, To0, Dom)
 1405        ;   From0 cis_gt n(S) -> narrow(Right0, From0, To0, Dom)
 1406        ;   narrow(Left0, From0, To0, Left1),
 1407            narrow(Right0, From0, To0, Right1),
 1408            (   Left1 == empty -> Dom = Right1
 1409            ;   Right1 == empty -> Dom = Left1
 1410            ;   Dom = split(S, Left1, Right1)
 1411            )
 1412        ).
 1413
 1414/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1415   Union of 2 domains.
 1416- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1417
 1418domains_union(D1, D2, Union) :-
 1419        domain_intervals(D1, Is1),
 1420        domain_intervals(D2, Is2),
 1421        append(Is1, Is2, IsU0),
 1422        merge_intervals(IsU0, IsU1),
 1423        intervals_to_domain(IsU1, Union).
 1424
 1425/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1426   Shift the domain by an offset.
 1427- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1428
 1429domain_shift(empty, _, empty).
 1430domain_shift(from_to(From0,To0), O, from_to(From,To)) :-
 1431        From cis From0 + n(O), To cis To0 + n(O).
 1432domain_shift(split(S0, Left0, Right0), O, split(S, Left, Right)) :-
 1433        S is S0 + O,
 1434        domain_shift(Left0, O, Left),
 1435        domain_shift(Right0, O, Right).
 1436
 1437/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1438   The new domain contains all values of the old domain,
 1439   multiplied by a constant multiplier.
 1440- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1441
 1442domain_expand(D0, M, D) :-
 1443        (   M < 0 ->
 1444            domain_negate(D0, D1),
 1445            M1 is abs(M),
 1446            domain_expand_(D1, M1, D)
 1447        ;   M =:= 1 -> D = D0
 1448        ;   domain_expand_(D0, M, D)
 1449        ).
 1450
 1451domain_expand_(empty, _, empty).
 1452domain_expand_(from_to(From0, To0), M, from_to(From,To)) :-
 1453        From cis From0*n(M),
 1454        To cis To0*n(M).
 1455domain_expand_(split(S0, Left0, Right0), M, split(S, Left, Right)) :-
 1456        S is M*S0,
 1457        domain_expand_(Left0, M, Left),
 1458        domain_expand_(Right0, M, Right).
 1459
 1460/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1461   similar to domain_expand/3, tailored for truncated division: an
 1462   interval [From,To] is extended to [From*M, ((To+1)*M - 1)], i.e.,
 1463   to all values that truncated integer-divided by M yield a value
 1464   from interval.
 1465- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1466
 1467domain_expand_more(D0, M, D) :-
 1468        %format("expanding ~w by ~w\n", [D0,M]),
 1469        (   M < 0 -> domain_negate(D0, D1), M1 is abs(M)
 1470        ;   D1 = D0, M1 = M
 1471        ),
 1472        domain_expand_more_(D1, M1, D).
 1473        %format("yield: ~w\n", [D]).
 1474
 1475domain_expand_more_(empty, _, empty).
 1476domain_expand_more_(from_to(From0, To0), M, from_to(From,To)) :-
 1477        (   From0 cis_leq n(0) ->
 1478            From cis (From0-n(1))*n(M) + n(1)
 1479        ;   From cis From0*n(M)
 1480        ),
 1481        (   To0 cis_lt n(0) ->
 1482            To cis To0*n(M)
 1483        ;   To cis (To0+n(1))*n(M) - n(1)
 1484        ).
 1485domain_expand_more_(split(S0, Left0, Right0), M, split(S, Left, Right)) :-
 1486        S is M*S0,
 1487        domain_expand_more_(Left0, M, Left),
 1488        domain_expand_more_(Right0, M, Right).
 1489
 1490/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1491   Scale a domain down by a constant multiplier. Assuming (//)/2.
 1492- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1493
 1494domain_contract(D0, M, D) :-
 1495        %format("contracting ~w by ~w\n", [D0,M]),
 1496        (   M < 0 -> domain_negate(D0, D1), M1 is abs(M)
 1497        ;   D1 = D0, M1 = M
 1498        ),
 1499        domain_contract_(D1, M1, D).
 1500
 1501domain_contract_(empty, _, empty).
 1502domain_contract_(from_to(From0, To0), M, from_to(From,To)) :-
 1503        (   From0 cis_geq n(0) ->
 1504            From cis (From0 + n(M) - n(1)) // n(M)
 1505        ;   From cis From0 // n(M)
 1506        ),
 1507        (   To0 cis_geq n(0) ->
 1508            To cis To0 // n(M)
 1509        ;   To cis (To0 - n(M) + n(1)) // n(M)
 1510        ).
 1511domain_contract_(split(_,Left0,Right0), M, D) :-
 1512        %  Scaled down domains do not necessarily retain any holes of
 1513        %  the original domain.
 1514        domain_contract_(Left0, M, Left),
 1515        domain_contract_(Right0, M, Right),
 1516        domains_union(Left, Right, D).
 1517
 1518/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1519   Similar to domain_contract, tailored for division, i.e.,
 1520   {21,23} contracted by 4 is 5. It contracts "less".
 1521- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1522
 1523domain_contract_less(D0, M, D) :-
 1524        (   M < 0 -> domain_negate(D0, D1), M1 is abs(M)
 1525        ;   D1 = D0, M1 = M
 1526        ),
 1527        domain_contract_less_(D1, M1, D).
 1528
 1529domain_contract_less_(empty, _, empty).
 1530domain_contract_less_(from_to(From0, To0), M, from_to(From,To)) :-
 1531        From cis From0 // n(M), To cis To0 // n(M).
 1532domain_contract_less_(split(_,Left0,Right0), M, D) :-
 1533        %  Scaled down domains do not necessarily retain any holes of
 1534        %  the original domain.
 1535        domain_contract_less_(Left0, M, Left),
 1536        domain_contract_less_(Right0, M, Right),
 1537        domains_union(Left, Right, D).
 1538
 1539/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1540   Negate the domain. Left and Right sub-domains and bounds switch sides.
 1541- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1542
 1543domain_negate(empty, empty).
 1544domain_negate(from_to(From0, To0), from_to(From, To)) :-
 1545        From cis -To0, To cis -From0.
 1546domain_negate(split(S0, Left0, Right0), split(S, Left, Right)) :-
 1547        S is -S0,
 1548        domain_negate(Left0, Right),
 1549        domain_negate(Right0, Left).
 1550
 1551/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1552   Construct a domain from a list of integers. Try to balance it.
 1553- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1554
 1555list_to_disjoint_intervals([], []).
 1556list_to_disjoint_intervals([N|Ns], Is) :-
 1557        list_to_disjoint_intervals(Ns, N, N, Is).
 1558
 1559list_to_disjoint_intervals([], M, N, [n(M)-n(N)]).
 1560list_to_disjoint_intervals([B|Bs], M, N, Is) :-
 1561        (   B =:= N + 1 ->
 1562            list_to_disjoint_intervals(Bs, M, B, Is)
 1563        ;   Is = [n(M)-n(N)|Rest],
 1564            list_to_disjoint_intervals(Bs, B, B, Rest)
 1565        ).
 1566
 1567list_to_domain(List0, D) :-
 1568        (   List0 == [] -> D = empty
 1569        ;   sort(List0, List),
 1570            list_to_disjoint_intervals(List, Is),
 1571            intervals_to_domain(Is, D)
 1572        ).
 1573
 1574intervals_to_domain([], empty) :- !.
 1575intervals_to_domain([M-N], from_to(M,N)) :- !.
 1576intervals_to_domain(Is, D) :-
 1577        length(Is, L),
 1578        FL is L // 2,
 1579        length(Front, FL),
 1580        append(Front, Tail, Is),
 1581        Tail = [n(Start)-_|_],
 1582        Hole is Start - 1,
 1583        intervals_to_domain(Front, Left),
 1584        intervals_to_domain(Tail, Right),
 1585        D = split(Hole, Left, Right).
 1586
 1587%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 in(?Var, +Domain)
Var is an element of Domain. Domain is one of:
Integer
Singleton set consisting only of Integer.
..(Lower, Upper)
All integers I such that Lower =< I =< Upper. Lower must be an integer or the atom inf, which denotes negative infinity. Upper must be an integer or the atom sup, which denotes positive infinity.
Domain1 \/ Domain2
The union of Domain1 and Domain2.
 1604Var in Dom :- clpfd_in(Var, Dom).
 1605
 1606clpfd_in(V, D) :-
 1607        fd_variable(V),
 1608        drep_to_domain(D, Dom),
 1609        domain(V, Dom).
 1610
 1611fd_variable(V) :-
 1612        (   var(V) -> true
 1613        ;   integer(V) -> true
 1614        ;   type_error(integer, V)
 1615        ).
 ins(+Vars, +Domain)
The variables in the list Vars are elements of Domain. See in/2 for the syntax of Domain.
 1622Vs ins D :-
 1623        fd_must_be_list(Vs),
 1624        maplist(fd_variable, Vs),
 1625        drep_to_domain(D, Dom),
 1626        domains(Vs, Dom).
 1627
 1628fd_must_be_list(Ls) :-
 1629        (   fd_var(Ls) -> type_error(list, Ls)
 1630        ;   must_be(list, Ls)
 1631        ).
 indomain(?Var)
Bind Var to all feasible values of its domain on backtracking. The domain of Var must be finite.
 1638indomain(Var) :- label([Var]).
 1639
 1640order_dom_next(up, Dom, Next)   :- domain_infimum(Dom, n(Next)).
 1641order_dom_next(down, Dom, Next) :- domain_supremum(Dom, n(Next)).
 1642order_dom_next(random_value(_), Dom, Next) :-
 1643        phrase(domain_to_intervals(Dom), Is),
 1644        length(Is, L),
 1645        R is random(L),
 1646        nth0(R, Is, From-To),
 1647        random_between(From, To, Next).
 1648
 1649domain_to_intervals(from_to(n(From),n(To))) --> [From-To].
 1650domain_to_intervals(split(_, Left, Right)) -->
 1651        domain_to_intervals(Left),
 1652        domain_to_intervals(Right).
 label(+Vars)
Equivalent to labeling([], Vars). See labeling/2.
 1658label(Vs) :- labeling([], Vs).
 labeling(+Options, +Vars)
Assign a value to each variable in Vars. Labeling means systematically trying out values for the finite domain variables Vars until all of them are ground. The domain of each variable in Vars must be finite. Options is a list of options that let you exhibit some control over the search process. Several categories of options exist:

The variable selection strategy lets you specify which variable of Vars is labeled next and is one of:

leftmost
Label the variables in the order they occur in Vars. This is the default.
ff
First fail. Label the leftmost variable with smallest domain next, in order to detect infeasibility early. This is often a good strategy.
ffc
Of the variables with smallest domains, the leftmost one participating in most constraints is labeled next.
min
Label the leftmost variable whose lower bound is the lowest next.
max
Label the leftmost variable whose upper bound is the highest next.

The value order is one of:

up
Try the elements of the chosen variable's domain in ascending order. This is the default.
down
Try the domain elements in descending order.

The branching strategy is one of:

step
For each variable X, a choice is made between X = V and X #\= V, where V is determined by the value ordering options. This is the default.
enum
For each variable X, a choice is made between X = V_1, X = V_2 etc., for all values V_i of the domain of X. The order is determined by the value ordering options.
bisect
For each variable X, a choice is made between X #=< M and X #> M, where M is the midpoint of the domain of X.

At most one option of each category can be specified, and an option must not occur repeatedly.

The order of solutions can be influenced with:

This generates solutions in ascending/descending order with respect to the evaluation of the arithmetic expression Expr. Labeling Vars must make Expr ground. If several such options are specified, they are interpreted from left to right, e.g.:

?- [X,Y] ins 10..20, labeling([max(X),min(Y)],[X,Y]).

This generates solutions in descending order of X, and for each binding of X, solutions are generated in ascending order of Y. To obtain the incomplete behaviour that other systems exhibit with "maximize(Expr)" and "minimize(Expr)", use once/1, e.g.:

once(labeling([max(Expr)], Vars))

Labeling is always complete, always terminates, and yields no redundant solutions. See core relations and search for usage advice.

 1745labeling(Options, Vars) :-
 1746        must_be(list, Options),
 1747        fd_must_be_list(Vars),
 1748        maplist(must_be_finite_fdvar, Vars),
 1749        label(Options, Options, default(leftmost), default(up), default(step), [], upto_ground, Vars).
 1750
 1751finite_domain(Dom) :-
 1752        domain_infimum(Dom, n(_)),
 1753        domain_supremum(Dom, n(_)).
 1754
 1755must_be_finite_fdvar(Var) :-
 1756        (   fd_get(Var, Dom, _) ->
 1757            (   finite_domain(Dom) -> true
 1758            ;   instantiation_error(Var)
 1759            )
 1760        ;   integer(Var) -> true
 1761        ;   must_be(integer, Var)
 1762        ).
 1763
 1764
 1765label([O|Os], Options, Selection, Order, Choice, Optim, Consistency, Vars) :-
 1766        (   var(O)-> instantiation_error(O)
 1767        ;   override(selection, Selection, O, Options, S1) ->
 1768            label(Os, Options, S1, Order, Choice, Optim, Consistency, Vars)
 1769        ;   override(order, Order, O, Options, O1) ->
 1770            label(Os, Options, Selection, O1, Choice, Optim, Consistency, Vars)
 1771        ;   override(choice, Choice, O, Options, C1) ->
 1772            label(Os, Options, Selection, Order, C1, Optim, Consistency, Vars)
 1773        ;   optimisation(O) ->
 1774            label(Os, Options, Selection, Order, Choice, [O|Optim], Consistency, Vars)
 1775        ;   consistency(O, O1) ->
 1776            label(Os, Options, Selection, Order, Choice, Optim, O1, Vars)
 1777        ;   domain_error(labeling_option, O)
 1778        ).
 1779label([], _, Selection, Order, Choice, Optim0, Consistency, Vars) :-
 1780        maplist(arg(1), [Selection,Order,Choice], [S,O,C]),
 1781        (   Optim0 == [] ->
 1782            label(Vars, S, O, C, Consistency)
 1783        ;   reverse(Optim0, Optim),
 1784            exprs_singlevars(Optim, SVs),
 1785            optimise(Vars, [S,O,C], SVs)
 1786        ).
 1787
 1788% Introduce new variables for each min/max expression to avoid
 1789% reparsing expressions during optimisation.
 1790
 1791exprs_singlevars([], []).
 1792exprs_singlevars([E|Es], [SV|SVs]) :-
 1793        E =.. [F,Expr],
 1794        ?(Single) #= Expr,
 1795        SV =.. [F,Single],
 1796        exprs_singlevars(Es, SVs).
 1797
 1798all_dead(fd_props(Bs,Gs,Os)) :-
 1799        all_dead_(Bs),
 1800        all_dead_(Gs),
 1801        all_dead_(Os).
 1802
 1803all_dead_([]).
 1804all_dead_([propagator(_, S)|Ps]) :- S == dead, all_dead_(Ps).
 1805
 1806label([], _, _, _, Consistency) :- !,
 1807        (   Consistency = upto_in(I0,I) -> I0 = I
 1808        ;   true
 1809        ).
 1810label(Vars, Selection, Order, Choice, Consistency) :-
 1811        (   Vars = [V|Vs], nonvar(V) -> label(Vs, Selection, Order, Choice, Consistency)
 1812        ;   select_var(Selection, Vars, Var, RVars),
 1813            (   var(Var) ->
 1814                (   Consistency = upto_in(I0,I), fd_get(Var, _, Ps), all_dead(Ps) ->
 1815                    fd_size(Var, Size),
 1816                    I1 is I0*Size,
 1817                    label(RVars, Selection, Order, Choice, upto_in(I1,I))
 1818                ;   Consistency = upto_in, fd_get(Var, _, Ps), all_dead(Ps) ->
 1819                    label(RVars, Selection, Order, Choice, Consistency)
 1820                ;   choice_order_variable(Choice, Order, Var, RVars, Vars, Selection, Consistency)
 1821                )
 1822            ;   label(RVars, Selection, Order, Choice, Consistency)
 1823            )
 1824        ).
 1825
 1826choice_order_variable(step, Order, Var, Vars, Vars0, Selection, Consistency) :-
 1827        fd_get(Var, Dom, _),
 1828        order_dom_next(Order, Dom, Next),
 1829        (   Var = Next,
 1830            label(Vars, Selection, Order, step, Consistency)
 1831        ;   neq_num(Var, Next),
 1832            do_queue,
 1833            label(Vars0, Selection, Order, step, Consistency)
 1834        ).
 1835choice_order_variable(enum, Order, Var, Vars, _, Selection, Consistency) :-
 1836        fd_get(Var, Dom0, _),
 1837        domain_direction_element(Dom0, Order, Var),
 1838        label(Vars, Selection, Order, enum, Consistency).
 1839choice_order_variable(bisect, Order, Var, _, Vars0, Selection, Consistency) :-
 1840        fd_get(Var, Dom, _),
 1841        domain_infimum(Dom, n(I)),
 1842        domain_supremum(Dom, n(S)),
 1843        Mid0 is (I + S) // 2,
 1844        (   Mid0 =:= S -> Mid is Mid0 - 1 ; Mid = Mid0 ),
 1845        (   Order == up -> ( Var #=< Mid ; Var #> Mid )
 1846        ;   Order == down -> ( Var #> Mid ; Var #=< Mid )
 1847        ;   domain_error(bisect_up_or_down, Order)
 1848        ),
 1849        label(Vars0, Selection, Order, bisect, Consistency).
 1850
 1851override(What, Prev, Value, Options, Result) :-
 1852        call(What, Value),
 1853        override_(Prev, Value, Options, Result).
 1854
 1855override_(default(_), Value, _, user(Value)).
 1856override_(user(Prev), Value, Options, _) :-
 1857        (   Value == Prev ->
 1858            domain_error(nonrepeating_labeling_options, Options)
 1859        ;   domain_error(consistent_labeling_options, Options)
 1860        ).
 1861
 1862selection(ff).
 1863selection(ffc).
 1864selection(min).
 1865selection(max).
 1866selection(leftmost).
 1867selection(random_variable(Seed)) :-
 1868        must_be(integer, Seed),
 1869        set_random(seed(Seed)).
 1870
 1871choice(step).
 1872choice(enum).
 1873choice(bisect).
 1874
 1875order(up).
 1876order(down).
 1877% TODO: random_variable and random_value currently both set the seed,
 1878% so exchanging the options can yield different results.
 1879order(random_value(Seed)) :-
 1880        must_be(integer, Seed),
 1881        set_random(seed(Seed)).
 1882
 1883consistency(upto_in(I), upto_in(1, I)).
 1884consistency(upto_in, upto_in).
 1885consistency(upto_ground, upto_ground).
 1886
 1887optimisation(min(_)).
 1888optimisation(max(_)).
 1889
 1890select_var(leftmost, [Var|Vars], Var, Vars).
 1891select_var(min, [V|Vs], Var, RVars) :-
 1892        find_min(Vs, V, Var),
 1893        delete_eq([V|Vs], Var, RVars).
 1894select_var(max, [V|Vs], Var, RVars) :-
 1895        find_max(Vs, V, Var),
 1896        delete_eq([V|Vs], Var, RVars).
 1897select_var(ff, [V|Vs], Var, RVars) :-
 1898        fd_size_(V, n(S)),
 1899        find_ff(Vs, V, S, Var),
 1900        delete_eq([V|Vs], Var, RVars).
 1901select_var(ffc, [V|Vs], Var, RVars) :-
 1902        find_ffc(Vs, V, Var),
 1903        delete_eq([V|Vs], Var, RVars).
 1904select_var(random_variable(_), Vars0, Var, Vars) :-
 1905        length(Vars0, L),
 1906        I is random(L),
 1907        nth0(I, Vars0, Var),
 1908        delete_eq(Vars0, Var, Vars).
 1909
 1910find_min([], Var, Var).
 1911find_min([V|Vs], CM, Min) :-
 1912        (   min_lt(V, CM) ->
 1913            find_min(Vs, V, Min)
 1914        ;   find_min(Vs, CM, Min)
 1915        ).
 1916
 1917find_max([], Var, Var).
 1918find_max([V|Vs], CM, Max) :-
 1919        (   max_gt(V, CM) ->
 1920            find_max(Vs, V, Max)
 1921        ;   find_max(Vs, CM, Max)
 1922        ).
 1923
 1924find_ff([], Var, _, Var).
 1925find_ff([V|Vs], CM, S0, FF) :-
 1926        (   nonvar(V) -> find_ff(Vs, CM, S0, FF)
 1927        ;   (   fd_size_(V, n(S1)), S1 < S0 ->
 1928                find_ff(Vs, V, S1, FF)
 1929            ;   find_ff(Vs, CM, S0, FF)
 1930            )
 1931        ).
 1932
 1933find_ffc([], Var, Var).
 1934find_ffc([V|Vs], Prev, FFC) :-
 1935        (   ffc_lt(V, Prev) ->
 1936            find_ffc(Vs, V, FFC)
 1937        ;   find_ffc(Vs, Prev, FFC)
 1938        ).
 1939
 1940
 1941ffc_lt(X, Y) :-
 1942        (   fd_get(X, XD, XPs) ->
 1943            domain_num_elements(XD, n(NXD))
 1944        ;   NXD = 1, XPs = []
 1945        ),
 1946        (   fd_get(Y, YD, YPs) ->
 1947            domain_num_elements(YD, n(NYD))
 1948        ;   NYD = 1, YPs = []
 1949        ),
 1950        (   NXD < NYD -> true
 1951        ;   NXD =:= NYD,
 1952            props_number(XPs, NXPs),
 1953            props_number(YPs, NYPs),
 1954            NXPs > NYPs
 1955        ).
 1956
 1957min_lt(X,Y) :- bounds(X,LX,_), bounds(Y,LY,_), LX < LY.
 1958
 1959max_gt(X,Y) :- bounds(X,_,UX), bounds(Y,_,UY), UX > UY.
 1960
 1961bounds(X, L, U) :-
 1962        (   fd_get(X, Dom, _) ->
 1963            domain_infimum(Dom, n(L)),
 1964            domain_supremum(Dom, n(U))
 1965        ;   L = X, U = L
 1966        ).
 1967
 1968delete_eq([], _, []).
 1969delete_eq([X|Xs], Y, List) :-
 1970        (   nonvar(X) -> delete_eq(Xs, Y, List)
 1971        ;   X == Y -> List = Xs
 1972        ;   List = [X|Tail],
 1973            delete_eq(Xs, Y, Tail)
 1974        ).
 1975
 1976/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1977   contracting/1 -- subject to change
 1978
 1979   This can remove additional domain elements from the boundaries.
 1980- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1981
 1982contracting(Vs) :-
 1983        must_be(list, Vs),
 1984        maplist(must_be_finite_fdvar, Vs),
 1985        contracting(Vs, false, Vs).
 1986
 1987contracting([], Repeat, Vars) :-
 1988        (   Repeat -> contracting(Vars, false, Vars)
 1989        ;   true
 1990        ).
 1991contracting([V|Vs], Repeat, Vars) :-
 1992        fd_inf(V, Min),
 1993        (   \+ \+ (V = Min) ->
 1994            fd_sup(V, Max),
 1995            (   \+ \+ (V = Max) ->
 1996                contracting(Vs, Repeat, Vars)
 1997            ;   V #\= Max,
 1998                contracting(Vs, true, Vars)
 1999            )
 2000        ;   V #\= Min,
 2001            contracting(Vs, true, Vars)
 2002        ).
 2003
 2004/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 2005   fds_sespsize(Vs, S).
 2006
 2007   S is an upper bound on the search space size with respect to finite
 2008   domain variables Vs.
 2009- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 2010
 2011fds_sespsize(Vs, S) :-
 2012        must_be(list, Vs),
 2013        maplist(fd_variable, Vs),
 2014        fds_sespsize(Vs, n(1), S1),
 2015        bound_portray(S1, S).
 2016
 2017fd_size_(V, S) :-
 2018        (   fd_get(V, D, _) ->
 2019            domain_num_elements(D, S)
 2020        ;   S = n(1)
 2021        ).
 2022
 2023fds_sespsize([], S, S).
 2024fds_sespsize([V|Vs], S0, S) :-
 2025        fd_size_(V, S1),
 2026        S2 cis S0*S1,
 2027        fds_sespsize(Vs, S2, S).
 2028
 2029/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 2030   Optimisation uses destructive assignment to save the computed
 2031   extremum over backtracking. Failure is used to get rid of copies of
 2032   attributed variables that are created in intermediate steps. At
 2033   least that's the intention - it currently doesn't work in SWI:
 2034
 2035   %?- X in 0..3, call_residue_vars(labeling([min(X)], [X]), Vs).
 2036   %@ X = 0,
 2037   %@ Vs = [_G6174, _G6177],
 2038   %@ _G6174 in 0..3
 2039
 2040- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 2041
 2042optimise(Vars, Options, Whats) :-
 2043        Whats = [What|WhatsRest],
 2044        Extremum = extremum(none),
 2045        (   catch(store_extremum(Vars, Options, What, Extremum),
 2046                  time_limit_exceeded,
 2047                  false)
 2048        ;   Extremum = extremum(n(Val)),
 2049            arg(1, What, Expr),
 2050            append(WhatsRest, Options, Options1),
 2051            (   Expr #= Val,
 2052                labeling(Options1, Vars)
 2053            ;   Expr #\= Val,
 2054                optimise(Vars, Options, Whats)
 2055            )
 2056        ).
 2057
 2058store_extremum(Vars, Options, What, Extremum) :-
 2059        catch((labeling(Options, Vars), throw(w(What))), w(What1), true),
 2060        functor(What, Direction, _),
 2061        maplist(arg(1), [What,What1], [Expr,Expr1]),
 2062        optimise(Direction, Options, Vars, Expr1, Expr, Extremum).
 2063
 2064optimise(Direction, Options, Vars, Expr0, Expr, Extremum) :-
 2065        must_be(ground, Expr0),
 2066        nb_setarg(1, Extremum, n(Expr0)),
 2067        catch((tighten(Direction, Expr, Expr0),
 2068               labeling(Options, Vars),
 2069               throw(v(Expr))), v(Expr1), true),
 2070        optimise(Direction, Options, Vars, Expr1, Expr, Extremum).
 2071
 2072tighten(min, E, V) :- E #< V.
 2073tighten(max, E, V) :- E #> V.
 2074
 2075%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 all_different(+Vars)
Like all_distinct/1, but with weaker propagation. Consider using all_distinct/1 instead, since all_distinct/1 is typically acceptably efficient and propagates much more strongly.
 2083all_different(Ls) :-
 2084        fd_must_be_list(Ls),
 2085        maplist(fd_variable, Ls),
 2086        Orig = original_goal(_, all_different(Ls)),
 2087        all_different(Ls, [], Orig),
 2088        do_queue.
 2089
 2090all_different([], _, _).
 2091all_different([X|Right], Left, Orig) :-
 2092        (   var(X) ->
 2093            make_propagator(pdifferent(Left,Right,X,Orig), Prop),
 2094            init_propagator(X, Prop),
 2095            trigger_prop(Prop)
 2096        ;   exclude_fire(Left, Right, X)
 2097        ),
 2098        all_different(Right, [X|Left], Orig).
 all_distinct(+Vars)
True iff Vars are pairwise distinct. For example, all_distinct/1 can detect that not all variables can assume distinct values given the following domains:
?- maplist(in, Vs,
           [1\/3..4, 1..2\/4, 1..2\/4, 1..3, 1..3, 1..6]),
   all_distinct(Vs).
false.
 2113all_distinct(Ls) :-
 2114        fd_must_be_list(Ls),
 2115        maplist(fd_variable, Ls),
 2116        make_propagator(pdistinct(Ls), Prop),
 2117        distinct_attach(Ls, Prop, []),
 2118        trigger_once(Prop).
 sum(+Vars, +Rel, ?Expr)
The sum of elements of the list Vars is in relation Rel to Expr. Rel is one of #=, #\=, #<, #>, #=< or #>=. For example:
?- [A,B,C] ins 0..sup, sum([A,B,C], #=, 100).
A in 0..100,
A+B+C#=100,
B in 0..100,
C in 0..100.
 2133sum(Vs, Op, Value) :-
 2134        must_be(list, Vs),
 2135        same_length(Vs, Ones),
 2136        maplist(=(1), Ones),
 2137        scalar_product(Ones, Vs, Op, Value).
 scalar_product(+Cs, +Vs, +Rel, ?Expr)
True iff the scalar product of Cs and Vs is in relation Rel to Expr. Cs is a list of integers, Vs is a list of variables and integers. Rel is #=, #\=, #<, #>, #=< or #>=.
 2145scalar_product(Cs, Vs, Op, Value) :-
 2146        must_be(list(integer), Cs),
 2147        must_be(list, Vs),
 2148        maplist(fd_variable, Vs),
 2149        (   Op = (#=), single_value(Value, Right), ground(Vs) ->
 2150            foldl(coeff_int_linsum, Cs, Vs, 0, Right)
 2151        ;   must_be(callable, Op),
 2152            (   memberchk(Op, [#=,#\=,#<,#>,#=<,#>=]) -> true
 2153            ;   domain_error(scalar_product_relation, Op)
 2154            ),
 2155            must_be(acyclic, Value),
 2156            foldl(coeff_var_plusterm, Cs, Vs, 0, Left),
 2157            (   left_right_linsum_const(Left, Value, Cs1, Vs1, Const) ->
 2158                scalar_product_(Op, Cs1, Vs1, Const)
 2159            ;   sum(Cs, Vs, 0, Op, Value)
 2160            )
 2161        ).
 2162
 2163single_value(V, V)    :- var(V), !, non_monotonic(V).
 2164single_value(V, V)    :- integer(V).
 2165single_value(?(V), V) :- fd_variable(V).
 2166
 2167coeff_var_plusterm(C, V, T0, T0+(C* ?(V))).
 2168
 2169coeff_int_linsum(C, I, S0, S) :- S is S0 + C*I.
 2170
 2171sum([], _, Sum, Op, Value) :- call(Op, Sum, Value).
 2172sum([C|Cs], [X|Xs], Acc, Op, Value) :-
 2173        ?(NAcc) #= Acc + C* ?(X),
 2174        sum(Cs, Xs, NAcc, Op, Value).
 2175
 2176multiples([], [], _).
 2177multiples([C|Cs], [V|Vs], Left) :-
 2178        (   (   Cs = [N|_] ; Left = [N|_] ) ->
 2179            (   N =\= 1, gcd(C,N) =:= 1 ->
 2180                gcd(Cs, N, GCD0),
 2181                gcd(Left, GCD0, GCD),
 2182                (   GCD > 1 -> ?(V) #= GCD * ?(_)
 2183                ;   true
 2184                )
 2185            ;   true
 2186            )
 2187        ;   true
 2188        ),
 2189        multiples(Cs, Vs, [C|Left]).
 2190
 2191abs(N, A) :- A is abs(N).
 2192
 2193divide(D, N, R) :- R is N // D.
 2194
 2195scalar_product_(#=, Cs0, Vs, S0) :-
 2196        (   Cs0 = [C|Rest] ->
 2197            gcd(Rest, C, GCD),
 2198            S0 mod GCD =:= 0,
 2199            maplist(divide(GCD), [S0|Cs0], [S|Cs])
 2200        ;   S0 =:= 0, S = S0, Cs = Cs0
 2201        ),
 2202        (   S0 =:= 0 ->
 2203            maplist(abs, Cs, As),
 2204            multiples(As, Vs, [])
 2205        ;   true
 2206        ),
 2207        propagator_init_trigger(Vs, scalar_product_eq(Cs, Vs, S)).
 2208scalar_product_(#\=, Cs, Vs, C) :-
 2209        propagator_init_trigger(Vs, scalar_product_neq(Cs, Vs, C)).
 2210scalar_product_(#=<, Cs, Vs, C) :-
 2211        propagator_init_trigger(Vs, scalar_product_leq(Cs, Vs, C)).
 2212scalar_product_(#<, Cs, Vs, C) :-
 2213        C1 is C - 1,
 2214        scalar_product_(#=<, Cs, Vs, C1).
 2215scalar_product_(#>, Cs, Vs, C) :-
 2216        C1 is C + 1,
 2217        scalar_product_(#>=, Cs, Vs, C1).
 2218scalar_product_(#>=, Cs, Vs, C) :-
 2219        maplist(negative, Cs, Cs1),
 2220        C1 is -C,
 2221        scalar_product_(#=<, Cs1, Vs, C1).
 2222
 2223negative(X0, X) :- X is -X0.
 2224
 2225coeffs_variables_const([], [], [], [], I, I).
 2226coeffs_variables_const([C|Cs], [V|Vs], Cs1, Vs1, I0, I) :-
 2227        (   var(V) ->
 2228            Cs1 = [C|CRest], Vs1 = [V|VRest], I1 = I0
 2229        ;   I1 is I0 + C*V,
 2230            Cs1 = CRest, Vs1 = VRest
 2231        ),
 2232        coeffs_variables_const(Cs, Vs, CRest, VRest, I1, I).
 2233
 2234sum_finite_domains([], [], [], [], Inf, Sup, Inf, Sup).
 2235sum_finite_domains([C|Cs], [V|Vs], Infs, Sups, Inf0, Sup0, Inf, Sup) :-
 2236        fd_get(V, _, Inf1, Sup1, _),
 2237        (   Inf1 = n(NInf) ->
 2238            (   C < 0 ->
 2239                Sup2 is Sup0 + C*NInf
 2240            ;   Inf2 is Inf0 + C*NInf
 2241            ),
 2242            Sups = Sups1,
 2243            Infs = Infs1
 2244        ;   (   C < 0 ->
 2245                Sup2 = Sup0,
 2246                Sups = [C*V|Sups1],
 2247                Infs = Infs1
 2248            ;   Inf2 = Inf0,
 2249                Infs = [C*V|Infs1],
 2250                Sups = Sups1
 2251            )
 2252        ),
 2253        (   Sup1 = n(NSup) ->
 2254            (   C < 0 ->
 2255                Inf2 is Inf0 + C*NSup
 2256            ;   Sup2 is Sup0 + C*NSup
 2257            ),
 2258            Sups1 = Sups2,
 2259            Infs1 = Infs2
 2260        ;   (   C < 0 ->
 2261                Inf2 = Inf0,
 2262                Infs1 = [C*V|Infs2],
 2263                Sups1 = Sups2
 2264            ;   Sup2 = Sup0,
 2265                Sups1 = [C*V|Sups2],
 2266                Infs1 = Infs2
 2267            )
 2268        ),
 2269        sum_finite_domains(Cs, Vs, Infs2, Sups2, Inf2, Sup2, Inf, Sup).
 2270
 2271remove_dist_upper_lower([], _, _, _).
 2272remove_dist_upper_lower([C|Cs], [V|Vs], D1, D2) :-
 2273        (   fd_get(V, VD, VPs) ->
 2274            (   C < 0 ->
 2275                domain_supremum(VD, n(Sup)),
 2276                L is Sup + D1//C,
 2277                domain_remove_smaller_than(VD, L, VD1),
 2278                domain_infimum(VD1, n(Inf)),
 2279                G is Inf - D2//C,
 2280                domain_remove_greater_than(VD1, G, VD2)
 2281            ;   domain_infimum(VD, n(Inf)),
 2282                G is Inf + D1//C,
 2283                domain_remove_greater_than(VD, G, VD1),
 2284                domain_supremum(VD1, n(Sup)),
 2285                L is Sup - D2//C,
 2286                domain_remove_smaller_than(VD1, L, VD2)
 2287            ),
 2288            fd_put(V, VD2, VPs)
 2289        ;   true
 2290        ),
 2291        remove_dist_upper_lower(Cs, Vs, D1, D2).
 2292
 2293
 2294remove_dist_upper_leq([], _, _).
 2295remove_dist_upper_leq([C|Cs], [V|Vs], D1) :-
 2296        (   fd_get(V, VD, VPs) ->
 2297            (   C < 0 ->
 2298                domain_supremum(VD, n(Sup)),
 2299                L is Sup + D1//C,
 2300                domain_remove_smaller_than(VD, L, VD1)
 2301            ;   domain_infimum(VD, n(Inf)),
 2302                G is Inf + D1//C,
 2303                domain_remove_greater_than(VD, G, VD1)
 2304            ),
 2305            fd_put(V, VD1, VPs)
 2306        ;   true
 2307        ),
 2308        remove_dist_upper_leq(Cs, Vs, D1).
 2309
 2310
 2311remove_dist_upper([], _).
 2312remove_dist_upper([C*V|CVs], D) :-
 2313        (   fd_get(V, VD, VPs) ->
 2314            (   C < 0 ->
 2315                (   domain_supremum(VD, n(Sup)) ->
 2316                    L is Sup + D//C,
 2317                    domain_remove_smaller_than(VD, L, VD1)
 2318                ;   VD1 = VD
 2319                )
 2320            ;   (   domain_infimum(VD, n(Inf)) ->
 2321                    G is Inf + D//C,
 2322                    domain_remove_greater_than(VD, G, VD1)
 2323                ;   VD1 = VD
 2324                )
 2325            ),
 2326            fd_put(V, VD1, VPs)
 2327        ;   true
 2328        ),
 2329        remove_dist_upper(CVs, D).
 2330
 2331remove_dist_lower([], _).
 2332remove_dist_lower([C*V|CVs], D) :-
 2333        (   fd_get(V, VD, VPs) ->
 2334            (   C < 0 ->
 2335                (   domain_infimum(VD, n(Inf)) ->
 2336                    G is Inf - D//C,
 2337                    domain_remove_greater_than(VD, G, VD1)
 2338                ;   VD1 = VD
 2339                )
 2340            ;   (   domain_supremum(VD, n(Sup)) ->
 2341                    L is Sup - D//C,
 2342                    domain_remove_smaller_than(VD, L, VD1)
 2343                ;   VD1 = VD
 2344                )
 2345            ),
 2346            fd_put(V, VD1, VPs)
 2347        ;   true
 2348        ),
 2349        remove_dist_lower(CVs, D).
 2350
 2351remove_upper([], _).
 2352remove_upper([C*X|CXs], Max) :-
 2353        (   fd_get(X, XD, XPs) ->
 2354            D is Max//C,
 2355            (   C < 0 ->
 2356                domain_remove_smaller_than(XD, D, XD1)
 2357            ;   domain_remove_greater_than(XD, D, XD1)
 2358            ),
 2359            fd_put(X, XD1, XPs)
 2360        ;   true
 2361        ),
 2362        remove_upper(CXs, Max).
 2363
 2364remove_lower([], _).
 2365remove_lower([C*X|CXs], Min) :-
 2366        (   fd_get(X, XD, XPs) ->
 2367            D is -Min//C,
 2368            (   C < 0 ->
 2369                domain_remove_greater_than(XD, D, XD1)
 2370            ;   domain_remove_smaller_than(XD, D, XD1)
 2371            ),
 2372            fd_put(X, XD1, XPs)
 2373        ;   true
 2374        ),
 2375        remove_lower(CXs, Min).
 2376
 2377%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 2378
 2379/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 2380   Constraint propagation proceeds as follows: Each CLP(FD) variable
 2381   has an attribute that stores its associated domain and constraints.
 2382   Constraints are triggered when the event they are registered for
 2383   occurs (for example: variable is instantiated, bounds change etc.).
 2384   do_queue/0 works off all triggered constraints, possibly triggering
 2385   new ones, until fixpoint.
 2386- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 2387
 2388% FIFO queue
 2389
 2390make_queue :- nb_setval('$clpfd_queue', fast_slow([], [])).
 2391
 2392push_queue(E, Which) :-
 2393        nb_getval('$clpfd_queue', Qs),
 2394        arg(Which, Qs, Q),
 2395        (   Q == [] ->
 2396            setarg(Which, Qs, [E|T]-T)
 2397        ;   Q = H-[E|T],
 2398            setarg(Which, Qs, H-T)
 2399        ).
 2400
 2401pop_queue(E) :-
 2402        nb_getval('$clpfd_queue', Qs),
 2403        (   pop_queue(E, Qs, 1) ->  true
 2404        ;   pop_queue(E, Qs, 2)
 2405        ).
 2406
 2407pop_queue(E, Qs, Which) :-
 2408        arg(Which, Qs, [E|NH]-T),
 2409        (   var(NH) ->
 2410            setarg(Which, Qs, [])
 2411        ;   setarg(Which, Qs, NH-T)
 2412        ).
 2413
 2414fetch_propagator(Prop) :-
 2415        pop_queue(P),
 2416        (   propagator_state(P, S), S == dead -> fetch_propagator(Prop)
 2417        ;   Prop = P
 2418        ).
 2419
 2420/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 2421   Parsing a CLP(FD) expression has two important side-effects: First,
 2422   it constrains the variables occurring in the expression to
 2423   integers. Second, it constrains some of them even more: For
 2424   example, in X/Y and X mod Y, Y is constrained to be #\= 0.
 2425- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 2426
 2427constrain_to_integer(Var) :-
 2428        (   integer(Var) -> true
 2429        ;   fd_get(Var, D, Ps),
 2430            fd_put(Var, D, Ps)
 2431        ).
 2432
 2433power_var_num(P, X, N) :-
 2434        (   var(P) -> X = P, N = 1
 2435        ;   P = Left*Right,
 2436            power_var_num(Left, XL, L),
 2437            power_var_num(Right, XR, R),
 2438            XL == XR,
 2439            X = XL,
 2440            N is L + R
 2441        ).
 2442
 2443/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 2444   Given expression E, we obtain the finite domain variable R by
 2445   interpreting a simple committed-choice language that is a list of
 2446   conditions and bodies. In conditions, g(Goal) means literally Goal,
 2447   and m(Match) means that E can be decomposed as stated. The
 2448   variables are to be understood as the result of parsing the
 2449   subexpressions recursively. In the body, g(Goal) means again Goal,
 2450   and p(Propagator) means to attach and trigger once a propagator.
 2451- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 2452
 2453:- op(800, xfx, =>). 2454
 2455parse_clpfd(E, R,
 2456            [g(cyclic_term(E)) => [g(domain_error(clpfd_expression, E))],
 2457             g(var(E))         => [g(non_monotonic(E)),
 2458                                   g(constrain_to_integer(E)), g(E = R)],
 2459             g(integer(E))     => [g(R = E)],
 2460             ?(E)              => [g(must_be_fd_integer(E)), g(R = E)],
 2461             #(E)              => [g(must_be_fd_integer(E)), g(R = E)],
 2462             m(A+B)            => [p(pplus(A, B, R))],
 2463             % power_var_num/3 must occur before */2 to be useful
 2464             g(power_var_num(E, V, N)) => [p(pexp(V, N, R))],
 2465             m(A*B)            => [p(ptimes(A, B, R))],
 2466             m(A-B)            => [p(pplus(R,B,A))],
 2467             m(-A)             => [p(ptimes(-1,A,R))],
 2468             m(max(A,B))       => [g(A #=< ?(R)), g(B #=< R), p(pmax(A, B, R))],
 2469             m(min(A,B))       => [g(A #>= ?(R)), g(B #>= R), p(pmin(A, B, R))],
 2470             m(A mod B)        => [g(B #\= 0), p(pmod(A, B, R))],
 2471             m(A rem B)        => [g(B #\= 0), p(prem(A, B, R))],
 2472             m(abs(A))         => [g(?(R) #>= 0), p(pabs(A, R))],
 2473%             m(A/B)            => [g(B #\= 0), p(ptzdiv(A, B, R))],
 2474             m(A//B)           => [g(B #\= 0), p(ptzdiv(A, B, R))],
 2475             m(A div B)        => [g(?(R) #= (A - (A mod B)) // B)],
 2476             m(A rdiv B)       => [g(B #\= 0), p(prdiv(A, B, R))],
 2477             m(A^B)            => [p(pexp(A, B, R))],
 2478             % bitwise operations
 2479             m(\A)             => [p(pfunction(\, A, R))],
 2480             m(msb(A))         => [p(pfunction(msb, A, R))],
 2481             m(lsb(A))         => [p(pfunction(lsb, A, R))],
 2482             m(popcount(A))    => [p(pfunction(popcount, A, R))],
 2483             m(A<<B)           => [p(pfunction(<<, A, B, R))],
 2484             m(A>>B)           => [p(pfunction(>>, A, B, R))],
 2485             m(A/\B)           => [p(pfunction(/\, A, B, R))],
 2486             m(A\/B)           => [p(pfunction(\/, A, B, R))],
 2487             m(A xor B)        => [p(pfunction(xor, A, B, R))],
 2488             g(true)           => [g(domain_error(clpfd_expression, E))]
 2489            ]).
 2490
 2491non_monotonic(X) :-
 2492        (   \+ fd_var(X), current_prolog_flag(clpfd_monotonic, true) ->
 2493            instantiation_error(X)
 2494        ;   true
 2495        ).
 2496
 2497% Here, we compile the committed choice language to a single
 2498% predicate, parse_clpfd/2.
 2499
 2500make_parse_clpfd(Clauses) :-
 2501        parse_clpfd_clauses(Clauses0),
 2502        maplist(goals_goal, Clauses0, Clauses).
 2503
 2504goals_goal((Head :- Goals), (Head :- Body)) :-
 2505        list_goal(Goals, Body).
 2506
 2507parse_clpfd_clauses(Clauses) :-
 2508        parse_clpfd(E, R, Matchers),
 2509        maplist(parse_matcher(E, R), Matchers, Clauses).
 2510
 2511parse_matcher(E, R, Matcher, Clause) :-
 2512        Matcher = (Condition0 => Goals0),
 2513        phrase((parse_condition(Condition0, E, Head),
 2514                parse_goals(Goals0)), Goals),
 2515        Clause = (parse_clpfd(Head, R) :- Goals).
 2516
 2517parse_condition(g(Goal), E, E)       --> [Goal, !].
 2518parse_condition(?(E), _, ?(E))       --> [!].
 2519parse_condition(#(E), _, #(E))       --> [!].
 2520parse_condition(m(Match), _, Match0) -->
 2521        [!],
 2522        { copy_term(Match, Match0),
 2523          term_variables(Match0, Vs0),
 2524          term_variables(Match, Vs)
 2525        },
 2526        parse_match_variables(Vs0, Vs).
 2527
 2528parse_match_variables([], []) --> [].
 2529parse_match_variables([V0|Vs0], [V|Vs]) -->
 2530        [parse_clpfd(V0, V)],
 2531        parse_match_variables(Vs0, Vs).
 2532
 2533parse_goals([]) --> [].
 2534parse_goals([G|Gs]) --> parse_goal(G), parse_goals(Gs).
 2535
 2536parse_goal(g(Goal)) --> [Goal].
 2537parse_goal(p(Prop)) -->
 2538        [make_propagator(Prop, P)],
 2539        { term_variables(Prop, Vs) },
 2540        parse_init(Vs, P),
 2541        [trigger_once(P)].
 2542
 2543parse_init([], _)     --> [].
 2544parse_init([V|Vs], P) --> [init_propagator(V, P)], parse_init(Vs, P).
 2545
 2546%?- set_prolog_flag(answer_write_options, [portray(true)]),
 2547%   clpfd:parse_clpfd_clauses(Clauses), maplist(portray_clause, Clauses).
 2548
 2549
 2550%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 2551%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 2552
 2553trigger_once(Prop) :- trigger_prop(Prop), do_queue.
 2554
 2555neq(A, B) :- propagator_init_trigger(pneq(A, B)).
 2556
 2557propagator_init_trigger(P) -->
 2558        { term_variables(P, Vs) },
 2559        propagator_init_trigger(Vs, P).
 2560
 2561propagator_init_trigger(Vs, P) -->
 2562        [p(Prop)],
 2563        { make_propagator(P, Prop),
 2564          maplist(prop_init(Prop), Vs),
 2565          trigger_once(Prop) }.
 2566
 2567propagator_init_trigger(P) :-
 2568        phrase(propagator_init_trigger(P), _).
 2569
 2570propagator_init_trigger(Vs, P) :-
 2571        phrase(propagator_init_trigger(Vs, P), _).
 2572
 2573prop_init(Prop, V) :- init_propagator(V, Prop).
 2574
 2575geq(A, B) :-
 2576        (   fd_get(A, AD, APs) ->
 2577            domain_infimum(AD, AI),
 2578            (   fd_get(B, BD, _) ->
 2579                domain_supremum(BD, BS),
 2580                (   AI cis_geq BS -> true
 2581                ;   propagator_init_trigger(pgeq(A,B))
 2582                )
 2583            ;   (   AI cis_geq n(B) -> true
 2584                ;   domain_remove_smaller_than(AD, B, AD1),
 2585                    fd_put(A, AD1, APs),
 2586                    do_queue
 2587                )
 2588            )
 2589        ;   fd_get(B, BD, BPs) ->
 2590            domain_remove_greater_than(BD, A, BD1),
 2591            fd_put(B, BD1, BPs),
 2592            do_queue
 2593        ;   A >= B
 2594        ).
 2595
 2596/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 2597   Naive parsing of inequalities and disequalities can result in a lot
 2598   of unnecessary work if expressions of non-trivial depth are
 2599   involved: Auxiliary variables are introduced for sub-expressions,
 2600   and propagation proceeds on them as if they were involved in a
 2601   tighter constraint (like equality), whereas eventually only very
 2602   little of the propagated information is actually used. For example,
 2603   only extremal values are of interest in inequalities. Introducing
 2604   auxiliary variables should be avoided when possible, and
 2605   specialised propagators should be used for common constraints.
 2606
 2607   We again use a simple committed-choice language for matching
 2608   special cases of constraints. m_c(M,C) means that M matches and C
 2609   holds. d(X, Y) means decomposition, i.e., it is short for
 2610   g(parse_clpfd(X, Y)). r(X, Y) means to rematch with X and Y.
 2611
 2612   Two things are important: First, although the actual constraint
 2613   functors (#\=2, #=/2 etc.) are used in the description, they must
 2614   expand to the respective auxiliary predicates (match_expand/2)
 2615   because the actual constraints are subject to goal expansion.
 2616   Second, when specialised constraints (like scalar product) post
 2617   simpler constraints on their own, these simpler versions must be
 2618   handled separately and must occur before.
 2619- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 2620
 2621match_expand(#>=, clpfd_geq_).
 2622match_expand(#=, clpfd_equal_).
 2623match_expand(#\=, clpfd_neq).
 2624
 2625symmetric(#=).
 2626symmetric(#\=).
 2627
 2628matches([
 2629         m_c(any(X) #>= any(Y), left_right_linsum_const(X, Y, Cs, Vs, Const)) =>
 2630            [g((   Cs = [1], Vs = [A] -> geq(A, Const)
 2631               ;   Cs = [-1], Vs = [A] -> Const1 is -Const, geq(Const1, A)
 2632               ;   Cs = [1,1], Vs = [A,B] -> ?(A) + ?(B) #= ?(S), geq(S, Const)
 2633               ;   Cs = [1,-1], Vs = [A,B] ->
 2634                   (   Const =:= 0 -> geq(A, B)
 2635                   ;   C1 is -Const,
 2636                       propagator_init_trigger(x_leq_y_plus_c(B, A, C1))
 2637                   )
 2638               ;   Cs = [-1,1], Vs = [A,B] ->
 2639                   (   Const =:= 0 -> geq(B, A)
 2640                   ;   C1 is -Const,
 2641                       propagator_init_trigger(x_leq_y_plus_c(A, B, C1))
 2642                   )
 2643               ;   Cs = [-1,-1], Vs = [A,B] ->
 2644                   ?(A) + ?(B) #= ?(S), Const1 is -Const, geq(Const1, S)
 2645               ;   scalar_product_(#>=, Cs, Vs, Const)
 2646               ))],
 2647         m(any(X) - any(Y) #>= integer(C))     => [d(X, X1), d(Y, Y1), g(C1 is -C), p(x_leq_y_plus_c(Y1, X1, C1))],
 2648         m(integer(X) #>= any(Z) + integer(A)) => [g(C is X - A), r(C, Z)],
 2649         m(abs(any(X)-any(Y)) #>= integer(I))  => [d(X, X1), d(Y, Y1), p(absdiff_geq(X1, Y1, I))],
 2650         m(abs(any(X)) #>= integer(I))         => [d(X, RX), g((I>0 -> I1 is -I, RX in inf..I1 \/ I..sup; true))],
 2651         m(integer(I) #>= abs(any(X)))         => [d(X, RX), g(I>=0), g(I1 is -I), g(RX in I1..I)],
 2652         m(any(X) #>= any(Y))                  => [d(X, RX), d(Y, RY), g(geq(RX, RY))],
 2653
 2654         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 2655         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 2656
 2657         m(var(X) #= var(Y))        => [g(constrain_to_integer(X)), g(X=Y)],
 2658         m(var(X) #= var(Y)+var(Z)) => [p(pplus(Y,Z,X))],
 2659         m(var(X) #= var(Y)-var(Z)) => [p(pplus(X,Z,Y))],
 2660         m(var(X) #= var(Y)*var(Z)) => [p(ptimes(Y,Z,X))],
 2661         m(var(X) #= -var(Z))       => [p(ptimes(-1, Z, X))],
 2662         m_c(any(X) #= any(Y), left_right_linsum_const(X, Y, Cs, Vs, S)) =>
 2663            [g(scalar_product_(#=, Cs, Vs, S))],
 2664         m(var(X) #= any(Y))       => [d(Y,X)],
 2665         m(any(X) #= any(Y))       => [d(X, RX), d(Y, RX)],
 2666
 2667         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 2668         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 2669
 2670         m(var(X) #\= integer(Y))             => [g(neq_num(X, Y))],
 2671         m(var(X) #\= var(Y))                 => [g(neq(X,Y))],
 2672         m(var(X) #\= var(Y) + var(Z))        => [p(x_neq_y_plus_z(X, Y, Z))],
 2673         m(var(X) #\= var(Y) - var(Z))        => [p(x_neq_y_plus_z(Y, X, Z))],
 2674         m(var(X) #\= var(Y)*var(Z))          => [p(ptimes(Y,Z,P)), g(neq(X,P))],
 2675         m(integer(X) #\= abs(any(Y)-any(Z))) => [d(Y, Y1), d(Z, Z1), p(absdiff_neq(Y1, Z1, X))],
 2676         m_c(any(X) #\= any(Y), left_right_linsum_const(X, Y, Cs, Vs, S)) =>
 2677            [g(scalar_product_(#\=, Cs, Vs, S))],
 2678         m(any(X) #\= any(Y) + any(Z))        => [d(X, X1), d(Y, Y1), d(Z, Z1), p(x_neq_y_plus_z(X1, Y1, Z1))],
 2679         m(any(X) #\= any(Y) - any(Z))        => [d(X, X1), d(Y, Y1), d(Z, Z1), p(x_neq_y_plus_z(Y1, X1, Z1))],
 2680         m(any(X) #\= any(Y)) => [d(X, RX), d(Y, RY), g(neq(RX, RY))]
 2681        ]).
 2682
 2683/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 2684   We again compile the committed-choice matching language to the
 2685   intended auxiliary predicates. We now must take care not to
 2686   unintentionally unify a variable with a complex term.
 2687- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 2688
 2689make_matches(Clauses) :-
 2690        matches(Ms),
 2691        findall(F, (member(M=>_, Ms), arg(1, M, M1), functor(M1, F, _)), Fs0),
 2692        sort(Fs0, Fs),
 2693        maplist(prevent_cyclic_argument, Fs, PrevCyclicClauses),
 2694        phrase(matchers(Ms), Clauses0),
 2695        maplist(goals_goal, Clauses0, MatcherClauses),
 2696        append(PrevCyclicClauses, MatcherClauses, Clauses1),
 2697        sort_by_predicate(Clauses1, Clauses).
 2698
 2699sort_by_predicate(Clauses, ByPred) :-
 2700        map_list_to_pairs(predname, Clauses, Keyed),
 2701        keysort(Keyed, KeyedByPred),
 2702        pairs_values(KeyedByPred, ByPred).
 2703
 2704predname((H:-_), Key)   :- !, predname(H, Key).
 2705predname(M:H, M:Key)    :- !, predname(H, Key).
 2706predname(H, Name/Arity) :- !, functor(H, Name, Arity).
 2707
 2708prevent_cyclic_argument(F0, Clause) :-
 2709        match_expand(F0, F),
 2710        Head =.. [F,X,Y],
 2711        Clause = (Head :- (   cyclic_term(X) ->
 2712                              domain_error(clpfd_expression, X)
 2713                          ;   cyclic_term(Y) ->
 2714                              domain_error(clpfd_expression, Y)
 2715                          ;   false
 2716                          )).
 2717
 2718matchers([]) --> [].
 2719matchers([Condition => Goals|Ms]) -->
 2720        matcher(Condition, Goals),
 2721        matchers(Ms).
 2722
 2723matcher(m(M), Gs) --> matcher(m_c(M,true), Gs).
 2724matcher(m_c(Matcher,Cond), Gs) -->
 2725        [(Head :- Goals0)],
 2726        { Matcher =.. [F,A,B],
 2727          match_expand(F, Expand),
 2728          Head =.. [Expand,X,Y],
 2729          phrase((match(A, X), match(B, Y)), Goals0, [Cond,!|Goals1]),
 2730          phrase(match_goals(Gs, Expand), Goals1) },
 2731        (   { symmetric(F), \+ (subsumes_term(A, B), subsumes_term(B, A)) } ->
 2732            { Head1 =.. [Expand,Y,X] },
 2733            [(Head1 :- Goals0)]
 2734        ;   []
 2735        ).
 2736
 2737match(any(A), T)     --> [A = T].
 2738match(var(V), T)     --> [( nonvar(T), ( T = ?(Var) ; T = #(Var) ) ->
 2739                            must_be_fd_integer(Var), V = Var
 2740                          ; v_or_i(T), V = T
 2741                          )].
 2742match(integer(I), T) --> [integer(T), I = T].
 2743match(-X, T)         --> [nonvar(T), T = -A], match(X, A).
 2744match(abs(X), T)     --> [nonvar(T), T = abs(A)], match(X, A).
 2745match(Binary, T)     -->
 2746        { Binary =.. [Op,X,Y], Term =.. [Op,A,B] },
 2747        [nonvar(T), T = Term],
 2748        match(X, A), match(Y, B).
 2749
 2750match_goals([], _)     --> [].
 2751match_goals([G|Gs], F) --> match_goal(G, F), match_goals(Gs, F).
 2752
 2753match_goal(r(X,Y), F)  --> { G =.. [F,X,Y] }, [G].
 2754match_goal(d(X,Y), _)  --> [parse_clpfd(X, Y)].
 2755match_goal(g(Goal), _) --> [Goal].
 2756match_goal(p(Prop), _) -->
 2757        [make_propagator(Prop, P)],
 2758        { term_variables(Prop, Vs) },
 2759        parse_init(Vs, P),
 2760        [trigger_once(P)].
 2761
 2762
 2763%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 #>=(?X, ?Y)
Same as Y #=< X. When reasoning over integers, replace (>=)/2 by #>=/2 to obtain more general relations. See declarative integer arithmetic.
 2773X #>= Y :- clpfd_geq(X, Y).
 2774
 2775clpfd_geq(X, Y) :- clpfd_geq_(X, Y), reinforce(X), reinforce(Y).
 #=<(?X, ?Y)
The arithmetic expression X is less than or equal to Y. When reasoning over integers, replace (=<)/2 by #=</2 to obtain more general relations. See declarative integer arithmetic.
 2784X #=< Y :- Y #>= X.
 #=(?X, ?Y)
The arithmetic expression X equals Y. This is the most important arithmetic constraint, subsuming and replacing both (is)/2 and (=:=)/2 over integers. See declarative integer arithmetic.
 2793X #= Y :- clpfd_equal(X, Y).
 2794
 2795clpfd_equal(X, Y) :- clpfd_equal_(X, Y), reinforce(X).
 2796
 2797/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 2798   Conditions under which an equality can be compiled to built-in
 2799   arithmetic. Their order is significant. (/)/2 becomes (//)/2.
 2800- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 2801
 2802expr_conds(E, E)                 --> [integer(E)],
 2803        { var(E), !, \+ current_prolog_flag(clpfd_monotonic, true) }.
 2804expr_conds(E, E)                 --> { integer(E) }.
 2805expr_conds(?(E), E)              --> [integer(E)].
 2806expr_conds(#(E), E)              --> [integer(E)].
 2807expr_conds(-E0, -E)              --> expr_conds(E0, E).
 2808expr_conds(abs(E0), abs(E))      --> expr_conds(E0, E).
 2809expr_conds(A0+B0, A+B)           --> expr_conds(A0, A), expr_conds(B0, B).
 2810expr_conds(A0*B0, A*B)           --> expr_conds(A0, A), expr_conds(B0, B).
 2811expr_conds(A0-B0, A-B)           --> expr_conds(A0, A), expr_conds(B0, B).
 2812expr_conds(A0//B0, A//B)         -->
 2813        expr_conds(A0, A), expr_conds(B0, B),
 2814        [B =\= 0].
 2815%expr_conds(A0/B0, AB)            --> expr_conds(A0//B0, AB).
 2816expr_conds(min(A0,B0), min(A,B)) --> expr_conds(A0, A), expr_conds(B0, B).
 2817expr_conds(max(A0,B0), max(A,B)) --> expr_conds(A0, A), expr_conds(B0, B).
 2818expr_conds(A0 mod B0, A mod B)   -->
 2819        expr_conds(A0, A), expr_conds(B0, B),
 2820        [B =\= 0].
 2821expr_conds(A0^B0, A^B)           -->
 2822        expr_conds(A0, A), expr_conds(B0, B),
 2823        [(B >= 0 ; A =:= -1)].
 2824% Bitwise operations, added to make CLP(FD) usable in more cases
 2825expr_conds(\ A0, \ A) --> expr_conds(A0, A).
 2826expr_conds(A0<<B0, A<<B) --> expr_conds(A0, A), expr_conds(B0, B).
 2827expr_conds(A0>>B0, A>>B) --> expr_conds(A0, A), expr_conds(B0, B).
 2828expr_conds(A0/\B0, A/\B) --> expr_conds(A0, A), expr_conds(B0, B).
 2829expr_conds(A0\/B0, A\/B) --> expr_conds(A0, A), expr_conds(B0, B).
 2830expr_conds(A0 xor B0, A xor B) --> expr_conds(A0, A), expr_conds(B0, B).
 2831expr_conds(lsb(A0), lsb(A)) --> expr_conds(A0, A).
 2832expr_conds(msb(A0), msb(A)) --> expr_conds(A0, A).
 2833expr_conds(popcount(A0), popcount(A)) --> expr_conds(A0, A).
 2834
 2835:- multifile
 2836        system:goal_expansion/2. 2837:- dynamic
 2838        system:goal_expansion/2. 2839
 2840system:goal_expansion(Goal, Expansion) :-
 2841        \+ current_prolog_flag(clpfd_goal_expansion, false),
 2842        clpfd_expandable(Goal),
 2843        prolog_load_context(module, M),
 2844	(   M == clpfd
 2845	->  true
 2846	;   predicate_property(M:Goal, imported_from(clpfd))
 2847	),
 2848        clpfd_expansion(Goal, Expansion).
 2849
 2850clpfd_expandable(_ in _).
 2851clpfd_expandable(_ #= _).
 2852clpfd_expandable(_ #>= _).
 2853clpfd_expandable(_ #=< _).
 2854clpfd_expandable(_ #> _).
 2855clpfd_expandable(_ #< _).
 2856
 2857clpfd_expansion(Var in Dom, In) :-
 2858        (   ground(Dom), Dom = L..U, integer(L), integer(U) ->
 2859            expansion_simpler(
 2860                (   integer(Var) ->
 2861                    between(L, U, Var)
 2862                ;   clpfd:clpfd_in(Var, Dom)
 2863                ), In)
 2864        ;   In = clpfd:clpfd_in(Var, Dom)
 2865        ).
 2866clpfd_expansion(X0 #= Y0, Equal) :-
 2867        phrase(expr_conds(X0, X), CsX),
 2868        phrase(expr_conds(Y0, Y), CsY),
 2869        list_goal(CsX, CondX),
 2870        list_goal(CsY, CondY),
 2871        expansion_simpler(
 2872                (   CondX ->
 2873                    (   var(Y) -> Y is X
 2874                    ;   CondY -> X =:= Y
 2875                    ;   T is X, clpfd:clpfd_equal(T, Y0)
 2876                    )
 2877                ;   CondY ->
 2878                    (   var(X) -> X is Y
 2879                    ;   T is Y, clpfd:clpfd_equal(X0, T)
 2880                    )
 2881                ;   clpfd:clpfd_equal(X0, Y0)
 2882                ), Equal).
 2883clpfd_expansion(X0 #>= Y0, Geq) :-
 2884        phrase(expr_conds(X0, X), CsX),
 2885        phrase(expr_conds(Y0, Y), CsY),
 2886        list_goal(CsX, CondX),
 2887        list_goal(CsY, CondY),
 2888        expansion_simpler(
 2889              (   CondX ->
 2890                  (   CondY -> X >= Y
 2891                  ;   T is X, clpfd:clpfd_geq(T, Y0)
 2892                  )
 2893              ;   CondY -> T is Y, clpfd:clpfd_geq(X0, T)
 2894              ;   clpfd:clpfd_geq(X0, Y0)
 2895              ), Geq).
 2896clpfd_expansion(X #=< Y,  Leq) :- clpfd_expansion(Y #>= X, Leq).
 2897clpfd_expansion(X #> Y, Gt)    :- clpfd_expansion(X #>= Y+1, Gt).
 2898clpfd_expansion(X #< Y, Lt)    :- clpfd_expansion(Y #> X, Lt).
 2899
 2900expansion_simpler((True->Then0;_), Then) :-
 2901        is_true(True), !,
 2902        expansion_simpler(Then0, Then).
 2903expansion_simpler((False->_;Else0), Else) :-
 2904        is_false(False), !,
 2905        expansion_simpler(Else0, Else).
 2906expansion_simpler((If->Then0;Else0), (If->Then;Else)) :- !,
 2907        expansion_simpler(Then0, Then),
 2908        expansion_simpler(Else0, Else).
 2909expansion_simpler((A0,B0), (A,B)) :-
 2910        expansion_simpler(A0, A),
 2911        expansion_simpler(B0, B).
 2912expansion_simpler(Var is Expr0, Goal) :-
 2913        ground(Expr0), !,
 2914        phrase(expr_conds(Expr0, Expr), Gs),
 2915        (   maplist(call, Gs) -> Value is Expr, Goal = (Var = Value)
 2916        ;   Goal = false
 2917        ).
 2918expansion_simpler(Var =:= Expr0, Goal) :-
 2919        ground(Expr0), !,
 2920        phrase(expr_conds(Expr0, Expr), Gs),
 2921        (   maplist(call, Gs) -> Value is Expr, Goal = (Var =:= Value)
 2922        ;   Goal = false
 2923        ).
 2924expansion_simpler(Var is Expr, Var = Expr) :- var(Expr), !.
 2925expansion_simpler(Var is Expr, Goal) :- !,
 2926        (   var(Var), nonvar(Expr),
 2927            Expr = E mod M, nonvar(E), E = A^B ->
 2928            Goal = ( ( integer(A), integer(B), integer(M),
 2929                       A >= 0, B >= 0, M > 0 ->
 2930                       Var is powm(A, B, M)
 2931                     ; Var is Expr
 2932                     ) )
 2933        ;   Goal = ( Var is Expr )
 2934        ).
 2935expansion_simpler(between(L,U,V), Goal) :- maplist(integer, [L,U,V]), !,
 2936        (   between(L,U,V) -> Goal = true
 2937        ;   Goal = false
 2938        ).
 2939expansion_simpler(Goal, Goal).
 2940
 2941is_true(true).
 2942is_true(integer(I))  :- integer(I).
 2943:- if(current_predicate(var_property/2)). 2944is_true(var(X))      :- var(X), var_property(X, fresh(true)).
 2945is_false(integer(X)) :- var(X), var_property(X, fresh(true)).
 2946is_false((A,B))      :- is_false(A) ; is_false(B).
 2947:- endif. 2948is_false(var(X)) :- nonvar(X).
 2949
 2950
 2951%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 2952
 2953linsum(X, S, S)    --> { var(X), !, non_monotonic(X) }, [vn(X,1)].
 2954linsum(I, S0, S)   --> { integer(I), S is S0 + I }.
 2955linsum(?(X), S, S) --> { must_be_fd_integer(X) }, [vn(X,1)].
 2956linsum(#(X), S, S) --> { must_be_fd_integer(X) }, [vn(X,1)].
 2957linsum(-A, S0, S)  --> mulsum(A, -1, S0, S).
 2958linsum(N*A, S0, S) --> { integer(N) }, !, mulsum(A, N, S0, S).
 2959linsum(A*N, S0, S) --> { integer(N) }, !, mulsum(A, N, S0, S).
 2960linsum(A+B, S0, S) --> linsum(A, S0, S1), linsum(B, S1, S).
 2961linsum(A-B, S0, S) --> linsum(A, S0, S1), mulsum(B, -1, S1, S).
 2962
 2963mulsum(A, M, S0, S) -->
 2964        { phrase(linsum(A, 0, CA), As), S is S0 + M*CA },
 2965        lin_mul(As, M).
 2966
 2967lin_mul([], _)             --> [].
 2968lin_mul([vn(X,N0)|VNs], M) --> { N is N0*M }, [vn(X,N)], lin_mul(VNs, M).
 2969
 2970v_or_i(V) :- var(V), !, non_monotonic(V).
 2971v_or_i(I) :- integer(I).
 2972
 2973must_be_fd_integer(X) :-
 2974        (   var(X) -> constrain_to_integer(X)
 2975        ;   must_be(integer, X)
 2976        ).
 2977
 2978left_right_linsum_const(Left, Right, Cs, Vs, Const) :-
 2979        phrase(linsum(Left, 0, CL), Lefts0, Rights),
 2980        phrase(linsum(Right, 0, CR), Rights0),
 2981        maplist(linterm_negate, Rights0, Rights),
 2982        msort(Lefts0, Lefts),
 2983        Lefts = [vn(First,N)|LeftsRest],
 2984        vns_coeffs_variables(LeftsRest, N, First, Cs0, Vs0),
 2985        filter_linsum(Cs0, Vs0, Cs, Vs),
 2986        Const is CR - CL.
 2987        %format("linear sum: ~w ~w ~w\n", [Cs,Vs,Const]).
 2988
 2989linterm_negate(vn(V,N0), vn(V,N)) :- N is -N0.
 2990
 2991vns_coeffs_variables([], N, V, [N], [V]).
 2992vns_coeffs_variables([vn(V,N)|VNs], N0, V0, Ns, Vs) :-
 2993        (   V == V0 ->
 2994            N1 is N0 + N,
 2995            vns_coeffs_variables(VNs, N1, V0, Ns, Vs)
 2996        ;   Ns = [N0|NRest],
 2997            Vs = [V0|VRest],
 2998            vns_coeffs_variables(VNs, N, V, NRest, VRest)
 2999        ).
 3000
 3001filter_linsum([], [], [], []).
 3002filter_linsum([C0|Cs0], [V0|Vs0], Cs, Vs) :-
 3003        (   C0 =:= 0 ->
 3004            constrain_to_integer(V0),
 3005            filter_linsum(Cs0, Vs0, Cs, Vs)
 3006        ;   Cs = [C0|Cs1], Vs = [V0|Vs1],
 3007            filter_linsum(Cs0, Vs0, Cs1, Vs1)
 3008        ).
 3009
 3010gcd([], G, G).
 3011gcd([N|Ns], G0, G) :-
 3012        G1 is gcd(N, G0),
 3013        gcd(Ns, G1, G).
 3014
 3015even(N) :- N mod 2 =:= 0.
 3016
 3017odd(N) :- \+ even(N).
 3018
 3019/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 3020   k-th root of N, if N is a k-th power.
 3021
 3022   TODO: Replace this when the GMP function becomes available.
 3023- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 3024
 3025integer_kth_root(N, K, R) :-
 3026        (   even(K) ->
 3027            N >= 0
 3028        ;   true
 3029        ),
 3030        (   N < 0 ->
 3031            odd(K),
 3032            integer_kroot(N, 0, N, K, R)
 3033        ;   integer_kroot(0, N, N, K, R)
 3034        ).
 3035
 3036integer_kroot(L, U, N, K, R) :-
 3037        (   L =:= U -> N =:= L^K, R = L
 3038        ;   L + 1 =:= U ->
 3039            (   L^K =:= N -> R = L
 3040            ;   U^K =:= N -> R = U
 3041            ;   false
 3042            )
 3043        ;   Mid is (L + U)//2,
 3044            (   Mid^K > N ->
 3045                integer_kroot(L, Mid, N, K, R)
 3046            ;   integer_kroot(Mid, U, N, K, R)
 3047            )
 3048        ).
 3049
 3050integer_log_b(N, B, Log0, Log) :-
 3051        T is B^Log0,
 3052        (   T =:= N -> Log = Log0
 3053        ;   T < N,
 3054            Log1 is Log0 + 1,
 3055            integer_log_b(N, B, Log1, Log)
 3056        ).
 3057
 3058floor_integer_log_b(N, B, Log0, Log) :-
 3059        T is B^Log0,
 3060        (   T > N -> Log is Log0 - 1
 3061        ;   T =:= N -> Log = Log0
 3062        ;   T < N,
 3063            Log1 is Log0 + 1,
 3064            floor_integer_log_b(N, B, Log1, Log)
 3065        ).
 3066
 3067
 3068/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 3069   Largest R such that R^K =< N.
 3070- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 3071
 3072:- if(current_predicate(nth_integer_root_and_remainder/4)). 3073
 3074% This currently only works for K >= 1, which is all that is needed for now.
 3075integer_kth_root_leq(N, K, R) :-
 3076        nth_integer_root_and_remainder(K, N, R, _).
 3077
 3078:- else. 3079
 3080integer_kth_root_leq(N, K, R) :-
 3081        (   even(K) ->
 3082            N >= 0
 3083        ;   true
 3084        ),
 3085        (   N < 0 ->
 3086            odd(K),
 3087            integer_kroot_leq(N, 0, N, K, R)
 3088        ;   integer_kroot_leq(0, N, N, K, R)
 3089        ).
 3090
 3091integer_kroot_leq(L, U, N, K, R) :-
 3092        (   L =:= U -> R = L
 3093        ;   L + 1 =:= U ->
 3094            (   U^K =< N -> R = U
 3095            ;   R = L
 3096            )
 3097        ;   Mid is (L + U)//2,
 3098            (   Mid^K > N ->
 3099                integer_kroot_leq(L, Mid, N, K, R)
 3100            ;   integer_kroot_leq(Mid, U, N, K, R)
 3101            )
 3102        ).
 3103
 3104:- endif.
 #\=(?X, ?Y)
The arithmetic expressions X and Y evaluate to distinct integers. When reasoning over integers, replace (=\=)/2 by #\=/2 to obtain more general relations. See declarative integer arithmetic.
 3113X #\= Y :- clpfd_neq(X, Y), do_queue.
 3114
 3115% X #\= Y + Z
 3116
 3117x_neq_y_plus_z(X, Y, Z) :-
 3118        propagator_init_trigger(x_neq_y_plus_z(X,Y,Z)).
 3119
 3120% X is distinct from the number N. This is used internally, and does
 3121% not reinforce other constraints.
 3122
 3123neq_num(X, N) :-
 3124        (   fd_get(X, XD, XPs) ->
 3125            domain_remove(XD, N, XD1),
 3126            fd_put(X, XD1, XPs)
 3127        ;   X =\= N
 3128        ).
 #>(?X, ?Y)
Same as Y #< X. When reasoning over integers, replace (>)/2 by #>/2 to obtain more general relations See declarative integer arithmetic.
 3136X #> Y  :- X #>= Y + 1.
 #<(?X, ?Y)
The arithmetic expression X is less than Y. When reasoning over integers, replace (<)/2 by #</2 to obtain more general relations. See declarative integer arithmetic.

In addition to its regular use in tasks that require it, this constraint can also be useful to eliminate uninteresting symmetries from a problem. For example, all possible matches between pairs built from four players in total:

?- Vs = [A,B,C,D], Vs ins 1..4,
        all_different(Vs),
        A #< B, C #< D, A #< C,
   findall(pair(A,B)-pair(C,D), label(Vs), Ms).
Ms = [ pair(1, 2)-pair(3, 4),
       pair(1, 3)-pair(2, 4),
       pair(1, 4)-pair(2, 3)].
 3159X #< Y  :- Y #> X.
 #\(+Q)
Q does not hold. See reification.

For example, to obtain the complement of a domain:

?- #\ X in -3..0\/10..80.
X in inf.. -4\/1..9\/81..sup.
 3172#\ Q       :- reify(Q, 0), do_queue.
 #<==>(?P, ?Q)
P and Q are equivalent. See reification.

For example:

?- X #= 4 #<==> B, X #\= 4.
B = 0,
X in inf..3\/5..sup.

The following example uses reified constraints to relate a list of finite domain variables to the number of occurrences of a given value:

vs_n_num(Vs, N, Num) :-
        maplist(eq_b(N), Vs, Bs),
        sum(Bs, #=, Num).

eq_b(X, Y, B) :- X #= Y #<==> B.

Sample queries and their results:

?- Vs = [X,Y,Z], Vs ins 0..1, vs_n_num(Vs, 4, Num).
Vs = [X, Y, Z],
Num = 0,
X in 0..1,
Y in 0..1,
Z in 0..1.

?- vs_n_num([X,Y,Z], 2, 3).
X = 2,
Y = 2,
Z = 2.
 3212L #<==> R  :- reify(L, B), reify(R, B), do_queue.
 #==>(?P, ?Q)
P implies Q. See reification.
 3218/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 3219   Implication is special in that created auxiliary constraints can be
 3220   retracted when the implication becomes entailed, for example:
 3221
 3222   %?- X + 1 #= Y #==> Z, Z #= 1.
 3223   %@ Z = 1,
 3224   %@ X in inf..sup,
 3225   %@ Y in inf..sup.
 3226
 3227   We cannot use propagator_init_trigger/1 here because the states of
 3228   auxiliary propagators are themselves part of the propagator.
 3229- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 3230
 3231L #==> R   :-
 3232        reify(L, LB, LPs),
 3233        reify(R, RB, RPs),
 3234        append(LPs, RPs, Ps),
 3235        propagator_init_trigger([LB,RB], pimpl(LB,RB,Ps)).
 #<==(?P, ?Q)
Q implies P. See reification.
 3241L #<== R   :- R #==> L.
 #/\(?P, ?Q)
P and Q hold. See reification.
 3247L #/\ R    :- reify(L, 1), reify(R, 1), do_queue.
 3248
 3249conjunctive_neqs_var_drep(Eqs, Var, Drep) :-
 3250        conjunctive_neqs_var(Eqs, Var),
 3251        phrase(conjunctive_neqs_vals(Eqs), Vals),
 3252        list_to_domain(Vals, Dom),
 3253        domain_complement(Dom, C),
 3254        domain_to_drep(C, Drep).
 3255
 3256conjunctive_neqs_var(V, _) :- var(V), !, false.
 3257conjunctive_neqs_var(L #\= R, Var) :-
 3258        (   var(L), integer(R) -> Var = L
 3259        ;   integer(L), var(R) -> Var = R
 3260        ;   false
 3261        ).
 3262conjunctive_neqs_var(A #/\ B, VA) :-
 3263        conjunctive_neqs_var(A, VA),
 3264        conjunctive_neqs_var(B, VB),
 3265        VA == VB.
 3266
 3267conjunctive_neqs_vals(L #\= R) --> ( { integer(L) } -> [L] ; [R] ).
 3268conjunctive_neqs_vals(A #/\ B) -->
 3269        conjunctive_neqs_vals(A),
 3270        conjunctive_neqs_vals(B).
 #\/(?P, ?Q)
P or Q holds. See reification.

For example, the sum of natural numbers below 1000 that are multiples of 3 or 5:

?- findall(N, (N mod 3 #= 0 #\/ N mod 5 #= 0, N in 0..999,
               indomain(N)),
           Ns),
   sum(Ns, #=, Sum).
Ns = [0, 3, 5, 6, 9, 10, 12, 15, 18|...],
Sum = 233168.
 3288L #\/ R :-
 3289        (   disjunctive_eqs_var_drep(L #\/ R, Var, Drep) -> Var in Drep
 3290        ;   reify(L, X, Ps1),
 3291            reify(R, Y, Ps2),
 3292            propagator_init_trigger([X,Y], reified_or(X,Ps1,Y,Ps2,1))
 3293        ).
 3294
 3295disjunctive_eqs_var_drep(Eqs, Var, Drep) :-
 3296        disjunctive_eqs_var(Eqs, Var),
 3297        phrase(disjunctive_eqs_vals(Eqs), Vals),
 3298        list_to_drep(Vals, Drep).
 3299
 3300disjunctive_eqs_var(V, _) :- var(V), !, false.
 3301disjunctive_eqs_var(V in I, V) :- var(V), integer(I).
 3302disjunctive_eqs_var(L #= R, Var) :-
 3303        (   var(L), integer(R) -> Var = L
 3304        ;   integer(L), var(R) -> Var = R
 3305        ;   false
 3306        ).
 3307disjunctive_eqs_var(A #\/ B, VA) :-
 3308        disjunctive_eqs_var(A, VA),
 3309        disjunctive_eqs_var(B, VB),
 3310        VA == VB.
 3311
 3312disjunctive_eqs_vals(L #= R)  --> ( { integer(L) } -> [L] ; [R] ).
 3313disjunctive_eqs_vals(_ in I)  --> [I].
 3314disjunctive_eqs_vals(A #\/ B) -->
 3315        disjunctive_eqs_vals(A),
 3316        disjunctive_eqs_vals(B).
 #\(?P, ?Q)
Either P holds or Q holds, but not both. See reification.
 3323L #\ R :- (L #\/ R) #/\ #\ (L #/\ R).
 3324
 3325/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 3326   A constraint that is being reified need not hold. Therefore, in
 3327   X/Y, Y can as well be 0, for example. Note that it is OK to
 3328   constrain the *result* of an expression (which does not appear
 3329   explicitly in the expression and is not visible to the outside),
 3330   but not the operands, except for requiring that they be integers.
 3331
 3332   In contrast to parse_clpfd/2, the result of an expression can now
 3333   also be undefined, in which case the constraint cannot hold.
 3334   Therefore, the committed-choice language is extended by an element
 3335   d(D) that states D is 1 iff all subexpressions are defined. a(V)
 3336   means that V is an auxiliary variable that was introduced while
 3337   parsing a compound expression. a(X,V) means V is auxiliary unless
 3338   it is ==/2 X, and a(X,Y,V) means V is auxiliary unless it is ==/2 X
 3339   or Y. l(L) means the literal L occurs in the described list.
 3340
 3341   When a constraint becomes entailed or subexpressions become
 3342   undefined, created auxiliary constraints are killed, and the
 3343   "clpfd" attribute is removed from auxiliary variables.
 3344
 3345   For (/)/2, mod/2 and rem/2, we create a skeleton propagator and
 3346   remember it as an auxiliary constraint. The pskeleton propagator
 3347   can use the skeleton when the constraint is defined.
 3348- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 3349
 3350parse_reified(E, R, D,
 3351              [g(cyclic_term(E)) => [g(domain_error(clpfd_expression, E))],
 3352               g(var(E))     => [g(non_monotonic(E)),
 3353                                 g(constrain_to_integer(E)), g(R = E), g(D=1)],
 3354               g(integer(E)) => [g(R=E), g(D=1)],
 3355               ?(E)          => [g(must_be_fd_integer(E)), g(R=E), g(D=1)],
 3356               #(E)          => [g(must_be_fd_integer(E)), g(R=E), g(D=1)],
 3357               m(A+B)        => [d(D), p(pplus(A,B,R)), a(A,B,R)],
 3358               m(A*B)        => [d(D), p(ptimes(A,B,R)), a(A,B,R)],
 3359               m(A-B)        => [d(D), p(pplus(R,B,A)), a(A,B,R)],
 3360               m(-A)         => [d(D), p(ptimes(-1,A,R)), a(R)],
 3361               m(max(A,B))   => [d(D), p(pgeq(R, A)), p(pgeq(R, B)), p(pmax(A,B,R)), a(A,B,R)],
 3362               m(min(A,B))   => [d(D), p(pgeq(A, R)), p(pgeq(B, R)), p(pmin(A,B,R)), a(A,B,R)],
 3363               m(abs(A))     => [g(?(R)#>=0), d(D), p(pabs(A, R)), a(A,R)],
 3364%               m(A/B)        => [skeleton(A,B,D,R,ptzdiv)],
 3365               m(A//B)       => [skeleton(A,B,D,R,ptzdiv)],
 3366               m(A div B)    => [skeleton(A,B,D,R,pdiv)],
 3367               m(A rdiv B)   => [skeleton(A,B,D,R,prdiv)],
 3368               m(A mod B)    => [skeleton(A,B,D,R,pmod)],
 3369               m(A rem B)    => [skeleton(A,B,D,R,prem)],
 3370               m(A^B)        => [d(D), p(pexp(A,B,R)), a(A,B,R)],
 3371               % bitwise operations
 3372               m(\A)         => [function(D,\,A,R)],
 3373               m(msb(A))     => [function(D,msb,A,R)],
 3374               m(lsb(A))     => [function(D,lsb,A,R)],
 3375               m(popcount(A)) => [function(D,popcount,A,R)],
 3376               m(A<<B)       => [function(D,<<,A,B,R)],
 3377               m(A>>B)       => [function(D,>>,A,B,R)],
 3378               m(A/\B)       => [function(D,/\,A,B,R)],
 3379               m(A\/B)       => [function(D,\/,A,B,R)],
 3380               m(A xor B)    => [function(D,xor,A,B,R)],
 3381               g(true)       => [g(domain_error(clpfd_expression, E))]]
 3382             ).
 3383
 3384% Again, we compile this to a predicate, parse_reified_clpfd//3. This
 3385% time, it is a DCG that describes the list of auxiliary variables and
 3386% propagators for the given expression, in addition to relating it to
 3387% its reified (Boolean) finite domain variable and its Boolean
 3388% definedness.
 3389
 3390make_parse_reified(Clauses) :-
 3391        parse_reified_clauses(Clauses0),
 3392        maplist(goals_goal_dcg, Clauses0, Clauses).
 3393
 3394goals_goal_dcg((Head --> Goals), Clause) :-
 3395        list_goal(Goals, Body),
 3396        expand_term((Head --> Body), Clause).
 3397
 3398parse_reified_clauses(Clauses) :-
 3399        parse_reified(E, R, D, Matchers),
 3400        maplist(parse_reified(E, R, D), Matchers, Clauses).
 3401
 3402parse_reified(E, R, D, Matcher, Clause) :-
 3403        Matcher = (Condition0 => Goals0),
 3404        phrase((reified_condition(Condition0, E, Head, Ds),
 3405                reified_goals(Goals0, Ds)), Goals, [a(D)]),
 3406        Clause = (parse_reified_clpfd(Head, R, D) --> Goals).
 3407
 3408reified_condition(g(Goal), E, E, []) --> [{Goal}, !].
 3409reified_condition(?(E), _, ?(E), []) --> [!].
 3410reified_condition(#(E), _, #(E), []) --> [!].
 3411reified_condition(m(Match), _, Match0, Ds) -->
 3412        [!],
 3413        { copy_term(Match, Match0),
 3414          term_variables(Match0, Vs0),
 3415          term_variables(Match, Vs)
 3416        },
 3417        reified_variables(Vs0, Vs, Ds).
 3418
 3419reified_variables([], [], []) --> [].
 3420reified_variables([V0|Vs0], [V|Vs], [D|Ds]) -->
 3421        [parse_reified_clpfd(V0, V, D)],
 3422        reified_variables(Vs0, Vs, Ds).
 3423
 3424reified_goals([], _) --> [].
 3425reified_goals([G|Gs], Ds) --> reified_goal(G, Ds), reified_goals(Gs, Ds).
 3426
 3427reified_goal(d(D), Ds) -->
 3428        (   { Ds = [X] } -> [{D=X}]
 3429        ;   { Ds = [X,Y] } ->
 3430            { phrase(reified_goal(p(reified_and(X,[],Y,[],D)), _), Gs),
 3431              list_goal(Gs, Goal) },
 3432            [( {X==1, Y==1} -> {D = 1} ; Goal )]
 3433        ;   { domain_error(one_or_two_element_list, Ds) }
 3434        ).
 3435reified_goal(g(Goal), _) --> [{Goal}].
 3436reified_goal(p(Vs, Prop), _) -->
 3437        [{make_propagator(Prop, P)}],
 3438        parse_init_dcg(Vs, P),
 3439        [{trigger_once(P)}],
 3440        [( { propagator_state(P, S), S == dead } -> [] ; [p(P)])].
 3441reified_goal(p(Prop), Ds) -->
 3442        { term_variables(Prop, Vs) },
 3443        reified_goal(p(Vs,Prop), Ds).
 3444reified_goal(function(D,Op,A,B,R), Ds) -->
 3445        reified_goals([d(D),p(pfunction(Op,A,B,R)),a(A,B,R)], Ds).
 3446reified_goal(function(D,Op,A,R), Ds) -->
 3447        reified_goals([d(D),p(pfunction(Op,A,R)),a(A,R)], Ds).
 3448reified_goal(skeleton(A,B,D,R,F), Ds) -->
 3449        { Prop =.. [F,X,Y,Z] },
 3450        reified_goals([d(D1),l(p(P)),g(make_propagator(Prop, P)),
 3451                       p([A,B,D2,R], pskeleton(A,B,D2,[X,Y,Z]-P,R,F)),
 3452                       p(reified_and(D1,[],D2,[],D)),a(D2),a(A,B,R)], Ds).
 3453reified_goal(a(V), _)     --> [a(V)].
 3454reified_goal(a(X,V), _)   --> [a(X,V)].
 3455reified_goal(a(X,Y,V), _) --> [a(X,Y,V)].
 3456reified_goal(l(L), _)     --> [[L]].
 3457
 3458parse_init_dcg([], _)     --> [].
 3459parse_init_dcg([V|Vs], P) --> [{init_propagator(V, P)}], parse_init_dcg(Vs, P).
 3460
 3461%?- set_prolog_flag(answer_write_options, [portray(true)]),
 3462%   clpfd:parse_reified_clauses(Cs), maplist(portray_clause, Cs).
 3463
 3464reify(E, B) :- reify(E, B, _).
 3465
 3466reify(Expr, B, Ps) :-
 3467        (   acyclic_term(Expr), reifiable(Expr) -> phrase(reify(Expr, B), Ps)
 3468        ;   domain_error(clpfd_reifiable_expression, Expr)
 3469        ).
 3470
 3471reifiable(E)      :- var(E), non_monotonic(E).
 3472reifiable(E)      :- integer(E), E in 0..1.
 3473reifiable(?(E))   :- must_be_fd_integer(E).
 3474reifiable(#(E))   :- must_be_fd_integer(E).
 3475reifiable(V in _) :- fd_variable(V).
 3476reifiable(Expr)   :-
 3477        Expr =.. [Op,Left,Right],
 3478        (   memberchk(Op, [#>=,#>,#=<,#<,#=,#\=])
 3479        ;   memberchk(Op, [#==>,#<==,#<==>,#/\,#\/,#\]),
 3480            reifiable(Left),
 3481            reifiable(Right)
 3482        ).
 3483reifiable(#\ E) :- reifiable(E).
 3484reifiable(tuples_in(Tuples, Relation)) :-
 3485        must_be(list(list), Tuples),
 3486        maplist(maplist(fd_variable), Tuples),
 3487        must_be(list(list(integer)), Relation).
 3488reifiable(finite_domain(V)) :- fd_variable(V).
 3489
 3490reify(E, B) --> { B in 0..1 }, reify_(E, B).
 3491
 3492reify_(E, B) --> { var(E), !, E = B }.
 3493reify_(E, B) --> { integer(E), E = B }.
 3494reify_(?(B), B) --> [].
 3495reify_(#(B), B) --> [].
 3496reify_(V in Drep, B) -->
 3497        { drep_to_domain(Drep, Dom) },
 3498        propagator_init_trigger(reified_in(V,Dom,B)),
 3499        a(B).
 3500reify_(tuples_in(Tuples, Relation), B) -->
 3501        { maplist(relation_tuple_b_prop(Relation), Tuples, Bs, Ps),
 3502          maplist(monotonic, Bs, Bs1),
 3503          fold_statement(conjunction, Bs1, And),
 3504          ?(B) #<==> And },
 3505        propagator_init_trigger([B], tuples_not_in(Tuples, Relation, B)),
 3506        kill_reified_tuples(Bs, Ps, Bs),
 3507        list(Ps),
 3508        as([B|Bs]).
 3509reify_(finite_domain(V), B) -->
 3510        propagator_init_trigger(reified_fd(V,B)),
 3511        a(B).
 3512reify_(L #>= R, B) --> arithmetic(L, R, B, reified_geq).
 3513reify_(L #= R, B)  --> arithmetic(L, R, B, reified_eq).
 3514reify_(L #\= R, B) --> arithmetic(L, R, B, reified_neq).
 3515reify_(L #> R, B)  --> reify_(L #>= (R+1), B).
 3516reify_(L #=< R, B) --> reify_(R #>= L, B).
 3517reify_(L #< R, B)  --> reify_(R #>= (L+1), B).
 3518reify_(L #==> R, B)  --> reify_((#\ L) #\/ R, B).
 3519reify_(L #<== R, B)  --> reify_(R #==> L, B).
 3520reify_(L #<==> R, B) --> reify_((L #==> R) #/\ (R #==> L), B).
 3521reify_(L #\ R, B) --> reify_((L #\/ R) #/\ #\ (L #/\ R), B).
 3522reify_(L #/\ R, B)   -->
 3523        (   { conjunctive_neqs_var_drep(L #/\ R, V, D) } -> reify_(V in D, B)
 3524        ;   boolean(L, R, B, reified_and)
 3525        ).
 3526reify_(L #\/ R, B) -->
 3527        (   { disjunctive_eqs_var_drep(L #\/ R, V, D) } -> reify_(V in D, B)
 3528        ;   boolean(L, R, B, reified_or)
 3529        ).
 3530reify_(#\ Q, B) -->
 3531        reify(Q, QR),
 3532        propagator_init_trigger(reified_not(QR,B)),
 3533        a(B).
 3534
 3535arithmetic(L, R, B, Functor) -->
 3536        { phrase((parse_reified_clpfd(L, LR, LD),
 3537                  parse_reified_clpfd(R, RR, RD)), Ps),
 3538          Prop =.. [Functor,LD,LR,RD,RR,Ps,B] },
 3539        list(Ps),
 3540        propagator_init_trigger([LD,LR,RD,RR,B], Prop),
 3541        a(B).
 3542
 3543boolean(L, R, B, Functor) -->
 3544        { reify(L, LR, Ps1), reify(R, RR, Ps2),
 3545          Prop =.. [Functor,LR,Ps1,RR,Ps2,B] },
 3546        list(Ps1), list(Ps2),
 3547        propagator_init_trigger([LR,RR,B], Prop),
 3548        a(LR, RR, B).
 3549
 3550list([])     --> [].
 3551list([L|Ls]) --> [L], list(Ls).
 3552
 3553a(X,Y,B) -->
 3554        (   { nonvar(X) } -> a(Y, B)
 3555        ;   { nonvar(Y) } -> a(X, B)
 3556        ;   [a(X,Y,B)]
 3557        ).
 3558
 3559a(X, B) -->
 3560        (   { var(X) } -> [a(X, B)]
 3561        ;   a(B)
 3562        ).
 3563
 3564a(B) -->
 3565        (   { var(B) } -> [a(B)]
 3566        ;   []
 3567        ).
 3568
 3569as([])     --> [].
 3570as([B|Bs]) --> a(B), as(Bs).
 3571
 3572kill_reified_tuples([], _, _) --> [].
 3573kill_reified_tuples([B|Bs], Ps, All) -->
 3574        propagator_init_trigger([B], kill_reified_tuples(B, Ps, All)),
 3575        kill_reified_tuples(Bs, Ps, All).
 3576
 3577relation_tuple_b_prop(Relation, Tuple, B, p(Prop)) :-
 3578        put_attr(R, clpfd_relation, Relation),
 3579        make_propagator(reified_tuple_in(Tuple, R, B), Prop),
 3580        tuple_freeze_(Tuple, Prop),
 3581        init_propagator(B, Prop).
 3582
 3583
 3584tuples_in_conjunction(Tuples, Relation, Conj) :-
 3585        maplist(tuple_in_disjunction(Relation), Tuples, Disjs),
 3586        fold_statement(conjunction, Disjs, Conj).
 3587
 3588tuple_in_disjunction(Relation, Tuple, Disj) :-
 3589        maplist(tuple_in_conjunction(Tuple), Relation, Conjs),
 3590        fold_statement(disjunction, Conjs, Disj).
 3591
 3592tuple_in_conjunction(Tuple, Element, Conj) :-
 3593        maplist(var_eq, Tuple, Element, Eqs),
 3594        fold_statement(conjunction, Eqs, Conj).
 3595
 3596fold_statement(Operation, List, Statement) :-
 3597        (   List = [] -> Statement = 1
 3598        ;   List = [First|Rest],
 3599            foldl(Operation, Rest, First, Statement)
 3600        ).
 3601
 3602conjunction(E, Conj, Conj #/\ E).
 3603
 3604disjunction(E, Disj, Disj #\/ E).
 3605
 3606var_eq(V, N, ?(V) #= N).
 3607
 3608% Match variables to created skeleton.
 3609
 3610skeleton(Vs, Vs-Prop) :-
 3611        maplist(prop_init(Prop), Vs),
 3612        trigger_once(Prop).
 3613
 3614/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 3615   A drep is a user-accessible and visible domain representation. N,
 3616   N..M, and D1 \/ D2 are dreps, if D1 and D2 are dreps.
 3617- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 3618
 3619is_drep(N)      :- integer(N).
 3620is_drep(N..M)   :- drep_bound(N), drep_bound(M), N \== sup, M \== inf.
 3621is_drep(D1\/D2) :- is_drep(D1), is_drep(D2).
 3622is_drep({AI})   :- is_and_integers(AI).
 3623is_drep(\D)     :- is_drep(D).
 3624
 3625is_and_integers(I)     :- integer(I).
 3626is_and_integers((A,B)) :- is_and_integers(A), is_and_integers(B).
 3627
 3628drep_bound(I)   :- integer(I).
 3629drep_bound(sup).
 3630drep_bound(inf).
 3631
 3632drep_to_intervals(I)        --> { integer(I) }, [n(I)-n(I)].
 3633drep_to_intervals(N..M)     -->
 3634        (   { defaulty_to_bound(N, N1), defaulty_to_bound(M, M1),
 3635              N1 cis_leq M1} -> [N1-M1]
 3636        ;   []
 3637        ).
 3638drep_to_intervals(D1 \/ D2) -->
 3639        drep_to_intervals(D1), drep_to_intervals(D2).
 3640drep_to_intervals(\D0) -->
 3641        { drep_to_domain(D0, D1),
 3642          domain_complement(D1, D),
 3643          domain_to_drep(D, Drep) },
 3644        drep_to_intervals(Drep).
 3645drep_to_intervals({AI}) -->
 3646        and_integers_(AI).
 3647
 3648and_integers_(I)     --> { integer(I) }, [n(I)-n(I)].
 3649and_integers_((A,B)) --> and_integers_(A), and_integers_(B).
 3650
 3651drep_to_domain(DR, D) :-
 3652        must_be(ground, DR),
 3653        (   is_drep(DR) -> true
 3654        ;   domain_error(clpfd_domain, DR)
 3655        ),
 3656        phrase(drep_to_intervals(DR), Is0),
 3657        merge_intervals(Is0, Is1),
 3658        intervals_to_domain(Is1, D).
 3659
 3660merge_intervals(Is0, Is) :-
 3661        keysort(Is0, Is1),
 3662        merge_overlapping(Is1, Is).
 3663
 3664merge_overlapping([], []).
 3665merge_overlapping([A-B0|ABs0], [A-B|ABs]) :-
 3666        merge_remaining(ABs0, B0, B, Rest),
 3667        merge_overlapping(Rest, ABs).
 3668
 3669merge_remaining([], B, B, []).
 3670merge_remaining([N-M|NMs], B0, B, Rest) :-
 3671        Next cis B0 + n(1),
 3672        (   N cis_gt Next -> B = B0, Rest = [N-M|NMs]
 3673        ;   B1 cis max(B0,M),
 3674            merge_remaining(NMs, B1, B, Rest)
 3675        ).
 3676
 3677domain(V, Dom) :-
 3678        (   fd_get(V, Dom0, VPs) ->
 3679            domains_intersection(Dom, Dom0, Dom1),
 3680            %format("intersected\n: ~w\n ~w\n==> ~w\n\n", [Dom,Dom0,Dom1]),
 3681            fd_put(V, Dom1, VPs),
 3682            do_queue,
 3683            reinforce(V)
 3684        ;   domain_contains(Dom, V)
 3685        ).
 3686
 3687domains([], _).
 3688domains([V|Vs], D) :- domain(V, D), domains(Vs, D).
 3689
 3690props_number(fd_props(Gs,Bs,Os), N) :-
 3691        length(Gs, N1),
 3692        length(Bs, N2),
 3693        length(Os, N3),
 3694        N is N1 + N2 + N3.
 3695
 3696fd_get(X, Dom, Ps) :-
 3697        (   get_attr(X, clpfd, Attr) -> Attr = clpfd_attr(_,_,_,Dom,Ps)
 3698        ;   var(X) -> default_domain(Dom), Ps = fd_props([],[],[])
 3699        ).
 3700
 3701fd_get(X, Dom, Inf, Sup, Ps) :-
 3702        fd_get(X, Dom, Ps),
 3703        domain_infimum(Dom, Inf),
 3704        domain_supremum(Dom, Sup).
 3705
 3706/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 3707   By default, propagation always terminates. Currently, this is
 3708   ensured by allowing the left and right boundaries, as well as the
 3709   distance between the smallest and largest number occurring in the
 3710   domain representation to be changed at most once after a constraint
 3711   is posted, unless the domain is bounded. Set the experimental
 3712   Prolog flag 'clpfd_propagation' to 'full' to make the solver
 3713   propagate as much as possible. This can make queries
 3714   non-terminating, like: X #> abs(X), or: X #> Y, Y #> X, X #> 0.
 3715   Importantly, it can also make labeling non-terminating, as in:
 3716
 3717   ?- B #==> X #> abs(X), indomain(B).
 3718- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 3719
 3720fd_put(X, Dom, Ps) :-
 3721        (   current_prolog_flag(clpfd_propagation, full) ->
 3722            put_full(X, Dom, Ps)
 3723        ;   put_terminating(X, Dom, Ps)
 3724        ).
 3725
 3726put_terminating(X, Dom, Ps) :-
 3727        Dom \== empty,
 3728        (   Dom = from_to(F, F) -> F = n(X)
 3729        ;   (   get_attr(X, clpfd, Attr) ->
 3730                Attr = clpfd_attr(Left,Right,Spread,OldDom, _OldPs),
 3731                put_attr(X, clpfd, clpfd_attr(Left,Right,Spread,Dom,Ps)),
 3732                (   OldDom == Dom -> true
 3733                ;   (   Left == (.) -> Bounded = yes
 3734                    ;   domain_infimum(Dom, Inf), domain_supremum(Dom, Sup),
 3735                        (   Inf = n(_), Sup = n(_) ->
 3736                            Bounded = yes
 3737                        ;   Bounded = no
 3738                        )
 3739                    ),
 3740                    (   Bounded == yes ->
 3741                        put_attr(X, clpfd, clpfd_attr(.,.,.,Dom,Ps)),
 3742                        trigger_props(Ps, X, OldDom, Dom)
 3743                    ;   % infinite domain; consider border and spread changes
 3744                        domain_infimum(OldDom, OldInf),
 3745                        (   Inf == OldInf -> LeftP = Left
 3746                        ;   LeftP = yes
 3747                        ),
 3748                        domain_supremum(OldDom, OldSup),
 3749                        (   Sup == OldSup -> RightP = Right
 3750                        ;   RightP = yes
 3751                        ),
 3752                        domain_spread(OldDom, OldSpread),
 3753                        domain_spread(Dom, NewSpread),
 3754                        (   NewSpread == OldSpread -> SpreadP = Spread
 3755                        ;   NewSpread cis_lt OldSpread -> SpreadP = no
 3756                        ;   SpreadP = yes
 3757                        ),
 3758                        put_attr(X, clpfd, clpfd_attr(LeftP,RightP,SpreadP,Dom,Ps)),
 3759                        (   RightP == yes, Right = yes -> true
 3760                        ;   LeftP == yes, Left = yes -> true
 3761                        ;   SpreadP == yes, Spread = yes -> true
 3762                        ;   trigger_props(Ps, X, OldDom, Dom)
 3763                        )
 3764                    )
 3765                )
 3766            ;   var(X) ->
 3767                put_attr(X, clpfd, clpfd_attr(no,no,no,Dom, Ps))
 3768            ;   true
 3769            )
 3770        ).
 3771
 3772domain_spread(Dom, Spread) :-
 3773        domain_smallest_finite(Dom, S),
 3774        domain_largest_finite(Dom, L),
 3775        Spread cis L - S.
 3776
 3777smallest_finite(inf, Y, Y).
 3778smallest_finite(n(N), _, n(N)).
 3779
 3780domain_smallest_finite(from_to(F,T), S)   :- smallest_finite(F, T, S).
 3781domain_smallest_finite(split(_, L, _), S) :- domain_smallest_finite(L, S).
 3782
 3783largest_finite(sup, Y, Y).
 3784largest_finite(n(N), _, n(N)).
 3785
 3786domain_largest_finite(from_to(F,T), L)   :- largest_finite(T, F, L).
 3787domain_largest_finite(split(_, _, R), L) :- domain_largest_finite(R, L).
 3788
 3789/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 3790   With terminating propagation, all relevant constraints get a
 3791   propagation opportunity whenever a new constraint is posted.
 3792- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 3793
 3794reinforce(X) :-
 3795        (   current_prolog_flag(clpfd_propagation, full) ->
 3796            % full propagation propagates everything in any case
 3797            true
 3798        ;   term_variables(X, Vs),
 3799            maplist(reinforce_, Vs),
 3800            do_queue
 3801        ).
 3802
 3803reinforce_(X) :-
 3804        (   fd_var(X), fd_get(X, Dom, Ps) ->
 3805            put_full(X, Dom, Ps)
 3806        ;   true
 3807        ).
 3808
 3809put_full(X, Dom, Ps) :-
 3810        Dom \== empty,
 3811        (   Dom = from_to(F, F) -> F = n(X)
 3812        ;   (   get_attr(X, clpfd, Attr) ->
 3813                Attr = clpfd_attr(_,_,_,OldDom, _OldPs),
 3814                put_attr(X, clpfd, clpfd_attr(no,no,no,Dom, Ps)),
 3815                %format("putting dom: ~w\n", [Dom]),
 3816                (   OldDom == Dom -> true
 3817                ;   trigger_props(Ps, X, OldDom, Dom)
 3818                )
 3819            ;   var(X) -> %format('\t~w in ~w .. ~w\n',[X,L,U]),
 3820                put_attr(X, clpfd, clpfd_attr(no,no,no,Dom, Ps))
 3821            ;   true
 3822            )
 3823        ).
 3824
 3825/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 3826   A propagator is a term of the form propagator(C, State), where C
 3827   represents a constraint, and State is a free variable that can be
 3828   used to destructively change the state of the propagator via
 3829   attributes. This can be used to avoid redundant invocation of the
 3830   same propagator, or to disable the propagator.
 3831- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 3832
 3833make_propagator(C, propagator(C, _)).
 3834
 3835propagator_state(propagator(_,S), S).
 3836
 3837trigger_props(fd_props(Gs,Bs,Os), X, D0, D) :-
 3838        (   ground(X) ->
 3839            trigger_props_(Gs),
 3840            trigger_props_(Bs)
 3841        ;   Bs \== [] ->
 3842            domain_infimum(D0, I0),
 3843            domain_infimum(D, I),
 3844            (   I == I0 ->
 3845                domain_supremum(D0, S0),
 3846                domain_supremum(D, S),
 3847                (   S == S0 -> true
 3848                ;   trigger_props_(Bs)
 3849                )
 3850            ;   trigger_props_(Bs)
 3851            )
 3852        ;   true
 3853        ),
 3854        trigger_props_(Os).
 3855
 3856trigger_props(fd_props(Gs,Bs,Os), X) :-
 3857        trigger_props_(Os),
 3858        trigger_props_(Bs),
 3859        (   ground(X) ->
 3860            trigger_props_(Gs)
 3861        ;   true
 3862        ).
 3863
 3864trigger_props(fd_props(Gs,Bs,Os)) :-
 3865        trigger_props_(Gs),
 3866        trigger_props_(Bs),
 3867        trigger_props_(Os).
 3868
 3869trigger_props_([]).
 3870trigger_props_([P|Ps]) :- trigger_prop(P), trigger_props_(Ps).
 3871
 3872trigger_prop(Propagator) :-
 3873        propagator_state(Propagator, State),
 3874        (   State == dead -> true
 3875        ;   get_attr(State, clpfd_aux, queued) -> true
 3876        ;   b_getval('$clpfd_current_propagator', C), C == State -> true
 3877        ;   % passive
 3878            % format("triggering: ~w\n", [Propagator]),
 3879            put_attr(State, clpfd_aux, queued),
 3880            (   arg(1, Propagator, C), functor(C, F, _), global_constraint(F) ->
 3881                push_queue(Propagator, 2)
 3882            ;   push_queue(Propagator, 1)
 3883            )
 3884        ).
 3885
 3886kill(State) :- del_attr(State, clpfd_aux), State = dead.
 3887
 3888kill(State, Ps) :-
 3889        kill(State),
 3890        maplist(kill_entailed, Ps).
 3891
 3892kill_entailed(p(Prop)) :-
 3893        propagator_state(Prop, State),
 3894        kill(State).
 3895kill_entailed(a(V)) :-
 3896        del_attr(V, clpfd).
 3897kill_entailed(a(X,B)) :-
 3898        (   X == B -> true
 3899        ;   del_attr(B, clpfd)
 3900        ).
 3901kill_entailed(a(X,Y,B)) :-
 3902        (   X == B -> true
 3903        ;   Y == B -> true
 3904        ;   del_attr(B, clpfd)
 3905        ).
 3906
 3907no_reactivation(rel_tuple(_,_)).
 3908no_reactivation(pdistinct(_)).
 3909no_reactivation(pgcc(_,_,_)).
 3910no_reactivation(pgcc_single(_,_)).
 3911%no_reactivation(scalar_product(_,_,_,_)).
 3912
 3913activate_propagator(propagator(P,State)) :-
 3914        % format("running: ~w\n", [P]),
 3915        del_attr(State, clpfd_aux),
 3916        (   no_reactivation(P) ->
 3917            b_setval('$clpfd_current_propagator', State),
 3918            run_propagator(P, State),
 3919            b_setval('$clpfd_current_propagator', [])
 3920        ;   run_propagator(P, State)
 3921        ).
 3922
 3923disable_queue :- b_setval('$clpfd_queue_status', disabled).
 3924enable_queue  :- b_setval('$clpfd_queue_status', enabled).
 3925
 3926portray_propagator(propagator(P,_), F) :- functor(P, F, _).
 3927
 3928portray_queue(V, []) :- var(V), !.
 3929portray_queue([P|Ps], [F|Fs]) :-
 3930        portray_propagator(P, F),
 3931        portray_queue(Ps, Fs).
 3932
 3933do_queue :-
 3934        % b_getval('$clpfd_queue', H-_),
 3935        % portray_queue(H, Port),
 3936        % format("queue: ~w\n", [Port]),
 3937        (   b_getval('$clpfd_queue_status', enabled) ->
 3938            (   fetch_propagator(Propagator) ->
 3939                activate_propagator(Propagator),
 3940                do_queue
 3941            ;   true
 3942            )
 3943        ;   true
 3944        ).
 3945
 3946init_propagator(Var, Prop) :-
 3947        (   fd_get(Var, Dom, Ps0) ->
 3948            insert_propagator(Prop, Ps0, Ps),
 3949            fd_put(Var, Dom, Ps)
 3950        ;   true
 3951        ).
 3952
 3953constraint_wake(pneq, ground).
 3954constraint_wake(x_neq_y_plus_z, ground).
 3955constraint_wake(absdiff_neq, ground).
 3956constraint_wake(pdifferent, ground).
 3957constraint_wake(pexclude, ground).
 3958constraint_wake(scalar_product_neq, ground).
 3959
 3960constraint_wake(x_leq_y_plus_c, bounds).
 3961constraint_wake(scalar_product_eq, bounds).
 3962constraint_wake(scalar_product_leq, bounds).
 3963constraint_wake(pplus, bounds).
 3964constraint_wake(pgeq, bounds).
 3965constraint_wake(pgcc_single, bounds).
 3966constraint_wake(pgcc_check_single, bounds).
 3967
 3968global_constraint(pdistinct).
 3969global_constraint(pgcc).
 3970global_constraint(pgcc_single).
 3971global_constraint(pcircuit).
 3972%global_constraint(rel_tuple).
 3973%global_constraint(scalar_product_eq).
 3974
 3975insert_propagator(Prop, Ps0, Ps) :-
 3976        Ps0 = fd_props(Gs,Bs,Os),
 3977        arg(1, Prop, Constraint),
 3978        functor(Constraint, F, _),
 3979        (   constraint_wake(F, ground) ->
 3980            Ps = fd_props([Prop|Gs], Bs, Os)
 3981        ;   constraint_wake(F, bounds) ->
 3982            Ps = fd_props(Gs, [Prop|Bs], Os)
 3983        ;   Ps = fd_props(Gs, Bs, [Prop|Os])
 3984        ).
 3985
 3986%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 lex_chain(+Lists)
Lists are lexicographically non-decreasing.
 3992lex_chain(Lss) :-
 3993        must_be(list(list), Lss),
 3994        maplist(maplist(fd_variable), Lss),
 3995        (   Lss == [] -> true
 3996        ;   Lss = [First|Rest],
 3997            make_propagator(presidual(lex_chain(Lss)), Prop),
 3998            foldl(lex_chain_(Prop), Rest, First, _)
 3999        ).
 4000
 4001lex_chain_(Prop, Ls, Prev, Ls) :-
 4002        maplist(prop_init(Prop), Ls),
 4003        lex_le(Prev, Ls).
 4004
 4005lex_le([], []).
 4006lex_le([V1|V1s], [V2|V2s]) :-
 4007        ?(V1) #=< ?(V2),
 4008        (   integer(V1) ->
 4009            (   integer(V2) ->
 4010                (   V1 =:= V2 -> lex_le(V1s, V2s) ;  true )
 4011            ;   freeze(V2, lex_le([V1|V1s], [V2|V2s]))
 4012            )
 4013        ;   freeze(V1, lex_le([V1|V1s], [V2|V2s]))
 4014        ).
 4015
 4016%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 tuples_in(+Tuples, +Relation)
True iff all Tuples are elements of Relation. Each element of the list Tuples is a list of integers or finite domain variables. Relation is a list of lists of integers. Arbitrary finite relations, such as compatibility tables, can be modeled in this way. For example, if 1 is compatible with 2 and 5, and 4 is compatible with 0 and 3:
?- tuples_in([[X,Y]], [[1,2],[1,5],[4,0],[4,3]]), X = 4.
X = 4,
Y in 0\/3.

As another example, consider a train schedule represented as a list of quadruples, denoting departure and arrival places and times for each train. In the following program, Ps is a feasible journey of length 3 from A to D via trains that are part of the given schedule.

trains([[1,2,0,1],
        [2,3,4,5],
        [2,3,0,1],
        [3,4,5,6],
        [3,4,2,3],
        [3,4,8,9]]).

threepath(A, D, Ps) :-
        Ps = [[A,B,_T0,T1],[B,C,T2,T3],[C,D,T4,_T5]],
        T2 #> T1,
        T4 #> T3,
        trains(Ts),
        tuples_in(Ps, Ts).

In this example, the unique solution is found without labeling:

?- threepath(1, 4, Ps).
Ps = [[1, 2, 0, 1], [2, 3, 4, 5], [3, 4, 8, 9]].
 4062tuples_in(Tuples, Relation) :-
 4063        must_be(list(list), Tuples),
 4064        maplist(maplist(fd_variable), Tuples),
 4065        must_be(list(list(integer)), Relation),
 4066        maplist(relation_tuple(Relation), Tuples),
 4067        do_queue.
 4068
 4069relation_tuple(Relation, Tuple) :-
 4070        relation_unifiable(Relation, Tuple, Us, _, _),
 4071        (   ground(Tuple) -> memberchk(Tuple, Relation)
 4072        ;   tuple_domain(Tuple, Us),
 4073            (   Tuple = [_,_|_] -> tuple_freeze(Tuple, Us)
 4074            ;   true
 4075            )
 4076        ).
 4077
 4078tuple_domain([], _).
 4079tuple_domain([T|Ts], Relation0) :-
 4080        maplist(list_first_rest, Relation0, Firsts, Relation1),
 4081        (   var(T) ->
 4082            (   Firsts = [Unique] -> T = Unique
 4083            ;   list_to_domain(Firsts, FDom),
 4084                fd_get(T, TDom, TPs),
 4085                domains_intersection(TDom, FDom, TDom1),
 4086                fd_put(T, TDom1, TPs)
 4087            )
 4088        ;   true
 4089        ),
 4090        tuple_domain(Ts, Relation1).
 4091
 4092tuple_freeze(Tuple, Relation) :-
 4093        put_attr(R, clpfd_relation, Relation),
 4094        make_propagator(rel_tuple(R, Tuple), Prop),
 4095        tuple_freeze_(Tuple, Prop).
 4096
 4097tuple_freeze_([], _).
 4098tuple_freeze_([T|Ts], Prop) :-
 4099        (   var(T) ->
 4100            init_propagator(T, Prop),
 4101            trigger_prop(Prop)
 4102        ;   true
 4103        ),
 4104        tuple_freeze_(Ts, Prop).
 4105
 4106relation_unifiable([], _, [], Changed, Changed).
 4107relation_unifiable([R|Rs], Tuple, Us, Changed0, Changed) :-
 4108        (   all_in_domain(R, Tuple) ->
 4109            Us = [R|Rest],
 4110            relation_unifiable(Rs, Tuple, Rest, Changed0, Changed)
 4111        ;   relation_unifiable(Rs, Tuple, Us, true, Changed)
 4112        ).
 4113
 4114all_in_domain([], []).
 4115all_in_domain([A|As], [T|Ts]) :-
 4116        (   fd_get(T, Dom, _) ->
 4117            domain_contains(Dom, A)
 4118        ;   T =:= A
 4119        ),
 4120        all_in_domain(As, Ts).
 4121
 4122%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 4123
 4124% trivial propagator, used only to remember pending constraints
 4125run_propagator(presidual(_), _).
 4126
 4127%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 4128run_propagator(pdifferent(Left,Right,X,_), MState) :-
 4129        run_propagator(pexclude(Left,Right,X), MState).
 4130
 4131run_propagator(weak_distinct(Left,Right,X,_), _MState) :-
 4132        (   ground(X) ->
 4133            disable_queue,
 4134            exclude_fire(Left, Right, X),
 4135            enable_queue
 4136        ;   outof_reducer(Left, Right, X)
 4137            %(   var(X) -> kill_if_isolated(Left, Right, X, MState)
 4138            %;   true
 4139            %)
 4140        ).
 4141
 4142run_propagator(pexclude(Left,Right,X), _) :-
 4143        (   ground(X) ->
 4144            disable_queue,
 4145            exclude_fire(Left, Right, X),
 4146            enable_queue
 4147        ;   true
 4148        ).
 4149
 4150run_propagator(pdistinct(Ls), _MState) :-
 4151        distinct(Ls).
 4152
 4153run_propagator(check_distinct(Left,Right,X), _) :-
 4154        \+ list_contains(Left, X),
 4155        \+ list_contains(Right, X).
 4156
 4157%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 4158
 4159run_propagator(pelement(N, Is, V), MState) :-
 4160        (   fd_get(N, NDom, _) ->
 4161            (   fd_get(V, VDom, VPs) ->
 4162                integers_remaining(Is, 1, NDom, empty, VDom1),
 4163                domains_intersection(VDom, VDom1, VDom2),
 4164                fd_put(V, VDom2, VPs)
 4165            ;   true
 4166            )
 4167        ;   kill(MState), nth1(N, Is, V)
 4168        ).
 4169
 4170%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 4171
 4172run_propagator(pgcc_single(Vs, Pairs), _) :- gcc_global(Vs, Pairs).
 4173
 4174run_propagator(pgcc_check_single(Pairs), _) :- gcc_check(Pairs).
 4175
 4176run_propagator(pgcc_check(Pairs), _) :- gcc_check(Pairs).
 4177
 4178run_propagator(pgcc(Vs, _, Pairs), _) :- gcc_global(Vs, Pairs).
 4179
 4180%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 4181
 4182run_propagator(pcircuit(Vs), _MState) :-
 4183        distinct(Vs),
 4184        propagate_circuit(Vs).
 4185
 4186%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 4187run_propagator(pneq(A, B), MState) :-
 4188        (   nonvar(A) ->
 4189            (   nonvar(B) -> A =\= B, kill(MState)
 4190            ;   fd_get(B, BD0, BExp0),
 4191                domain_remove(BD0, A, BD1),
 4192                kill(MState),
 4193                fd_put(B, BD1, BExp0)
 4194            )
 4195        ;   nonvar(B) -> run_propagator(pneq(B, A), MState)
 4196        ;   A \== B,
 4197            fd_get(A, _, AI, AS, _), fd_get(B, _, BI, BS, _),
 4198            (   AS cis_lt BI -> kill(MState)
 4199            ;   AI cis_gt BS -> kill(MState)
 4200            ;   true
 4201            )
 4202        ).
 4203
 4204%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 4205run_propagator(pgeq(A,B), MState) :-
 4206        (   A == B -> kill(MState)
 4207        ;   nonvar(A) ->
 4208            (   nonvar(B) -> kill(MState), A >= B
 4209            ;   fd_get(B, BD, BPs),
 4210                domain_remove_greater_than(BD, A, BD1),
 4211                kill(MState),
 4212                fd_put(B, BD1, BPs)
 4213            )
 4214        ;   nonvar(B) ->
 4215            fd_get(A, AD, APs),
 4216            domain_remove_smaller_than(AD, B, AD1),
 4217            kill(MState),
 4218            fd_put(A, AD1, APs)
 4219        ;   fd_get(A, AD, AL, AU, APs),
 4220            fd_get(B, _, BL, BU, _),
 4221            AU cis_geq BL,
 4222            (   AL cis_geq BU -> kill(MState)
 4223            ;   AU == BL -> kill(MState), A = B
 4224            ;   NAL cis max(AL,BL),
 4225                domains_intersection(AD, from_to(NAL,AU), NAD),
 4226                fd_put(A, NAD, APs),
 4227                (   fd_get(B, BD2, BL2, BU2, BPs2) ->
 4228                    NBU cis min(BU2, AU),
 4229                    domains_intersection(BD2, from_to(BL2,NBU), NBD),
 4230                    fd_put(B, NBD, BPs2)
 4231                ;   true
 4232                )
 4233            )
 4234        ).
 4235
 4236%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 4237
 4238run_propagator(rel_tuple(R, Tuple), MState) :-
 4239        get_attr(R, clpfd_relation, Relation),
 4240        (   ground(Tuple) -> kill(MState), memberchk(Tuple, Relation)
 4241        ;   relation_unifiable(Relation, Tuple, Us, false, Changed),
 4242            Us = [_|_],
 4243            (   Tuple = [First,Second], ( ground(First) ; ground(Second) ) ->
 4244                kill(MState)
 4245            ;   true
 4246            ),
 4247            (   Us = [Single] -> kill(MState), Single = Tuple
 4248            ;   Changed ->
 4249                put_attr(R, clpfd_relation, Us),
 4250                disable_queue,
 4251                tuple_domain(Tuple, Us),
 4252                enable_queue
 4253            ;   true
 4254            )
 4255        ).
 4256
 4257%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 4258
 4259run_propagator(pserialized(S_I, D_I, S_J, D_J, _), MState) :-
 4260        (   nonvar(S_I), nonvar(S_J) ->
 4261            kill(MState),
 4262            (   S_I + D_I =< S_J -> true
 4263            ;   S_J + D_J =< S_I -> true
 4264            ;   false
 4265            )
 4266        ;   serialize_lower_upper(S_I, D_I, S_J, D_J, MState),
 4267            serialize_lower_upper(S_J, D_J, S_I, D_I, MState)
 4268        ).
 4269
 4270%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 4271
 4272% abs(X-Y) #\= C
 4273run_propagator(absdiff_neq(X,Y,C), MState) :-
 4274        (   C < 0 -> kill(MState)
 4275        ;   nonvar(X) ->
 4276            kill(MState),
 4277            (   nonvar(Y) -> abs(X - Y) =\= C
 4278            ;   V1 is X - C, neq_num(Y, V1),
 4279                V2 is C + X, neq_num(Y, V2)
 4280            )
 4281        ;   nonvar(Y) -> kill(MState),
 4282            V1 is C + Y, neq_num(X, V1),
 4283            V2 is Y - C, neq_num(X, V2)
 4284        ;   true
 4285        ).
 4286
 4287% abs(X-Y) #>= C
 4288run_propagator(absdiff_geq(X,Y,C), MState) :-
 4289        (   C =< 0 -> kill(MState)
 4290        ;   nonvar(X) ->
 4291            kill(MState),
 4292            (   nonvar(Y) -> abs(X-Y) >= C
 4293            ;   P1 is X - C, P2 is X + C,
 4294                Y in inf..P1 \/ P2..sup
 4295            )
 4296        ;   nonvar(Y) ->
 4297            kill(MState),
 4298            P1 is Y - C, P2 is Y + C,
 4299            X in inf..P1 \/ P2..sup
 4300        ;   true
 4301        ).
 4302
 4303% X #\= Y + Z
 4304run_propagator(x_neq_y_plus_z(X,Y,Z), MState) :-
 4305        (   nonvar(X) ->
 4306            (   nonvar(Y) ->
 4307                (   nonvar(Z) -> kill(MState), X =\= Y + Z
 4308                ;   kill(MState), XY is X - Y, neq_num(Z, XY)
 4309                )
 4310            ;   nonvar(Z) -> kill(MState), XZ is X - Z, neq_num(Y, XZ)
 4311            ;   true
 4312            )
 4313        ;   nonvar(Y) ->
 4314            (   nonvar(Z) ->
 4315                kill(MState), YZ is Y + Z, neq_num(X, YZ)
 4316            ;   Y =:= 0 -> kill(MState), neq(X, Z)
 4317            ;   true
 4318            )
 4319        ;   Z == 0 -> kill(MState), neq(X, Y)
 4320        ;   true
 4321        ).
 4322
 4323% X #=< Y + C
 4324run_propagator(x_leq_y_plus_c(X,Y,C), MState) :-
 4325        (   nonvar(X) ->
 4326            (   nonvar(Y) -> kill(MState), X =< Y + C
 4327            ;   kill(MState),
 4328                R is X - C,
 4329                fd_get(Y, YD, YPs),
 4330                domain_remove_smaller_than(YD, R, YD1),
 4331                fd_put(Y, YD1, YPs)
 4332            )
 4333        ;   nonvar(Y) ->
 4334            kill(MState),
 4335            R is Y + C,
 4336            fd_get(X, XD, XPs),
 4337            domain_remove_greater_than(XD, R, XD1),
 4338            fd_put(X, XD1, XPs)
 4339        ;   (   X == Y -> C >= 0, kill(MState)
 4340            ;   fd_get(Y, YD, _),
 4341                (   domain_supremum(YD, n(YSup)) ->
 4342                    YS1 is YSup + C,
 4343                    fd_get(X, XD, XPs),
 4344                    domain_remove_greater_than(XD, YS1, XD1),
 4345                    fd_put(X, XD1, XPs)
 4346                ;   true
 4347                ),
 4348                (   fd_get(X, XD2, _), domain_infimum(XD2, n(XInf)) ->
 4349                    XI1 is XInf - C,
 4350                    (   fd_get(Y, YD1, YPs1) ->
 4351                        domain_remove_smaller_than(YD1, XI1, YD2),
 4352                        (   domain_infimum(YD2, n(YInf)),
 4353                            domain_supremum(XD2, n(XSup)),
 4354                            XSup =< YInf + C ->
 4355                            kill(MState)
 4356                        ;   true
 4357                        ),
 4358                        fd_put(Y, YD2, YPs1)
 4359                    ;   true
 4360                    )
 4361                ;   true
 4362                )
 4363            )
 4364        ).
 4365
 4366run_propagator(scalar_product_neq(Cs0,Vs0,P0), MState) :-
 4367        coeffs_variables_const(Cs0, Vs0, Cs, Vs, 0, I),
 4368        P is P0 - I,
 4369        (   Vs = [] -> kill(MState), P =\= 0
 4370        ;   Vs = [V], Cs = [C] ->
 4371            kill(MState),
 4372            (   C =:= 1 -> neq_num(V, P)
 4373            ;   C*V #\= P
 4374            )
 4375        ;   Cs == [1,-1] -> kill(MState), Vs = [A,B], x_neq_y_plus_z(A, B, P)
 4376        ;   Cs == [-1,1] -> kill(MState), Vs = [A,B], x_neq_y_plus_z(B, A, P)
 4377        ;   P =:= 0, Cs = [1,1,-1] ->
 4378            kill(MState), Vs = [A,B,C], x_neq_y_plus_z(C, A, B)
 4379        ;   P =:= 0, Cs = [1,-1,1] ->
 4380            kill(MState), Vs = [A,B,C], x_neq_y_plus_z(B, A, C)
 4381        ;   P =:= 0, Cs = [-1,1,1] ->
 4382            kill(MState), Vs = [A,B,C], x_neq_y_plus_z(A, B, C)
 4383        ;   true
 4384        ).
 4385
 4386run_propagator(scalar_product_leq(Cs0,Vs0,P0), MState) :-
 4387        coeffs_variables_const(Cs0, Vs0, Cs, Vs, 0, I),
 4388        P is P0 - I,
 4389        (   Vs = [] -> kill(MState), P >= 0
 4390        ;   sum_finite_domains(Cs, Vs, Infs, Sups, 0, 0, Inf, Sup),
 4391            D1 is P - Inf,
 4392            disable_queue,
 4393            (   Infs == [], Sups == [] ->
 4394                Inf =< P,
 4395                (   Sup =< P -> kill(MState)
 4396                ;   remove_dist_upper_leq(Cs, Vs, D1)
 4397                )
 4398            ;   Infs == [] -> Inf =< P, remove_dist_upper(Sups, D1)
 4399            ;   Sups = [_], Infs = [_] ->
 4400                remove_upper(Infs, D1)
 4401            ;   Infs = [_] -> remove_upper(Infs, D1)
 4402            ;   true
 4403            ),
 4404            enable_queue
 4405        ).
 4406
 4407run_propagator(scalar_product_eq(Cs0,Vs0,P0), MState) :-
 4408        coeffs_variables_const(Cs0, Vs0, Cs, Vs, 0, I),
 4409        P is P0 - I,
 4410        (   Vs = [] -> kill(MState), P =:= 0
 4411        ;   Vs = [V], Cs = [C] -> kill(MState), P mod C =:= 0, V is P // C
 4412        ;   Cs == [1,1] -> kill(MState), Vs = [A,B], A + B #= P
 4413        ;   Cs == [1,-1] -> kill(MState), Vs = [A,B], A #= P + B
 4414        ;   Cs == [-1,1] -> kill(MState), Vs = [A,B], B #= P + A
 4415        ;   Cs == [-1,-1] -> kill(MState), Vs = [A,B], P1 is -P, A + B #= P1
 4416        ;   P =:= 0, Cs == [1,1,-1] -> kill(MState), Vs = [A,B,C], A + B #= C
 4417        ;   P =:= 0, Cs == [1,-1,1] -> kill(MState), Vs = [A,B,C], A + C #= B
 4418        ;   P =:= 0, Cs == [-1,1,1] -> kill(MState), Vs = [A,B,C], B + C #= A
 4419        ;   sum_finite_domains(Cs, Vs, Infs, Sups, 0, 0, Inf, Sup),
 4420            % nl, writeln(Infs-Sups-Inf-Sup),
 4421            D1 is P - Inf,
 4422            D2 is Sup - P,
 4423            disable_queue,
 4424            (   Infs == [], Sups == [] ->
 4425                between(Inf, Sup, P),
 4426                remove_dist_upper_lower(Cs, Vs, D1, D2)
 4427            ;   Sups = [] -> P =< Sup, remove_dist_lower(Infs, D2)
 4428            ;   Infs = [] -> Inf =< P, remove_dist_upper(Sups, D1)
 4429            ;   Sups = [_], Infs = [_] ->
 4430                remove_lower(Sups, D2),
 4431                remove_upper(Infs, D1)
 4432            ;   Infs = [_] -> remove_upper(Infs, D1)
 4433            ;   Sups = [_] -> remove_lower(Sups, D2)
 4434            ;   true
 4435            ),
 4436            enable_queue
 4437        ).
 4438
 4439% X + Y = Z
 4440run_propagator(pplus(X,Y,Z), MState) :-
 4441        (   nonvar(X) ->
 4442            (   X =:= 0 -> kill(MState), Y = Z
 4443            ;   Y == Z -> kill(MState), X =:= 0
 4444            ;   nonvar(Y) -> kill(MState), Z is X + Y
 4445            ;   nonvar(Z) -> kill(MState), Y is Z - X
 4446            ;   fd_get(Z, ZD, ZPs),
 4447                fd_get(Y, YD, _),
 4448                domain_shift(YD, X, Shifted_YD),
 4449                domains_intersection(ZD, Shifted_YD, ZD1),
 4450                fd_put(Z, ZD1, ZPs),
 4451                (   fd_get(Y, YD1, YPs) ->
 4452                    O is -X,
 4453                    domain_shift(ZD1, O, YD2),
 4454                    domains_intersection(YD1, YD2, YD3),
 4455                    fd_put(Y, YD3, YPs)
 4456                ;   true
 4457                )
 4458            )
 4459        ;   nonvar(Y) -> run_propagator(pplus(Y,X,Z), MState)
 4460        ;   nonvar(Z) ->
 4461            (   X == Y -> kill(MState), even(Z), X is Z // 2
 4462            ;   fd_get(X, XD, _),
 4463                fd_get(Y, YD, YPs),
 4464                domain_negate(XD, XDN),
 4465                domain_shift(XDN, Z, YD1),
 4466                domains_intersection(YD, YD1, YD2),
 4467                fd_put(Y, YD2, YPs),
 4468                (   fd_get(X, XD1, XPs) ->
 4469                    domain_negate(YD2, YD2N),
 4470                    domain_shift(YD2N, Z, XD2),
 4471                    domains_intersection(XD1, XD2, XD3),
 4472                    fd_put(X, XD3, XPs)
 4473                ;   true
 4474                )
 4475            )
 4476        ;   (   X == Y -> kill(MState), 2*X #= Z
 4477            ;   X == Z -> kill(MState), Y = 0
 4478            ;   Y == Z -> kill(MState), X = 0
 4479            ;   fd_get(X, XD, XL, XU, XPs),
 4480                fd_get(Y, _, YL, YU, _),
 4481                fd_get(Z, _, ZL, ZU, _),
 4482                NXL cis max(XL, ZL-YU),
 4483                NXU cis min(XU, ZU-YL),
 4484                update_bounds(X, XD, XPs, XL, XU, NXL, NXU),
 4485                (   fd_get(Y, YD2, YL2, YU2, YPs2) ->
 4486                    NYL cis max(YL2, ZL-NXU),
 4487                    NYU cis min(YU2, ZU-NXL),
 4488                    update_bounds(Y, YD2, YPs2, YL2, YU2, NYL, NYU)
 4489                ;   NYL = n(Y), NYU = n(Y)
 4490                ),
 4491                (   fd_get(Z, ZD2, ZL2, ZU2, ZPs2) ->
 4492                    NZL cis max(ZL2,NXL+NYL),
 4493                    NZU cis min(ZU2,NXU+NYU),
 4494                    update_bounds(Z, ZD2, ZPs2, ZL2, ZU2, NZL, NZU)
 4495                ;   true
 4496                )
 4497            )
 4498        ).
 4499
 4500%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 4501
 4502run_propagator(ptimes(X,Y,Z), MState) :-
 4503        (   nonvar(X) ->
 4504            (   nonvar(Y) -> kill(MState), Z is X * Y
 4505            ;   X =:= 0 -> kill(MState), Z = 0
 4506            ;   X =:= 1 -> kill(MState), Z = Y
 4507            ;   nonvar(Z) -> kill(MState), 0 =:= Z mod X, Y is Z // X
 4508            ;   (   Y == Z -> kill(MState), Y = 0
 4509                ;   fd_get(Y, YD, _),
 4510                    fd_get(Z, ZD, ZPs),
 4511                    domain_expand(YD, X, Scaled_YD),
 4512                    domains_intersection(ZD, Scaled_YD, ZD1),
 4513                    fd_put(Z, ZD1, ZPs),
 4514                    (   fd_get(Y, YDom2, YPs2) ->
 4515                        domain_contract(ZD1, X, Contract),
 4516                        domains_intersection(YDom2, Contract, NYDom),
 4517                        fd_put(Y, NYDom, YPs2)
 4518                    ;   kill(MState), Z is X * Y
 4519                    )
 4520                )
 4521            )
 4522        ;   nonvar(Y) -> run_propagator(ptimes(Y,X,Z), MState)
 4523        ;   nonvar(Z) ->
 4524            (   X == Y ->
 4525                kill(MState),
 4526                integer_kth_root(Z, 2, R),
 4527                NR is -R,
 4528                X in NR \/ R
 4529            ;   fd_get(X, XD, XL, XU, XPs),
 4530                fd_get(Y, YD, YL, YU, _),
 4531                min_max_factor(n(Z), n(Z), YL, YU, XL, XU, NXL, NXU),
 4532                update_bounds(X, XD, XPs, XL, XU, NXL, NXU),
 4533                (   fd_get(Y, YD2, YL2, YU2, YPs2) ->
 4534                    min_max_factor(n(Z), n(Z), NXL, NXU, YL2, YU2, NYL, NYU),
 4535                    update_bounds(Y, YD2, YPs2, YL2, YU2, NYL, NYU)
 4536                ;   (   Y =\= 0 -> 0 =:= Z mod Y, kill(MState), X is Z // Y
 4537                    ;   kill(MState), Z = 0
 4538                    )
 4539                ),
 4540                (   Z =:= 0 ->
 4541                    (   \+ domain_contains(XD, 0) -> kill(MState), Y = 0
 4542                    ;   \+ domain_contains(YD, 0) -> kill(MState), X = 0
 4543                    ;   true
 4544                    )
 4545                ;  neq_num(X, 0), neq_num(Y, 0)
 4546                )
 4547            )
 4548        ;   (   X == Y -> kill(MState), X^2 #= Z
 4549            ;   fd_get(X, XD, XL, XU, XPs),
 4550                fd_get(Y, _, YL, YU, _),
 4551                fd_get(Z, ZD, ZL, ZU, _),
 4552                (   Y == Z, \+ domain_contains(ZD, 0) -> kill(MState), X = 1
 4553                ;   X == Z, \+ domain_contains(ZD, 0) -> kill(MState), Y = 1
 4554                ;   min_max_factor(ZL, ZU, YL, YU, XL, XU, NXL, NXU),
 4555                    update_bounds(X, XD, XPs, XL, XU, NXL, NXU),
 4556                    (   fd_get(Y, YD2, YL2, YU2, YPs2) ->
 4557                        min_max_factor(ZL, ZU, NXL, NXU, YL2, YU2, NYL, NYU),
 4558                        update_bounds(Y, YD2, YPs2, YL2, YU2, NYL, NYU)
 4559                    ;   NYL = n(Y), NYU = n(Y)
 4560                    ),
 4561                    (   fd_get(Z, ZD2, ZL2, ZU2, ZPs2) ->
 4562                        min_product(NXL, NXU, NYL, NYU, NZL),
 4563                        max_product(NXL, NXU, NYL, NYU, NZU),
 4564                        (   NZL cis_leq ZL2, NZU cis_geq ZU2 -> ZD3 = ZD2
 4565                        ;   domains_intersection(ZD2, from_to(NZL,NZU), ZD3),
 4566                            fd_put(Z, ZD3, ZPs2)
 4567                        ),
 4568                        (   domain_contains(ZD3, 0) ->  true
 4569                        ;   neq_num(X, 0), neq_num(Y, 0)
 4570                        )
 4571                    ;   true
 4572                    )
 4573                )
 4574            )
 4575        ).
 4576
 4577%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 4578
 4579% X div Y = Z
 4580run_propagator(pdiv(X,Y,Z), MState) :- kill(MState), Z #= (X-(X mod Y)) // Y.
 4581
 4582% X rdiv Y = Z
 4583run_propagator(prdiv(X,Y,Z), MState) :- kill(MState), Z*Y #= X.
 4584
 4585% X // Y = Z (round towards zero)
 4586run_propagator(ptzdiv(X,Y,Z), MState) :-
 4587        (   nonvar(X) ->
 4588            (   nonvar(Y) -> kill(MState), Y =\= 0, Z is X // Y
 4589            ;   fd_get(Y, YD, YL, YU, YPs),
 4590                (   nonvar(Z) ->
 4591                    (   Z =:= 0 ->
 4592                        NYL is -abs(X) - 1,
 4593                        NYU is abs(X) + 1,
 4594                        domains_intersection(YD, split(0, from_to(inf,n(NYL)),
 4595                                                       from_to(n(NYU), sup)),
 4596                                             NYD),
 4597                        fd_put(Y, NYD, YPs)
 4598                    ;   (   sign(X) =:= sign(Z) ->
 4599                            NYL cis max(n(X) // (n(Z)+sign(n(Z))) + n(1), YL),
 4600                            NYU cis min(n(X) // n(Z), YU)
 4601                        ;   NYL cis max(n(X) // n(Z), YL),
 4602                            NYU cis min(n(X) // (n(Z)+sign(n(Z))) - n(1), YU)
 4603                        ),
 4604                        update_bounds(Y, YD, YPs, YL, YU, NYL, NYU)
 4605                    )
 4606                ;   fd_get(Z, ZD, ZL, ZU, ZPs),
 4607                    (   X >= 0, ( YL cis_gt n(0) ; YU cis_lt n(0) )->
 4608                        NZL cis max(n(X)//YU, ZL),
 4609                        NZU cis min(n(X)//YL, ZU)
 4610                    ;   X < 0, ( YL cis_gt n(0) ; YU cis_lt n(0) ) ->
 4611                        NZL cis max(n(X)//YL, ZL),
 4612                        NZU cis min(n(X)//YU, ZU)
 4613                    ;   % TODO: more stringent bounds, cover Y
 4614                        NZL cis max(-abs(n(X)), ZL),
 4615                        NZU cis min(abs(n(X)), ZU)
 4616                    ),
 4617                    update_bounds(Z, ZD, ZPs, ZL, ZU, NZL, NZU),
 4618                    (   X >= 0, NZL cis_gt n(0), fd_get(Y, YD1, YPs1) ->
 4619                        NYL cis n(X) // (NZU + n(1)) + n(1),
 4620                        NYU cis n(X) // NZL,
 4621                        domains_intersection(YD1, from_to(NYL, NYU), NYD1),
 4622                        fd_put(Y, NYD1, YPs1)
 4623                    ;   true
 4624                    )
 4625                )
 4626            )
 4627        ;   nonvar(Y) ->
 4628            Y =\= 0,
 4629            (   Y =:= 1 -> kill(MState), X = Z
 4630            ;   Y =:= -1 -> kill(MState), Z #= -X
 4631            ;   fd_get(X, XD, XL, XU, XPs),
 4632                (   nonvar(Z) ->
 4633                    kill(MState),
 4634                    (   sign(Z) =:= sign(Y) ->
 4635                        NXL cis max(n(Z)*n(Y), XL),
 4636                        NXU cis min((abs(n(Z))+n(1))*abs(n(Y))-n(1), XU)
 4637                    ;   Z =:= 0 ->
 4638                        NXL cis max(-abs(n(Y)) + n(1), XL),
 4639                        NXU cis min(abs(n(Y)) - n(1), XU)
 4640                    ;   NXL cis max((n(Z)+sign(n(Z)))*n(Y)+n(1), XL),
 4641                        NXU cis min(n(Z)*n(Y), XU)
 4642                    ),
 4643                    update_bounds(X, XD, XPs, XL, XU, NXL, NXU)
 4644                ;   fd_get(Z, ZD, ZPs),
 4645                    domain_contract_less(XD, Y, Contracted),
 4646                    domains_intersection(ZD, Contracted, NZD),
 4647                    fd_put(Z, NZD, ZPs),
 4648                    (   fd_get(X, XD2, XPs2) ->
 4649                        domain_expand_more(NZD, Y, Expanded),
 4650                        domains_intersection(XD2, Expanded, NXD2),
 4651                        fd_put(X, NXD2, XPs2)
 4652                    ;   true
 4653                    )
 4654                )
 4655            )
 4656        ;   nonvar(Z) ->
 4657            fd_get(X, XD, XL, XU, XPs),
 4658            fd_get(Y, _, YL, YU, _),
 4659            (   YL cis_geq n(0), XL cis_geq n(0) ->
 4660                NXL cis max(YL*n(Z), XL),
 4661                NXU cis min(YU*(n(Z)+n(1))-n(1), XU)
 4662            ;   %TODO: cover more cases
 4663                NXL = XL, NXU = XU
 4664            ),
 4665            update_bounds(X, XD, XPs, XL, XU, NXL, NXU)
 4666        ;   (   X == Y -> kill(MState), Z = 1
 4667            ;   fd_get(X, _, XL, XU, _),
 4668                fd_get(Y, _, YL, _, _),
 4669                fd_get(Z, ZD, ZPs),
 4670                NZU cis max(abs(XL), XU),
 4671                NZL cis -NZU,
 4672                domains_intersection(ZD, from_to(NZL,NZU), NZD0),
 4673                (   XL cis_geq n(0), YL cis_geq n(0) ->
 4674                    domain_remove_smaller_than(NZD0, 0, NZD1)
 4675                ;   % TODO: cover more cases
 4676                    NZD1 = NZD0
 4677                ),
 4678                fd_put(Z, NZD1, ZPs)
 4679            )
 4680        ).
 4681
 4682
 4683%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 4684% Y = abs(X)
 4685
 4686run_propagator(pabs(X,Y), MState) :-
 4687        (   nonvar(X) -> kill(MState), Y is abs(X)
 4688        ;   nonvar(Y) ->
 4689            kill(MState),
 4690            Y >= 0,
 4691            YN is -Y,
 4692            X in YN \/ Y
 4693        ;   fd_get(X, XD, XPs),
 4694            fd_get(Y, YD, _),
 4695            domain_negate(YD, YDNegative),
 4696            domains_union(YD, YDNegative, XD1),
 4697            domains_intersection(XD, XD1, XD2),
 4698            fd_put(X, XD2, XPs),
 4699            (   fd_get(Y, YD1, YPs1) ->
 4700                domain_negate(XD2, XD2Neg),
 4701                domains_union(XD2, XD2Neg, YD2),
 4702                domain_remove_smaller_than(YD2, 0, YD3),
 4703                domains_intersection(YD1, YD3, YD4),
 4704                fd_put(Y, YD4, YPs1)
 4705            ;   true
 4706            )
 4707        ).
 4708%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 4709% Z = X mod Y
 4710
 4711run_propagator(pmod(X,Y,Z), MState) :-
 4712        (   nonvar(X) ->
 4713            (   nonvar(Y) -> kill(MState), Y =\= 0, Z is X mod Y
 4714            ;   true
 4715            )
 4716        ;   nonvar(Y) ->
 4717            Y =\= 0,
 4718            (   abs(Y) =:= 1 -> kill(MState), Z = 0
 4719            ;   var(Z) ->
 4720                YP is abs(Y) - 1,
 4721                (   Y > 0, fd_get(X, _, n(XL), n(XU), _) ->
 4722                    (   XL >= 0, XU < Y ->
 4723                        kill(MState), Z = X, ZL = XL, ZU = XU
 4724                    ;   ZL = 0, ZU = YP
 4725                    )
 4726                ;   Y > 0 -> ZL = 0, ZU = YP
 4727                ;   YN is -YP, ZL = YN, ZU = 0
 4728                ),
 4729                (   fd_get(Z, ZD, ZPs) ->
 4730                    domains_intersection(ZD, from_to(n(ZL), n(ZU)), ZD1),
 4731                    domain_infimum(ZD1, n(ZMin)),
 4732                    domain_supremum(ZD1, n(ZMax)),
 4733                    fd_put(Z, ZD1, ZPs)
 4734                ;   ZMin = Z, ZMax = Z
 4735                ),
 4736                (   fd_get(X, XD, XPs), domain_infimum(XD, n(XMin)) ->
 4737                    Z1 is XMin mod Y,
 4738                    (   between(ZMin, ZMax, Z1) -> true
 4739                    ;   Y > 0 ->
 4740                        Next is ((XMin - ZMin + Y - 1) div Y)*Y + ZMin,
 4741                        domain_remove_smaller_than(XD, Next, XD1),
 4742                        fd_put(X, XD1, XPs)
 4743                    ;   neq_num(X, XMin)
 4744                    )
 4745                ;   true
 4746                ),
 4747                (   fd_get(X, XD2, XPs2), domain_supremum(XD2, n(XMax)) ->
 4748                    Z2 is XMax mod Y,
 4749                    (   between(ZMin, ZMax, Z2) -> true
 4750                    ;   Y > 0 ->
 4751                        Prev is ((XMax - ZMin) div Y)*Y + ZMax,
 4752                        domain_remove_greater_than(XD2, Prev, XD3),
 4753                        fd_put(X, XD3, XPs2)
 4754                    ;   neq_num(X, XMax)
 4755                    )
 4756                ;   true
 4757                )
 4758            ;   fd_get(X, XD, XPs),
 4759                % if possible, propagate at the boundaries
 4760                (   domain_infimum(XD, n(Min)) ->
 4761                    (   Min mod Y =:= Z -> true
 4762                    ;   Y > 0 ->
 4763                        Next is ((Min - Z + Y - 1) div Y)*Y + Z,
 4764                        domain_remove_smaller_than(XD, Next, XD1),
 4765                        fd_put(X, XD1, XPs)
 4766                    ;   neq_num(X, Min)
 4767                    )
 4768                ;   true
 4769                ),
 4770                (   fd_get(X, XD2, XPs2) ->
 4771                    (   domain_supremum(XD2, n(Max)) ->
 4772                        (   Max mod Y =:= Z -> true
 4773                        ;   Y > 0 ->
 4774                            Prev is ((Max - Z) div Y)*Y + Z,
 4775                            domain_remove_greater_than(XD2, Prev, XD3),
 4776                            fd_put(X, XD3, XPs2)
 4777                        ;   neq_num(X, Max)
 4778                        )
 4779                    ;   true
 4780                    )
 4781                ;   true
 4782                )
 4783            )
 4784        ;   X == Y -> kill(MState), Z = 0
 4785        ;   true % TODO: propagate more
 4786        ).
 4787
 4788%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 4789% Z = X rem Y
 4790
 4791run_propagator(prem(X,Y,Z), MState) :-
 4792        (   nonvar(X) ->
 4793            (   nonvar(Y) -> kill(MState), Y =\= 0, Z is X rem Y
 4794            ;   U is abs(X),
 4795                fd_get(Y, YD, _),
 4796                (   X >=0, domain_infimum(YD, n(Min)), Min >= 0 -> L = 0
 4797                ;   L is -U
 4798                ),
 4799                Z in L..U
 4800            )
 4801        ;   nonvar(Y) ->
 4802            Y =\= 0,
 4803            (   abs(Y) =:= 1 -> kill(MState), Z = 0
 4804            ;   var(Z) ->
 4805                YP is abs(Y) - 1,
 4806                YN is -YP,
 4807                (   Y > 0, fd_get(X, _, n(XL), n(XU), _) ->
 4808                    (   abs(XL) < Y, XU < Y -> kill(MState), Z = X, ZL = XL
 4809                    ;   XL < 0, abs(XL) < Y -> ZL = XL
 4810                    ;   XL >= 0 -> ZL = 0
 4811                    ;   ZL = YN
 4812                    ),
 4813                    (   XU > 0, XU < Y -> ZU = XU
 4814                    ;   XU < 0 -> ZU = 0
 4815                    ;   ZU = YP
 4816                    )
 4817                ;   ZL = YN, ZU = YP
 4818                ),
 4819                (   fd_get(Z, ZD, ZPs) ->
 4820                    domains_intersection(ZD, from_to(n(ZL), n(ZU)), ZD1),
 4821                    fd_put(Z, ZD1, ZPs)
 4822                ;   ZD1 = from_to(n(Z), n(Z))
 4823                ),
 4824                (   fd_get(X, XD, _), domain_infimum(XD, n(Min)) ->
 4825                    Z1 is Min rem Y,
 4826                    (   domain_contains(ZD1, Z1) -> true
 4827                    ;   neq_num(X, Min)
 4828                    )
 4829                ;   true
 4830                ),
 4831                (   fd_get(X, XD1, _), domain_supremum(XD1, n(Max)) ->
 4832                    Z2 is Max rem Y,
 4833                    (   domain_contains(ZD1, Z2) -> true
 4834                    ;   neq_num(X, Max)
 4835                    )
 4836                ;   true
 4837                )
 4838            ;   fd_get(X, XD1, XPs1),
 4839                % if possible, propagate at the boundaries
 4840                (   domain_infimum(XD1, n(Min)) ->
 4841                    (   Min rem Y =:= Z -> true
 4842                    ;   Y > 0, Min > 0 ->
 4843                        Next is ((Min - Z + Y - 1) div Y)*Y + Z,
 4844                        domain_remove_smaller_than(XD1, Next, XD2),
 4845                        fd_put(X, XD2, XPs1)
 4846                    ;   % TODO: bigger steps in other cases as well
 4847                        neq_num(X, Min)
 4848                    )
 4849                ;   true
 4850                ),
 4851                (   fd_get(X, XD3, XPs3) ->
 4852                    (   domain_supremum(XD3, n(Max)) ->
 4853                        (   Max rem Y =:= Z -> true
 4854                        ;   Y > 0, Max > 0  ->
 4855                            Prev is ((Max - Z) div Y)*Y + Z,
 4856                            domain_remove_greater_than(XD3, Prev, XD4),
 4857                            fd_put(X, XD4, XPs3)
 4858                        ;   % TODO: bigger steps in other cases as well
 4859                            neq_num(X, Max)
 4860                        )
 4861                    ;   true
 4862                    )
 4863                ;   true
 4864                )
 4865            )
 4866        ;   X == Y -> kill(MState), Z = 0
 4867        ;   fd_get(Z, ZD, ZPs) ->
 4868            fd_get(Y, _, YInf, YSup, _),
 4869            fd_get(X, _, XInf, XSup, _),
 4870            M cis max(abs(YInf),YSup),
 4871            (   XInf cis_geq n(0) -> Inf0 = n(0)
 4872            ;   Inf0 = XInf
 4873            ),
 4874            (   XSup cis_leq n(0) -> Sup0 = n(0)
 4875            ;   Sup0 = XSup
 4876            ),
 4877            NInf cis max(max(Inf0, -M + n(1)), min(XInf,-XSup)),
 4878            NSup cis min(min(Sup0, M - n(1)), max(abs(XInf),XSup)),
 4879            domains_intersection(ZD, from_to(NInf,NSup), ZD1),
 4880            fd_put(Z, ZD1, ZPs)
 4881        ;   true % TODO: propagate more
 4882        ).
 4883
 4884%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 4885% Z = max(X,Y)
 4886
 4887run_propagator(pmax(X,Y,Z), MState) :-
 4888        (   nonvar(X) ->
 4889            (   nonvar(Y) -> kill(MState), Z is max(X,Y)
 4890            ;   nonvar(Z) ->
 4891                (   Z =:= X -> kill(MState), X #>= Y
 4892                ;   Z > X -> Z = Y
 4893                ;   false % Z < X
 4894                )
 4895            ;   fd_get(Y, _, YInf, YSup, _),
 4896                (   YInf cis_gt n(X) -> Z = Y
 4897                ;   YSup cis_lt n(X) -> Z = X
 4898                ;   YSup = n(M) ->
 4899                    fd_get(Z, ZD, ZPs),
 4900                    domain_remove_greater_than(ZD, M, ZD1),
 4901                    fd_put(Z, ZD1, ZPs)
 4902                ;   true
 4903                )
 4904            )
 4905        ;   nonvar(Y) -> run_propagator(pmax(Y,X,Z), MState)
 4906        ;   fd_get(Z, ZD, ZPs) ->
 4907            fd_get(X, _, XInf, XSup, _),
 4908            fd_get(Y, _, YInf, YSup, _),
 4909            (   YInf cis_gt YSup -> kill(MState), Z = Y
 4910            ;   YSup cis_lt XInf -> kill(MState), Z = X
 4911            ;   n(M) cis max(XSup, YSup) ->
 4912                domain_remove_greater_than(ZD, M, ZD1),
 4913                fd_put(Z, ZD1, ZPs)
 4914            ;   true
 4915            )
 4916        ;   true
 4917        ).
 4918
 4919%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 4920% Z = min(X,Y)
 4921
 4922run_propagator(pmin(X,Y,Z), MState) :-
 4923        (   nonvar(X) ->
 4924            (   nonvar(Y) -> kill(MState), Z is min(X,Y)
 4925            ;   nonvar(Z) ->
 4926                (   Z =:= X -> kill(MState), X #=< Y
 4927                ;   Z < X -> Z = Y
 4928                ;   false % Z > X
 4929                )
 4930            ;   fd_get(Y, _, YInf, YSup, _),
 4931                (   YSup cis_lt n(X) -> Z = Y
 4932                ;   YInf cis_gt n(X) -> Z = X
 4933                ;   YInf = n(M) ->
 4934                    fd_get(Z, ZD, ZPs),
 4935                    domain_remove_smaller_than(ZD, M, ZD1),
 4936                    fd_put(Z, ZD1, ZPs)
 4937                ;   true
 4938                )
 4939            )
 4940        ;   nonvar(Y) -> run_propagator(pmin(Y,X,Z), MState)
 4941        ;   fd_get(Z, ZD, ZPs) ->
 4942            fd_get(X, _, XInf, XSup, _),
 4943            fd_get(Y, _, YInf, YSup, _),
 4944            (   YSup cis_lt YInf -> kill(MState), Z = Y
 4945            ;   YInf cis_gt XSup -> kill(MState), Z = X
 4946            ;   n(M) cis min(XInf, YInf) ->
 4947                domain_remove_smaller_than(ZD, M, ZD1),
 4948                fd_put(Z, ZD1, ZPs)
 4949            ;   true
 4950            )
 4951        ;   true
 4952        ).
 4953%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 4954% Z = X ^ Y
 4955
 4956run_propagator(pexp(X,Y,Z), MState) :-
 4957        (   X == 1 -> kill(MState), Z = 1
 4958        ;   X == 0 -> kill(MState), Z in 0..1, Z #<==> Y #= 0
 4959        ;   Y == 0 -> kill(MState), Z = 1
 4960        ;   Y == 1 -> kill(MState), Z = X
 4961        ;   nonvar(X) ->
 4962            (   nonvar(Y) ->
 4963                (   Y >= 0 -> true ; X =:= -1 ),
 4964                kill(MState),
 4965                Z is X^Y
 4966            ;   nonvar(Z) ->
 4967                (   Z > 1 ->
 4968                    kill(MState),
 4969                    integer_log_b(Z, X, 1, Y)
 4970                ;   true
 4971                )
 4972            ;   fd_get(Y, _, YL, YU, _),
 4973                fd_get(Z, ZD, ZPs),
 4974                (   X > 0, YL cis_geq n(0) ->
 4975                    NZL cis n(X)^YL,
 4976                    NZU cis n(X)^YU,
 4977                    domains_intersection(ZD, from_to(NZL,NZU), NZD),
 4978                    fd_put(Z, NZD, ZPs)
 4979                ;   true
 4980                ),
 4981                (   X > 0,
 4982                    fd_get(Z, _, _, n(ZMax), _),
 4983                    ZMax > 0 ->
 4984                    floor_integer_log_b(ZMax, X, 1, YCeil),
 4985                    Y in inf..YCeil
 4986                ;   true
 4987                )
 4988            )
 4989        ;   nonvar(Z) ->
 4990            (   nonvar(Y) ->
 4991                integer_kth_root(Z, Y, R),
 4992                kill(MState),
 4993                (   even(Y) ->
 4994                    N is -R,
 4995                    X in N \/ R
 4996                ;   X = R
 4997                )
 4998            ;   fd_get(X, _, n(NXL), _, _), NXL > 1 ->
 4999                (   Z > 1, between(NXL, Z, Exp), NXL^Exp > Z ->
 5000                    Exp1 is Exp - 1,
 5001                    fd_get(Y, YD, YPs),
 5002                    domains_intersection(YD, from_to(n(1),n(Exp1)), YD1),
 5003                    fd_put(Y, YD1, YPs),
 5004                    (   fd_get(X, XD, XPs) ->
 5005                        domain_infimum(YD1, n(YL)),
 5006                        integer_kth_root_leq(Z, YL, RU),
 5007                        domains_intersection(XD, from_to(n(NXL),n(RU)), XD1),
 5008                        fd_put(X, XD1, XPs)
 5009                    ;   true
 5010                    )
 5011                ;   true
 5012                )
 5013            ;   true
 5014            )
 5015        ;   nonvar(Y), Y > 0 ->
 5016            (   even(Y) ->
 5017                geq(Z, 0)
 5018            ;   true
 5019            ),
 5020            (   fd_get(X, XD, XL, XU, _), fd_get(Z, ZD, ZL, ZU, ZPs) ->
 5021                (   domain_contains(ZD, 0) -> XD1 = XD
 5022                ;   domain_remove(XD, 0, XD1)
 5023                ),
 5024                (   domain_contains(XD, 0) -> ZD1 = ZD
 5025                ;   domain_remove(ZD, 0, ZD1)
 5026                ),
 5027                (   even(Y) ->
 5028                    (   XL cis_geq n(0) ->
 5029                        NZL cis XL^n(Y)
 5030                    ;   XU cis_leq n(0) ->
 5031                        NZL cis XU^n(Y)
 5032                    ;   NZL = n(0)
 5033                    ),
 5034                    NZU cis max(abs(XL),abs(XU))^n(Y),
 5035                    domains_intersection(ZD1, from_to(NZL,NZU), ZD2)
 5036                ;   (   finite(XL) ->
 5037                        NZL cis XL^n(Y),
 5038                        NZU cis XU^n(Y),
 5039                        domains_intersection(ZD1, from_to(NZL,NZU), ZD2)
 5040                    ;   ZD2 = ZD1
 5041                    )
 5042                ),
 5043                fd_put(Z, ZD2, ZPs),
 5044                (   even(Y), ZU = n(Num) ->
 5045                    integer_kth_root_leq(Num, Y, RU),
 5046                    (   XL cis_geq n(0), ZL = n(Num1) ->
 5047                        integer_kth_root_leq(Num1, Y, RL0),
 5048                        (   RL0^Y < Num1 -> RL is RL0 + 1
 5049                        ;   RL = RL0
 5050                        )
 5051                    ;   RL is -RU
 5052                    ),
 5053                    RL =< RU,
 5054                    NXD = from_to(n(RL),n(RU))
 5055                ;   odd(Y), ZL cis_geq n(0), ZU = n(Num) ->
 5056                    integer_kth_root_leq(Num, Y, RU),
 5057                    ZL = n(Num1),
 5058                    integer_kth_root_leq(Num1, Y, RL0),
 5059                    (   RL0^Y < Num1 -> RL is RL0 + 1
 5060                    ;   RL = RL0
 5061                    ),
 5062                    RL =< RU,
 5063                    NXD = from_to(n(RL),n(RU))
 5064                ;   NXD = XD1   % TODO: propagate more
 5065                ),
 5066                (   fd_get(X, XD2, XPs) ->
 5067                    domains_intersection(XD2, XD1, XD3),
 5068                    domains_intersection(XD3, NXD, XD4),
 5069                    fd_put(X, XD4, XPs)
 5070                ;   true
 5071                )
 5072            ;   true
 5073            )
 5074        ;   fd_get(X, _, XL, _, _),
 5075            XL cis_gt n(0),
 5076            fd_get(Y, _, YL, _, _),
 5077            YL cis_gt n(0),
 5078            fd_get(Z, ZD, ZPs) ->
 5079            n(NZL) cis XL^YL,
 5080            domain_remove_smaller_than(ZD, NZL, ZD1),
 5081            fd_put(Z, ZD1, ZPs)
 5082        ;   true
 5083        ).
 5084
 5085%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 5086run_propagator(pzcompare(Order, A, B), MState) :-
 5087        (   A == B -> kill(MState), Order = (=)
 5088        ;   (   nonvar(A) ->
 5089                (   nonvar(B) ->
 5090                    kill(MState),
 5091                    (   A > B -> Order = (>)
 5092                    ;   Order = (<)
 5093                    )
 5094                ;   fd_get(B, _, BL, BU, _),
 5095                    (   BL cis_gt n(A) -> kill(MState), Order = (<)
 5096                    ;   BU cis_lt n(A) -> kill(MState), Order = (>)
 5097                    ;   true
 5098                    )
 5099                )
 5100            ;   nonvar(B) ->
 5101                fd_get(A, _, AL, AU, _),
 5102                (   AL cis_gt n(B) -> kill(MState), Order = (>)
 5103                ;   AU cis_lt n(B) -> kill(MState), Order = (<)
 5104                ;   true
 5105                )
 5106            ;   fd_get(A, _, AL, AU, _),
 5107                fd_get(B, _, BL, BU, _),
 5108                (   AL cis_gt BU -> kill(MState), Order = (>)
 5109                ;   AU cis_lt BL -> kill(MState), Order = (<)
 5110                ;   true
 5111                )
 5112            )
 5113        ).
 5114
 5115%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 5116
 5117% reified constraints
 5118
 5119%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 5120
 5121run_propagator(reified_in(V,Dom,B), MState) :-
 5122        (   integer(V) ->
 5123            kill(MState),
 5124            (   domain_contains(Dom, V) -> B = 1
 5125            ;   B = 0
 5126            )
 5127        ;   B == 1 -> kill(MState), domain(V, Dom)
 5128        ;   B == 0 -> kill(MState), domain_complement(Dom, C), domain(V, C)
 5129        ;   fd_get(V, VD, _),
 5130            (   domains_intersection(VD, Dom, I) ->
 5131                (   I == VD -> kill(MState), B = 1
 5132                ;   true
 5133                )
 5134            ;   kill(MState), B = 0
 5135            )
 5136        ).
 5137
 5138run_propagator(reified_tuple_in(Tuple, R, B), MState) :-
 5139        get_attr(R, clpfd_relation, Relation),
 5140        (   B == 1 -> kill(MState), tuples_in([Tuple], Relation)
 5141        ;   (   ground(Tuple) ->
 5142                kill(MState),
 5143                (   memberchk(Tuple, Relation) -> B = 1
 5144                ;   B = 0
 5145                )
 5146            ;   relation_unifiable(Relation, Tuple, Us, _, _),
 5147                (   Us = [] -> kill(MState), B = 0
 5148                ;   true
 5149                )
 5150            )
 5151        ).
 5152
 5153run_propagator(tuples_not_in(Tuples, Relation, B), MState) :-
 5154        (   B == 0 ->
 5155            kill(MState),
 5156            tuples_in_conjunction(Tuples, Relation, Conj),
 5157            #\ Conj
 5158        ;   true
 5159        ).
 5160
 5161run_propagator(kill_reified_tuples(B, Ps, Bs), _) :-
 5162        (   B == 0 ->
 5163            maplist(kill_entailed, Ps),
 5164            phrase(as(Bs), As),
 5165            maplist(kill_entailed, As)
 5166        ;   true
 5167        ).
 5168
 5169run_propagator(reified_fd(V,B), MState) :-
 5170        (   fd_inf(V, I), I \== inf, fd_sup(V, S), S \== sup ->
 5171            kill(MState),
 5172            B = 1
 5173        ;   B == 0 ->
 5174            (   fd_inf(V, inf) -> true
 5175            ;   fd_sup(V, sup) -> true
 5176            ;   false
 5177            )
 5178        ;   true
 5179        ).
 5180
 5181% The result of X/Y, X mod Y, and X rem Y is undefined iff Y is 0.
 5182
 5183run_propagator(pskeleton(X,Y,D,Skel,Z,_), MState) :-
 5184        (   Y == 0 -> kill(MState), D = 0
 5185        ;   D == 1 -> kill(MState), neq_num(Y, 0), skeleton([X,Y,Z], Skel)
 5186        ;   integer(Y), Y =\= 0 -> kill(MState), D = 1, skeleton([X,Y,Z], Skel)
 5187        ;   fd_get(Y, YD, _), \+ domain_contains(YD, 0) ->
 5188            kill(MState),
 5189            D = 1, skeleton([X,Y,Z], Skel)
 5190        ;   true
 5191        ).
 5192
 5193/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 5194   Propagators for arithmetic functions that only propagate
 5195   functionally. These are currently the bitwise operations.
 5196- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 5197
 5198run_propagator(pfunction(Op,A,B,R), MState) :-
 5199        (   integer(A), integer(B) ->
 5200            kill(MState),
 5201            Expr =.. [Op,A,B],
 5202            R is Expr
 5203        ;   true
 5204        ).
 5205run_propagator(pfunction(Op,A,R), MState) :-
 5206        (   integer(A) ->
 5207            kill(MState),
 5208            Expr =.. [Op,A],
 5209            R is Expr
 5210        ;   true
 5211        ).
 5212
 5213%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 5214
 5215run_propagator(reified_geq(DX,X,DY,Y,Ps,B), MState) :-
 5216        (   DX == 0 -> kill(MState, Ps), B = 0
 5217        ;   DY == 0 -> kill(MState, Ps), B = 0
 5218        ;   B == 1 -> kill(MState), DX = 1, DY = 1, geq(X, Y)
 5219        ;   DX == 1, DY == 1 ->
 5220            (   var(B) ->
 5221                (   nonvar(X) ->
 5222                    (   nonvar(Y) ->
 5223                        kill(MState),
 5224                        (   X >= Y -> B = 1 ; B = 0 )
 5225                    ;   fd_get(Y, _, YL, YU, _),
 5226                        (   n(X) cis_geq YU -> kill(MState, Ps), B = 1
 5227                        ;   n(X) cis_lt YL -> kill(MState, Ps), B = 0
 5228                        ;   true
 5229                        )
 5230                    )
 5231                ;   nonvar(Y) ->
 5232                    fd_get(X, _, XL, XU, _),
 5233                    (   XL cis_geq n(Y) -> kill(MState, Ps), B = 1
 5234                    ;   XU cis_lt n(Y) -> kill(MState, Ps), B = 0
 5235                    ;   true
 5236                    )
 5237                ;   X == Y -> kill(MState, Ps), B = 1
 5238                ;   fd_get(X, _, XL, XU, _),
 5239                    fd_get(Y, _, YL, YU, _),
 5240                    (   XL cis_geq YU -> kill(MState, Ps), B = 1
 5241                    ;   XU cis_lt YL -> kill(MState, Ps), B = 0
 5242                    ;   true
 5243                    )
 5244                )
 5245            ;   B =:= 0 -> kill(MState), X #< Y
 5246            ;   true
 5247            )
 5248        ;   true
 5249</