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