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