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

:- use_module(library(clpfd)).

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

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

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

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

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

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

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

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

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

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

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

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

is executed as if it were written as:

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

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

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

Example: Factorial relation

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

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

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

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

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

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

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

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

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

?- n_factorial(N, 3).
false.

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

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

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

Combinatorial constraints

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

Domains

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

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

Example: Sudoku

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

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

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

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

Sample query:

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

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

Residual goals

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Example: Eight queens puzzle

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Optimisation

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

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

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

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

Reification

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

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