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