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

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

 1743labeling(Options, Vars) :-
 1744        must_be(list, Options),
 1745        fd_must_be_list(Vars),
 1746        maplist(must_be_finite_fdvar, Vars),
 1747        label(Options, Options, default(leftmost), default(up), default(step), [], upto_ground, Vars).
 1748
 1749finite_domain(Dom) :-
 1750        domain_infimum(Dom, n(_)),
 1751        domain_supremum(Dom, n(_)).
 1752
 1753must_be_finite_fdvar(Var) :-
 1754        (   fd_get(Var, Dom, _) ->
 1755            (   finite_domain(Dom) -> true
 1756            ;   instantiation_error(Var)
 1757            )
 1758        ;   integer(Var) -> true
 1759        ;   must_be(integer, Var)
 1760        ).
 1761
 1762
 1763label([O|Os], Options, Selection, Order, Choice, Optim, Consistency, Vars) :-
 1764        (   var(O)-> instantiation_error(O)
 1765        ;   override(selection, Selection, O, Options, S1) ->
 1766            label(Os, Options, S1, Order, Choice, Optim, Consistency, Vars)
 1767        ;   override(order, Order, O, Options, O1) ->
 1768            label(Os, Options, Selection, O1, Choice, Optim, Consistency, Vars)
 1769        ;   override(choice, Choice, O, Options, C1) ->
 1770            label(Os, Options, Selection, Order, C1, Optim, Consistency, Vars)
 1771        ;   optimisation(O) ->
 1772            label(Os, Options, Selection, Order, Choice, [O|Optim], Consistency, Vars)
 1773        ;   consistency(O, O1) ->
 1774            label(Os, Options, Selection, Order, Choice, Optim, O1, Vars)
 1775        ;   domain_error(labeling_option, O)
 1776        ).
 1777label([], _, Selection, Order, Choice, Optim0, Consistency, Vars) :-
 1778        maplist(arg(1), [Selection,Order,Choice], [S,O,C]),
 1779        (   Optim0 == [] ->
 1780            label(Vars, S, O, C, Consistency)
 1781        ;   reverse(Optim0, Optim),
 1782            exprs_singlevars(Optim, SVs),
 1783            optimise(Vars, [S,O,C], SVs)
 1784        ).
 1785
 1786% Introduce new variables for each min/max expression to avoid
 1787% reparsing expressions during optimisation.
 1788
 1789exprs_singlevars([], []).
 1790exprs_singlevars([E|Es], [SV|SVs]) :-
 1791        E =.. [F,Expr],
 1792        ?(Single) #= Expr,
 1793        SV =.. [F,Single],
 1794        exprs_singlevars(Es, SVs).
 1795
 1796all_dead(fd_props(Bs,Gs,Os)) :-
 1797        all_dead_(Bs),
 1798        all_dead_(Gs),
 1799        all_dead_(Os).
 1800
 1801all_dead_([]).
 1802all_dead_([propagator(_, S)|Ps]) :- S == dead, all_dead_(Ps).
 1803
 1804label([], _, _, _, Consistency) :- !,
 1805        (   Consistency = upto_in(I0,I) -> I0 = I
 1806        ;   true
 1807        ).
 1808label(Vars, Selection, Order, Choice, Consistency) :-
 1809        (   Vars = [V|Vs], nonvar(V) -> label(Vs, Selection, Order, Choice, Consistency)
 1810        ;   select_var(Selection, Vars, Var, RVars),
 1811            (   var(Var) ->
 1812                (   Consistency = upto_in(I0,I), fd_get(Var, _, Ps), all_dead(Ps) ->
 1813                    fd_size(Var, Size),
 1814                    I1 is I0*Size,
 1815                    label(RVars, Selection, Order, Choice, upto_in(I1,I))
 1816                ;   Consistency = upto_in, fd_get(Var, _, Ps), all_dead(Ps) ->
 1817                    label(RVars, Selection, Order, Choice, Consistency)
 1818                ;   choice_order_variable(Choice, Order, Var, RVars, Vars, Selection, Consistency)
 1819                )
 1820            ;   label(RVars, Selection, Order, Choice, Consistency)
 1821            )
 1822        ).
 1823
 1824choice_order_variable(step, Order, Var, Vars, Vars0, Selection, Consistency) :-
 1825        fd_get(Var, Dom, _),
 1826        order_dom_next(Order, Dom, Next),
 1827        (   Var = Next,
 1828            label(Vars, Selection, Order, step, Consistency)
 1829        ;   neq_num(Var, Next),
 1830            do_queue,
 1831            label(Vars0, Selection, Order, step, Consistency)
 1832        ).
 1833choice_order_variable(enum, Order, Var, Vars, _, Selection, Consistency) :-
 1834        fd_get(Var, Dom0, _),
 1835        domain_direction_element(Dom0, Order, Var),
 1836        label(Vars, Selection, Order, enum, Consistency).
 1837choice_order_variable(bisect, Order, Var, _, Vars0, Selection, Consistency) :-
 1838        fd_get(Var, Dom, _),
 1839        domain_infimum(Dom, n(I)),
 1840        domain_supremum(Dom, n(S)),
 1841        Mid0 is (I + S) // 2,
 1842        (   Mid0 =:= S -> Mid is Mid0 - 1 ; Mid = Mid0 ),
 1843        (   Order == up -> ( Var #=< Mid ; Var #> Mid )
 1844        ;   Order == down -> ( Var #> Mid ; Var #=< Mid )
 1845        ;   domain_error(bisect_up_or_down, Order)
 1846        ),
 1847        label(Vars0, Selection, Order, bisect, Consistency).
 1848
 1849override(What, Prev, Value, Options, Result) :-
 1850        call(What, Value),
 1851        override_(Prev, Value, Options, Result).
 1852
 1853override_(default(_), Value, _, user(Value)).
 1854override_(user(Prev), Value, Options, _) :-
 1855        (   Value == Prev ->
 1856            domain_error(nonrepeating_labeling_options, Options)
 1857        ;   domain_error(consistent_labeling_options, Options)
 1858        ).
 1859
 1860selection(ff).
 1861selection(ffc).
 1862selection(min).
 1863selection(max).
 1864selection(leftmost).
 1865selection(random_variable(Seed)) :-
 1866        must_be(integer, Seed),
 1867        set_random(seed(Seed)).
 1868
 1869choice(step).
 1870choice(enum).
 1871choice(bisect).
 1872
 1873order(up).
 1874order(down).
 1875% TODO: random_variable and random_value currently both set the seed,
 1876% so exchanging the options can yield different results.
 1877order(random_value(Seed)) :-
 1878        must_be(integer, Seed),
 1879        set_random(seed(Seed)).
 1880
 1881consistency(upto_in(I), upto_in(1, I)).
 1882consistency(upto_in, upto_in).
 1883consistency(upto_ground, upto_ground).
 1884
 1885optimisation(min(_)).
 1886optimisation(max(_)).
 1887
 1888select_var(leftmost, [Var|Vars], Var, Vars).
 1889select_var(min, [V|Vs], Var, RVars) :-
 1890        find_min(Vs, V, Var),
 1891        delete_eq([V|Vs], Var, RVars).
 1892select_var(max, [V|Vs], Var, RVars) :-
 1893        find_max(Vs, V, Var),
 1894        delete_eq([V|Vs], Var, RVars).
 1895select_var(ff, [V|Vs], Var, RVars) :-
 1896        fd_size_(V, n(S)),
 1897        find_ff(Vs, V, S, Var),
 1898        delete_eq([V|Vs], Var, RVars).
 1899select_var(ffc, [V|Vs], Var, RVars) :-
 1900        find_ffc(Vs, V, Var),
 1901        delete_eq([V|Vs], Var, RVars).
 1902select_var(random_variable(_), Vars0, Var, Vars) :-
 1903        length(Vars0, L),
 1904        I is random(L),
 1905        nth0(I, Vars0, Var),
 1906        delete_eq(Vars0, Var, Vars).
 1907
 1908find_min([], Var, Var).
 1909find_min([V|Vs], CM, Min) :-
 1910        (   min_lt(V, CM) ->
 1911            find_min(Vs, V, Min)
 1912        ;   find_min(Vs, CM, Min)
 1913        ).
 1914
 1915find_max([], Var, Var).
 1916find_max([V|Vs], CM, Max) :-
 1917        (   max_gt(V, CM) ->
 1918            find_max(Vs, V, Max)
 1919        ;   find_max(Vs, CM, Max)
 1920        ).
 1921
 1922find_ff([], Var, _, Var).
 1923find_ff([V|Vs], CM, S0, FF) :-
 1924        (   nonvar(V) -> find_ff(Vs, CM, S0, FF)
 1925        ;   (   fd_size_(V, n(S1)), S1 < S0 ->
 1926                find_ff(Vs, V, S1, FF)
 1927            ;   find_ff(Vs, CM, S0, FF)
 1928            )
 1929        ).
 1930
 1931find_ffc([], Var, Var).
 1932find_ffc([V|Vs], Prev, FFC) :-
 1933        (   ffc_lt(V, Prev) ->
 1934            find_ffc(Vs, V, FFC)
 1935        ;   find_ffc(Vs, Prev, FFC)
 1936        ).
 1937
 1938
 1939ffc_lt(X, Y) :-
 1940        (   fd_get(X, XD, XPs) ->
 1941            domain_num_elements(XD, n(NXD))
 1942        ;   NXD = 1, XPs = []
 1943        ),
 1944        (   fd_get(Y, YD, YPs) ->
 1945            domain_num_elements(YD, n(NYD))
 1946        ;   NYD = 1, YPs = []
 1947        ),
 1948        (   NXD < NYD -> true
 1949        ;   NXD =:= NYD,
 1950            props_number(XPs, NXPs),
 1951            props_number(YPs, NYPs),
 1952            NXPs > NYPs
 1953        ).
 1954
 1955min_lt(X,Y) :- bounds(X,LX,_), bounds(Y,LY,_), LX < LY.
 1956
 1957max_gt(X,Y) :- bounds(X,_,UX), bounds(Y,_,UY), UX > UY.
 1958
 1959bounds(X, L, U) :-
 1960        (   fd_get(X, Dom, _) ->
 1961            domain_infimum(Dom, n(L)),
 1962            domain_supremum(Dom, n(U))
 1963        ;   L = X, U = L
 1964        ).
 1965
 1966delete_eq([], _, []).
 1967delete_eq([X|Xs], Y, List) :-
 1968        (   nonvar(X) -> delete_eq(Xs, Y, List)
 1969        ;   X == Y -> List = Xs
 1970        ;   List = [X|Tail],
 1971            delete_eq(Xs, Y, Tail)
 1972        ).
 1973
 1974/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1975   contracting/1 -- subject to change
 1976
 1977   This can remove additional domain elements from the boundaries.
 1978- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1979
 1980contracting(Vs) :-
 1981        must_be(list, Vs),
 1982        maplist(must_be_finite_fdvar, Vs),
 1983        contracting(Vs, false, Vs).
 1984
 1985contracting([], Repeat, Vars) :-
 1986        (   Repeat -> contracting(Vars, false, Vars)
 1987        ;   true
 1988        ).
 1989contracting([V|Vs], Repeat, Vars) :-
 1990        fd_inf(V, Min),
 1991        (   \+ \+ (V = Min) ->
 1992            fd_sup(V, Max),
 1993            (   \+ \+ (V = Max) ->
 1994                contracting(Vs, Repeat, Vars)
 1995            ;   V #\= Max,
 1996                contracting(Vs, true, Vars)
 1997            )
 1998        ;   V #\= Min,
 1999            contracting(Vs, true, Vars)
 2000        ).
 2001
 2002/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 2003   fds_sespsize(Vs, S).
 2004
 2005   S is an upper bound on the search space size with respect to finite
 2006   domain variables Vs.
 2007- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 2008
 2009fds_sespsize(Vs, S) :-
 2010        must_be(list, Vs),
 2011        maplist(fd_variable, Vs),
 2012        fds_sespsize(Vs, n(1), S1),
 2013        bound_portray(S1, S).
 2014
 2015fd_size_(V, S) :-
 2016        (   fd_get(V, D, _) ->
 2017            domain_num_elements(D, S)
 2018        ;   S = n(1)
 2019        ).
 2020
 2021fds_sespsize([], S, S).
 2022fds_sespsize([V|Vs], S0, S) :-
 2023        fd_size_(V, S1),
 2024        S2 cis S0*S1,
 2025        fds_sespsize(Vs, S2, S).
 2026
 2027/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 2028   Optimisation uses destructive assignment to save the computed
 2029   extremum over backtracking. Failure is used to get rid of copies of
 2030   attributed variables that are created in intermediate steps. At
 2031   least that's the intention - it currently doesn't work in SWI:
 2032
 2033   %?- X in 0..3, call_residue_vars(labeling([min(X)], [X]), Vs).
 2034   %@ X = 0,
 2035   %@ Vs = [_G6174, _G6177],
 2036   %@ _G6174 in 0..3
 2037
 2038- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 2039
 2040optimise(Vars, Options, Whats) :-
 2041        Whats = [What|WhatsRest],
 2042        Extremum = extremum(none),
 2043        (   catch(store_extremum(Vars, Options, What, Extremum),
 2044                  time_limit_exceeded,
 2045                  false)
 2046        ;   Extremum = extremum(n(Val)),
 2047            arg(1, What, Expr),
 2048            append(WhatsRest, Options, Options1),
 2049            (   Expr #= Val,
 2050                labeling(Options1, Vars)
 2051            ;   Expr #\= Val,
 2052                optimise(Vars, Options, Whats)
 2053            )
 2054        ).
 2055
 2056store_extremum(Vars, Options, What, Extremum) :-
 2057        catch((labeling(Options, Vars), throw(w(What))), w(What1), true),
 2058        functor(What, Direction, _),
 2059        maplist(arg(1), [What,What1], [Expr,Expr1]),
 2060        optimise(Direction, Options, Vars, Expr1, Expr, Extremum).
 2061
 2062optimise(Direction, Options, Vars, Expr0, Expr, Extremum) :-
 2063        must_be(ground, Expr0),
 2064        nb_setarg(1, Extremum, n(Expr0)),
 2065        catch((tighten(Direction, Expr, Expr0),
 2066               labeling(Options, Vars),
 2067               throw(v(Expr))), v(Expr1), true),
 2068        optimise(Direction, Options, Vars, Expr1, Expr, Extremum).
 2069
 2070tighten(min, E, V) :- E #< V.
 2071tighten(max, E, V) :- E #> V.
 2072
 2073%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 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.
 2081all_different(Ls) :-
 2082        fd_must_be_list(Ls),
 2083        maplist(fd_variable, Ls),
 2084        Orig = original_goal(_, all_different(Ls)),
 2085        all_different(Ls, [], Orig),
 2086        do_queue.
 2087
 2088all_different([], _, _).
 2089all_different([X|Right], Left, Orig) :-
 2090        (   var(X) ->
 2091            make_propagator(pdifferent(Left,Right,X,Orig), Prop),
 2092            init_propagator(X, Prop),
 2093            trigger_prop(Prop)
 2094        ;   exclude_fire(Left, Right, X)
 2095        ),
 2096        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.
 2111all_distinct(Ls) :-
 2112        fd_must_be_list(Ls),
 2113        maplist(fd_variable, Ls),
 2114        make_propagator(pdistinct(Ls), Prop),
 2115        distinct_attach(Ls, Prop, []),
 2116        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.
 2131sum(Vs, Op, Value) :-
 2132        must_be(list, Vs),
 2133        same_length(Vs, Ones),
 2134        maplist(=(1), Ones),
 2135        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 #>=.
 2143scalar_product(Cs, Vs, Op, Value) :-
 2144        must_be(list(integer), Cs),
 2145        must_be(list, Vs),
 2146        maplist(fd_variable, Vs),
 2147        (   Op = (#=), single_value(Value, Right), ground(Vs) ->
 2148            foldl(coeff_int_linsum, Cs, Vs, 0, Right)
 2149        ;   must_be(callable, Op),
 2150            (   memberchk(Op, [#=,#\=,#<,#>,#=<,#>=]) -> true
 2151            ;   domain_error(scalar_product_relation, Op)
 2152            ),
 2153            must_be(acyclic, Value),
 2154            foldl(coeff_var_plusterm, Cs, Vs, 0, Left),
 2155            (   left_right_linsum_const(Left, Value, Cs1, Vs1, Const) ->
 2156                scalar_product_(Op, Cs1, Vs1, Const)
 2157            ;   sum(Cs, Vs, 0, Op, Value)
 2158            )
 2159        ).
 2160
 2161single_value(V, V)    :- var(V), !, non_monotonic(V).
 2162single_value(V, V)    :- integer(V).
 2163single_value(?(V), V) :- fd_variable(V).
 2164
 2165coeff_var_plusterm(C, V, T0, T0+(C* ?(V))).
 2166
 2167coeff_int_linsum(C, I, S0, S) :- S is S0 + C*I.
 2168
 2169sum([], _, Sum, Op, Value) :- call(Op, Sum, Value).
 2170sum([C|Cs], [X|Xs], Acc, Op, Value) :-
 2171        ?(NAcc) #= Acc + C* ?(X),
 2172        sum(Cs, Xs, NAcc, Op, Value).
 2173
 2174multiples([], [], _).
 2175multiples([C|Cs], [V|Vs], Left) :-
 2176        (   (   Cs = [N|_] ; Left = [N|_] ) ->
 2177            (   N =\= 1, gcd(C,N) =:= 1 ->
 2178                gcd(Cs, N, GCD0),
 2179                gcd(Left, GCD0, GCD),
 2180                (   GCD > 1 -> ?(V) #= GCD * ?(_)
 2181                ;   true
 2182                )
 2183            ;   true
 2184            )
 2185        ;   true
 2186        ),
 2187        multiples(Cs, Vs, [C|Left]).
 2188
 2189abs(N, A) :- A is abs(N).
 2190
 2191divide(D, N, R) :- R is N // D.
 2192
 2193scalar_product_(#=, Cs0, Vs, S0) :-
 2194        (   Cs0 = [C|Rest] ->
 2195            gcd(Rest, C, GCD),
 2196            S0 mod GCD =:= 0,
 2197            maplist(divide(GCD), [S0|Cs0], [S|Cs])
 2198        ;   S0 =:= 0, S = S0, Cs = Cs0
 2199        ),
 2200        (   S0 =:= 0 ->
 2201            maplist(abs, Cs, As),
 2202            multiples(As, Vs, [])
 2203        ;   true
 2204        ),
 2205        propagator_init_trigger(Vs, scalar_product_eq(Cs, Vs, S)).
 2206scalar_product_(#\=, Cs, Vs, C) :-
 2207        propagator_init_trigger(Vs, scalar_product_neq(Cs, Vs, C)).
 2208scalar_product_(#=<, Cs, Vs, C) :-
 2209        propagator_init_trigger(Vs, scalar_product_leq(Cs, Vs, C)).
 2210scalar_product_(#<, Cs, Vs, C) :-
 2211        C1 is C - 1,
 2212        scalar_product_(#=<, Cs, Vs, C1).
 2213scalar_product_(#>, Cs, Vs, C) :-
 2214        C1 is C + 1,
 2215        scalar_product_(#>=, Cs, Vs, C1).
 2216scalar_product_(#>=, Cs, Vs, C) :-
 2217        maplist(negative, Cs, Cs1),
 2218        C1 is -C,
 2219        scalar_product_(#=<, Cs1, Vs, C1).
 2220
 2221negative(X0, X) :- X is -X0.
 2222
 2223coeffs_variables_const([], [], [], [], I, I).
 2224coeffs_variables_const([C|Cs], [V|Vs], Cs1, Vs1, I0, I) :-
 2225        (   var(V) ->
 2226            Cs1 = [C|CRest], Vs1 = [V|VRest], I1 = I0
 2227        ;   I1 is I0 + C*V,
 2228            Cs1 = CRest, Vs1 = VRest
 2229        ),
 2230        coeffs_variables_const(Cs, Vs, CRest, VRest, I1, I).
 2231
 2232sum_finite_domains([], [], [], [], Inf, Sup, Inf, Sup).
 2233sum_finite_domains([C|Cs], [V|Vs], Infs, Sups, Inf0, Sup0, Inf, Sup) :-
 2234        fd_get(V, _, Inf1, Sup1, _),
 2235        (   Inf1 = n(NInf) ->
 2236            (   C < 0 ->
 2237                Sup2 is Sup0 + C*NInf
 2238            ;   Inf2 is Inf0 + C*NInf
 2239            ),
 2240            Sups = Sups1,
 2241            Infs = Infs1
 2242        ;   (   C < 0 ->
 2243                Sup2 = Sup0,
 2244                Sups = [C*V|Sups1],
 2245                Infs = Infs1
 2246            ;   Inf2 = Inf0,
 2247                Infs = [C*V|Infs1],
 2248                Sups = Sups1
 2249            )
 2250        ),
 2251        (   Sup1 = n(NSup) ->
 2252            (   C < 0 ->
 2253                Inf2 is Inf0 + C*NSup
 2254            ;   Sup2 is Sup0 + C*NSup
 2255            ),
 2256            Sups1 = Sups2,
 2257            Infs1 = Infs2
 2258        ;   (   C < 0 ->
 2259                Inf2 = Inf0,
 2260                Infs1 = [C*V|Infs2],
 2261                Sups1 = Sups2
 2262            ;   Sup2 = Sup0,
 2263                Sups1 = [C*V|Sups2],
 2264                Infs1 = Infs2
 2265            )
 2266        ),
 2267        sum_finite_domains(Cs, Vs, Infs2, Sups2, Inf2, Sup2, Inf, Sup).
 2268
 2269remove_dist_upper_lower([], _, _, _).
 2270remove_dist_upper_lower([C|Cs], [V|Vs], D1, D2) :-
 2271        (   fd_get(V, VD, VPs) ->
 2272            (   C < 0 ->
 2273                domain_supremum(VD, n(Sup)),
 2274                L is Sup + D1//C,
 2275                domain_remove_smaller_than(VD, L, VD1),
 2276                domain_infimum(VD1, n(Inf)),
 2277                G is Inf - D2//C,
 2278                domain_remove_greater_than(VD1, G, VD2)
 2279            ;   domain_infimum(VD, n(Inf)),
 2280                G is Inf + D1//C,
 2281                domain_remove_greater_than(VD, G, VD1),
 2282                domain_supremum(VD1, n(Sup)),
 2283                L is Sup - D2//C,
 2284                domain_remove_smaller_than(VD1, L, VD2)
 2285            ),
 2286            fd_put(V, VD2, VPs)
 2287        ;   true
 2288        ),
 2289        remove_dist_upper_lower(Cs, Vs, D1, D2).
 2290
 2291
 2292remove_dist_upper_leq([], _, _).
 2293remove_dist_upper_leq([C|Cs], [V|Vs], D1) :-
 2294        (   fd_get(V, VD, VPs) ->
 2295            (   C < 0 ->
 2296                domain_supremum(VD, n(Sup)),
 2297                L is Sup + D1//C,
 2298                domain_remove_smaller_than(VD, L, VD1)
 2299            ;   domain_infimum(VD, n(Inf)),
 2300                G is Inf + D1//C,
 2301                domain_remove_greater_than(VD, G, VD1)
 2302            ),
 2303            fd_put(V, VD1, VPs)
 2304        ;   true
 2305        ),
 2306        remove_dist_upper_leq(Cs, Vs, D1).
 2307
 2308
 2309remove_dist_upper([], _).
 2310remove_dist_upper([C*V|CVs], D) :-
 2311        (   fd_get(V, VD, VPs) ->
 2312            (   C < 0 ->
 2313                (   domain_supremum(VD, n(Sup)) ->
 2314                    L is Sup + D//C,
 2315                    domain_remove_smaller_than(VD, L, VD1)
 2316                ;   VD1 = VD
 2317                )
 2318            ;   (   domain_infimum(VD, n(Inf)) ->
 2319                    G is Inf + D//C,
 2320                    domain_remove_greater_than(VD, G, VD1)
 2321                ;   VD1 = VD
 2322                )
 2323            ),
 2324            fd_put(V, VD1, VPs)
 2325        ;   true
 2326        ),
 2327        remove_dist_upper(CVs, D).
 2328
 2329remove_dist_lower([], _).
 2330remove_dist_lower([C*V|CVs], D) :-
 2331        (   fd_get(V, VD, VPs) ->
 2332            (   C < 0 ->
 2333                (   domain_infimum(VD, n(Inf)) ->
 2334                    G is Inf - D//C,
 2335                    domain_remove_greater_than(VD, G, VD1)
 2336                ;   VD1 = VD
 2337                )
 2338            ;   (   domain_supremum(VD, n(Sup)) ->
 2339                    L is Sup - D//C,
 2340                    domain_remove_smaller_than(VD, L, VD1)
 2341                ;   VD1 = VD
 2342                )
 2343            ),
 2344            fd_put(V, VD1, VPs)
 2345        ;   true
 2346        ),
 2347        remove_dist_lower(CVs, D).
 2348
 2349remove_upper([], _).
 2350remove_upper([C*X|CXs], Max) :-
 2351        (   fd_get(X, XD, XPs) ->
 2352            D is Max//C,
 2353            (   C < 0 ->
 2354                domain_remove_smaller_than(XD, D, XD1)
 2355            ;   domain_remove_greater_than(XD, D, XD1)
 2356            ),
 2357            fd_put(X, XD1, XPs)
 2358        ;   true
 2359        ),
 2360        remove_upper(CXs, Max).
 2361
 2362remove_lower([], _).
 2363remove_lower([C*X|CXs], Min) :-
 2364        (   fd_get(X, XD, XPs) ->
 2365            D is -Min//C,
 2366            (   C < 0 ->
 2367                domain_remove_greater_than(XD, D, XD1)
 2368            ;   domain_remove_smaller_than(XD, D, XD1)
 2369            ),
 2370            fd_put(X, XD1, XPs)
 2371        ;   true
 2372        ),
 2373        remove_lower(CXs, Min).
 2374
 2375%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 2376
 2377/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 2378   Constraint propagation proceeds as follows: Each CLP(FD) variable
 2379   has an attribute that stores its associated domain and constraints.
 2380   Constraints are triggered when the event they are registered for
 2381   occurs (for example: variable is instantiated, bounds change etc.).
 2382   do_queue/0 works off all triggered constraints, possibly triggering
 2383   new ones, until fixpoint.
 2384- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 2385
 2386% FIFO queue
 2387
 2388make_queue :- nb_setval('$clpfd_queue', fast_slow([], [])).
 2389
 2390push_queue(E, Which) :-
 2391        nb_getval('$clpfd_queue', Qs),
 2392        arg(Which, Qs, Q),
 2393        (   Q == [] ->
 2394            setarg(Which, Qs, [E|T]-T)
 2395        ;   Q = H-[E|T],
 2396            setarg(Which, Qs, H-T)
 2397        ).
 2398
 2399pop_queue(E) :-
 2400        nb_getval('$clpfd_queue', Qs),
 2401        (   pop_queue(E, Qs, 1) ->  true
 2402        ;   pop_queue(E, Qs, 2)
 2403        ).
 2404
 2405pop_queue(E, Qs, Which) :-
 2406        arg(Which, Qs, [E|NH]-T),
 2407        (   var(NH) ->
 2408            setarg(Which, Qs, [])
 2409        ;   setarg(Which, Qs, NH-T)
 2410        ).
 2411
 2412fetch_propagator(Prop) :-
 2413        pop_queue(P),
 2414        (   propagator_state(P, S), S == dead -> fetch_propagator(Prop)
 2415        ;   Prop = P
 2416        ).
 2417
 2418/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 2419   Parsing a CLP(FD) expression has two important side-effects: First,
 2420   it constrains the variables occurring in the expression to
 2421   integers. Second, it constrains some of them even more: For
 2422   example, in X/Y and X mod Y, Y is constrained to be #\= 0.
 2423- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 2424
 2425constrain_to_integer(Var) :-
 2426        (   integer(Var) -> true
 2427        ;   fd_get(Var, D, Ps),
 2428            fd_put(Var, D, Ps)
 2429        ).
 2430
 2431power_var_num(P, X, N) :-
 2432        (   var(P) -> X = P, N = 1
 2433        ;   P = Left*Right,
 2434            power_var_num(Left, XL, L),
 2435            power_var_num(Right, XR, R),
 2436            XL == XR,
 2437            X = XL,
 2438            N is L + R
 2439        ).
 2440
 2441/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 2442   Given expression E, we obtain the finite domain variable R by
 2443   interpreting a simple committed-choice language that is a list of
 2444   conditions and bodies. In conditions, g(Goal) means literally Goal,
 2445   and m(Match) means that E can be decomposed as stated. The
 2446   variables are to be understood as the result of parsing the
 2447   subexpressions recursively. In the body, g(Goal) means again Goal,
 2448   and p(Propagator) means to attach and trigger once a propagator.
 2449- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 2450
 2451:- op(800, xfx, =>). 2452
 2453parse_clpfd(E, R,
 2454            [g(cyclic_term(E)) => [g(domain_error(clpfd_expression, E))],
 2455             g(var(E))         => [g(non_monotonic(E)),
 2456                                   g(constrain_to_integer(E)), g(E = R)],
 2457             g(integer(E))     => [g(R = E)],
 2458             ?(E)              => [g(must_be_fd_integer(E)), g(R = E)],
 2459             #(E)              => [g(must_be_fd_integer(E)), g(R = E)],
 2460             m(A+B)            => [p(pplus(A, B, R))],
 2461             % power_var_num/3 must occur before */2 to be useful
 2462             g(power_var_num(E, V, N)) => [p(pexp(V, N, R))],
 2463             m(A*B)            => [p(ptimes(A, B, R))],
 2464             m(A-B)            => [p(pplus(R,B,A))],
 2465             m(-A)             => [p(ptimes(-1,A,R))],
 2466             m(max(A,B))       => [g(A #=< ?(R)), g(B #=< R), p(pmax(A, B, R))],
 2467             m(min(A,B))       => [g(A #>= ?(R)), g(B #>= R), p(pmin(A, B, R))],
 2468             m(A mod B)        => [g(B #\= 0), p(pmod(A, B, R))],
 2469             m(A rem B)        => [g(B #\= 0), p(prem(A, B, R))],
 2470             m(abs(A))         => [g(?(R) #>= 0), p(pabs(A, R))],
 2471%             m(A/B)            => [g(B #\= 0), p(ptzdiv(A, B, R))],
 2472             m(A//B)           => [g(B #\= 0), p(ptzdiv(A, B, R))],
 2473             m(A div B)        => [g(?(R) #= (A - (A mod B)) // B)],
 2474             m(A rdiv B)       => [g(B #\= 0), p(prdiv(A, B, R))],
 2475             m(A^B)            => [p(pexp(A, B, R))],
 2476             % bitwise operations
 2477             m(\A)             => [p(pfunction(\, A, R))],
 2478             m(msb(A))         => [p(pfunction(msb, A, R))],
 2479             m(lsb(A))         => [p(pfunction(lsb, A, R))],
 2480             m(popcount(A))    => [p(pfunction(popcount, A, R))],
 2481             m(A<<B)           => [p(pfunction(<<, A, B, R))],
 2482             m(A>>B)           => [p(pfunction(>>, A, B, R))],
 2483             m(A/\B)           => [p(pfunction(/\, A, B, R))],
 2484             m(A\/B)           => [p(pfunction(\/, A, B, R))],
 2485             m(A xor B)        => [p(pfunction(xor, A, B, R))],
 2486             g(true)           => [g(domain_error(clpfd_expression, E))]
 2487            ]).
 2488
 2489non_monotonic(X) :-
 2490        (   \+ fd_var(X), current_prolog_flag(clpfd_monotonic, true) ->
 2491            instantiation_error(X)
 2492        ;   true
 2493        ).
 2494
 2495% Here, we compile the committed choice language to a single
 2496% predicate, parse_clpfd/2.
 2497
 2498make_parse_clpfd(Clauses) :-
 2499        parse_clpfd_clauses(Clauses0),
 2500        maplist(goals_goal, Clauses0, Clauses).
 2501
 2502goals_goal((Head :- Goals), (Head :- Body)) :-
 2503        list_goal(Goals, Body).
 2504
 2505parse_clpfd_clauses(Clauses) :-
 2506        parse_clpfd(E, R, Matchers),
 2507        maplist(parse_matcher(E, R), Matchers, Clauses).
 2508
 2509parse_matcher(E, R, Matcher, Clause) :-
 2510        Matcher = (Condition0 => Goals0),
 2511        phrase((parse_condition(Condition0, E, Head),
 2512                parse_goals(Goals0)), Goals),
 2513        Clause = (parse_clpfd(Head, R) :- Goals).
 2514
 2515parse_condition(g(Goal), E, E)       --> [Goal, !].
 2516parse_condition(?(E), _, ?(E))       --> [!].
 2517parse_condition(#(E), _, #(E))       --> [!].
 2518parse_condition(m(Match), _, Match0) -->
 2519        [!],
 2520        { copy_term(Match, Match0),
 2521          term_variables(Match0, Vs0),
 2522          term_variables(Match, Vs)
 2523        },
 2524        parse_match_variables(Vs0, Vs).
 2525
 2526parse_match_variables([], []) --> [].
 2527parse_match_variables([V0|Vs0], [V|Vs]) -->
 2528        [parse_clpfd(V0, V)],
 2529        parse_match_variables(Vs0, Vs).
 2530
 2531parse_goals([]) --> [].
 2532parse_goals([G|Gs]) --> parse_goal(G), parse_goals(Gs).
 2533
 2534parse_goal(g(Goal)) --> [Goal].
 2535parse_goal(p(Prop)) -->
 2536        [make_propagator(Prop, P)],
 2537        { term_variables(Prop, Vs) },
 2538        parse_init(Vs, P),
 2539        [trigger_once(P)].
 2540
 2541parse_init([], _)     --> [].
 2542parse_init([V|Vs], P) --> [init_propagator(V, P)], parse_init(Vs, P).
 2543
 2544%?- set_prolog_flag(answer_write_options, [portray(true)]),
 2545%   clpfd:parse_clpfd_clauses(Clauses), maplist(portray_clause, Clauses).
 2546
 2547
 2548%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 2549%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 2550
 2551trigger_once(Prop) :- trigger_prop(Prop), do_queue.
 2552
 2553neq(A, B) :- propagator_init_trigger(pneq(A, B)).
 2554
 2555propagator_init_trigger(P) -->
 2556        { term_variables(P, Vs) },
 2557        propagator_init_trigger(Vs, P).
 2558
 2559propagator_init_trigger(Vs, P) -->
 2560        [p(Prop)],
 2561        { make_propagator(P, Prop),
 2562          maplist(prop_init(Prop), Vs),
 2563          trigger_once(Prop) }.
 2564
 2565propagator_init_trigger(P) :-
 2566        phrase(propagator_init_trigger(P), _).
 2567
 2568propagator_init_trigger(Vs, P) :-
 2569        phrase(propagator_init_trigger(Vs, P), _).
 2570
 2571prop_init(Prop, V) :- init_propagator(V, Prop).
 2572
 2573geq(A, B) :-
 2574        (   fd_get(A, AD, APs) ->
 2575            domain_infimum(AD, AI),
 2576            (   fd_get(B, BD, _) ->
 2577                domain_supremum(BD, BS),
 2578                (   AI cis_geq BS -> true
 2579                ;   propagator_init_trigger(pgeq(A,B))
 2580                )
 2581            ;   (   AI cis_geq n(B) -> true
 2582                ;   domain_remove_smaller_than(AD, B, AD1),
 2583                    fd_put(A, AD1, APs),
 2584                    do_queue
 2585                )
 2586            )
 2587        ;   fd_get(B, BD, BPs) ->
 2588            domain_remove_greater_than(BD, A, BD1),
 2589            fd_put(B, BD1, BPs),
 2590            do_queue
 2591        ;   A >= B
 2592        ).
 2593
 2594/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 2595   Naive parsing of inequalities and disequalities can result in a lot
 2596   of unnecessary work if expressions of non-trivial depth are
 2597   involved: Auxiliary variables are introduced for sub-expressions,
 2598   and propagation proceeds on them as if they were involved in a
 2599   tighter constraint (like equality), whereas eventually only very
 2600   little of the propagated information is actually used. For example,
 2601   only extremal values are of interest in inequalities. Introducing
 2602   auxiliary variables should be avoided when possible, and
 2603   specialised propagators should be used for common constraints.
 2604
 2605   We again use a simple committed-choice language for matching
 2606   special cases of constraints. m_c(M,C) means that M matches and C
 2607   holds. d(X, Y) means decomposition, i.e., it is short for
 2608   g(parse_clpfd(X, Y)). r(X, Y) means to rematch with X and Y.
 2609
 2610   Two things are important: First, although the actual constraint
 2611   functors (#\=2, #=/2 etc.) are used in the description, they must
 2612   expand to the respective auxiliary predicates (match_expand/2)
 2613   because the actual constraints are subject to goal expansion.
 2614   Second, when specialised constraints (like scalar product) post
 2615   simpler constraints on their own, these simpler versions must be
 2616   handled separately and must occur before.
 2617- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 2618
 2619match_expand(#>=, clpfd_geq_).
 2620match_expand(#=, clpfd_equal_).
 2621match_expand(#\=, clpfd_neq).
 2622
 2623symmetric(#=).
 2624symmetric(#\=).
 2625
 2626matches([
 2627         m_c(any(X) #>= any(Y), left_right_linsum_const(X, Y, Cs, Vs, Const)) =>
 2628            [g((   Cs = [1], Vs = [A] -> geq(A, Const)
 2629               ;   Cs = [-1], Vs = [A] -> Const1 is -Const, geq(Const1, A)
 2630               ;   Cs = [1,1], Vs = [A,B] -> ?(A) + ?(B) #= ?(S), geq(S, Const)
 2631               ;   Cs = [1,-1], Vs = [A,B] ->
 2632                   (   Const =:= 0 -> geq(A, B)
 2633                   ;   C1 is -Const,
 2634                       propagator_init_trigger(x_leq_y_plus_c(B, A, C1))
 2635                   )
 2636               ;   Cs = [-1,1], Vs = [A,B] ->
 2637                   (   Const =:= 0 -> geq(B, A)
 2638                   ;   C1 is -Const,
 2639                       propagator_init_trigger(x_leq_y_plus_c(A, B, C1))
 2640                   )
 2641               ;   Cs = [-1,-1], Vs = [A,B] ->
 2642                   ?(A) + ?(B) #= ?(S), Const1 is -Const, geq(Const1, S)
 2643               ;   scalar_product_(#>=, Cs, Vs, Const)
 2644               ))],
 2645         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))],
 2646         m(integer(X) #>= any(Z) + integer(A)) => [g(C is X - A), r(C, Z)],
 2647         m(abs(any(X)-any(Y)) #>= integer(I))  => [d(X, X1), d(Y, Y1), p(absdiff_geq(X1, Y1, I))],
 2648         m(abs(any(X)) #>= integer(I))         => [d(X, RX), g((I>0 -> I1 is -I, RX in inf..I1 \/ I..sup; true))],
 2649         m(integer(I) #>= abs(any(X)))         => [d(X, RX), g(I>=0), g(I1 is -I), g(RX in I1..I)],
 2650         m(any(X) #>= any(Y))                  => [d(X, RX), d(Y, RY), g(geq(RX, RY))],
 2651
 2652         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 2653         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 2654
 2655         m(var(X) #= var(Y))        => [g(constrain_to_integer(X)), g(X=Y)],
 2656         m(var(X) #= var(Y)+var(Z)) => [p(pplus(Y,Z,X))],
 2657         m(var(X) #= var(Y)-var(Z)) => [p(pplus(X,Z,Y))],
 2658         m(var(X) #= var(Y)*var(Z)) => [p(ptimes(Y,Z,X))],
 2659         m(var(X) #= -var(Z))       => [p(ptimes(-1, Z, X))],
 2660         m_c(any(X) #= any(Y), left_right_linsum_const(X, Y, Cs, Vs, S)) =>
 2661            [g(scalar_product_(#=, Cs, Vs, S))],
 2662         m(var(X) #= any(Y))       => [d(Y,X)],
 2663         m(any(X) #= any(Y))       => [d(X, RX), d(Y, RX)],
 2664
 2665         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 2666         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 2667
 2668         m(var(X) #\= integer(Y))             => [g(neq_num(X, Y))],
 2669         m(var(X) #\= var(Y))                 => [g(neq(X,Y))],
 2670         m(var(X) #\= var(Y) + var(Z))        => [p(x_neq_y_plus_z(X, Y, Z))],
 2671         m(var(X) #\= var(Y) - var(Z))        => [p(x_neq_y_plus_z(Y, X, Z))],
 2672         m(var(X) #\= var(Y)*var(Z))          => [p(ptimes(Y,Z,P)), g(neq(X,P))],
 2673         m(integer(X) #\= abs(any(Y)-any(Z))) => [d(Y, Y1), d(Z, Z1), p(absdiff_neq(Y1, Z1, X))],
 2674         m_c(any(X) #\= any(Y), left_right_linsum_const(X, Y, Cs, Vs, S)) =>
 2675            [g(scalar_product_(#\=, Cs, Vs, S))],
 2676         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))],
 2677         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))],
 2678         m(any(X) #\= any(Y)) => [d(X, RX), d(Y, RY), g(neq(RX, RY))]
 2679        ]).
 2680
 2681/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 2682   We again compile the committed-choice matching language to the
 2683   intended auxiliary predicates. We now must take care not to
 2684   unintentionally unify a variable with a complex term.
 2685- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 2686
 2687make_matches(Clauses) :-
 2688        matches(Ms),
 2689        findall(F, (member(M=>_, Ms), arg(1, M, M1), functor(M1, F, _)), Fs0),
 2690        sort(Fs0, Fs),
 2691        maplist(prevent_cyclic_argument, Fs, PrevCyclicClauses),
 2692        phrase(matchers(Ms), Clauses0),
 2693        maplist(goals_goal, Clauses0, MatcherClauses),
 2694        append(PrevCyclicClauses, MatcherClauses, Clauses1),
 2695        sort_by_predicate(Clauses1, Clauses).
 2696
 2697sort_by_predicate(Clauses, ByPred) :-
 2698        map_list_to_pairs(predname, Clauses, Keyed),
 2699        keysort(Keyed, KeyedByPred),
 2700        pairs_values(KeyedByPred, ByPred).
 2701
 2702predname((H:-_), Key)   :- !, predname(H, Key).
 2703predname(M:H, M:Key)    :- !, predname(H, Key).
 2704predname(H, Name/Arity) :- !, functor(H, Name, Arity).
 2705
 2706prevent_cyclic_argument(F0, Clause) :-
 2707        match_expand(F0, F),
 2708        Head =.. [F,X,Y],
 2709        Clause = (Head :- (   cyclic_term(X) ->
 2710                              domain_error(clpfd_expression, X)
 2711                          ;   cyclic_term(Y) ->
 2712                              domain_error(clpfd_expression, Y)
 2713                          ;   false
 2714                          )).
 2715
 2716matchers([]) --> [].
 2717matchers([Condition => Goals|Ms]) -->
 2718        matcher(Condition, Goals),
 2719        matchers(Ms).
 2720
 2721matcher(m(M), Gs) --> matcher(m_c(M,true), Gs).
 2722matcher(m_c(Matcher,Cond), Gs) -->
 2723        [(Head :- Goals0)],
 2724        { Matcher =.. [F,A,B],
 2725          match_expand(F, Expand),
 2726          Head =.. [Expand,X,Y],
 2727          phrase((match(A, X), match(B, Y)), Goals0, [Cond,!|Goals1]),
 2728          phrase(match_goals(Gs, Expand), Goals1) },
 2729        (   { symmetric(F), \+ (subsumes_term(A, B), subsumes_term(B, A)) } ->
 2730            { Head1 =.. [Expand,Y,X] },
 2731            [(Head1 :- Goals0)]
 2732        ;   []
 2733        ).
 2734
 2735match(any(A), T)     --> [A = T].
 2736match(var(V), T)     --> [( nonvar(T), ( T = ?(Var) ; T = #(Var) ) ->
 2737                            must_be_fd_integer(Var), V = Var
 2738                          ; v_or_i(T), V = T
 2739                          )].
 2740match(integer(I), T) --> [integer(T), I = T].
 2741match(-X, T)         --> [nonvar(T), T = -A], match(X, A).
 2742match(abs(X), T)     --> [nonvar(T), T = abs(A)], match(X, A).
 2743match(Binary, T)     -->
 2744        { Binary =.. [Op,X,Y], Term =.. [Op,A,B] },
 2745        [nonvar(T), T = Term],
 2746        match(X, A), match(Y, B).
 2747
 2748match_goals([], _)     --> [].
 2749match_goals([G|Gs], F) --> match_goal(G, F), match_goals(Gs, F).
 2750
 2751match_goal(r(X,Y), F)  --> { G =.. [F,X,Y] }, [G].
 2752match_goal(d(X,Y), _)  --> [parse_clpfd(X, Y)].
 2753match_goal(g(Goal), _) --> [Goal].
 2754match_goal(p(Prop), _) -->
 2755        [make_propagator(Prop, P)],
 2756        { term_variables(Prop, Vs) },
 2757        parse_init(Vs, P),
 2758        [trigger_once(P)].
 2759
 2760
 2761%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 #>=(?X, ?Y)
Same as Y #=< X. When reasoning over integers, replace (>=)/2 by #>=/2 to obtain more general relations. See declarative integer arithmetic.
 2771X #>= Y :- clpfd_geq(X, Y).
 2772
 2773clpfd_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.
 2782X #=< 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.
 2791X #= Y :- clpfd_equal(X, Y).
 2792
 2793clpfd_equal(X, Y) :- clpfd_equal_(X, Y), reinforce(X).
 2794
 2795/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 2796   Conditions under which an equality can be compiled to built-in
 2797   arithmetic. Their order is significant. (/)/2 becomes (//)/2.
 2798- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 2799
 2800expr_conds(E, E)                 --> [integer(E)],
 2801        { var(E), !, \+ current_prolog_flag(clpfd_monotonic, true) }.
 2802expr_conds(E, E)                 --> { integer(E) }.
 2803expr_conds(?(E), E)              --> [integer(E)].
 2804expr_conds(#(E), E)              --> [integer(E)].
 2805expr_conds(-E0, -E)              --> expr_conds(E0, E).
 2806expr_conds(abs(E0), abs(E))      --> expr_conds(E0, E).
 2807expr_conds(A0+B0, A+B)           --> expr_conds(A0, A), expr_conds(B0, B).
 2808expr_conds(A0*B0, A*B)           --> expr_conds(A0, A), expr_conds(B0, B).
 2809expr_conds(A0-B0, A-B)           --> expr_conds(A0, A), expr_conds(B0, B).
 2810expr_conds(A0//B0, A//B)         -->
 2811        expr_conds(A0, A), expr_conds(B0, B),
 2812        [B =\= 0].
 2813%expr_conds(A0/B0, AB)            --> expr_conds(A0//B0, AB).
 2814expr_conds(min(A0,B0), min(A,B)) --> expr_conds(A0, A), expr_conds(B0, B).
 2815expr_conds(max(A0,B0), max(A,B)) --> expr_conds(A0, A), expr_conds(B0, B).
 2816expr_conds(A0 mod B0, A mod B)   -->
 2817        expr_conds(A0, A), expr_conds(B0, B),
 2818        [B =\= 0].
 2819expr_conds(A0^B0, A^B)           -->
 2820        expr_conds(A0, A), expr_conds(B0, B),
 2821        [(B >= 0 ; A =:= -1)].
 2822% Bitwise operations, added to make CLP(FD) usable in more cases
 2823expr_conds(\ A0, \ A) --> expr_conds(A0, A).
 2824expr_conds(A0<<B0, A<<B) --> expr_conds(A0, A), expr_conds(B0, B).
 2825expr_conds(A0>>B0, A>>B) --> expr_conds(A0, A), expr_conds(B0, B).
 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 xor B0, A xor B) --> expr_conds(A0, A), expr_conds(B0, B).
 2829expr_conds(lsb(A0), lsb(A)) --> expr_conds(A0, A).
 2830expr_conds(msb(A0), msb(A)) --> expr_conds(A0, A).
 2831expr_conds(popcount(A0), popcount(A)) --> expr_conds(A0, A).
 2832
 2833:- multifile
 2834        system:goal_expansion/2. 2835:- dynamic
 2836        system:goal_expansion/2. 2837
 2838system:goal_expansion(Goal, Expansion) :-
 2839        \+ current_prolog_flag(clpfd_goal_expansion, false),
 2840        clpfd_expandable(Goal),
 2841        prolog_load_context(module, M),
 2842	(   M == clpfd
 2843	->  true
 2844	;   predicate_property(M:Goal, imported_from(clpfd))
 2845	),
 2846        clpfd_expansion(Goal, Expansion).
 2847
 2848clpfd_expandable(_ in _).
 2849clpfd_expandable(_ #= _).
 2850clpfd_expandable(_ #>= _).
 2851clpfd_expandable(_ #=< _).
 2852clpfd_expandable(_ #> _).
 2853clpfd_expandable(_ #< _).
 2854
 2855clpfd_expansion(Var in Dom, In) :-
 2856        (   ground(Dom), Dom = L..U, integer(L), integer(U) ->
 2857            expansion_simpler(
 2858                (   integer(Var) ->
 2859                    between(L, U, Var)
 2860                ;   clpfd:clpfd_in(Var, Dom)
 2861                ), In)
 2862        ;   In = clpfd:clpfd_in(Var, Dom)
 2863        ).
 2864clpfd_expansion(X0 #= Y0, Equal) :-
 2865        phrase(expr_conds(X0, X), CsX),
 2866        phrase(expr_conds(Y0, Y), CsY),
 2867        list_goal(CsX, CondX),
 2868        list_goal(CsY, CondY),
 2869        expansion_simpler(
 2870                (   CondX ->
 2871                    (   var(Y) -> Y is X
 2872                    ;   CondY -> X =:= Y
 2873                    ;   T is X, clpfd:clpfd_equal(T, Y0)
 2874                    )
 2875                ;   CondY ->
 2876                    (   var(X) -> X is Y
 2877                    ;   T is Y, clpfd:clpfd_equal(X0, T)
 2878                    )
 2879                ;   clpfd:clpfd_equal(X0, Y0)
 2880                ), Equal).
 2881clpfd_expansion(X0 #>= Y0, Geq) :-
 2882        phrase(expr_conds(X0, X), CsX),
 2883        phrase(expr_conds(Y0, Y), CsY),
 2884        list_goal(CsX, CondX),
 2885        list_goal(CsY, CondY),
 2886        expansion_simpler(
 2887              (   CondX ->
 2888                  (   CondY -> X >= Y
 2889                  ;   T is X, clpfd:clpfd_geq(T, Y0)
 2890                  )
 2891              ;   CondY -> T is Y, clpfd:clpfd_geq(X0, T)
 2892              ;   clpfd:clpfd_geq(X0, Y0)
 2893              ), Geq).
 2894clpfd_expansion(X #=< Y,  Leq) :- clpfd_expansion(Y #>= X, Leq).
 2895clpfd_expansion(X #> Y, Gt)    :- clpfd_expansion(X #>= Y+1, Gt).
 2896clpfd_expansion(X #< Y, Lt)    :- clpfd_expansion(Y #> X, Lt).
 2897
 2898expansion_simpler((True->Then0;_), Then) :-
 2899        is_true(True), !,
 2900        expansion_simpler(Then0, Then).
 2901expansion_simpler((False->_;Else0), Else) :-
 2902        is_false(False), !,
 2903        expansion_simpler(Else0, Else).
 2904expansion_simpler((If->Then0;Else0), (If->Then;Else)) :- !,
 2905        expansion_simpler(Then0, Then),
 2906        expansion_simpler(Else0, Else).
 2907expansion_simpler((A0,B0), (A,B)) :-
 2908        expansion_simpler(A0, A),
 2909        expansion_simpler(B0, B).
 2910expansion_simpler(Var is Expr0, Goal) :-
 2911        ground(Expr0), !,
 2912        phrase(expr_conds(Expr0, Expr), Gs),
 2913        (   maplist(call, Gs) -> Value is Expr, Goal = (Var = Value)
 2914        ;   Goal = false
 2915        ).
 2916expansion_simpler(Var =:= Expr0, Goal) :-
 2917        ground(Expr0), !,
 2918        phrase(expr_conds(Expr0, Expr), Gs),
 2919        (   maplist(call, Gs) -> Value is Expr, Goal = (Var =:= Value)
 2920        ;   Goal = false
 2921        ).
 2922expansion_simpler(Var is Expr, Var = Expr) :- var(Expr), !.
 2923expansion_simpler(Var is Expr, Goal) :- !,
 2924        (   var(Var), nonvar(Expr),
 2925            Expr = E mod M, nonvar(E), E = A^B ->
 2926            Goal = ( ( integer(A), integer(B), integer(M),
 2927                       A >= 0, B >= 0, M > 0 ->
 2928                       Var is powm(A, B, M)
 2929                     ; Var is Expr
 2930                     ) )
 2931        ;   Goal = ( Var is Expr )
 2932        ).
 2933expansion_simpler(between(L,U,V), Goal) :- maplist(integer, [L,U,V]), !,
 2934        (   between(L,U,V) -> Goal = true
 2935        ;   Goal = false
 2936        ).
 2937expansion_simpler(Goal, Goal).
 2938
 2939is_true(true).
 2940is_true(integer(I))  :- integer(I).
 2941:- if(current_predicate(var_property/2)). 2942is_true(var(X))      :- var(X), var_property(X, fresh(true)).
 2943is_false(integer(X)) :- var(X), var_property(X, fresh(true)).
 2944is_false((A,B))      :- is_false(A) ; is_false(B).
 2945:- endif. 2946is_false(var(X)) :- nonvar(X).
 2947
 2948
 2949%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 2950
 2951linsum(X, S, S)    --> { var(X), !, non_monotonic(X) }, [vn(X,1)].
 2952linsum(I, S0, S)   --> { integer(I), S is S0 + I }.
 2953linsum(?(X), S, S) --> { must_be_fd_integer(X) }, [vn(X,1)].
 2954linsum(#(X), S, S) --> { must_be_fd_integer(X) }, [vn(X,1)].
 2955linsum(-A, S0, S)  --> mulsum(A, -1, S0, S).
 2956linsum(N*A, S0, S) --> { integer(N) }, !, mulsum(A, N, S0, S).
 2957linsum(A*N, S0, S) --> { integer(N) }, !, mulsum(A, N, S0, S).
 2958linsum(A+B, S0, S) --> linsum(A, S0, S1), linsum(B, S1, S).
 2959linsum(A-B, S0, S) --> linsum(A, S0, S1), mulsum(B, -1, S1, S).
 2960
 2961mulsum(A, M, S0, S) -->
 2962        { phrase(linsum(A, 0, CA), As), S is S0 + M*CA },
 2963        lin_mul(As, M).
 2964
 2965lin_mul([], _)             --> [].
 2966lin_mul([vn(X,N0)|VNs], M) --> { N is N0*M }, [vn(X,N)], lin_mul(VNs, M).
 2967
 2968v_or_i(V) :- var(V), !, non_monotonic(V).
 2969v_or_i(I) :- integer(I).
 2970
 2971must_be_fd_integer(X) :-
 2972        (   var(X) -> constrain_to_integer(X)
 2973        ;   must_be(integer, X)
 2974        ).
 2975
 2976left_right_linsum_const(Left, Right, Cs, Vs, Const) :-
 2977        phrase(linsum(Left, 0, CL), Lefts0, Rights),
 2978        phrase(linsum(Right, 0, CR), Rights0),
 2979        maplist(linterm_negate, Rights0, Rights),
 2980        msort(Lefts0, Lefts),
 2981        Lefts = [vn(First,N)|LeftsRest],
 2982        vns_coeffs_variables(LeftsRest, N, First, Cs0, Vs0),
 2983        filter_linsum(Cs0, Vs0, Cs, Vs),
 2984        Const is CR - CL.
 2985        %format("linear sum: ~w ~w ~w\n", [Cs,Vs,Const]).
 2986
 2987linterm_negate(vn(V,N0), vn(V,N)) :- N is -N0.
 2988
 2989vns_coeffs_variables([], N, V, [N], [V]).
 2990vns_coeffs_variables([vn(V,N)|VNs], N0, V0, Ns, Vs) :-
 2991        (   V == V0 ->
 2992            N1 is N0 + N,
 2993            vns_coeffs_variables(VNs, N1, V0, Ns, Vs)
 2994        ;   Ns = [N0|NRest],
 2995            Vs = [V0|VRest],
 2996            vns_coeffs_variables(VNs, N, V, NRest, VRest)
 2997        ).
 2998
 2999filter_linsum([], [], [], []).
 3000filter_linsum([C0|Cs0], [V0|Vs0], Cs, Vs) :-
 3001        (   C0 =:= 0 ->
 3002            constrain_to_integer(V0),
 3003            filter_linsum(Cs0, Vs0, Cs, Vs)
 3004        ;   Cs = [C0|Cs1], Vs = [V0|Vs1],
 3005            filter_linsum(Cs0, Vs0, Cs1, Vs1)
 3006        ).
 3007
 3008gcd([], G, G).
 3009gcd([N|Ns], G0, G) :-
 3010        G1 is gcd(N, G0),
 3011        gcd(Ns, G1, G).
 3012
 3013even(N) :- N mod 2 =:= 0.
 3014
 3015odd(N) :- \+ even(N).
 3016
 3017/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 3018   k-th root of N, if N is a k-th power.
 3019
 3020   TODO: Replace this when the GMP function becomes available.
 3021- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 3022
 3023integer_kth_root(N, K, R) :-
 3024        (   even(K) ->
 3025            N >= 0
 3026        ;   true
 3027        ),
 3028        (   N < 0 ->
 3029            odd(K),
 3030            integer_kroot(N, 0, N, K, R)
 3031        ;   integer_kroot(0, N, N, K, R)
 3032        ).
 3033
 3034integer_kroot(L, U, N, K, R) :-
 3035        (   L =:= U -> N =:= L^K, R = L
 3036        ;   L + 1 =:= U ->
 3037            (   L^K =:= N -> R = L
 3038            ;   U^K =:= N -> R = U
 3039            ;   false
 3040            )
 3041        ;   Mid is (L + U)//2,
 3042            (   Mid^K > N ->
 3043                integer_kroot(L, Mid, N, K, R)
 3044            ;   integer_kroot(Mid, U, N, K, R)
 3045            )
 3046        ).
 3047
 3048integer_log_b(N, B, Log0, Log) :-
 3049        T is B^Log0,
 3050        (   T =:= N -> Log = Log0
 3051        ;   T < N,
 3052            Log1 is Log0 + 1,
 3053            integer_log_b(N, B, Log1, Log)
 3054        ).
 3055
 3056floor_integer_log_b(N, B, Log0, Log) :-
 3057        T is B^Log0,
 3058        (   T > N -> Log is Log0 - 1
 3059        ;   T =:= N -> Log = Log0
 3060        ;   T < N,
 3061            Log1 is Log0 + 1,
 3062            floor_integer_log_b(N, B, Log1, Log)
 3063        ).
 3064
 3065
 3066/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 3067   Largest R such that R^K =< N.
 3068- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 3069
 3070:- if(current_predicate(nth_integer_root_and_remainder/4)). 3071
 3072% This currently only works for K >= 1, which is all that is needed for now.
 3073integer_kth_root_leq(N, K, R) :-
 3074        nth_integer_root_and_remainder(K, N, R, _).
 3075
 3076:- else. 3077
 3078integer_kth_root_leq(N, K, R) :-
 3079        (   even(K) ->
 3080            N >= 0
 3081        ;   true
 3082        ),
 3083        (   N < 0 ->
 3084            odd(K),
 3085            integer_kroot_leq(N, 0, N, K, R)
 3086        ;   integer_kroot_leq(0, N, N, K, R)
 3087        ).
 3088
 3089integer_kroot_leq(L, U, N, K, R) :-
 3090        (   L =:= U -> R = L
 3091        ;   L + 1 =:= U ->
 3092            (   U^K =< N -> R = U
 3093            ;   R = L
 3094            )
 3095        ;   Mid is (L + U)//2,
 3096            (   Mid^K > N ->
 3097                integer_kroot_leq(L, Mid, N, K, R)
 3098            ;   integer_kroot_leq(Mid, U, N, K, R)
 3099            )
 3100        ).
 3101
 3102:- 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.
 3111X #\= Y :- clpfd_neq(X, Y), do_queue.
 3112
 3113% X #\= Y + Z
 3114
 3115x_neq_y_plus_z(X, Y, Z) :-
 3116        propagator_init_trigger(x_neq_y_plus_z(X,Y,Z)).
 3117
 3118% X is distinct from the number N. This is used internally, and does
 3119% not reinforce other constraints.
 3120
 3121neq_num(X, N) :-
 3122        (   fd_get(X, XD, XPs) ->
 3123            domain_remove(XD, N, XD1),
 3124            fd_put(X, XD1, XPs)
 3125        ;   X =\= N
 3126        ).
 #>(?X, ?Y)
Same as Y #< X. When reasoning over integers, replace (>)/2 by #>/2 to obtain more general relations See declarative integer arithmetic.
 3134X #> 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)].
 3157X #< 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.
 3170#\ 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.
 3210L #<==> R  :- reify(L, B), reify(R, B), do_queue.
 #==>(?P, ?Q)
P implies Q. See reification.
 3216/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 3217   Implication is special in that created auxiliary constraints can be
 3218   retracted when the implication becomes entailed, for example:
 3219
 3220   %?- X + 1 #= Y #==> Z, Z #= 1.
 3221   %@ Z = 1,
 3222   %@ X in inf..sup,
 3223   %@ Y in inf..sup.
 3224
 3225   We cannot use propagator_init_trigger/1 here because the states of
 3226   auxiliary propagators are themselves part of the propagator.
 3227- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 3228
 3229L #==> R   :-
 3230        reify(L, LB, LPs),
 3231        reify(R, RB, RPs),
 3232        append(LPs, RPs, Ps),
 3233        propagator_init_trigger([LB,RB], pimpl(LB,RB,Ps)).
 #<==(?P, ?Q)
Q implies P. See reification.
 3239L #<== R   :- R #==> L.
 #/\(?P, ?Q)
P and Q hold. See reification.
 3245L #/\ R    :- reify(L, 1), reify(R, 1), do_queue.
 3246
 3247conjunctive_neqs_var_drep(Eqs, Var, Drep) :-
 3248        conjunctive_neqs_var(Eqs, Var),
 3249        phrase(conjunctive_neqs_vals(Eqs), Vals),
 3250        list_to_domain(Vals, Dom),
 3251        domain_complement(Dom, C),
 3252        domain_to_drep(C, Drep).
 3253
 3254conjunctive_neqs_var(V, _) :- var(V), !, false.
 3255conjunctive_neqs_var(L #\= R, Var) :-
 3256        (   var(L), integer(R) -> Var = L
 3257        ;   integer(L), var(R) -> Var = R
 3258        ;   false
 3259        ).
 3260conjunctive_neqs_var(A #/\ B, VA) :-
 3261        conjunctive_neqs_var(A, VA),
 3262        conjunctive_neqs_var(B, VB),
 3263        VA == VB.
 3264
 3265conjunctive_neqs_vals(L #\= R) --> ( { integer(L) } -> [L] ; [R] ).
 3266conjunctive_neqs_vals(A #/\ B) -->
 3267        conjunctive_neqs_vals(A),
 3268        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.
 3286L #\/ R :-
 3287        (   disjunctive_eqs_var_drep(L #\/ R, Var, Drep) -> Var in Drep
 3288        ;   reify(L, X, Ps1),
 3289            reify(R, Y, Ps2),
 3290            propagator_init_trigger([X,Y], reified_or(X,Ps1,Y,Ps2,1))
 3291        ).
 3292
 3293disjunctive_eqs_var_drep(Eqs, Var, Drep) :-
 3294        disjunctive_eqs_var(Eqs, Var),
 3295        phrase(disjunctive_eqs_vals(Eqs), Vals),
 3296        list_to_drep(Vals, Drep).
 3297
 3298disjunctive_eqs_var(V, _) :- var(V), !, false.
 3299disjunctive_eqs_var(V in I, V) :- var(V), integer(I).
 3300disjunctive_eqs_var(L #= R, Var) :-
 3301        (   var(L), integer(R) -> Var = L
 3302        ;   integer(L), var(R) -> Var = R
 3303        ;   false
 3304        ).
 3305disjunctive_eqs_var(A #\/ B, VA) :-
 3306        disjunctive_eqs_var(A, VA),
 3307        disjunctive_eqs_var(B, VB),
 3308        VA == VB.
 3309
 3310disjunctive_eqs_vals(L #= R)  --> ( { integer(L) } -> [L] ; [R] ).
 3311disjunctive_eqs_vals(_ in I)  --> [I].
 3312disjunctive_eqs_vals(A #\/ B) -->
 3313        disjunctive_eqs_vals(A),
 3314        disjunctive_eqs_vals(B).
 #\(?P, ?Q)
Either P holds or Q holds, but not both. See reification.
 3321L #\ R :- (L #\/ R) #/\ #\ (L #/\ R).
 3322
 3323/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 3324   A constraint that is being reified need not hold. Therefore, in
 3325   X/Y, Y can as well be 0, for example. Note that it is OK to
 3326   constrain the *result* of an expression (which does not appear
 3327   explicitly in the expression and is not visible to the outside),
 3328   but not the operands, except for requiring that they be integers.
 3329
 3330   In contrast to parse_clpfd/2, the result of an expression can now
 3331   also be undefined, in which case the constraint cannot hold.
 3332   Therefore, the committed-choice language is extended by an element
 3333   d(D) that states D is 1 iff all subexpressions are defined. a(V)
 3334   means that V is an auxiliary variable that was introduced while
 3335   parsing a compound expression. a(X,V) means V is auxiliary unless
 3336   it is ==/2 X, and a(X,Y,V) means V is auxiliary unless it is ==/2 X
 3337   or Y. l(L) means the literal L occurs in the described list.
 3338
 3339   When a constraint becomes entailed or subexpressions become
 3340   undefined, created auxiliary constraints are killed, and the
 3341   "clpfd" attribute is removed from auxiliary variables.
 3342
 3343   For (/)/2, mod/2 and rem/2, we create a skeleton propagator and
 3344   remember it as an auxiliary constraint. The pskeleton propagator
 3345   can use the skeleton when the constraint is defined.
 3346- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 3347
 3348parse_reified(E, R, D,
 3349              [g(cyclic_term(E)) => [g(domain_error(clpfd_expression, E))],
 3350               g(var(E))     => [g(non_monotonic(E)),
 3351                                 g(constrain_to_integer(E)), g(R = E), g(D=1)],
 3352               g(integer(E)) => [g(R=E), g(D=1)],
 3353               ?(E)          => [g(must_be_fd_integer(E)), g(R=E), g(D=1)],
 3354               #(E)          => [g(must_be_fd_integer(E)), g(R=E), g(D=1)],
 3355               m(A+B)        => [d(D), p(pplus(A,B,R)), a(A,B,R)],
 3356               m(A*B)        => [d(D), p(ptimes(A,B,R)), a(A,B,R)],
 3357               m(A-B)        => [d(D), p(pplus(R,B,A)), a(A,B,R)],
 3358               m(-A)         => [d(D), p(ptimes(-1,A,R)), a(R)],
 3359               m(max(A,B))   => [d(D), p(pgeq(R, A)), p(pgeq(R, B)), p(pmax(A,B,R)), a(A,B,R)],
 3360               m(min(A,B))   => [d(D), p(pgeq(A, R)), p(pgeq(B, R)), p(pmin(A,B,R)), a(A,B,R)],
 3361               m(abs(A))     => [g(?(R)#>=0), d(D), p(pabs(A, R)), a(A,R)],
 3362%               m(A/B)        => [skeleton(A,B,D,R,ptzdiv)],
 3363               m(A//B)       => [skeleton(A,B,D,R,ptzdiv)],
 3364               m(A div B)    => [skeleton(A,B,D,R,pdiv)],
 3365               m(A rdiv B)   => [skeleton(A,B,D,R,prdiv)],
 3366               m(A mod B)    => [skeleton(A,B,D,R,pmod)],
 3367               m(A rem B)    => [skeleton(A,B,D,R,prem)],
 3368               m(A^B)        => [d(D), p(pexp(A,B,R)), a(A,B,R)],
 3369               % bitwise operations
 3370               m(\A)         => [function(D,\,A,R)],
 3371               m(msb(A))     => [function(D,msb,A,R)],
 3372               m(lsb(A))     => [function(D,lsb,A,R)],
 3373               m(popcount(A)) => [function(D,popcount,A,R)],
 3374               m(A<<B)       => [function(D,<<,A,B,R)],
 3375               m(A>>B)       => [function(D,>>,A,B,R)],
 3376               m(A/\B)       => [function(D,/\,A,B,R)],
 3377               m(A\/B)       => [function(D,\/,A,B,R)],
 3378               m(A xor B)    => [function(D,xor,A,B,R)],
 3379               g(true)       => [g(domain_error(clpfd_expression, E))]]
 3380             ).
 3381
 3382% Again, we compile this to a predicate, parse_reified_clpfd//3. This
 3383% time, it is a DCG that describes the list of auxiliary variables and
 3384% propagators for the given expression, in addition to relating it to
 3385% its reified (Boolean) finite domain variable and its Boolean
 3386% definedness.
 3387
 3388make_parse_reified(Clauses) :-
 3389        parse_reified_clauses(Clauses0),
 3390        maplist(goals_goal_dcg, Clauses0, Clauses).
 3391
 3392goals_goal_dcg((Head --> Goals), Clause) :-
 3393        list_goal(Goals, Body),
 3394        expand_term((Head --> Body), Clause).
 3395
 3396parse_reified_clauses(Clauses) :-
 3397        parse_reified(E, R, D, Matchers),
 3398        maplist(parse_reified(E, R, D), Matchers, Clauses).
 3399
 3400parse_reified(E, R, D, Matcher, Clause) :-
 3401        Matcher = (Condition0 => Goals0),
 3402        phrase((reified_condition(Condition0, E, Head, Ds),
 3403                reified_goals(Goals0, Ds)), Goals, [a(D)]),
 3404        Clause = (parse_reified_clpfd(Head, R, D) --> Goals).
 3405
 3406reified_condition(g(Goal), E, E, []) --> [{Goal}, !].
 3407reified_condition(?(E), _, ?(E), []) --> [!].
 3408reified_condition(#(E), _, #(E), []) --> [!].
 3409reified_condition(m(Match), _, Match0, Ds) -->
 3410        [!],
 3411        { copy_term(Match, Match0),
 3412          term_variables(Match0, Vs0),
 3413          term_variables(Match, Vs)
 3414        },
 3415        reified_variables(Vs0, Vs, Ds).
 3416
 3417reified_variables([], [], []) --> [].
 3418reified_variables([V0|Vs0], [V|Vs], [D|Ds]) -->
 3419        [parse_reified_clpfd(V0, V, D)],
 3420        reified_variables(Vs0, Vs, Ds).
 3421
 3422reified_goals([], _) --> [].
 3423reified_goals([G|Gs], Ds) --> reified_goal(G, Ds), reified_goals(Gs, Ds).
 3424
 3425reified_goal(d(D), Ds) -->
 3426        (   { Ds = [X] } -> [{D=X}]
 3427        ;   { Ds = [X,Y] } ->
 3428            { phrase(reified_goal(p(reified_and(X,[],Y,[],D)), _), Gs),
 3429              list_goal(Gs, Goal) },
 3430            [( {X==1, Y==1} -> {D = 1} ; Goal )]
 3431        ;   { domain_error(one_or_two_element_list, Ds) }
 3432        ).
 3433reified_goal(g(Goal), _) --> [{Goal}].
 3434reified_goal(p(Vs, Prop), _) -->
 3435        [{make_propagator(Prop, P)}],
 3436        parse_init_dcg(Vs, P),
 3437        [{trigger_once(P)}],
 3438        [( { propagator_state(P, S), S == dead } -> [] ; [p(P)])].
 3439reified_goal(p(Prop), Ds) -->
 3440        { term_variables(Prop, Vs) },
 3441        reified_goal(p(Vs,Prop), Ds).
 3442reified_goal(function(D,Op,A,B,R), Ds) -->
 3443        reified_goals([d(D),p(pfunction(Op,A,B,R)),a(A,B,R)], Ds).
 3444reified_goal(function(D,Op,A,R), Ds) -->
 3445        reified_goals([d(D),p(pfunction(Op,A,R)),a(A,R)], Ds).
 3446reified_goal(skeleton(A,B,D,R,F), Ds) -->
 3447        { Prop =.. [F,X,Y,Z] },
 3448        reified_goals([d(D1),l(p(P)),g(make_propagator(Prop, P)),
 3449                       p([A,B,D2,R], pskeleton(A,B,D2,[X,Y,Z]-P,R,F)),
 3450                       p(reified_and(D1,[],D2,[],D)),a(D2),a(A,B,R)], Ds).
 3451reified_goal(a(V), _)     --> [a(V)].
 3452reified_goal(a(X,V), _)   --> [a(X,V)].
 3453reified_goal(a(X,Y,V), _) --> [a(X,Y,V)].
 3454reified_goal(l(L), _)     --> [[L]].
 3455
 3456parse_init_dcg([], _)     --> [].
 3457parse_init_dcg([V|Vs], P) --> [{init_propagator(V, P)}], parse_init_dcg(Vs, P).
 3458
 3459%?- set_prolog_flag(answer_write_options, [portray(true)]),
 3460%   clpfd:parse_reified_clauses(Cs), maplist(portray_clause, Cs).
 3461
 3462reify(E, B) :- reify(E, B, _).
 3463
 3464reify(Expr, B, Ps) :-
 3465        (   acyclic_term(Expr), reifiable(Expr) -> phrase(reify(Expr, B), Ps)
 3466        ;   domain_error(clpfd_reifiable_expression, Expr)
 3467        ).
 3468
 3469reifiable(E)      :- var(E), non_monotonic(E).
 3470reifiable(E)      :- integer(E), E in 0..1.
 3471reifiable(?(E))   :- must_be_fd_integer(E).
 3472reifiable(#(E))   :- must_be_fd_integer(E).
 3473reifiable(V in _) :- fd_variable(V).
 3474reifiable(Expr)   :-
 3475        Expr =.. [Op,Left,Right],
 3476        (   memberchk(Op, [#>=,#>,#=<,#<,#=,#\=])
 3477        ;   memberchk(Op, [#==>,#<==,#<==>,#/\,#\/,#\]),
 3478            reifiable(Left),
 3479            reifiable(Right)
 3480        ).
 3481reifiable(#\ E) :- reifiable(E).
 3482reifiable(tuples_in(Tuples, Relation)) :-
 3483        must_be(list(list), Tuples),
 3484        maplist(maplist(fd_variable), Tuples),
 3485        must_be(list(list(integer)), Relation).
 3486reifiable(finite_domain(V)) :- fd_variable(V).
 3487
 3488reify(E, B) --> { B in 0..1 }, reify_(E, B).
 3489
 3490reify_(E, B) --> { var(E), !, E = B }.
 3491reify_(E, B) --> { integer(E), E = B }.
 3492reify_(?(B), B) --> [].
 3493reify_(#(B), B) --> [].
 3494reify_(V in Drep, B) -->
 3495        { drep_to_domain(Drep, Dom) },
 3496        propagator_init_trigger(reified_in(V,Dom,B)),
 3497        a(B).
 3498reify_(tuples_in(Tuples, Relation), B) -->
 3499        { maplist(relation_tuple_b_prop(Relation), Tuples, Bs, Ps),
 3500          maplist(monotonic, Bs, Bs1),
 3501          fold_statement(conjunction, Bs1, And),
 3502          ?(B) #<==> And },
 3503        propagator_init_trigger([B], tuples_not_in(Tuples, Relation, B)),
 3504        kill_reified_tuples(Bs, Ps, Bs),
 3505        list(Ps),
 3506        as([B|Bs]).
 3507reify_(finite_domain(V), B) -->
 3508        propagator_init_trigger(reified_fd(V,B)),
 3509        a(B).
 3510reify_(L #>= R, B) --> arithmetic(L, R, B, reified_geq).
 3511reify_(L #= R, B)  --> arithmetic(L, R, B, reified_eq).
 3512reify_(L #\= R, B) --> arithmetic(L, R, B, reified_neq).
 3513reify_(L #> R, B)  --> reify_(L #>= (R+1), B).
 3514reify_(L #=< R, B) --> reify_(R #>= L, B).
 3515reify_(L #< R, B)  --> reify_(R #>= (L+1), B).
 3516reify_(L #==> R, B)  --> reify_((#\ L) #\/ R, B).
 3517reify_(L #<== R, B)  --> reify_(R #==> L, B).
 3518reify_(L #<==> R, B) --> reify_((L #==> R) #/\ (R #==> L), B).
 3519reify_(L #\ R, B) --> reify_((L #\/ R) #/\ #\ (L #/\ R), B).
 3520reify_(L #/\ R, B)   -->
 3521        (   { conjunctive_neqs_var_drep(L #/\ R, V, D) } -> reify_(V in D, B)
 3522        ;   boolean(L, R, B, reified_and)
 3523        ).
 3524reify_(L #\/ R, B) -->
 3525        (   { disjunctive_eqs_var_drep(L #\/ R, V, D) } -> reify_(V in D, B)
 3526        ;   boolean(L, R, B, reified_or)
 3527        ).
 3528reify_(#\ Q, B) -->
 3529        reify(Q, QR),
 3530        propagator_init_trigger(reified_not(QR,B)),
 3531        a(B).
 3532
 3533arithmetic(L, R, B, Functor) -->
 3534        { phrase((parse_reified_clpfd(L, LR, LD),
 3535                  parse_reified_clpfd(R, RR, RD)), Ps),
 3536          Prop =.. [Functor,LD,LR,RD,RR,Ps,B] },
 3537        list(Ps),
 3538        propagator_init_trigger([LD,LR,RD,RR,B], Prop),
 3539        a(B).
 3540
 3541boolean(L, R, B, Functor) -->
 3542        { reify(L, LR, Ps1), reify(R, RR, Ps2),
 3543          Prop =.. [Functor,LR,Ps1,RR,Ps2,B] },
 3544        list(Ps1), list(Ps2),
 3545        propagator_init_trigger([LR,RR,B], Prop),
 3546        a(LR, RR, B).
 3547
 3548list([])     --> [].
 3549list([L|Ls]) --> [L], list(Ls).
 3550
 3551a(X,Y,B) -->
 3552        (   { nonvar(X) } -> a(Y, B)
 3553        ;   { nonvar(Y) } -> a(X, B)
 3554        ;   [a(X,Y,B)]
 3555        ).
 3556
 3557a(X, B) -->
 3558        (   { var(X) } -> [a(X, B)]
 3559        ;   a(B)
 3560        ).
 3561
 3562a(B) -->
 3563        (   { var(B) } -> [a(B)]
 3564        ;   []
 3565        ).
 3566
 3567as([])     --> [].
 3568as([B|Bs]) --> a(B), as(Bs).
 3569
 3570kill_reified_tuples([], _, _) --> [].
 3571kill_reified_tuples([B|Bs], Ps, All) -->
 3572        propagator_init_trigger([B], kill_reified_tuples(B, Ps, All)),
 3573        kill_reified_tuples(Bs, Ps, All).
 3574
 3575relation_tuple_b_prop(Relation, Tuple, B, p(Prop)) :-
 3576        put_attr(R, clpfd_relation, Relation),
 3577        make_propagator(reified_tuple_in(Tuple, R, B), Prop),
 3578        tuple_freeze_(Tuple, Prop),
 3579        init_propagator(B, Prop).
 3580
 3581
 3582tuples_in_conjunction(Tuples, Relation, Conj) :-
 3583        maplist(tuple_in_disjunction(Relation), Tuples, Disjs),
 3584        fold_statement(conjunction, Disjs, Conj).
 3585
 3586tuple_in_disjunction(Relation, Tuple, Disj) :-
 3587        maplist(tuple_in_conjunction(Tuple), Relation, Conjs),
 3588        fold_statement(disjunction, Conjs, Disj).
 3589
 3590tuple_in_conjunction(Tuple, Element, Conj) :-
 3591        maplist(var_eq, Tuple, Element, Eqs),
 3592        fold_statement(conjunction, Eqs, Conj).
 3593
 3594fold_statement(Operation, List, Statement) :-
 3595        (   List = [] -> Statement = 1
 3596        ;   List = [First|Rest],
 3597            foldl(Operation, Rest, First, Statement)
 3598        ).
 3599
 3600conjunction(E, Conj, Conj #/\ E).
 3601
 3602disjunction(E, Disj, Disj #\/ E).
 3603
 3604var_eq(V, N, ?(V) #= N).
 3605
 3606% Match variables to created skeleton.
 3607
 3608skeleton(Vs, Vs-Prop) :-
 3609        maplist(prop_init(Prop), Vs),
 3610        trigger_once(Prop).
 3611
 3612/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 3613   A drep is a user-accessible and visible domain representation. N,
 3614   N..M, and D1 \/ D2 are dreps, if D1 and D2 are dreps.
 3615- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 3616
 3617is_drep(N)      :- integer(N).
 3618is_drep(N..M)   :- drep_bound(N), drep_bound(M), N \== sup, M \== inf.
 3619is_drep(D1\/D2) :- is_drep(D1), is_drep(D2).
 3620is_drep({AI})   :- is_and_integers(AI).
 3621is_drep(\D)     :- is_drep(D).
 3622
 3623is_and_integers(I)     :- integer(I).
 3624is_and_integers((A,B)) :- is_and_integers(A), is_and_integers(B).
 3625
 3626drep_bound(I)   :- integer(I).
 3627drep_bound(sup).
 3628drep_bound(inf).
 3629
 3630drep_to_intervals(I)        --> { integer(I) }, [n(I)-n(I)].
 3631drep_to_intervals(N..M)     -->
 3632        (   { defaulty_to_bound(N, N1), defaulty_to_bound(M, M1),
 3633              N1 cis_leq M1} -> [N1-M1]
 3634        ;   []
 3635        ).
 3636drep_to_intervals(D1 \/ D2) -->
 3637        drep_to_intervals(D1), drep_to_intervals(D2).
 3638drep_to_intervals(\D0) -->
 3639        { drep_to_domain(D0, D1),
 3640          domain_complement(D1, D),
 3641          domain_to_drep(D, Drep) },
 3642        drep_to_intervals(Drep).
 3643drep_to_intervals({AI}) -->
 3644        and_integers_(AI).
 3645
 3646and_integers_(I)     --> { integer(I) }, [n(I)-n(I)].
 3647and_integers_((A,B)) --> and_integers_(A), and_integers_(B).
 3648
 3649drep_to_domain(DR, D) :-
 3650        must_be(ground, DR),
 3651        (   is_drep(DR) -> true
 3652        ;   domain_error(clpfd_domain, DR)
 3653        ),
 3654        phrase(drep_to_intervals(DR), Is0),
 3655        merge_intervals(Is0, Is1),
 3656        intervals_to_domain(Is1, D).
 3657
 3658merge_intervals(Is0, Is) :-
 3659        keysort(Is0, Is1),
 3660        merge_overlapping(Is1, Is).
 3661
 3662merge_overlapping([], []).
 3663merge_overlapping([A-B0|ABs0], [A-B|ABs]) :-
 3664        merge_remaining(ABs0, B0, B, Rest),
 3665        merge_overlapping(Rest, ABs).
 3666
 3667merge_remaining([], B, B, []).
 3668merge_remaining([N-M|NMs], B0, B, Rest) :-
 3669        Next cis B0 + n(1),
 3670        (   N cis_gt Next -> B = B0, Rest = [N-M|NMs]
 3671        ;   B1 cis max(B0,M),
 3672            merge_remaining(NMs, B1, B, Rest)
 3673        ).
 3674
 3675domain(V, Dom) :-
 3676        (   fd_get(V, Dom0, VPs) ->
 3677            domains_intersection(Dom, Dom0, Dom1),
 3678            %format("intersected\n: ~w\n ~w\n==> ~w\n\n", [Dom,Dom0,Dom1]),
 3679            fd_put(V, Dom1, VPs),
 3680            do_queue,
 3681            reinforce(V)
 3682        ;   domain_contains(Dom, V)
 3683        ).
 3684
 3685domains([], _).
 3686domains([V|Vs], D) :- domain(V, D), domains(Vs, D).
 3687
 3688props_number(fd_props(Gs,Bs,Os), N) :-
 3689        length(Gs, N1),
 3690        length(Bs, N2),
 3691        length(Os, N3),
 3692        N is N1 + N2 + N3.
 3693
 3694fd_get(X, Dom, Ps) :-
 3695        (   get_attr(X, clpfd, Attr) -> Attr = clpfd_attr(_,_,_,Dom,Ps)
 3696        ;   var(X) -> default_domain(Dom), Ps = fd_props([],[],[])
 3697        ).
 3698
 3699fd_get(X, Dom, Inf, Sup, Ps) :-
 3700        fd_get(X, Dom, Ps),
 3701        domain_infimum(Dom, Inf),
 3702        domain_supremum(Dom, Sup).
 3703
 3704/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 3705   By default, propagation always terminates. Currently, this is
 3706   ensured by allowing the left and right boundaries, as well as the
 3707   distance between the smallest and largest number occurring in the
 3708   domain representation to be changed at most once after a constraint
 3709   is posted, unless the domain is bounded. Set the experimental
 3710   Prolog flag 'clpfd_propagation' to 'full' to make the solver
 3711   propagate as much as possible. This can make queries
 3712   non-terminating, like: X #> abs(X), or: X #> Y, Y #> X, X #> 0.
 3713   Importantly, it can also make labeling non-terminating, as in:
 3714
 3715   ?- B #==> X #> abs(X), indomain(B).
 3716- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 3717
 3718fd_put(X, Dom, Ps) :-
 3719        (   current_prolog_flag(clpfd_propagation, full) ->
 3720            put_full(X, Dom, Ps)
 3721        ;   put_terminating(X, Dom, Ps)
 3722        ).
 3723
 3724put_terminating(X, Dom, Ps) :-
 3725        Dom \== empty,
 3726        (   Dom = from_to(F, F) -> F = n(X)
 3727        ;   (   get_attr(X, clpfd, Attr) ->
 3728                Attr = clpfd_attr(Left,Right,Spread,OldDom, _OldPs),
 3729                put_attr(X, clpfd, clpfd_attr(Left,Right,Spread,Dom,Ps)),
 3730                (   OldDom == Dom -> true
 3731                ;   (   Left == (.) -> Bounded = yes
 3732                    ;   domain_infimum(Dom, Inf), domain_supremum(Dom, Sup),
 3733                        (   Inf = n(_), Sup = n(_) ->
 3734                            Bounded = yes
 3735                        ;   Bounded = no
 3736                        )
 3737                    ),
 3738                    (   Bounded == yes ->
 3739                        put_attr(X, clpfd, clpfd_attr(.,.,.,Dom,Ps)),
 3740                        trigger_props(Ps, X, OldDom, Dom)
 3741                    ;   % infinite domain; consider border and spread changes
 3742                        domain_infimum(OldDom, OldInf),
 3743                        (   Inf == OldInf -> LeftP = Left
 3744                        ;   LeftP = yes
 3745                        ),
 3746                        domain_supremum(OldDom, OldSup),
 3747                        (   Sup == OldSup -> RightP = Right
 3748                        ;   RightP = yes
 3749                        ),
 3750                        domain_spread(OldDom, OldSpread),
 3751                        domain_spread(Dom, NewSpread),
 3752                        (   NewSpread == OldSpread -> SpreadP = Spread
 3753                        ;   NewSpread cis_lt OldSpread -> SpreadP = no
 3754                        ;   SpreadP = yes
 3755                        ),
 3756                        put_attr(X, clpfd, clpfd_attr(LeftP,RightP,SpreadP,Dom,Ps)),
 3757                        (   RightP == yes, Right = yes -> true
 3758                        ;   LeftP == yes, Left = yes -> true
 3759                        ;   SpreadP == yes, Spread = yes -> true
 3760                        ;   trigger_props(Ps, X, OldDom, Dom)
 3761                        )
 3762                    )
 3763                )
 3764            ;   var(X) ->
 3765                put_attr(X, clpfd, clpfd_attr(no,no,no,Dom, Ps))
 3766            ;   true
 3767            )
 3768        ).
 3769
 3770domain_spread(Dom, Spread) :-
 3771        domain_smallest_finite(Dom, S),
 3772        domain_largest_finite(Dom, L),
 3773        Spread cis L - S.
 3774
 3775smallest_finite(inf, Y, Y).
 3776smallest_finite(n(N), _, n(N)).
 3777
 3778domain_smallest_finite(from_to(F,T), S)   :- smallest_finite(F, T, S).
 3779domain_smallest_finite(split(_, L, _), S) :- domain_smallest_finite(L, S).
 3780
 3781largest_finite(sup, Y, Y).
 3782largest_finite(n(N), _, n(N)).
 3783
 3784domain_largest_finite(from_to(F,T), L)   :- largest_finite(T, F, L).
 3785domain_largest_finite(split(_, _, R), L) :- domain_largest_finite(R, L).
 3786
 3787/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 3788   With terminating propagation, all relevant constraints get a
 3789   propagation opportunity whenever a new constraint is posted.
 3790- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 3791
 3792reinforce(X) :-
 3793        (   current_prolog_flag(clpfd_propagation, full) ->
 3794            % full propagation propagates everything in any case
 3795            true
 3796        ;   term_variables(X, Vs),
 3797            maplist(reinforce_, Vs),
 3798            do_queue
 3799        ).
 3800
 3801reinforce_(X) :-
 3802        (   fd_var(X), fd_get(X, Dom, Ps) ->
 3803            put_full(X, Dom, Ps)
 3804        ;   true
 3805        ).
 3806
 3807put_full(X, Dom, Ps) :-
 3808        Dom \== empty,
 3809        (   Dom = from_to(F, F) -> F = n(X)
 3810        ;   (   get_attr(X, clpfd, Attr) ->
 3811                Attr = clpfd_attr(_,_,_,OldDom, _OldPs),
 3812                put_attr(X, clpfd, clpfd_attr(no,no,no,Dom, Ps)),
 3813                %format("putting dom: ~w\n", [Dom]),
 3814                (   OldDom == Dom -> true
 3815                ;   trigger_props(Ps, X, OldDom, Dom)
 3816                )
 3817            ;   var(X) -> %format('\t~w in ~w .. ~w\n',[X,L,U]),
 3818                put_attr(X, clpfd, clpfd_attr(no,no,no,Dom, Ps))
 3819            ;   true
 3820            )
 3821        ).
 3822
 3823/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 3824   A propagator is a term of the form propagator(C, State), where C
 3825   represents a constraint, and State is a free variable that can be
 3826   used to destructively change the state of the propagator via
 3827   attributes. This can be used to avoid redundant invocation of the
 3828   same propagator, or to disable the propagator.
 3829- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 3830
 3831make_propagator(C, propagator(C, _)).
 3832
 3833propagator_state(propagator(_,S), S).
 3834
 3835trigger_props(fd_props(Gs,Bs,Os), X, D0, D) :-
 3836        (   ground(X) ->
 3837            trigger_props_(Gs),
 3838            trigger_props_(Bs)
 3839        ;   Bs \== [] ->
 3840            domain_infimum(D0, I0),
 3841            domain_infimum(D, I),
 3842            (   I == I0 ->
 3843                domain_supremum(D0, S0),
 3844                domain_supremum(D, S),
 3845                (   S == S0 -> true
 3846                ;   trigger_props_(Bs)
 3847                )
 3848            ;   trigger_props_(Bs)
 3849            )
 3850        ;   true
 3851        ),
 3852        trigger_props_(Os).
 3853
 3854trigger_props(fd_props(Gs,Bs,Os), X) :-
 3855        trigger_props_(Os),
 3856        trigger_props_(Bs),
 3857        (   ground(X) ->
 3858            trigger_props_(Gs)
 3859        ;   true
 3860        ).
 3861
 3862trigger_props(fd_props(Gs,Bs,Os)) :-
 3863        trigger_props_(Gs),
 3864        trigger_props_(Bs),
 3865        trigger_props_(Os).
 3866
 3867trigger_props_([]).
 3868trigger_props_([P|Ps]) :- trigger_prop(P), trigger_props_(Ps).
 3869
 3870trigger_prop(Propagator) :-
 3871        propagator_state(Propagator, State),
 3872        (   State == dead -> true
 3873        ;   get_attr(State, clpfd_aux, queued) -> true
 3874        ;   b_getval('$clpfd_current_propagator', C), C == State -> true
 3875        ;   % passive
 3876            % format("triggering: ~w\n", [Propagator]),
 3877            put_attr(State, clpfd_aux, queued),
 3878            (   arg(1, Propagator, C), functor(C, F, _), global_constraint(F) ->
 3879                push_queue(Propagator, 2)
 3880            ;   push_queue(Propagator, 1)
 3881            )
 3882        ).
 3883
 3884kill(State) :- del_attr(State, clpfd_aux), State = dead.
 3885
 3886kill(State, Ps) :-
 3887        kill(State),
 3888        maplist(kill_entailed, Ps).
 3889
 3890kill_entailed(p(Prop)) :-
 3891        propagator_state(Prop, State),
 3892        kill(State).
 3893kill_entailed(a(V)) :-
 3894        del_attr(V, clpfd).
 3895kill_entailed(a(X,B)) :-
 3896        (   X == B -> true
 3897        ;   del_attr(B, clpfd)
 3898        ).
 3899kill_entailed(a(X,Y,B)) :-
 3900        (   X == B -> true
 3901        ;   Y == B -> true
 3902        ;   del_attr(B, clpfd)
 3903        ).
 3904
 3905no_reactivation(rel_tuple(_,_)).
 3906no_reactivation(pdistinct(_)).
 3907no_reactivation(pgcc(_,_,_)).
 3908no_reactivation(pgcc_single(_,_)).
 3909%no_reactivation(scalar_product(_,_,_,_)).
 3910
 3911activate_propagator(propagator(P,State)) :-
 3912        % format("running: ~w\n", [P]),
 3913        del_attr(State, clpfd_aux),
 3914        (   no_reactivation(P) ->
 3915            b_setval('$clpfd_current_propagator', State),
 3916            run_propagator(P, State),
 3917            b_setval('$clpfd_current_propagator', [])
 3918        ;   run_propagator(P, State)
 3919        ).
 3920
 3921disable_queue :- b_setval('$clpfd_queue_status', disabled).
 3922enable_queue  :- b_setval('$clpfd_queue_status', enabled).
 3923
 3924portray_propagator(propagator(P,_), F) :- functor(P, F, _).
 3925
 3926portray_queue(V, []) :- var(V), !.
 3927portray_queue([P|Ps], [F|Fs]) :-
 3928        portray_propagator(P, F),
 3929        portray_queue(Ps, Fs).
 3930
 3931do_queue :-
 3932        % b_getval('$clpfd_queue', H-_),
 3933        % portray_queue(H, Port),
 3934        % format("queue: ~w\n", [Port]),
 3935        (   b_getval('$clpfd_queue_status', enabled) ->
 3936            (   fetch_propagator(Propagator) ->
 3937                activate_propagator(Propagator),
 3938                do_queue
 3939            ;   true
 3940            )
 3941        ;   true
 3942        ).
 3943
 3944init_propagator(Var, Prop) :-
 3945        (   fd_get(Var, Dom, Ps0) ->
 3946            insert_propagator(Prop, Ps0, Ps),
 3947            fd_put(Var, Dom, Ps)
 3948        ;   true
 3949        ).
 3950
 3951constraint_wake(pneq, ground).
 3952constraint_wake(x_neq_y_plus_z, ground).
 3953constraint_wake(absdiff_neq, ground).
 3954constraint_wake(pdifferent, ground).
 3955constraint_wake(pexclude, ground).
 3956constraint_wake(scalar_product_neq, ground).
 3957
 3958constraint_wake(x_leq_y_plus_c, bounds).
 3959constraint_wake(scalar_product_eq, bounds).
 3960constraint_wake(scalar_product_leq, bounds).
 3961constraint_wake(pplus, bounds).
 3962constraint_wake(pgeq, bounds).
 3963constraint_wake(pgcc_single, bounds).
 3964constraint_wake(pgcc_check_single, bounds).
 3965
 3966global_constraint(pdistinct).
 3967global_constraint(pgcc).
 3968global_constraint(pgcc_single).
 3969global_constraint(pcircuit).
 3970%global_constraint(rel_tuple).
 3971%global_constraint(scalar_product_eq).
 3972
 3973insert_propagator(Prop, Ps0, Ps) :-
 3974        Ps0 = fd_props(Gs,Bs,Os),
 3975        arg(1, Prop, Constraint),
 3976        functor(Constraint, F, _),
 3977        (   constraint_wake(F, ground) ->
 3978            Ps = fd_props([Prop|Gs], Bs, Os)
 3979        ;   constraint_wake(F, bounds) ->
 3980            Ps = fd_props(Gs, [Prop|Bs], Os)
 3981        ;   Ps = fd_props(Gs, Bs, [Prop|Os])
 3982        ).
 3983
 3984%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 lex_chain(+Lists)
Lists are lexicographically non-decreasing.
 3990lex_chain(Lss) :-
 3991        must_be(list(list), Lss),
 3992        maplist(maplist(fd_variable), Lss),
 3993        (   Lss == [] -> true
 3994        ;   Lss = [First|Rest],
 3995            make_propagator(presidual(lex_chain(Lss)), Prop),
 3996            foldl(lex_chain_(Prop), Rest, First, _)
 3997        ).
 3998
 3999lex_chain_(Prop, Ls, Prev, Ls) :-
 4000        maplist(prop_init(Prop), Ls),
 4001        lex_le(Prev, Ls).
 4002
 4003lex_le([], []).
 4004lex_le([V1|V1s], [V2|V2s]) :-
 4005        ?(V1) #=< ?(V2),
 4006        (   integer(V1) ->
 4007            (   integer(V2) ->
 4008                (   V1 =:= V2 -> lex_le(V1s, V2s) ;  true )
 4009            ;   freeze(V2, lex_le([V1|V1s], [V2|V2s]))
 4010            )
 4011        ;   freeze(V1, lex_le([V1|V1s], [V2|V2s]))
 4012        ).
 4013
 4014%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 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]].
 4060tuples_in(Tuples, Relation) :-
 4061        must_be(list(list), Tuples),
 4062        maplist(maplist(fd_variable), Tuples),
 4063        must_be(list(list(integer)), Relation),
 4064        maplist(relation_tuple(Relation), Tuples),
 4065        do_queue.
 4066
 4067relation_tuple(Relation, Tuple) :-
 4068        relation_unifiable(Relation, Tuple, Us, _, _),
 4069        (   ground(Tuple) -> memberchk(Tuple, Relation)
 4070        ;   tuple_domain(Tuple, Us),
 4071            (   Tuple = [_,_|_] -> tuple_freeze(Tuple, Us)
 4072            ;   true
 4073            )
 4074        ).
 4075
 4076tuple_domain([], _).
 4077tuple_domain([T|Ts], Relation0) :-
 4078        maplist(list_first_rest, Relation0, Firsts, Relation1),
 4079        (   var(T) ->
 4080            (   Firsts = [Unique] -> T = Unique
 4081            ;   list_to_domain(Firsts, FDom),
 4082                fd_get(T, TDom, TPs),
 4083                domains_intersection(TDom, FDom, TDom1),
 4084                fd_put(T, TDom1, TPs)
 4085            )
 4086        ;   true
 4087        ),
 4088        tuple_domain(Ts, Relation1).
 4089
 4090tuple_freeze(Tuple, Relation) :-
 4091        put_attr(R, clpfd_relation, Relation),
 4092        make_propagator(rel_tuple(R, Tuple), Prop),
 4093        tuple_freeze_(Tuple, Prop).
 4094
 4095tuple_freeze_([], _).
 4096tuple_freeze_([T|Ts], Prop) :-
 4097        (   var(T) ->
 4098            init_propagator(T, Prop),
 4099            trigger_prop(Prop)
 4100        ;   true
 4101        ),
 4102        tuple_freeze_(Ts, Prop).
 4103
 4104relation_unifiable([], _, [], Changed, Changed).
 4105relation_unifiable([R|Rs], Tuple, Us, Changed0, Changed) :-
 4106        (   all_in_domain(R, Tuple) ->
 4107            Us = [R|Rest],
 4108            relation_unifiable(Rs, Tuple, Rest, Changed0, Changed)
 4109        ;   relation_unifiable(Rs, Tuple, Us, true, Changed)
 4110        ).
 4111
 4112all_in_domain([], []).
 4113all_in_domain([A|As], [T|Ts]) :-
 4114        (   fd_get(T, Dom, _) ->
 4115            domain_contains(Dom, A)
 4116        ;   T =:= A
 4117        ),
 4118        all_in_domain(As, Ts).
 4119
 4120%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 4121
 4122% trivial propagator, used only to remember pending constraints
 4123run_propagator(presidual(_), _).
 4124
 4125%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 4126run_propagator(pdifferent(Left,Right,X,_), MState) :-
 4127        run_propagator(pexclude(Left,Right,X), MState).
 4128
 4129run_propagator(weak_distinct(Left,Right,X,_), _MState) :-
 4130        (   ground(X) ->
 4131            disable_queue,
 4132            exclude_fire(Left, Right, X),
 4133            enable_queue
 4134        ;   outof_reducer(Left, Right, X)
 4135            %(   var(X) -> kill_if_isolated(Left, Right, X, MState)
 4136            %;   true
 4137            %)
 4138        ).
 4139
 4140run_propagator(pexclude(Left,Right,X), _) :-
 4141        (   ground(X) ->
 4142            disable_queue,
 4143            exclude_fire(Left, Right, X),
 4144            enable_queue
 4145        ;   true
 4146        ).
 4147
 4148run_propagator(pdistinct(Ls), _MState) :-
 4149        distinct(Ls).
 4150
 4151run_propagator(check_distinct(Left,Right,X), _) :-
 4152        \+ list_contains(Left, X),
 4153        \+ list_contains(Right, X).
 4154
 4155%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 4156
 4157run_propagator(pelement(N, Is, V), MState) :-
 4158        (   fd_get(N, NDom, _) ->
 4159            (   fd_get(V, VDom, VPs) ->
 4160                integers_remaining(Is, 1, NDom, empty, VDom1),
 4161                domains_intersection(VDom, VDom1, VDom2),
 4162                fd_put(V, VDom2, VPs)
 4163            ;   true
 4164            )
 4165        ;   kill(MState), nth1(N, Is, V)
 4166        ).
 4167
 4168%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 4169
 4170run_propagator(pgcc_single(Vs, Pairs), _) :- gcc_global(Vs, Pairs).
 4171
 4172run_propagator(pgcc_check_single(Pairs), _) :- gcc_check(Pairs).
 4173
 4174run_propagator(pgcc_check(Pairs), _) :- gcc_check(Pairs).
 4175
 4176run_propagator(pgcc(Vs, _, Pairs), _) :- gcc_global(Vs, Pairs).
 4177
 4178%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 4179
 4180run_propagator(pcircuit(Vs), _MState) :-
 4181        distinct(Vs),
 4182        propagate_circuit(Vs).
 4183
 4184%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 4185run_propagator(pneq(A, B), MState) :-
 4186        (   nonvar(A) ->
 4187            (   nonvar(B) -> A =\= B, kill(MState)
 4188            ;   fd_get(B, BD0, BExp0),
 4189                domain_remove(BD0, A, BD1),
 4190                kill(MState),
 4191                fd_put(B, BD1, BExp0)
 4192            )
 4193        ;   nonvar(B) -> run_propagator(pneq(B, A), MState)
 4194        ;   A \== B,
 4195            fd_get(A, _, AI, AS, _), fd_get(B, _, BI, BS, _),
 4196            (   AS cis_lt BI -> kill(MState)
 4197            ;   AI cis_gt BS -> kill(MState)
 4198            ;   true
 4199            )
 4200        ).
 4201
 4202%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 4203run_propagator(pgeq(A,B), MState) :-
 4204        (   A == B -> kill(MState)
 4205        ;   nonvar(A) ->
 4206            (   nonvar(B) -> kill(MState), A >= B
 4207            ;   fd_get(B, BD, BPs),
 4208                domain_remove_greater_than(BD, A, BD1),
 4209                kill(MState),
 4210                fd_put(B, BD1, BPs)
 4211            )
 4212        ;   nonvar(B) ->
 4213            fd_get(A, AD, APs),
 4214            domain_remove_smaller_than(AD, B, AD1),
 4215            kill(MState),
 4216            fd_put(A, AD1, APs)
 4217        ;   fd_get(A, AD, AL, AU, APs),
 4218            fd_get(B, _, BL, BU, _),
 4219            AU cis_geq BL,
 4220            (   AL cis_geq BU -> kill(MState)
 4221            ;   AU == BL -> kill(MState), A = B
 4222            ;   NAL cis max(AL,BL),
 4223                domains_intersection(AD, from_to(NAL,AU), NAD),
 4224                fd_put(A, NAD, APs),
 4225                (   fd_get(B, BD2, BL2, BU2, BPs2) ->
 4226                    NBU cis min(BU2, AU),
 4227                    domains_intersection(BD2, from_to(BL2,NBU), NBD),
 4228                    fd_put(B, NBD, BPs2)
 4229                ;   true
 4230                )
 4231            )
 4232        ).
 4233
 4234%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 4235
 4236run_propagator(rel_tuple(R, Tuple), MState) :-
 4237        get_attr(R, clpfd_relation, Relation),
 4238        (   ground(Tuple) -> kill(MState), memberchk(Tuple, Relation)
 4239        ;   relation_unifiable(Relation, Tuple, Us, false, Changed),
 4240            Us = [_|_],
 4241            (   Tuple = [First,Second], ( ground(First) ; ground(Second) ) ->
 4242                kill(MState)
 4243            ;   true
 4244            ),
 4245            (   Us = [Single] -> kill(MState), Single = Tuple
 4246            ;   Changed ->
 4247                put_attr(R, clpfd_relation, Us),
 4248                disable_queue,
 4249                tuple_domain(Tuple, Us),
 4250                enable_queue
 4251            ;   true
 4252            )
 4253        ).
 4254
 4255%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 4256
 4257run_propagator(pserialized(S_I, D_I, S_J, D_J, _), MState) :-
 4258        (   nonvar(S_I), nonvar(S_J) ->
 4259            kill(MState),
 4260            (   S_I + D_I =< S_J -> true
 4261            ;   S_J + D_J =< S_I -> true
 4262            ;   false
 4263            )
 4264        ;   serialize_lower_upper(S_I, D_I, S_J, D_J, MState),
 4265            serialize_lower_upper(S_J, D_J, S_I, D_I, MState)
 4266        ).
 4267
 4268%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 4269
 4270% abs(X-Y) #\= C
 4271run_propagator(absdiff_neq(X,Y,C), MState) :-
 4272        (   C < 0 -> kill(MState)
 4273        ;   nonvar(X) ->
 4274            kill(MState),
 4275            (   nonvar(Y) -> abs(X - Y) =\= C
 4276            ;   V1 is X - C, neq_num(Y, V1),
 4277                V2 is C + X, neq_num(Y, V2)
 4278            )
 4279        ;   nonvar(Y) -> kill(MState),
 4280            V1 is C + Y, neq_num(X, V1),
 4281            V2 is Y - C, neq_num(X, V2)
 4282        ;   true
 4283        ).
 4284
 4285% abs(X-Y) #>= C
 4286run_propagator(absdiff_geq(X,Y,C), MState) :-
 4287        (   C =< 0 -> kill(MState)
 4288        ;   nonvar(X) ->
 4289            kill(MState),
 4290            (   nonvar(Y) -> abs(X-Y) >= C
 4291            ;   P1 is X - C, P2 is X + C,
 4292                Y in inf..P1 \/ P2..sup
 4293            )
 4294        ;   nonvar(Y) ->
 4295            kill(MState),
 4296            P1 is Y - C, P2 is Y + C,
 4297            X in inf..P1 \/ P2..sup
 4298        ;   true
 4299        ).
 4300
 4301% X #\= Y + Z
 4302run_propagator(x_neq_y_plus_z(X,Y,Z), MState) :-
 4303        (   nonvar(X) ->
 4304            (   nonvar(Y) ->
 4305                (   nonvar(Z) -> kill(MState), X =\= Y + Z
 4306                ;   kill(MState), XY is X - Y, neq_num(Z, XY)
 4307                )
 4308            ;   nonvar(Z) -> kill(MState), XZ is X - Z, neq_num(Y, XZ)
 4309            ;   true
 4310            )
 4311        ;   nonvar(Y) ->
 4312            (   nonvar(Z) ->
 4313                kill(MState), YZ is Y + Z, neq_num(X, YZ)
 4314            ;   Y =:= 0 -> kill(MState), neq(X, Z)
 4315            ;   true
 4316            )
 4317        ;   Z == 0 -> kill(MState), neq(X, Y)
 4318        ;   true
 4319        ).
 4320
 4321% X #=< Y + C
 4322run_propagator(x_leq_y_plus_c(X,Y,C), MState) :-
 4323        (   nonvar(X) ->
 4324            (   nonvar(Y) -> kill(MState), X =< Y + C
 4325            ;   kill(MState),
 4326                R is X - C,
 4327                fd_get(Y, YD, YPs),
 4328                domain_remove_smaller_than(YD, R, YD1),
 4329                fd_put(Y, YD1, YPs)
 4330            )
 4331        ;   nonvar(Y) ->
 4332            kill(MState),
 4333            R is Y + C,
 4334            fd_get(X, XD, XPs),
 4335            domain_remove_greater_than(XD, R, XD1),
 4336            fd_put(X, XD1, XPs)
 4337        ;   (   X == Y -> C >= 0, kill(MState)
 4338            ;   fd_get(Y, YD, _),
 4339                (   domain_supremum(YD, n(YSup)) ->
 4340                    YS1 is YSup + C,
 4341                    fd_get(X, XD, XPs),
 4342                    domain_remove_greater_than(XD, YS1, XD1),
 4343                    fd_put(X, XD1, XPs)
 4344                ;   true
 4345                ),
 4346                (   fd_get(X, XD2, _), domain_infimum(XD2, n(XInf)) ->
 4347                    XI1 is XInf - C,
 4348                    (   fd_get(Y, YD1, YPs1) ->
 4349                        domain_remove_smaller_than(YD1, XI1, YD2),
 4350                        (   domain_infimum(YD2, n(YInf)),
 4351                            domain_supremum(XD2, n(XSup)),
 4352                            XSup =< YInf + C ->
 4353                            kill(MState)
 4354                        ;   true
 4355                        ),
 4356                        fd_put(Y, YD2, YPs1)
 4357                    ;   true
 4358                    )
 4359                ;   true
 4360                )
 4361            )
 4362        ).
 4363
 4364run_propagator(scalar_product_neq(Cs0,Vs0,P0), MState) :-
 4365        coeffs_variables_const(Cs0, Vs0, Cs, Vs, 0, I),
 4366        P is P0 - I,
 4367        (   Vs = [] -> kill(MState), P =\= 0
 4368        ;   Vs = [V], Cs = [C] ->
 4369            kill(MState),
 4370            (   C =:= 1 -> neq_num(V, P)
 4371            ;   C*V #\= P
 4372            )
 4373        ;   Cs == [1,-1] -> kill(MState), Vs = [A,B], x_neq_y_plus_z(A, B, P)
 4374        ;   Cs == [-1,1] -> kill(MState), Vs = [A,B], x_neq_y_plus_z(B, A, P)
 4375        ;   P =:= 0, Cs = [1,1,-1] ->
 4376            kill(MState), Vs = [A,B,C], x_neq_y_plus_z(C, A, B)
 4377        ;   P =:= 0, Cs = [1,-1,1] ->
 4378            kill(MState), Vs = [A,B,C], x_neq_y_plus_z(B, A, C)
 4379        ;   P =:= 0, Cs = [-1,1,1] ->
 4380            kill(MState), Vs = [A,B,C], x_neq_y_plus_z(A, B, C)
 4381        ;   true
 4382        ).
 4383
 4384run_propagator(scalar_product_leq(Cs0,Vs0,P0), MState) :-
 4385        coeffs_variables_const(Cs0, Vs0, Cs, Vs, 0, I),
 4386        P is P0 - I,
 4387        (   Vs = [] -> kill(MState), P >= 0
 4388        ;   sum_finite_domains(Cs, Vs, Infs, Sups, 0, 0, Inf, Sup),
 4389            D1 is P - Inf,
 4390            disable_queue,
 4391            (   Infs == [], Sups == [] ->
 4392                Inf =< P,
 4393                (   Sup =< P -> kill(MState)
 4394                ;   remove_dist_upper_leq(Cs, Vs, D1)
 4395                )
 4396            ;   Infs == [] -> Inf =< P, remove_dist_upper(Sups, D1)
 4397            ;   Sups = [_], Infs = [_] ->
 4398                remove_upper(Infs, D1)
 4399            ;   Infs = [_] -> remove_upper(Infs, D1)
 4400            ;   true
 4401            ),
 4402            enable_queue
 4403        ).
 4404
 4405run_propagator(scalar_product_eq(Cs0,Vs0,P0), MState) :-
 4406        coeffs_variables_const(Cs0, Vs0, Cs, Vs, 0, I),
 4407        P is P0 - I,
 4408        (   Vs = [] -> kill(MState), P =:= 0
 4409        ;   Vs = [V], Cs = [C] -> kill(MState), P mod C =:= 0, V is P // C
 4410        ;   Cs == [1,1] -> kill(MState), Vs = [A,B], A + B #= P
 4411        ;   Cs == [1,-1] -> kill(MState), Vs = [A,B], A #= P + B
 4412        ;   Cs == [-1,1] -> kill(MState), Vs = [A,B], B #= P + A
 4413        ;   Cs == [-1,-1] -> kill(MState), Vs = [A,B], P1 is -P, A + B #= P1
 4414        ;   P =:= 0, Cs == [1,1,-1] -> kill(MState), Vs = [A,B,C], A + B #= C
 4415        ;   P =:= 0, Cs == [1,-1,1] -> kill(MState), Vs = [A,B,C], A + C #= B
 4416        ;   P =:= 0, Cs == [-1,1,1] -> kill(MState), Vs = [A,B,C], B + C #= A
 4417        ;   sum_finite_domains(Cs, Vs, Infs, Sups, 0, 0, Inf, Sup),
 4418            % nl, writeln(Infs-Sups-Inf-Sup),
 4419            D1 is P - Inf,
 4420            D2 is Sup - P,
 4421            disable_queue,
 4422            (   Infs == [], Sups == [] ->
 4423                between(Inf, Sup, P),
 4424                remove_dist_upper_lower(Cs, Vs, D1, D2)
 4425            ;   Sups = [] -> P =< Sup, remove_dist_lower(Infs, D2)
 4426            ;   Infs = [] -> Inf =< P, remove_dist_upper(Sups, D1)
 4427            ;   Sups = [_], Infs = [_] ->
 4428                remove_lower(Sups, D2),
 4429                remove_upper(Infs, D1)
 4430            ;   Infs = [_] -> remove_upper(Infs, D1)
 4431            ;   Sups = [_] -> remove_lower(Sups, D2)
 4432            ;   true
 4433            ),
 4434            enable_queue
 4435        ).
 4436
 4437% X + Y = Z
 4438run_propagator(pplus(X,Y,Z), MState) :-
 4439        (   nonvar(X) ->
 4440            (   X =:= 0 -> kill(MState), Y = Z
 4441            ;   Y == Z -> kill(MState), X =:= 0
 4442            ;   nonvar(Y) -> kill(MState), Z is X + Y
 4443            ;   nonvar(Z) -> kill(MState), Y is Z - X
 4444            ;   fd_get(Z, ZD, ZPs),
 4445                fd_get(Y, YD, _),
 4446                domain_shift(YD, X, Shifted_YD),
 4447                domains_intersection(ZD, Shifted_YD, ZD1),
 4448                fd_put(Z, ZD1, ZPs),
 4449                (   fd_get(Y, YD1, YPs) ->
 4450                    O is -X,
 4451                    domain_shift(ZD1, O, YD2),
 4452                    domains_intersection(YD1, YD2, YD3),
 4453                    fd_put(Y, YD3, YPs)
 4454                ;   true
 4455                )
 4456            )
 4457        ;   nonvar(Y) -> run_propagator(pplus(Y,X,Z), MState)
 4458        ;   nonvar(Z) ->
 4459            (   X == Y -> kill(MState), even(Z), X is Z // 2
 4460            ;   fd_get(X, XD, _),
 4461                fd_get(Y, YD, YPs),
 4462                domain_negate(XD, XDN),
 4463                domain_shift(XDN, Z, YD1),
 4464                domains_intersection(YD, YD1, YD2),
 4465                fd_put(Y, YD2, YPs),
 4466                (   fd_get(X, XD1, XPs) ->
 4467                    domain_negate(YD2, YD2N),
 4468                    domain_shift(YD2N, Z, XD2),
 4469                    domains_intersection(XD1, XD2, XD3),
 4470                    fd_put(X, XD3, XPs)
 4471                ;   true
 4472                )
 4473            )
 4474        ;   (   X == Y -> kill(MState), 2*X #= Z
 4475            ;   X == Z -> kill(MState), Y = 0
 4476            ;   Y == Z -> kill(MState), X = 0
 4477            ;   fd_get(X, XD, XL, XU, XPs),
 4478                fd_get(Y, _, YL, YU, _),
 4479                fd_get(Z, _, ZL, ZU, _),
 4480                NXL cis max(XL, ZL-YU),
 4481                NXU cis min(XU, ZU-YL),
 4482                update_bounds(X, XD, XPs, XL, XU, NXL, NXU),
 4483                (   fd_get(Y, YD2, YL2, YU2, YPs2) ->
 4484                    NYL cis max(YL2, ZL-NXU),
 4485                    NYU cis min(YU2, ZU-NXL),
 4486                    update_bounds(Y, YD2, YPs2, YL2, YU2, NYL, NYU)
 4487                ;   NYL = n(Y), NYU = n(Y)
 4488                ),
 4489                (   fd_get(Z, ZD2, ZL2, ZU2, ZPs2) ->
 4490                    NZL cis max(ZL2,NXL+NYL),
 4491                    NZU cis min(ZU2,NXU+NYU),
 4492                    update_bounds(Z, ZD2, ZPs2, ZL2, ZU2, NZL, NZU)
 4493                ;   true
 4494                )
 4495            )
 4496        ).
 4497
 4498%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 4499
 4500run_propagator(ptimes(X,Y,Z), MState) :-
 4501        (   nonvar(X) ->
 4502            (   nonvar(Y) -> kill(MState), Z is X * Y
 4503            ;   X =:= 0 -> kill(MState), Z = 0
 4504            ;   X =:= 1 -> kill(MState), Z = Y
 4505            ;   nonvar(Z) -> kill(MState), 0 =:= Z mod X, Y is Z // X
 4506            ;   (   Y == Z -> kill(MState), Y = 0
 4507                ;   fd_get(Y, YD, _),
 4508                    fd_get(Z, ZD, ZPs),
 4509                    domain_expand(YD, X, Scaled_YD),
 4510                    domains_intersection(ZD, Scaled_YD, ZD1),
 4511                    fd_put(Z, ZD1, ZPs),
 4512                    (   fd_get(Y, YDom2, YPs2) ->
 4513                        domain_contract(ZD1, X, Contract),
 4514                        domains_intersection(YDom2, Contract, NYDom),
 4515                        fd_put(Y, NYDom, YPs2)
 4516                    ;   kill(MState), Z is X * Y
 4517                    )
 4518                )
 4519            )
 4520        ;   nonvar(Y) -> run_propagator(ptimes(Y,X,Z), MState)
 4521        ;   nonvar(Z) ->
 4522            (   X == Y ->
 4523                kill(MState),
 4524                integer_kth_root(Z, 2, R),
 4525                NR is -R,
 4526                X in NR \/ R
 4527            ;   fd_get(X, XD, XL, XU, XPs),
 4528                fd_get(Y, YD, YL, YU, _),
 4529                min_max_factor(n(Z), n(Z), YL, YU, XL, XU, NXL, NXU),
 4530                update_bounds(X, XD, XPs, XL, XU, NXL, NXU),
 4531                (   fd_get(Y, YD2, YL2, YU2, YPs2) ->
 4532                    min_max_factor(n(Z), n(Z), NXL, NXU, YL2, YU2, NYL, NYU),
 4533                    update_bounds(Y, YD2, YPs2, YL2, YU2, NYL, NYU)
 4534                ;   (   Y =\= 0 -> 0 =:= Z mod Y, kill(MState), X is Z // Y
 4535                    ;   kill(MState), Z = 0
 4536                    )
 4537                ),
 4538                (   Z =:= 0 ->
 4539                    (   \+ domain_contains(XD, 0) -> kill(MState), Y = 0
 4540                    ;   \+ domain_contains(YD, 0) -> kill(MState), X = 0
 4541                    ;   true
 4542                    )
 4543                ;  neq_num(X, 0), neq_num(Y, 0)
 4544                )
 4545            )
 4546        ;   (   X == Y -> kill(MState), X^2 #= Z
 4547            ;   fd_get(X, XD, XL, XU, XPs),
 4548                fd_get(Y, _, YL, YU, _),
 4549                fd_get(Z, ZD, ZL, ZU, _),
 4550                (   Y == Z, \+ domain_contains(ZD, 0) -> kill(MState), X = 1
 4551                ;   X == Z, \+ domain_contains(ZD, 0) -> kill(MState), Y = 1
 4552                ;   min_max_factor(ZL, ZU, YL, YU, XL, XU, NXL, NXU),
 4553                    update_bounds(X, XD, XPs, XL, XU, NXL, NXU),
 4554                    (   fd_get(Y, YD2, YL2, YU2, YPs2) ->
 4555                        min_max_factor(ZL, ZU, NXL, NXU, YL2, YU2, NYL, NYU),
 4556                        update_bounds(Y, YD2, YPs2, YL2, YU2, NYL, NYU)
 4557                    ;   NYL = n(Y), NYU = n(Y)
 4558                    ),
 4559                    (   fd_get(Z, ZD2, ZL2, ZU2, ZPs2) ->
 4560                        min_product(NXL, NXU, NYL, NYU, NZL),
 4561                        max_product(NXL, NXU, NYL, NYU, NZU),
 4562                        (   NZL cis_leq ZL2, NZU cis_geq ZU2 -> ZD3 = ZD2
 4563                        ;   domains_intersection(ZD2, from_to(NZL,NZU), ZD3),
 4564                            fd_put(Z, ZD3, ZPs2)
 4565                        ),
 4566                        (   domain_contains(ZD3, 0) ->  true
 4567                        ;   neq_num(X, 0), neq_num(Y, 0)
 4568                        )
 4569                    ;   true
 4570                    )
 4571                )
 4572            )
 4573        ).
 4574
 4575%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 4576
 4577% X div Y = Z
 4578run_propagator(pdiv(X,Y,Z), MState) :- kill(MState), Z #= (X-(X mod Y)) // Y.
 4579
 4580% X rdiv Y = Z
 4581run_propagator(prdiv(X,Y,Z), MState) :- kill(MState), Z*Y #= X.
 4582
 4583% X // Y = Z (round towards zero)
 4584run_propagator(ptzdiv(X,Y,Z), MState) :-
 4585        (   nonvar(X) ->
 4586            (   nonvar(Y) -> kill(MState), Y =\= 0, Z is X // Y
 4587            ;   fd_get(Y, YD, YL, YU, YPs),
 4588                (   nonvar(Z) ->
 4589                    (   Z =:= 0 ->
 4590                        NYL is -abs(X) - 1,
 4591                        NYU is abs(X) + 1,
 4592                        domains_intersection(YD, split(0, from_to(inf,n(NYL)),
 4593                                                       from_to(n(NYU), sup)),
 4594                                             NYD),
 4595                        fd_put(Y, NYD, YPs)
 4596                    ;   (   sign(X) =:= sign(Z) ->
 4597                            NYL cis max(n(X) // (n(Z)+sign(n(Z))) + n(1), YL),
 4598                            NYU cis min(n(X) // n(Z), YU)
 4599                        ;   NYL cis max(n(X) // n(Z), YL),
 4600                            NYU cis min(n(X) // (n(Z)+sign(n(Z))) - n(1), YU)
 4601                        ),
 4602                        update_bounds(Y, YD, YPs, YL, YU, NYL, NYU)
 4603                    )
 4604                ;   fd_get(Z, ZD, ZL, ZU, ZPs),
 4605                    (   X >= 0, ( YL cis_gt n(0) ; YU cis_lt n(0) )->
 4606                        NZL cis max(n(X)//YU, ZL),
 4607                        NZU cis min(n(X)//YL, ZU)
 4608                    ;   X < 0, ( YL cis_gt n(0) ; YU cis_lt n(0) ) ->
 4609                        NZL cis max(n(X)//YL, ZL),
 4610                        NZU cis min(n(X)//YU, ZU)
 4611                    ;   % TODO: more stringent bounds, cover Y
 4612                        NZL cis max(-abs(n(X)), ZL),
 4613                        NZU cis min(abs(n(X)), ZU)
 4614                    ),
 4615                    update_bounds(Z, ZD, ZPs, ZL, ZU, NZL, NZU),
 4616                    (   X >= 0, NZL cis_gt n(0), fd_get(Y, YD1, YPs1) ->
 4617                        NYL cis n(X) // (NZU + n(1)) + n(1),
 4618                        NYU cis n(X) // NZL,
 4619                        domains_intersection(YD1, from_to(NYL, NYU), NYD1),
 4620                        fd_put(Y, NYD1, YPs1)
 4621                    ;   true
 4622                    )
 4623                )
 4624            )
 4625        ;   nonvar(Y) ->
 4626            Y =\= 0,
 4627            (   Y =:= 1 -> kill(MState), X = Z
 4628            ;   Y =:= -1 -> kill(MState), Z #= -X
 4629            ;   fd_get(X, XD, XL, XU, XPs),
 4630                (   nonvar(Z) ->
 4631                    kill(MState),
 4632                    (   sign(Z) =:= sign(Y) ->
 4633                        NXL cis max(n(Z)*n(Y), XL),
 4634                        NXU cis min((abs(n(Z))+n(1))*abs(n(Y))-n(1), XU)
 4635                    ;   Z =:= 0 ->
 4636                        NXL cis max(-abs(n(Y)) + n(1), XL),
 4637                        NXU cis min(abs(n(Y)) - n(1), XU)
 4638                    ;   NXL cis max((n(Z)+sign(n(Z)))*n(Y)+n(1), XL),
 4639                        NXU cis min(n(Z)*n(Y), XU)
 4640                    ),
 4641                    update_bounds(X, XD, XPs, XL, XU, NXL, NXU)
 4642                ;   fd_get(Z, ZD, ZPs),
 4643                    domain_contract_less(XD, Y, Contracted),
 4644                    domains_intersection(ZD, Contracted, NZD),
 4645                    fd_put(Z, NZD, ZPs),
 4646                    (   fd_get(X, XD2, XPs2) ->
 4647                        domain_expand_more(NZD, Y, Expanded),
 4648                        domains_intersection(XD2, Expanded, NXD2),
 4649                        fd_put(X, NXD2, XPs2)
 4650                    ;   true
 4651                    )
 4652                )
 4653            )
 4654        ;   nonvar(Z) ->
 4655            fd_get(X, XD, XL, XU, XPs),
 4656            fd_get(Y, _, YL, YU, _),
 4657            (   YL cis_geq n(0), XL cis_geq n(0) ->
 4658                NXL cis max(YL*n(Z), XL),
 4659                NXU cis min(YU*(n(Z)+n(1))-n(1), XU)
 4660            ;   %TODO: cover more cases
 4661                NXL = XL, NXU = XU
 4662            ),
 4663            update_bounds(X, XD, XPs, XL, XU, NXL, NXU)
 4664        ;   (   X == Y -> kill(MState), Z = 1
 4665            ;   fd_get(X, _, XL, XU, _),
 4666                fd_get(Y, _, YL, _, _),
 4667                fd_get(Z, ZD, ZPs),
 4668                NZU cis max(abs(XL), XU),
 4669                NZL cis -NZU,
 4670                domains_intersection(ZD, from_to(NZL,NZU), NZD0),
 4671                (   XL cis_geq n(0), YL cis_geq n(0) ->
 4672                    domain_remove_smaller_than(NZD0, 0, NZD1)
 4673                ;   % TODO: cover more cases
 4674                    NZD1 = NZD0
 4675                ),
 4676                fd_put(Z, NZD1, ZPs)
 4677            )
 4678        ).
 4679
 4680
 4681%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 4682% Y = abs(X)
 4683
 4684run_propagator(pabs(X,Y), MState) :-
 4685        (   nonvar(X) -> kill(MState), Y is abs(X)
 4686        ;   nonvar(Y) ->
 4687            kill(MState),
 4688            Y >= 0,
 4689            YN is -Y,
 4690            X in YN \/ Y
 4691        ;   fd_get(X, XD, XPs),
 4692            fd_get(Y, YD, _),
 4693            domain_negate(YD, YDNegative),
 4694            domains_union(YD, YDNegative, XD1),
 4695            domains_intersection(XD, XD1, XD2),
 4696            fd_put(X, XD2, XPs),
 4697            (   fd_get(Y, YD1, YPs1) ->
 4698                domain_negate(XD2, XD2Neg),
 4699                domains_union(XD2, XD2Neg, YD2),
 4700                domain_remove_smaller_than(YD2, 0, YD3),
 4701                domains_intersection(YD1, YD3, YD4),
 4702                fd_put(Y, YD4, YPs1)
 4703            ;   true
 4704            )
 4705        ).
 4706%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 4707% Z = X mod Y
 4708
 4709run_propagator(pmod(X,Y,Z), MState) :-
 4710        (   nonvar(X) ->
 4711            (   nonvar(Y) -> kill(MState), Y =\= 0, Z is X mod Y
 4712            ;   true
 4713            )
 4714        ;   nonvar(Y) ->
 4715            Y =\= 0,
 4716            (   abs(Y) =:= 1 -> kill(MState), Z = 0
 4717            ;   var(Z) ->
 4718                YP is abs(Y) - 1,
 4719                (   Y > 0, fd_get(X, _, n(XL), n(XU), _) ->
 4720                    (   XL >= 0, XU < Y ->
 4721                        kill(MState), Z = X, ZL = XL, ZU = XU
 4722                    ;   ZL = 0, ZU = YP
 4723                    )
 4724                ;   Y > 0 -> ZL = 0, ZU = YP
 4725                ;   YN is -YP, ZL = YN, ZU = 0
 4726                ),
 4727                (   fd_get(Z, ZD, ZPs) ->
 4728                    domains_intersection(ZD, from_to(n(ZL), n(ZU)), ZD1),
 4729                    domain_infimum(ZD1, n(ZMin)),
 4730                    domain_supremum(ZD1, n(ZMax)),
 4731                    fd_put(Z, ZD1, ZPs)
 4732                ;   ZMin = Z, ZMax = Z
 4733                ),
 4734                (   fd_get(X, XD, XPs), domain_infimum(XD, n(XMin)) ->
 4735                    Z1 is XMin mod Y,
 4736                    (   between(ZMin, ZMax, Z1) -> true
 4737                    ;   Y > 0 ->
 4738                        Next is ((XMin - ZMin + Y - 1) div Y)*Y + ZMin,
 4739                        domain_remove_smaller_than(XD, Next, XD1),
 4740                        fd_put(X, XD1, XPs)
 4741                    ;   neq_num(X, XMin)
 4742                    )
 4743                ;   true
 4744                ),
 4745                (   fd_get(X, XD2, XPs2), domain_supremum(XD2, n(XMax)) ->
 4746                    Z2 is XMax mod Y,
 4747                    (   between(ZMin, ZMax, Z2) -> true
 4748                    ;   Y > 0 ->
 4749                        Prev is ((XMax - ZMin) div Y)*Y + ZMax,
 4750                        domain_remove_greater_than(XD2, Prev, XD3),
 4751                        fd_put(X, XD3, XPs2)
 4752                    ;   neq_num(X, XMax)
 4753                    )
 4754                ;   true
 4755                )
 4756            ;   fd_get(X, XD, XPs),
 4757                % if possible, propagate at the boundaries
 4758                (   domain_infimum(XD, n(Min)) ->
 4759                    (   Min mod Y =:= Z -> true
 4760                    ;   Y > 0 ->
 4761                        Next is ((Min - Z + Y - 1) div Y)*Y + Z,
 4762                        domain_remove_smaller_than(XD, Next, XD1),
 4763                        fd_put(X, XD1, XPs)
 4764                    ;   neq_num(X, Min)
 4765                    )
 4766                ;   true
 4767                ),
 4768                (   fd_get(X, XD2, XPs2) ->
 4769                    (   domain_supremum(XD2, n(Max)) ->
 4770                        (   Max mod Y =:= Z -> true
 4771                        ;   Y > 0 ->
 4772                            Prev is ((Max - Z) div Y)*Y + Z,
 4773                            domain_remove_greater_than(XD2, Prev, XD3),
 4774                            fd_put(X, XD3, XPs2)
 4775                        ;   neq_num(X, Max)
 4776                        )
 4777                    ;   true
 4778                    )
 4779                ;   true
 4780                )
 4781            )
 4782        ;   X == Y -> kill(MState), Z = 0
 4783        ;   true % TODO: propagate more
 4784        ).
 4785
 4786%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 4787% Z = X rem Y
 4788
 4789run_propagator(prem(X,Y,Z), MState) :-
 4790        (   nonvar(X) ->
 4791            (   nonvar(Y) -> kill(MState), Y =\= 0, Z is X rem Y
 4792            ;   U is abs(X),
 4793                fd_get(Y, YD, _),
 4794                (   X >=0, domain_infimum(YD, n(Min)), Min >= 0 -> L = 0
 4795                ;   L is -U
 4796                ),
 4797                Z in L..U
 4798            )
 4799        ;   nonvar(Y) ->
 4800            Y =\= 0,
 4801            (   abs(Y) =:= 1 -> kill(MState), Z = 0
 4802            ;   var(Z) ->
 4803                YP is abs(Y) - 1,
 4804                YN is -YP,
 4805                (   Y > 0, fd_get(X, _, n(XL), n(XU), _) ->
 4806                    (   abs(XL) < Y, XU < Y -> kill(MState), Z = X, ZL = XL
 4807                    ;   XL < 0, abs(XL) < Y -> ZL = XL
 4808                    ;   XL >= 0 -> ZL = 0
 4809                    ;   ZL = YN
 4810                    ),
 4811                    (   XU > 0, XU < Y -> ZU = XU
 4812                    ;   XU < 0 -> ZU = 0
 4813                    ;   ZU = YP
 4814                    )
 4815                ;   ZL = YN, ZU = YP
 4816                ),
 4817                (   fd_get(Z, ZD, ZPs) ->
 4818                    domains_intersection(ZD, from_to(n(ZL), n(ZU)), ZD1),
 4819                    fd_put(Z, ZD1, ZPs)
 4820                ;   ZD1 = from_to(n(Z), n(Z))
 4821                ),
 4822                (   fd_get(X, XD, _), domain_infimum(XD, n(Min)) ->
 4823                    Z1 is Min rem Y,
 4824                    (   domain_contains(ZD1, Z1) -> true
 4825                    ;   neq_num(X, Min)
 4826                    )
 4827                ;   true
 4828                ),
 4829                (   fd_get(X, XD1, _), domain_supremum(XD1, n(Max)) ->
 4830                    Z2 is Max rem Y,
 4831                    (   domain_contains(ZD1, Z2) -> true
 4832                    ;   neq_num(X, Max)
 4833                    )
 4834                ;   true
 4835                )
 4836            ;   fd_get(X, XD1, XPs1),
 4837                % if possible, propagate at the boundaries
 4838                (   domain_infimum(XD1, n(Min)) ->
 4839                    (   Min rem Y =:= Z -> true
 4840                    ;   Y > 0, Min > 0 ->
 4841                        Next is ((Min - Z + Y - 1) div Y)*Y + Z,
 4842                        domain_remove_smaller_than(XD1, Next, XD2),
 4843                        fd_put(X, XD2, XPs1)
 4844                    ;   % TODO: bigger steps in other cases as well
 4845                        neq_num(X, Min)
 4846                    )
 4847                ;   true
 4848                ),
 4849                (   fd_get(X, XD3, XPs3) ->
 4850                    (   domain_supremum(XD3, n(Max)) ->
 4851                        (   Max rem Y =:= Z -> true
 4852                        ;   Y > 0, Max > 0  ->
 4853                            Prev is ((Max - Z) div Y)*Y + Z,
 4854                            domain_remove_greater_than(XD3, Prev, XD4),
 4855                            fd_put(X, XD4, XPs3)
 4856                        ;   % TODO: bigger steps in other cases as well
 4857                            neq_num(X, Max)
 4858                        )
 4859                    ;   true
 4860                    )
 4861                ;   true
 4862                )
 4863            )
 4864        ;   X == Y -> kill(MState), Z = 0
 4865        ;   fd_get(Z, ZD, ZPs) ->
 4866            fd_get(Y, _, YInf, YSup, _),
 4867            fd_get(X, _, XInf, XSup, _),
 4868            M cis max(abs(YInf),YSup),
 4869            (   XInf cis_geq n(0) -> Inf0 = n(0)
 4870            ;   Inf0 = XInf
 4871            ),
 4872            (   XSup cis_leq n(0) -> Sup0 = n(0)
 4873            ;   Sup0 = XSup
 4874            ),
 4875            NInf cis max(max(Inf0, -M + n(1)), min(XInf,-XSup)),
 4876            NSup cis min(min(Sup0, M - n(1)), max(abs(XInf),XSup)),
 4877            domains_intersection(ZD, from_to(NInf,NSup), ZD1),
 4878            fd_put(Z, ZD1, ZPs)
 4879        ;   true % TODO: propagate more
 4880        ).
 4881
 4882%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 4883% Z = max(X,Y)
 4884
 4885run_propagator(pmax(X,Y,Z), MState) :-
 4886        (   nonvar(X) ->
 4887            (   nonvar(Y) -> kill(MState), Z is max(X,Y)
 4888            ;   nonvar(Z) ->
 4889                (   Z =:= X -> kill(MState), X #>= Y
 4890                ;   Z > X -> Z = Y
 4891                ;   false % Z < X
 4892                )
 4893            ;   fd_get(Y, _, YInf, YSup, _),
 4894                (   YInf cis_gt n(X) -> Z = Y
 4895                ;   YSup cis_lt n(X) -> Z = X
 4896                ;   YSup = n(M) ->
 4897                    fd_get(Z, ZD, ZPs),
 4898                    domain_remove_greater_than(ZD, M, ZD1),
 4899                    fd_put(Z, ZD1, ZPs)
 4900                ;   true
 4901                )
 4902            )
 4903        ;   nonvar(Y) -> run_propagator(pmax(Y,X,Z), MState)
 4904        ;   fd_get(Z, ZD, ZPs) ->
 4905            fd_get(X, _, XInf, XSup, _),
 4906            fd_get(Y, _, YInf, YSup, _),
 4907            (   YInf cis_gt YSup -> kill(MState), Z = Y
 4908            ;   YSup cis_lt XInf -> kill(MState), Z = X
 4909            ;   n(M) cis max(XSup, YSup) ->
 4910                domain_remove_greater_than(ZD, M, ZD1),
 4911                fd_put(Z, ZD1, ZPs)
 4912            ;   true
 4913            )
 4914        ;   true
 4915        ).
 4916
 4917%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 4918% Z = min(X,Y)
 4919
 4920run_propagator(pmin(X,Y,Z), MState) :-
 4921        (   nonvar(X) ->
 4922            (   nonvar(Y) -> kill(MState), Z is min(X,Y)
 4923            ;   nonvar(Z) ->
 4924                (   Z =:= X -> kill(MState), X #=< Y
 4925                ;   Z < X -> Z = Y
 4926                ;   false % Z > X
 4927                )
 4928            ;   fd_get(Y, _, YInf, YSup, _),
 4929                (   YSup cis_lt n(X) -> Z = Y
 4930                ;   YInf cis_gt n(X) -> Z = X
 4931                ;   YInf = n(M) ->
 4932                    fd_get(Z, ZD, ZPs),
 4933                    domain_remove_smaller_than(ZD, M, ZD1),
 4934                    fd_put(Z, ZD1, ZPs)
 4935                ;   true
 4936                )
 4937            )
 4938        ;   nonvar(Y) -> run_propagator(pmin(Y,X,Z), MState)
 4939        ;   fd_get(Z, ZD, ZPs) ->
 4940            fd_get(X, _, XInf, XSup, _),
 4941            fd_get(Y, _, YInf, YSup, _),
 4942            (   YSup cis_lt YInf -> kill(MState), Z = Y
 4943            ;   YInf cis_gt XSup -> kill(MState), Z = X
 4944            ;   n(M) cis min(XInf, YInf) ->
 4945                domain_remove_smaller_than(ZD, M, ZD1),
 4946                fd_put(Z, ZD1, ZPs)
 4947            ;   true
 4948            )
 4949        ;   true
 4950        ).
 4951%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 4952% Z = X ^ Y
 4953
 4954run_propagator(pexp(X,Y,Z), MState) :-
 4955        (   X == 1 -> kill(MState), Z = 1
 4956        ;   X == 0 -> kill(MState), Z in 0..1, Z #<==> Y #= 0
 4957        ;   Y == 0 -> kill(MState), Z = 1
 4958        ;   Y == 1 -> kill(MState), Z = X
 4959        ;   nonvar(X) ->
 4960            (   nonvar(Y) ->
 4961                (   Y >= 0 -> true ; X =:= -1 ),
 4962                kill(MState),
 4963                Z is X^Y
 4964            ;   nonvar(Z) ->
 4965                (   Z > 1 ->
 4966                    kill(MState),
 4967                    integer_log_b(Z, X, 1, Y)
 4968                ;   true
 4969                )
 4970            ;   fd_get(Y, _, YL, YU, _),
 4971                fd_get(Z, ZD, ZPs),
 4972                (   X > 0, YL cis_geq n(0) ->
 4973                    NZL cis n(X)^YL,
 4974                    NZU cis n(X)^YU,
 4975                    domains_intersection(ZD, from_to(NZL,NZU), NZD),
 4976                    fd_put(Z, NZD, ZPs)
 4977                ;   true
 4978                ),
 4979                (   X > 0,
 4980                    fd_get(Z, _, _, n(ZMax), _),
 4981                    ZMax > 0 ->
 4982                    floor_integer_log_b(ZMax, X, 1, YCeil),
 4983                    Y in inf..YCeil
 4984                ;   true
 4985                )
 4986            )
 4987        ;   nonvar(Z) ->
 4988            (   nonvar(Y) ->
 4989                integer_kth_root(Z, Y, R),
 4990                kill(MState),
 4991                (   even(Y) ->
 4992                    N is -R,
 4993                    X in N \/ R
 4994                ;   X = R
 4995                )
 4996            ;   fd_get(X, _, n(NXL), _, _), NXL > 1 ->
 4997                (   Z > 1, between(NXL, Z, Exp), NXL^Exp > Z ->
 4998                    Exp1 is Exp - 1,
 4999                    fd_get(Y, YD, YPs),
 5000                    domains_intersection(YD, from_to(n(1),n(Exp1)), YD1),
 5001                    fd_put(Y, YD1, YPs),
 5002                    (   fd_get(X, XD, XPs) ->
 5003                        domain_infimum(YD1, n(YL)),
 5004                        integer_kth_root_leq(Z, YL, RU),
 5005                        domains_intersection(XD, from_to(n(NXL),n(RU)), XD1),
 5006                        fd_put(X, XD1, XPs)
 5007                    ;   true
 5008                    )
 5009                ;   true
 5010                )
 5011            ;   true
 5012            )
 5013        ;   nonvar(Y), Y > 0 ->
 5014            (   even(Y) ->
 5015                geq(Z, 0)
 5016            ;   true
 5017            ),
 5018            (   fd_get(X, XD, XL, XU, _), fd_get(Z, ZD, ZL, ZU, ZPs) ->
 5019                (   domain_contains(ZD, 0) -> XD1 = XD
 5020                ;   domain_remove(XD, 0, XD1)
 5021                ),
 5022                (   domain_contains(XD, 0) -> ZD1 = ZD
 5023                ;   domain_remove(ZD, 0, ZD1)
 5024                ),
 5025                (   even(Y) ->
 5026                    (   XL cis_geq n(0) ->
 5027                        NZL cis XL^n(Y)
 5028                    ;   XU cis_leq n(0) ->
 5029                        NZL cis XU^n(Y)
 5030                    ;   NZL = n(0)
 5031                    ),
 5032                    NZU cis max(abs(XL),abs(XU))^n(Y),
 5033                    domains_intersection(ZD1, from_to(NZL,NZU), ZD2)
 5034                ;   (   finite(XL) ->
 5035                        NZL cis XL^n(Y),
 5036                        NZU cis XU^n(Y),
 5037                        domains_intersection(ZD1, from_to(NZL,NZU), ZD2)
 5038                    ;   ZD2 = ZD1
 5039                    )
 5040                ),
 5041                fd_put(Z, ZD2, ZPs),
 5042                (   even(Y), ZU = n(Num) ->
 5043                    integer_kth_root_leq(Num, Y, RU),
 5044                    (   XL cis_geq n(0), ZL = n(Num1) ->
 5045                        integer_kth_root_leq(Num1, Y, RL0),
 5046                        (   RL0^Y < Num1 -> RL is RL0 + 1
 5047                        ;   RL = RL0
 5048                        )
 5049                    ;   RL is -RU
 5050                    ),
 5051                    RL =< RU,
 5052                    NXD = from_to(n(RL),n(RU))
 5053                ;   odd(Y), ZL cis_geq n(0), ZU = n(Num) ->
 5054                    integer_kth_root_leq(Num, Y, RU),
 5055                    ZL = n(Num1),
 5056                    integer_kth_root_leq(Num1, Y, RL0),
 5057                    (   RL0^Y < Num1 -> RL is RL0 + 1
 5058                    ;   RL = RL0
 5059                    ),
 5060                    RL =< RU,
 5061                    NXD = from_to(n(RL),n(RU))
 5062                ;   NXD = XD1   % TODO: propagate more
 5063                ),
 5064                (   fd_get(X, XD2, XPs) ->
 5065                    domains_intersection(XD2, XD1, XD3),
 5066                    domains_intersection(XD3, NXD, XD4),
 5067                    fd_put(X, XD4, XPs)
 5068                ;   true
 5069                )
 5070            ;   true
 5071            )
 5072        ;   fd_get(X, _, XL, _, _),
 5073            XL cis_gt n(0),
 5074            fd_get(Y, _, YL, _, _),
 5075            YL cis_gt n(0),
 5076            fd_get(Z, ZD, ZPs) ->
 5077            n(NZL) cis XL^YL,
 5078            domain_remove_smaller_than(ZD, NZL, ZD1),
 5079            fd_put(Z, ZD1, ZPs)
 5080        ;   true
 5081        ).
 5082
 5083%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 5084run_propagator(pzcompare(Order, A, B), MState) :-
 5085        (   A == B -> kill(MState), Order = (=)
 5086        ;   (   nonvar(A) ->
 5087                (   nonvar(B) ->
 5088                    kill(MState),
 5089                    (   A > B -> Order = (>)
 5090                    ;   Order = (<)
 5091                    )
 5092                ;   fd_get(B, _, BL, BU, _),
 5093                    (   BL cis_gt n(A) -> kill(MState), Order = (<)
 5094                    ;   BU cis_lt n(A) -> kill(MState), Order = (>)
 5095                    ;   true
 5096                    )
 5097                )
 5098            ;   nonvar(B) ->
 5099                fd_get(A, _, AL, AU, _),
 5100                (   AL cis_gt n(B) -> kill(MState), Order = (>)
 5101                ;   AU cis_lt n(B) -> kill(MState), Order = (<)
 5102                ;   true
 5103                )
 5104            ;   fd_get(A, _, AL, AU, _),
 5105                fd_get(B, _, BL, BU, _),
 5106                (   AL cis_gt BU -> kill(MState), Order = (>)
 5107                ;   AU cis_lt BL -> kill(MState), Order = (<)
 5108                ;   true
 5109                )
 5110            )
 5111        ).
 5112
 5113%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 5114
 5115% reified constraints
 5116
 5117%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 5118
 5119run_propagator(reified_in(V,Dom,B), MState) :-
 5120        (   integer(V) ->
 5121            kill(MState),
 5122            (   domain_contains(Dom, V) -> B = 1
 5123            ;   B = 0
 5124            )
 5125        ;   B == 1 -> kill(MState), domain(V, Dom)
 5126        ;   B == 0 -> kill(MState), domain_complement(Dom, C), domain(V, C)
 5127        ;   fd_get(V, VD, _),
 5128            (   domains_intersection(VD, Dom, I) ->
 5129                (   I == VD -> kill(MState), B = 1
 5130                ;   true
 5131                )
 5132            ;   kill(MState), B = 0
 5133            )
 5134        ).
 5135
 5136run_propagator(reified_tuple_in(Tuple, R, B), MState) :-
 5137        get_attr(R, clpfd_relation, Relation),
 5138        (   B == 1 -> kill(MState), tuples_in([Tuple], Relation)
 5139        ;   (   ground(Tuple) ->
 5140                kill(MState),
 5141                (   memberchk(Tuple, Relation) -> B = 1
 5142                ;   B = 0
 5143                )
 5144            ;   relation_unifiable(Relation, Tuple, Us, _, _),
 5145                (   Us = [] -> kill(MState), B = 0
 5146                ;   true
 5147                )
 5148            )
 5149        ).
 5150
 5151run_propagator(tuples_not_in(Tuples, Relation, B), MState) :-
 5152        (   B == 0 ->
 5153            kill(MState),
 5154            tuples_in_conjunction(Tuples, Relation, Conj),
 5155            #\ Conj
 5156        ;   true
 5157        ).
 5158
 5159run_propagator(kill_reified_tuples(B, Ps, Bs), _) :-
 5160        (   B == 0 ->
 5161            maplist(kill_entailed, Ps),
 5162            phrase(as(Bs), As),
 5163            maplist(kill_entailed, As)
 5164        ;   true
 5165        ).
 5166
 5167run_propagator(reified_fd(V,B), MState) :-
 5168        (   fd_inf(V, I), I \== inf, fd_sup(V, S), S \== sup ->
 5169            kill(MState),
 5170            B = 1
 5171        ;   B == 0 ->
 5172            (   fd_inf(V, inf) -> true
 5173            ;   fd_sup(V, sup) -> true
 5174            ;   false
 5175            )
 5176        ;   true
 5177        ).
 5178
 5179% The result of X/Y, X mod Y, and X rem Y is undefined iff Y is 0.
 5180
 5181run_propagator(pskeleton(X,Y,D,Skel,Z,_), MState) :-
 5182        (   Y == 0 -> kill(MState), D = 0
 5183        ;   D == 1 -> kill(MState), neq_num(Y, 0), skeleton([X,Y,Z], Skel)
 5184        ;   integer(Y), Y =\= 0 -> kill(MState), D = 1, skeleton([X,Y,Z], Skel)
 5185        ;   fd_get(Y, YD, _), \+ domain_contains(YD, 0) ->
 5186            kill(MState),
 5187            D = 1, skeleton([X,Y,Z], Skel)
 5188        ;   true
 5189        ).
 5190
 5191/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 5192   Propagators for arithmetic functions that only propagate
 5193   functionally. These are currently the bitwise operations.
 5194- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 5195
 5196run_propagator(pfunction(Op,A,B,R), MState) :-
 5197        (   integer(A), integer(B) ->
 5198            kill(MState),
 5199            Expr =.. [Op,A,B],
 5200            R is Expr
 5201        ;   true
 5202        ).
 5203run_propagator(pfunction(Op,A,R), MState) :-
 5204        (   integer(A) ->
 5205            kill(MState),
 5206            Expr =.. [Op,A],
 5207            R is Expr
 5208        ;   true
 5209        ).
 5210
 5211%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 5212
 5213run_propagator(reified_geq(DX,X,DY,Y,Ps,B), MState) :-
 5214        (   DX == 0 -> kill(MState, Ps), B = 0
 5215        ;   DY == 0 -> kill(MState, Ps), B = 0
 5216        ;   B == 1 -> kill(MState), DX = 1, DY = 1, geq(X, Y)
 5217        ;   DX == 1, DY == 1 ->
 5218            (   var(B) ->
 5219                (   nonvar(X) ->
 5220                    (   nonvar(Y) ->
 5221                        kill(MState),
 5222                        (   X >= Y -> B = 1 ; B = 0 )
 5223                    ;   fd_get(Y, _, YL, YU, _),
 5224                        (   n(X) cis_geq YU -> kill(MState, Ps), B = 1
 5225                        ;   n(X) cis_lt YL -> kill(MState, Ps), B = 0
 5226                        ;   true
 5227                        )
 5228                    )
 5229                ;   nonvar(Y) ->
 5230                    fd_get(X, _, XL, XU, _),
 5231                    (   XL cis_geq n(Y) -> kill(MState, Ps), B = 1
 5232                    ;   XU cis_lt n(Y) -> kill(MState, Ps), B = 0
 5233                    ;   true
 5234                    )
 5235                ;   X == Y -> kill(MState, Ps), B = 1
 5236                ;   fd_get(X, _, XL, XU, _),
 5237                    fd_get(Y, _, YL, YU, _),
 5238                    (   XL cis_geq YU -> kill(MState, Ps), B = 1
 5239                    ;   XU cis_lt YL -> kill(MState, Ps), B = 0
 5240                    ;   true
 5241                    )
 5242                )
 5243            ;   B =:= 0 -> kill(MState), X #< Y
 5244            ;   true
 5245            )
 5246        ;   true
 5247        ).
 5248
 5249%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
<