/*
 * lavabdd.h
 *
 *  Created on: 25 janv. 2010
 *      Author: nmaquet
 */

/*!
 *  \mainpage LaVaBDD : Lattice-Valued Binary Decision Diagrams Template Library
 *  \version 0.3
 *  \date March 28th, 2010
 *  \author Nicolas Maquet
 *  \section intro_sec Introduction
 *
 *  LaVaBDD is a C++ template library for Lattice-Valued Binary Decision Diagrams.
 *  The goal of this implementation is to be as simple as possible and reasonably efficient.
 *
 *  The LaVaBDD library is composed of a single C++ header file : lavabdd.h
 *
 *  \section install_sec Installation
 *
 *  The only dependency of the LaVaBDD library is Google's SparseHash.
 *  SparseHash is an efficient hash map implementation.
 *
 *  You can download it here : http://code.google.com/p/google-sparsehash/
 *
 *  \section howto_sec How to use the library
 *
 *  \htmlonly
 *
 *  Lattice-Valued Binary Decision Diagrams (LVBDD for short) encode functions of the
 *  form V(p<sub>1</sub>,...,p<sub>n</sub>) &rarr; L where:
 *  <li> p<sub>1</sub>,...,p<sub>n</sub> is a finite set of <em>propositional variables</em> </li>
 *  <li> V(p<sub>1</sub>,...,p<sub>n</sub>) is the set of <em>valuations</em> over p<sub>1</sub>,...,p<sub>n</sub> </li>
 *  <li> L is a <em>finite distributive lattice</em> </li>
 *
 *  \endhtmlonly
 *
 *  \subsection define_lattice Step 1 : define your lattice class
 *
 *  The first step to use LaVaBDD is to define your own lattice class.
 *  This class must follow the lava::Lattice interface defined in the package.
 *
 *  The lava::Lattice class is a class template. To define your lattice class,
 *  you must subclass this template as in the following example:
 *
 *
 *  \code

 #include <stdexcept> // for logic_error
 #include <bdd.h> // the BuDDy ROBDD library defines the "bdd" type
 #include "lavabdd.h"

 using namespace std;

 class BDDLattice : public lava::Lattice<bdd> {

 public:

     static bool less(const bdd& x, const bdd& y) { return (x & bdd_not(y)) == bddfalse; }
     static bool equal(const bdd& x, const bdd& y) { return x == y; }
     static void join(const bdd& x, const bdd& y, bdd& z) { z = x | y; }
     static void meet(const bdd& x, const bdd& y, bdd& z) { z = x & y; }
     static void rpc(const bdd& x, const bdd& y, bdd& z) { z = bdd_not(x) | y; }
     static void drpc(const bdd& x, const bdd& y, bdd& z) { z = y & bdd_not(x); }
     static void top(bdd& x) { x = bddtrue; }
     static void bot(bdd& x) { x = bddfalse; }
     static size_t hash(const bdd& x) { return lava::hash_int(x.id()); }
     static void print(const bdd& x, ostream& out) { out << x; }
     static int size(const bdd& x) { return bdd_nodecount(x); }

     // not implemented here for the sake of brevity
     static int list_size(const list<bdd>& dds) { throw logic_error("not implemented"); }

 };

 *  \endcode
 *
 *  Here is another, simpler example.
 *
 *  \code

 #include <algorithm> // for min() and max()
 #include "lavabdd.h"

 using namespace std;

 typedef unsigned char byte;

 class ByteLattice : public lava::Lattice<byte> {

 public:

     static bool less(const byte& x, const byte& y) { x < y; }
     static bool equal(const byte& x, const byte& y) { return x == y; }
     static void join(const byte& x, const byte& y, byte& z) { z = max(x, y); }
     static void meet(const byte& x, const byte& y, byte& z) { z = min(x, y); }
     static void rpc(const byte& x, const byte& y, byte& z) { z = y; }
     static void drpc(const byte& x, const byte& y, byte& z) { z = y; }
     static void top(byte& x) { x = 255; }
     static void bot(byte& x) { x = 0; }
     static size_t hash(const byte& x) { return lava::hash_int((int) x); }
     static void print(const byte& x, ostream& out) { out << x; }
     static int size(const byte& x) { return 1; }
     static int list_size(const list<byte>& xs) { return xs.size(); }

 };

 *  \endcode
 *
 *  \subsection chose_normal_form Step 2 : Choose an LVBDD normal form
 *
 *  There are three available normal forms for LVBDD defined in the package.
 *  \li the Unshared Normal Form (lava::UNF)
 *  \li the \e meet-semantics Shared Normal Form (lava::SNF_MEET)
 *  \li the \e join-semantics Shared Normal Form (lava::SNF_JOIN)
 *
 *  For a full comparison of each, see the paper. Here is a quick bottom line:
 *  \li both SNF are potentially exponentially more compact than UNF
 *  \li the meet and join are both quadratic-time for UNF
 *  \li the meet for \e meet-semantics SNF is polynomial, join is exponential
 *  \li the join for \e join-semantics SNF is polynomial, meet is exponential
 *
 *  \subsection define_lvbdd_type Step 3 : Define your LVBDD type
 *
 *  Since the LVBDD class template takes two parameters, it is convenient to
 *  define your own LVBDD type. This is done simply as follows :
 *
 *  \code
 *
 *  class BDDLattice; // as defined above
 *
 *  typedef lava::LVBDD<BDDLattice, lava::SNF_MEET> LVBDD;
 *
 *  \endcode
 *
 *  \subsection init_LVBDD Step 4 : Initialize your LVBDD type
 *
 *  Every LVBDD type \b must be initialized before use. See lava::LVBDD::init.
 *
 *  \subsection use_LVBDD Step 5 : You're good to go !
 *
 *  You can now use your LVBDD type. See lava::LVBDD for a complete reference on the
 *  available methods.
 *
 *  A good starting point is lava::LVBDD::terminal and lava::LVBDD::literal.
 *
 *  To help you, here is a short example.
 *
 *  \code

 #include <bdd.h> // BuDDy
 #include "lavabdd.h"

 using namespace std;

 class BDDLattice; // as defined above

 typedef lava::LVBDD<BDDLattice, lava::SNF_MEET> LVBDD;

 int main() {
     LVBDD::init();              // initialize the LVBDD type
     bdd_init(10000, 10000);     // initialize the BuDDy library
     bdd_setvarnum(4);

     LVBDD p1 = LVBDD::literal(1, true);
     LVBDD p2 = LVBDD::literal(2, true);
     LVBDD p3 = LVBDD::literal(3, true);
     LVBDD c1 = LVBDD::terminal(bdd_ithvar(1));
     LVBDD c2 = LVBDD::terminal(bdd_ithvar(2));
     LVBDD c3 = LVBDD::terminal(bdd_ithvar(3));

     LVBDD dd = (p1 | c1) & (p2 | c2) & (p3 | c3);

     dd.print_graphviz();
 }

 *  \endcode
 *
 *  \section sect_changelog Changelog
 *
 *  \subsection version_0_dot_1 Version 0.1 (Released February 23rd 2010)
 *
 *  Initial version.
 *
 *  \subsection version_0_dot_2 Version 0.2 (Released February 26th 2010)
 *
 *  \li Added a changelog to the documentation ;-)
 *  \li Lots of bugfixes / typos
 *  \li lava::hash_int is now a global function and not a static member of lava:LVBDD
 *  \li Added a definition for operator<< for writing LVBDD on output streams
 *  \li Added a definition for operator< in lava::NodePtr to allow std::set<NodePtr>
 *
 *  \subsection version_0_dot_3 Version 0.3 (Released March 28th 2010)
 *  \li a few bugfixes
 *  \li a few minor optimizations
 *  \li the private member Node<L,F>::_refcount is now "mutable", implementation is cleaner as a result
 *  \li memoization tables are now to be cleared manually; this means that, by default,
 *      <b>everything</b> is memoized and kept in memory until manually cleared
 *  \li added the static method lava::LVBDD::clear_memoization_tables
 *  \li added the static method lava::LVBDD::memoization_tables_size
 */

#ifndef __LAVABDD_H__
#define __LAVABDD_H__

/*!
 *  \brief Define this to use Google's *dense* hashmap (fast but memory hungry).
 *  Don't define it if you want to use the *sparse* hashmap (slower but memory-efficient).
 */
#define USING_DENSE_HASH_MAP

/*!
 *  \brief Initial reference count for LVBDD nodes. Should be zero.
 *  Increase the value to track reference counting leaks.
 */
#define INITIAL_REFCOUNT 0

#include <list>
#include <stdexcept>
#include <iomanip>
#include <iostream>
#include <fstream>
#include <map>
#include <set>
#include <vector>

#ifdef USING_DENSE_HASH_MAP
#include <google/dense_hash_map>
using google::dense_hash_map;
#else
#include <google/sparse_hash_map>
using google::sparse_hash_map;
#endif

using namespace std;

/*!
 *  \brief Everything in the LaVaBDD package is defined in the "lava" namespace.
 */
namespace lava {

/*!
 *  \brief Pass this value to the LVBDD template to obtain a Unshared Normal Form LVBDD.
 */
const int UNF = 0;

/*!
 *  \brief Pass this value to the LVBDD template to obtain a meet-semantics
 *         Shared Normal Form LVBDD.
 */
const int SNF_MEET = 1;

/*!
 *  \brief Pass this value to the LVBDD template to obtain a join-semantics
 *         Shared Normal Form LVBDD.
 */
const int SNF_JOIN = 2;

/* =====================================================================================*/

/*!
 *   \brief Simple hash function for int values.
 *          It is used internally to hash pointers and indexes.
 */
inline int hash_int(int key) {
    key = ~key + (key << 15);
    key = key ^ (key >> 12);
    key = key + (key << 2);
    key = key ^ (key >> 4);
    key = key * 2057;
    key = key ^ (key >> 16);
    return key;
}

/* =====================================================================================*/

/*!
 *  \brief The Lattice class template represents a finite distributive lattice.
 */
template<typename V>
class Lattice {

public:

    typedef V ValueType; // This type is propagated as V in Node, NodePtr, and LVBDD.

    /*!
     *  \brief The partial order of the lattice.
     *         It must be antisymmetric, reflexive and transitive.
     */
    static bool less(const V& x, const V& y) {
        throw logic_error("Lattice::less not implemented.");
    }

    /*!
     *  \brief Is equivalent to "less(x,y) and less(y,x)" but ideally more efficient.
     */
    static bool equal(const V& x, const V& y) {
        throw logic_error("Lattice::equal not implemented.");
    }

    /*!
     *  \brief Least upper bound operator of the lattice (or "join").
     */
    static void join(const V& x, const V& y, V& z) {
        throw logic_error("Lattice::join not implemented.");
    }

    /*!
     *  \brief Greatest lower bound operator of the lattice (or "meet").
     */
    static void meet(const V& x, const V& y, V& z) {
        throw logic_error("Lattice::meet not implemented.");
    }

    /*!
     *  \brief Relative pseudocomplementation. Assumes that less(y, x) is true.
     *         Computes the greatest value z such that "x meet z = y"
     */
    static void rpc(const V& x, const V& y, V& z) {
        throw logic_error("Lattice::rpc not implemented.");
    }

    /*!
     *  \brief Dual relative pseudocomplementation. Assumes that less(x, y) is true.
     *         Computes the smallest value z such that "x join z = y"
     */
    static void drpc(const V& x, const V& y, V& z) {
        throw logic_error("Lattice::drpc not implemented.");
    }

    /*!
     *  \brief Returns the "top" element of the lattice.
     */
    static void top(V& x) {
        throw logic_error("Lattice::top not implemented.");
    }

    /*!
     *  \brief Returns the "bottom" element of the lattice.
     */
    static void bot(V& x) {
        throw logic_error("Lattice::bot not implemented.");
    }

    /*!
     *  \brief Computes an approriate hash value for a lattice.
     *         It must be the case that "equal(x, y)" implies "hash(x) == hash(y)".
     *         Note that you can first compute an integer and then use LVBDD::hash_int.
     *         A poor implementation of this function will lead to bad performance
     *         due to hash map collisions.
     *         See http://www.sgi.com/tech/stl/HashFunction.html for more information.
     */
    static size_t hash(const V& x) {
        throw logic_error("Lattice::hash not implemented.");
    }

    /*!
     *   \brief Prints a textual representation of a lattice value.
     */
    static void print(const V& x, ostream& out) {
        throw logic_error("Lattice::print not implemented.");
    }

    /*!
     *   \brief Returns a reasonable measure of the size of the encoding of a lattice value.
     *          For instance, a reasonable measure for the size of a lattice value encoded
     *          as an ROBDD is the number of its nodes.
     */
    static int size(const V& x) {
        throw logic_error("Lattice::size not implemented.");
    }

    /*!
     *   \brief Returns a reasonable measure of the size of a list of lattice values.
     *          In many cases, it is valid to simply return the sum of Lattice::size
     *          for each value in the list. If memory is shared among lattice values, this
     *          is not accurate (e.g. ROBDD), hence the need for this function.
     */
    static int list_size(const list<V>& l) {
        throw logic_error("Lattice::list_size not implemented.");
    }

};

/* ==================================================================================== */

template<class L, int F> class Node;

/*!
 *  \brief The NodePtr class provides automatic reference-counting for Node pointers.
 *
 *         It encapsulates a Node* variable and can be used exactly as such, thanks to
 *         the overloaded operator definitions.
 */
template<class L, int F> class NodePtr {

private:

    typedef typename L::ValueType V;
    typedef Node<L, F> Node;

    Node* _node;

public:

    NodePtr() {
        _node = NULL;
    }

    /*!
     *  \brief Constructor. It is not explicit so Node* values (such as NULL) will be
     *         converted automatically.
     */
    NodePtr(Node* node) {
        _node = node;
        if (_node)
            _node->_ref();
    }

    /*!
     *  \brief Destructor.
     */
    ~NodePtr() {
        if (_node)
            _node->_unref();
    }

    /*!
     *  \brief Copy constructor.
     */
    NodePtr(const NodePtr& other) {
        _node = other._node;
        if (_node)
            _node->_ref();
    }

    /*!
     *  \brief Assignment operator.
     */
    NodePtr& operator=(const NodePtr& other) {
        if (this != &other) {
            if (_node)
                _node->_unref();
            _node = other._node;
            if (_node)
                _node->_ref();
        }
        return *this;
    }

    /*!
     *  \brief Pointer dereference. Throws std::logic_error on NULL access.
     */
    Node& operator*() const {
        if (_node == NULL)
            throw logic_error("NodePtr : NULL pointer access");
        return *_node;
    }

    /*!
     *  \brief Pointer member access operator. Throws std::logic_error on NULL access.
     */
    Node* operator->() const {
        if (_node == NULL)
            throw logic_error("NodePtr : NULL pointer access");
        return _node;
    }

    /*!
     *  \brief Returns the encapsulated pointer.
     *         You should not use this unless you know what you're doing.
     *         Wrong use of this will mess up reference counting.
     */
    Node* get_ptr() const {
        return _node;
    }

    /*!
     *  \brief Pointer equality.
     */
    bool operator==(const NodePtr& other) const {
        return _node == other._node;
    }

    /*!
     *  \brief Pointer inequality.
     */
    bool operator!=(const NodePtr& other) const {
        return _node != other._node;
    }

    /*!
     *  \brief <em>Less than</em> comparison on pointers. Allows the use of set<NodePtr>.
     */
    bool operator<(const NodePtr& other) const {
        return _node < other._node;
    }

    /*!
     *  \brief Conversion to bool. Returns true iff encapsulated pointer is non-NULL.
     */
    operator bool() const {
        return _node ? true : false;
    }

    /*!
     *  \brief Releases the encapsulated pointer and sets it to NULL.
     */
    void nullify() {
        if (_node)
            _node->_unref();
        _node = NULL;
    }

    /*!
     *  \brief Releases the current encapsulated pointer and grabs a new one.
     */
    void grab(Node *node) {
        if (_node)
            _node->_unref();
        _node = node;
        if (_node)
            _node->_ref();
    }

};

template<class L, int F> class LVBDD;

/*!
 *   \brief The Node class represents an LVBDD node.
 *
 *          To avoid memory leaks and corruption, this package will not let you access
 *          this class directly; see NodePtr.
 */
template<class L, int F> class Node {

private:

    typedef typename L::ValueType V;

    friend class NodePtr<L, F> ;
    typedef NodePtr<L, F> NodePtr;
    typedef LVBDD<L, F> LVBDD;

    mutable int _refcount; /* you can change the RC of const Node references */
    int _index;
    V _value;
    size_t _hash;

    NodePtr _lo;
    NodePtr _hi;

    void _unref() {
        if (_refcount <= 0) {
            throw runtime_error("unref on orphan Node");
        }
        _refcount--;
        if (_refcount == 0) {
            LVBDD::_add_orphan(this);
            _lo.nullify();
            _hi.nullify();
            _index = -1;
        }
    }

    void _ref() {
        _refcount++;
    }

public:

    /*
     *  Default constructor. Do not use.
     */
    Node() {
    }

    /*!
     *  \brief Constructor.
     *         Do not use this unless you know what you're doing.
     *         Wrong use of this constructor will lead to normal form violation.
     *  \param [in] index : An index value strictly larger than zero for non-terminal nodes and
     *                 equal to zero for terminal nodes.
     *  \param [in] val : The value labeling the node.
     *  \param [in] lo : The node's low-child (must be NULL for terminal nodes, non-NULL for
     *              nonterminal nodes)
     *  \param [in] hi : The node's high-child (must be NULL for terminal nodes, non-NULL for
     *              nonterminal nodes)
     */
    Node(int index, const V& val, const NodePtr& lo, const NodePtr& hi) :
        _refcount(INITIAL_REFCOUNT), _index(index), _value(val), _lo(lo), _hi(hi) {
        _hash = L::hash(_value) ^ hash_int(_index);
        if (lo)
            _hash ^= hash_int(((size_t) lo.get_ptr()) * 17);
        if (hi)
            _hash ^= hash_int(((size_t) hi.get_ptr()));
    }

    /*!
     *  \brief Returns the nodes's index
     *  (zero for terminal nodes, and larger for nonterminal nodes).
     */
    int index() const {
        return _index;
    }

    /*!
     *  \brief Returns a hash value for the node. This is the XOR of four hashes :
     *  the index, the value and the lo and hi child pointers.
     */
    size_t hash() const {
        return _hash;
    }

    /*!
     *  \brief Returns the node's reference count.
     */
    int refcount() const {
        return _refcount;
    }

    /*!
     *  \brief Returns the node's label value.
     */
    const V& value() const {
        return _value;
    }

    /*!
     *  \brief Returns the node's label value.
     */
    void value(V& value) const {
        value = _value;
    }

    /*!
     *  \brief Returns the node's low-child.
     */
    NodePtr lo() const {
        return _lo;
    }

    /*!
     *  \brief Returns the node's high-child.
     */
    NodePtr hi() const {
        return _hi;
    }

    /*!
     *  \brief Node equality operator. Returns true iff the compared Node has the same
     *  index, the same value and the same low- and high-child pointers. Runs in O(1).
     */
    bool operator==(const Node& other) const {
        return _hash == other._hash && _index == other._index && _lo == other._lo && _hi
                == other._hi && L::equal(_value, other._value);
    }

    /*!
     *  \brief Node inequality operator. Same remark as Node::operator==.
     */
    bool operator!=(const Node& other) const {
        return not (*this == other);
    }

    /*
     * Needless to say, be careful with this one ;-)
     */
    void set_refcount(int refcount) {
        _refcount = refcount;
    }

    /*
     *  Do not use. This is used as a hashmap sentinel key.
     */
    static Node deleted_key() {
        Node key;
        key._refcount = -1;
        key._index = -1;
        key._hash = -1;
        return key;
    }

    /*
     * Do not use. This is used as a hashmap sentinel key.
     */
    static Node empty_key() {
        Node key;
        key._refcount = -2;
        key._index = -2;
        key._hash = -2;
        return key;
    }

    /*!
     *  \brief Outputs a DOT diagram representation of all descendants of a node.
     *  \param [in] out : The output stream on which to write.
     *  \param [in] root : The root node of which we want to print the descendants.
     *                     If NULL, all nodes are printed.
     */
    static void print_graphviz(ostream& out = cout, NodePtr root = NULL) {
        set<Node*> my_nodes;
        if (not root) {
            typename LVBDD::HashMap::const_iterator it;
            for (it = LVBDD::_hash_map.begin(); it != LVBDD::_hash_map.end(); it++) {
                my_nodes.insert(it->second);
            }
        }
        else {
            root->_descendants(my_nodes);
        }
        typename set<Node*>::iterator it;
        out << "digraph G {" << endl;
        for (it = my_nodes.begin(); it != my_nodes.end(); it++) {
            if ((*it)->index() == 0) {
                out << "\"" << *it << "\" [label=\"";
                L::print((*it)->value(), out);
                if (*it == root.get_ptr())
                    out << " (" << (*it)->refcount() - 1 << ")";
                else
                    out << " (" << (*it)->refcount() << ")";
                out << "\", shape=box];" << endl;
            }
            else {
                out << "\"" << *it << "\" [label=\"" << (*it)->index() << " : ";
                L::print((*it)->value(), out);
                if (*it == root.get_ptr())
                    out << " (" << (*it)->refcount() - 1 << ")";
                else
                    out << " (" << (*it)->refcount() << ")";
                out << "\"];" << endl;
            }
        }
        for (it = my_nodes.begin(); it != my_nodes.end(); it++) {
            if ((*it)->lo()) {
                out << "\"" << *it << "\" -> " << "\"" << ((*it)->lo().get_ptr());
                out << "\" [style=dotted, arrowhead=none];" << endl;
            }
            if ((*it)->lo()) {
                out << "\"" << *it << "\" -> " << "\"" << ((*it)->hi().get_ptr());
                out << "\" [arrowhead=none];" << endl;
            }
        }
        out << "}" << endl;
    }

    void _descendants(set<Node*>& nodes) {
        nodes.clear();
        set<Node*> to_visit;
        typename set<Node*>::iterator it;
        to_visit.insert(this);
        Node* n;
        while (to_visit.size() > 0) {
            it = to_visit.begin();
            n = *it;
            to_visit.erase(it);
            if (nodes.find(n) == nodes.end()) {
                nodes.insert(n);
                if (n->index() > 0) {
                    to_visit.insert(n->lo().get_ptr());
                    to_visit.insert(n->hi().get_ptr());
                }
            }
        }
    }

};

/* ==================================================================================== */

/*!
 *  \brief The LVBDD class template represents a Lattice-Valued Binary Decision Diagram.
 *
 *  See lava::UNF, lava::SNF_JOIN and  and laval::SNF_MEET.
 */
template<class L, int F>
class LVBDD {

private:

    typedef typename L::ValueType V;

    /* ================================================================================ */

    struct NodeHashFunction {
        size_t operator()(Node<L, F>* const & n) const {
            return n->hash();
        }
    };

    struct NodeEqualFunction {
        bool operator()(Node<L, F>* const & n1, Node<L, F>* const & n2) const {
            return (*n1) == (*n2);
        }
    };

    /* ================================================================================ */

    struct Quadruple {
        NodePtr<L, F> n1, n2;
        V d1, d2;
    };

    struct QuadrupleLT {
        bool operator()(const Quadruple& a, const Quadruple& b) const {
            if (a.n1->hash() != b.n1->hash())
                return a.n1->hash() < b.n1->hash();
            if (a.n2->hash() != b.n2->hash())
                return a.n2->hash() < b.n2->hash();
            if (L::hash(a.d1) != L::hash(b.d1))
                return L::hash(a.d1) < L::hash(b.d1);
            if (L::hash(a.d2) != L::hash(b.d2))
                return L::hash(a.d2) < L::hash(b.d2);
            return false;
        }
    };

    /* ================================================================================ */

    class QuadrupleMemoizer {

        typedef Node<L, F> Node;
        typedef NodePtr<L, F> NodePtr;

        typedef map<const Quadruple, NodePtr, QuadrupleLT> MAP;
        MAP _map;

    public:

        NodePtr get(const NodePtr& n1, const NodePtr& n2, const V& d1, const V& d2) const {
            Quadruple q1, q2;
            q1.n1 = n1;
            q2.n1 = n2;
            q1.n2 = n2;
            q2.n2 = n1;
            q1.d1 = d1;
            q2.d1 = d2;
            q1.d2 = d2;
            q2.d2 = d1;
            typename MAP::const_iterator it;
            it = _map.find(q1);
            if (it != _map.end())
                return it->second;
            it = _map.find(q2);
            if (it != _map.end())
                return it->second;
            return NULL;
        }

        void set(const NodePtr& n1, const NodePtr& n2, const V& d1, const V& d2,
                const NodePtr& memo) {
            Quadruple q;
            q.n1 = n1;
            q.n2 = n2;
            q.d1 = d1;
            q.d2 = d2;
            _map[q] = memo;
        }

        void clear() {
            _map.clear();
        }

    };

    /* ================================================================================ */

    class NodeNodeMemoizer {

        typedef Node<L, F> Node;
        typedef NodePtr<L, F> NodePtr;

        map<pair<NodePtr, NodePtr> , NodePtr> _map;

    public:

        NodePtr get(const NodePtr& n1, const NodePtr& n2) const {
            pair<NodePtr, NodePtr> p1(n1, n2);
            pair<NodePtr, NodePtr> p2(n2, n1);
            typename map<pair<NodePtr, NodePtr> , NodePtr>::const_iterator it;
            NodePtr result;
            it = _map.find(p1);
            if (it != _map.end()) {
                result = it->second;
            }
            else {
                it = _map.find(p2);
                if (it != _map.end())
                    result = it->second;
                else
                    result = NULL;
            }
            return result;
        }

        void set(const NodePtr& n1, const NodePtr& n2, const NodePtr& memo) {
            pair<NodePtr, NodePtr> p(n1, n2);
            _map[p] = memo;
        }

        void clear() {
            _map.clear();
        }

    };

    /* ================================================================================ */

    struct NodeValueLT {
        typedef Node<L, F> Node;
        typedef NodePtr<L, F> NodePtr;

        bool operator()(const pair<NodePtr, V>& a, const pair<NodePtr, V>& b) const {
            if (a.first == b.first)
                return L::hash(a.second) < L::hash(b.second);
            else
                return a.first.get_ptr() < b.first.get_ptr();
        }
    };

    /* ================================================================================ */

    class NodeValueMemoizer {

        typedef Node<L, F> Node;
        typedef NodePtr<L, F> NodePtr;
        map<pair<NodePtr, V> , NodePtr, NodeValueLT> _map;

    public:

        NodePtr get(const NodePtr& n, const V& value) const {
            pair<NodePtr, V> p(n, value);
            typename map<pair<NodePtr, V> , NodePtr, NodeValueLT>::const_iterator it;
            it = _map.find(p);
            if (it != _map.end())
                return it->second;
            else
                return NULL;
        }

        void set(const NodePtr& n, const V& value, const NodePtr& memo) {
            pair<NodePtr, V> p(n, value);
            _map[p] = memo;
        }

        void clear() {
            _map.clear();
        }

    };

    friend class Node<L, F> ;

    typedef Node<L, F> Node;
    typedef NodeHashFunction HashNode;
    typedef NodeEqualFunction EqualNode;
    typedef NodePtr<L, F> NodePtr;

#ifdef USING_DENSE_HASH_MAP
    typedef dense_hash_map<Node*, Node*, HashNode, EqualNode> HashMap;
#else
    typedef sparse_hash_map<Node*, Node*, HashNode, EqualNode> HashMap;
#endif

    NodePtr _root;

    static HashMap _hash_map;
    static list<Node*> _orphans;

    static int _mk_call_count;
    static int _orphan_fetches;
    static int _node_allocations;

    static NodeNodeMemoizer _memo_nn;
    static NodeValueMemoizer _memo_nv;
    static QuadrupleMemoizer _memo_quad;

    static NodePtr _unf_apply_rec(const NodePtr& n1, const NodePtr& n2, bool meet) {
        if (not n1 or not n2)
            abort();
        NodePtr result = _memo_nn.get(n1, n2);
        V value;
        L::top(value);
        if (not result) {
            if (n1->index() == 0 and n2->index() == 0) {
                V value = n1->value();
                if (meet)
                    L::meet(value, n2->value(), value);
                else
                    L::join(value, n2->value(), value);
                result = _mk(0, value, NULL, NULL);
            }
            else if (n1->index() == n2->index()) {
                NodePtr lo = _unf_apply_rec(n1->lo(), n2->lo(), meet);
                NodePtr hi = _unf_apply_rec(n1->hi(), n2->hi(), meet);
                if (lo == hi) {
                    result = lo;
                }
                else {
                    result = _mk(n1->index(), value, lo, hi);
                }
            }
            else {
                const NodePtr *new_n1 = &n1, *new_n2 = &n2;
                if (n1->index() < n2->index()) {
                    new_n1 = &n2;
                    new_n2 = &n1;
                }
                NodePtr lo = _unf_apply_rec((*new_n1)->lo(), *new_n2, meet);
                NodePtr hi = _unf_apply_rec((*new_n1)->hi(), *new_n2, meet);
                if (lo == hi) {
                    result = lo;
                }
                else {
                    result = _mk((*new_n1)->index(), value, lo, hi);
                }
            }
            _memo_nn.set(n1, n2, result);
        }
        return result;
    }

    static NodePtr _snf_rpc(const NodePtr& n, const V& d) {
        V value = n->value();
        if (F == SNF_MEET) {
            L::rpc(d, value, value);
        }
        else {
            L::drpc(d, value, value);
        }
        if (n->index() == 0)
            return _mk(0, value, NULL, NULL);
        V x = n->lo()->value();
        if (F == SNF_MEET) {
            L::join(x, n->hi()->value(), x);
            L::meet(value, x, value);
        }
        else {
            L::meet(x, n->hi()->value(), x);
            L::join(value, x, value);
        }
        return _mk(n->index(), value, n->lo(), n->hi());
    }

    static NodePtr _snf_easy_cst(const NodePtr& n, const V& value) {
        NodePtr result = _memo_nv.get(n, value);
        if (result) {
            return result;
        }
        if ((F == SNF_MEET and L::less(n->value(), value)) or (F == SNF_JOIN and L::less(value,
                n->value()))) {
            result = n;
        }
        else {
            V new_value = value;
            if (F == SNF_MEET)
                L::meet(new_value, n->value(), new_value);
            else
                L::join(new_value, n->value(), new_value);
            if (n->index() == 0) {
                result = _mk(0, new_value, NULL, NULL);
            }
            else {
                NodePtr lo, hi, new_lo, new_hi;
                lo = _snf_easy_cst(n->lo(), value);
                hi = _snf_easy_cst(n->hi(), value);
                if (lo == hi) {
                    V label = lo->value();
                    if (F == SNF_MEET)
                        L::meet(label, value, label);
                    else
                        L::join(label, value, label);
                    result = _mk(lo->index(), label, lo->lo(), lo->hi());
                }
                else {
                    new_lo = _snf_rpc(lo, value);
                    new_hi = _snf_rpc(hi, value);
                    if (new_lo == new_hi)
                        abort();
                    result = _mk(n->index(), new_value, new_lo, new_hi);
                }
            }
        }
        _memo_nv.set(n, value, result);
        return result;
    }

    static NodePtr _snf_easy_rec(const NodePtr& n1, const NodePtr& n2) {
        if (not n1 or not n2)
            abort();
        NodePtr result = _memo_nn.get(n1, n2);
        if (result) {
            return result;
        }
        if (n1->index() == 0)
            result = _snf_easy_cst(n2, n1->value());
        else if (n2->index() == 0)
            result = _snf_easy_cst(n1, n2->value());
        else if (n1->index() == n2->index()) {
            NodePtr lo = _snf_easy_rec(n1->lo(), n2->lo());
            NodePtr hi = _snf_easy_rec(n1->hi(), n2->hi());
            V cst_value = n1->value();
            if (F == SNF_MEET)
                L::meet(cst_value, n2->value(), cst_value);
            else
                L::join(cst_value, n2->value(), cst_value);
            NodePtr new_lo = _snf_easy_cst(lo, cst_value);
            NodePtr new_hi = _snf_easy_cst(hi, cst_value);
            if (new_lo == new_hi) {
                result = new_lo;
            }
            else {
                cst_value = new_lo->value();
                if (F == SNF_MEET) {
                    L::join(cst_value, new_hi->value(), cst_value);
                }
                else {
                    L::meet(cst_value, new_hi->value(), cst_value);
                }
                lo = _snf_rpc(new_lo, cst_value);
                hi = _snf_rpc(new_hi, cst_value);
                if (lo == hi) {
                    result = lo;
                }
                else {
                    result = _mk(n1->index(), cst_value, lo, hi);
                }
            }
        }
        else {
            const NodePtr *new_n1 = &n1, *new_n2 = &n2;
            if (n1->index() < n2->index()) {
                new_n1 = &n2;
                new_n2 = &n1;
            }
            NodePtr lo = _snf_easy_rec((*new_n1)->lo(), *new_n2);
            NodePtr hi = _snf_easy_rec((*new_n1)->hi(), *new_n2);
            V cst_value = (*new_n1)->value();
            NodePtr new_lo = _snf_easy_cst(lo, cst_value);
            NodePtr new_hi = _snf_easy_cst(hi, cst_value);
            if (new_lo == new_hi) {
                result = new_lo;
            }
            else {
                cst_value = new_lo->value();
                if (F == SNF_MEET) {
                    L::join(cst_value, new_hi->value(), cst_value);
                }
                else {
                    L::meet(cst_value, new_hi->value(), cst_value);
                }
                lo = _snf_rpc(new_lo, cst_value);
                hi = _snf_rpc(new_hi, cst_value);
                if (lo == hi) {
                    result = lo;
                }
                else {
                    result = _mk((*new_n1)->index(), cst_value, lo, hi);
                }
            }
        }
        _memo_nn.set(n1, n2, result);
        return result;
    }

    static NodePtr _snf_hard_rec(const NodePtr& n1, const NodePtr& n2, const V& d1, const V& d2) {
        NodePtr result = _memo_quad.get(n1, n2, d1, d2);
        if (result) {
            return result;
        }
        V new_d1 = d1, new_d2 = d2;
        if (F == SNF_MEET) {
            L::meet(new_d1, n1->value(), new_d1);
            L::meet(new_d2, n2->value(), new_d2);
        }
        else {
            L::join(new_d1, n1->value(), new_d1);
            L::join(new_d2, n2->value(), new_d2);
        }
        if (n1->index() == 0 and n2->index() == 0) {
            V value = new_d1;
            if (F == SNF_MEET)
                L::join(value, new_d2, value);
            else
                L::meet(value, new_d2, value);
            result = _mk(0, value, NULL, NULL);
            _memo_quad.set(n1, n2, d1, d2, result);
            return result;
        }
        else if (n1->index() == n2->index()) {
            NodePtr lo = _snf_hard_rec(n1->lo(), n2->lo(), new_d1, new_d2);
            NodePtr hi = _snf_hard_rec(n1->hi(), n2->hi(), new_d1, new_d2);
            if (lo == hi) {
                result = lo;
            }
            else {
                V value = lo->value();
                if (F == SNF_MEET) {
                    L::join(value, hi->value(), value);
                }
                else {
                    L::meet(value, hi->value(), value);
                }
                NodePtr new_lo = _snf_rpc(lo, value);
                NodePtr new_hi = _snf_rpc(hi, value);
                result = _mk(n1->index(), value, new_lo, new_hi);
            }
            _memo_quad.set(n1, n2, d1, d2, result);
            return result;
        }
        else {
            const NodePtr *new_n1 = &n1, *new_n2 = &n2;
            bool swapped = false;
            if (n1->index() < n2->index()) {
                new_n1 = &n2;
                new_n2 = &n1;
                swapped = true;
                swap(new_d1, new_d2);
            }
            NodePtr lo = _snf_hard_rec((*new_n1)->lo(), *new_n2, new_d1, new_d2);
            NodePtr hi = _snf_hard_rec((*new_n1)->hi(), *new_n2, new_d1, new_d2);
            if (lo == hi) {
                result = lo;
            }
            else {
                V value = lo->value();
                if (F == SNF_MEET) {
                    L::join(value, hi->value(), value);
                }
                else {
                    L::meet(value, hi->value(), value);
                }
                NodePtr new_lo = _snf_rpc(lo, value);
                NodePtr new_hi = _snf_rpc(hi, value);
                result = _mk((*new_n1)->index(), value, new_lo, new_hi);
            }
            if (swapped)
                _memo_quad.set(n2, n1, d1, d2, result);
            else
                _memo_quad.set(n1, n2, d1, d2, result);
            return result;
        }
    }

    static NodePtr _fetch_orphan() {
        Node* result;
        typename list<Node*>::iterator it;
        for (it = _orphans.begin(); it != _orphans.end();) {
            if ((*it)->refcount() == 0)
                break;
            it = _orphans.erase(it);
        }
        if (it == _orphans.end()) {
            result = NULL;
        }
        else {
            Node* n = *it;
            _orphans.erase(it);
            _hash_map.erase(n);
            result = n;
        }
        return NodePtr(result);
    }

    static NodePtr _mkbot() {
        V v;
        L::bot(v);
        return _mk(0, v, NULL, NULL);
    }

    static NodePtr _mktop() {
        V v;
        L::top(v);
        return _mk(0, v, NULL, NULL);
    }

    static void _add_orphan(Node* n) {
        if (n == NULL or n->refcount() != 0)
            throw runtime_error("invalid orphan node");
        _hash_map.erase(n);
        _orphans.push_back(n);
    }

    static NodePtr _mk(int index, const V& value, const NodePtr& lo, const NodePtr& hi) {

        typename HashMap::const_iterator it;

        _mk_call_count++;

        NodePtr result;
        if (lo == hi and lo) {
            abort();
        }
        else {
            Node key(index, value, lo, hi);
            it = _hash_map.find(&key);
            if (it == _hash_map.end()) {
                result = _fetch_orphan();
                if (not result) {
                    _node_allocations++;
                    result.grab(new Node(index, value, lo, hi));
                }
                else {
                    if (result->refcount() != 1)
                        abort();
                    _orphan_fetches++;
                    *result = key;
                    result->set_refcount(1);
                }
                _hash_map[result.get_ptr()] = result.get_ptr();
            }
            else {
                result.grab(it->second);
            }
        }

        return result;
    }

    void _quantify_rec(const NodePtr& n, map<NodePtr, V>& memo, V& result, bool exists) const {
        typename map<NodePtr, V>::iterator it;
        it = memo.find(n);
        if (it != memo.end())
            result = it->second;
        else if (n->index() == 0) {
            result = n->value();
            memo[n] = result;
        }
        else {
            V tmp;
            _quantify_rec(n->lo(), memo, tmp, exists);
            _quantify_rec(n->hi(), memo, result, exists);
            if (exists) {
                L::join(result, tmp, result);
            }
            else {
                L::meet(result, tmp, result);
            }
            if (F == SNF_MEET) {
                L::meet(result, n->value(), result);
            }
            else {
                L::join(result, n->value(), result);
            }
            memo[n] = result;
        }
    }

    static void _eval_rec(const vector<bool>& valuation, const NodePtr& node, V& result) {
        if (node->index() == 0) {
            node->value(result);
            return;
        }
        if (node->index() > (int) valuation.size()) {
            throw logic_error("eval : bad valuation");
        }
        NodePtr next = valuation[node->index() - 1] ? node->hi() : node->lo();
        if (F == SNF_JOIN) {
            V tmp;
            _eval_rec(valuation, next, tmp);
            L::join(tmp, node->value(), result);
        }
        else {
            V tmp;
            _eval_rec(valuation, next, tmp);
            L::meet(tmp, node->value(), result);
        }
    }

public:

    /*!
     *   \brief Default constructor. Is equal to L::bot().
     */
    LVBDD() {
        _root = _mkbot();
    }

    /*!
     *   \brief Destructor.
     */
    ~LVBDD() {
        // nothing to do (NodePtr's destructor is called automatically).
    }

    /*!
     *   \brief Copy constructor.
     */
    LVBDD(const LVBDD& other) {
        _root = other._root;
    }

    /*!
     *   \brief Assignment operator.
     */
    LVBDD& operator=(const LVBDD& other) {
        if (this != &other) {
            _root = other._root;
        }
        return *this;
    }

    /*!
     *   \brief Equality operator. This tests only for pointer equality of the roots
     *          and runs thus in O(1).
     */
    bool operator==(const LVBDD& other) const {
        return _root == other._root;
    }

    /*!
     *   \brief Inequality operator. Same remark as LVBDD::operator==.
     */
    bool operator!=(const LVBDD& other) const {
        return _root != other._root;
    }

    /*!
     *   \brief Creates a terminal LVBDD labeled by a chosen value.
     */
    static LVBDD terminal(const V& value) {
        LVBDD res;
        res._root = _mk(0, value, NULL, NULL);
        return res;
    }

    /*!
     *   \brief Creates a three-node LVBDD representing a positive or negative literal.
     *   \param [in] index : The proposition index of the literal. Must be strictly positive.
     *   \param [in] positive : The polarity of the literal (positive or negative).
     */
    static LVBDD literal(int index, bool positive) {
        if (index <= 0)
            throw invalid_argument("bad index");
        LVBDD res;
        NodePtr lo = positive ? _mkbot() : _mktop();
        NodePtr hi = positive ? _mktop() : _mkbot();
        V value;
        if (F == SNF_JOIN)
            L::bot(value);
        else
            L::top(value);
        res._root = _mk(index, value, lo, hi);
        return res;
    }

    /*!
     *  \brief Creates an LVBDD representing the "top" lattice constant.
     */
    static LVBDD top() {
        LVBDD result;
        V value;
        L::top(value);
        result._root = _mk(0, value, NULL, NULL);
        return result;
    }

    /*!
     *  \brief Creates an LVBDD representing the "bottom" lattice constant.
     */
    static LVBDD bot() {
        LVBDD result;
        V value;
        L::bot(value);
        result._root = _mk(0, value, NULL, NULL);
        return result;
    }

    /*!
     *  \brief Creates an *non-terminal* LVBDD "by hand".
     *         Do not use unless you know what you're doing.
     *         Wrong use of this fonction leads to normal form violation.
     *  \param [in] index : The proposition index of the root node.
     *                      Must be strictly positive.
     *  \param [in] value : The lattice value labeling the root node.
     *  \param [in] lo    : The node's low-child.
     *  \param [in] hi    : The node's high-child.
     */
    static LVBDD _mk_unsafe(int index, V value, LVBDD lo, LVBDD hi) {
        if (index < 1)
            throw logic_error("bad index");
        if (lo == hi)
            throw logic_error("lo and hi are the same");
        LVBDD result;
        result._root = _mk(index, value, lo._root, hi._root);
        return result;
    }

    /*!
     *  \brief Returns the number of orphaned nodes (reference count = 0) in the hash map.
     */
    static int orphan_count() {
        int count = 0;
        for (typename list<Node*>::iterator it = _orphans.begin(); it != _orphans.end(); it++) {
            if ((*it)->refcount() == 0)
                count++;
        }
        return count;
    }

    /*!
     *  \brief Prints a complete dump of the hash map. Use for debugging.
     */
    static void print_hashmap_debug(ostream& out = cout) {
        out << endl;
        typename HashMap::iterator it;
        for (it = _hash_map.begin(); it != _hash_map.end(); it++) {
            Node* n = (*it).second;
            out << "@" << n << "  index=" << n->index();
            out << "  lo@" << setw(10) << n->lo();
            out << "  hi@" << setw(10) << n->hi();
            out << "  rc=" << n->refcount();
            out << "  value=" << n->value();
            out << "  hash=" << n->hash();
            out << endl;
        }
        print_hashmap_stats();
        out << endl;
    }

    /*!
     *  \brief Prints various satistics about the hash map.
     */
    static void print_hashmap_stats(ostream& out = cout) {
        set<size_t, less<size_t> > hashes;
        map<int, int, less<int> > hash_count;
        for (typename HashMap::iterator it = _hash_map.begin(); it != _hash_map.end(); it++) {
            hashes.insert(it->second->hash());
            if (hash_count.find(it->second->hash()) == hash_count.end()) {
                hash_count[it->second->hash()] = 1;
            }
            else {
                hash_count[it->second->hash()] += 1;
            }
        }
        int max = -1;
        int max_hash = -1;
        map<int, int, less<int> >::const_iterator itt;
        for (itt = hash_count.begin(); itt != hash_count.end(); itt++) {
            if (itt->second > max) {
                max = itt->second;
                max_hash = itt->first;
            }
        }
        float lf = float(_hash_map.size()) / _hash_map.bucket_count();
        out << "HashMap size           = " << _hash_map.size() << endl;
        out << "HashMap load factor    = " << lf << endl;
        out << "HashMap bucket count   = " << _hash_map.bucket_count() << endl;
        out << "# of different hashes  = " << hashes.size() << endl;
        out << "Most usage of hash key = " << max << " (" << max_hash << ")";
        out << endl;
        out << "# of calls to MK       = " << _mk_call_count << endl;
        float pct = 100.0 - ((100.0 * (_orphan_fetches + _node_allocations)) / _mk_call_count);
        out << "# of MK cache hits     = " << _mk_call_count - _orphan_fetches - _node_allocations;
        out << " (" << pct << "%)" << endl;
        pct = (100.0 * _orphan_fetches) / _mk_call_count;
        out << "# of orphans reused    = " << _orphan_fetches << " (";
        out << pct << "%)" << endl;
        pct = (100.0 * _node_allocations) / _mk_call_count;
        out << "# of node allocations  = " << _node_allocations << " (";
        out << pct << "%)" << endl;
        out << "# of orphan nodes      = " << orphan_count() << endl;
    }

    /*!
     *   \brief Deallocates all memory used by orphaned nodes.
     */
    static void free_orphaned_nodes() {
        typename list<Node*>::iterator it;
        Node* n;
        for (it = _orphans.begin(); it != _orphans.end(); it++) {
            if ((*it)->refcount())
                continue;
            n = *it;
            delete n;
        }
        _orphans.clear();
    }

    /*!
     *   \brief Clears the memoization tables of join and meet operations.
     */
    static void clear_memoization_tables() {
        _memo_quad.clear();
        _memo_nn.clear();
        _memo_nv.clear();
    }

    /*!
     *   \brief Returns the number of entries in the memoization tables.
     */
    static int memoization_tables_size() {
        return _memo_quad.size() + _memo_nn.size() + _memo_nv.size();
    }

    /*!
     *  \brief Returns the total number of allocated nodes (orphans included).
     *         Always equal to orphan_count() + node_count().
     */
    static size_t total_node_count() {
        return _hash_map.size();
    }

    /*!
     *  \brief Returns the number of nodes in use (non orphaned) in the hash map.
     */
    static size_t node_count() {
        return total_node_count() - orphan_count();
    }

    /*!
     *  \brief Initialises the hash map.
     *  You MUST call this function before calling any other in the package.
     */
    static void init() {
        Node *dk = new Node();
        *dk = Node::deleted_key();
        _hash_map.set_deleted_key(dk);
#ifdef USING_DENSE_HASH_MAP
        Node *ek = new Node();
        *ek = Node::empty_key();
        _hash_map.set_empty_key(ek);
#endif
    }

    /*!
     *   \brief Computes the join of two LVBDD.
     *          Runtime complexity depends on the chosen normal form.
     */
    LVBDD join(const LVBDD& other) const {
        LVBDD result;
        if (F == SNF_MEET) {
            V value;
            L::top(value);
            result._root = _snf_hard_rec(_root, other._root, value, value);
        }
        else if (F == SNF_JOIN) {
            result._root = _snf_easy_rec(_root, other._root);
        }
        else { // UNF
            result._root = _unf_apply_rec(_root, other._root, false);
        }
        return result;
    }

    /*!
     *   \brief Syntactic sugar for LVBDD::join.
     */
    LVBDD operator|(const LVBDD& other) const {
        return join(other);
    }

    /*!
     *   \brief Computes the meet of two LVBDD.
     *          Runtime complexity depends on the chosen normal form.
     */
    LVBDD meet(const LVBDD& other) const {
        LVBDD result;
        if (F == SNF_JOIN) {
            V value;
            L::bot(value);
            result._root = _snf_hard_rec(_root, other._root, value, value);
        }
        else if (F == SNF_MEET) {
            result._root = _snf_easy_rec(_root, other._root);
        }
        else { // UNF
            result._root = _unf_apply_rec(_root, other._root, true);
        }
        return result;
    }

    /*!
     *   \brief Syntactic sugar for LVBDD::meet.
     */
    LVBDD operator&(const LVBDD& other) const {
        return meet(other);
    }

    /*!
     *   \brief Returns the number of nodes in the LVBDD.
     */
    int size() const {
        set<Node*> my_nodes;
        _root->_descendants(my_nodes);
        return my_nodes.size();
    }

    /*!
     *   \brief Returns the number of non-terminal nodes in the LVBDD.
     */
    int nonterminal_size() const {
        set<Node*> my_nodes;
        _root->_descendants(my_nodes);
        typename set<Node*>::iterator it;
        int count = 0;
        for (it = my_nodes.begin(); it != my_nodes.end(); it++) {
            if ((*it)->index() != 0) {
                count++;
            }
        }
        return count;
    }

    /*!
     *   \brief Returns the number of terminal nodes in the LVBDD.
     */
    int terminal_size() const {
        set<Node*> my_nodes;
        _root->_descendants(my_nodes);
        typename set<Node*>::iterator it;
        int count = 0;
        for (it = my_nodes.begin(); it != my_nodes.end(); it++) {
            if ((*it)->index() == 0) {
                count++;
            }
        }
        return count;
    }

    /*!
     *  \brief Returns Lattice<V>::list_size() on the list of all values
     *         labeling the LVBDD.
     */
    int values_size() const {
        set<Node*> my_nodes;
        _root->_descendants(my_nodes);
        list<V> l;
        typename set<Node*>::iterator it;
        for (it = my_nodes.begin(); it != my_nodes.end(); it++) {
            V val = (*it)->value();
            l.push_back(val);
        }
        return L::list_size(l);
    }

    /*!
     *   \brief Prints a DOT diagram of the LVBDD.
     */
    void print_graphviz(ostream& out = cout) const {
        Node::print_graphviz(out, _root);
    }

    /*!
     *   \brief Prints a DOT diagram of the entire hash map. Use for debugging.
     */
    static void print_graphviz_debug(ostream& out = cout) {
        Node::print_graphviz(out, NULL);
    }

    /*!
     *   \brief Quantifies all decision variables existentially.
     */
    void exists(V& result) const {
        if (F == SNF_MEET) {
            result = _root->value();
        }
        else {
            map<NodePtr, V> memo;
            _quantify_rec(_root, memo, result, true);
        }
    }

    /*!
     *   \brief Quantifies all decision variables existentially.
     */
    V exists() const {
        V result;
        exists(result);
        return result;
    }

    /*!
     *   \brief Quantifies all decision variables universally.
     */
    void forall(V& result) const {
        if (F == SNF_JOIN) {
            result = _root->value();
        }
        else {
            map<NodePtr, V> memo;
            _quantify_rec(_root, memo, result, false);
        }
    }

    /*!
     *   \brief Quantifies all decision variables universally.
     */
    V forall() const {
        V result;
        forall(result);
        return result;
    }

    /*!
     *   \brief Returns the low-child of the LVBDD. T
     *          Throws std::logic_error on terminal LVBDD.
     */
    LVBDD lo() const {
        if (_root->index() == 0)
            throw logic_error("called lo() on terminal LVBDD");
        LVBDD result;
        result._root = _root->lo();
        return result;
    }

    /*!
     *   \brief Returns the high-child of the LVBDD. T
     *          Throws std::logic_error on terminal LVBDD.
     */
    LVBDD hi() const {
        if (_root->index() == 0)
            throw logic_error("called hi() on terminal LVBDD");
        LVBDD result;
        result._root = _root->hi();
        return result;
    }

    /*!
     *   \brief Returns the index (strictly positive) of the node's proposition
     *          for nonterminal LVBDD; returns 0 for terminal LVBDD.
     */
    int index() const {
        return _root->index();
    }

    /*!
     *   \brief Returns the lattice value labeling the root.
     */
    void root_value(V& result) const {
        _root->value(result);
    }

    /*!
     *   \brief Returns the lattice value labeling the root.
     */
    V root_value() const {
        return _root->value();
    }

    /*!
     *   \brief Evaluate the LVBDD with a single propositional valuation.
     *          Throws std::logic_error if the valuation is incomplete.
     *   \param [in] valuation : The valuation. Proposition i is indexed at valuation[i-1].
     *   \param [out] result : The resulting lattice value.
     */
    void evaluate(const vector<bool>& valuation, V& result) const {
        result = _eval_rec(valuation, _root, result);
    }

    /*!
     *   \brief Evaluate the LVBDD with a single propositional valuation.
     *          Throws std::logic_error if the valuation is incomplete.
     *   \param [in] valuation : The valuation. Proposition i is indexed at valuation[i-1].
     *   \return The resulting lattice value.
     */
    V evaluate(const vector<bool>& valuation) const {
        V result;
        _eval_rec(valuation, _root, result);
        return result;
    }

    /*!
     *   \brief Prints the memory address of an LVBDD's root node on an output stream.
     */
    friend ostream& operator<<(ostream& out, const LVBDD& dd) {
        out << dd._root.get_ptr();
        return out;
    }

};

template<class L, int F>
typename LVBDD<L, F>::HashMap LVBDD<L, F>::_hash_map;

template<class L, int F>
list<Node<L, F>*> LVBDD<L, F>::_orphans;

template<class L, int F>
int LVBDD<L, F>::_mk_call_count;

template<class L, int F>
int LVBDD<L, F>::_orphan_fetches;

template<class L, int F>
int LVBDD<L, F>::_node_allocations;

template<class L, int F>
typename LVBDD<L, F>::NodeNodeMemoizer LVBDD<L, F>::_memo_nn;

template<class L, int F>
typename LVBDD<L, F>::NodeValueMemoizer LVBDD<L, F>::_memo_nv;

template<class L, int F>
typename LVBDD<L, F>::QuadrupleMemoizer LVBDD<L, F>::_memo_quad;

} // namespace lava


#endif /* __LAVABDD_H__ */