1/*  Part of SWI-Prolog
    2
    3    Author:        Markus Triska
    4    E-mail:        triska@metalevel.at
    5    WWW:           http://www.swi-prolog.org
    6    Copyright (C): 2007-2017 Markus Triska
    7    All rights reserved.
    8
    9    Redistribution and use in source and binary forms, with or without
   10    modification, are permitted provided that the following conditions
   11    are met:
   12
   13    1. Redistributions of source code must retain the above copyright
   14       notice, this list of conditions and the following disclaimer.
   15
   16    2. Redistributions in binary form must reproduce the above copyright
   17       notice, this list of conditions and the following disclaimer in
   18       the documentation and/or other materials provided with the
   19       distribution.
   20
   21    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
   22    "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
   23    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
   24    FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
   25    COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
   26    INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
   27    BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
   28    LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
   29    CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
   30    LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
   31    ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
   32    POSSIBILITY OF SUCH DAMAGE.
   33*/
   34
   35/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
   36   Thanks to Tom Schrijvers for his "bounds.pl", the first finite
   37   domain constraint solver included with SWI-Prolog. I've learned a
   38   lot from it and could even use some of the code for this solver.
   39   The propagation queue idea is taken from "prop.pl", a prototype
   40   solver also written by Tom. Highlights of the present solver:
   41
   42   Symbolic constants for infinities
   43   ---------------------------------
   44
   45   ?- X #>= 0, Y #=< 0.
   46   %@ X in 0..sup,
   47   %@ Y in inf..0.
   48
   49   No artificial limits (using GMP)
   50   ---------------------------------
   51
   52   ?- N #= 2^66, X #\= N.
   53   %@ N = 73786976294838206464,
   54   %@ X in inf..73786976294838206463\/73786976294838206465..sup.
   55
   56   Often stronger propagation
   57   ---------------------------------
   58
   59   ?- Y #= abs(X), Y #\= 3, Z * Z #= 4.
   60   %@ Y in 0..2\/4..sup,
   61   %@ Y#=abs(X),
   62   %@ X in inf.. -4\/ -2..2\/4..sup,
   63   %@ Z in -2\/2.
   64
   65   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
   66
   67   Development of this library has moved to SICStus Prolog. If you
   68   need any additional features or want to help, please file an issue at:
   69
   70                    https://github.com/triska/clpz
   71                    ==============================
   72
   73- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
   74
   75:- module(clpfd, [
   76                  op(760, yfx, #<==>),
   77                  op(750, xfy, #==>),
   78                  op(750, yfx, #<==),
   79                  op(740, yfx, #\/),
   80                  op(730, yfx, #\),
   81                  op(720, yfx, #/\),
   82                  op(710,  fy, #\),
   83                  op(700, xfx, #>),
   84                  op(700, xfx, #<),
   85                  op(700, xfx, #>=),
   86                  op(700, xfx, #=<),
   87                  op(700, xfx, #=),
   88                  op(700, xfx, #\=),
   89                  op(700, xfx, in),
   90                  op(700, xfx, ins),
   91                  op(450, xfx, ..), % should bind more tightly than \/
   92                  (#>)/2,
   93                  (#<)/2,
   94                  (#>=)/2,
   95                  (#=<)/2,
   96                  (#=)/2,
   97                  (#\=)/2,
   98                  (#\)/1,
   99                  (#<==>)/2,
  100                  (#==>)/2,
  101                  (#<==)/2,
  102                  (#\/)/2,
  103                  (#\)/2,
  104                  (#/\)/2,
  105                  (in)/2,
  106                  (ins)/2,
  107                  all_different/1,
  108                  all_distinct/1,
  109                  sum/3,
  110                  scalar_product/4,
  111                  tuples_in/2,
  112                  labeling/2,
  113                  label/1,
  114                  indomain/1,
  115                  lex_chain/1,
  116                  serialized/2,
  117                  global_cardinality/2,
  118                  global_cardinality/3,
  119                  circuit/1,
  120                  cumulative/1,
  121                  cumulative/2,
  122                  disjoint2/1,
  123                  element/3,
  124                  automaton/3,
  125                  automaton/8,
  126                  transpose/2,
  127                  zcompare/3,
  128                  chain/2,
  129                  fd_var/1,
  130                  fd_inf/2,
  131                  fd_sup/2,
  132                  fd_size/2,
  133                  fd_dom/2
  134                 ]).  135
  136:- public                               % called from goal_expansion
  137        clpfd_equal/2,
  138        clpfd_geq/2.  139
  140:- use_module(library(apply)).  141:- use_module(library(apply_macros)).  142:- use_module(library(assoc)).  143:- use_module(library(error)).  144:- use_module(library(lists)).  145:- use_module(library(pairs)).  146
  147:- op(700, xfx, cis).  148:- op(700, xfx, cis_geq).  149:- op(700, xfx, cis_gt).  150:- op(700, xfx, cis_leq).  151:- op(700, xfx, cis_lt).  152
  153/** <module> CLP(FD): Constraint Logic Programming over Finite Domains
  154
  155**Development of this library has moved to SICStus Prolog.**
  156
  157Please see [**CLP(Z)**](https://github.com/triska/clpz) for more
  158information.
  159
  160## Introduction			{#clpfd-intro}
  161
  162This library provides CLP(FD): Constraint Logic Programming over
  163Finite Domains. This is an instance of the general [CLP(_X_)
  164scheme](<#clp>), extending logic programming with reasoning over
  165specialised domains.
  166
  167CLP(FD) lets us reason about **integers** in a way that honors the
  168relational nature of Prolog. Read [**The Power of
  169Prolog**](https://www.metalevel.at/prolog) to understand how this
  170library is meant to be used in practice.
  171
  172There are two major use cases of CLP(FD) constraints:
  173
  174    1. [**declarative integer arithmetic**](<#clpfd-integer-arith>)
  175    2. solving **combinatorial problems** such as planning, scheduling
  176       and allocation tasks.
  177
  178The predicates of this library can be classified as:
  179
  180    * _arithmetic_ constraints like #=/2, #>/2 and #\=/2 [](<#clpfd-arithmetic>)
  181    * the _membership_ constraints in/2 and ins/2 [](<#clpfd-membership>)
  182    * the _enumeration_ predicates indomain/1, label/1 and labeling/2 [](<#clpfd-enumeration>)
  183    * _combinatorial_ constraints like all_distinct/1 and global_cardinality/2 [](<#clpfd-global>)
  184    * _reification_ predicates such as #<==>/2 [](<#clpfd-reification-predicates>)
  185    * _reflection_ predicates such as fd_dom/2 [](<#clpfd-reflection-predicates>)
  186
  187In most cases, [_arithmetic constraints_](<#clpfd-arith-constraints>)
  188are the only predicates you will ever need from this library. When
  189reasoning over integers, simply replace low-level arithmetic
  190predicates like `(is)/2` and `(>)/2` by the corresponding CLP(FD)
  191constraints like #=/2 and #>/2 to honor and preserve declarative
  192properties of your programs. For satisfactory performance, arithmetic
  193constraints are implicitly rewritten at compilation time so that
  194low-level fallback predicates are automatically used whenever
  195possible.
  196
  197Almost all Prolog programs also reason about integers. Therefore, it
  198is highly advisable that you make CLP(FD) constraints available in all
  199your programs. One way to do this is to put the following directive in
  200your =|~/.swiplrc|= initialisation file:
  201
  202==
  203:- use_module(library(clpfd)).
  204==
  205
  206All example programs that appear in the CLP(FD) documentation assume
  207that you have done this.
  208
  209Important concepts and principles of this library are illustrated by
  210means of usage examples that are available in a public git repository:
  211[**github.com/triska/clpfd**](https://github.com/triska/clpfd)
  212
  213If you are used to the complicated operational considerations that
  214low-level arithmetic primitives necessitate, then moving to CLP(FD)
  215constraints may, due to their power and convenience, at first feel to
  216you excessive and almost like cheating. It _isn't_. Constraints are an
  217integral part of all popular Prolog systems, and they are designed
  218to help you eliminate and avoid the use of low-level and less general
  219primitives by providing declarative alternatives that are meant to be
  220used instead.
  221
  222When teaching Prolog, CLP(FD) constraints should be introduced
  223_before_ explaining low-level arithmetic predicates and their
  224procedural idiosyncrasies. This is because constraints are easy to
  225explain, understand and use due to their purely relational nature. In
  226contrast, the modedness and directionality of low-level arithmetic
  227primitives are impure limitations that are better deferred to more
  228advanced lectures.
  229
  230We recommend the following reference (PDF:
  231[metalevel.at/swiclpfd.pdf](https://www.metalevel.at/swiclpfd.pdf)) for
  232citing this library in scientific publications:
  233
  234==
  235@inproceedings{Triska12,
  236  author    = {Markus Triska},
  237  title     = {The Finite Domain Constraint Solver of {SWI-Prolog}},
  238  booktitle = {FLOPS},
  239  series    = {LNCS},
  240  volume    = {7294},
  241  year      = {2012},
  242  pages     = {307-316}
  243}
  244==
  245
  246More information about CLP(FD) constraints and their implementation is
  247contained in: [**metalevel.at/drt.pdf**](https://www.metalevel.at/drt.pdf)
  248
  249The best way to discuss applying, improving and extending CLP(FD)
  250constraints is to use the dedicated `clpfd` tag on
  251[stackoverflow.com](http://stackoverflow.com). Several of the world's
  252foremost CLP(FD) experts regularly participate in these discussions
  253and will help you for free on this platform.
  254
  255## Arithmetic constraints		{#clpfd-arith-constraints}
  256
  257In modern Prolog systems, *arithmetic constraints* subsume and
  258supersede low-level predicates over integers. The main advantage of
  259arithmetic constraints is that they are true _relations_ and can be
  260used in all directions. For most programs, arithmetic constraints are
  261the only predicates you will ever need from this library.
  262
  263The most important arithmetic constraint is #=/2, which subsumes both
  264`(is)/2` and `(=:=)/2` over integers. Use #=/2 to make your programs
  265more general.
  266
  267In total, the arithmetic constraints are:
  268
  269    | Expr1 `#=`  Expr2  | Expr1 equals Expr2                       |
  270    | Expr1 `#\=` Expr2  | Expr1 is not equal to Expr2              |
  271    | Expr1 `#>=` Expr2  | Expr1 is greater than or equal to Expr2  |
  272    | Expr1 `#=<` Expr2  | Expr1 is less than or equal to Expr2     |
  273    | Expr1 `#>` Expr2   | Expr1 is greater than Expr2              |
  274    | Expr1 `#<` Expr2   | Expr1 is less than Expr2                 |
  275
  276`Expr1` and `Expr2` denote *arithmetic expressions*, which are:
  277
  278    | _integer_          | Given value                          |
  279    | _variable_         | Unknown integer                      |
  280    | ?(_variable_)      | Unknown integer                      |
  281    | -Expr              | Unary minus                          |
  282    | Expr + Expr        | Addition                             |
  283    | Expr * Expr        | Multiplication                       |
  284    | Expr - Expr        | Subtraction                          |
  285    | Expr ^ Expr        | Exponentiation                       |
  286    | min(Expr,Expr)     | Minimum of two expressions           |
  287    | max(Expr,Expr)     | Maximum of two expressions           |
  288    | Expr `mod` Expr    | Modulo induced by floored division   |
  289    | Expr `rem` Expr    | Modulo induced by truncated division |
  290    | abs(Expr)          | Absolute value                       |
  291    | Expr // Expr       | Truncated integer division           |
  292    | Expr div Expr      | Floored integer division             |
  293
  294where `Expr` again denotes an arithmetic expression.
  295
  296The bitwise operations `(\)/1`, `(/\)/2`, `(\/)/2`, `(>>)/2`,
  297`(<<)/2`, `lsb/1`, `msb/1`, `popcount/1` and `(xor)/2` are also
  298supported.
  299
  300## Declarative integer arithmetic		{#clpfd-integer-arith}
  301
  302The [_arithmetic constraints_](<#clpfd-arith-constraints>) #=/2, #>/2
  303etc. are meant to be used _instead_ of the primitives `(is)/2`,
  304`(=:=)/2`, `(>)/2` etc. over integers. Almost all Prolog programs also
  305reason about integers. Therefore, it is recommended that you put the
  306following directive in your =|~/.swiplrc|= initialisation file to make
  307CLP(FD) constraints available in all your programs:
  308
  309==
  310:- use_module(library(clpfd)).
  311==
  312
  313Throughout the following, it is assumed that you have done this.
  314
  315The most basic use of CLP(FD) constraints is _evaluation_ of
  316arithmetic expressions involving integers. For example:
  317
  318==
  319?- X #= 1+2.
  320X = 3.
  321==
  322
  323This could in principle also be achieved with the lower-level
  324predicate `(is)/2`. However, an important advantage of arithmetic
  325constraints is their purely relational nature: Constraints can be used
  326in _all directions_, also if one or more of their arguments are only
  327partially instantiated. For example:
  328
  329==
  330?- 3 #= Y+2.
  331Y = 1.
  332==
  333
  334This relational nature makes CLP(FD) constraints easy to explain and
  335use, and well suited for beginners and experienced Prolog programmers
  336alike. In contrast, when using low-level integer arithmetic, we get:
  337
  338==
  339?- 3 is Y+2.
  340ERROR: is/2: Arguments are not sufficiently instantiated
  341
  342?- 3 =:= Y+2.
  343ERROR: =:=/2: Arguments are not sufficiently instantiated
  344==
  345
  346Due to the necessary operational considerations, the use of these
  347low-level arithmetic predicates is considerably harder to understand
  348and should therefore be deferred to more advanced lectures.
  349
  350For supported expressions, CLP(FD) constraints are drop-in
  351replacements of these low-level arithmetic predicates, often yielding
  352more general programs. See [`n_factorial/2`](<#clpfd-factorial>) for an
  353example.
  354
  355This library uses goal_expansion/2 to automatically rewrite
  356constraints at compilation time so that low-level arithmetic
  357predicates are _automatically_ used whenever possible. For example,
  358the predicate:
  359
  360==
  361positive_integer(N) :- N #>= 1.
  362==
  363
  364is executed as if it were written as:
  365
  366==
  367positive_integer(N) :-
  368        (   integer(N)
  369        ->  N >= 1
  370        ;   N #>= 1
  371        ).
  372==
  373
  374This illustrates why the performance of CLP(FD) constraints is almost
  375always completely satisfactory when they are used in modes that can be
  376handled by low-level arithmetic. To disable the automatic rewriting,
  377set the Prolog flag `clpfd_goal_expansion` to `false`.
  378
  379If you are used to the complicated operational considerations that
  380low-level arithmetic primitives necessitate, then moving to CLP(FD)
  381constraints may, due to their power and convenience, at first feel to
  382you excessive and almost like cheating. It _isn't_. Constraints are an
  383integral part of all popular Prolog systems, and they are designed
  384to help you eliminate and avoid the use of low-level and less general
  385primitives by providing declarative alternatives that are meant to be
  386used instead.
  387
  388
  389## Example: Factorial relation {#clpfd-factorial}
  390
  391We illustrate the benefit of using #=/2 for more generality with a
  392simple example.
  393
  394Consider first a rather conventional definition of `n_factorial/2`,
  395relating each natural number _N_ to its factorial _F_:
  396
  397==
  398n_factorial(0, 1).
  399n_factorial(N, F) :-
  400        N #> 0,
  401        N1 #= N - 1,
  402        n_factorial(N1, F1),
  403        F #= N * F1.
  404==
  405
  406This program uses CLP(FD) constraints _instead_ of low-level
  407arithmetic throughout, and everything that _would have worked_ with
  408low-level arithmetic _also_ works with CLP(FD) constraints, retaining
  409roughly the same performance. For example:
  410
  411==
  412?- n_factorial(47, F).
  413F = 258623241511168180642964355153611979969197632389120000000000 ;
  414false.
  415==
  416
  417Now the point: Due to the increased flexibility and generality of
  418CLP(FD) constraints, we are free to _reorder_ the goals as follows:
  419
  420==
  421n_factorial(0, 1).
  422n_factorial(N, F) :-
  423        N #> 0,
  424        N1 #= N - 1,
  425        F #= N * F1,
  426        n_factorial(N1, F1).
  427==
  428
  429In this concrete case, _termination_ properties of the predicate are
  430improved. For example, the following queries now both terminate:
  431
  432==
  433?- n_factorial(N, 1).
  434N = 0 ;
  435N = 1 ;
  436false.
  437
  438?- n_factorial(N, 3).
  439false.
  440==
  441
  442To make the predicate terminate if _any_ argument is instantiated, add
  443the (implied) constraint `F #\= 0` before the recursive call.
  444Otherwise, the query `n_factorial(N, 0)` is the only non-terminating
  445case of this kind.
  446
  447The value of CLP(FD) constraints does _not_ lie in completely freeing
  448us from _all_ procedural phenomena. For example, the two programs do
  449not even have the same _termination properties_ in all cases.
  450Instead, the primary benefit of CLP(FD) constraints is that they allow
  451you to try different execution orders and apply [**declarative
  452debugging**](https://www.metalevel.at/prolog/debugging)
  453techniques _at all_!  Reordering goals (and clauses) can significantly
  454impact the performance of Prolog programs, and you are free to try
  455different variants if you use declarative approaches. Moreover, since
  456all CLP(FD) constraints _always terminate_, placing them earlier can
  457at most _improve_, never worsen, the termination properties of your
  458programs. An additional benefit of CLP(FD) constraints is that they
  459eliminate the complexity of introducing `(is)/2` and `(=:=)/2` to
  460beginners, since _both_ predicates are subsumed by #=/2 when reasoning
  461over integers.
  462
  463In the case above, the clauses are mutually exclusive _if_ the first
  464argument is sufficiently instantiated. To make the predicate
  465deterministic in such cases while retaining its generality, you can
  466use zcompare/3 to _reify_ a comparison, making the different cases
  467distinguishable by pattern matching. For example, in this concrete
  468case and others like it, you can use `zcompare(Comp, 0, N)` to obtain as
  469`Comp` the symbolic outcome (`<`, `=`, `>`) of 0 compared to N.
  470
  471## Combinatorial constraints  {#clpfd-combinatorial}
  472
  473In addition to subsuming and replacing low-level arithmetic
  474predicates, CLP(FD) constraints are often used to solve combinatorial
  475problems such as planning, scheduling and allocation tasks. Among the
  476most frequently used *combinatorial constraints* are all_distinct/1,
  477global_cardinality/2 and cumulative/2. This library also provides
  478several other constraints like disjoint2/1 and automaton/8, which are
  479useful in more specialized applications.
  480
  481## Domains                             {#clpfd-domains}
  482
  483Each CLP(FD) variable has an associated set of admissible integers,
  484which we call the variable's *domain*. Initially, the domain of each
  485CLP(FD) variable is the set of _all_ integers. CLP(FD) constraints
  486like #=/2, #>/2 and #\=/2 can at most reduce, and never extend, the
  487domains of their arguments. The constraints in/2 and ins/2 let us
  488explicitly state domains of CLP(FD) variables. The process of
  489determining and adjusting domains of variables is called constraint
  490*propagation*, and it is performed automatically by this library. When
  491the domain of a variable contains only one element, then the variable
  492is automatically unified to that element.
  493
  494Domains are taken into account when further constraints are stated,
  495and by enumeration predicates like labeling/2.
  496
  497## Example: Sudoku {#clpfd-sudoku}
  498
  499As another example, consider _Sudoku_: It is a popular puzzle
  500over integers that can be easily solved with CLP(FD) constraints.
  501
  502==
  503sudoku(Rows) :-
  504        length(Rows, 9), maplist(same_length(Rows), Rows),
  505        append(Rows, Vs), Vs ins 1..9,
  506        maplist(all_distinct, Rows),
  507        transpose(Rows, Columns),
  508        maplist(all_distinct, Columns),
  509        Rows = [As,Bs,Cs,Ds,Es,Fs,Gs,Hs,Is],
  510        blocks(As, Bs, Cs),
  511        blocks(Ds, Es, Fs),
  512        blocks(Gs, Hs, Is).
  513
  514blocks([], [], []).
  515blocks([N1,N2,N3|Ns1], [N4,N5,N6|Ns2], [N7,N8,N9|Ns3]) :-
  516        all_distinct([N1,N2,N3,N4,N5,N6,N7,N8,N9]),
  517        blocks(Ns1, Ns2, Ns3).
  518
  519problem(1, [[_,_,_,_,_,_,_,_,_],
  520            [_,_,_,_,_,3,_,8,5],
  521            [_,_,1,_,2,_,_,_,_],
  522            [_,_,_,5,_,7,_,_,_],
  523            [_,_,4,_,_,_,1,_,_],
  524            [_,9,_,_,_,_,_,_,_],
  525            [5,_,_,_,_,_,_,7,3],
  526            [_,_,2,_,1,_,_,_,_],
  527            [_,_,_,_,4,_,_,_,9]]).
  528==
  529
  530Sample query:
  531
  532==
  533?- problem(1, Rows), sudoku(Rows), maplist(writeln, Rows).
  534[9,8,7,6,5,4,3,2,1]
  535[2,4,6,1,7,3,9,8,5]
  536[3,5,1,9,2,8,7,4,6]
  537[1,2,8,5,3,7,6,9,4]
  538[6,3,4,8,9,2,1,5,7]
  539[7,9,5,4,6,1,8,3,2]
  540[5,1,9,2,8,6,4,7,3]
  541[4,7,2,3,1,9,5,6,8]
  542[8,6,3,7,4,5,2,1,9]
  543Rows = [[9, 8, 7, 6, 5, 4, 3, 2|...], ... , [...|...]].
  544==
  545
  546In this concrete case, the constraint solver is strong enough to find
  547the unique solution without any search.
  548
  549
  550## Residual goals				{#clpfd-residual-goals}
  551
  552Here is an example session with a few queries and their answers:
  553
  554==
  555?- X #> 3.
  556X in 4..sup.
  557
  558?- X #\= 20.
  559X in inf..19\/21..sup.
  560
  561?- 2*X #= 10.
  562X = 5.
  563
  564?- X*X #= 144.
  565X in -12\/12.
  566
  567?- 4*X + 2*Y #= 24, X + Y #= 9, [X,Y] ins 0..sup.
  568X = 3,
  569Y = 6.
  570
  571?- X #= Y #<==> B, X in 0..3, Y in 4..5.
  572B = 0,
  573X in 0..3,
  574Y in 4..5.
  575==
  576
  577The answers emitted by the toplevel are called _residual programs_,
  578and the goals that comprise each answer are called **residual goals**.
  579In each case above, and as for all pure programs, the residual program
  580is declaratively equivalent to the original query. From the residual
  581goals, it is clear that the constraint solver has deduced additional
  582domain restrictions in many cases.
  583
  584To inspect residual goals, it is best to let the toplevel display them
  585for us. Wrap the call of your predicate into call_residue_vars/2 to
  586make sure that all constrained variables are displayed. To make the
  587constraints a variable is involved in available as a Prolog term for
  588further reasoning within your program, use copy_term/3. For example:
  589
  590==
  591?- X #= Y + Z, X in 0..5, copy_term([X,Y,Z], [X,Y,Z], Gs).
  592Gs = [clpfd: (X in 0..5), clpfd: (Y+Z#=X)],
  593X in 0..5,
  594Y+Z#=X.
  595==
  596
  597This library also provides _reflection_ predicates (like fd_dom/2,
  598fd_size/2 etc.) with which we can inspect a variable's current
  599domain. These predicates can be useful if you want to implement your
  600own labeling strategies.
  601
  602## Core relations and search    {#clpfd-search}
  603
  604Using CLP(FD) constraints to solve combinatorial tasks typically
  605consists of two phases:
  606
  607    1. First, all relevant constraints are stated.
  608    2. Second, if the domain of each involved variable is _finite_,
  609       then _enumeration predicates_ can be used to search for
  610       concrete solutions.
  611
  612It is good practice to keep the modeling part, via a dedicated
  613predicate called the *core relation*, separate from the actual
  614search for solutions. This lets us observe termination and
  615determinism properties of the core relation in isolation from the
  616search, and more easily try different search strategies.
  617
  618As an example of a constraint satisfaction problem, consider the
  619cryptoarithmetic puzzle SEND + MORE = MONEY, where different letters
  620denote distinct integers between 0 and 9. It can be modeled in CLP(FD)
  621as follows:
  622
  623==
  624puzzle([S,E,N,D] + [M,O,R,E] = [M,O,N,E,Y]) :-
  625        Vars = [S,E,N,D,M,O,R,Y],
  626        Vars ins 0..9,
  627        all_different(Vars),
  628                  S*1000 + E*100 + N*10 + D +
  629                  M*1000 + O*100 + R*10 + E #=
  630        M*10000 + O*1000 + N*100 + E*10 + Y,
  631        M #\= 0, S #\= 0.
  632==
  633
  634Notice that we are _not_ using labeling/2 in this predicate, so that
  635we can first execute and observe the modeling part in isolation.
  636Sample query and its result (actual variables replaced for
  637readability):
  638
  639==
  640?- puzzle(As+Bs=Cs).
  641As = [9, A2, A3, A4],
  642Bs = [1, 0, B3, A2],
  643Cs = [1, 0, A3, A2, C5],
  644A2 in 4..7,
  645all_different([9, A2, A3, A4, 1, 0, B3, C5]),
  64691*A2+A4+10*B3#=90*A3+C5,
  647A3 in 5..8,
  648A4 in 2..8,
  649B3 in 2..8,
  650C5 in 2..8.
  651==
  652
  653From this answer, we see that this core relation _terminates_ and is in
  654fact _deterministic_. Moreover, we see from the residual goals that
  655the constraint solver has deduced more stringent bounds for all
  656variables. Such observations are only possible if modeling and search
  657parts are cleanly separated.
  658
  659Labeling can then be used to search for solutions in a separate
  660predicate or goal:
  661
  662==
  663?- puzzle(As+Bs=Cs), label(As).
  664As = [9, 5, 6, 7],
  665Bs = [1, 0, 8, 5],
  666Cs = [1, 0, 6, 5, 2] ;
  667false.
  668==
  669
  670In this case, it suffices to label a subset of variables to find the
  671puzzle's unique solution, since the constraint solver is strong enough
  672to reduce the domains of remaining variables to singleton sets. In
  673general though, it is necessary to label all variables to obtain
  674ground solutions.
  675
  676## Example: Eight queens puzzle {#clpfd-n-queens}
  677
  678We illustrate the concepts of the preceding sections by means of the
  679so-called _eight queens puzzle_. The task is to place 8 queens on an
  6808x8 chessboard such that none of the queens is under attack. This
  681means that no two queens share the same row, column or diagonal.
  682
  683To express this puzzle via CLP(FD) constraints, we must first pick a
  684suitable representation. Since CLP(FD) constraints reason over
  685_integers_, we must find a way to map the positions of queens to
  686integers. Several such mappings are conceivable, and it is not
  687immediately obvious which we should use. On top of that, different
  688constraints can be used to express the desired relations. For such
  689reasons, _modeling_ combinatorial problems via CLP(FD) constraints
  690often necessitates some creativity and has been described as more of
  691an art than a science.
  692
  693In our concrete case, we observe that there must be exactly one queen
  694per column. The following representation therefore suggests itself: We
  695are looking for 8 integers, one for each column, where each integer
  696denotes the _row_ of the queen that is placed in the respective
  697column, and which are subject to certain constraints.
  698
  699In fact, let us now generalize the task to the so-called _N queens
  700puzzle_, which is obtained by replacing 8 by _N_ everywhere it occurs
  701in the above description. We implement the above considerations in the
  702**core relation** `n_queens/2`, where the first argument is the number
  703of queens (which is identical to the number of rows and columns of the
  704generalized chessboard), and the second argument is a list of _N_
  705integers that represents a solution in the form described above.
  706
  707==
  708n_queens(N, Qs) :-
  709        length(Qs, N),
  710        Qs ins 1..N,
  711        safe_queens(Qs).
  712
  713safe_queens([]).
  714safe_queens([Q|Qs]) :- safe_queens(Qs, Q, 1), safe_queens(Qs).
  715
  716safe_queens([], _, _).
  717safe_queens([Q|Qs], Q0, D0) :-
  718        Q0 #\= Q,
  719        abs(Q0 - Q) #\= D0,
  720        D1 #= D0 + 1,
  721        safe_queens(Qs, Q0, D1).
  722==
  723
  724Note that all these predicates can be used in _all directions_: We
  725can use them to _find_ solutions, _test_ solutions and _complete_
  726partially instantiated solutions.
  727
  728The original task can be readily solved with the following query:
  729
  730==
  731?- n_queens(8, Qs), label(Qs).
  732Qs = [1, 5, 8, 6, 3, 7, 2, 4] .
  733==
  734
  735Using suitable labeling strategies, we can easily find solutions with
  73680 queens and more:
  737
  738==
  739?- n_queens(80, Qs), labeling([ff], Qs).
  740Qs = [1, 3, 5, 44, 42, 4, 50, 7, 68|...] .
  741
  742?- time((n_queens(90, Qs), labeling([ff], Qs))).
  743% 5,904,401 inferences, 0.722 CPU in 0.737 seconds (98% CPU)
  744Qs = [1, 3, 5, 50, 42, 4, 49, 7, 59|...] .
  745==
  746
  747Experimenting with different search strategies is easy because we have
  748separated the core relation from the actual search.
  749
  750
  751
  752## Optimisation    {#clpfd-optimisation}
  753
  754We can use labeling/2 to minimize or maximize the value of a CLP(FD)
  755expression, and generate solutions in increasing or decreasing order
  756of the value. See the labeling options `min(Expr)` and `max(Expr)`,
  757respectively.
  758
  759Again, to easily try different labeling options in connection with
  760optimisation, we recommend to introduce a dedicated predicate for
  761posting constraints, and to use `labeling/2` in a separate goal. This
  762way, we can observe properties of the core relation in isolation,
  763and try different labeling options without recompiling our code.
  764
  765If necessary, we can use `once/1` to commit to the first optimal
  766solution. However, it is often very valuable to see alternative
  767solutions that are _also_ optimal, so that we can choose among optimal
  768solutions by other criteria. For the sake of
  769[**purity**](https://www.metalevel.at/prolog/purity) and
  770completeness, we recommend to avoid `once/1` and other constructs that
  771lead to impurities in CLP(FD) programs.
  772
  773Related to optimisation with CLP(FD) constraints are
  774[`library(simplex)`](http://eu.swi-prolog.org/man/simplex.html) and
  775CLP(Q) which reason about _linear_ constraints over rational numbers.
  776
  777## Reification				{#clpfd-reification}
  778
  779The constraints in/2, #=/2, #\=/2, #</2, #>/2, #=</2, and #>=/2 can be
  780_reified_, which means reflecting their truth values into Boolean
  781values represented by the integers 0 and 1. Let P and Q denote
  782reifiable constraints or Boolean variables, then:
  783
  784    | #\ Q      | True iff Q is false                  |
  785    | P #\/ Q   | True iff either P or Q               |
  786    | P #/\ Q   | True iff both P and Q                |
  787    | P #\ Q    | True iff either P or Q, but not both |
  788    | P #<==> Q | True iff P and Q are equivalent      |
  789    | P #==> Q  | True iff P implies Q                 |
  790    | P #<== Q  | True iff Q implies P                 |
  791
  792The constraints of this table are reifiable as well.
  793
  794When reasoning over Boolean variables, also consider using
  795CLP(B) constraints as provided by
  796[`library(clpb)`](http://eu.swi-prolog.org/man/clpb.html).
  797
  798## Enabling monotonic CLP(FD)		{#clpfd-monotonicity}
  799
  800In the default execution mode, CLP(FD) constraints still exhibit some
  801non-relational properties. For example, _adding_ constraints can yield
  802new solutions:
  803
  804==
  805?-          X #= 2, X = 1+1.
  806false.
  807
  808?- X = 1+1, X #= 2, X = 1+1.
  809X = 1+1.
  810==
  811
  812This behaviour is highly problematic from a logical point of view, and
  813it may render declarative debugging techniques inapplicable.
  814
  815Set the Prolog flag `clpfd_monotonic` to `true` to make CLP(FD)
  816**monotonic**: This means that _adding_ new constraints _cannot_ yield
  817new solutions. When this flag is `true`, we must wrap variables that
  818occur in arithmetic expressions with the functor `(?)/1` or `(#)/1`. For
  819example:
  820
  821==
  822?- set_prolog_flag(clpfd_monotonic, true).
  823true.
  824
  825?- #(X) #= #(Y) + #(Z).
  826#(Y)+ #(Z)#= #(X).
  827
  828?-          X #= 2, X = 1+1.
  829ERROR: Arguments are not sufficiently instantiated
  830==
  831
  832The wrapper can be omitted for variables that are already constrained
  833to integers.
  834
  835## Custom constraints			{#clpfd-custom-constraints}
  836
  837We can define custom constraints. The mechanism to do this is not yet
  838finalised, and we welcome suggestions and descriptions of use cases
  839that are important to you.
  840
  841As an example of how it can be done currently, let us define a new
  842custom constraint `oneground(X,Y,Z)`, where Z shall be 1 if at least
  843one of X and Y is instantiated:
  844
  845==
  846:- multifile clpfd:run_propagator/2.
  847
  848oneground(X, Y, Z) :-
  849        clpfd:make_propagator(oneground(X, Y, Z), Prop),
  850        clpfd:init_propagator(X, Prop),
  851        clpfd:init_propagator(Y, Prop),
  852        clpfd:trigger_once(Prop).
  853
  854clpfd:run_propagator(oneground(X, Y, Z), MState) :-
  855        (   integer(X) -> clpfd:kill(MState), Z = 1
  856        ;   integer(Y) -> clpfd:kill(MState), Z = 1
  857        ;   true
  858        ).
  859==
  860
  861First, clpfd:make_propagator/2 is used to transform a user-defined
  862representation of the new constraint to an internal form. With
  863clpfd:init_propagator/2, this internal form is then attached to X and
  864Y. From now on, the propagator will be invoked whenever the domains of
  865X or Y are changed. Then, clpfd:trigger_once/1 is used to give the
  866propagator its first chance for propagation even though the variables'
  867domains have not yet changed. Finally, clpfd:run_propagator/2 is
  868extended to define the actual propagator. As explained, this predicate
  869is automatically called by the constraint solver. The first argument
  870is the user-defined representation of the constraint as used in
  871clpfd:make_propagator/2, and the second argument is a mutable state
  872that can be used to prevent further invocations of the propagator when
  873the constraint has become entailed, by using clpfd:kill/1. An example
  874of using the new constraint:
  875
  876==
  877?- oneground(X, Y, Z), Y = 5.
  878Y = 5,
  879Z = 1,
  880X in inf..sup.
  881==
  882
  883## Applications   {#clpfd-applications}
  884
  885CLP(FD) applications that we find particularly impressive and worth
  886studying include:
  887
  888  * Michael Hendricks uses CLP(FD) constraints for flexible reasoning
  889    about _dates_ and _times_ in the
  890    [`julian`](http://www.swi-prolog.org/pack/list?p=julian) package.
  891  * Julien Cumin uses CLP(FD) constraints for integer arithmetic in
  892    [=Brachylog=](https://github.com/JCumin/Brachylog).
  893
  894## Acknowledgments {#clpfd-acknowledgments}
  895
  896This library gives you a glimpse of what [**SICStus
  897Prolog**](https://sicstus.sics.se/) can do. The API is intentionally
  898mostly compatible with that of SICStus Prolog, so that you can easily
  899switch to a much more feature-rich and much faster CLP(FD) system when
  900you need it. I thank [Mats Carlsson](https://www.sics.se/~matsc/), the
  901designer and main implementor of SICStus Prolog, for his elegant
  902example. I first encountered his system as part of the excellent
  903[**GUPU**](http://www.complang.tuwien.ac.at/ulrich/gupu/) teaching
  904environment by [Ulrich
  905Neumerkel](http://www.complang.tuwien.ac.at/ulrich/). Ulrich was also
  906the first and most determined tester of the present system, filing
  907hundreds of comments and suggestions for improvement. [Tom
  908Schrijvers](https://people.cs.kuleuven.be/~tom.schrijvers/) has
  909contributed several constraint libraries to SWI-Prolog, and I learned
  910a lot from his coding style and implementation examples. [Bart
  911Demoen](https://people.cs.kuleuven.be/~bart.demoen/) was a driving
  912force behind the implementation of attributed variables in SWI-Prolog,
  913and this library could not even have started without his prior work
  914and contributions. Thank you all!
  915
  916## CLP(FD) predicate index			{#clpfd-predicate-index}
  917
  918In the following, each CLP(FD) predicate is described in more detail.
  919
  920We recommend the following link to refer to this manual:
  921
  922http://eu.swi-prolog.org/man/clpfd.html
  923
  924@author [Markus Triska](https://www.metalevel.at)
  925*/
  926
  927:- create_prolog_flag(clpfd_monotonic, false, []).  928
  929/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  930   A bound is either:
  931
  932   n(N):    integer N
  933   inf:     infimum of Z (= negative infinity)
  934   sup:     supremum of Z (= positive infinity)
  935- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  936
  937is_bound(n(N)) :- integer(N).
  938is_bound(inf).
  939is_bound(sup).
  940
  941defaulty_to_bound(D, P) :- ( integer(D) -> P = n(D) ; P = D ).
  942
  943/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  944   Compactified is/2 and predicates for several arithmetic expressions
  945   with infinities, tailored for the modes needed by this solver.
  946- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  947
  948goal_expansion(A cis B, Expansion) :-
  949        phrase(cis_goals(B, A), Goals),
  950        list_goal(Goals, Expansion).
  951goal_expansion(A cis_lt B, B cis_gt A).
  952goal_expansion(A cis_leq B, B cis_geq A).
  953goal_expansion(A cis_geq B, cis_leq_numeric(B, N)) :- nonvar(A), A = n(N).
  954goal_expansion(A cis_geq B, cis_geq_numeric(A, N)) :- nonvar(B), B = n(N).
  955goal_expansion(A cis_gt B, cis_lt_numeric(B, N))   :- nonvar(A), A = n(N).
  956goal_expansion(A cis_gt B, cis_gt_numeric(A, N))   :- nonvar(B), B = n(N).
  957
  958% cis_gt only works for terms of depth 0 on both sides
  959cis_gt(sup, B0) :- B0 \== sup.
  960cis_gt(n(N), B) :- cis_lt_numeric(B, N).
  961
  962cis_lt_numeric(inf, _).
  963cis_lt_numeric(n(B), A) :- B < A.
  964
  965cis_gt_numeric(sup, _).
  966cis_gt_numeric(n(B), A) :- B > A.
  967
  968cis_geq(inf, inf).
  969cis_geq(sup, _).
  970cis_geq(n(N), B) :- cis_leq_numeric(B, N).
  971
  972cis_leq_numeric(inf, _).
  973cis_leq_numeric(n(B), A) :- B =< A.
  974
  975cis_geq_numeric(sup, _).
  976cis_geq_numeric(n(B), A) :- B >= A.
  977
  978cis_min(inf, _, inf).
  979cis_min(sup, B, B).
  980cis_min(n(N), B, Min) :- cis_min_(B, N, Min).
  981
  982cis_min_(inf, _, inf).
  983cis_min_(sup, N, n(N)).
  984cis_min_(n(B), A, n(M)) :- M is min(A,B).
  985
  986cis_max(sup, _, sup).
  987cis_max(inf, B, B).
  988cis_max(n(N), B, Max) :- cis_max_(B, N, Max).
  989
  990cis_max_(inf, N, n(N)).
  991cis_max_(sup, _, sup).
  992cis_max_(n(B), A, n(M)) :- M is max(A,B).
  993
  994cis_plus(inf, _, inf).
  995cis_plus(sup, _, sup).
  996cis_plus(n(A), B, Plus) :- cis_plus_(B, A, Plus).
  997
  998cis_plus_(sup, _, sup).
  999cis_plus_(inf, _, inf).
 1000cis_plus_(n(B), A, n(S)) :- S is A + B.
 1001
 1002cis_minus(inf, _, inf).
 1003cis_minus(sup, _, sup).
 1004cis_minus(n(A), B, M) :- cis_minus_(B, A, M).
 1005
 1006cis_minus_(inf, _, sup).
 1007cis_minus_(sup, _, inf).
 1008cis_minus_(n(B), A, n(M)) :- M is A - B.
 1009
 1010cis_uminus(inf, sup).
 1011cis_uminus(sup, inf).
 1012cis_uminus(n(A), n(B)) :- B is -A.
 1013
 1014cis_abs(inf, sup).
 1015cis_abs(sup, sup).
 1016cis_abs(n(A), n(B)) :- B is abs(A).
 1017
 1018cis_times(inf, B, P) :-
 1019        (   B cis_lt n(0) -> P = sup
 1020        ;   B cis_gt n(0) -> P = inf
 1021        ;   P = n(0)
 1022        ).
 1023cis_times(sup, B, P) :-
 1024        (   B cis_gt n(0) -> P = sup
 1025        ;   B cis_lt n(0) -> P = inf
 1026        ;   P = n(0)
 1027        ).
 1028cis_times(n(N), B, P) :- cis_times_(B, N, P).
 1029
 1030cis_times_(inf, A, P)     :- cis_times(inf, n(A), P).
 1031cis_times_(sup, A, P)     :- cis_times(sup, n(A), P).
 1032cis_times_(n(B), A, n(P)) :- P is A * B.
 1033
 1034cis_exp(inf, n(Y), R) :-
 1035        (   even(Y) -> R = sup
 1036        ;   R = inf
 1037        ).
 1038cis_exp(sup, _, sup).
 1039cis_exp(n(N), Y, R) :- cis_exp_(Y, N, R).
 1040
 1041cis_exp_(n(Y), N, n(R)) :- R is N^Y.
 1042cis_exp_(sup, _, sup).
 1043cis_exp_(inf, _, inf).
 1044
 1045cis_goals(V, V)          --> { var(V) }, !.
 1046cis_goals(n(N), n(N))    --> [].
 1047cis_goals(inf, inf)      --> [].
 1048cis_goals(sup, sup)      --> [].
 1049cis_goals(sign(A0), R)   --> cis_goals(A0, A), [cis_sign(A, R)].
 1050cis_goals(abs(A0), R)    --> cis_goals(A0, A), [cis_abs(A, R)].
 1051cis_goals(-A0, R)        --> cis_goals(A0, A), [cis_uminus(A, R)].
 1052cis_goals(A0+B0, R)      -->
 1053        cis_goals(A0, A),
 1054        cis_goals(B0, B),
 1055        [cis_plus(A, B, R)].
 1056cis_goals(A0-B0, R)      -->
 1057        cis_goals(A0, A),
 1058        cis_goals(B0, B),
 1059        [cis_minus(A, B, R)].
 1060cis_goals(min(A0,B0), R) -->
 1061        cis_goals(A0, A),
 1062        cis_goals(B0, B),
 1063        [cis_min(A, B, R)].
 1064cis_goals(max(A0,B0), R) -->
 1065        cis_goals(A0, A),
 1066        cis_goals(B0, B),
 1067        [cis_max(A, B, R)].
 1068cis_goals(A0*B0, R)      -->
 1069        cis_goals(A0, A),
 1070        cis_goals(B0, B),
 1071        [cis_times(A, B, R)].
 1072cis_goals(div(A0,B0), R) -->
 1073        cis_goals(A0, A),
 1074        cis_goals(B0, B),
 1075        [cis_div(A, B, R)].
 1076cis_goals(A0//B0, R)     -->
 1077        cis_goals(A0, A),
 1078        cis_goals(B0, B),
 1079        [cis_slash(A, B, R)].
 1080cis_goals(A0^B0, R)      -->
 1081        cis_goals(A0, A),
 1082        cis_goals(B0, B),
 1083        [cis_exp(A, B, R)].
 1084
 1085list_goal([], true).
 1086list_goal([G|Gs], Goal) :- foldl(list_goal_, Gs, G, Goal).
 1087
 1088list_goal_(G, G0, (G0,G)).
 1089
 1090cis_sign(sup, n(1)).
 1091cis_sign(inf, n(-1)).
 1092cis_sign(n(N), n(S)) :- S is sign(N).
 1093
 1094cis_div(sup, Y, Z)  :- ( Y cis_geq n(0) -> Z = sup ; Z = inf ).
 1095cis_div(inf, Y, Z)  :- ( Y cis_geq n(0) -> Z = inf ; Z = sup ).
 1096cis_div(n(X), Y, Z) :- cis_div_(Y, X, Z).
 1097
 1098cis_div_(sup, _, n(0)).
 1099cis_div_(inf, _, n(0)).
 1100cis_div_(n(Y), X, Z) :-
 1101        (   Y =:= 0 -> (  X >= 0 -> Z = sup ; Z = inf )
 1102        ;   Z0 is X // Y, Z = n(Z0)
 1103        ).
 1104
 1105cis_slash(sup, _, sup).
 1106cis_slash(inf, _, inf).
 1107cis_slash(n(N), B, S) :- cis_slash_(B, N, S).
 1108
 1109cis_slash_(sup, _, n(0)).
 1110cis_slash_(inf, _, n(0)).
 1111cis_slash_(n(B), A, n(S)) :- S is A // B.
 1112
 1113
 1114/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1115   A domain is a finite set of disjoint intervals. Internally, domains
 1116   are represented as trees. Each node is one of:
 1117
 1118   empty: empty domain.
 1119
 1120   split(N, Left, Right)
 1121      - split on integer N, with Left and Right domains whose elements are
 1122        all less than and greater than N, respectively. The domain is the
 1123        union of Left and Right, i.e., N is a hole.
 1124
 1125   from_to(From, To)
 1126      - interval (From-1, To+1); From and To are bounds
 1127
 1128   Desiderata: rebalance domains; singleton intervals.
 1129- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1130
 1131/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1132   Type definition and inspection of domains.
 1133- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1134
 1135check_domain(D) :-
 1136        (   var(D) -> instantiation_error(D)
 1137        ;   is_domain(D) -> true
 1138        ;   domain_error(clpfd_domain, D)
 1139        ).
 1140
 1141is_domain(empty).
 1142is_domain(from_to(From,To)) :-
 1143        is_bound(From), is_bound(To),
 1144        From cis_leq To.
 1145is_domain(split(S, Left, Right)) :-
 1146        integer(S),
 1147        is_domain(Left), is_domain(Right),
 1148        all_less_than(Left, S),
 1149        all_greater_than(Right, S).
 1150
 1151all_less_than(empty, _).
 1152all_less_than(from_to(From,To), S) :-
 1153        From cis_lt n(S), To cis_lt n(S).
 1154all_less_than(split(S0,Left,Right), S) :-
 1155        S0 < S,
 1156        all_less_than(Left, S),
 1157        all_less_than(Right, S).
 1158
 1159all_greater_than(empty, _).
 1160all_greater_than(from_to(From,To), S) :-
 1161        From cis_gt n(S), To cis_gt n(S).
 1162all_greater_than(split(S0,Left,Right), S) :-
 1163        S0 > S,
 1164        all_greater_than(Left, S),
 1165        all_greater_than(Right, S).
 1166
 1167default_domain(from_to(inf,sup)).
 1168
 1169domain_infimum(from_to(I, _), I).
 1170domain_infimum(split(_, Left, _), I) :- domain_infimum(Left, I).
 1171
 1172domain_supremum(from_to(_, S), S).
 1173domain_supremum(split(_, _, Right), S) :- domain_supremum(Right, S).
 1174
 1175domain_num_elements(empty, n(0)).
 1176domain_num_elements(from_to(From,To), Num) :- Num cis To - From + n(1).
 1177domain_num_elements(split(_, Left, Right), Num) :-
 1178        domain_num_elements(Left, NL),
 1179        domain_num_elements(Right, NR),
 1180        Num cis NL + NR.
 1181
 1182domain_direction_element(from_to(n(From), n(To)), Dir, E) :-
 1183        (   Dir == up -> between(From, To, E)
 1184        ;   between(From, To, E0),
 1185            E is To - (E0 - From)
 1186        ).
 1187domain_direction_element(split(_, D1, D2), Dir, E) :-
 1188        (   Dir == up ->
 1189            (   domain_direction_element(D1, Dir, E)
 1190            ;   domain_direction_element(D2, Dir, E)
 1191            )
 1192        ;   (   domain_direction_element(D2, Dir, E)
 1193            ;   domain_direction_element(D1, Dir, E)
 1194            )
 1195        ).
 1196
 1197/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1198   Test whether domain contains a given integer.
 1199- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1200
 1201domain_contains(from_to(From,To), I) :- From cis_leq n(I), n(I) cis_leq To.
 1202domain_contains(split(S, Left, Right), I) :-
 1203        (   I < S -> domain_contains(Left, I)
 1204        ;   I > S -> domain_contains(Right, I)
 1205        ).
 1206
 1207/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1208   Test whether a domain contains another domain.
 1209- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1210
 1211domain_subdomain(Dom, Sub) :- domain_subdomain(Dom, Dom, Sub).
 1212
 1213domain_subdomain(from_to(_,_), Dom, Sub) :-
 1214        domain_subdomain_fromto(Sub, Dom).
 1215domain_subdomain(split(_, _, _), Dom, Sub) :-
 1216        domain_subdomain_split(Sub, Dom, Sub).
 1217
 1218domain_subdomain_split(empty, _, _).
 1219domain_subdomain_split(from_to(From,To), split(S,Left0,Right0), Sub) :-
 1220        (   To cis_lt n(S) -> domain_subdomain(Left0, Left0, Sub)
 1221        ;   From cis_gt n(S) -> domain_subdomain(Right0, Right0, Sub)
 1222        ).
 1223domain_subdomain_split(split(_,Left,Right), Dom, _) :-
 1224        domain_subdomain(Dom, Dom, Left),
 1225        domain_subdomain(Dom, Dom, Right).
 1226
 1227domain_subdomain_fromto(empty, _).
 1228domain_subdomain_fromto(from_to(From,To), from_to(From0,To0)) :-
 1229        From0 cis_leq From, To0 cis_geq To.
 1230domain_subdomain_fromto(split(_,Left,Right), Dom) :-
 1231        domain_subdomain_fromto(Left, Dom),
 1232        domain_subdomain_fromto(Right, Dom).
 1233
 1234/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1235   Remove an integer from a domain. The domain is traversed until an
 1236   interval is reached from which the element can be removed, or until
 1237   it is clear that no such interval exists.
 1238- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1239
 1240domain_remove(empty, _, empty).
 1241domain_remove(from_to(L0, U0), X, D) :- domain_remove_(L0, U0, X, D).
 1242domain_remove(split(S, Left0, Right0), X, D) :-
 1243        (   X =:= S -> D = split(S, Left0, Right0)
 1244        ;   X < S ->
 1245            domain_remove(Left0, X, Left1),
 1246            (   Left1 == empty -> D = Right0
 1247            ;   D = split(S, Left1, Right0)
 1248            )
 1249        ;   domain_remove(Right0, X, Right1),
 1250            (   Right1 == empty -> D = Left0
 1251            ;   D = split(S, Left0, Right1)
 1252            )
 1253        ).
 1254
 1255%?- domain_remove(from_to(n(0),n(5)), 3, D).
 1256
 1257domain_remove_(inf, U0, X, D) :-
 1258        (   U0 == n(X) -> U1 is X - 1, D = from_to(inf, n(U1))
 1259        ;   U0 cis_lt n(X) -> D = from_to(inf,U0)
 1260        ;   L1 is X + 1, U1 is X - 1,
 1261            D = split(X, from_to(inf, n(U1)), from_to(n(L1),U0))
 1262        ).
 1263domain_remove_(n(N), U0, X, D) :- domain_remove_upper(U0, N, X, D).
 1264
 1265domain_remove_upper(sup, L0, X, D) :-
 1266        (   L0 =:= X -> L1 is X + 1, D = from_to(n(L1),sup)
 1267        ;   L0 > X -> D = from_to(n(L0),sup)
 1268        ;   L1 is X + 1, U1 is X - 1,
 1269            D = split(X, from_to(n(L0),n(U1)), from_to(n(L1),sup))
 1270        ).
 1271domain_remove_upper(n(U0), L0, X, D) :-
 1272        (   L0 =:= U0, X =:= L0 -> D = empty
 1273        ;   L0 =:= X -> L1 is X + 1, D = from_to(n(L1), n(U0))
 1274        ;   U0 =:= X -> U1 is X - 1, D = from_to(n(L0), n(U1))
 1275        ;   between(L0, U0, X) ->
 1276            U1 is X - 1, L1 is X + 1,
 1277            D = split(X, from_to(n(L0), n(U1)), from_to(n(L1), n(U0)))
 1278        ;   D = from_to(n(L0),n(U0))
 1279        ).
 1280
 1281/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1282   Remove all elements greater than / less than a constant.
 1283- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1284
 1285domain_remove_greater_than(empty, _, empty).
 1286domain_remove_greater_than(from_to(From0,To0), G, D) :-
 1287        (   From0 cis_gt n(G) -> D = empty
 1288        ;   To cis min(To0,n(G)), D = from_to(From0,To)
 1289        ).
 1290domain_remove_greater_than(split(S,Left0,Right0), G, D) :-
 1291        (   S =< G ->
 1292            domain_remove_greater_than(Right0, G, Right),
 1293            (   Right == empty -> D = Left0
 1294            ;   D = split(S, Left0, Right)
 1295            )
 1296        ;   domain_remove_greater_than(Left0, G, D)
 1297        ).
 1298
 1299domain_remove_smaller_than(empty, _, empty).
 1300domain_remove_smaller_than(from_to(From0,To0), V, D) :-
 1301        (   To0 cis_lt n(V) -> D = empty
 1302        ;   From cis max(From0,n(V)), D = from_to(From,To0)
 1303        ).
 1304domain_remove_smaller_than(split(S,Left0,Right0), V, D) :-
 1305        (   S >= V ->
 1306            domain_remove_smaller_than(Left0, V, Left),
 1307            (   Left == empty -> D = Right0
 1308            ;   D = split(S, Left, Right0)
 1309            )
 1310        ;   domain_remove_smaller_than(Right0, V, D)
 1311        ).
 1312
 1313
 1314/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1315   Remove a whole domain from another domain. (Set difference.)
 1316- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1317
 1318domain_subtract(Dom0, Sub, Dom) :- domain_subtract(Dom0, Dom0, Sub, Dom).
 1319
 1320domain_subtract(empty, _, _, empty).
 1321domain_subtract(from_to(From0,To0), Dom, Sub, D) :-
 1322        (   Sub == empty -> D = Dom
 1323        ;   Sub = from_to(From,To) ->
 1324            (   From == To -> From = n(X), domain_remove(Dom, X, D)
 1325            ;   From cis_gt To0 -> D = Dom
 1326            ;   To cis_lt From0 -> D = Dom
 1327            ;   From cis_leq From0 ->
 1328                (   To cis_geq To0 -> D = empty
 1329                ;   From1 cis To + n(1),
 1330                    D = from_to(From1, To0)
 1331                )
 1332            ;   To1 cis From - n(1),
 1333                (   To cis_lt To0 ->
 1334                    From = n(S),
 1335                    From2 cis To + n(1),
 1336                    D = split(S,from_to(From0,To1),from_to(From2,To0))
 1337                ;   D = from_to(From0,To1)
 1338                )
 1339            )
 1340        ;   Sub = split(S, Left, Right) ->
 1341            (   n(S) cis_gt To0 -> domain_subtract(Dom, Dom, Left, D)
 1342            ;   n(S) cis_lt From0 -> domain_subtract(Dom, Dom, Right, D)
 1343            ;   domain_subtract(Dom, Dom, Left, D1),
 1344                domain_subtract(D1, D1, Right, D)
 1345            )
 1346        ).
 1347domain_subtract(split(S, Left0, Right0), _, Sub, D) :-
 1348        domain_subtract(Left0, Left0, Sub, Left),
 1349        domain_subtract(Right0, Right0, Sub, Right),
 1350        (   Left == empty -> D = Right
 1351        ;   Right == empty -> D = Left
 1352        ;   D = split(S, Left, Right)
 1353        ).
 1354
 1355/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1356   Complement of a domain
 1357- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1358
 1359domain_complement(D, C) :-
 1360        default_domain(Default),
 1361        domain_subtract(Default, D, C).
 1362
 1363/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1364   Convert domain to a list of disjoint intervals From-To.
 1365- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1366
 1367domain_intervals(D, Is) :- phrase(domain_intervals(D), Is).
 1368
 1369domain_intervals(split(_, Left, Right)) -->
 1370        domain_intervals(Left), domain_intervals(Right).
 1371domain_intervals(empty)                 --> [].
 1372domain_intervals(from_to(From,To))      --> [From-To].
 1373
 1374/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1375   To compute the intersection of two domains D1 and D2, we choose D1
 1376   as the reference domain. For each interval of D1, we compute how
 1377   far and to which values D2 lets us extend it.
 1378- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1379
 1380domains_intersection(D1, D2, Intersection) :-
 1381        domains_intersection_(D1, D2, Intersection),
 1382        Intersection \== empty.
 1383
 1384domains_intersection_(empty, _, empty).
 1385domains_intersection_(from_to(L0,U0), D2, Dom) :-
 1386        narrow(D2, L0, U0, Dom).
 1387domains_intersection_(split(S,Left0,Right0), D2, Dom) :-
 1388        domains_intersection_(Left0, D2, Left1),
 1389        domains_intersection_(Right0, D2, Right1),
 1390        (   Left1 == empty -> Dom = Right1
 1391        ;   Right1 == empty -> Dom = Left1
 1392        ;   Dom = split(S, Left1, Right1)
 1393        ).
 1394
 1395narrow(empty, _, _, empty).
 1396narrow(from_to(L0,U0), From0, To0, Dom) :-
 1397        From1 cis max(From0,L0), To1 cis min(To0,U0),
 1398        (   From1 cis_gt To1 -> Dom = empty
 1399        ;   Dom = from_to(From1,To1)
 1400        ).
 1401narrow(split(S, Left0, Right0), From0, To0, Dom) :-
 1402        (   To0 cis_lt n(S) -> narrow(Left0, From0, To0, Dom)
 1403        ;   From0 cis_gt n(S) -> narrow(Right0, From0, To0, Dom)
 1404        ;   narrow(Left0, From0, To0, Left1),
 1405            narrow(Right0, From0, To0, Right1),
 1406            (   Left1 == empty -> Dom = Right1
 1407            ;   Right1 == empty -> Dom = Left1
 1408            ;   Dom = split(S, Left1, Right1)
 1409            )
 1410        ).
 1411
 1412/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1413   Union of 2 domains.
 1414- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1415
 1416domains_union(D1, D2, Union) :-
 1417        domain_intervals(D1, Is1),
 1418        domain_intervals(D2, Is2),
 1419        append(Is1, Is2, IsU0),
 1420        merge_intervals(IsU0, IsU1),
 1421        intervals_to_domain(IsU1, Union).
 1422
 1423/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1424   Shift the domain by an offset.
 1425- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1426
 1427domain_shift(empty, _, empty).
 1428domain_shift(from_to(From0,To0), O, from_to(From,To)) :-
 1429        From cis From0 + n(O), To cis To0 + n(O).
 1430domain_shift(split(S0, Left0, Right0), O, split(S, Left, Right)) :-
 1431        S is S0 + O,
 1432        domain_shift(Left0, O, Left),
 1433        domain_shift(Right0, O, Right).
 1434
 1435/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1436   The new domain contains all values of the old domain,
 1437   multiplied by a constant multiplier.
 1438- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1439
 1440domain_expand(D0, M, D) :-
 1441        (   M < 0 ->
 1442            domain_negate(D0, D1),
 1443            M1 is abs(M),
 1444            domain_expand_(D1, M1, D)
 1445        ;   M =:= 1 -> D = D0
 1446        ;   domain_expand_(D0, M, D)
 1447        ).
 1448
 1449domain_expand_(empty, _, empty).
 1450domain_expand_(from_to(From0, To0), M, from_to(From,To)) :-
 1451        From cis From0*n(M),
 1452        To cis To0*n(M).
 1453domain_expand_(split(S0, Left0, Right0), M, split(S, Left, Right)) :-
 1454        S is M*S0,
 1455        domain_expand_(Left0, M, Left),
 1456        domain_expand_(Right0, M, Right).
 1457
 1458/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1459   similar to domain_expand/3, tailored for truncated division: an
 1460   interval [From,To] is extended to [From*M, ((To+1)*M - 1)], i.e.,
 1461   to all values that truncated integer-divided by M yield a value
 1462   from interval.
 1463- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1464
 1465domain_expand_more(D0, M, D) :-
 1466        %format("expanding ~w by ~w\n", [D0,M]),
 1467        (   M < 0 -> domain_negate(D0, D1), M1 is abs(M)
 1468        ;   D1 = D0, M1 = M
 1469        ),
 1470        domain_expand_more_(D1, M1, D).
 1471        %format("yield: ~w\n", [D]).
 1472
 1473domain_expand_more_(empty, _, empty).
 1474domain_expand_more_(from_to(From0, To0), M, from_to(From,To)) :-
 1475        (   From0 cis_leq n(0) ->
 1476            From cis (From0-n(1))*n(M) + n(1)
 1477        ;   From cis From0*n(M)
 1478        ),
 1479        (   To0 cis_lt n(0) ->
 1480            To cis To0*n(M)
 1481        ;   To cis (To0+n(1))*n(M) - n(1)
 1482        ).
 1483domain_expand_more_(split(S0, Left0, Right0), M, split(S, Left, Right)) :-
 1484        S is M*S0,
 1485        domain_expand_more_(Left0, M, Left),
 1486        domain_expand_more_(Right0, M, Right).
 1487
 1488/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1489   Scale a domain down by a constant multiplier. Assuming (//)/2.
 1490- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1491
 1492domain_contract(D0, M, D) :-
 1493        %format("contracting ~w by ~w\n", [D0,M]),
 1494        (   M < 0 -> domain_negate(D0, D1), M1 is abs(M)
 1495        ;   D1 = D0, M1 = M
 1496        ),
 1497        domain_contract_(D1, M1, D).
 1498
 1499domain_contract_(empty, _, empty).
 1500domain_contract_(from_to(From0, To0), M, from_to(From,To)) :-
 1501        (   From0 cis_geq n(0) ->
 1502            From cis (From0 + n(M) - n(1)) // n(M)
 1503        ;   From cis From0 // n(M)
 1504        ),
 1505        (   To0 cis_geq n(0) ->
 1506            To cis To0 // n(M)
 1507        ;   To cis (To0 - n(M) + n(1)) // n(M)
 1508        ).
 1509domain_contract_(split(_,Left0,Right0), M, D) :-
 1510        %  Scaled down domains do not necessarily retain any holes of
 1511        %  the original domain.
 1512        domain_contract_(Left0, M, Left),
 1513        domain_contract_(Right0, M, Right),
 1514        domains_union(Left, Right, D).
 1515
 1516/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1517   Similar to domain_contract, tailored for division, i.e.,
 1518   {21,23} contracted by 4 is 5. It contracts "less".
 1519- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1520
 1521domain_contract_less(D0, M, D) :-
 1522        (   M < 0 -> domain_negate(D0, D1), M1 is abs(M)
 1523        ;   D1 = D0, M1 = M
 1524        ),
 1525        domain_contract_less_(D1, M1, D).
 1526
 1527domain_contract_less_(empty, _, empty).
 1528domain_contract_less_(from_to(From0, To0), M, from_to(From,To)) :-
 1529        From cis From0 // n(M), To cis To0 // n(M).
 1530domain_contract_less_(split(_,Left0,Right0), M, D) :-
 1531        %  Scaled down domains do not necessarily retain any holes of
 1532        %  the original domain.
 1533        domain_contract_less_(Left0, M, Left),
 1534        domain_contract_less_(Right0, M, Right),
 1535        domains_union(Left, Right, D).
 1536
 1537/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1538   Negate the domain. Left and Right sub-domains and bounds switch sides.
 1539- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1540
 1541domain_negate(empty, empty).
 1542domain_negate(from_to(From0, To0), from_to(From, To)) :-
 1543        From cis -To0, To cis -From0.
 1544domain_negate(split(S0, Left0, Right0), split(S, Left, Right)) :-
 1545        S is -S0,
 1546        domain_negate(Left0, Right),
 1547        domain_negate(Right0, Left).
 1548
 1549/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1550   Construct a domain from a list of integers. Try to balance it.
 1551- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1552
 1553list_to_disjoint_intervals([], []).
 1554list_to_disjoint_intervals([N|Ns], Is) :-
 1555        list_to_disjoint_intervals(Ns, N, N, Is).
 1556
 1557list_to_disjoint_intervals([], M, N, [n(M)-n(N)]).
 1558list_to_disjoint_intervals([B|Bs], M, N, Is) :-
 1559        (   B =:= N + 1 ->
 1560            list_to_disjoint_intervals(Bs, M, B, Is)
 1561        ;   Is = [n(M)-n(N)|Rest],
 1562            list_to_disjoint_intervals(Bs, B, B, Rest)
 1563        ).
 1564
 1565list_to_domain(List0, D) :-
 1566        (   List0 == [] -> D = empty
 1567        ;   sort(List0, List),
 1568            list_to_disjoint_intervals(List, Is),
 1569            intervals_to_domain(Is, D)
 1570        ).
 1571
 1572intervals_to_domain([], empty) :- !.
 1573intervals_to_domain([M-N], from_to(M,N)) :- !.
 1574intervals_to_domain(Is, D) :-
 1575        length(Is, L),
 1576        FL is L // 2,
 1577        length(Front, FL),
 1578        append(Front, Tail, Is),
 1579        Tail = [n(Start)-_|_],
 1580        Hole is Start - 1,
 1581        intervals_to_domain(Front, Left),
 1582        intervals_to_domain(Tail, Right),
 1583        D = split(Hole, Left, Right).
 1584
 1585%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 1586
 1587
 1588%% ?Var in +Domain
 1589%
 1590%  Var is an element of Domain. Domain is one of:
 1591%
 1592%         * Integer
 1593%           Singleton set consisting only of _Integer_.
 1594%         * Lower..Upper
 1595%           All integers _I_ such that _Lower_ =< _I_ =< _Upper_.
 1596%           _Lower_ must be an integer or the atom *inf*, which
 1597%           denotes negative infinity. _Upper_ must be an integer or
 1598%           the atom *sup*, which denotes positive infinity.
 1599%         * Domain1 \/ Domain2
 1600%           The union of Domain1 and Domain2.
 1601
 1602Var in Dom :- clpfd_in(Var, Dom).
 1603
 1604clpfd_in(V, D) :-
 1605        fd_variable(V),
 1606        drep_to_domain(D, Dom),
 1607        domain(V, Dom).
 1608
 1609fd_variable(V) :-
 1610        (   var(V) -> true
 1611        ;   integer(V) -> true
 1612        ;   type_error(integer, V)
 1613        ).
 1614
 1615%% +Vars ins +Domain
 1616%
 1617%  The variables in the list Vars are elements of Domain. See in/2 for
 1618%  the syntax of Domain.
 1619
 1620Vs ins D :-
 1621        fd_must_be_list(Vs),
 1622        maplist(fd_variable, Vs),
 1623        drep_to_domain(D, Dom),
 1624        domains(Vs, Dom).
 1625
 1626fd_must_be_list(Ls) :-
 1627        (   fd_var(Ls) -> type_error(list, Ls)
 1628        ;   must_be(list, Ls)
 1629        ).
 1630
 1631%% indomain(?Var)
 1632%
 1633% Bind Var to all feasible values of its domain on backtracking. The
 1634% domain of Var must be finite.
 1635
 1636indomain(Var) :- label([Var]).
 1637
 1638order_dom_next(up, Dom, Next)   :- domain_infimum(Dom, n(Next)).
 1639order_dom_next(down, Dom, Next) :- domain_supremum(Dom, n(Next)).
 1640order_dom_next(random_value(_), Dom, Next) :-
 1641        phrase(domain_to_intervals(Dom), Is),
 1642        length(Is, L),
 1643        R is random(L),
 1644        nth0(R, Is, From-To),
 1645        random_between(From, To, Next).
 1646
 1647domain_to_intervals(from_to(n(From),n(To))) --> [From-To].
 1648domain_to_intervals(split(_, Left, Right)) -->
 1649        domain_to_intervals(Left),
 1650        domain_to_intervals(Right).
 1651
 1652%% label(+Vars)
 1653%
 1654% Equivalent to labeling([], Vars). See labeling/2.
 1655
 1656label(Vs) :- labeling([], Vs).
 1657
 1658%% labeling(+Options, +Vars)
 1659%
 1660% Assign a value to each variable in Vars. Labeling means systematically
 1661% trying out values for the finite domain   variables  Vars until all of
 1662% them are ground. The domain of each   variable in Vars must be finite.
 1663% Options is a list of options that   let  you exhibit some control over
 1664% the search process. Several categories of options exist:
 1665%
 1666% The variable selection strategy lets you specify which variable of
 1667% Vars is labeled next and is one of:
 1668%
 1669%   * leftmost
 1670%   Label the variables in the order they occur in Vars. This is the
 1671%   default.
 1672%
 1673%   * ff
 1674%   _|First fail|_. Label the leftmost variable with smallest domain next,
 1675%   in order to detect infeasibility early. This is often a good
 1676%   strategy.
 1677%
 1678%   * ffc
 1679%   Of the variables with smallest domains, the leftmost one
 1680%   participating in most constraints is labeled next.
 1681%
 1682%   * min
 1683%   Label the leftmost variable whose lower bound is the lowest next.
 1684%
 1685%   * max
 1686%   Label the leftmost variable whose upper bound is the highest next.
 1687%
 1688% The value order is one of:
 1689%
 1690%   * up
 1691%   Try the elements of the chosen variable's domain in ascending order.
 1692%   This is the default.
 1693%
 1694%   * down
 1695%   Try the domain elements in descending order.
 1696%
 1697% The branching strategy is one of:
 1698%
 1699%   * step
 1700%   For each variable X, a choice is made between X = V and X #\= V,
 1701%   where V is determined by the value ordering options. This is the
 1702%   default.
 1703%
 1704%   * enum
 1705%   For each variable X, a choice is made between X = V_1, X = V_2
 1706%   etc., for all values V_i of the domain of X. The order is
 1707%   determined by the value ordering options.
 1708%
 1709%   * bisect
 1710%   For each variable X, a choice is made between X #=< M and X #> M,
 1711%   where M is the midpoint of the domain of X.
 1712%
 1713% At most one option of each category can be specified, and an option
 1714% must not occur repeatedly.
 1715%
 1716% The order of solutions can be influenced with:
 1717%
 1718%   * min(Expr)
 1719%   * max(Expr)
 1720%
 1721% This generates solutions in ascending/descending order with respect
 1722% to the evaluation of the arithmetic expression Expr. Labeling Vars
 1723% must make Expr ground. If several such options are specified, they
 1724% are interpreted from left to right, e.g.:
 1725%
 1726% ==
 1727% ?- [X,Y] ins 10..20, labeling([max(X),min(Y)],[X,Y]).
 1728% ==
 1729%
 1730% This generates solutions in descending order of X, and for each
 1731% binding of X, solutions are generated in ascending order of Y. To
 1732% obtain the incomplete behaviour that other systems exhibit with
 1733% "maximize(Expr)" and "minimize(Expr)", use once/1, e.g.:
 1734%
 1735% ==
 1736% once(labeling([max(Expr)], Vars))
 1737% ==
 1738%
 1739% Labeling is always complete, always terminates, and yields no
 1740% redundant solutions. See [core relations and
 1741% search](<#clpfd-search>) for usage advice.
 1742
 1743labeling(Options, Vars) :-
 1744        must_be(list, Options),
 1745        fd_must_be_list(Vars),
 1746        maplist(must_be_finite_fdvar, Vars),
 1747        label(Options, Options, default(leftmost), default(up), default(step), [], upto_ground, Vars).
 1748
 1749finite_domain(Dom) :-
 1750        domain_infimum(Dom, n(_)),
 1751        domain_supremum(Dom, n(_)).
 1752
 1753must_be_finite_fdvar(Var) :-
 1754        (   fd_get(Var, Dom, _) ->
 1755            (   finite_domain(Dom) -> true
 1756            ;   instantiation_error(Var)
 1757            )
 1758        ;   integer(Var) -> true
 1759        ;   must_be(integer, Var)
 1760        ).
 1761
 1762
 1763label([O|Os], Options, Selection, Order, Choice, Optim, Consistency, Vars) :-
 1764        (   var(O)-> instantiation_error(O)
 1765        ;   override(selection, Selection, O, Options, S1) ->
 1766            label(Os, Options, S1, Order, Choice, Optim, Consistency, Vars)
 1767        ;   override(order, Order, O, Options, O1) ->
 1768            label(Os, Options, Selection, O1, Choice, Optim, Consistency, Vars)
 1769        ;   override(choice, Choice, O, Options, C1) ->
 1770            label(Os, Options, Selection, Order, C1, Optim, Consistency, Vars)
 1771        ;   optimisation(O) ->
 1772            label(Os, Options, Selection, Order, Choice, [O|Optim], Consistency, Vars)
 1773        ;   consistency(O, O1) ->
 1774            label(Os, Options, Selection, Order, Choice, Optim, O1, Vars)
 1775        ;   domain_error(labeling_option, O)
 1776        ).
 1777label([], _, Selection, Order, Choice, Optim0, Consistency, Vars) :-
 1778        maplist(arg(1), [Selection,Order,Choice], [S,O,C]),
 1779        (   Optim0 == [] ->
 1780            label(Vars, S, O, C, Consistency)
 1781        ;   reverse(Optim0, Optim),
 1782            exprs_singlevars(Optim, SVs),
 1783            optimise(Vars, [S,O,C], SVs)
 1784        ).
 1785
 1786% Introduce new variables for each min/max expression to avoid
 1787% reparsing expressions during optimisation.
 1788
 1789exprs_singlevars([], []).
 1790exprs_singlevars([E|Es], [SV|SVs]) :-
 1791        E =.. [F,Expr],
 1792        ?(Single) #= Expr,
 1793        SV =.. [F,Single],
 1794        exprs_singlevars(Es, SVs).
 1795
 1796all_dead(fd_props(Bs,Gs,Os)) :-
 1797        all_dead_(Bs),
 1798        all_dead_(Gs),
 1799        all_dead_(Os).
 1800
 1801all_dead_([]).
 1802all_dead_([propagator(_, S)|Ps]) :- S == dead, all_dead_(Ps).
 1803
 1804label([], _, _, _, Consistency) :- !,
 1805        (   Consistency = upto_in(I0,I) -> I0 = I
 1806        ;   true
 1807        ).
 1808label(Vars, Selection, Order, Choice, Consistency) :-
 1809        (   Vars = [V|Vs], nonvar(V) -> label(Vs, Selection, Order, Choice, Consistency)
 1810        ;   select_var(Selection, Vars, Var, RVars),
 1811            (   var(Var) ->
 1812                (   Consistency = upto_in(I0,I), fd_get(Var, _, Ps), all_dead(Ps) ->
 1813                    fd_size(Var, Size),
 1814                    I1 is I0*Size,
 1815                    label(RVars, Selection, Order, Choice, upto_in(I1,I))
 1816                ;   Consistency = upto_in, fd_get(Var, _, Ps), all_dead(Ps) ->
 1817                    label(RVars, Selection, Order, Choice, Consistency)
 1818                ;   choice_order_variable(Choice, Order, Var, RVars, Vars, Selection, Consistency)
 1819                )
 1820            ;   label(RVars, Selection, Order, Choice, Consistency)
 1821            )
 1822        ).
 1823
 1824choice_order_variable(step, Order, Var, Vars, Vars0, Selection, Consistency) :-
 1825        fd_get(Var, Dom, _),
 1826        order_dom_next(Order, Dom, Next),
 1827        (   Var = Next,
 1828            label(Vars, Selection, Order, step, Consistency)
 1829        ;   neq_num(Var, Next),
 1830            do_queue,
 1831            label(Vars0, Selection, Order, step, Consistency)
 1832        ).
 1833choice_order_variable(enum, Order, Var, Vars, _, Selection, Consistency) :-
 1834        fd_get(Var, Dom0, _),
 1835        domain_direction_element(Dom0, Order, Var),
 1836        label(Vars, Selection, Order, enum, Consistency).
 1837choice_order_variable(bisect, Order, Var, _, Vars0, Selection, Consistency) :-
 1838        fd_get(Var, Dom, _),
 1839        domain_infimum(Dom, n(I)),
 1840        domain_supremum(Dom, n(S)),
 1841        Mid0 is (I + S) // 2,
 1842        (   Mid0 =:= S -> Mid is Mid0 - 1 ; Mid = Mid0 ),
 1843        (   Order == up -> ( Var #=< Mid ; Var #> Mid )
 1844        ;   Order == down -> ( Var #> Mid ; Var #=< Mid )
 1845        ;   domain_error(bisect_up_or_down, Order)
 1846        ),
 1847        label(Vars0, Selection, Order, bisect, Consistency).
 1848
 1849override(What, Prev, Value, Options, Result) :-
 1850        call(What, Value),
 1851        override_(Prev, Value, Options, Result).
 1852
 1853override_(default(_), Value, _, user(Value)).
 1854override_(user(Prev), Value, Options, _) :-
 1855        (   Value == Prev ->
 1856            domain_error(nonrepeating_labeling_options, Options)
 1857        ;   domain_error(consistent_labeling_options, Options)
 1858        ).
 1859
 1860selection(ff).
 1861selection(ffc).
 1862selection(min).
 1863selection(max).
 1864selection(leftmost).
 1865selection(random_variable(Seed)) :-
 1866        must_be(integer, Seed),
 1867        set_random(seed(Seed)).
 1868
 1869choice(step).
 1870choice(enum).
 1871choice(bisect).
 1872
 1873order(up).
 1874order(down).
 1875% TODO: random_variable and random_value currently both set the seed,
 1876% so exchanging the options can yield different results.
 1877order(random_value(Seed)) :-
 1878        must_be(integer, Seed),
 1879        set_random(seed(Seed)).
 1880
 1881consistency(upto_in(I), upto_in(1, I)).
 1882consistency(upto_in, upto_in).
 1883consistency(upto_ground, upto_ground).
 1884
 1885optimisation(min(_)).
 1886optimisation(max(_)).
 1887
 1888select_var(leftmost, [Var|Vars], Var, Vars).
 1889select_var(min, [V|Vs], Var, RVars) :-
 1890        find_min(Vs, V, Var),
 1891        delete_eq([V|Vs], Var, RVars).
 1892select_var(max, [V|Vs], Var, RVars) :-
 1893        find_max(Vs, V, Var),
 1894        delete_eq([V|Vs], Var, RVars).
 1895select_var(ff, [V|Vs], Var, RVars) :-
 1896        fd_size_(V, n(S)),
 1897        find_ff(Vs, V, S, Var),
 1898        delete_eq([V|Vs], Var, RVars).
 1899select_var(ffc, [V|Vs], Var, RVars) :-
 1900        find_ffc(Vs, V, Var),
 1901        delete_eq([V|Vs], Var, RVars).
 1902select_var(random_variable(_), Vars0, Var, Vars) :-
 1903        length(Vars0, L),
 1904        I is random(L),
 1905        nth0(I, Vars0, Var),
 1906        delete_eq(Vars0, Var, Vars).
 1907
 1908find_min([], Var, Var).
 1909find_min([V|Vs], CM, Min) :-
 1910        (   min_lt(V, CM) ->
 1911            find_min(Vs, V, Min)
 1912        ;   find_min(Vs, CM, Min)
 1913        ).
 1914
 1915find_max([], Var, Var).
 1916find_max([V|Vs], CM, Max) :-
 1917        (   max_gt(V, CM) ->
 1918            find_max(Vs, V, Max)
 1919        ;   find_max(Vs, CM, Max)
 1920        ).
 1921
 1922find_ff([], Var, _, Var).
 1923find_ff([V|Vs], CM, S0, FF) :-
 1924        (   nonvar(V) -> find_ff(Vs, CM, S0, FF)
 1925        ;   (   fd_size_(V, n(S1)), S1 < S0 ->
 1926                find_ff(Vs, V, S1, FF)
 1927            ;   find_ff(Vs, CM, S0, FF)
 1928            )
 1929        ).
 1930
 1931find_ffc([], Var, Var).
 1932find_ffc([V|Vs], Prev, FFC) :-
 1933        (   ffc_lt(V, Prev) ->
 1934            find_ffc(Vs, V, FFC)
 1935        ;   find_ffc(Vs, Prev, FFC)
 1936        ).
 1937
 1938
 1939ffc_lt(X, Y) :-
 1940        (   fd_get(X, XD, XPs) ->
 1941            domain_num_elements(XD, n(NXD))
 1942        ;   NXD = 1, XPs = []
 1943        ),
 1944        (   fd_get(Y, YD, YPs) ->
 1945            domain_num_elements(YD, n(NYD))
 1946        ;   NYD = 1, YPs = []
 1947        ),
 1948        (   NXD < NYD -> true
 1949        ;   NXD =:= NYD,
 1950            props_number(XPs, NXPs),
 1951            props_number(YPs, NYPs),
 1952            NXPs > NYPs
 1953        ).
 1954
 1955min_lt(X,Y) :- bounds(X,LX,_), bounds(Y,LY,_), LX < LY.
 1956
 1957max_gt(X,Y) :- bounds(X,_,UX), bounds(Y,_,UY), UX > UY.
 1958
 1959bounds(X, L, U) :-
 1960        (   fd_get(X, Dom, _) ->
 1961            domain_infimum(Dom, n(L)),
 1962            domain_supremum(Dom, n(U))
 1963        ;   L = X, U = L
 1964        ).
 1965
 1966delete_eq([], _, []).
 1967delete_eq([X|Xs], Y, List) :-
 1968        (   nonvar(X) -> delete_eq(Xs, Y, List)
 1969        ;   X == Y -> List = Xs
 1970        ;   List = [X|Tail],
 1971            delete_eq(Xs, Y, Tail)
 1972        ).
 1973
 1974/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1975   contracting/1 -- subject to change
 1976
 1977   This can remove additional domain elements from the boundaries.
 1978- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1979
 1980contracting(Vs) :-
 1981        must_be(list, Vs),
 1982        maplist(must_be_finite_fdvar, Vs),
 1983        contracting(Vs, false, Vs).
 1984
 1985contracting([], Repeat, Vars) :-
 1986        (   Repeat -> contracting(Vars, false, Vars)
 1987        ;   true
 1988        ).
 1989contracting([V|Vs], Repeat, Vars) :-
 1990        fd_inf(V, Min),
 1991        (   \+ \+ (V = Min) ->
 1992            fd_sup(V, Max),
 1993            (   \+ \+ (V = Max) ->
 1994                contracting(Vs, Repeat, Vars)
 1995            ;   V #\= Max,
 1996                contracting(Vs, true, Vars)
 1997            )
 1998        ;   V #\= Min,
 1999            contracting(Vs, true, Vars)
 2000        ).
 2001
 2002/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 2003   fds_sespsize(Vs, S).
 2004
 2005   S is an upper bound on the search space size with respect to finite
 2006   domain variables Vs.
 2007- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 2008
 2009fds_sespsize(Vs, S) :-
 2010        must_be(list, Vs),
 2011        maplist(fd_variable, Vs),
 2012        fds_sespsize(Vs, n(1), S1),
 2013        bound_portray(S1, S).
 2014
 2015fd_size_(V, S) :-
 2016        (   fd_get(V, D, _) ->
 2017            domain_num_elements(D, S)
 2018        ;   S = n(1)
 2019        ).
 2020
 2021fds_sespsize([], S, S).
 2022fds_sespsize([V|Vs], S0, S) :-
 2023        fd_size_(V, S1),
 2024        S2 cis S0*S1,
 2025        fds_sespsize(Vs, S2, S).
 2026
 2027/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 2028   Optimisation uses destructive assignment to save the computed
 2029   extremum over backtracking. Failure is used to get rid of copies of
 2030   attributed variables that are created in intermediate steps. At
 2031   least that's the intention - it currently doesn't work in SWI:
 2032
 2033   %?- X in 0..3, call_residue_vars(labeling([min(X)], [X]), Vs).
 2034   %@ X = 0,
 2035   %@ Vs = [_G6174, _G6177],
 2036   %@ _G6174 in 0..3
 2037
 2038- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 2039
 2040optimise(Vars, Options, Whats) :-
 2041        Whats = [What|WhatsRest],
 2042        Extremum = extremum(none),
 2043        (   catch(store_extremum(Vars, Options, What, Extremum),
 2044                  time_limit_exceeded,
 2045                  false)
 2046        ;   Extremum = extremum(n(Val)),
 2047            arg(1, What, Expr),
 2048            append(WhatsRest, Options, Options1),
 2049            (   Expr #= Val,
 2050                labeling(Options1, Vars)
 2051            ;   Expr #\= Val,
 2052                optimise(Vars, Options, Whats)
 2053            )
 2054        ).
 2055
 2056store_extremum(Vars, Options, What, Extremum) :-
 2057        catch((labeling(Options, Vars), throw(w(What))), w(What1), true),
 2058        functor(What, Direction, _),
 2059        maplist(arg(1), [What,What1], [Expr,Expr1]),
 2060        optimise(Direction, Options, Vars, Expr1, Expr, Extremum).
 2061
 2062optimise(Direction, Options, Vars, Expr0, Expr, Extremum) :-
 2063        must_be(ground, Expr0),
 2064        nb_setarg(1, Extremum, n(Expr0)),
 2065        catch((tighten(Direction, Expr, Expr0),
 2066               labeling(Options, Vars),
 2067               throw(v(Expr))), v(Expr1), true),
 2068        optimise(Direction, Options, Vars, Expr1, Expr, Extremum).
 2069
 2070tighten(min, E, V) :- E #< V.
 2071tighten(max, E, V) :- E #> V.
 2072
 2073%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 2074
 2075%% all_different(+Vars)
 2076%
 2077% Like all_distinct/1, but with weaker propagation. Consider using
 2078% all_distinct/1 instead, since all_distinct/1 is typically acceptably
 2079% efficient and propagates much more strongly.
 2080
 2081all_different(Ls) :-
 2082        fd_must_be_list(Ls),
 2083        maplist(fd_variable, Ls),
 2084        Orig = original_goal(_, all_different(Ls)),
 2085        all_different(Ls, [], Orig),
 2086        do_queue.
 2087
 2088all_different([], _, _).
 2089all_different([X|Right], Left, Orig) :-
 2090        (   var(X) ->
 2091            make_propagator(pdifferent(Left,Right,X,Orig), Prop),
 2092            init_propagator(X, Prop),
 2093            trigger_prop(Prop)
 2094        ;   exclude_fire(Left, Right, X)
 2095        ),
 2096        all_different(Right, [X|Left], Orig).
 2097
 2098%% all_distinct(+Vars).
 2099%
 2100%  True iff Vars are pairwise distinct. For example, all_distinct/1
 2101%  can detect that not all variables can assume distinct values given
 2102%  the following domains:
 2103%
 2104%  ==
 2105%  ?- maplist(in, Vs,
 2106%             [1\/3..4, 1..2\/4, 1..2\/4, 1..3, 1..3, 1..6]),
 2107%     all_distinct(Vs).
 2108%  false.
 2109%  ==
 2110
 2111all_distinct(Ls) :-
 2112        fd_must_be_list(Ls),
 2113        maplist(fd_variable, Ls),
 2114        make_propagator(pdistinct(Ls), Prop),
 2115        distinct_attach(Ls, Prop, []),
 2116        trigger_once(Prop).
 2117
 2118%% sum(+Vars, +Rel, ?Expr)
 2119%
 2120% The sum of elements of the list Vars is in relation Rel to Expr.
 2121% Rel is one of #=, #\=, #<, #>, #=< or #>=. For example:
 2122%
 2123% ==
 2124% ?- [A,B,C] ins 0..sup, sum([A,B,C], #=, 100).
 2125% A in 0..100,
 2126% A+B+C#=100,
 2127% B in 0..100,
 2128% C in 0..100.
 2129% ==
 2130
 2131sum(Vs, Op, Value) :-
 2132        must_be(list, Vs),
 2133        same_length(Vs, Ones),
 2134        maplist(=(1), Ones),
 2135        scalar_product(Ones, Vs, Op, Value).
 2136
 2137%% scalar_product(+Cs, +Vs, +Rel, ?Expr)
 2138%
 2139% True iff the scalar product of Cs and Vs is in relation Rel to Expr.
 2140% Cs is a list of integers, Vs is a list of variables and integers.
 2141% Rel is #=, #\=, #<, #>, #=< or #>=.
 2142
 2143scalar_product(Cs, Vs, Op, Value) :-
 2144        must_be(list(integer), Cs),
 2145        must_be(list, Vs),
 2146        maplist(fd_variable, Vs),
 2147        (   Op = (#=), single_value(Value, Right), ground(Vs) ->
 2148            foldl(coeff_int_linsum, Cs, Vs, 0, Right)
 2149        ;   must_be(callable, Op),
 2150            (   memberchk(Op, [#=,#\=,#<,#>,#=<,#>=]) -> true
 2151            ;   domain_error(scalar_product_relation, Op)
 2152            ),
 2153            must_be(acyclic, Value),
 2154            foldl(coeff_var_plusterm, Cs, Vs, 0, Left),
 2155            (   left_right_linsum_const(Left, Value, Cs1, Vs1, Const) ->
 2156                scalar_product_(Op, Cs1, Vs1, Const)
 2157            ;   sum(Cs, Vs, 0, Op, Value)
 2158            )
 2159        ).
 2160
 2161single_value(V, V)    :- var(V), !, non_monotonic(V).
 2162single_value(V, V)    :- integer(V).
 2163single_value(?(V), V) :- fd_variable(V).
 2164
 2165coeff_var_plusterm(C, V, T0, T0+(C* ?(V))).
 2166
 2167coeff_int_linsum(C, I, S0, S) :- S is S0 + C*I.
 2168
 2169sum([], _, Sum, Op, Value) :- call(Op, Sum, Value).
 2170sum([C|Cs], [X|Xs], Acc, Op, Value) :-
 2171        ?(NAcc) #= Acc + C* ?(X),
 2172        sum(Cs, Xs, NAcc, Op, Value).
 2173
 2174multiples([], [], _).
 2175multiples([C|Cs], [V|Vs], Left) :-
 2176        (   (   Cs = [N|_] ; Left = [N|_] ) ->
 2177            (   N =\= 1, gcd(C,N) =:= 1 ->
 2178                gcd(Cs, N, GCD0),
 2179                gcd(Left, GCD0, GCD),
 2180                (   GCD > 1 -> ?(V) #= GCD * ?(_)
 2181                ;   true
 2182                )
 2183            ;   true
 2184            )
 2185        ;   true
 2186        ),
 2187        multiples(Cs, Vs, [C|Left]).
 2188
 2189abs(N, A) :- A is abs(N).
 2190
 2191divide(D, N, R) :- R is N // D.
 2192
 2193scalar_product_(#=, Cs0, Vs, S0) :-
 2194        (   Cs0 = [C|Rest] ->
 2195            gcd(Rest, C, GCD),
 2196            S0 mod GCD =:= 0,
 2197            maplist(divide(GCD), [S0|Cs0], [S|Cs])
 2198        ;   S0 =:= 0, S = S0, Cs = Cs0
 2199        ),
 2200        (   S0 =:= 0 ->
 2201            maplist(abs, Cs, As),
 2202            multiples(As, Vs, [])
 2203        ;   true
 2204        ),
 2205        propagator_init_trigger(Vs, scalar_product_eq(Cs, Vs, S)).
 2206scalar_product_(#\=, Cs, Vs, C) :-
 2207        propagator_init_trigger(Vs, scalar_product_neq(Cs, Vs, C)).
 2208scalar_product_(#=<, Cs, Vs, C) :-
 2209        propagator_init_trigger(Vs, scalar_product_leq(Cs, Vs, C)).
 2210scalar_product_(#<, Cs, Vs, C) :-
 2211        C1 is C - 1,
 2212        scalar_product_(#=<, Cs, Vs, C1).
 2213scalar_product_(#>, Cs, Vs, C) :-
 2214        C1 is C + 1,
 2215        scalar_product_(#>=, Cs, Vs, C1).
 2216scalar_product_(#>=, Cs, Vs, C) :-
 2217        maplist(negative, Cs, Cs1),
 2218        C1 is -C,
 2219        scalar_product_(#=<, Cs1, Vs, C1).
 2220
 2221negative(X0, X) :- X is -X0.
 2222
 2223coeffs_variables_const([], [], [], [], I, I).
 2224coeffs_variables_const([C|Cs], [V|Vs], Cs1, Vs1, I0, I) :-
 2225        (   var(V) ->
 2226            Cs1 = [C|CRest], Vs1 = [V|VRest], I1 = I0
 2227        ;   I1 is I0 + C*V,
 2228            Cs1 = CRest, Vs1 = VRest
 2229        ),
 2230        coeffs_variables_const(Cs, Vs, CRest, VRest, I1, I).
 2231
 2232sum_finite_domains([], [], [], [], Inf, Sup, Inf, Sup).
 2233sum_finite_domains([C|Cs], [V|Vs], Infs, Sups, Inf0, Sup0, Inf, Sup) :-
 2234        fd_get(V, _, Inf1, Sup1, _),
 2235        (   Inf1 = n(NInf) ->
 2236            (   C < 0 ->
 2237                Sup2 is Sup0 + C*NInf
 2238            ;   Inf2 is Inf0 + C*NInf
 2239            ),
 2240            Sups = Sups1,
 2241            Infs = Infs1
 2242        ;   (   C < 0 ->
 2243                Sup2 = Sup0,
 2244                Sups = [C*V|Sups1],
 2245                Infs = Infs1
 2246            ;   Inf2 = Inf0,
 2247                Infs = [C*V|Infs1],
 2248                Sups = Sups1
 2249            )
 2250        ),
 2251        (   Sup1 = n(NSup) ->
 2252            (   C < 0 ->
 2253                Inf2 is Inf0 + C*NSup
 2254            ;   Sup2 is Sup0 + C*NSup
 2255            ),
 2256            Sups1 = Sups2,
 2257            Infs1 = Infs2
 2258        ;   (   C < 0 ->
 2259                Inf2 = Inf0,
 2260                Infs1 = [C*V|Infs2],
 2261                Sups1 = Sups2
 2262            ;   Sup2 = Sup0,
 2263                Sups1 = [C*V|Sups2],
 2264                Infs1 = Infs2
 2265            )
 2266        ),
 2267        sum_finite_domains(Cs, Vs, Infs2, Sups2, Inf2, Sup2, Inf, Sup).
 2268
 2269remove_dist_upper_lower([], _, _, _).
 2270remove_dist_upper_lower([C|Cs], [V|Vs], D1, D2) :-
 2271        (   fd_get(V, VD, VPs) ->
 2272            (   C < 0 ->
 2273                domain_supremum(VD, n(Sup)),
 2274                L is Sup + D1//C,
 2275                domain_remove_smaller_than(VD, L, VD1),
 2276                domain_infimum(VD1, n(Inf)),
 2277                G is Inf - D2//C,
 2278                domain_remove_greater_than(VD1, G, VD2)
 2279            ;   domain_infimum(VD, n(Inf)),
 2280                G is Inf + D1//C,
 2281                domain_remove_greater_than(VD, G, VD1),
 2282                domain_supremum(VD1, n(Sup)),
 2283                L is Sup - D2//C,
 2284                domain_remove_smaller_than(VD1, L, VD2)
 2285            ),
 2286            fd_put(V, VD2, VPs)
 2287        ;   true
 2288        ),
 2289        remove_dist_upper_lower(Cs, Vs, D1, D2).
 2290
 2291
 2292remove_dist_upper_leq([], _, _).
 2293remove_dist_upper_leq([C|Cs], [V|Vs], D1) :-
 2294        (   fd_get(V, VD, VPs) ->
 2295            (   C < 0 ->
 2296                domain_supremum(VD, n(Sup)),
 2297                L is Sup + D1//C,
 2298                domain_remove_smaller_than(VD, L, VD1)
 2299            ;   domain_infimum(VD, n(Inf)),
 2300                G is Inf + D1//C,
 2301                domain_remove_greater_than(VD, G, VD1)
 2302            ),
 2303            fd_put(V, VD1, VPs)
 2304        ;   true
 2305        ),
 2306        remove_dist_upper_leq(Cs, Vs, D1).
 2307
 2308
 2309remove_dist_upper([], _).
 2310remove_dist_upper([C*V|CVs], D) :-
 2311        (   fd_get(V, VD, VPs) ->
 2312            (   C < 0 ->
 2313                (   domain_supremum(VD, n(Sup)) ->
 2314                    L is Sup + D//C,
 2315                    domain_remove_smaller_than(VD, L, VD1)
 2316                ;   VD1 = VD
 2317                )
 2318            ;   (   domain_infimum(VD, n(Inf)) ->
 2319                    G is Inf + D//C,
 2320                    domain_remove_greater_than(VD, G, VD1)
 2321                ;   VD1 = VD
 2322                )
 2323            ),
 2324            fd_put(V, VD1, VPs)
 2325        ;   true
 2326        ),
 2327        remove_dist_upper(CVs, D).
 2328
 2329remove_dist_lower([], _).
 2330remove_dist_lower([C*V|CVs], D) :-
 2331        (   fd_get(V, VD, VPs) ->
 2332            (   C < 0 ->
 2333                (   domain_infimum(VD, n(Inf)) ->
 2334                    G is Inf - D//C,
 2335                    domain_remove_greater_than(VD, G, VD1)
 2336                ;   VD1 = VD
 2337                )
 2338            ;   (   domain_supremum(VD, n(Sup)) ->
 2339                    L is Sup - D//C,
 2340                    domain_remove_smaller_than(VD, L, VD1)
 2341                ;   VD1 = VD
 2342                )
 2343            ),
 2344            fd_put(V, VD1, VPs)
 2345        ;   true
 2346        ),
 2347        remove_dist_lower(CVs, D).
 2348
 2349remove_upper([], _).
 2350remove_upper([C*X|CXs], Max) :-
 2351        (   fd_get(X, XD, XPs) ->
 2352            D is Max//C,
 2353            (   C < 0 ->
 2354                domain_remove_smaller_than(XD, D, XD1)
 2355            ;   domain_remove_greater_than(XD, D, XD1)
 2356            ),
 2357            fd_put(X, XD1, XPs)
 2358        ;   true
 2359        ),
 2360        remove_upper(CXs, Max).
 2361
 2362remove_lower([], _).
 2363remove_lower([C*X|CXs], Min) :-
 2364        (   fd_get(X, XD, XPs) ->
 2365            D is -Min//C,
 2366            (   C < 0 ->
 2367                domain_remove_greater_than(XD, D, XD1)
 2368            ;   domain_remove_smaller_than(XD, D, XD1)
 2369            ),
 2370            fd_put(X, XD1, XPs)
 2371        ;   true
 2372        ),
 2373        remove_lower(CXs, Min).
 2374
 2375%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 2376
 2377/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 2378   Constraint propagation proceeds as follows: Each CLP(FD) variable
 2379   has an attribute that stores its associated domain and constraints.
 2380   Constraints are triggered when the event they are registered for
 2381   occurs (for example: variable is instantiated, bounds change etc.).
 2382   do_queue/0 works off all triggered constraints, possibly triggering
 2383   new ones, until fixpoint.
 2384- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 2385
 2386% FIFO queue
 2387
 2388make_queue :- nb_setval('$clpfd_queue', fast_slow([], [])).
 2389
 2390push_queue(E, Which) :-
 2391        nb_getval('$clpfd_queue', Qs),
 2392        arg(Which, Qs, Q),
 2393        (   Q == [] ->
 2394            setarg(Which, Qs, [E|T]-T)
 2395        ;   Q = H-[E|T],
 2396            setarg(Which, Qs, H-T)
 2397        ).
 2398
 2399pop_queue(E) :-
 2400        nb_getval('$clpfd_queue', Qs),
 2401        (   pop_queue(E, Qs, 1) ->  true
 2402        ;   pop_queue(E, Qs, 2)
 2403        ).
 2404
 2405pop_queue(E, Qs, Which) :-
 2406        arg(Which, Qs, [E|NH]-T),
 2407        (   var(NH) ->
 2408            setarg(Which, Qs, [])
 2409        ;   setarg(Which, Qs, NH-T)
 2410        ).
 2411
 2412fetch_propagator(Prop) :-
 2413        pop_queue(P),
 2414        (   propagator_state(P, S), S == dead -> fetch_propagator(Prop)
 2415        ;   Prop = P
 2416        ).
 2417
 2418/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 2419   Parsing a CLP(FD) expression has two important side-effects: First,
 2420   it constrains the variables occurring in the expression to
 2421   integers. Second, it constrains some of them even more: For
 2422   example, in X/Y and X mod Y, Y is constrained to be #\= 0.
 2423- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 2424
 2425constrain_to_integer(Var) :-
 2426        (   integer(Var) -> true
 2427        ;   fd_get(Var, D, Ps),
 2428            fd_put(Var, D, Ps)
 2429        ).
 2430
 2431power_var_num(P, X, N) :-
 2432        (   var(P) -> X = P, N = 1
 2433        ;   P = Left*Right,
 2434            power_var_num(Left, XL, L),
 2435            power_var_num(Right, XR, R),
 2436            XL == XR,
 2437            X = XL,
 2438            N is L + R
 2439        ).
 2440
 2441/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 2442   Given expression E, we obtain the finite domain variable R by
 2443   interpreting a simple committed-choice language that is a list of
 2444   conditions and bodies. In conditions, g(Goal) means literally Goal,
 2445   and m(Match) means that E can be decomposed as stated. The
 2446   variables are to be understood as the result of parsing the
 2447   subexpressions recursively. In the body, g(Goal) means again Goal,
 2448   and p(Propagator) means to attach and trigger once a propagator.
 2449- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 2450
 2451:- op(800, xfx, =>). 2452
 2453parse_clpfd(E, R,
 2454            [g(cyclic_term(E)) => [g(domain_error(clpfd_expression, E))],
 2455             g(var(E))         => [g(non_monotonic(E)),
 2456                                   g(constrain_to_integer(E)), g(E = R)],
 2457             g(integer(E))     => [g(R = E)],
 2458             ?(E)              => [g(must_be_fd_integer(E)), g(R = E)],
 2459             #(E)              => [g(must_be_fd_integer(E)), g(R = E)],
 2460             m(A+B)            => [p(pplus(A, B, R))],
 2461             % power_var_num/3 must occur before */2 to be useful
 2462             g(power_var_num(E, V, N)) => [p(pexp(V, N, R))],
 2463             m(A*B)            => [p(ptimes(A, B, R))],
 2464             m(A-B)            => [p(pplus(R,B,A))],
 2465             m(-A)             => [p(ptimes(-1,A,R))],
 2466             m(max(A,B))       => [g(A #=< ?(R)), g(B #=< R), p(pmax(A, B, R))],
 2467             m(min(A,B))       => [g(A #>= ?(R)), g(B #>= R), p(pmin(A, B, R))],
 2468             m(A mod B)        => [g(B #\= 0), p(pmod(A, B, R))],
 2469             m(A rem B)        => [g(B #\= 0), p(prem(A, B, R))],
 2470             m(abs(A))         => [g(?(R) #>= 0), p(pabs(A, R))],
 2471%             m(A/B)            => [g(B #\= 0), p(ptzdiv(A, B, R))],
 2472             m(A//B)           => [g(B #\= 0), p(ptzdiv(A, B, R))],
 2473             m(A div B)        => [g(?(R) #= (A - (A mod B)) // B)],
 2474             m(A rdiv B)       => [g(B #\= 0), p(prdiv(A, B, R))],
 2475             m(A^B)            => [p(pexp(A, B, R))],
 2476             % bitwise operations
 2477             m(\A)             => [p(pfunction(\, A, R))],
 2478             m(msb(A))         => [p(pfunction(msb, A, R))],
 2479             m(lsb(A))         => [p(pfunction(lsb, A, R))],
 2480             m(popcount(A))    => [p(pfunction(popcount, A, R))],
 2481             m(A<<B)           => [p(pfunction(<<, A, B, R))],
 2482             m(A>>B)           => [p(pfunction(>>, A, B, R))],
 2483             m(A/\B)           => [p(pfunction(/\, A, B, R))],
 2484             m(A\/B)           => [p(pfunction(\/, A, B, R))],
 2485             m(A xor B)        => [p(pfunction(xor, A, B, R))],
 2486             g(true)           => [g(domain_error(clpfd_expression, E))]
 2487            ]).
 2488
 2489non_monotonic(X) :-
 2490        (   \+ fd_var(X), current_prolog_flag(clpfd_monotonic, true) ->
 2491            instantiation_error(X)
 2492        ;   true
 2493        ).
 2494
 2495% Here, we compile the committed choice language to a single
 2496% predicate, parse_clpfd/2.
 2497
 2498make_parse_clpfd(Clauses) :-
 2499        parse_clpfd_clauses(Clauses0),
 2500        maplist(goals_goal, Clauses0, Clauses).
 2501
 2502goals_goal((Head :- Goals), (Head :- Body)) :-
 2503        list_goal(Goals, Body).
 2504
 2505parse_clpfd_clauses(Clauses) :-
 2506        parse_clpfd(E, R, Matchers),
 2507        maplist(parse_matcher(E, R), Matchers, Clauses).
 2508
 2509parse_matcher(E, R, Matcher, Clause) :-
 2510        Matcher = (Condition0 => Goals0),
 2511        phrase((parse_condition(Condition0, E, Head),
 2512                parse_goals(Goals0)), Goals),
 2513        Clause = (parse_clpfd(Head, R) :- Goals).
 2514
 2515parse_condition(g(Goal), E, E)       --> [Goal, !].
 2516parse_condition(?(E), _, ?(E))       --> [!].
 2517parse_condition(#(E), _, #(E))       --> [!].
 2518parse_condition(m(Match), _, Match0) -->
 2519        [!],
 2520        { copy_term(Match, Match0),
 2521          term_variables(Match0, Vs0),
 2522          term_variables(Match, Vs)
 2523        },
 2524        parse_match_variables(Vs0, Vs).
 2525
 2526parse_match_variables([], []) --> [].
 2527parse_match_variables([V0|Vs0], [V|Vs]) -->
 2528        [parse_clpfd(V0, V)],
 2529        parse_match_variables(Vs0, Vs).
 2530
 2531parse_goals([]) --> [].
 2532parse_goals([G|Gs]) --> parse_goal(G), parse_goals(Gs).
 2533
 2534parse_goal(g(Goal)) --> [Goal].
 2535parse_goal(p(Prop)) -->
 2536        [make_propagator(Prop, P)],
 2537        { term_variables(Prop, Vs) },
 2538        parse_init(Vs, P),
 2539        [trigger_once(P)].
 2540
 2541parse_init([], _)     --> [].
 2542parse_init([V|Vs], P) --> [init_propagator(V, P)], parse_init(Vs, P).
 2543
 2544%?- set_prolog_flag(answer_write_options, [portray(true)]),
 2545%   clpfd:parse_clpfd_clauses(Clauses), maplist(portray_clause, Clauses).
 2546
 2547
 2548%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 2549%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 2550
 2551trigger_once(Prop) :- trigger_prop(Prop), do_queue.
 2552
 2553neq(A, B) :- propagator_init_trigger(pneq(A, B)).
 2554
 2555propagator_init_trigger(P) -->
 2556        { term_variables(P, Vs) },
 2557        propagator_init_trigger(Vs, P).
 2558
 2559propagator_init_trigger(Vs, P) -->
 2560        [p(Prop)],
 2561        { make_propagator(P, Prop),
 2562          maplist(prop_init(Prop), Vs),
 2563          trigger_once(Prop) }.
 2564
 2565propagator_init_trigger(P) :-
 2566        phrase(propagator_init_trigger(P), _).
 2567
 2568propagator_init_trigger(Vs, P) :-
 2569        phrase(propagator_init_trigger(Vs, P), _).
 2570
 2571prop_init(Prop, V) :- init_propagator(V, Prop).
 2572
 2573geq(A, B) :-
 2574        (   fd_get(A, AD, APs) ->
 2575            domain_infimum(AD, AI),
 2576            (   fd_get(B, BD, _) ->
 2577                domain_supremum(BD, BS),
 2578                (   AI cis_geq BS -> true
 2579                ;   propagator_init_trigger(pgeq(A,B))
 2580                )
 2581            ;   (   AI cis_geq n(B) -> true
 2582                ;   domain_remove_smaller_than(AD, B, AD1),
 2583                    fd_put(A, AD1, APs),
 2584                    do_queue
 2585                )
 2586            )
 2587        ;   fd_get(B, BD, BPs) ->
 2588            domain_remove_greater_than(BD, A, BD1),
 2589            fd_put(B, BD1, BPs),
 2590            do_queue
 2591        ;   A >= B
 2592        ).
 2593
 2594/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 2595   Naive parsing of inequalities and disequalities can result in a lot
 2596   of unnecessary work if expressions of non-trivial depth are
 2597   involved: Auxiliary variables are introduced for sub-expressions,
 2598   and propagation proceeds on them as if they were involved in a
 2599   tighter constraint (like equality), whereas eventually only very
 2600   little of the propagated information is actually used. For example,
 2601   only extremal values are of interest in inequalities. Introducing
 2602   auxiliary variables should be avoided when possible, and
 2603   specialised propagators should be used for common constraints.
 2604
 2605   We again use a simple committed-choice language for matching
 2606   special cases of constraints. m_c(M,C) means that M matches and C
 2607   holds. d(X, Y) means decomposition, i.e., it is short for
 2608   g(parse_clpfd(X, Y)). r(X, Y) means to rematch with X and Y.
 2609
 2610   Two things are important: First, although the actual constraint
 2611   functors (#\=2, #=/2 etc.) are used in the description, they must
 2612   expand to the respective auxiliary predicates (match_expand/2)
 2613   because the actual constraints are subject to goal expansion.
 2614   Second, when specialised constraints (like scalar product) post
 2615   simpler constraints on their own, these simpler versions must be
 2616   handled separately and must occur before.
 2617- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 2618
 2619match_expand(#>=, clpfd_geq_).
 2620match_expand(#=, clpfd_equal_).
 2621match_expand(#\=, clpfd_neq).
 2622
 2623symmetric(#=).
 2624symmetric(#\=).
 2625
 2626matches([
 2627         m_c(any(X) #>= any(Y), left_right_linsum_const(X, Y, Cs, Vs, Const)) =>
 2628            [g((   Cs = [1], Vs = [A] -> geq(A, Const)
 2629               ;   Cs = [-1], Vs = [A] -> Const1 is -Const, geq(Const1, A)
 2630               ;   Cs = [1,1], Vs = [A,B] -> ?(A) + ?(B) #= ?(S), geq(S, Const)
 2631               ;   Cs = [1,-1], Vs = [A,B] ->
 2632                   (   Const =:= 0 -> geq(A, B)
 2633                   ;   C1 is -Const,
 2634                       propagator_init_trigger(x_leq_y_plus_c(B, A, C1))
 2635                   )
 2636               ;   Cs = [-1,1], Vs = [A,B] ->
 2637                   (   Const =:= 0 -> geq(B, A)
 2638                   ;   C1 is -Const,
 2639                       propagator_init_trigger(x_leq_y_plus_c(A, B, C1))
 2640                   )
 2641               ;   Cs = [-1,-1], Vs = [A,B] ->
 2642                   ?(A) + ?(B) #= ?(S), Const1 is -Const, geq(Const1, S)
 2643               ;   scalar_product_(#>=, Cs, Vs, Const)
 2644               ))],
 2645         m(any(X) - any(Y) #>= integer(C))     => [d(X, X1), d(Y, Y1), g(C1 is -C), p(x_leq_y_plus_c(Y1, X1, C1))],
 2646         m(integer(X) #>= any(Z) + integer(A)) => [g(C is X - A), r(C, Z)],
 2647         m(abs(any(X)-any(Y)) #>= integer(I))  => [d(X, X1), d(Y, Y1), p(absdiff_geq(X1, Y1, I))],
 2648         m(abs(any(X)) #>= integer(I))         => [d(X, RX), g((I>0 -> I1 is -I, RX in inf..I1 \/ I..sup; true))],
 2649         m(integer(I) #>= abs(any(X)))         => [d(X, RX), g(I>=0), g(I1 is -I), g(RX in I1..I)],
 2650         m(any(X) #>= any(Y))                  => [d(X, RX), d(Y, RY), g(geq(RX, RY))],
 2651
 2652         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 2653         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 2654
 2655         m(var(X) #= var(Y))        => [g(constrain_to_integer(X)), g(X=Y)],
 2656         m(var(X) #= var(Y)+var(Z)) => [p(pplus(Y,Z,X))],
 2657         m(var(X) #= var(Y)-var(Z)) => [p(pplus(X,Z,Y))],
 2658         m(var(X) #= var(Y)*var(Z)) => [p(ptimes(Y,Z,X))],
 2659         m(var(X) #= -var(Z))       => [p(ptimes(-1, Z, X))],
 2660         m_c(any(X) #= any(Y), left_right_linsum_const(X, Y, Cs, Vs, S)) =>
 2661            [g(scalar_product_(#=, Cs, Vs, S))],
 2662         m(var(X) #= any(Y))       => [d(Y,X)],
 2663         m(any(X) #= any(Y))       => [d(X, RX), d(Y, RX)],
 2664
 2665         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 2666         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 2667
 2668         m(var(X) #\= integer(Y))             => [g(neq_num(X, Y))],
 2669         m(var(X) #\= var(Y))                 => [g(neq(X,Y))],
 2670         m(var(X) #\= var(Y) + var(Z))        => [p(x_neq_y_plus_z(X, Y, Z))],
 2671         m(var(X) #\= var(Y) - var(Z))        => [p(x_neq_y_plus_z(Y, X, Z))],
 2672         m(var(X) #\= var(Y)*var(Z))          => [p(ptimes(Y,Z,P)), g(neq(X,P))],
 2673         m(integer(X) #\= abs(any(Y)-any(Z))) => [d(Y, Y1), d(Z, Z1), p(absdiff_neq(Y1, Z1, X))],
 2674         m_c(any(X) #\= any(Y), left_right_linsum_const(X, Y, Cs, Vs, S)) =>
 2675            [g(scalar_product_(#\=, Cs, Vs, S))],
 2676         m(any(X) #\= any(Y) + any(Z))        => [d(X, X1), d(Y, Y1), d(Z, Z1), p(x_neq_y_plus_z(X1, Y1, Z1))],
 2677         m(any(X) #\= any(Y) - any(Z))        => [d(X, X1), d(Y, Y1), d(Z, Z1), p(x_neq_y_plus_z(Y1, X1, Z1))],
 2678         m(any(X) #\= any(Y)) => [d(X, RX), d(Y, RY), g(neq(RX, RY))]
 2679        ]).
 2680
 2681/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 2682   We again compile the committed-choice matching language to the
 2683   intended auxiliary predicates. We now must take care not to
 2684   unintentionally unify a variable with a complex term.
 2685- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 2686
 2687make_matches(Clauses) :-
 2688        matches(Ms),
 2689        findall(F, (member(M=>_, Ms), arg(1, M, M1), functor(M1, F, _)), Fs0),
 2690        sort(Fs0, Fs),
 2691        maplist(prevent_cyclic_argument, Fs, PrevCyclicClauses),
 2692        phrase(matchers(Ms), Clauses0),
 2693        maplist(goals_goal, Clauses0, MatcherClauses),
 2694        append(PrevCyclicClauses, MatcherClauses, Clauses1),
 2695        sort_by_predicate(Clauses1, Clauses).
 2696
 2697sort_by_predicate(Clauses, ByPred) :-
 2698        map_list_to_pairs(predname, Clauses, Keyed),
 2699        keysort(Keyed, KeyedByPred),
 2700        pairs_values(KeyedByPred, ByPred).
 2701
 2702predname((H:-_), Key)   :- !, predname(H, Key).
 2703predname(M:H, M:Key)    :- !, predname(H, Key).
 2704predname(H, Name/Arity) :- !, functor(H, Name, Arity).
 2705
 2706prevent_cyclic_argument(F0, Clause) :-
 2707        match_expand(F0, F),
 2708        Head =.. [F,X,Y],
 2709        Clause = (Head :- (   cyclic_term(X) ->
 2710                              domain_error(clpfd_expression, X)
 2711                          ;   cyclic_term(Y) ->
 2712                              domain_error(clpfd_expression, Y)
 2713                          ;   false
 2714                          )).
 2715
 2716matchers([]) --> [].
 2717matchers([Condition => Goals|Ms]) -->
 2718        matcher(Condition, Goals),
 2719        matchers(Ms).
 2720
 2721matcher(m(M), Gs) --> matcher(m_c(M,true), Gs).
 2722matcher(m_c(Matcher,Cond), Gs) -->
 2723        [(Head :- Goals0)],
 2724        { Matcher =.. [F,A,B],
 2725          match_expand(F, Expand),
 2726          Head =.. [Expand,X,Y],
 2727          phrase((match(A, X), match(B, Y)), Goals0, [Cond,!|Goals1]),
 2728          phrase(match_goals(Gs, Expand), Goals1) },
 2729        (   { symmetric(F), \+ (subsumes_term(A, B), subsumes_term(B, A)) } ->
 2730            { Head1 =.. [Expand,Y,X] },
 2731            [(Head1 :- Goals0)]
 2732        ;   []
 2733        ).
 2734
 2735match(any(A), T)     --> [A = T].
 2736match(var(V), T)     --> [( nonvar(T), ( T = ?(Var) ; T = #(Var) ) ->
 2737                            must_be_fd_integer(Var), V = Var
 2738                          ; v_or_i(T), V = T
 2739                          )].
 2740match(integer(I), T) --> [integer(T), I = T].
 2741match(-X, T)         --> [nonvar(T), T = -A], match(X, A).
 2742match(abs(X), T)     --> [nonvar(T), T = abs(A)], match(X, A).
 2743match(Binary, T)     -->
 2744        { Binary =.. [Op,X,Y], Term =.. [Op,A,B] },
 2745        [nonvar(T), T = Term],
 2746        match(X, A), match(Y, B).
 2747
 2748match_goals([], _)     --> [].
 2749match_goals([G|Gs], F) --> match_goal(G, F), match_goals(Gs, F).
 2750
 2751match_goal(r(X,Y), F)  --> { G =.. [F,X,Y] }, [G].
 2752match_goal(d(X,Y), _)  --> [parse_clpfd(X, Y)].
 2753match_goal(g(Goal), _) --> [Goal].
 2754match_goal(p(Prop), _) -->
 2755        [make_propagator(Prop, P)],
 2756        { term_variables(Prop, Vs) },
 2757        parse_init(Vs, P),
 2758        [trigger_once(P)].
 2759
 2760
 2761%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 2762
 2763
 2764
 2765%% ?X #>= ?Y
 2766%
 2767% Same as Y #=< X. When reasoning over integers, replace `(>=)/2` by
 2768% #>=/2 to obtain more general relations. See [declarative integer
 2769% arithmetic](<#clpfd-integer-arith>).
 2770
 2771X #>= Y :- clpfd_geq(X, Y).
 2772
 2773clpfd_geq(X, Y) :- clpfd_geq_(X, Y), reinforce(X), reinforce(Y).
 2774
 2775%% ?X #=< ?Y
 2776%
 2777% The arithmetic expression X is less than or equal to Y. When
 2778% reasoning over integers, replace `(=<)/2` by #=</2 to obtain more
 2779% general relations. See [declarative integer
 2780% arithmetic](<#clpfd-integer-arith>).
 2781
 2782X #=< Y :- Y #>= X.
 2783
 2784%% ?X #= ?Y
 2785%
 2786% The arithmetic expression X equals Y. This is the most important
 2787% [arithmetic constraint](<#clpfd-arith-constraints>), subsuming and
 2788% replacing both `(is)/2` _and_ `(=:=)/2` over integers. See
 2789% [declarative integer arithmetic](<#clpfd-integer-arith>).
 2790
 2791X #= Y :- clpfd_equal(X, Y).
 2792
 2793clpfd_equal(X, Y) :- clpfd_equal_(X, Y), reinforce(X).
 2794
 2795/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 2796   Conditions under which an equality can be compiled to built-in
 2797   arithmetic. Their order is significant. (/)/2 becomes (//)/2.
 2798- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 2799
 2800expr_conds(E, E)                 --> [integer(E)],
 2801        { var(E), !, \+ current_prolog_flag(clpfd_monotonic, true) }.
 2802expr_conds(E, E)                 --> { integer(E) }.
 2803expr_conds(?(E), E)              --> [integer(E)].
 2804expr_conds(#(E), E)              --> [integer(E)].
 2805expr_conds(-E0, -E)              --> expr_conds(E0, E).
 2806expr_conds(abs(E0), abs(E))      --> expr_conds(E0, E).
 2807expr_conds(A0+B0, A+B)           --> expr_conds(A0, A), expr_conds(B0, B).
 2808expr_conds(A0*B0, A*B)           --> expr_conds(A0, A), expr_conds(B0, B).
 2809expr_conds(A0-B0, A-B)           --> expr_conds(A0, A), expr_conds(B0, B).
 2810expr_conds(A0//B0, A//B)         -->
 2811        expr_conds(A0, A), expr_conds(B0, B),
 2812        [B =\= 0].
 2813%expr_conds(A0/B0, AB)            --> expr_conds(A0//B0, AB).
 2814expr_conds(min(A0,B0), min(A,B)) --> expr_conds(A0, A), expr_conds(B0, B).
 2815expr_conds(max(A0,B0), max(A,B)) --> expr_conds(A0, A), expr_conds(B0, B).
 2816expr_conds(A0 mod B0, A mod B)   -->
 2817        expr_conds(A0, A), expr_conds(B0, B),
 2818        [B =\= 0].
 2819expr_conds(A0^B0, A^B)           -->
 2820        expr_conds(A0, A), expr_conds(B0, B),
 2821        [(B >= 0 ; A =:= -1)].
 2822% Bitwise operations, added to make CLP(FD) usable in more cases
 2823expr_conds(\ A0, \ A) --> expr_conds(A0, A).
 2824expr_conds(A0<<B0, A<<B) --> expr_conds(A0, A), expr_conds(B0, B).
 2825expr_conds(A0>>B0, A>>B) --> expr_conds(A0, A), expr_conds(B0, B).
 2826expr_conds(A0/\B0, A/\B) --> expr_conds(A0, A), expr_conds(B0, B).
 2827expr_conds(A0\/B0, A\/B) --> expr_conds(A0, A), expr_conds(B0, B).
 2828expr_conds(A0 xor B0, A xor B) --> expr_conds(A0, A), expr_conds(B0, B).
 2829expr_conds(lsb(A0), lsb(A)) --> expr_conds(A0, A).
 2830expr_conds(msb(A0), msb(A)) --> expr_conds(A0, A).
 2831expr_conds(popcount(A0), popcount(A)) --> expr_conds(A0, A).
 2832
 2833:- multifile
 2834        system:goal_expansion/2. 2835:- dynamic
 2836        system:goal_expansion/2. 2837
 2838system:goal_expansion(Goal, Expansion) :-
 2839        \+ current_prolog_flag(clpfd_goal_expansion, false),
 2840        clpfd_expandable(Goal),
 2841        prolog_load_context(module, M),
 2842	(   M == clpfd
 2843	->  true
 2844	;   predicate_property(M:Goal, imported_from(clpfd))
 2845	),
 2846        clpfd_expansion(Goal, Expansion).
 2847
 2848clpfd_expandable(_ in _).
 2849clpfd_expandable(_ #= _).
 2850clpfd_expandable(_ #>= _).
 2851clpfd_expandable(_ #=< _).
 2852clpfd_expandable(_ #> _).
 2853clpfd_expandable(_ #< _).
 2854
 2855clpfd_expansion(Var in Dom, In) :-
 2856        (   ground(Dom), Dom = L..U, integer(L), integer(U) ->
 2857            expansion_simpler(
 2858                (   integer(Var) ->
 2859                    between(L, U, Var)
 2860                ;   clpfd:clpfd_in(Var, Dom)
 2861                ), In)
 2862        ;   In = clpfd:clpfd_in(Var, Dom)
 2863        ).
 2864clpfd_expansion(X0 #= Y0, Equal) :-
 2865        phrase(expr_conds(X0, X), CsX),
 2866        phrase(expr_conds(Y0, Y), CsY),
 2867        list_goal(CsX, CondX),
 2868        list_goal(CsY, CondY),
 2869        expansion_simpler(
 2870                (   CondX ->
 2871                    (   var(Y) -> Y is X
 2872                    ;   CondY -> X =:= Y
 2873                    ;   T is X, clpfd:clpfd_equal(T, Y0)
 2874                    )
 2875                ;   CondY ->
 2876                    (   var(X) -> X is Y
 2877                    ;   T is Y, clpfd:clpfd_equal(X0, T)
 2878                    )
 2879                ;   clpfd:clpfd_equal(X0, Y0)
 2880                ), Equal).
 2881clpfd_expansion(X0 #>= Y0, Geq) :-
 2882        phrase(expr_conds(X0, X), CsX),
 2883        phrase(expr_conds(Y0, Y), CsY),
 2884        list_goal(CsX, CondX),
 2885        list_goal(CsY, CondY),
 2886        expansion_simpler(
 2887              (   CondX ->
 2888                  (   CondY -> X >= Y
 2889                  ;   T is X, clpfd:clpfd_geq(T, Y0)
 2890                  )
 2891              ;   CondY -> T is Y, clpfd:clpfd_geq(X0, T)
 2892              ;   clpfd:clpfd_geq(X0, Y0)
 2893              ), Geq).
 2894clpfd_expansion(X #=< Y,  Leq) :- clpfd_expansion(Y #>= X, Leq).
 2895clpfd_expansion(X #> Y, Gt)    :- clpfd_expansion(X #>= Y+1, Gt).
 2896clpfd_expansion(X #< Y, Lt)    :- clpfd_expansion(Y #> X, Lt).
 2897
 2898expansion_simpler((True->Then0;_), Then) :-
 2899        is_true(True), !,
 2900        expansion_simpler(Then0, Then).
 2901expansion_simpler((False->_;Else0), Else) :-
 2902        is_false(False), !,
 2903        expansion_simpler(Else0, Else).
 2904expansion_simpler((If->Then0;Else0), (If->Then;Else)) :- !,
 2905        expansion_simpler(Then0, Then),
 2906        expansion_simpler(Else0, Else).
 2907expansion_simpler((A0,B0), (A,B)) :-
 2908        expansion_simpler(A0, A),
 2909        expansion_simpler(B0, B).
 2910expansion_simpler(Var is Expr0, Goal) :-
 2911        ground(Expr0), !,
 2912        phrase(expr_conds(Expr0, Expr), Gs),
 2913        (   maplist(call, Gs) -> Value is Expr, Goal = (Var = Value)
 2914        ;   Goal = false
 2915        ).
 2916expansion_simpler(Var =:= Expr0, Goal) :-
 2917        ground(Expr0), !,
 2918        phrase(expr_conds(Expr0, Expr), Gs),
 2919        (   maplist(call, Gs) -> Value is Expr, Goal = (Var =:= Value)
 2920        ;   Goal = false
 2921        ).
 2922expansion_simpler(Var is Expr, Var = Expr) :- var(Expr), !.
 2923expansion_simpler(Var is Expr, Goal) :- !,
 2924        (   var(Var), nonvar(Expr),
 2925            Expr = E mod M, nonvar(E), E = A^B ->
 2926            Goal = ( ( integer(A), integer(B), integer(M),
 2927                       A >= 0, B >= 0, M > 0 ->
 2928                       Var is powm(A, B, M)
 2929                     ; Var is Expr
 2930                     ) )
 2931        ;   Goal = ( Var is Expr )
 2932        ).
 2933expansion_simpler(between(L,U,V), Goal) :- maplist(integer, [L,U,V]), !,
 2934        (   between(L,U,V) -> Goal = true
 2935        ;   Goal = false
 2936        ).
 2937expansion_simpler(Goal, Goal).
 2938
 2939is_true(true).
 2940is_true(integer(I))  :- integer(I).
 2941:- if(current_predicate(var_property/2)). 2942is_true(var(X))      :- var(X), var_property(X, fresh(true)).
 2943is_false(integer(X)) :- var(X), var_property(X, fresh(true)).
 2944is_false((A,B))      :- is_false(A) ; is_false(B).
 2945:- endif. 2946is_false(var(X)) :- nonvar(X).
 2947
 2948
 2949%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 2950
 2951linsum(X, S, S)    --> { var(X), !, non_monotonic(X) }, [vn(X,1)].
 2952linsum(I, S0, S)   --> { integer(I), S is S0 + I }.
 2953linsum(?(X), S, S) --> { must_be_fd_integer(X) }, [vn(X,1)].
 2954linsum(#(X), S, S) --> { must_be_fd_integer(X) }, [vn(X,1)].
 2955linsum(-A, S0, S)  --> mulsum(A, -1, S0, S).
 2956linsum(N*A, S0, S) --> { integer(N) }, !, mulsum(A, N, S0, S).
 2957linsum(A*N, S0, S) --> { integer(N) }, !, mulsum(A, N, S0, S).
 2958linsum(A+B, S0, S) --> linsum(A, S0, S1), linsum(B, S1, S).
 2959linsum(A-B, S0, S) --> linsum(A, S0, S1), mulsum(B, -1, S1, S).
 2960
 2961mulsum(A, M, S0, S) -->
 2962        { phrase(linsum(A, 0, CA), As), S is S0 + M*CA },
 2963        lin_mul(As, M).
 2964
 2965lin_mul([], _)             --> [].
 2966lin_mul([vn(X,N0)|VNs], M) --> { N is N0*M }, [vn(X,N)], lin_mul(VNs, M).
 2967
 2968v_or_i(V) :- var(V), !, non_monotonic(V).
 2969v_or_i(I) :- integer(I).
 2970
 2971must_be_fd_integer(X) :-
 2972        (   var(X) -> constrain_to_integer(X)
 2973        ;   must_be(integer, X)
 2974        ).
 2975
 2976left_right_linsum_const(Left, Right, Cs, Vs, Const) :-
 2977        phrase(linsum(Left, 0, CL), Lefts0, Rights),
 2978        phrase(linsum(Right, 0, CR), Rights0),
 2979        maplist(linterm_negate, Rights0, Rights),
 2980        msort(Lefts0, Lefts),
 2981        Lefts = [vn(First,N)|LeftsRest],
 2982        vns_coeffs_variables(LeftsRest, N, First, Cs0, Vs0),
 2983        filter_linsum(Cs0, Vs0, Cs, Vs),
 2984        Const is CR - CL.
 2985        %format("linear sum: ~w ~w ~w\n", [Cs,Vs,Const]).
 2986
 2987linterm_negate(vn(V,N0), vn(V,N)) :- N is -N0.
 2988
 2989vns_coeffs_variables([], N, V, [N], [V]).
 2990vns_coeffs_variables([vn(V,N)|VNs], N0, V0, Ns, Vs) :-
 2991        (   V == V0 ->
 2992            N1 is N0 + N,
 2993            vns_coeffs_variables(VNs, N1, V0, Ns, Vs)
 2994        ;   Ns = [N0|NRest],
 2995            Vs = [V0|VRest],
 2996            vns_coeffs_variables(VNs, N, V, NRest, VRest)
 2997        ).
 2998
 2999filter_linsum([], [], [], []).
 3000filter_linsum([C0|Cs0], [V0|Vs0], Cs, Vs) :-
 3001        (   C0 =:= 0 ->
 3002            constrain_to_integer(V0),
 3003            filter_linsum(Cs0, Vs0, Cs, Vs)
 3004        ;   Cs = [C0|Cs1], Vs = [V0|Vs1],
 3005            filter_linsum(Cs0, Vs0, Cs1, Vs1)
 3006        ).
 3007
 3008gcd([], G, G).
 3009gcd([N|Ns], G0, G) :-
 3010        G1 is gcd(N, G0),
 3011        gcd(Ns, G1, G).
 3012
 3013even(N) :- N mod 2 =:= 0.
 3014
 3015odd(N) :- \+ even(N).
 3016
 3017/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 3018   k-th root of N, if N is a k-th power.
 3019
 3020   TODO: Replace this when the GMP function becomes available.
 3021- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 3022
 3023integer_kth_root(N, K, R) :-
 3024        (   even(K) ->
 3025            N >= 0
 3026        ;   true
 3027        ),
 3028        (   N < 0 ->
 3029            odd(K),
 3030            integer_kroot(N, 0, N, K, R)
 3031        ;   integer_kroot(0, N, N, K, R)
 3032        ).
 3033
 3034integer_kroot(L, U, N, K, R) :-
 3035        (   L =:= U -> N =:= L^K, R = L
 3036        ;   L + 1 =:= U ->
 3037            (   L^K =:= N -> R = L
 3038            ;   U^K =:= N -> R = U
 3039            ;   false
 3040            )
 3041        ;   Mid is (L + U)//2,
 3042            (   Mid^K > N ->
 3043                integer_kroot(L, Mid, N, K, R)
 3044            ;   integer_kroot(Mid, U, N, K, R)
 3045            )
 3046        ).
 3047
 3048integer_log_b(N, B, Log0, Log) :-
 3049        T is B^Log0,
 3050        (   T =:= N -> Log = Log0
 3051        ;   T < N,
 3052            Log1 is Log0 + 1,
 3053            integer_log_b(N, B, Log1, Log)
 3054        ).
 3055
 3056floor_integer_log_b(N, B, Log0, Log) :-
 3057        T is B^Log0,
 3058        (   T > N -> Log is Log0 - 1
 3059        ;   T =:= N -> Log = Log0
 3060        ;   T < N,
 3061            Log1 is Log0 + 1,
 3062            floor_integer_log_b(N, B, Log1, Log)
 3063        ).
 3064
 3065
 3066/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 3067   Largest R such that R^K =< N.
 3068- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 3069
 3070:- if(current_predicate(nth_integer_root_and_remainder/4)). 3071
 3072% This currently only works for K >= 1, which is all that is needed for now.
 3073integer_kth_root_leq(N, K, R) :-
 3074        nth_integer_root_and_remainder(K, N, R, _).
 3075
 3076:- else. 3077
 3078integer_kth_root_leq(N, K, R) :-
 3079        (   even(K) ->
 3080            N >= 0
 3081        ;   true
 3082        ),
 3083        (   N < 0 ->
 3084            odd(K),
 3085            integer_kroot_leq(N, 0, N, K, R)
 3086        ;   integer_kroot_leq(0, N, N, K, R)
 3087        ).
 3088
 3089integer_kroot_leq(L, U, N, K, R) :-
 3090        (   L =:= U -> R = L
 3091        ;   L + 1 =:= U ->
 3092            (   U^K =< N -> R = U
 3093            ;   R = L
 3094            )
 3095        ;   Mid is (L + U)//2,
 3096            (   Mid^K > N ->
 3097                integer_kroot_leq(L, Mid, N, K, R)
 3098            ;   integer_kroot_leq(Mid, U, N, K, R)
 3099            )
 3100        ).
 3101
 3102:- endif. 3103
 3104%% ?X #\= ?Y
 3105%
 3106% The arithmetic expressions X and Y evaluate to distinct integers.
 3107% When reasoning over integers, replace `(=\=)/2` by #\=/2 to obtain
 3108% more general relations. See [declarative integer
 3109% arithmetic](<#clpfd-integer-arith>).
 3110
 3111X #\= Y :- clpfd_neq(X, Y), do_queue.
 3112
 3113% X #\= Y + Z
 3114
 3115x_neq_y_plus_z(X, Y, Z) :-
 3116        propagator_init_trigger(x_neq_y_plus_z(X,Y,Z)).
 3117
 3118% X is distinct from the number N. This is used internally, and does
 3119% not reinforce other constraints.
 3120
 3121neq_num(X, N) :-
 3122        (   fd_get(X, XD, XPs) ->
 3123            domain_remove(XD, N, XD1),
 3124            fd_put(X, XD1, XPs)
 3125        ;   X =\= N
 3126        ).
 3127
 3128%% ?X #> ?Y
 3129%
 3130% Same as Y #< X. When reasoning over integers, replace `(>)/2` by
 3131% #>/2 to obtain more general relations See [declarative integer
 3132% arithmetic](<#clpfd-integer-arith>).
 3133
 3134X #> Y  :- X #>= Y + 1.
 3135
 3136%% #<(?X, ?Y)
 3137%
 3138% The arithmetic expression X is less than Y. When reasoning over
 3139% integers, replace `(<)/2` by #</2 to obtain more general relations. See
 3140% [declarative integer arithmetic](<#clpfd-integer-arith>).
 3141%
 3142% In addition to its regular use in tasks that require it, this
 3143% constraint can also be useful to eliminate uninteresting symmetries
 3144% from a problem. For example, all possible matches between pairs
 3145% built from four players in total:
 3146%
 3147% ==
 3148% ?- Vs = [A,B,C,D], Vs ins 1..4,
 3149%         all_different(Vs),
 3150%         A #< B, C #< D, A #< C,
 3151%    findall(pair(A,B)-pair(C,D), label(Vs), Ms).
 3152% Ms = [ pair(1, 2)-pair(3, 4),
 3153%        pair(1, 3)-pair(2, 4),
 3154%        pair(1, 4)-pair(2, 3)].
 3155% ==
 3156
 3157X #< Y  :- Y #> X.
 3158
 3159%% #\ +Q
 3160%
 3161% Q does _not_ hold. See [reification](<#clpfd-reification>).
 3162%
 3163% For example, to obtain the complement of a domain:
 3164%
 3165% ==
 3166% ?- #\ X in -3..0\/10..80.
 3167% X in inf.. -4\/1..9\/81..sup.
 3168% ==
 3169
 3170#\ Q       :- reify(Q, 0), do_queue.
 3171
 3172%% ?P #<==> ?Q
 3173%
 3174% P and Q are equivalent. See [reification](<#clpfd-reification>).
 3175%
 3176% For example:
 3177%
 3178% ==
 3179% ?- X #= 4 #<==> B, X #\= 4.
 3180% B = 0,
 3181% X in inf..3\/5..sup.
 3182% ==
 3183% The following example uses reified constraints to relate a list of
 3184% finite domain variables to the number of occurrences of a given value:
 3185%
 3186% ==
 3187% vs_n_num(Vs, N, Num) :-
 3188%         maplist(eq_b(N), Vs, Bs),
 3189%         sum(Bs, #=, Num).
 3190%
 3191% eq_b(X, Y, B) :- X #= Y #<==> B.
 3192% ==
 3193%
 3194% Sample queries and their results:
 3195%
 3196% ==
 3197% ?- Vs = [X,Y,Z], Vs ins 0..1, vs_n_num(Vs, 4, Num).
 3198% Vs = [X, Y, Z],
 3199% Num = 0,
 3200% X in 0..1,
 3201% Y in 0..1,
 3202% Z in 0..1.
 3203%
 3204% ?- vs_n_num([X,Y,Z], 2, 3).
 3205% X = 2,
 3206% Y = 2,
 3207% Z = 2.
 3208% ==
 3209
 3210L #<==> R  :- reify(L, B), reify(R, B), do_queue.
 3211
 3212%% ?P #==> ?Q
 3213%
 3214% P implies Q. See [reification](<#clpfd-reification>).
 3215
 3216/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 3217   Implication is special in that created auxiliary constraints can be
 3218   retracted when the implication becomes entailed, for example:
 3219
 3220   %?- X + 1 #= Y #==> Z, Z #= 1.
 3221   %@ Z = 1,
 3222   %@ X in inf..sup,
 3223   %@ Y in inf..sup.
 3224
 3225   We cannot use propagator_init_trigger/1 here because the states of
 3226   auxiliary propagators are themselves part of the propagator.
 3227- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 3228
 3229L #==> R   :-
 3230        reify(L, LB, LPs),
 3231        reify(R, RB, RPs),
 3232        append(LPs, RPs, Ps),
 3233        propagator_init_trigger([LB,RB], pimpl(LB,RB,Ps)).
 3234
 3235%% ?P #<== ?Q
 3236%
 3237% Q implies P. See [reification](<#clpfd-reification>).
 3238
 3239L #<== R   :- R #==> L.
 3240
 3241%% ?P #/\ ?Q
 3242%
 3243% P and Q hold. See [reification](<#clpfd-reification>).
 3244
 3245L #/\ R    :- reify(L, 1), reify(R, 1), do_queue.
 3246
 3247conjunctive_neqs_var_drep(Eqs, Var, Drep) :-
 3248        conjunctive_neqs_var(Eqs, Var),
 3249        phrase(conjunctive_neqs_vals(Eqs), Vals),
 3250        list_to_domain(Vals, Dom),
 3251        domain_complement(Dom, C),
 3252        domain_to_drep(C, Drep).
 3253
 3254conjunctive_neqs_var(V, _) :- var(V), !, false.
 3255conjunctive_neqs_var(L #\= R, Var) :-
 3256        (   var(L), integer(R) -> Var = L
 3257        ;   integer(L), var(R) -> Var = R
 3258        ;   false
 3259        ).
 3260conjunctive_neqs_var(A #/\ B, VA) :-
 3261        conjunctive_neqs_var(A, VA),
 3262        conjunctive_neqs_var(B, VB),
 3263        VA == VB.
 3264
 3265conjunctive_neqs_vals(L #\= R) --> ( { integer(L) } -> [L] ; [R] ).
 3266conjunctive_neqs_vals(A #/\ B) -->
 3267        conjunctive_neqs_vals(A),
 3268        conjunctive_neqs_vals(B).
 3269
 3270%% ?P #\/ ?Q
 3271%
 3272% P or Q holds. See [reification](<#clpfd-reification>).
 3273%
 3274% For example, the sum of natural numbers below 1000 that are
 3275% multiples of 3 or 5:
 3276%
 3277% ==
 3278% ?- findall(N, (N mod 3 #= 0 #\/ N mod 5 #= 0, N in 0..999,
 3279%                indomain(N)),
 3280%            Ns),
 3281%    sum(Ns, #=, Sum).
 3282% Ns = [0, 3, 5, 6, 9, 10, 12, 15, 18|...],
 3283% Sum = 233168.
 3284% ==
 3285
 3286L #\/ R :-
 3287        (   disjunctive_eqs_var_drep(L #\/ R, Var, Drep) -> Var in Drep
 3288        ;   reify(L, X, Ps1),
 3289            reify(R, Y, Ps2),
 3290            propagator_init_trigger([X,Y], reified_or(X,Ps1,Y,Ps2,1))
 3291        ).
 3292
 3293disjunctive_eqs_var_drep(Eqs, Var, Drep) :-
 3294        disjunctive_eqs_var(Eqs, Var),
 3295        phrase(disjunctive_eqs_vals(Eqs), Vals),
 3296        list_to_drep(Vals, Drep).
 3297
 3298disjunctive_eqs_var(V, _) :- var(V), !, false.
 3299disjunctive_eqs_var(V in I, V) :- var(V), integer(I).
 3300disjunctive_eqs_var(L #= R, Var) :-
 3301        (   var(L), integer(R) -> Var = L
 3302        ;   integer(L), var(R) -> Var = R
 3303        ;   false
 3304        ).
 3305disjunctive_eqs_var(A #\/ B, VA) :-
 3306        disjunctive_eqs_var(A, VA),
 3307        disjunctive_eqs_var(B, VB),
 3308        VA == VB.
 3309
 3310disjunctive_eqs_vals(L #= R)  --> ( { integer(L) } -> [L] ; [R] ).
 3311disjunctive_eqs_vals(_ in I)  --> [I].
 3312disjunctive_eqs_vals(A #\/ B) -->
 3313        disjunctive_eqs_vals(A),
 3314        disjunctive_eqs_vals(B).
 3315
 3316%% ?P #\ ?Q
 3317%
 3318% Either P holds or Q holds, but not both. See
 3319% [reification](<#clpfd-reification>).
 3320
 3321L #\ R :- (L #\/ R) #/\ #\ (L #/\ R).
 3322
 3323/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 3324   A constraint that is being reified need not hold. Therefore, in
 3325   X/Y, Y can as well be 0, for example. Note that it is OK to
 3326   constrain the *result* of an expression (which does not appear
 3327   explicitly in the expression and is not visible to the outside),
 3328   but not the operands, except for requiring that they be integers.
 3329
 3330   In contrast to parse_clpfd/2, the result of an expression can now
 3331   also be undefined, in which case the constraint cannot hold.
 3332   Therefore, the committed-choice language is extended by an element
 3333   d(D) that states D is 1 iff all subexpressions are defined. a(V)
 3334   means that V is an auxiliary variable that was introduced while
 3335   parsing a compound expression. a(X,V) means V is auxiliary unless
 3336   it is ==/2 X, and a(X,Y,V) means V is auxiliary unless it is ==/2 X
 3337   or Y. l(L) means the literal L occurs in the described list.
 3338
 3339   When a constraint becomes entailed or subexpressions become
 3340   undefined, created auxiliary constraints are killed, and the
 3341   "clpfd" attribute is removed from auxiliary variables.
 3342
 3343   For (/)/2, mod/2 and rem/2, we create a skeleton propagator and
 3344   remember it as an auxiliary constraint. The pskeleton propagator
 3345   can use the skeleton when the constraint is defined.
 3346- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 3347
 3348parse_reified(E, R, D,
 3349              [g(cyclic_term(E)) => [g(domain_error(clpfd_expression, E))],
 3350               g(var(E))     => [g(non_monotonic(E)),
 3351                                 g(constrain_to_integer(E)), g(R = E), g(D=1)],
 3352               g(integer(E)) => [g(R=E), g(D=1)],
 3353               ?(E)          => [g(must_be_fd_integer(E)), g(R=E), g(D=1)],
 3354               #(E)          => [g(must_be_fd_integer(E)), g(R=E), g(D=1)],
 3355               m(A+B)        => [d(D), p(pplus(A,B,R)), a(A,B,R)],
 3356               m(A*B)        => [d(D), p(ptimes(A,B,R)), a(A,B,R)],
 3357               m(A-B)        => [d(D), p(pplus(R,B,A)), a(A,B,R)],
 3358               m(-A)         => [d(D), p(ptimes(-1,A,R)), a(R)],
 3359               m(max(A,B))   => [d(D), p(pgeq(R, A)), p(pgeq(R, B)), p(pmax(A,B,R)), a(A,B,R)],
 3360               m(min(A,B))   => [d(D), p(pgeq(A, R)), p(pgeq(B, R)), p(pmin(A,B,R)), a(A,B,R)],
 3361               m(abs(A))     => [g(?(R)#>=0), d(D), p(pabs(A, R)), a(A,R)],
 3362%               m(A/B)        => [skeleton(A,B,D,R,ptzdiv)],
 3363               m(A//B)       => [skeleton(A,B,D,R,ptzdiv)],
 3364               m(A div B)    => [skeleton(A,B,D,R,pdiv)],
 3365               m(A rdiv B)   => [skeleton(A,B,D,R,prdiv)],
 3366               m(A mod B)    => [skeleton(A,B,D,R,pmod)],
 3367               m(A rem B)    => [skeleton(A,B,D,R,prem)],
 3368               m(A^B)        => [d(D), p(pexp(A,B,R)), a(A,B,R)],
 3369               % bitwise operations
 3370               m(\A)         => [function(D,\,A,R)],
 3371               m(msb(A))     => [function(D,msb,A,R)],
 3372               m(lsb(A))     => [function(D,lsb,A,R)],
 3373               m(popcount(A)) => [function(D,popcount,A,R)],
 3374               m(A<<B)       => [function(D,<<,A,B,R)],
 3375               m(A>>B)       => [function(D,>>,A,B,R)],
 3376               m(A/\B)       => [function(D,/\,A,B,R)],
 3377               m(A\/B)       => [function(D,\/,A,B,R)],
 3378               m(A xor B)    => [function(D,xor,A,B,R)],
 3379               g(true)       => [g(domain_error(clpfd_expression, E))]]
 3380             ).
 3381
 3382% Again, we compile this to a predicate, parse_reified_clpfd//3. This
 3383% time, it is a DCG that describes the list of auxiliary variables and
 3384% propagators for the given expression, in addition to relating it to
 3385% its reified (Boolean) finite domain variable and its Boolean
 3386% definedness.
 3387
 3388make_parse_reified(Clauses) :-
 3389        parse_reified_clauses(Clauses0),
 3390        maplist(goals_goal_dcg, Clauses0, Clauses).
 3391
 3392goals_goal_dcg((Head --> Goals), Clause) :-
 3393        list_goal(Goals, Body),
 3394        expand_term((Head --> Body), Clause).
 3395
 3396parse_reified_clauses(Clauses) :-
 3397        parse_reified(E, R, D, Matchers),
 3398        maplist(parse_reified(E, R, D), Matchers, Clauses).
 3399
 3400parse_reified(E, R, D, Matcher, Clause) :-
 3401        Matcher = (Condition0 => Goals0),
 3402        phrase((reified_condition(Condition0, E, Head, Ds),
 3403                reified_goals(Goals0, Ds)), Goals, [a(D)]),
 3404        Clause = (parse_reified_clpfd(Head, R, D) --> Goals).
 3405
 3406reified_condition(g(Goal), E, E, []) --> [{Goal}, !].
 3407reified_condition(?(E), _, ?(E), []) --> [!].
 3408reified_condition(#(E), _, #(E), []) --> [!].
 3409reified_condition(m(Match), _, Match0, Ds) -->
 3410        [!],
 3411        { copy_term(Match, Match0),
 3412          term_variables(Match0, Vs0),
 3413          term_variables(Match, Vs)
 3414        },
 3415        reified_variables(Vs0, Vs, Ds).
 3416
 3417reified_variables([], [], []) --> [].
 3418reified_variables([V0|Vs0], [V|Vs], [D|Ds]) -->
 3419        [parse_reified_clpfd(V0, V, D)],
 3420        reified_variables(Vs0, Vs, Ds).
 3421
 3422reified_goals([], _) --> [].
 3423reified_goals([G|Gs], Ds) --> reified_goal(G, Ds), reified_goals(Gs, Ds).
 3424
 3425reified_goal(d(D), Ds) -->
 3426        (   { Ds = [X] } -> [{D=X}]
 3427        ;   { Ds = [X,Y] } ->
 3428            { phrase(reified_goal(p(reified_and(X,[],Y,[],D)), _), Gs),
 3429              list_goal(Gs, Goal) },
 3430            [( {X==1, Y==1} -> {D = 1} ; Goal )]
 3431        ;   { domain_error(one_or_two_element_list, Ds) }
 3432        ).
 3433reified_goal(g(Goal), _) --> [{Goal}].
 3434reified_goal(p(Vs, Prop), _) -->
 3435        [{make_propagator(Prop, P)}],
 3436        parse_init_dcg(Vs, P),
 3437        [{trigger_once(P)}],
 3438        [( { propagator_state(P, S), S == dead } -> [] ; [p(P)])].
 3439reified_goal(p(Prop), Ds) -->
 3440        { term_variables(Prop, Vs) },
 3441        reified_goal(p(Vs,Prop), Ds).
 3442reified_goal(function(D,Op,A,B,R), Ds) -->
 3443        reified_goals([d(D),p(pfunction(Op,A,B,R)),a(A,B,R)], Ds).
 3444reified_goal(function(D,Op,A,R), Ds) -->
 3445        reified_goals([d(D),p(pfunction(Op,A,R)),a(A,R)], Ds).
 3446reified_goal(skeleton(A,B,D,R,F), Ds) -->
 3447        { Prop =.. [F,X,Y,Z] },
 3448        reified_goals([d(D1),l(p(P)),g(make_propagator(Prop, P)),
 3449                       p([A,B,D2,R], pskeleton(A,B,D2,[X,Y,Z]-P,R,F)),
 3450                       p(reified_and(D1,[],D2,[],D)),a(D2),a(A,B,R)], Ds).
 3451reified_goal(a(V), _)     --> [a(V)].
 3452reified_goal(a(X,V), _)   --> [a(X,V)].
 3453reified_goal(a(X,Y,V), _) --> [a(X,Y,V)].
 3454reified_goal(l(L), _)     --> [[L]].
 3455
 3456parse_init_dcg([], _)     --> [].
 3457parse_init_dcg([V|Vs], P) --> [{init_propagator(V, P)}], parse_init_dcg(Vs, P).
 3458
 3459%?- set_prolog_flag(answer_write_options, [portray(true)]),
 3460%   clpfd:parse_reified_clauses(Cs), maplist(portray_clause, Cs).
 3461
 3462reify(E, B) :- reify(E, B, _).
 3463
 3464reify(Expr, B, Ps) :-
 3465        (   acyclic_term(Expr), reifiable(Expr) -> phrase(reify(Expr, B), Ps)
 3466        ;   domain_error(clpfd_reifiable_expression, Expr)
 3467        ).
 3468
 3469reifiable(E)      :- var(E), non_monotonic(E).
 3470reifiable(E)      :- integer(E), E in 0..1.
 3471reifiable(?(E))   :- must_be_fd_integer(E).
 3472reifiable(#(E))   :- must_be_fd_integer(E).
 3473reifiable(V in _) :- fd_variable(V).
 3474reifiable(Expr)   :-
 3475        Expr =.. [Op,Left,Right],
 3476        (   memberchk(Op, [#>=,#>,#=<,#<,#=,#\=])
 3477        ;   memberchk(Op, [#==>,#<==,#<==>,#/\,#\/,#\]),
 3478            reifiable(Left),
 3479            reifiable(Right)
 3480        ).
 3481reifiable(#\ E) :- reifiable(E).
 3482reifiable(tuples_in(Tuples, Relation)) :-
 3483        must_be(list(list), Tuples),
 3484        maplist(maplist(fd_variable), Tuples),
 3485        must_be(list(list(integer)), Relation).
 3486reifiable(finite_domain(V)) :- fd_variable(V).
 3487
 3488reify(E, B) --> { B in 0..1 }, reify_(E, B).
 3489
 3490reify_(E, B) --> { var(E), !, E = B }.
 3491reify_(E, B) --> { integer(E), E = B }.
 3492reify_(?(B), B) --> [].
 3493reify_(#(B), B) --> [].
 3494reify_(V in Drep, B) -->
 3495        { drep_to_domain(Drep, Dom) },
 3496        propagator_init_trigger(reified_in(V,Dom,B)),
 3497        a(B).
 3498reify_(tuples_in(Tuples, Relation), B) -->
 3499        { maplist(relation_tuple_b_prop(Relation), Tuples, Bs, Ps),
 3500          maplist(monotonic, Bs, Bs1),
 3501          fold_statement(conjunction, Bs1, And),
 3502          ?(B) #<==> And },
 3503        propagator_init_trigger([B], tuples_not_in(Tuples, Relation, B)),
 3504        kill_reified_tuples(Bs, Ps, Bs),
 3505        list(Ps),
 3506        as([B|Bs]).
 3507reify_(finite_domain(V), B) -->
 3508        propagator_init_trigger(reified_fd(V,B)),
 3509        a(B).
 3510reify_(L #>= R, B) --> arithmetic(L, R, B, reified_geq).
 3511reify_(L #= R, B)  --> arithmetic(L, R, B, reified_eq).
 3512reify_(L #\= R, B) --> arithmetic(L, R, B, reified_neq).
 3513reify_(L #> R, B)  --> reify_(L #>= (R+1), B).
 3514reify_(L #=< R, B) --> reify_(R #>= L, B).
 3515reify_(L #< R, B)  --> reify_(R #>= (L+1), B).
 3516reify_(L #==> R, B)  --> reify_((#\ L) #\/ R, B).
 3517reify_(L #<== R, B)  --> reify_(R #==> L, B).
 3518reify_(L #<==> R, B) --> reify_((L #==> R) #/\ (R #==> L), B).
 3519reify_(L #\ R, B) --> reify_((L #\/ R) #/\ #\ (L #/\ R), B).
 3520reify_(L #/\ R, B)   -->
 3521        (   { conjunctive_neqs_var_drep(L #/\ R, V, D) } -> reify_(V in D, B)
 3522        ;   boolean(L, R, B, reified_and)
 3523        ).
 3524reify_(L #\/ R, B) -->
 3525        (   { disjunctive_eqs_var_drep(L #\/ R, V, D) } -> reify_(V in D, B)
 3526        ;   boolean(L, R, B, reified_or)
 3527        ).
 3528reify_(#\ Q, B) -->
 3529        reify(Q, QR),
 3530        propagator_init_trigger(reified_not(QR,B)),
 3531        a(B).
 3532
 3533arithmetic(L, R, B, Functor) -->
 3534        { phrase((parse_reified_clpfd(L, LR, LD),
 3535                  parse_reified_clpfd(R, RR, RD)), Ps),
 3536          Prop =.. [Functor,LD,LR,RD,RR,Ps,B] },
 3537        list(Ps),
 3538        propagator_init_trigger([LD,LR,RD,RR,B], Prop),
 3539        a(B).
 3540
 3541boolean(L, R, B, Functor) -->
 3542        { reify(L, LR, Ps1), reify(R, RR, Ps2),
 3543          Prop =.. [Functor,LR,Ps1,RR,Ps2,B] },
 3544        list(Ps1), list(Ps2),
 3545        propagator_init_trigger([LR,RR,B], Prop),
 3546        a(LR, RR, B).
 3547
 3548list([])     --> [].
 3549list([L|Ls]) --> [L], list(Ls).
 3550
 3551a(X,Y,B) -->
 3552        (   { nonvar(X) } -> a(Y, B)
 3553        ;   { nonvar(Y) } -> a(X, B)
 3554        ;   [a(X,Y,B)]
 3555        ).
 3556
 3557a(X, B) -->
 3558        (   { var(X) } -> [a(X, B)]
 3559        ;   a(B)
 3560        ).
 3561
 3562a(B) -->
 3563        (   { var(B) } -> [a(B)]
 3564        ;   []
 3565        ).
 3566
 3567as([])     --> [].
 3568as([B|Bs]) --> a(B), as(Bs).
 3569
 3570kill_reified_tuples([], _, _) --> [].
 3571kill_reified_tuples([B|Bs], Ps, All) -->
 3572        propagator_init_trigger([B], kill_reified_tuples(B, Ps, All)),
 3573        kill_reified_tuples(Bs, Ps, All).
 3574
 3575relation_tuple_b_prop(Relation, Tuple, B, p(Prop)) :-
 3576        put_attr(R, clpfd_relation, Relation),
 3577        make_propagator(reified_tuple_in(Tuple, R, B), Prop),
 3578        tuple_freeze_(Tuple, Prop),
 3579        init_propagator(B, Prop).
 3580
 3581
 3582tuples_in_conjunction(Tuples, Relation, Conj) :-
 3583        maplist(tuple_in_disjunction(Relation), Tuples, Disjs),
 3584        fold_statement(conjunction, Disjs, Conj).
 3585
 3586tuple_in_disjunction(Relation, Tuple, Disj) :-
 3587        maplist(tuple_in_conjunction(Tuple), Relation, Conjs),
 3588        fold_statement(disjunction, Conjs, Disj).
 3589
 3590tuple_in_conjunction(Tuple, Element, Conj) :-
 3591        maplist(var_eq, Tuple, Element, Eqs),
 3592        fold_statement(conjunction, Eqs, Conj).
 3593
 3594fold_statement(Operation, List, Statement) :-
 3595        (   List = [] -> Statement = 1
 3596        ;   List = [First|Rest],
 3597            foldl(Operation, Rest, First, Statement)
 3598        ).
 3599
 3600conjunction(E, Conj, Conj #/\ E).
 3601
 3602disjunction(E, Disj, Disj #\/ E).
 3603
 3604var_eq(V, N, ?(V) #= N).
 3605
 3606% Match variables to created skeleton.
 3607
 3608skeleton(Vs, Vs-Prop) :-
 3609        maplist(prop_init(Prop), Vs),
 3610        trigger_once(Prop).
 3611
 3612/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 3613   A drep is a user-accessible and visible domain representation. N,
 3614   N..M, and D1 \/ D2 are dreps, if D1 and D2 are dreps.
 3615- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 3616
 3617is_drep(N)      :- integer(N).
 3618is_drep(N..M)   :- drep_bound(N), drep_bound(M), N \== sup, M \== inf.
 3619is_drep(D1\/D2) :- is_drep(D1), is_drep(D2).
 3620is_drep({AI})   :- is_and_integers(AI).
 3621is_drep(\D)     :- is_drep(D).
 3622
 3623is_and_integers(I)     :- integer(I).
 3624is_and_integers((A,B)) :- is_and_integers(A), is_and_integers(B).
 3625
 3626drep_bound(I)   :- integer(I).
 3627drep_bound(sup).
 3628drep_bound(inf).
 3629
 3630drep_to_intervals(I)        --> { integer(I) }, [n(I)-n(I)].
 3631drep_to_intervals(N..M)     -->
 3632        (   { defaulty_to_bound(N, N1), defaulty_to_bound(M, M1),
 3633              N1 cis_leq M1} -> [N1-M1]
 3634        ;   []
 3635        ).
 3636drep_to_intervals(D1 \/ D2) -->
 3637        drep_to_intervals(D1), drep_to_intervals(D2).
 3638drep_to_intervals(\D0) -->
 3639        { drep_to_domain(D0, D1),
 3640          domain_complement(D1, D),
 3641          domain_to_drep(D, Drep) },
 3642        drep_to_intervals(Drep).
 3643drep_to_intervals({AI}) -->
 3644        and_integers_(AI).
 3645
 3646and_integers_(I)     --> { integer(I) }, [n(I)-n(I)].
 3647and_integers_((A,B)) --> and_integers_(A), and_integers_(B).
 3648
 3649drep_to_domain(DR, D) :-
 3650        must_be(ground, DR),
 3651        (   is_drep(DR) -> true
 3652        ;   domain_error(clpfd_domain, DR)
 3653        ),
 3654        phrase(drep_to_intervals(DR), Is0),
 3655        merge_intervals(Is0, Is1),
 3656        intervals_to_domain(Is1, D).
 3657
 3658merge_intervals(Is0, Is) :-
 3659        keysort(Is0, Is1),
 3660        merge_overlapping(Is1, Is).
 3661
 3662merge_overlapping([], []).
 3663merge_overlapping([A-B0|ABs0], [A-B|ABs]) :-
 3664        merge_remaining(ABs0, B0, B, Rest),
 3665        merge_overlapping(Rest, ABs).
 3666
 3667merge_remaining([], B, B, []).
 3668merge_remaining([N-M|NMs], B0, B, Rest) :-
 3669        Next cis B0 + n(1),
 3670        (   N cis_gt Next -> B = B0, Rest = [N-M|NMs]
 3671        ;   B1 cis max(B0,M),
 3672            merge_remaining(NMs, B1, B, Rest)
 3673        ).
 3674
 3675domain(V, Dom) :-
 3676        (   fd_get(V, Dom0, VPs) ->
 3677            domains_intersection(Dom, Dom0, Dom1),
 3678            %format("intersected\n: ~w\n ~w\n==> ~w\n\n", [Dom,Dom0,Dom1]),
 3679            fd_put(V, Dom1, VPs),
 3680            do_queue,
 3681            reinforce(V)
 3682        ;   domain_contains(Dom, V)
 3683        ).
 3684
 3685domains([], _).
 3686domains([V|Vs], D) :- domain(V, D), domains(Vs, D).
 3687
 3688props_number(fd_props(Gs,Bs,Os), N) :-
 3689        length(Gs, N1),
 3690        length(Bs, N2),
 3691        length(Os, N3),
 3692        N is N1 + N2 + N3.
 3693
 3694fd_get(X, Dom, Ps) :-
 3695        (   get_attr(X, clpfd, Attr) -> Attr = clpfd_attr(_,_,_,Dom,Ps)
 3696        ;   var(X) -> default_domain(Dom), Ps = fd_props([],[],[])
 3697        ).
 3698
 3699fd_get(X, Dom, Inf, Sup, Ps) :-
 3700        fd_get(X, Dom, Ps),
 3701        domain_infimum(Dom, Inf),
 3702        domain_supremum(Dom, Sup).
 3703
 3704/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 3705   By default, propagation always terminates. Currently, this is
 3706   ensured by allowing the left and right boundaries, as well as the
 3707   distance between the smallest and largest number occurring in the
 3708   domain representation to be changed at most once after a constraint
 3709   is posted, unless the domain is bounded. Set the experimental
 3710   Prolog flag 'clpfd_propagation' to 'full' to make the solver
 3711   propagate as much as possible. This can make queries
 3712   non-terminating, like: X #> abs(X), or: X #> Y, Y #> X, X #> 0.
 3713   Importantly, it can also make labeling non-terminating, as in:
 3714
 3715   ?- B #==> X #> abs(X), indomain(B).
 3716- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 3717
 3718fd_put(X, Dom, Ps) :-
 3719        (   current_prolog_flag(clpfd_propagation, full) ->
 3720            put_full(X, Dom, Ps)
 3721        ;   put_terminating(X, Dom, Ps)
 3722        ).
 3723
 3724put_terminating(X, Dom, Ps) :-
 3725        Dom \== empty,
 3726        (   Dom = from_to(F, F) -> F = n(X)
 3727        ;   (   get_attr(X, clpfd, Attr) ->
 3728                Attr = clpfd_attr(Left,Right,Spread,OldDom, _OldPs),
 3729                put_attr(X, clpfd, clpfd_attr(Left,Right,Spread,Dom,Ps)),
 3730                (   OldDom == Dom -> true
 3731                ;   (   Left == (.) -> Bounded = yes
 3732                    ;   domain_infimum(Dom, Inf), domain_supremum(Dom, Sup),
 3733                        (   Inf = n(_), Sup = n(_) ->
 3734                            Bounded = yes
 3735                        ;   Bounded = no
 3736                        )
 3737                    ),
 3738                    (   Bounded == yes ->
 3739                        put_attr(X, clpfd, clpfd_attr(.,.,.,Dom,Ps)),
 3740                        trigger_props(Ps, X, OldDom, Dom)
 3741                    ;   % infinite domain; consider border and spread changes
 3742                        domain_infimum(OldDom, OldInf),
 3743                        (   Inf == OldInf -> LeftP = Left
 3744                        ;   LeftP = yes
 3745                        ),
 3746                        domain_supremum(OldDom, OldSup),
 3747                        (   Sup == OldSup -> RightP = Right
 3748                        ;   RightP = yes
 3749                        ),
 3750                        domain_spread(OldDom, OldSpread),
 3751                        domain_spread(Dom, NewSpread),
 3752                        (   NewSpread == OldSpread -> SpreadP = Spread
 3753                        ;   NewSpread cis_lt OldSpread -> SpreadP = no
 3754                        ;   SpreadP = yes
 3755                        ),
 3756                        put_attr(X, clpfd, clpfd_attr(LeftP,RightP,SpreadP,Dom,Ps)),
 3757                        (   RightP == yes, Right = yes -> true
 3758                        ;   LeftP == yes, Left = yes -> true
 3759                        ;   SpreadP == yes, Spread = yes -> true
 3760                        ;   trigger_props(Ps, X, OldDom, Dom)
 3761                        )
 3762                    )
 3763                )
 3764            ;   var(X) ->
 3765                put_attr(X, clpfd, clpfd_attr(no,no,no,Dom, Ps))
 3766            ;   true
 3767            )
 3768        ).
 3769
 3770domain_spread(Dom, Spread) :-
 3771        domain_smallest_finite(Dom, S),
 3772        domain_largest_finite(Dom, L),
 3773        Spread cis L - S.
 3774
 3775smallest_finite(inf, Y, Y).
 3776smallest_finite(n(N), _, n(N)).
 3777
 3778domain_smallest_finite(from_to(F,T), S)   :- smallest_finite(F, T, S).
 3779domain_smallest_finite(split(_, L, _), S) :- domain_smallest_finite(L, S).
 3780
 3781largest_finite(sup, Y, Y).
 3782largest_finite(n(N), _, n(N)).
 3783
 3784domain_largest_finite(from_to(F,T), L)   :- largest_finite(T, F, L).
 3785domain_largest_finite(split(_, _, R), L) :- domain_largest_finite(R, L).
 3786
 3787/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 3788   With terminating propagation, all relevant constraints get a
 3789   propagation opportunity whenever a new constraint is posted.
 3790- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 3791
 3792reinforce(X) :-
 3793        (   current_prolog_flag(clpfd_propagation, full) ->
 3794            % full propagation propagates everything in any case
 3795            true
 3796        ;   term_variables(X, Vs),
 3797            maplist(reinforce_, Vs),
 3798            do_queue
 3799        ).
 3800
 3801reinforce_(X) :-
 3802        (   fd_var(X), fd_get(X, Dom, Ps) ->
 3803            put_full(X, Dom, Ps)
 3804        ;   true
 3805        ).
 3806
 3807put_full(X, Dom, Ps) :-
 3808        Dom \== empty,
 3809        (   Dom = from_to(F, F) -> F = n(X)
 3810        ;   (   get_attr(X, clpfd, Attr) ->
 3811                Attr = clpfd_attr(_,_,_,OldDom, _OldPs),
 3812                put_attr(X, clpfd, clpfd_attr(no,no,no,Dom, Ps)),
 3813                %format("putting dom: ~w\n", [Dom]),
 3814                (   OldDom == Dom -> true
 3815                ;   trigger_props(Ps, X, OldDom, Dom)
 3816                )
 3817            ;   var(X) -> %format('\t~w in ~w .. ~w\n',[X,L,U]),
 3818                put_attr(X, clpfd, clpfd_attr(no,no,no,Dom, Ps))
 3819            ;   true
 3820            )
 3821        ).
 3822
 3823/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 3824   A propagator is a term of the form propagator(C, State), where C
 3825   represents a constraint, and State is a free variable that can be
 3826   used to destructively change the state of the propagator via
 3827   attributes. This can be used to avoid redundant invocation of the
 3828   same propagator, or to disable the propagator.
 3829- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 3830
 3831make_propagator(C, propagator(C, _)).
 3832
 3833propagator_state(propagator(_,S), S).
 3834
 3835trigger_props(fd_props(Gs,Bs,Os), X, D0, D) :-
 3836        (   ground(X) ->
 3837            trigger_props_(Gs),
 3838            trigger_props_(Bs)
 3839        ;   Bs \== [] ->
 3840            domain_infimum(D0, I0),
 3841            domain_infimum(D, I),
 3842            (   I == I0 ->
 3843                domain_supremum(D0, S0),
 3844                domain_supremum(D, S),
 3845                (   S == S0 -> true
 3846                ;   trigger_props_(Bs)
 3847                )
 3848            ;   trigger_props_(Bs)
 3849            )
 3850        ;   true
 3851        ),
 3852        trigger_props_(Os).
 3853
 3854trigger_props(fd_props(Gs,Bs,Os), X) :-
 3855        trigger_props_(Os),
 3856        trigger_props_(Bs),
 3857        (   ground(X) ->
 3858            trigger_props_(Gs)
 3859        ;   true
 3860        ).
 3861
 3862trigger_props(fd_props(Gs,Bs,Os)) :-
 3863        trigger_props_(Gs),
 3864        trigger_props_(Bs),
 3865        trigger_props_(Os).
 3866
 3867trigger_props_([]).
 3868trigger_props_([P|Ps]) :- trigger_prop(P), trigger_props_(Ps).
 3869
 3870trigger_prop(Propagator) :-
 3871        propagator_state(Propagator, State),
 3872        (   State == dead -> true
 3873        ;   get_attr(State, clpfd_aux, queued) -> true
 3874        ;   b_getval('$clpfd_current_propagator', C), C == State -> true
 3875        ;   % passive
 3876            % format("triggering: ~w\n", [Propagator]),
 3877            put_attr(State, clpfd_aux, queued),
 3878            (   arg(1, Propagator, C), functor(C, F, _), global_constraint(F) ->
 3879                push_queue(Propagator, 2)
 3880            ;   push_queue(Propagator, 1)
 3881            )
 3882        ).
 3883
 3884kill(State) :- del_attr(State, clpfd_aux), State = dead.
 3885
 3886kill(State, Ps) :-
 3887        kill(State),
 3888        maplist(kill_entailed, Ps).
 3889
 3890kill_entailed(p(Prop)) :-
 3891        propagator_state(Prop, State),
 3892        kill(State).
 3893kill_entailed(a(V)) :-
 3894        del_attr(V, clpfd).
 3895kill_entailed(a(X,B)) :-
 3896        (   X == B -> true
 3897        ;   del_attr(B, clpfd)
 3898        ).
 3899kill_entailed(a(X,Y,B)) :-
 3900        (   X == B -> true
 3901        ;   Y == B -> true
 3902        ;   del_attr(B, clpfd)
 3903        ).
 3904
 3905no_reactivation(rel_tuple(_,_)).
 3906no_reactivation(pdistinct(_)).
 3907no_reactivation(pgcc(_,_,_)).
 3908no_reactivation(pgcc_single(_,_)).
 3909%no_reactivation(scalar_product(_,_,_,_)).
 3910
 3911activate_propagator(propagator(P,State)) :-
 3912        % format("running: ~w\n", [P]),
 3913        del_attr(State, clpfd_aux),
 3914        (   no_reactivation(P) ->
 3915            b_setval('$clpfd_current_propagator', State),
 3916            run_propagator(P, State),
 3917            b_setval('$clpfd_current_propagator', [])
 3918        ;   run_propagator(P, State)
 3919        ).
 3920
 3921disable_queue :- b_setval('$clpfd_queue_status', disabled).
 3922enable_queue  :- b_setval('$clpfd_queue_status', enabled).
 3923
 3924portray_propagator(propagator(P,_), F) :- functor(P, F, _).
 3925
 3926portray_queue(V, []) :- var(V), !.
 3927portray_queue([P|Ps], [F|Fs]) :-
 3928        portray_propagator(P, F),
 3929        portray_queue(Ps, Fs).
 3930
 3931do_queue :-
 3932        % b_getval('$clpfd_queue', H-_),
 3933        % portray_queue(H, Port),
 3934        % format("queue: ~w\n", [Port]),
 3935        (   b_getval('$clpfd_queue_status', enabled) ->
 3936            (   fetch_propagator(Propagator) ->
 3937                activate_propagator(Propagator),
 3938                do_queue
 3939            ;   true
 3940            )
 3941        ;   true
 3942        ).
 3943
 3944init_propagator(Var, Prop) :-
 3945        (   fd_get(Var, Dom, Ps0) ->
 3946            insert_propagator(Prop, Ps0, Ps),
 3947            fd_put(Var, Dom, Ps)
 3948        ;   true
 3949        ).
 3950
 3951constraint_wake(pneq, ground).
 3952constraint_wake(x_neq_y_plus_z, ground).
 3953constraint_wake(absdiff_neq, ground).
 3954constraint_wake(pdifferent, ground).
 3955constraint_wake(pexclude, ground).
 3956constraint_wake(scalar_product_neq, ground).
 3957
 3958constraint_wake(x_leq_y_plus_c, bounds).
 3959constraint_wake(scalar_product_eq, bounds).
 3960constraint_wake(scalar_product_leq, bounds).
 3961constraint_wake(pplus, bounds).
 3962constraint_wake(pgeq, bounds).
 3963constraint_wake(pgcc_single, bounds).
 3964constraint_wake(pgcc_check_single, bounds).
 3965
 3966global_constraint(pdistinct).
 3967global_constraint(pgcc).
 3968global_constraint(pgcc_single).
 3969global_constraint(pcircuit).
 3970%global_constraint(rel_tuple).
 3971%global_constraint(scalar_product_eq).
 3972
 3973insert_propagator(Prop, Ps0, Ps) :-
 3974        Ps0 = fd_props(Gs,Bs,Os),
 3975        arg(1, Prop, Constraint),
 3976        functor(Constraint, F, _),
 3977        (   constraint_wake(F, ground) ->
 3978            Ps = fd_props([Prop|Gs], Bs, Os)
 3979        ;   constraint_wake(F, bounds) ->
 3980            Ps = fd_props(Gs, [Prop|Bs], Os)
 3981        ;   Ps = fd_props(Gs, Bs, [Prop|Os])
 3982        ).
 3983
 3984%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 3985
 3986%% lex_chain(+Lists)
 3987%
 3988% Lists are lexicographically non-decreasing.
 3989
 3990lex_chain(Lss) :-
 3991        must_be(list(list), Lss),
 3992        maplist(maplist(fd_variable), Lss),
 3993        (   Lss == [] -> true
 3994        ;   Lss = [First|Rest],
 3995            make_propagator(presidual(lex_chain(Lss)), Prop),
 3996            foldl(lex_chain_(Prop), Rest, First, _)
 3997        ).
 3998
 3999lex_chain_(Prop, Ls, Prev, Ls) :-
 4000        maplist(prop_init(Prop), Ls),
 4001        lex_le(Prev, Ls).
 4002
 4003lex_le([], []).
 4004lex_le([V1|V1s], [V2|V2s]) :-
 4005        ?(V1) #=< ?(V2),
 4006        (   integer(V1) ->
 4007            (   integer(V2) ->
 4008                (   V1 =:= V2 -> lex_le(V1s, V2s) ;  true )
 4009            ;   freeze(V2, lex_le([V1|V1s], [V2|V2s]))
 4010            )
 4011        ;   freeze(V1, lex_le([V1|V1s], [V2|V2s]))
 4012        ).
 4013
 4014%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 4015
 4016
 4017%% tuples_in(+Tuples, +Relation).
 4018%
 4019% True iff all Tuples are elements of Relation. Each element of the
 4020% list Tuples is a list of integers or finite domain variables.
 4021% Relation is a list of lists of integers. Arbitrary finite relations,
 4022% such as compatibility tables, can be modeled in this way. For
 4023% example, if 1 is compatible with 2 and 5, and 4 is compatible with 0
 4024% and 3:
 4025%
 4026% ==
 4027% ?- tuples_in([[X,Y]], [[1,2],[1,5],[4,0],[4,3]]), X = 4.
 4028% X = 4,
 4029% Y in 0\/3.
 4030% ==
 4031%
 4032% As another example, consider a train schedule represented as a list
 4033% of quadruples, denoting departure and arrival places and times for
 4034% each train. In the following program, Ps is a feasible journey of
 4035% length 3 from A to D via trains that are part of the given schedule.
 4036%
 4037% ==
 4038% trains([[1,2,0,1],
 4039%         [2,3,4,5],
 4040%         [2,3,0,1],
 4041%         [3,4,5,6],
 4042%         [3,4,2,3],
 4043%         [3,4,8,9]]).
 4044%
 4045% threepath(A, D, Ps) :-
 4046%         Ps = [[A,B,_T0,T1],[B,C,T2,T3],[C,D,T4,_T5]],
 4047%         T2 #> T1,
 4048%         T4 #> T3,
 4049%         trains(Ts),
 4050%         tuples_in(Ps, Ts).
 4051% ==
 4052%
 4053% In this example, the unique solution is found without labeling:
 4054%
 4055% ==
 4056% ?- threepath(1, 4, Ps).
 4057% Ps = [[1, 2, 0, 1], [2, 3, 4, 5], [3, 4, 8, 9]].
 4058% ==
 4059
 4060tuples_in(Tuples, Relation) :-
 4061        must_be(list(list), Tuples),
 4062        maplist(maplist(fd_variable), Tuples),
 4063        must_be(list(list(integer)), Relation),
 4064        maplist(relation_tuple(Relation), Tuples),
 4065        do_queue.
 4066
 4067relation_tuple(Relation, Tuple) :-
 4068        relation_unifiable(Relation, Tuple, Us, _, _),
 4069        (   ground(Tuple) -> memberchk(Tuple, Relation)
 4070        ;   tuple_domain(Tuple, Us),
 4071            (   Tuple = [_,_|_] -> tuple_freeze(Tuple, Us)
 4072            ;   true
 4073            )
 4074        ).
 4075
 4076tuple_domain([], _).
 4077tuple_domain([T|Ts], Relation0) :-
 4078        maplist(list_first_rest, Relation0, Firsts, Relation1),
 4079        (   var(T) ->
 4080            (   Firsts = [Unique] -> T = Unique
 4081            ;   list_to_domain(Firsts, FDom),
 4082                fd_get(T, TDom, TPs),
 4083                domains_intersection(TDom, FDom, TDom1),
 4084                fd_put(T, TDom1, TPs)
 4085            )
 4086        ;   true
 4087        ),
 4088        tuple_domain(Ts, Relation1).
 4089
 4090tuple_freeze(Tuple, Relation) :-
 4091        put_attr(R, clpfd_relation, Relation),
 4092        make_propagator(rel_tuple(R, Tuple), Prop),
 4093        tuple_freeze_(Tuple, Prop).
 4094
 4095tuple_freeze_([], _).
 4096tuple_freeze_([T|Ts], Prop) :-
 4097        (   var(T) ->
 4098            init_propagator(T, Prop),
 4099            trigger_prop(Prop)
 4100        ;   true
 4101        ),
 4102        tuple_freeze_(Ts, Prop).
 4103
 4104relation_unifiable([], _, [], Changed, Changed).
 4105relation_unifiable([R|Rs], Tuple, Us, Changed0, Changed) :-
 4106        (   all_in_domain(R, Tuple) ->
 4107            Us = [R|Rest],
 4108            relation_unifiable(Rs, Tuple, Rest, Changed0, Changed)
 4109        ;   relation_unifiable(Rs, Tuple, Us, true, Changed)
 4110        ).
 4111
 4112all_in_domain([], []).
 4113all_in_domain([A|As], [T|Ts]) :-
 4114        (   fd_get(T, Dom, _) ->
 4115            domain_contains(Dom, A)
 4116        ;   T =:= A
 4117        ),
 4118        all_in_domain(As, Ts).
 4119
 4120%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 4121
 4122% trivial propagator, used only to remember pending constraints
 4123run_propagator(presidual(_), _).
 4124
 4125%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 4126run_propagator(pdifferent(Left,Right,X,_), MState) :-
 4127        run_propagator(pexclude(Left,Right,X), MState).
 4128
 4129run_propagator(weak_distinct(Left,Right,X,_), _MState) :-
 4130        (   ground(X) ->
 4131            disable_queue,
 4132            exclude_fire(Left, Right, X),
 4133            enable_queue
 4134        ;   outof_reducer(Left, Right, X)
 4135            %(   var(X) -> kill_if_isolated(Left, Right, X, MState)
 4136            %;   true
 4137            %)
 4138        ).
 4139
 4140run_propagator(pexclude(Left,Right,X), _) :-
 4141        (   ground(X) ->
 4142            disable_queue,
 4143            exclude_fire(Left, Right, X),
 4144            enable_queue
 4145        ;   true
 4146        ).
 4147
 4148run_propagator(pdistinct(Ls), _MState) :-
 4149        distinct(Ls).
 4150
 4151run_propagator(check_distinct(Left,Right,X), _) :-
 4152        \+ list_contains(Left, X),
 4153        \+ list_contains(Right, X).
 4154
 4155%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 4156
 4157run_propagator(pelement(N, Is, V), MState) :-
 4158        (   fd_get(N, NDom, _) ->
 4159            (   fd_get(V, VDom, VPs) ->
 4160                integers_remaining(Is, 1, NDom, empty, VDom1),
 4161                domains_intersection(VDom, VDom1, VDom2),
 4162                fd_put(V, VDom2, VPs)
 4163            ;   true
 4164            )
 4165        ;   kill(MState), nth1(N, Is, V)
 4166        ).
 4167
 4168%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 4169
 4170run_propagator(pgcc_single(Vs, Pairs), _) :- gcc_global(Vs, Pairs).
 4171
 4172run_propagator(pgcc_check_single(Pairs), _) :- gcc_check(Pairs).
 4173
 4174run_propagator(pgcc_check(Pairs), _) :- gcc_check(Pairs).
 4175
 4176run_propagator(pgcc(Vs, _, Pairs), _) :- gcc_global(Vs, Pairs).
 4177
 4178%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 4179
 4180run_propagator(pcircuit(Vs), _MState) :-
 4181        distinct(Vs),
 4182        propagate_circuit(Vs).
 4183
 4184%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 4185run_propagator(pneq(A, B), MState) :-
 4186        (   nonvar(A) ->
 4187            (   nonvar(B) -> A =\= B, kill(MState)
 4188            ;   fd_get(B, BD0, BExp0),
 4189                domain_remove(BD0, A, BD1),
 4190                kill(MState),
 4191                fd_put(B, BD1, BExp0)
 4192            )
 4193        ;   nonvar(B) -> run_propagator(pneq(B, A), MState)
 4194        ;   A \== B,
 4195            fd_get(A, _, AI, AS, _), fd_get(B, _, BI, BS, _),
 4196            (   AS cis_lt BI -> kill(MState)
 4197            ;   AI cis_gt BS -> kill(MState)
 4198            ;   true
 4199            )
 4200        ).
 4201
 4202%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 4203run_propagator(pgeq(A,B), MState) :-
 4204        (   A == B -> kill(MState)
 4205        ;   nonvar(A) ->
 4206            (   nonvar(B) -> kill(MState), A >= B
 4207            ;   fd_get(B, BD, BPs),
 4208                domain_remove_greater_than(BD, A, BD1),
 4209                kill(MState),
 4210                fd_put(B, BD1, BPs)
 4211            )
 4212        ;   nonvar(B) ->
 4213            fd_get(A, AD, APs),
 4214            domain_remove_smaller_than(AD, B, AD1),
 4215            kill(MState),
 4216            fd_put(A, AD1, APs)
 4217        ;   fd_get(A, AD, AL, AU, APs),
 4218            fd_get(B, _, BL, BU, _),
 4219            AU cis_geq BL,
 4220            (   AL cis_geq BU -> kill(MState)
 4221            ;   AU == BL -> kill(MState), A = B
 4222            ;   NAL cis max(AL,BL),
 4223                domains_intersection(AD, from_to(NAL,AU), NAD),
 4224                fd_put(A, NAD, APs),
 4225                (   fd_get(B, BD2, BL2, BU2, BPs2) ->
 4226                    NBU cis min(BU2, AU),
 4227                    domains_intersection(BD2, from_to(BL2,NBU), NBD),
 4228                    fd_put(B, NBD, BPs2)
 4229                ;   true
 4230                )
 4231            )
 4232        ).
 4233
 4234%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 4235
 4236run_propagator(rel_tuple(R, Tuple), MState) :-
 4237        get_attr(R, clpfd_relation, Relation),
 4238        (   ground(Tuple) -> kill(MState), memberchk(Tuple, Relation)
 4239        ;   relation_unifiable(Relation, Tuple, Us, false, Changed),
 4240            Us = [_|_],
 4241            (   Tuple = [First,Second], ( ground(First) ; ground(Second) ) ->
 4242                kill(MState)
 4243            ;   true
 4244            ),
 4245            (   Us = [Single] -> kill(MState), Single = Tuple
 4246            ;   Changed ->
 4247                put_attr(R, clpfd_relation, Us),
 4248                disable_queue,
 4249                tuple_domain(Tuple, Us),
 4250                enable_queue
 4251            ;   true
 4252            )
 4253        ).
 4254
 4255%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 4256
 4257run_propagator(pserialized(S_I, D_I, S_J, D_J, _), MState) :-
 4258        (   nonvar(S_I), nonvar(S_J) ->
 4259            kill(MState),
 4260            (   S_I + D_I =< S_J -> true
 4261            ;   S_J + D_J =< S_I -> true
 4262            ;   false
 4263            )
 4264        ;   serialize_lower_upper(S_I, D_I, S_J, D_J, MState),
 4265            serialize_lower_upper(S_J, D_J, S_I, D_I, MState)
 4266        ).
 4267
 4268%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 4269
 4270% abs(X-Y) #\= C
 4271run_propagator(absdiff_neq(X,Y,C), MState) :-
 4272        (   C < 0 -> kill(MState)
 4273        ;   nonvar(X) ->
 4274            kill(MState),
 4275            (   nonvar(Y) -> abs(X - Y) =\= C
 4276            ;   V1 is X - C, neq_num(Y, V1),
 4277                V2 is C + X, neq_num(Y, V2)
 4278            )
 4279        ;   nonvar(Y) -> kill(MState),
 4280            V1 is C + Y, neq_num(X, V1),
 4281            V2 is Y - C, neq_num(X, V2)
 4282        ;   true
 4283        ).
 4284
 4285% abs(X-Y) #>= C
 4286run_propagator(absdiff_geq(X,Y,C), MState) :-
 4287        (   C =< 0 -> kill(MState)
 4288        ;   nonvar(X) ->
 4289            kill(MState),
 4290            (   nonvar(Y) -> abs(X-Y) >= C
 4291            ;   P1 is X - C, P2 is X + C,
 4292                Y in inf..P1 \/ P2..sup
 4293            )
 4294        ;   nonvar(Y) ->
 4295            kill(MState),
 4296            P1 is Y - C, P2 is Y + C,
 4297            X in inf..P1 \/ P2..sup
 4298        ;   true
 4299        ).
 4300
 4301% X #\= Y + Z
 4302run_propagator(x_neq_y_plus_z(X,Y,Z), MState) :-
 4303        (   nonvar(X) ->
 4304            (   nonvar(Y) ->
 4305                (   nonvar(Z) -> kill(MState), X =\= Y + Z
 4306                ;   kill(MState), XY is X - Y, neq_num(Z, XY)
 4307                )
 4308            ;   nonvar(Z) -> kill(MState), XZ is X - Z, neq_num(Y, XZ)
 4309            ;   true
 4310            )
 4311        ;   nonvar(Y) ->
 4312            (   nonvar(Z) ->
 4313                kill(MState), YZ is Y + Z, neq_num(X, YZ)
 4314            ;   Y =:= 0 -> kill(MState), neq(X, Z)
 4315            ;   true
 4316            )
 4317        ;   Z == 0 -> kill(MState), neq(X, Y)
 4318        ;   true
 4319        ).
 4320
 4321% X #=< Y + C
 4322run_propagator(x_leq_y_plus_c(X,Y,C), MState) :-
 4323        (   nonvar(X) ->
 4324            (   nonvar(Y) -> kill(MState), X =< Y + C
 4325            ;   kill(MState),
 4326                R is X - C,
 4327                fd_get(Y, YD, YPs),
 4328                domain_remove_smaller_than(YD, R, YD1),
 4329                fd_put(Y, YD1, YPs)
 4330            )
 4331        ;   nonvar(Y) ->
 4332            kill(MState),
 4333            R is Y + C,
 4334            fd_get(X, XD, XPs),
 4335            domain_remove_greater_than(XD, R, XD1),
 4336            fd_put(X, XD1, XPs)
 4337        ;   (   X == Y -> C >= 0, kill(MState)
 4338            ;   fd_get(Y, YD, _),
 4339                (   domain_supremum(YD, n(YSup)) ->
 4340                    YS1 is YSup + C,
 4341                    fd_get(X, XD, XPs),
 4342                    domain_remove_greater_than(XD, YS1, XD1),
 4343                    fd_put(X, XD1, XPs)
 4344                ;   true
 4345                ),
 4346                (   fd_get(X, XD2, _), domain_infimum(XD2, n(XInf)) ->
 4347                    XI1 is XInf - C,
 4348                    (   fd_get(Y, YD1, YPs1) ->
 4349                        domain_remove_smaller_than(YD1, XI1, YD2),
 4350                        (   domain_infimum(YD2, n(YInf)),
 4351                            domain_supremum(XD2, n(XSup)),
 4352                            XSup =< YInf + C ->
 4353                            kill(MState)
 4354                        ;   true
 4355                        ),
 4356                        fd_put(Y, YD2, YPs1)
 4357                    ;   true
 4358                    )
 4359                ;   true
 4360                )
 4361            )
 4362        ).
 4363
 4364run_propagator(scalar_product_neq(Cs0,Vs0,P0), MState) :-
 4365        coeffs_variables_const(Cs0, Vs0, Cs, Vs, 0, I),
 4366        P is P0 - I,
 4367        (   Vs = [] -> kill(MState), P =\= 0
 4368        ;   Vs = [V], Cs = [C] ->
 4369            kill(MState),
 4370            (   C =:= 1 -> neq_num(V, P)
 4371            ;   C*V #\= P
 4372            )
 4373        ;   Cs == [1,-1] -> kill(MState), Vs = [A,B], x_neq_y_plus_z(A, B, P)
 4374        ;   Cs == [-1,1] -> kill(MState), Vs = [A,B], x_neq_y_plus_z(B, A, P)
 4375        ;   P =:= 0, Cs = [1,1,-1] ->
 4376            kill(MState), Vs = [A,B,C], x_neq_y_plus_z(C, A, B)
 4377        ;   P =:= 0, Cs = [1,-1,1] ->
 4378            kill(MState), Vs = [A,B,C], x_neq_y_plus_z(B, A, C)
 4379        ;   P =:= 0, Cs = [-1,1,1] ->
 4380            kill(MState), Vs = [A,B,C], x_neq_y_plus_z(A, B, C)
 4381        ;   true
 4382        ).
 4383
 4384run_propagator(scalar_product_leq(Cs0,Vs0,P0), MState) :-
 4385        coeffs_variables_const(Cs0, Vs0, Cs, Vs, 0, I),
 4386        P is P0 - I,
 4387        (   Vs = [] -> kill(MState), P >= 0
 4388        ;   sum_finite_domains(Cs, Vs, Infs, Sups, 0, 0, Inf, Sup),
 4389            D1 is P - Inf,
 4390            disable_queue,
 4391            (   Infs == [], Sups == [] ->
 4392                Inf =< P,
 4393                (   Sup =< P -> kill(MState)
 4394                ;   remove_dist_upper_leq(Cs, Vs, D1)
 4395                )
 4396            ;   Infs == [] -> Inf =< P, remove_dist_upper(Sups, D1)
 4397            ;   Sups = [_], Infs = [_] ->
 4398                remove_upper(Infs, D1)
 4399            ;   Infs = [_] -> remove_upper(Infs, D1)
 4400            ;   true
 4401            ),
 4402            enable_queue
 4403        ).
 4404
 4405run_propagator(scalar_product_eq(Cs0,Vs0,P0), MState) :-
 4406        coeffs_variables_const(Cs0, Vs0, Cs, Vs, 0, I),
 4407        P is P0 - I,
 4408        (   Vs = [] -> kill(MState), P =:= 0
 4409        ;   Vs = [V], Cs = [C] -> kill(MState), P mod C =:= 0, V is P // C
 4410        ;   Cs == [1,1] -> kill(MState), Vs = [A,B], A + B #= P
 4411        ;   Cs == [1,-1] -> kill(MState), Vs = [A,B], A #= P + B
 4412        ;   Cs == [-1,1] -> kill(MState), Vs = [A,B], B #= P + A
 4413        ;   Cs == [-1,-1] -> kill(MState), Vs = [A,B], P1 is -P, A + B #= P1
 4414        ;   P =:= 0, Cs == [1,1,-1] -> kill(MState), Vs = [A,B,C], A + B #= C
 4415        ;   P =:= 0, Cs == [1,-1,1] -> kill(MState), Vs = [A,B,C], A + C #= B
 4416        ;   P =:= 0, Cs == [-1,1,1] -> kill(MState), Vs = [A,B,C], B + C #= A
 4417        ;   sum_finite_domains(Cs, Vs, Infs, Sups, 0, 0, Inf, Sup),
 4418            % nl, writeln(Infs-Sups-Inf-Sup),
 4419            D1 is P - Inf,
 4420            D2 is Sup - P,
 4421            disable_queue,
 4422            (   Infs == [], Sups == [] ->
 4423                between(Inf, Sup, P),
 4424                remove_dist_upper_lower(Cs, Vs, D1, D2)
 4425            ;   Sups = [] -> P =< Sup, remove_dist_lower(Infs, D2)
 4426            ;   Infs = [] -> Inf =< P, remove_dist_upper(Sups, D1)
 4427            ;   Sups = [_], Infs = [_] ->
 4428                remove_lower(Sups, D2),
 4429                remove_upper(Infs, D1)
 4430            ;   Infs = [_] -> remove_upper(Infs, D1)
 4431            ;   Sups = [_] -> remove_lower(Sups, D2)
 4432            ;   true
 4433            ),
 4434            enable_queue
 4435        ).
 4436
 4437% X + Y = Z
 4438run_propagator(pplus(X,Y,Z), MState) :-
 4439        (   nonvar(X) ->
 4440            (   X =:= 0 -> kill(MState), Y = Z
 4441            ;   Y == Z -> kill(MState), X =:= 0
 4442            ;   nonvar(Y) -> kill(MState), Z is X + Y
 4443            ;   nonvar(Z) -> kill(MState), Y is Z - X
 4444            ;   fd_get(Z, ZD, ZPs),
 4445                fd_get(Y, YD, _),
 4446                domain_shift(YD, X, Shifted_YD),
 4447                domains_intersection(ZD, Shifted_YD, ZD1),
 4448                fd_put(Z, ZD1, ZPs),
 4449                (   fd_get(Y, YD1, YPs) ->
 4450                    O is -X,
 4451                    domain_shift(ZD1, O, YD2),
 4452                    domains_intersection(YD1, YD2, YD3),
 4453                    fd_put(Y, YD3, YPs)
 4454                ;   true
 4455                )
 4456            )
 4457        ;   nonvar(Y) -> run_propagator(pplus(Y,X,Z), MState)
 4458        ;   nonvar(Z) ->
 4459            (   X == Y -> kill(MState), even(Z), X is Z // 2
 4460            ;   fd_get(X, XD, _),
 4461                fd_get(Y, YD, YPs),
 4462                domain_negate(XD, XDN),
 4463                domain_shift(XDN, Z, YD1),
 4464                domains_intersection(YD, YD1, YD2),
 4465                fd_put(Y, YD2, YPs),
 4466                (   fd_get(X, XD1, XPs) ->
 4467                    domain_negate(YD2, YD2N),
 4468                    domain_shift(YD2N, Z, XD2),
 4469                    domains_intersection(XD1, XD2, XD3),
 4470                    fd_put(X, XD3, XPs)
 4471                ;   true
 4472                )
 4473            )
 4474        ;   (   X == Y -> kill(MState), 2*X #= Z
 4475            ;   X == Z -> kill(MState), Y = 0
 4476            ;   Y == Z -> kill(MState), X = 0
 4477            ;   fd_get(X, XD, XL, XU, XPs),
 4478                fd_get(Y, _, YL, YU, _),
 4479                fd_get(Z, _, ZL, ZU, _),
 4480                NXL cis max(XL, ZL-YU),
 4481                NXU cis min(XU, ZU-YL),
 4482                update_bounds(X, XD, XPs, XL, XU, NXL, NXU),
 4483                (   fd_get(Y, YD2, YL2, YU2, YPs2) ->
 4484                    NYL cis max(YL2, ZL-NXU),
 4485                    NYU cis min(YU2, ZU-NXL),
 4486                    update_bounds(Y, YD2, YPs2, YL2, YU2, NYL, NYU)
 4487                ;   NYL = n(Y), NYU = n(Y)
 4488                ),
 4489                (   fd_get(Z, ZD2, ZL2, ZU2, ZPs2) ->
 4490                    NZL cis max(ZL2,NXL+NYL),
 4491                    NZU cis min(ZU2,NXU+NYU),
 4492                    update_bounds(Z, ZD2, ZPs2, ZL2, ZU2, NZL, NZU)
 4493                ;   true
 4494                )
 4495            )
 4496        ).
 4497
 4498%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 4499
 4500run_propagator(ptimes(X,Y,Z), MState) :-
 4501        (   nonvar(X) ->
 4502            (   nonvar(Y) -> kill(MState), Z is X * Y
 4503            ;   X =:= 0 -> kill(MState), Z = 0
 4504            ;   X =:= 1 -> kill(MState), Z = Y
 4505            ;   nonvar(Z) -> kill(MState), 0 =:= Z mod X, Y is Z // X
 4506            ;   (   Y == Z -> kill(MState), Y = 0
 4507                ;   fd_get(Y, YD, _),
 4508                    fd_get(Z, ZD, ZPs),
 4509                    domain_expand(YD, X, Scaled_YD),
 4510                    domains_intersection(ZD, Scaled_YD, ZD1),
 4511                    fd_put(Z, ZD1, ZPs),
 4512                    (   fd_get(Y, YDom2, YPs2) ->
 4513                        domain_contract(ZD1, X, Contract),
 4514                        domains_intersection(YDom2, Contract, NYDom),
 4515                        fd_put(Y, NYDom, YPs2)
 4516                    ;   kill(MState), Z is X * Y
 4517                    )
 4518                )
 4519            )
 4520        ;   nonvar(Y) -> run_propagator(ptimes(Y,X,Z), MState)
 4521        ;   nonvar(Z) ->
 4522            (   X == Y ->
 4523                kill(MState),
 4524                integer_kth_root(Z, 2, R),
 4525                NR is -R,
 4526                X in NR \/ R
 4527            ;   fd_get(X, XD, XL, XU, XPs),
 4528                fd_get(Y, YD, YL, YU, _),
 4529                min_max_factor(n(Z), n(Z), YL, YU, XL, XU, NXL, NXU),
 4530                update_bounds(X, XD, XPs, XL, XU, NXL, NXU),
 4531                (   fd_get(Y, YD2, YL2, YU2, YPs2) ->
 4532                    min_max_factor(n(Z), n(Z), NXL, NXU, YL2, YU2, NYL, NYU),
 4533                    update_bounds(Y, YD2, YPs2, YL2, YU2, NYL, NYU)
 4534                ;   (   Y =\= 0 -> 0 =:= Z mod Y, kill(MState), X is Z // Y
 4535                    ;   kill(MState), Z = 0
 4536                    )
 4537                ),
 4538                (   Z =:= 0 ->
 4539                    (   \+ domain_contains(XD, 0) -> kill(MState), Y = 0
 4540                    ;   \+ domain_contains(YD, 0) -> kill(MState), X = 0
 4541                    ;   true
 4542                    )
 4543                ;  neq_num(X, 0), neq_num(Y, 0)
 4544                )
 4545            )
 4546        ;   (   X == Y -> kill(MState), X^2 #= Z
 4547            ;   fd_get(X, XD, XL, XU, XPs),
 4548                fd_get(Y, _, YL, YU, _),
 4549                fd_get(Z, ZD, ZL, ZU, _),
 4550                (   Y == Z, \+ domain_contains(ZD, 0) -> kill(MState), X = 1
 4551                ;   X == Z, \+ domain_contains(ZD, 0) -> kill(MState), Y = 1
 4552                ;   min_max_factor(ZL, ZU, YL, YU, XL, XU, NXL, NXU),
 4553                    update_bounds(X, XD, XPs, XL, XU, NXL, NXU),
 4554                    (   fd_get(Y, YD2, YL2, YU2, YPs2) ->
 4555                        min_max_factor(ZL, ZU, NXL, NXU, YL2, YU2, NYL, NYU),
 4556                        update_bounds(Y, YD2, YPs2, YL2, YU2, NYL, NYU)
 4557                    ;   NYL = n(Y), NYU = n(Y)
 4558                    ),
 4559                    (   fd_get(Z, ZD2, ZL2, ZU2, ZPs2) ->
 4560                        min_product(NXL, NXU, NYL, NYU, NZL),
 4561                        max_product(NXL, NXU, NYL, NYU, NZU),
 4562                        (   NZL cis_leq ZL2, NZU cis_geq ZU2 -> ZD3 = ZD2
 4563                        ;   domains_intersection(ZD2, from_to(NZL,NZU), ZD3),
 4564                            fd_put(Z, ZD3, ZPs2)
 4565                        ),
 4566                        (   domain_contains(ZD3, 0) ->  true
 4567                        ;   neq_num(X, 0), neq_num(Y, 0)
 4568                        )
 4569                    ;   true
 4570                    )
 4571                )
 4572            )
 4573        ).
 4574
 4575%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 4576
 4577% X div Y = Z
 4578run_propagator(pdiv(X,Y,Z), MState) :- kill(MState), Z #= (X-(X mod Y)) // Y.
 4579
 4580% X rdiv Y = Z
 4581run_propagator(prdiv(X,Y,Z), MState) :- kill(MState), Z*Y #= X.
 4582
 4583% X // Y = Z (round towards zero)
 4584run_propagator(ptzdiv(X,Y,Z), MState) :-
 4585        (   nonvar(X) ->
 4586            (   nonvar(Y) -> kill(MState), Y =\= 0, Z is X // Y
 4587            ;   fd_get(Y, YD, YL, YU, YPs),
 4588                (   nonvar(Z) ->
 4589                    (   Z =:= 0 ->
 4590                        NYL is -abs(X) - 1,
 4591                        NYU is abs(X) + 1,
 4592                        domains_intersection(YD, split(0, from_to(inf,n(NYL)),
 4593                                                       from_to(n(NYU), sup)),
 4594                                             NYD),
 4595                        fd_put(Y, NYD, YPs)
 4596                    ;   (   sign(X) =:= sign(Z) ->
 4597                            NYL cis max(n(X) // (n(Z)+sign(n(Z))) + n(1), YL),
 4598                            NYU cis min(n(X) // n(Z), YU)
 4599                        ;   NYL cis max(n(X) // n(Z), YL),
 4600                            NYU cis min(n(X) // (n(Z)+sign(n(Z))) - n(1), YU)
 4601                        ),
 4602                        update_bounds(Y, YD, YPs, YL, YU, NYL, NYU)
 4603                    )
 4604                ;   fd_get(Z, ZD, ZL, ZU, ZPs),
 4605                    (   X >= 0, ( YL cis_gt n(0) ; YU cis_lt n(0) )->
 4606                        NZL cis max(n(X)//YU, ZL),
 4607                        NZU cis min(n(X)//YL, ZU)
 4608                    ;   X < 0, ( YL cis_gt n(0) ; YU cis_lt n(0) ) ->
 4609                        NZL cis max(n(X)//YL, ZL),
 4610                        NZU cis min(n(X)//YU, ZU)
 4611                    ;   % TODO: more stringent bounds, cover Y
 4612                        NZL cis max(-abs(n(X)), ZL),
 4613                        NZU cis min(abs(n(X)), ZU)
 4614                    ),
 4615                    update_bounds(Z, ZD, ZPs, ZL, ZU, NZL, NZU),
 4616                    (   X >= 0, NZL cis_gt n(0), fd_get(Y, YD1, YPs1) ->
 4617                        NYL cis n(X) // (NZU + n(1)) + n(1),
 4618                        NYU cis n(X) // NZL,
 4619                        domains_intersection(YD1, from_to(NYL, NYU), NYD1),
 4620                        fd_put(Y, NYD1, YPs1)
 4621                    ;   true
 4622                    )
 4623                )
 4624            )
 4625        ;   nonvar(Y) ->
 4626            Y =\= 0,
 4627            (   Y =:= 1 -> kill(MState), X = Z
 4628            ;   Y =:= -1 -> kill(MState), Z #= -X
 4629            ;   fd_get(X, XD, XL, XU, XPs),
 4630                (   nonvar(Z) ->
 4631                    kill(MState),
 4632                    (   sign(Z) =:= sign(Y) ->
 4633                        NXL cis max(n(Z)*n(Y), XL),
 4634                        NXU cis min((abs(n(Z))+n(1))*abs(n(Y))-n(1), XU)
 4635                    ;   Z =:= 0 ->
 4636                        NXL cis max(-abs(n(Y)) + n(1), XL),
 4637                        NXU cis min(abs(n(Y)) - n(1), XU)
 4638                    ;   NXL cis max((n(Z)+sign(n(Z)))*n(Y)+n(1), XL),
 4639                        NXU cis min(n(Z)*n(Y), XU)
 4640                    ),
 4641                    update_bounds(X, XD, XPs, XL, XU, NXL, NXU)
 4642                ;   fd_get(Z, ZD, ZPs),
 4643                    domain_contract_less(XD, Y, Contracted),
 4644                    domains_intersection(ZD, Contracted, NZD),
 4645                    fd_put(Z, NZD, ZPs),
 4646                    (   fd_get(X, XD2, XPs2) ->
 4647                        domain_expand_more(NZD, Y, Expanded),
 4648                        domains_intersection(XD2, Expanded, NXD2),
 4649                        fd_put(X, NXD2, XPs2)
 4650                    ;   true
 4651                    )
 4652                )
 4653            )
 4654        ;   nonvar(Z) ->
 4655            fd_get(X, XD, XL, XU, XPs),
 4656            fd_get(Y, _, YL, YU, _),
 4657            (   YL cis_geq n(0), XL cis_geq n(0) ->
 4658                NXL cis max(YL*n(Z), XL),
 4659                NXU cis min(YU*(n(Z)+n(1))-n(1), XU)
 4660            ;   %TODO: cover more cases
 4661                NXL = XL, NXU = XU
 4662            ),
 4663            update_bounds(X, XD, XPs, XL, XU, NXL, NXU)
 4664        ;   (   X == Y -> kill(MState), Z = 1
 4665            ;   fd_get(X, _, XL, XU, _),
 4666                fd_get(Y, _, YL, _, _),
 4667                fd_get(Z, ZD, ZPs),
 4668                NZU cis max(abs(XL), XU),
 4669                NZL cis -NZU,
 4670                domains_intersection(ZD, from_to(NZL,NZU), NZD0),
 4671                (   XL cis_geq n(0), YL cis_geq n(0) ->
 4672                    domain_remove_smaller_than(NZD0, 0, NZD1)
 4673                ;   % TODO: cover more cases
 4674                    NZD1 = NZD0
 4675                ),
 4676                fd_put(Z, NZD1, ZPs)
 4677            )
 4678        ).
 4679
 4680
 4681%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 4682% Y = abs(X)
 4683
 4684run_propagator(pabs(X,Y), MState) :-
 4685        (   nonvar(X) -> kill(MState), Y is abs(X)
 4686        ;   nonvar(Y) ->
 4687            kill(MState),
 4688            Y >= 0,
 4689            YN is -Y,
 4690            X in YN \/ Y
 4691        ;   fd_get(X, XD, XPs),
 4692            fd_get(Y, YD, _),
 4693            domain_negate(YD, YDNegative),
 4694            domains_union(YD, YDNegative, XD1),
 4695            domains_intersection(XD, XD1, XD2),
 4696            fd_put(X, XD2, XPs),
 4697            (   fd_get(Y, YD1, YPs1) ->
 4698                domain_negate(XD2, XD2Neg),
 4699                domains_union(XD2, XD2Neg, YD2),
 4700                domain_remove_smaller_than(YD2, 0, YD3),
 4701                domains_intersection(YD1, YD3, YD4),
 4702                fd_put(Y, YD4, YPs1)
 4703            ;   true
 4704            )
 4705        ).
 4706%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 4707% Z = X mod Y
 4708
 4709run_propagator(pmod(X,Y,Z), MState) :-
 4710        (   nonvar(X) ->
 4711            (   nonvar(Y) -> kill(MState), Y =\= 0, Z is X mod Y
 4712            ;   true
 4713            )
 4714        ;   nonvar(Y) ->
 4715            Y =\= 0,
 4716            (   abs(Y) =:= 1 -> kill(MState), Z = 0
 4717            ;   var(Z) ->
 4718                YP is abs(Y) - 1,
 4719                (   Y > 0, fd_get(X, _, n(XL), n(XU), _) ->
 4720                    (   XL >= 0, XU < Y ->
 4721                        kill(MState), Z = X, ZL = XL, ZU = XU
 4722                    ;   ZL = 0, ZU = YP
 4723                    )
 4724                ;   Y > 0 -> ZL = 0, ZU = YP
 4725                ;   YN is -YP, ZL = YN, ZU = 0
 4726                ),
 4727                (   fd_get(Z, ZD, ZPs) ->
 4728                    domains_intersection(ZD, from_to(n(ZL), n(ZU)), ZD1),
 4729                    domain_infimum(ZD1, n(ZMin)),
 4730                    domain_supremum(ZD1, n(ZMax)),
 4731                    fd_put(Z, ZD1, ZPs)
 4732                ;   ZMin = Z, ZMax = Z
 4733                ),
 4734                (   fd_get(X, XD, XPs), domain_infimum(XD, n(XMin)) ->
 4735                    Z1 is XMin mod Y,
 4736                    (   between(ZMin, ZMax, Z1) -> true
 4737                    ;   Y > 0 ->
 4738                        Next is ((XMin - ZMin + Y - 1) div Y)*Y + ZMin,
 4739                        domain_remove_smaller_than(XD, Next, XD1),
 4740                        fd_put(X, XD1, XPs)
 4741                    ;   neq_num(X, XMin)
 4742                    )
 4743                ;   true
 4744                ),
 4745                (   fd_get(X, XD2, XPs2), domain_supremum(XD2, n(XMax)) ->
 4746                    Z2 is XMax mod Y,
 4747                    (   between(ZMin, ZMax, Z2) -> true
 4748                    ;   Y > 0 ->
 4749                        Prev is ((XMax - ZMin) div Y)*Y + ZMax,
 4750                        domain_remove_greater_than(XD2, Prev, XD3),
 4751                        fd_put(X, XD3, XPs2)
 4752                    ;   neq_num(X, XMax)
 4753                    )
 4754                ;   true
 4755                )
 4756            ;   fd_get(X, XD, XPs),
 4757                % if possible, propagate at the boundaries
 4758                (   domain_infimum(XD, n(Min)) ->
 4759                    (   Min mod Y =:= Z -> true
 4760                    ;   Y > 0 ->
 4761                        Next is ((Min - Z + Y - 1) div Y)*Y + Z,
 4762                        domain_remove_smaller_than(XD, Next, XD1),
 4763                        fd_put(X, XD1, XPs)
 4764                    ;   neq_num(X, Min)
 4765                    )
 4766                ;   true
 4767                ),
 4768                (   fd_get(X, XD2, XPs2) ->
 4769                    (   domain_supremum(XD2, n(Max)) ->
 4770                        (   Max mod Y =:= Z -> true
 4771                        ;   Y > 0 ->
 4772                            Prev is ((Max - Z) div Y)*Y + Z,
 4773                            domain_remove_greater_than(XD2, Prev, XD3),
 4774                            fd_put(X, XD3, XPs2)
 4775                        ;   neq_num(X, Max)
 4776                        )
 4777                    ;   true
 4778                    )
 4779                ;   true
 4780                )
 4781            )
 4782        ;   X == Y -> kill(MState), Z = 0
 4783        ;   true % TODO: propagate more
 4784        ).
 4785
 4786%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 4787% Z = X rem Y
 4788
 4789run_propagator(prem(X,Y,Z), MState) :-
 4790        (   nonvar(X) ->
 4791            (   nonvar(Y) -> kill(MState), Y =\= 0, Z is X rem Y
 4792            ;   U is abs(X),
 4793                fd_get(Y, YD, _),
 4794                (   X >=0, domain_infimum(YD, n(Min)), Min >= 0 -> L = 0
 4795                ;   L is -U
 4796                ),
 4797                Z in L..U
 4798            )
 4799        ;   nonvar(Y) ->
 4800            Y =\= 0,
 4801            (   abs(Y) =:= 1 -> kill(MState), Z = 0
 4802            ;   var(Z) ->
 4803                YP is abs(Y) - 1,
 4804                YN is -YP,
 4805                (   Y > 0, fd_get(X, _, n(XL), n(XU), _) ->
 4806                    (   abs(XL) < Y, XU < Y -> kill(MState), Z = X, ZL = XL
 4807                    ;   XL < 0, abs(XL) < Y -> ZL = XL
 4808                    ;   XL >= 0 -> ZL = 0
 4809                    ;   ZL = YN
 4810                    ),
 4811                    (   XU > 0, XU < Y -> ZU = XU
 4812                    ;   XU < 0 -> ZU = 0
 4813                    ;   ZU = YP
 4814                    )
 4815                ;   ZL = YN, ZU = YP
 4816                ),
 4817                (   fd_get(Z, ZD, ZPs) ->
 4818                    domains_intersection(ZD, from_to(n(ZL), n(ZU)), ZD1),
 4819                    fd_put(Z, ZD1, ZPs)
 4820                ;   ZD1 = from_to(n(Z), n(Z))
 4821                ),
 4822                (   fd_get(X, XD, _), domain_infimum(XD, n(Min)) ->
 4823                    Z1 is Min rem Y,
 4824                    (   domain_contains(ZD1, Z1) -> true
 4825                    ;   neq_num(X, Min)
 4826                    )
 4827                ;   true
 4828                ),
 4829                (   fd_get(X, XD1, _), domain_supremum(XD1, n(Max)) ->
 4830                    Z2 is Max rem Y,
 4831                    (   domain_contains(ZD1, Z2) -> true
 4832                    ;   neq_num(X, Max)
 4833                    )
 4834                ;   true
 4835                )
 4836            ;   fd_get(X, XD1, XPs1),
 4837                % if possible, propagate at the boundaries
 4838                (   domain_infimum(XD1, n(Min)) ->
 4839                    (   Min rem Y =:= Z -> true
 4840                    ;   Y > 0, Min > 0 ->
 4841                        Next is ((Min - Z + Y - 1) div Y)*Y + Z,
 4842                        domain_remove_smaller_than(XD1, Next, XD2),
 4843                        fd_put(X, XD2, XPs1)
 4844                    ;   % TODO: bigger steps in other cases as well
 4845                        neq_num(X, Min)
 4846                    )
 4847                ;   true
 4848                ),
 4849                (   fd_get(X, XD3, XPs3) ->
 4850                    (   domain_supremum(XD3, n(Max)) ->
 4851                        (   Max rem Y =:= Z -> true
 4852                        ;   Y > 0, Max > 0  ->
 4853                            Prev is ((Max - Z) div Y)*Y + Z,
 4854                            domain_remove_greater_than(XD3, Prev, XD4),
 4855                            fd_put(X, XD4, XPs3)
 4856                        ;   % TODO: bigger steps in other cases as well
 4857                            neq_num(X, Max)
 4858                        )
 4859                    ;   true
 4860                    )
 4861                ;   true
 4862                )
 4863            )
 4864        ;   X == Y -> kill(MState), Z = 0
 4865        ;   fd_get(Z, ZD, ZPs) ->
 4866            fd_get(Y, _, YInf, YSup, _),
 4867            fd_get(X, _, XInf, XSup, _),
 4868            M cis max(abs(YInf),YSup),
 4869            (   XInf cis_geq n(0) -> Inf0 = n(0)
 4870            ;   Inf0 = XInf
 4871            ),
 4872            (   XSup cis_leq n(0) -> Sup0 = n(0)
 4873            ;   Sup0 = XSup
 4874            ),
 4875            NInf cis max(max(Inf0, -M + n(1)), min(XInf,-XSup)),
 4876            NSup cis min(min(Sup0, M - n(1)), max(abs(XInf),XSup)),
 4877            domains_intersection(ZD, from_to(NInf,NSup), ZD1),
 4878            fd_put(Z, ZD1, ZPs)
 4879        ;   true % TODO: propagate more
 4880        ).
 4881
 4882%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 4883% Z = max(X,Y)
 4884
 4885run_propagator(pmax(X,Y,Z), MState) :-
 4886        (   nonvar(X) ->
 4887            (   nonvar(Y) -> kill(MState), Z is max(X,Y)
 4888            ;   nonvar(Z) ->
 4889                (   Z =:= X -> kill(MState), X #>= Y
 4890                ;   Z > X -> Z = Y
 4891                ;   false % Z < X
 4892                )
 4893            ;   fd_get(Y, _, YInf, YSup, _),
 4894                (   YInf cis_gt n(X) -> Z = Y
 4895                ;   YSup cis_lt n(X) -> Z = X
 4896                ;   YSup = n(M) ->
 4897                    fd_get(Z, ZD, ZPs),
 4898                    domain_remove_greater_than(ZD, M, ZD1),
 4899                    fd_put(Z, ZD1, ZPs)
 4900                ;   true
 4901                )
 4902            )
 4903        ;   nonvar(Y) -> run_propagator(pmax(Y,X,Z), MState)
 4904        ;   fd_get(Z, ZD, ZPs) ->
 4905            fd_get(X, _, XInf, XSup, _),
 4906            fd_get(Y, _, YInf, YSup, _),
 4907            (   YInf cis_gt YSup -> kill(MState), Z = Y
 4908            ;   YSup cis_lt XInf -> kill(MState), Z = X
 4909            ;   n(M) cis max(XSup, YSup) ->
 4910                domain_remove_greater_than(ZD, M, ZD1),
 4911                fd_put(Z, ZD1, ZPs)
 4912            ;   true
 4913            )
 4914        ;   true
 4915        ).
 4916
 4917%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 4918% Z = min(X,Y)
 4919
 4920run_propagator(pmin(X,Y,Z), MState) :-
 4921        (   nonvar(X) ->
 4922            (   nonvar(Y) -> kill(MState), Z is min(X,Y)
 4923            ;   nonvar(Z) ->
 4924                (   Z =:= X -> kill(MState), X #=< Y
 4925                ;   Z < X -> Z = Y
 4926                ;   false % Z > X
 4927                )
 4928            ;   fd_get(Y, _, YInf, YSup, _),
 4929                (   YSup cis_lt n(X) -> Z = Y
 4930                ;   YInf cis_gt n(X) -> Z = X
 4931                ;   YInf = n(M) ->
 4932                    fd_get(Z, ZD, ZPs),
 4933                    domain_remove_smaller_than(ZD, M, ZD1),
 4934                    fd_put(Z, ZD1, ZPs)
 4935                ;   true
 4936                )
 4937            )
 4938        ;   nonvar(Y) -> run_propagator(pmin(Y,X,Z), MState)
 4939        ;   fd_get(Z, ZD, ZPs) ->
 4940            fd_get(X, _, XInf, XSup, _),
 4941            fd_get(Y, _, YInf, YSup, _),
 4942            (   YSup cis_lt YInf -> kill(MState), Z = Y
 4943            ;   YInf cis_gt XSup -> kill(MState), Z = X
 4944            ;   n(M) cis min(XInf, YInf) ->
 4945                domain_remove_smaller_than(ZD, M, ZD1),
 4946                fd_put(Z, ZD1, ZPs)
 4947            ;   true
 4948            )
 4949        ;   true
 4950        ).
 4951%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 4952% Z = X ^ Y
 4953
 4954run_propagator(pexp(X,Y,Z), MState) :-
 4955        (   X == 1 -> kill(MState), Z = 1
 4956        ;   X == 0 -> kill(MState), Z in 0..1, Z #<==> Y #= 0
 4957        ;   Y == 0 -> kill(MState), Z = 1
 4958        ;   Y == 1 -> kill(MState), Z = X
 4959        ;   nonvar(X) ->
 4960            (   nonvar(Y) ->
 4961                (   Y >= 0 -> true ; X =:= -1 ),
 4962                kill(MState),
 4963                Z is X^Y
 4964            ;   nonvar(Z) ->
 4965                (   Z > 1 ->
 4966                    kill(MState),
 4967                    integer_log_b(Z, X, 1, Y)
 4968                ;   true
 4969                )
 4970            ;   fd_get(Y, _, YL, YU, _),
 4971                fd_get(Z, ZD, ZPs),
 4972                (   X > 0, YL cis_geq n(0) ->
 4973                    NZL cis n(X)^YL,
 4974                    NZU cis n(X)^YU,
 4975                    domains_intersection(ZD, from_to(NZL,NZU), NZD),
 4976                    fd_put(Z, NZD, ZPs)
 4977                ;   true
 4978                ),
 4979                (   X > 0,
 4980                    fd_get(Z, _, _, n(ZMax), _),
 4981                    ZMax > 0 ->
 4982                    floor_integer_log_b(ZMax, X, 1, YCeil),
 4983                    Y in inf..YCeil
 4984                ;   true
 4985                )
 4986            )
 4987        ;   nonvar(Z) ->
 4988            (   nonvar(Y) ->
 4989                integer_kth_root(Z, Y, R),
 4990                kill(MState),
 4991                (   even(Y) ->
 4992                    N is -R,
 4993                    X in N \/ R
 4994                ;   X = R
 4995                )
 4996            ;   fd_get(X, _, n(NXL), _, _), NXL > 1 ->
 4997                (   Z > 1, between(NXL, Z, Exp), NXL^Exp > Z ->
 4998                    Exp1 is Exp - 1,
 4999                    fd_get(Y, YD, YPs),
 5000                    domains_intersection(YD, from_to(n(1),n(Exp1)), YD1),
 5001                    fd_put(Y, YD1, YPs),
 5002                    (   fd_get(X, XD, XPs) ->
 5003                        domain_infimum(YD1, n(YL)),
 5004                        integer_kth_root_leq(Z, YL, RU),
 5005                        domains_intersection(XD, from_to(n(NXL),n(RU)), XD1),
 5006                        fd_put(X, XD1, XPs)
 5007                    ;   true
 5008                    )
 5009                ;   true
 5010                )
 5011            ;   true
 5012            )
 5013        ;   nonvar(Y), Y > 0 ->
 5014            (   even(Y) ->
 5015                geq(Z, 0)
 5016            ;   true
 5017            ),
 5018            (   fd_get(X, XD, XL, XU, _), fd_get(Z, ZD, ZL, ZU, ZPs) ->
 5019                (   domain_contains(ZD, 0) -> XD1 = XD
 5020                ;   domain_remove(XD, 0, XD1)
 5021                ),
 5022                (   domain_contains(XD, 0) -> ZD1 = ZD
 5023                ;   domain_remove(ZD, 0, ZD1)
 5024                ),
 5025                (   even(Y) ->
 5026                    (   XL cis_geq n(0) ->
 5027                        NZL cis XL^n(Y)
 5028                    ;   XU cis_leq n(0) ->
 5029                        NZL cis XU^n(Y)
 5030                    ;   NZL = n(0)
 5031                    ),
 5032                    NZU cis max(abs(XL),abs(XU))^n(Y),
 5033                    domains_intersection(ZD1, from_to(NZL,NZU), ZD2)
 5034                ;   (   finite(XL) ->
 5035                        NZL cis XL^n(Y),
 5036                        NZU cis XU^n(Y),
 5037                        domains_intersection(ZD1, from_to(NZL,NZU), ZD2)
 5038                    ;   ZD2 = ZD1
 5039                    )
 5040                ),
 5041                fd_put(Z, ZD2, ZPs),
 5042                (   even(Y), ZU = n(Num) ->
 5043                    integer_kth_root_leq(Num, Y, RU),
 5044                    (   XL cis_geq n(0), ZL = n(Num1) ->
 5045                        integer_kth_root_leq(Num1, Y, RL0),
 5046                        (   RL0^Y < Num1 -> RL is RL0 + 1
 5047                        ;   RL = RL0
 5048                        )
 5049                    ;   RL is -RU
 5050                    ),
 5051                    RL =< RU,
 5052                    NXD = from_to(n(RL),n(RU))
 5053                ;   odd(Y), ZL cis_geq n(0), ZU = n(Num) ->
 5054                    integer_kth_root_leq(Num, Y, RU),
 5055                    ZL = n(Num1),
 5056                    integer_kth_root_leq(Num1, Y, RL0),
 5057                    (   RL0^Y < Num1 -> RL is RL0 + 1
 5058                    ;   RL = RL0
 5059                    ),
 5060                    RL =< RU,
 5061                    NXD = from_to(n(RL),n(RU))
 5062                ;   NXD = XD1   % TODO: propagate more
 5063                ),
 5064                (   fd_get(X, XD2, XPs) ->
 5065                    domains_intersection(XD2, XD1, XD3),
 5066                    domains_intersection(XD3, NXD, XD4),
 5067                    fd_put(X, XD4, XPs)
 5068                ;   true
 5069                )
 5070            ;   true
 5071            )
 5072        ;   fd_get(X, _, XL, _, _),
 5073            XL cis_gt n(0),
 5074            fd_get(Y, _, YL, _, _),
 5075            YL cis_gt n(0),
 5076            fd_get(Z, ZD, ZPs) ->
 5077            n(NZL) cis XL^YL,
 5078            domain_remove_smaller_than(ZD, NZL, ZD1),
 5079            fd_put(Z, ZD1, ZPs)
 5080        ;   true
 5081        ).
 5082
 5083%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 5084run_propagator(pzcompare(Order, A, B), MState) :-
 5085        (   A == B -> kill(MState), Order = (=)
 5086        ;   (   nonvar(A) ->
 5087                (   nonvar(B) ->
 5088                    kill(MState),
 5089                    (   A > B -> Order = (>)
 5090                    ;   Order = (<)
 5091                    )
 5092                ;   fd_get(B, _, BL, BU, _),
 5093                    (   BL cis_gt n(A) -> kill(MState), Order = (<)
 5094                    ;   BU cis_lt n(A) -> kill(MState), Order = (>)
 5095                    ;   true
 5096                    )
 5097                )
 5098            ;   nonvar(B) ->
 5099                fd_get(A, _, AL, AU, _),
 5100                (   AL cis_gt n(B) -> kill(MState), Order = (>)
 5101                ;   AU cis_lt n(B) -> kill(MState), Order = (<)
 5102                ;   true
 5103                )
 5104            ;   fd_get(A, _, AL, AU, _),
 5105                fd_get(B, _, BL, BU, _),
 5106                (   AL cis_gt BU -> kill(MState), Order = (>)
 5107                ;   AU cis_lt BL -> kill(MState), Order = (<)
 5108                ;   true
 5109                )
 5110            )
 5111        ).
 5112
 5113%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 5114
 5115% reified constraints
 5116
 5117%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 5118
 5119run_propagator(reified_in(V,Dom,B), MState) :-
 5120        (   integer(V) ->
 5121            kill(MState),
 5122            (   domain_contains(Dom, V) -> B = 1
 5123            ;   B = 0
 5124            )
 5125        ;   B == 1 -> kill(MState), domain(V, Dom)
 5126        ;   B == 0 -> kill(MState), domain_complement(Dom, C), domain(V, C)
 5127        ;   fd_get(V, VD, _),
 5128            (   domains_intersection(VD, Dom, I) ->
 5129                (   I == VD -> kill(MState), B = 1
 5130                ;   true
 5131                )
 5132            ;   kill(MState), B = 0
 5133            )
 5134        ).
 5135
 5136run_propagator(reified_tuple_in(Tuple, R, B), MState) :-
 5137        get_attr(R, clpfd_relation, Relation),
 5138        (   B == 1 -> kill(MState), tuples_in([Tuple], Relation)
 5139        ;   (   ground(Tuple) ->
 5140                kill(MState),
 5141                (   memberchk(Tuple, Relation) -> B = 1
 5142                ;   B = 0
 5143                )
 5144            ;