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