1/*  Part of SWI-Prolog
    2
    3    Author:        Markus Triska
    4    E-mail:        triska@metalevel.at
    5    WWW:           http://www.swi-prolog.org
    6    Copyright (C): 2014-2018 Markus Triska
    7    All rights reserved.
    8
    9    Redistribution and use in source and binary forms, with or without
   10    modification, are permitted provided that the following conditions
   11    are met:
   12
   13    1. Redistributions of source code must retain the above copyright
   14       notice, this list of conditions and the following disclaimer.
   15
   16    2. Redistributions in binary form must reproduce the above copyright
   17       notice, this list of conditions and the following disclaimer in
   18       the documentation and/or other materials provided with the
   19       distribution.
   20
   21    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
   22    "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
   23    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
   24    FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
   25    COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
   26    INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
   27    BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
   28    LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
   29    CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
   30    LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
   31    ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
   32    POSSIBILITY OF SUCH DAMAGE.
   33*/
   34
   35/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
   36   CLP(B): Constraint Logic Programming over Boolean variables.
   37- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
   38
   39:- module(clpb, [
   40                 op(300, fy, ~),
   41                 op(500, yfx, #),
   42                 sat/1,
   43                 taut/2,
   44                 labeling/1,
   45                 sat_count/2,
   46                 weighted_maximum/3,
   47                 random_labeling/2
   48                ]).   49
   50:- use_module(library(error)).   51:- use_module(library(assoc)).   52:- use_module(library(apply_macros)).   53
   54:- create_prolog_flag(clpb_monotonic, false, []).   55:- create_prolog_flag(clpb_residuals, default, []).   56
   57/** <module> CLP(B): Constraint Logic Programming over Boolean Variables
   58
   59## Introduction                        {#clpb-intro}
   60
   61This library provides CLP(B), Constraint Logic Programming over
   62Boolean variables. It can be used to model and solve combinatorial
   63problems such as verification, allocation and covering tasks.
   64
   65CLP(B) is an instance of the general [CLP(_X_) scheme](<#clp>),
   66extending logic programming with reasoning over specialised domains.
   67
   68The implementation is based on reduced and ordered Binary Decision
   69Diagrams (BDDs).
   70
   71Benchmarks and usage examples of this library are available from:
   72[__https://www.metalevel.at/clpb/__](https://www.metalevel.at/clpb/)
   73
   74We recommend the following references for citing this library in
   75scientific publications:
   76
   77==
   78@inproceedings{Triska2016,
   79  author    = "Markus Triska",
   80  title     = "The {Boolean} Constraint Solver of {SWI-Prolog}:
   81               System Description",
   82  booktitle = "FLOPS",
   83  series    = "LNCS",
   84  volume    = 9613,
   85  year      = 2016,
   86  pages     = "45--61"
   87}
   88
   89@article{Triska2018,
   90  title = "Boolean constraints in {SWI-Prolog}:
   91           A comprehensive system description",
   92  journal = "Science of Computer Programming",
   93  volume = "164",
   94  pages = "98 - 115",
   95  year = "2018",
   96  note = "Special issue of selected papers from FLOPS 2016",
   97  issn = "0167-6423",
   98  doi = "https://doi.org/10.1016/j.scico.2018.02.001",
   99  url = "http://www.sciencedirect.com/science/article/pii/S0167642318300273",
  100  author = "Markus Triska",
  101  keywords = "CLP(B), Boolean unification, Decision diagrams, BDD"
  102}
  103==
  104
  105These papers are available from
  106[https://www.metalevel.at/swiclpb.pdf](https://www.metalevel.at/swiclpb.pdf)
  107and
  108[https://www.metalevel.at/boolean.pdf](https://www.metalevel.at/boolean.pdf)
  109respectively.
  110
  111## Boolean expressions {#clpb-exprs}
  112
  113A _Boolean expression_ is one of:
  114
  115    | `0`                | false                                |
  116    | `1`                | true                                 |
  117    | _variable_         | unknown truth value                  |
  118    | _atom_             | universally quantified variable      |
  119    | ~ _Expr_           | logical NOT                          |
  120    | _Expr_ + _Expr_    | logical OR                           |
  121    | _Expr_ * _Expr_    | logical AND                          |
  122    | _Expr_ # _Expr_    | exclusive OR                         |
  123    | _Var_ ^ _Expr_     | existential quantification           |
  124    | _Expr_ =:= _Expr_  | equality                             |
  125    | _Expr_ =\= _Expr_  | disequality (same as #)              |
  126    | _Expr_ =< _Expr_   | less or equal (implication)          |
  127    | _Expr_ >= _Expr_   | greater or equal                     |
  128    | _Expr_ < _Expr_    | less than                            |
  129    | _Expr_ > _Expr_    | greater than                         |
  130    | card(Is,Exprs)     | cardinality constraint (_see below_) |
  131    | `+(Exprs)`         | n-fold disjunction (_see below_)     |
  132    | `*(Exprs)`         | n-fold conjunction (_see below_)     |
  133
  134where _Expr_ again denotes a Boolean expression.
  135
  136The Boolean expression card(Is,Exprs) is true iff the number of true
  137expressions in the list `Exprs` is a member of the list `Is` of
  138integers and integer ranges of the form `From-To`. For example, to
  139state that precisely two of the three variables `X`, `Y` and `Z` are
  140`true`, you can use `sat(card([2],[X,Y,Z]))`.
  141
  142`+(Exprs)` and `*(Exprs)` denote, respectively, the disjunction and
  143conjunction of all elements in the list `Exprs` of Boolean
  144expressions.
  145
  146Atoms denote parametric values that are universally quantified. All
  147universal quantifiers appear implicitly in front of the entire
  148expression. In residual goals, universally quantified variables always
  149appear on the right-hand side of equations. Therefore, they can be
  150used to express functional dependencies on input variables.
  151
  152## Interface predicates   {#clpb-interface}
  153
  154The most frequently used CLP(B) predicates are:
  155
  156    * sat(+Expr)
  157      True iff the Boolean expression Expr is satisfiable.
  158
  159    * taut(+Expr, -T)
  160      If Expr is a tautology with respect to the posted constraints, succeeds
  161      with *T = 1*. If Expr cannot be satisfied, succeeds with *T = 0*.
  162      Otherwise, it fails.
  163
  164    * labeling(+Vs)
  165      Assigns truth values to the variables Vs such that all constraints
  166      are satisfied.
  167
  168The unification of a CLP(B) variable _X_ with a term _T_ is equivalent
  169to posting the constraint sat(X=:=T).
  170
  171## Examples                            {#clpb-examples}
  172
  173Here is an example session with a few queries and their answers:
  174
  175==
  176?- use_module(library(clpb)).
  177true.
  178
  179?- sat(X*Y).
  180X = Y, Y = 1.
  181
  182?- sat(X * ~X).
  183false.
  184
  185?- taut(X * ~X, T).
  186T = 0,
  187sat(X=:=X).
  188
  189?- sat(X^Y^(X+Y)).
  190sat(X=:=X),
  191sat(Y=:=Y).
  192
  193?- sat(X*Y + X*Z), labeling([X,Y,Z]).
  194X = Z, Z = 1, Y = 0 ;
  195X = Y, Y = 1, Z = 0 ;
  196X = Y, Y = Z, Z = 1.
  197
  198?- sat(X =< Y), sat(Y =< Z), taut(X =< Z, T).
  199T = 1,
  200sat(X=:=X*Y),
  201sat(Y=:=Y*Z).
  202
  203?- sat(1#X#a#b).
  204sat(X=:=a#b).
  205==
  206
  207The pending residual goals constrain remaining variables to Boolean
  208expressions and are declaratively equivalent to the original query.
  209The last example illustrates that when applicable, remaining variables
  210are expressed as functions of universally quantified variables.
  211
  212## Obtaining BDDs {#clpb-residual-goals}
  213
  214By default, CLP(B) residual goals appear in (approximately) algebraic
  215normal form (ANF). This projection is often computationally expensive.
  216You can set the Prolog flag `clpb_residuals` to the value `bdd` to see
  217the BDD representation of all constraints. This results in faster
  218projection to residual goals, and is also useful for learning more
  219about BDDs. For example:
  220
  221==
  222?- set_prolog_flag(clpb_residuals, bdd).
  223true.
  224
  225?- sat(X#Y).
  226node(3)- (v(X, 0)->node(2);node(1)),
  227node(1)- (v(Y, 1)->true;false),
  228node(2)- (v(Y, 1)->false;true).
  229==
  230
  231Note that this representation cannot be pasted back on the toplevel,
  232and its details are subject to change. Use copy_term/3 to obtain
  233such answers as Prolog terms.
  234
  235The variable order of the BDD is determined by the order in which the
  236variables first appear in constraints. To obtain different orders,
  237you can for example use:
  238
  239==
  240?- sat(+[1,Y,X]), sat(X#Y).
  241node(3)- (v(Y, 0)->node(2);node(1)),
  242node(1)- (v(X, 1)->true;false),
  243node(2)- (v(X, 1)->false;true).
  244==
  245
  246## Enabling monotonic CLP(B) {#clpb-monotonic}
  247
  248In the default execution mode, CLP(B) constraints are _not_ monotonic.
  249This means that _adding_ constraints can yield new solutions. For
  250example:
  251
  252==
  253?-          sat(X=:=1), X = 1+0.
  254false.
  255
  256?- X = 1+0, sat(X=:=1), X = 1+0.
  257X = 1+0.
  258==
  259
  260This behaviour is highly problematic from a logical point of view, and
  261it may render [**declarative
  262debugging**](https://www.metalevel.at/prolog/debugging)
  263techniques inapplicable.
  264
  265Set the flag `clpb_monotonic` to `true` to make CLP(B) *monotonic*. If
  266this mode is enabled, then you must wrap CLP(B) variables with the
  267functor `v/1`. For example:
  268
  269==
  270?- set_prolog_flag(clpb_monotonic, true).
  271true.
  272
  273?- sat(v(X)=:=1#1).
  274X = 0.
  275==
  276
  277@author [Markus Triska](https://www.metalevel.at)
  278*/
  279
  280
  281/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  282   Each CLP(B) variable belongs to exactly one BDD. Each CLP(B)
  283   variable gets an attribute (in module "clpb") of the form:
  284
  285        index_root(Index,Root)
  286
  287   where Index is the variable's unique integer index, and Root is the
  288   root of the BDD that the variable belongs to.
  289
  290   Each CLP(B) variable also gets an attribute in module clpb_hash: an
  291   association table node(LID,HID) -> Node, to keep the BDD reduced.
  292   The association table of each variable must be rebuilt on occasion
  293   to remove nodes that are no longer reachable. We rebuild the
  294   association tables of involved variables after BDDs are merged to
  295   build a new root. This only serves to reclaim memory: Keeping a
  296   node in a local table even when it no longer occurs in any BDD does
  297   not affect the solver's correctness. However, apply_shortcut/4
  298   relies on the invariant that every node that occurs in the relevant
  299   BDDs is also registered in the table of its branching variable.
  300
  301   A root is a logical variable with a single attribute ("clpb_bdd")
  302   of the form:
  303
  304        Sat-BDD
  305
  306   where Sat is the SAT formula (in original form) that corresponds to
  307   BDD. Sat is necessary to rebuild the BDD after variable aliasing,
  308   and to project all remaining constraints to a list of sat/1 goals.
  309
  310   Finally, a BDD is either:
  311
  312      *)  The integers 0 or 1, denoting false and true, respectively, or
  313      *)  A node of the form
  314
  315           node(ID, Var, Low, High, Aux)
  316               Where ID is the node's unique integer ID, Var is the
  317               node's branching variable, and Low and High are the
  318               node's low (Var = 0) and high (Var = 1) children. Aux
  319               is a free variable, one for each node, that can be used
  320               to attach attributes and store intermediate results.
  321
  322   Variable aliasing is treated as a conjunction of corresponding SAT
  323   formulae.
  324
  325   You should think of CLP(B) as a potentially vast collection of BDDs
  326   that can range from small to gigantic in size, and which can merge.
  327- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  328
  329/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  330   Type checking.
  331- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  332
  333is_sat(V)     :- var(V), !, non_monotonic(V).
  334is_sat(v(V))  :- var(V), !.
  335is_sat(v(I))  :- integer(I), between(0, 1, I).
  336is_sat(I)     :- integer(I), between(0, 1, I).
  337is_sat(A)     :- atom(A).
  338is_sat(~A)    :- is_sat(A).
  339is_sat(A*B)   :- is_sat(A), is_sat(B).
  340is_sat(A+B)   :- is_sat(A), is_sat(B).
  341is_sat(A#B)   :- is_sat(A), is_sat(B).
  342is_sat(A=:=B) :- is_sat(A), is_sat(B).
  343is_sat(A=\=B) :- is_sat(A), is_sat(B).
  344is_sat(A=<B)  :- is_sat(A), is_sat(B).
  345is_sat(A>=B)  :- is_sat(A), is_sat(B).
  346is_sat(A<B)   :- is_sat(A), is_sat(B).
  347is_sat(A>B)   :- is_sat(A), is_sat(B).
  348is_sat(+(Ls)) :- must_be(list, Ls), maplist(is_sat, Ls).
  349is_sat(*(Ls)) :- must_be(list, Ls), maplist(is_sat, Ls).
  350is_sat(X^F)   :- var(X), is_sat(F).
  351is_sat(card(Is,Fs)) :-
  352        must_be(list(ground), Is),
  353        must_be(list, Fs),
  354        maplist(is_sat, Fs).
  355
  356non_monotonic(X) :-
  357        (   var_index(X, _) ->
  358            % OK: already constrained to a CLP(B) variable
  359            true
  360        ;   current_prolog_flag(clpb_monotonic, true) ->
  361            instantiation_error(X)
  362        ;   true
  363        ).
  364
  365/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  366   Rewriting to canonical expressions.
  367   Atoms are converted to variables with a special attribute.
  368   A global lookup table maintains the correspondence between atoms and
  369   their variables throughout different sat/1 goals.
  370- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  371
  372% elementary
  373sat_rewrite(V, V)       :- var(V), !.
  374sat_rewrite(I, I)       :- integer(I), !.
  375sat_rewrite(A, V)       :- atom(A), !, clpb_atom_var(A, V).
  376sat_rewrite(v(V), V).
  377sat_rewrite(P0*Q0, P*Q) :- sat_rewrite(P0, P), sat_rewrite(Q0, Q).
  378sat_rewrite(P0+Q0, P+Q) :- sat_rewrite(P0, P), sat_rewrite(Q0, Q).
  379sat_rewrite(P0#Q0, P#Q) :- sat_rewrite(P0, P), sat_rewrite(Q0, Q).
  380sat_rewrite(X^F0, X^F)  :- sat_rewrite(F0, F).
  381sat_rewrite(card(Is,Fs0), card(Is,Fs)) :-
  382        maplist(sat_rewrite, Fs0, Fs).
  383% synonyms
  384sat_rewrite(~P, R)      :- sat_rewrite(1 # P, R).
  385sat_rewrite(P =:= Q, R) :- sat_rewrite(~P # Q, R).
  386sat_rewrite(P =\= Q, R) :- sat_rewrite(P # Q, R).
  387sat_rewrite(P =< Q, R)  :- sat_rewrite(~P + Q, R).
  388sat_rewrite(P >= Q, R)  :- sat_rewrite(Q =< P, R).
  389sat_rewrite(P < Q, R)   :- sat_rewrite(~P * Q, R).
  390sat_rewrite(P > Q, R)   :- sat_rewrite(Q < P, R).
  391sat_rewrite(+(Ls), R)   :- foldl(or, Ls, 0, F), sat_rewrite(F, R).
  392sat_rewrite(*(Ls), R)   :- foldl(and, Ls, 1, F), sat_rewrite(F, R).
  393
  394or(A, B, B + A).
  395
  396and(A, B, B * A).
  397
  398must_be_sat(Sat) :-
  399        must_be(acyclic, Sat),
  400        (   is_sat(Sat) -> true
  401        ;   no_truth_value(Sat)
  402        ).
  403
  404no_truth_value(Term) :- domain_error(clpb_expr, Term).
  405
  406parse_sat(Sat0, Sat) :-
  407        must_be_sat(Sat0),
  408        sat_rewrite(Sat0, Sat),
  409        term_variables(Sat, Vs),
  410        maplist(enumerate_variable, Vs).
  411
  412enumerate_variable(V) :-
  413        (   var_index_root(V, _, _) -> true
  414        ;   clpb_next_id('$clpb_next_var', Index),
  415            put_attr(V, clpb, index_root(Index,_)),
  416            put_empty_hash(V)
  417        ).
  418
  419var_index(V, I) :- var_index_root(V, I, _).
  420
  421var_index_root(V, I, Root) :- get_attr(V, clpb, index_root(I,Root)).
  422
  423put_empty_hash(V) :-
  424        empty_assoc(H0),
  425        put_attr(V, clpb_hash, H0).
  426
  427sat_roots(Sat, Roots) :-
  428        term_variables(Sat, Vs),
  429        maplist(var_index_root, Vs, _, Roots0),
  430        term_variables(Roots0, Roots).
  431
  432%% sat(+Expr) is semidet.
  433%
  434% True iff Expr is a satisfiable Boolean expression.
  435
  436sat(Sat0) :-
  437        (   phrase(sat_ands(Sat0), Ands), Ands = [_,_|_] ->
  438            maplist(sat, Ands)
  439        ;   parse_sat(Sat0, Sat),
  440            sat_bdd(Sat, BDD),
  441            sat_roots(Sat, Roots),
  442            roots_and(Roots, Sat0-BDD, And-BDD1),
  443            maplist(del_bdd, Roots),
  444            maplist(=(Root), Roots),
  445            root_put_formula_bdd(Root, And, BDD1),
  446            is_bdd(BDD1),
  447            satisfiable_bdd(BDD1)
  448        ).
  449
  450/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  451   Posting many small sat/1 constraints is better than posting a huge
  452   conjunction (or negated disjunction), because unneeded nodes are
  453   removed from node tables after BDDs are merged. This is not
  454   possible in sat_bdd/2 because the nodes may occur in other BDDs. A
  455   better version of sat_bdd/2 or a proper implementation of a unique
  456   table including garbage collection would make this obsolete and
  457   also improve taut/2 and sat_count/2 in such cases.
  458- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  459
  460sat_ands(X) -->
  461        (   { var(X) } -> [X]
  462        ;   { X = (A*B) } -> sat_ands(A), sat_ands(B)
  463        ;   { X = *(Ls) } -> sat_ands_(Ls)
  464        ;   { X = ~Y } -> not_ors(Y)
  465        ;   [X]
  466        ).
  467
  468sat_ands_([]) --> [].
  469sat_ands_([L|Ls]) --> [L], sat_ands_(Ls).
  470
  471not_ors(X) -->
  472        (   { var(X) } -> [~X]
  473        ;   { X = (A+B) } -> not_ors(A), not_ors(B)
  474        ;   { X = +(Ls) } -> not_ors_(Ls)
  475        ;   [~X]
  476        ).
  477
  478not_ors_([]) --> [].
  479not_ors_([L|Ls]) --> [~L], not_ors_(Ls).
  480
  481del_bdd(Root) :- del_attr(Root, clpb_bdd).
  482
  483root_get_formula_bdd(Root, F, BDD) :- get_attr(Root, clpb_bdd, F-BDD).
  484
  485root_put_formula_bdd(Root, F, BDD) :- put_attr(Root, clpb_bdd, F-BDD).
  486
  487roots_and(Roots, Sat0-BDD0, Sat-BDD) :-
  488        foldl(root_and, Roots, Sat0-BDD0, Sat-BDD),
  489        rebuild_hashes(BDD).
  490
  491root_and(Root, Sat0-BDD0, Sat-BDD) :-
  492        (   root_get_formula_bdd(Root, F, B) ->
  493            Sat = F*Sat0,
  494            bdd_and(B, BDD0, BDD)
  495        ;   Sat = Sat0,
  496            BDD = BDD0
  497        ).
  498
  499bdd_and(NA, NB, And) :-
  500        apply(*, NA, NB, And),
  501        is_bdd(And).
  502
  503%% taut(+Expr, -T) is semidet
  504%
  505% Tautology check. Succeeds with T = 0 if the Boolean expression Expr
  506% cannot be satisfied, and with T = 1 if Expr is always true with
  507% respect to the current constraints. Fails otherwise.
  508
  509taut(Sat0, T) :-
  510        parse_sat(Sat0, Sat),
  511        (   T = 0, \+ sat(Sat) -> true
  512        ;   T = 1, tautology(Sat) -> true
  513        ;   false
  514        ).
  515
  516/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  517   The algebraic equivalence: tautology(F) <=> \+ sat(~F) does NOT
  518   hold in CLP(B) because the quantifiers of universally quantified
  519   variables always implicitly appear in front of the *entire*
  520   expression. Thus we have for example: X+a is not a tautology, but
  521   ~(X+a), meaning forall(a, ~(X+a)), is unsatisfiable:
  522
  523      sat(~(X+a)) = sat(~X * ~a) = sat(~X), sat(~a) = X=0, false
  524
  525   The actual negation of X+a, namely ~forall(A,X+A), in terms of
  526   CLP(B): ~ ~exists(A, ~(X+A)), is of course satisfiable:
  527
  528      ?- sat(~ ~A^ ~(X+A)).
  529      %@ X = 0,
  530      %@ sat(A=:=A).
  531
  532   Instead, of such rewriting, we test whether the BDD of the negated
  533   formula is 0. Critically, this avoids constraint propagation.
  534- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  535
  536tautology(Sat) :-
  537        (   phrase(sat_ands(Sat), Ands), Ands = [_,_|_] ->
  538            maplist(tautology, Ands)
  539        ;   catch((sat_roots(Sat, Roots),
  540                   roots_and(Roots, _-1, _-Ands),
  541                   sat_bdd(1#Sat, BDD),
  542                   bdd_and(BDD, Ands, B),
  543                   B == 0,
  544                   % reset all attributes
  545                   throw(tautology)),
  546                  tautology,
  547                  true)
  548        ).
  549
  550satisfiable_bdd(BDD) :-
  551        (   BDD == 0 -> false
  552        ;   BDD == 1 -> true
  553        ;   (   bdd_nodes(var_unbound, BDD, Nodes) ->
  554                bdd_variables_classification(BDD, Nodes, Classes),
  555                partition(var_class, Classes, Eqs, Bs, Ds),
  556                domain_consistency(Eqs, Goal),
  557                aliasing_consistency(Bs, Ds, Goals),
  558                maplist(unification, [Goal|Goals])
  559            ;   % if any variable is instantiated, we do not perform
  560                % any propagation for now
  561                true
  562            )
  563        ).
  564
  565var_class(_=_, <).
  566var_class(further_branching(_,_), =).
  567var_class(negative_decisive(_), >).
  568
  569unification(true).
  570unification(A=B) :- A = B.      % safe_goal/1 detects safety of this call
  571
  572var_unbound(Node) :-
  573        node_var_low_high(Node, Var, _, _),
  574        var(Var).
  575
  576universal_var(Var) :- get_attr(Var, clpb_atom, _).
  577
  578/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  579   By aliasing consistency, we mean that all unifications X=Y, where
  580   taut(X=:=Y, 1) holds, are posted.
  581
  582   To detect this, we distinguish two kinds of variables among those
  583   variables that are not skipped in any branch: further-branching and
  584   negative-decisive. X is negative-decisive iff every node where X
  585   appears as a branching variable has 0 as one of its children. X is
  586   further-branching iff 1 is not a direct child of any node where X
  587   appears as a branching variable.
  588
  589   Any potential aliasing must involve one further-branching, and one
  590   negative-decisive variable. X=Y must hold if, for each low branch
  591   of nodes with X as branching variable, Y has high branch 0, and for
  592   each high branch of nodes involving X, Y has low branch 0.
  593- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  594
  595aliasing_consistency(Bs, Ds, Goals) :-
  596        phrase(aliasings(Bs, Ds), Goals).
  597
  598aliasings([], _) --> [].
  599aliasings([further_branching(B,Nodes)|Bs], Ds) -->
  600        { var_index(B, BI) },
  601        aliasings_(Ds, B, BI, Nodes),
  602        aliasings(Bs, Ds).
  603
  604aliasings_([], _, _, _) --> [].
  605aliasings_([negative_decisive(D)|Ds], B, BI, Nodes) -->
  606        { var_index(D, DI) },
  607        (   { DI > BI,
  608              always_false(high, DI, Nodes),
  609              always_false(low, DI, Nodes),
  610              var_or_atom(D, DA), var_or_atom(B, BA) } ->
  611            [DA=BA]
  612        ;   []
  613        ),
  614        aliasings_(Ds, B, BI, Nodes).
  615
  616var_or_atom(Var, VA) :-
  617        (   get_attr(Var, clpb_atom, VA) -> true
  618        ;   VA = Var
  619        ).
  620
  621always_false(Which, DI, Nodes) :-
  622        phrase(nodes_always_false(Nodes, Which, DI), Opposites),
  623        maplist(with_aux(unvisit), Opposites).
  624
  625nodes_always_false([], _, _) --> [].
  626nodes_always_false([Node|Nodes], Which, DI) -->
  627        { which_node_child(Which, Node, Child),
  628          opposite(Which, Opposite) },
  629        opposite_always_false(Opposite, DI, Child),
  630        nodes_always_false(Nodes, Which, DI).
  631
  632which_node_child(low, Node, Child) :-
  633        node_var_low_high(Node, _, Child, _).
  634which_node_child(high, Node, Child) :-
  635        node_var_low_high(Node, _, _, Child).
  636
  637opposite(low, high).
  638opposite(high, low).
  639
  640opposite_always_false(Opposite, DI, Node) -->
  641        (   { node_visited(Node) } -> []
  642        ;   { node_var_low_high(Node, Var, Low, High),
  643              with_aux(put_visited, Node),
  644              var_index(Var, VI) },
  645            [Node],
  646            (   { VI =:= DI } ->
  647                { which_node_child(Opposite, Node, Child),
  648                  Child == 0 }
  649            ;   opposite_always_false(Opposite, DI, Low),
  650                opposite_always_false(Opposite, DI, High)
  651            )
  652        ).
  653
  654further_branching(Node) :-
  655        node_var_low_high(Node, _, Low, High),
  656        Low \== 1,
  657        High \== 1.
  658
  659negative_decisive(Node) :-
  660        node_var_low_high(Node, _, Low, High),
  661        (   Low == 0 -> true
  662        ;   High == 0 -> true
  663        ;   false
  664        ).
  665
  666/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  667   Instantiate all variables that only admit a single Boolean value.
  668
  669   This is the case if: The variable is not skipped in any branch
  670   leading to 1 (its being skipped means that it may be assigned
  671   either 0 or 1 and can thus not be fixed yet), and all nodes where
  672   it occurs as a branching variable have either lower or upper child
  673   fixed to 0 consistently.
  674- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  675
  676domain_consistency(Eqs, Goal) :-
  677        maplist(eq_a_b, Eqs, Vs, Values),
  678        Goal = (Vs = Values). % propagate all assignments at once
  679
  680eq_a_b(A=B, A, B).
  681
  682consistently_false_(Which, Node) :-
  683        which_node_child(Which, Node, Child),
  684        Child == 0.
  685
  686/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  687   In essentially one sweep of the BDD, all variables can be classified:
  688   Unification with 0 or 1, further branching and/or negative decisive.
  689
  690   Strategy: Breadth-first traversal of the BDD, failing (and thus
  691   clearing all attributes) if the variable is skipped in some branch,
  692   and moving the frontier along each time.
  693
  694   A formula is only satisfiable if it is a tautology after all (also
  695   implicitly) existentially quantified variables are projected away.
  696   However, we only need to check this explicitly if at least one
  697   universally quantified variable appears. Otherwise, we know that
  698   the formula is satisfiable at this point, because its BDD is not 0.
  699- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  700
  701bdd_variables_classification(BDD, Nodes, Classes) :-
  702        nodes_variables(Nodes, Vs0),
  703        variables_in_index_order(Vs0, Vs),
  704        (   partition(universal_var, Vs, [_|_], Es) ->
  705            foldl(existential, Es, BDD, 1)
  706        ;   true
  707        ),
  708        phrase(variables_classification(Vs, [BDD]), Classes),
  709        maplist(with_aux(unvisit), Nodes).
  710
  711variables_classification([], _) --> [].
  712variables_classification([V|Vs], Nodes0) -->
  713        { var_index(V, Index) },
  714        (   { phrase(nodes_with_variable(Nodes0, Index), Nodes) } ->
  715            (   { maplist(consistently_false_(low), Nodes) } -> [V=1]
  716            ;   { maplist(consistently_false_(high), Nodes) } -> [V=0]
  717            ;   []
  718            ),
  719            (   { maplist(further_branching, Nodes) } ->
  720                [further_branching(V, Nodes)]
  721            ;   []
  722            ),
  723            (   { maplist(negative_decisive, Nodes) } ->
  724                [negative_decisive(V)]
  725            ;   []
  726            ),
  727            { maplist(with_aux(unvisit), Nodes) },
  728            variables_classification(Vs, Nodes)
  729        ;   variables_classification(Vs, Nodes0)
  730        ).
  731
  732nodes_with_variable([], _) --> [].
  733nodes_with_variable([Node|Nodes], VI) -->
  734        { Node \== 1 },
  735        (   { node_visited(Node) } -> nodes_with_variable(Nodes, VI)
  736        ;   { with_aux(put_visited, Node),
  737              node_var_low_high(Node, OVar, Low, High),
  738              var_index(OVar, OVI) },
  739            { OVI =< VI },
  740            (   { OVI =:= VI } -> [Node]
  741            ;   nodes_with_variable([Low,High], VI)
  742            ),
  743            nodes_with_variable(Nodes, VI)
  744        ).
  745
  746/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  747   Node management. Always use an existing node, if there is one.
  748- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  749
  750make_node(Var, Low, High, Node) :-
  751        (   Low == High -> Node = Low
  752        ;   low_high_key(Low, High, Key),
  753            (   lookup_node(Var, Key, Node) -> true
  754            ;   clpb_next_id('$clpb_next_node', ID),
  755                Node = node(ID,Var,Low,High,_Aux),
  756                register_node(Var, Key, Node)
  757            )
  758        ).
  759
  760make_node(Var, Low, High, Node) -->
  761        % make it conveniently usable within DCGs
  762        { make_node(Var, Low, High, Node) }.
  763
  764
  765% The key of a node for hashing is determined by the IDs of its
  766% children.
  767
  768low_high_key(Low, High, node(LID,HID)) :-
  769        node_id(Low, LID),
  770        node_id(High, HID).
  771
  772
  773rebuild_hashes(BDD) :-
  774        bdd_nodes(nodevar_put_empty_hash, BDD, Nodes),
  775        maplist(re_register_node, Nodes).
  776
  777nodevar_put_empty_hash(Node) :-
  778        node_var_low_high(Node, Var, _, _),
  779        empty_assoc(H0),
  780        put_attr(Var, clpb_hash, H0).
  781
  782re_register_node(Node) :-
  783        node_var_low_high(Node, Var, Low, High),
  784        low_high_key(Low, High, Key),
  785        register_node(Var, Key, Node).
  786
  787register_node(Var, Key, Node) :-
  788        get_attr(Var, clpb_hash, H0),
  789        put_assoc(Key, H0, Node, H),
  790        put_attr(Var, clpb_hash, H).
  791
  792lookup_node(Var, Key, Node) :-
  793        get_attr(Var, clpb_hash, H0),
  794        get_assoc(Key, H0, Node).
  795
  796
  797node_id(0, false).
  798node_id(1, true).
  799node_id(node(ID,_,_,_,_), ID).
  800
  801node_aux(Node, Aux) :- arg(5, Node, Aux).
  802
  803node_var_low_high(Node, Var, Low, High) :-
  804        Node = node(_,Var,Low,High,_).
  805
  806
  807/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  808   sat_bdd/2 converts a SAT formula in canonical form to an ordered
  809   and reduced BDD.
  810- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  811
  812sat_bdd(V, Node)           :- var(V), !, make_node(V, 0, 1, Node).
  813sat_bdd(I, I)              :- integer(I), !.
  814sat_bdd(V^Sat, Node)       :- !, sat_bdd(Sat, BDD), existential(V, BDD, Node).
  815sat_bdd(card(Is,Fs), Node) :- !, counter_network(Is, Fs, Node).
  816sat_bdd(Sat, Node)         :- !,
  817        Sat =.. [F,A,B],
  818        sat_bdd(A, NA),
  819        sat_bdd(B, NB),
  820        apply(F, NA, NB, Node).
  821
  822existential(V, BDD, Node) :-
  823        var_index(V, Index),
  824        bdd_restriction(BDD, Index, 0, NA),
  825        bdd_restriction(BDD, Index, 1, NB),
  826        apply(+, NA, NB, Node).
  827
  828/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  829   Counter network for card(Is,Fs).
  830- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  831
  832counter_network(Cs, Fs, Node) :-
  833        same_length([_|Fs], Indicators),
  834        fill_indicators(Indicators, 0, Cs),
  835        phrase(formulas_variables(Fs, Vars0), ExBDDs),
  836        maplist(unvisit, Vars0),
  837        % The counter network is built bottom-up, so variables with
  838        % highest index must be processed first.
  839        variables_in_index_order(Vars0, Vars1),
  840        reverse(Vars1, Vars),
  841        counter_network_(Vars, Indicators, Node0),
  842        foldl(existential_and, ExBDDs, Node0, Node).
  843
  844% Introduce fresh variables for expressions that are not variables.
  845% These variables are later existentially quantified to remove them.
  846% Also, new variables are introduced for variables that are used more
  847% than once, as in card([0,1],[X,X,Y]), to keep the BDD ordered.
  848
  849formulas_variables([], []) --> [].
  850formulas_variables([F|Fs], [V|Vs]) -->
  851        (   { var(F), \+ is_visited(F) } ->
  852            { V = F,
  853              put_visited(F) }
  854        ;   { enumerate_variable(V),
  855              sat_rewrite(V =:= F, Sat),
  856              sat_bdd(Sat, BDD) },
  857            [V-BDD]
  858        ),
  859        formulas_variables(Fs, Vs).
  860
  861counter_network_([], [Node], Node).
  862counter_network_([Var|Vars], [I|Is0], Node) :-
  863        foldl(indicators_pairing(Var), Is0, Is, I, _),
  864        counter_network_(Vars, Is, Node).
  865
  866indicators_pairing(Var, I, Node, Prev, I) :- make_node(Var, Prev, I, Node).
  867
  868fill_indicators([], _, _).
  869fill_indicators([I|Is], Index0, Cs) :-
  870        (   memberchk(Index0, Cs) -> I = 1
  871        ;   member(A-B, Cs), between(A, B, Index0) -> I = 1
  872        ;   I = 0
  873        ),
  874        Index1 is Index0 + 1,
  875        fill_indicators(Is, Index1, Cs).
  876
  877existential_and(Ex-BDD, Node0, Node) :-
  878        bdd_and(BDD, Node0, Node1),
  879        existential(Ex, Node1, Node),
  880        % remove attributes to avoid residual goals for variables that
  881        % are only used as substitutes for formulas
  882        del_attrs(Ex).
  883
  884/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  885   Compute F(NA, NB).
  886
  887   We use a DCG to thread through an implicit argument G0, an
  888   association table F(IDA,IDB) -> Node, used for memoization.
  889- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  890
  891apply(F, NA, NB, Node) :-
  892        empty_assoc(G0),
  893        phrase(apply(F, NA, NB, Node), [G0], _).
  894
  895apply(F, NA, NB, Node) -->
  896        (   { integer(NA), integer(NB) } -> { once(bool_op(F, NA, NB, Node)) }
  897        ;   { apply_shortcut(F, NA, NB, Node) } -> []
  898        ;   { node_id(NA, IDA), node_id(NB, IDB), Key =.. [F,IDA,IDB] },
  899            (   state(G0), { get_assoc(Key, G0, Node) } -> []
  900            ;   apply_(F, NA, NB, Node),
  901                state(G0, G),
  902                { put_assoc(Key, G0, Node, G) }
  903            )
  904        ).
  905
  906apply_shortcut(+, NA, NB, Node) :-
  907        (   NA == 0 -> Node = NB
  908        ;   NA == 1 -> Node = 1
  909        ;   NB == 0 -> Node = NA
  910        ;   NB == 1 -> Node = 1
  911        ;   false
  912        ).
  913
  914apply_shortcut(*, NA, NB, Node) :-
  915        (   NA == 0 -> Node = 0
  916        ;   NA == 1 -> Node = NB
  917        ;   NB == 0 -> Node = 0
  918        ;   NB == 1 -> Node = NA
  919        ;   false
  920        ).
  921
  922
  923apply_(F, NA, NB, Node) -->
  924        { var_less_than(NA, NB),
  925          !,
  926          node_var_low_high(NA, VA, LA, HA) },
  927        apply(F, LA, NB, Low),
  928        apply(F, HA, NB, High),
  929        make_node(VA, Low, High, Node).
  930apply_(F, NA, NB, Node) -->
  931        { node_var_low_high(NA, VA, LA, HA),
  932          node_var_low_high(NB, VB, LB, HB),
  933          VA == VB },
  934        !,
  935        apply(F, LA, LB, Low),
  936        apply(F, HA, HB, High),
  937        make_node(VA, Low, High, Node).
  938apply_(F, NA, NB, Node) --> % NB < NA
  939        { node_var_low_high(NB, VB, LB, HB) },
  940        apply(F, NA, LB, Low),
  941        apply(F, NA, HB, High),
  942        make_node(VB, Low, High, Node).
  943
  944
  945node_varindex(Node, VI) :-
  946        node_var_low_high(Node, V, _, _),
  947        var_index(V, VI).
  948
  949var_less_than(NA, NB) :-
  950        (   integer(NB) -> true
  951        ;   node_varindex(NA, VAI),
  952            node_varindex(NB, VBI),
  953            VAI < VBI
  954        ).
  955
  956bool_op(+, 0, 0, 0).
  957bool_op(+, 0, 1, 1).
  958bool_op(+, 1, 0, 1).
  959bool_op(+, 1, 1, 1).
  960
  961bool_op(*, 0, 0, 0).
  962bool_op(*, 0, 1, 0).
  963bool_op(*, 1, 0, 0).
  964bool_op(*, 1, 1, 1).
  965
  966bool_op(#, 0, 0, 0).
  967bool_op(#, 0, 1, 1).
  968bool_op(#, 1, 0, 1).
  969bool_op(#, 1, 1, 0).
  970
  971/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  972   Access implicit state in DCGs.
  973- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  974
  975state(S) --> state(S, S).
  976
  977state(S0, S), [S] --> [S0].
  978
  979/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  980   Unification. X = Expr is equivalent to sat(X =:= Expr).
  981
  982   Current limitation:
  983   ===================
  984
  985   The current interface of attributed variables is not general enough
  986   to express what we need. For example,
  987
  988       ?- sat(A + B), A = A + 1.
  989
  990   should be equivalent to
  991
  992       ?- sat(A + B), sat(A =:= A + 1).
  993
  994   However, attr_unify_hook/2 is only called *after* the unification
  995   of A with A + 1 has already taken place and turned A into a cyclic
  996   ground term, raised an error or failed (depending on the flag
  997   occurs_check), making it impossible to reason about the variable A
  998   in the unification hook. Therefore, a more general interface for
  999   attributed variables should replace the current one. In particular,
 1000   unification filters should be able to reason about terms before
 1001   they are unified with anything.
 1002- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1003
 1004attr_unify_hook(index_root(I,Root), Other) :-
 1005        (   integer(Other) ->
 1006            (   between(0, 1, Other) ->
 1007                root_get_formula_bdd(Root, Sat, BDD0),
 1008                bdd_restriction(BDD0, I, Other, BDD),
 1009                root_put_formula_bdd(Root, Sat, BDD),
 1010                satisfiable_bdd(BDD)
 1011            ;   no_truth_value(Other)
 1012            )
 1013        ;   atom(Other) ->
 1014            root_get_formula_bdd(Root, Sat0, _),
 1015            parse_sat(Sat0, Sat),
 1016            sat_bdd(Sat, BDD),
 1017            root_put_formula_bdd(Root, Sat0, BDD),
 1018            is_bdd(BDD),
 1019            satisfiable_bdd(BDD)
 1020        ;   % due to variable aliasing, any BDDs may now be unordered,
 1021            % so we need to rebuild the new BDD from the conjunction.
 1022            root_get_formula_bdd(Root, Sat0, _),
 1023            Sat = Sat0*OtherSat,
 1024            (   var(Other), var_index_root(Other, _, OtherRoot),
 1025                OtherRoot \== Root ->
 1026                root_get_formula_bdd(OtherRoot, OtherSat, _),
 1027                parse_sat(Sat, Sat1),
 1028                sat_bdd(Sat1, BDD1),
 1029                And = Sat,
 1030                sat_roots(Sat, Roots)
 1031            ;   parse_sat(Other, OtherSat),
 1032                sat_roots(Sat, Roots),
 1033                maplist(root_rebuild_bdd, Roots),
 1034                roots_and(Roots, 1-1, And-BDD1)
 1035            ),
 1036            maplist(del_bdd, Roots),
 1037            maplist(=(NewRoot), Roots),
 1038            root_put_formula_bdd(NewRoot, And, BDD1),
 1039            is_bdd(BDD1),
 1040            satisfiable_bdd(BDD1)
 1041        ).
 1042
 1043root_rebuild_bdd(Root) :-
 1044        (   root_get_formula_bdd(Root, F0, _) ->
 1045            parse_sat(F0, Sat),
 1046            sat_bdd(Sat, BDD),
 1047            root_put_formula_bdd(Root, F0, BDD)
 1048        ;   true
 1049        ).
 1050
 1051/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1052   Support for project_attributes/2.
 1053
 1054   This is called by the toplevel as
 1055
 1056      project_attributes(+QueryVars, +AttrVars)
 1057
 1058   in order to project all remaining constraints onto QueryVars.
 1059
 1060   All CLP(B) variables that do not occur in QueryVars or AttrVars
 1061   need to be existentially quantified, so that they do not occur in
 1062   residual goals. This is very easy to do in the case of CLP(B).
 1063- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1064
 1065project_attributes(QueryVars0, AttrVars) :-
 1066        append(QueryVars0, AttrVars, QueryVars1),
 1067        include(clpb_variable, QueryVars1, QueryVars),
 1068        maplist(var_index_root, QueryVars, _, Roots0),
 1069        sort(Roots0, Roots),
 1070        maplist(remove_hidden_variables(QueryVars), Roots).
 1071
 1072clpb_variable(Var) :- var_index(Var, _).
 1073
 1074/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1075   All CLP(B) variables occurring in BDDs but not in query variables
 1076   become existentially quantified. This must also be reflected in the
 1077   formula. In addition, an attribute is attached to these variables
 1078   to suppress superfluous sat(V=:=V) goals.
 1079- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1080
 1081remove_hidden_variables(QueryVars, Root) :-
 1082        root_get_formula_bdd(Root, Formula, BDD0),
 1083        maplist(put_visited, QueryVars),
 1084        bdd_variables(BDD0, HiddenVars0),
 1085        exclude(universal_var, HiddenVars0, HiddenVars),
 1086        maplist(unvisit, QueryVars),
 1087        foldl(existential, HiddenVars, BDD0, BDD),
 1088        foldl(quantify_existantially, HiddenVars, Formula, ExFormula),
 1089        root_put_formula_bdd(Root, ExFormula, BDD).
 1090
 1091quantify_existantially(E, E0, E^E0) :- put_attr(E, clpb_omit_boolean, true).
 1092
 1093/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1094   BDD restriction.
 1095- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1096
 1097bdd_restriction(Node, VI, Value, Res) :-
 1098        empty_assoc(G0),
 1099        phrase(bdd_restriction_(Node, VI, Value, Res), [G0], _),
 1100        is_bdd(Res).
 1101
 1102bdd_restriction_(Node, VI, Value, Res) -->
 1103        (   { integer(Node) } -> { Res = Node }
 1104        ;   { node_var_low_high(Node, Var, Low, High) } ->
 1105            (   { integer(Var) } ->
 1106                (   { Var =:= 0 } -> bdd_restriction_(Low, VI, Value, Res)
 1107                ;   { Var =:= 1 } -> bdd_restriction_(High, VI, Value, Res)
 1108                ;   { no_truth_value(Var) }
 1109                )
 1110            ;   { var_index(Var, I0),
 1111                  node_id(Node, ID) },
 1112                (   { I0 =:= VI } ->
 1113                    (   { Value =:= 0 } -> { Res = Low }
 1114                    ;   { Value =:= 1 } -> { Res = High }
 1115                    )
 1116                ;   { I0 > VI } -> { Res = Node }
 1117                ;   state(G0), { get_assoc(ID, G0, Res) } -> []
 1118                ;   bdd_restriction_(Low, VI, Value, LRes),
 1119                    bdd_restriction_(High, VI, Value, HRes),
 1120                    make_node(Var, LRes, HRes, Res),
 1121                    state(G0, G),
 1122                    { put_assoc(ID, G0, Res, G) }
 1123                )
 1124            )
 1125        ;   { domain_error(node, Node) }
 1126        ).
 1127
 1128/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1129   Relating a BDD to its elements (nodes and variables).
 1130
 1131   Note that BDDs can become quite big (easily millions of nodes), and
 1132   memory space is a major bottleneck for many problems. If possible,
 1133   we therefore do not duplicate the entire BDD in memory (as in
 1134   bdd_ites/2), but only extract its features as needed.
 1135- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1136
 1137bdd_nodes(BDD, Ns) :- bdd_nodes(ignore_node, BDD, Ns).
 1138
 1139ignore_node(_).
 1140
 1141% VPred is a unary predicate that is called for each node that has a
 1142% branching variable (= each inner node).
 1143
 1144bdd_nodes(VPred, BDD, Ns) :-
 1145        phrase(bdd_nodes_(VPred, BDD), Ns),
 1146        maplist(with_aux(unvisit), Ns).
 1147
 1148bdd_nodes_(VPred, Node) -->
 1149        (   { node_visited(Node) } -> []
 1150        ;   { call(VPred, Node),
 1151              with_aux(put_visited, Node),
 1152              node_var_low_high(Node, _, Low, High) },
 1153            [Node],
 1154            bdd_nodes_(VPred, Low),
 1155            bdd_nodes_(VPred, High)
 1156        ).
 1157
 1158node_visited(Node) :- integer(Node).
 1159node_visited(Node) :- with_aux(is_visited, Node).
 1160
 1161bdd_variables(BDD, Vs) :-
 1162        bdd_nodes(BDD, Nodes),
 1163        nodes_variables(Nodes, Vs).
 1164
 1165nodes_variables(Nodes, Vs) :-
 1166        phrase(nodes_variables_(Nodes), Vs),
 1167        maplist(unvisit, Vs).
 1168
 1169nodes_variables_([]) --> [].
 1170nodes_variables_([Node|Nodes]) -->
 1171        { node_var_low_high(Node, Var, _, _) },
 1172        (   { integer(Var) } -> []
 1173        ;   { is_visited(Var) } -> []
 1174        ;   { put_visited(Var) },
 1175            [Var]
 1176        ),
 1177        nodes_variables_(Nodes).
 1178
 1179unvisit(V) :- del_attr(V, clpb_visited).
 1180
 1181is_visited(V) :- get_attr(V, clpb_visited, true).
 1182
 1183put_visited(V) :- put_attr(V, clpb_visited, true).
 1184
 1185with_aux(Pred, Node) :-
 1186        node_aux(Node, Aux),
 1187        call(Pred, Aux).
 1188
 1189/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1190   Internal consistency checks.
 1191
 1192   To enable these checks, set the flag clpb_validation to true.
 1193- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1194
 1195is_bdd(BDD) :-
 1196        (   current_prolog_flag(clpb_validation, true) ->
 1197            bdd_ites(BDD, ITEs),
 1198            pairs_values(ITEs, Ls0),
 1199            sort(Ls0, Ls1),
 1200            (   same_length(Ls0, Ls1) -> true
 1201            ;   domain_error(reduced_ites, (ITEs,Ls0,Ls1))
 1202            ),
 1203            (   member(ITE, ITEs), \+ registered_node(ITE) ->
 1204                domain_error(registered_node, ITE)
 1205            ;   true
 1206            ),
 1207            (   member(I, ITEs), \+ ordered(I) ->
 1208                domain_error(ordered_node, I)
 1209            ;   true
 1210            )
 1211        ;   true
 1212        ).
 1213
 1214ordered(_-ite(Var,High,Low)) :-
 1215        (   var_index(Var, VI) ->
 1216            greater_varindex_than(High, VI),
 1217            greater_varindex_than(Low, VI)
 1218        ;   true
 1219        ).
 1220
 1221greater_varindex_than(Node, VI) :-
 1222        (   integer(Node) -> true
 1223        ;   node_var_low_high(Node, Var, _, _),
 1224            (   var_index(Var, OI) ->
 1225                OI > VI
 1226            ;   true
 1227            )
 1228        ).
 1229
 1230registered_node(Node-ite(Var,High,Low)) :-
 1231        (   var(Var) ->
 1232            low_high_key(Low, High, Key),
 1233            lookup_node(Var, Key, Node0),
 1234            Node == Node0
 1235        ;   true
 1236        ).
 1237
 1238bdd_ites(BDD, ITEs) :-
 1239        bdd_nodes(BDD, Nodes),
 1240        maplist(node_ite, Nodes, ITEs).
 1241
 1242node_ite(Node, Node-ite(Var,High,Low)) :-
 1243        node_var_low_high(Node, Var, Low, High).
 1244
 1245%% labeling(+Vs) is multi.
 1246%
 1247% Enumerate concrete solutions. Assigns truth values to the Boolean
 1248% variables Vs such that all stated constraints are satisfied.
 1249
 1250labeling(Vs0) :-
 1251        must_be(list, Vs0),
 1252        maplist(labeling_var, Vs0),
 1253        variables_in_index_order(Vs0, Vs),
 1254        maplist(indomain, Vs).
 1255
 1256labeling_var(V) :- var(V), !.
 1257labeling_var(V) :- V == 0, !.
 1258labeling_var(V) :- V == 1, !.
 1259labeling_var(V) :- domain_error(clpb_variable, V).
 1260
 1261variables_in_index_order(Vs0, Vs) :-
 1262        maplist(var_with_index, Vs0, IVs0),
 1263        keysort(IVs0, IVs),
 1264        pairs_values(IVs, Vs).
 1265
 1266var_with_index(V, I-V) :-
 1267        (   var_index_root(V, I, _) -> true
 1268        ;   I = 0
 1269        ).
 1270
 1271indomain(0).
 1272indomain(1).
 1273
 1274
 1275%% sat_count(+Expr, -Count) is det.
 1276%
 1277% Count the number of admissible assignments. Count is the number of
 1278% different assignments of truth values to the variables in the
 1279% Boolean expression Expr, such that Expr is true and all posted
 1280% constraints are satisfiable.
 1281%
 1282% A common form of invocation is `sat_count(+[1|Vs], Count)`: This
 1283% counts the number of admissible assignments to `Vs` without imposing
 1284% any further constraints.
 1285%
 1286% Examples:
 1287%
 1288% ==
 1289% ?- sat(A =< B), Vs = [A,B], sat_count(+[1|Vs], Count).
 1290% Vs = [A, B],
 1291% Count = 3,
 1292% sat(A=:=A*B).
 1293%
 1294% ?- length(Vs, 120),
 1295%    sat_count(+Vs, CountOr),
 1296%    sat_count(*(Vs), CountAnd).
 1297% Vs = [...],
 1298% CountOr = 1329227995784915872903807060280344575,
 1299% CountAnd = 1.
 1300% ==
 1301
 1302sat_count(Sat0, N) :-
 1303        catch((parse_sat(Sat0, Sat),
 1304               sat_bdd(Sat, BDD),
 1305               sat_roots(Sat, Roots),
 1306               roots_and(Roots, _-BDD, _-BDD1),
 1307               % we mark variables that occur in Sat0 as visited ...
 1308               term_variables(Sat0, Vs),
 1309               maplist(put_visited, Vs),
 1310               % ... so that they do not appear in Vs1 ...
 1311               bdd_variables(BDD1, Vs1),
 1312               partition(universal_var, Vs1, Univs, Exis),
 1313               % ... and then remove remaining variables:
 1314               foldl(universal, Univs, BDD1, BDD2),
 1315               foldl(existential, Exis, BDD2, BDD3),
 1316               variables_in_index_order(Vs, IVs),
 1317               foldl(renumber_variable, IVs, 1, VNum),
 1318               bdd_count(BDD3, VNum, Count0),
 1319               var_u(BDD3, VNum, P),
 1320               % Do not unify N directly, because we are not prepared
 1321               % for propagation here in case N is a CLP(B) variable.
 1322               N0 is 2^(P - 1)*Count0,
 1323               % reset all attributes and Aux variables
 1324               throw(count(N0))),
 1325              count(N0),
 1326              N = N0).
 1327
 1328universal(V, BDD, Node) :-
 1329        var_index(V, Index),
 1330        bdd_restriction(BDD, Index, 0, NA),
 1331        bdd_restriction(BDD, Index, 1, NB),
 1332        apply(*, NA, NB, Node).
 1333
 1334renumber_variable(V, I0, I) :-
 1335        put_attr(V, clpb, index_root(I0,_)),
 1336        I is I0 + 1.
 1337
 1338bdd_count(Node, VNum, Count) :-
 1339        (   integer(Node) -> Count = Node
 1340        ;   node_aux(Node, Count),
 1341            (   integer(Count) -> true
 1342            ;   node_var_low_high(Node, V, Low, High),
 1343                bdd_count(Low, VNum, LCount),
 1344                bdd_count(High, VNum, HCount),
 1345                bdd_pow(Low, V, VNum, LPow),
 1346                bdd_pow(High, V, VNum, HPow),
 1347                Count is LPow*LCount + HPow*HCount
 1348            )
 1349        ).
 1350
 1351
 1352bdd_pow(Node, V, VNum, Pow) :-
 1353        var_index(V, Index),
 1354        var_u(Node, VNum, P),
 1355        Pow is 2^(P - Index - 1).
 1356
 1357var_u(Node, VNum, Index) :-
 1358        (   integer(Node) -> Index = VNum
 1359        ;   node_varindex(Node, Index)
 1360        ).
 1361
 1362/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1363   Pick a solution in such a way that each solution is equally likely.
 1364- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1365
 1366%% random_labeling(+Seed, +Vs) is det.
 1367%
 1368% Select a single random solution. An admissible assignment of truth
 1369% values to the Boolean variables in Vs is chosen in such a way that
 1370% each admissible assignment is equally likely. Seed is an integer,
 1371% used as the initial seed for the random number generator.
 1372
 1373single_bdd(Vars0) :-
 1374        maplist(monotonic_variable, Vars0, Vars),
 1375        % capture all variables with a single BDD
 1376        sat(+[1|Vars]).
 1377
 1378random_labeling(Seed, Vars) :-
 1379        must_be(list, Vars),
 1380        set_random(seed(Seed)),
 1381        (   ground(Vars) -> true
 1382        ;   catch((single_bdd(Vars),
 1383                   once((member(Var, Vars),var(Var))),
 1384                   var_index_root(Var, _, Root),
 1385                   root_get_formula_bdd(Root, _, BDD),
 1386                   bdd_variables(BDD, Vs),
 1387                   variables_in_index_order(Vs, IVs),
 1388                   foldl(renumber_variable, IVs, 1, VNum),
 1389                   phrase(random_bindings(VNum, BDD), Bs),
 1390                   maplist(del_attrs, Vs),
 1391                   % reset all attribute modifications
 1392                   throw(randsol(Vars, Bs))),
 1393                  randsol(Vars, Bs),
 1394                  true),
 1395            maplist(call, Bs),
 1396            % set remaining variables to 0 or 1 with equal probability
 1397            include(var, Vars, Remaining),
 1398            maplist(maybe_zero, Remaining)
 1399        ).
 1400
 1401maybe_zero(Var) :-
 1402        (   maybe -> Var = 0
 1403        ;   Var = 1
 1404        ).
 1405
 1406random_bindings(_, Node) --> { Node == 1 }, !.
 1407random_bindings(VNum, Node) -->
 1408        { node_var_low_high(Node, Var, Low, High),
 1409          bdd_count(Node, VNum, Total),
 1410          bdd_count(Low, VNum, LCount) },
 1411        (   { maybe(LCount, Total) } ->
 1412            [Var=0], random_bindings(VNum, Low)
 1413        ;   [Var=1], random_bindings(VNum, High)
 1414        ).
 1415
 1416/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1417   Find solutions with maximum weight.
 1418- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1419
 1420%% weighted_maximum(+Weights, +Vs, -Maximum) is multi.
 1421%
 1422% Enumerate weighted optima over admissible assignments. Maximize a
 1423% linear objective function over Boolean variables Vs with integer
 1424% coefficients Weights. This predicate assigns 0 and 1 to the
 1425% variables in Vs such that all stated constraints are satisfied, and
 1426% Maximum is the maximum of sum(Weight_i*V_i) over all admissible
 1427% assignments.  On backtracking, all admissible assignments that
 1428% attain the optimum are generated.
 1429%
 1430% This predicate can also be used to _minimize_ a linear Boolean
 1431% program, since negative integers can appear in Weights.
 1432%
 1433% Example:
 1434%
 1435% ==
 1436% ?- sat(A#B), weighted_maximum([1,2,1], [A,B,C], Maximum).
 1437% A = 0, B = 1, C = 1, Maximum = 3.
 1438% ==
 1439
 1440weighted_maximum(Ws, Vars, Max) :-
 1441        must_be(list(integer), Ws),
 1442        must_be(list(var), Vars),
 1443        single_bdd(Vars),
 1444        Vars = [Var|_],
 1445        var_index_root(Var, _,  Root),
 1446        root_get_formula_bdd(Root, _, BDD0),
 1447        bdd_variables(BDD0, Vs),
 1448        % existentially quantify variables that are not considered
 1449        maplist(put_visited, Vars),
 1450        exclude(is_visited, Vs, Unvisited),
 1451        maplist(unvisit, Vars),
 1452        foldl(existential, Unvisited, BDD0, BDD),
 1453        maplist(var_with_index, Vars, IVs),
 1454        pairs_keys_values(Pairs0, IVs, Ws),
 1455        keysort(Pairs0, Pairs1),
 1456        pairs_keys_values(Pairs1, IVs1, WeightsIndexOrder),
 1457        pairs_values(IVs1, VarsIndexOrder),
 1458        % Pairs is a list of Var-Weight terms, in index order of Vars
 1459        pairs_keys_values(Pairs, VarsIndexOrder, WeightsIndexOrder),
 1460        bdd_maximum(BDD, Pairs, Max),
 1461        max_labeling(BDD, Pairs).
 1462
 1463max_labeling(1, Pairs) :- max_upto(Pairs, _, _).
 1464max_labeling(node(_,Var,Low,High,Aux), Pairs0) :-
 1465        max_upto(Pairs0, Var, Pairs),
 1466        get_attr(Aux, clpb_max, max(_,Dir)),
 1467        direction_labeling(Dir, Var, Low, High, Pairs).
 1468
 1469max_upto([], _, _).
 1470max_upto([Var0-Weight|VWs0], Var, VWs) :-
 1471        (   Var == Var0 -> VWs = VWs0
 1472        ;   Weight =:= 0 ->
 1473            (   Var0 = 0 ; Var0 = 1 ),
 1474            max_upto(VWs0, Var, VWs)
 1475        ;   Weight < 0 -> Var0 = 0, max_upto(VWs0, Var, VWs)
 1476        ;   Var0 = 1, max_upto(VWs0, Var, VWs)
 1477        ).
 1478
 1479direction_labeling(low, 0, Low, _, Pairs)   :- max_labeling(Low, Pairs).
 1480direction_labeling(high, 1, _, High, Pairs) :- max_labeling(High, Pairs).
 1481
 1482bdd_maximum(1, Pairs, Max) :-
 1483        pairs_values(Pairs, Weights0),
 1484        include(<(0), Weights0, Weights),
 1485        sum_list(Weights, Max).
 1486bdd_maximum(node(_,Var,Low,High,Aux), Pairs0, Max) :-
 1487        (   get_attr(Aux, clpb_max, max(Max,_)) -> true
 1488        ;   (   skip_to_var(Var, Weight, Pairs0, Pairs),
 1489                (   Low == 0 ->
 1490                    bdd_maximum_(High, Pairs, MaxHigh, MaxToHigh),
 1491                    Max is MaxToHigh + MaxHigh + Weight,
 1492                    Dir = high
 1493                ;   High == 0 ->
 1494                    bdd_maximum_(Low, Pairs, MaxLow, MaxToLow),
 1495                    Max is MaxToLow + MaxLow,
 1496                    Dir = low
 1497                ;   bdd_maximum_(Low, Pairs, MaxLow, MaxToLow),
 1498                    bdd_maximum_(High, Pairs, MaxHigh, MaxToHigh),
 1499                    Max0 is MaxToLow + MaxLow,
 1500                    Max1 is MaxToHigh + MaxHigh + Weight,
 1501                    Max is max(Max0,Max1),
 1502                    (   Max0 =:= Max1 -> Dir = _Any
 1503                    ;   Max0 < Max1 -> Dir = high
 1504                    ;   Dir = low
 1505                    )
 1506                ),
 1507                store_maximum(Aux, Max, Dir)
 1508            )
 1509        ).
 1510
 1511bdd_maximum_(Node, Pairs, Max, MaxTo) :-
 1512	bdd_maximum(Node, Pairs, Max),
 1513	between_weights(Node, Pairs, MaxTo).
 1514
 1515store_maximum(Aux, Max, Dir) :- put_attr(Aux, clpb_max, max(Max,Dir)).
 1516
 1517between_weights(Node, Pairs0, MaxTo) :-
 1518        (   Node == 1 -> MaxTo = 0
 1519        ;   node_var_low_high(Node, Var, _, _),
 1520            phrase(skip_to_var_(Var, _, Pairs0, _), Weights0),
 1521            include(<(0), Weights0, Weights),
 1522            sum_list(Weights, MaxTo)
 1523        ).
 1524
 1525skip_to_var(Var, Weight, Pairs0, Pairs) :-
 1526        phrase(skip_to_var_(Var, Weight, Pairs0, Pairs), _).
 1527
 1528skip_to_var_(Var, Weight, [Var0-Weight0|VWs0], VWs) -->
 1529        (   { Var == Var0 } ->
 1530            { Weight = Weight0, VWs0 = VWs }
 1531        ;   (   { Weight0 =< 0 } -> []
 1532            ;   [Weight0]
 1533            ),
 1534            skip_to_var_(Var, Weight, VWs0, VWs)
 1535        ).
 1536
 1537/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1538   Projection to residual goals.
 1539- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1540
 1541attribute_goals(Var) -->
 1542        { var_index_root(Var, _, Root) },
 1543        (   { root_get_formula_bdd(Root, Formula, BDD) } ->
 1544            { del_bdd(Root) },
 1545            (   { current_prolog_flag(clpb_residuals, bdd) } ->
 1546                { bdd_nodes(BDD, Nodes),
 1547                  phrase(nodes(Nodes), Ns) },
 1548                [clpb:'$clpb_bdd'(Ns)]
 1549            ;   { prepare_global_variables(BDD),
 1550                  phrase(sat_ands(Formula), Ands0),
 1551                  ands_fusion(Ands0, Ands),
 1552                  maplist(formula_anf, Ands, ANFs0),
 1553                  sort(ANFs0, ANFs1),
 1554                  exclude(eq_1, ANFs1, ANFs2),
 1555                  variables_separation(ANFs2, ANFs) },
 1556                sats(ANFs)
 1557            ),
 1558            (   { get_attr(Var, clpb_atom, Atom) } ->
 1559                [clpb:sat(Var=:=Atom)]
 1560            ;   []
 1561            ),
 1562            % formula variables not occurring in the BDD should be booleans
 1563            { bdd_variables(BDD, Vs),
 1564              maplist(del_clpb, Vs),
 1565              term_variables(Formula, RestVs0),
 1566              include(clpb_variable, RestVs0, RestVs) },
 1567            booleans(RestVs)
 1568        ;   boolean(Var)  % the variable may have occurred only in taut/2
 1569        ).
 1570
 1571del_clpb(Var) :- del_attr(Var, clpb).
 1572
 1573/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1574   To make residual projection work with recorded constraints, the
 1575   global counters must be adjusted so that new variables and nodes
 1576   also get new IDs. Also, clpb_next_id/2 is used to actually create
 1577   these counters, because creating them with b_setval/2 would make
 1578   them [] on backtracking, which is quite unfortunate in itself.
 1579- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1580
 1581prepare_global_variables(BDD) :-
 1582        clpb_next_id('$clpb_next_var', V0),
 1583        clpb_next_id('$clpb_next_node', N0),
 1584        bdd_nodes(BDD, Nodes),
 1585        foldl(max_variable_node, Nodes, V0-N0, MaxV0-MaxN0),
 1586        MaxV is MaxV0 + 1,
 1587        MaxN is MaxN0 + 1,
 1588        b_setval('$clpb_next_var', MaxV),
 1589        b_setval('$clpb_next_node', MaxN).
 1590
 1591max_variable_node(Node, V0-N0, V-N) :-
 1592        node_id(Node, N1),
 1593        node_varindex(Node, V1),
 1594        N is max(N0,N1),
 1595        V is max(V0,V1).
 1596
 1597/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1598   Fuse formulas that share the same variables into single conjunctions.
 1599- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1600
 1601ands_fusion(Ands0, Ands) :-
 1602        maplist(with_variables, Ands0, Pairs0),
 1603        keysort(Pairs0, Pairs),
 1604        group_pairs_by_key(Pairs, Groups),
 1605        pairs_values(Groups, Andss),
 1606        maplist(list_to_conjunction, Andss, Ands).
 1607
 1608with_variables(F, Vs-F) :-
 1609        term_variables(F, Vs0),
 1610        variables_in_index_order(Vs0, Vs).
 1611
 1612/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1613   If possible, separate variables into different sat/1 goals.
 1614   A formula F can be split in two if for two of its variables A and B,
 1615   taut((A^F)*(B^F) =:= F, 1) holds. In the first conjunct, A does not
 1616   occur, and in the second, B does not occur. We separate variables
 1617   until that is no longer possible. There may be a better way to do this.
 1618- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1619
 1620variables_separation(Fs0, Fs) :- separation_fixpoint(Fs0, [], Fs).
 1621
 1622separation_fixpoint(Fs0, Ds0, Fs) :-
 1623        phrase(variables_separation_(Fs0, Ds0, Rest), Fs1),
 1624        partition(anf_done, Fs1, Ds1, Fs2),
 1625        maplist(arg(1), Ds1, Ds2),
 1626        maplist(arg(1), Fs2, Fs3),
 1627        append(Ds0, Ds2, Ds3),
 1628        append(Rest, Fs3, Fs4),
 1629        sort(Fs4, Fs5),
 1630        sort(Ds3, Ds4),
 1631        (   Fs5 == [] -> Fs = Ds4
 1632        ;   separation_fixpoint(Fs5, Ds4, Fs)
 1633        ).
 1634
 1635anf_done(done(_)).
 1636
 1637variables_separation_([], _, []) --> [].
 1638variables_separation_([F0|Fs0], Ds, Rest) -->
 1639        (   { member(Done, Ds), F0 == Done } ->
 1640            variables_separation_(Fs0, Ds, Rest)
 1641        ;   { sat_rewrite(F0, F),
 1642              sat_bdd(F, BDD),
 1643              bdd_variables(BDD, Vs0),
 1644              exclude(universal_var, Vs0, Vs),
 1645              maplist(existential_(BDD), Vs, Nodes),
 1646              phrase(pairs(Nodes), Pairs),
 1647              group_pairs_by_key(Pairs, Groups),
 1648              phrase(groups_separation(Groups, BDD), ANFs) },
 1649            (   { ANFs = [_|_] } ->
 1650                list(ANFs),
 1651                { Rest = Fs0 }
 1652            ;   [done(F0)],
 1653                variables_separation_(Fs0, Ds, Rest)
 1654            )
 1655        ).
 1656
 1657
 1658existential_(BDD, V, Node) :- existential(V, BDD, Node).
 1659
 1660groups_separation([], _) --> [].
 1661groups_separation([BDD1-BDDs|Groups], OrigBDD) -->
 1662        { phrase(separate_pairs(BDDs, BDD1, OrigBDD), Nodes) },
 1663        (   { Nodes = [_|_] } ->
 1664            nodes_anfs([BDD1|Nodes])
 1665        ;   []
 1666        ),
 1667        groups_separation(Groups, OrigBDD).
 1668
 1669separate_pairs([], _, _) --> [].
 1670separate_pairs([BDD2|Ps], BDD1, OrigBDD) -->
 1671        (   { apply(*, BDD1, BDD2, And),
 1672              And == OrigBDD } ->
 1673            [BDD2]
 1674        ;   []
 1675        ),
 1676        separate_pairs(Ps, BDD1, OrigBDD).
 1677
 1678nodes_anfs([]) --> [].
 1679nodes_anfs([N|Ns]) --> { node_anf(N, ANF) }, [anf(ANF)], nodes_anfs(Ns).
 1680
 1681pairs([]) --> [].
 1682pairs([V|Vs]) --> pairs_(Vs, V), pairs(Vs).
 1683
 1684pairs_([], _) --> [].
 1685pairs_([B|Bs], A) --> [A-B], pairs_(Bs, A).
 1686
 1687/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1688   Set the Prolog flag clpb_residuals to bdd to obtain the BDD nodes
 1689   as residuals. Note that they cannot be used as regular goals.
 1690- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1691
 1692nodes([]) --> [].
 1693nodes([Node|Nodes]) -->
 1694        { node_var_low_high(Node, Var0, Low, High),
 1695          var_or_atom(Var0, Var),
 1696          maplist(node_projection, [Node,High,Low], [ID,HID,LID]),
 1697          var_index(Var0, VI) },
 1698        [ID-(v(Var,VI) -> HID ; LID)],
 1699        nodes(Nodes).
 1700
 1701
 1702node_projection(Node, Projection) :-
 1703        node_id(Node, ID),
 1704        (   integer(ID) -> Projection = node(ID)
 1705        ;   Projection = ID
 1706        ).
 1707
 1708/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1709   By default, residual goals are sat/1 calls of the remaining formulas,
 1710   using (mostly) algebraic normal form.
 1711- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1712
 1713sats([]) --> [].
 1714sats([A|As]) --> [clpb:sat(A)], sats(As).
 1715
 1716booleans([]) --> [].
 1717booleans([B|Bs]) --> boolean(B), { del_clpb(B) }, booleans(Bs).
 1718
 1719boolean(Var) -->
 1720        (   { get_attr(Var, clpb_omit_boolean, true) } -> []
 1721        ;   [clpb:sat(Var =:= Var)]
 1722        ).
 1723
 1724/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1725   Relate a formula to its algebraic normal form (ANF).
 1726- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1727
 1728formula_anf(Formula0, ANF) :-
 1729        parse_sat(Formula0, Formula),
 1730        sat_bdd(Formula, Node),
 1731        node_anf(Node, ANF).
 1732
 1733node_anf(Node, ANF) :-
 1734        node_xors(Node, Xors0),
 1735        maplist(maplist(monotonic_variable), Xors0, Xors),
 1736        maplist(list_to_conjunction, Xors, Conjs),
 1737        (   Conjs = [Var,C|Rest], clpb_var(Var) ->
 1738            foldl(xor, Rest, C, RANF),
 1739            ANF = (Var =\= RANF)
 1740        ;   Conjs = [One,Var,C|Rest], One == 1, clpb_var(Var) ->
 1741            foldl(xor, Rest, C, RANF),
 1742            ANF = (Var =:= RANF)
 1743        ;   Conjs = [C|Cs],
 1744            foldl(xor, Cs, C, ANF)
 1745        ).
 1746
 1747monotonic_variable(Var0, Var) :-
 1748        (   var(Var0), current_prolog_flag(clpb_monotonic, true) ->
 1749            Var = v(Var0)
 1750        ;   Var = Var0
 1751        ).
 1752
 1753clpb_var(Var) :- var(Var), !.
 1754clpb_var(v(_)).
 1755
 1756list_to_conjunction([], 1).
 1757list_to_conjunction([L|Ls], Conj) :- foldl(and, Ls, L, Conj).
 1758
 1759xor(A, B, B # A).
 1760
 1761eq_1(V) :- V == 1.
 1762
 1763node_xors(Node, Xors) :-
 1764        phrase(xors(Node), Xors0),
 1765        % we remove elements that occur an even number of times (A#A --> 0)
 1766        maplist(sort, Xors0, Xors1),
 1767        pairs_keys_values(Pairs0, Xors1, _),
 1768        keysort(Pairs0, Pairs),
 1769        group_pairs_by_key(Pairs, Groups),
 1770        exclude(even_occurrences, Groups, Odds),
 1771        pairs_keys(Odds, Xors2),
 1772        maplist(exclude(eq_1), Xors2, Xors).
 1773
 1774even_occurrences(_-Ls) :- length(Ls, L), L mod 2 =:= 0.
 1775
 1776xors(Node) -->
 1777        (   { Node == 0 } -> []
 1778        ;   { Node == 1 } -> [[1]]
 1779        ;   { node_var_low_high(Node, Var0, Low, High),
 1780              var_or_atom(Var0, Var),
 1781              node_xors(Low, Ls0),
 1782              node_xors(High, Hs0),
 1783              maplist(with_var(Var), Ls0, Ls),
 1784              maplist(with_var(Var), Hs0, Hs) },
 1785            list(Ls0),
 1786            list(Ls),
 1787            list(Hs)
 1788        ).
 1789
 1790list([]) --> [].
 1791list([L|Ls]) --> [L], list(Ls).
 1792
 1793with_var(Var, Ls, [Var|Ls]).
 1794
 1795/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1796   Global variables for unique node and variable IDs and atoms.
 1797- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1798
 1799make_clpb_var('$clpb_next_var') :- nb_setval('$clpb_next_var', 0).
 1800
 1801make_clpb_var('$clpb_next_node') :- nb_setval('$clpb_next_node', 0).
 1802
 1803make_clpb_var('$clpb_atoms') :-
 1804        empty_assoc(E),
 1805        nb_setval('$clpb_atoms', E).
 1806
 1807:- multifile user:exception/3. 1808
 1809user:exception(undefined_global_variable, Name, retry) :-
 1810        make_clpb_var(Name), !.
 1811
 1812clpb_next_id(Var, ID) :-
 1813        b_getval(Var, ID),
 1814        Next is ID + 1,
 1815        b_setval(Var, Next).
 1816
 1817clpb_atom_var(Atom, Var) :-
 1818        b_getval('$clpb_atoms', A0),
 1819        (   get_assoc(Atom, A0, Var) -> true
 1820        ;   put_attr(Var, clpb_atom, Atom),
 1821            put_attr(Var, clpb_omit_boolean, true),
 1822            put_assoc(Atom, A0, Var, A),
 1823            b_setval('$clpb_atoms', A)
 1824        ).
 1825
 1826/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1827   The variable attributes below are not used as constraints by this
 1828   library. Project remaining attributes to empty lists of residuals.
 1829
 1830   Because accessing these hooks is basically a cross-module call, we
 1831   must declare them public.
 1832- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1833
 1834:- public
 1835        clpb_hash:attr_unify_hook/2,
 1836        clpb_bdd:attribute_goals//1,
 1837        clpb_hash:attribute_goals//1,
 1838        clpb_omit_boolean:attr_unify_hook/2,
 1839        clpb_omit_boolean:attribute_goals//1,
 1840        clpb_atom:attr_unify_hook/2,
 1841        clpb_atom:attribute_goals//1. 1842
 1843clpb_hash:attr_unify_hook(_,_).  % this unification is always admissible
 1844
 1845/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1846   If a universally quantified variable is unified to a Boolean value,
 1847   it indicates that the formula does not hold for the other value, so
 1848   it is false.
 1849- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1850
 1851clpb_atom:attr_unify_hook(_, _) :- false.
 1852
 1853clpb_omit_boolean:attr_unify_hook(_,_).
 1854
 1855clpb_bdd:attribute_goals(_)          --> [].
 1856clpb_hash:attribute_goals(_)         --> [].
 1857clpb_omit_boolean:attribute_goals(_) --> [].
 1858clpb_atom:attribute_goals(_)         --> [].
 1859
 1860% clpb_hash:attribute_goals(Var) -->
 1861%         { get_attr(Var, clpb_hash, Assoc),
 1862%           assoc_to_list(Assoc, List0),
 1863%           maplist(node_portray, List0, List) }, [Var-List].
 1864
 1865% node_portray(Key-Node, Key-Node-ite(Var,High,Low)) :-
 1866%         node_var_low_high(Node, Var, Low, High).
 1867
 1868/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1869   Messages
 1870- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1871
 1872:- multifile prolog:message//1. 1873
 1874prolog:message(clpb(bounded)) -->
 1875        ['Using CLP(B) with bounded arithmetic may yield wrong results.'-[]].
 1876
 1877warn_if_bounded_arithmetic :-
 1878        (   current_prolog_flag(bounded, true) ->
 1879            print_message(warning, clpb(bounded))
 1880        ;   true
 1881        ).
 1882
 1883:- initialization(warn_if_bounded_arithmetic). 1884
 1885/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 1886   Sanbox declarations
 1887- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 1888
 1889:- multifile
 1890        sandbox:safe_global_variable/1,
 1891        sandbox:safe_primitive/1. 1892
 1893sandbox:safe_global_variable('$clpb_next_var').
 1894sandbox:safe_global_variable('$clpb_next_node').
 1895sandbox:safe_global_variable('$clpb_atoms').
 1896sandbox:safe_primitive(set_prolog_flag(clpb_residuals, _))