diff --git a/Cargo.lock b/Cargo.lock
index 18df324..fc27e76 100644
--- a/Cargo.lock
+++ b/Cargo.lock
@@ -4,7 +4,7 @@ version = 3
 
 [[package]]
 name = "adf_bdd"
-version = "0.2.1"
+version = "0.2.2"
 dependencies = [
  "biodivine-lib-bdd",
  "derivative",
diff --git a/bin/src/main.rs b/bin/src/main.rs
index 3ba4acf..1dbcf20 100644
--- a/bin/src/main.rs
+++ b/bin/src/main.rs
@@ -40,6 +40,7 @@ OPTIONS:
         --rust_log <RUST_LOG>     Sets the verbosity to 'warn', 'info', 'debug' or 'trace' if -v and
                                   -q are not use [env: RUST_LOG=debug]
         --stm                     Compute the stable models
+        --stmc                    Compute the stable models with the help of modelcounting
         --stmpre                  Compute the stable models with a pre-filter (only hybrid lib-mode)
         --stmrew                  Compute the stable models with a single-formula rewriting (only
                                   hybrid lib-mode)
@@ -104,6 +105,9 @@ struct App {
     /// Compute the stable models
     #[clap(long = "stm")]
     stable: bool,
+    /// Compute the stable models with the help of modelcounting
+    #[clap(long = "stmc")]
+    stable_counting: bool,
     /// Compute the stable models with a pre-filter (only hybrid lib-mode)
     #[clap(long = "stmpre")]
     stable_pre: bool,
@@ -201,6 +205,7 @@ impl App {
                 }
 
                 let printer = naive_adf.print_dictionary();
+
                 if self.complete {
                     for model in naive_adf.complete() {
                         print!("{}", printer.print_interpretation(&model));
@@ -213,6 +218,12 @@ impl App {
                     }
                 }
 
+                if self.stable_counting {
+                    for model in naive_adf.stable_count_optimisation() {
+                        print!("{}", printer.print_interpretation(&model));
+                    }
+                }
+
                 if self.stable_pre {
                     for model in naive_adf.stable_with_prefilter() {
                         print!("{}", printer.print_interpretation(&model));
diff --git a/lib/Cargo.toml b/lib/Cargo.toml
index 4a8365a..2c4c86c 100644
--- a/lib/Cargo.toml
+++ b/lib/Cargo.toml
@@ -1,6 +1,6 @@
 [package]
 name = "adf_bdd"
-version = "0.2.1"
+version = "0.2.2"
 authors = ["Stefan Ellmauthaler <stefan.ellmauthaler@tu-dresden.de>"]
 edition = "2021"
 repository = "https://github.com/ellmau/adf-obdd/"
diff --git a/lib/src/adf.rs b/lib/src/adf.rs
index 195d7c5..dac5c57 100644
--- a/lib/src/adf.rs
+++ b/lib/src/adf.rs
@@ -1,7 +1,9 @@
-//! This module describes the abstract dialectical framework
-//!
-//!  - computing interpretations
-//!  - computing fixpoints
+/*!
+This module describes the abstract dialectical framework
+
+ - computing interpretations
+ - computing fixpoints
+*/
 
 use serde::{Deserialize, Serialize};
 
@@ -18,7 +20,7 @@ use crate::{
 };
 
 #[derive(Serialize, Deserialize, Debug)]
-/// Representation of an ADF, with an ordering and dictionary of statement <-> number relations, a binary decision diagram, and a list of acceptance functions in Term representation
+/// Representation of an ADF, with an ordering and dictionary of statement <-> number relations, a binary decision diagram, and a list of acceptance functions in Term representation.
 ///
 /// Please note that due to the nature of the underlying reduced and ordered Bdd the concept of a [`Term`][crate::datatypes::Term] represents one (sub) formula as well as truth-values.
 pub struct Adf {
@@ -226,8 +228,8 @@ impl Adf {
         new_interpretation
     }
 
-    /// Computes the stable models
-    /// Returns an Iterator which contains all stable models
+    /// Computes the stable models.
+    /// Returns an Iterator which contains all stable models.
     pub fn stable<'a, 'c>(&'a mut self) -> impl Iterator<Item = Vec<Term>> + 'c
     where
         'a: 'c,
@@ -264,9 +266,9 @@ impl Adf {
             .map(|(int, _grd)| int)
     }
 
-    /// Computes the stable models
-    /// Returns a vector with all stable models, using a single-formula representation in biodivine to enumerate the possible models
-    /// Note that the biodivine adf needs to be the one which instantiated the adf (if applicable)
+    /// Computes the stable models.
+    /// Returns a vector with all stable models, using a single-formula representation in biodivine to enumerate the possible models.
+    /// Note that the biodivine adf needs to be the one which instantiated the adf (if applicable).
     pub fn stable_bdd_representation(
         &mut self,
         biodivine: &crate::adfbiodivine::Adf,
@@ -294,8 +296,8 @@ impl Adf {
             .collect::<Vec<Vec<Term>>>()
     }
 
-    /// Computes the stable models
-    /// Returns an Iterator which contains all stable models
+    /// Computes the stable models.
+    /// Returns an Iterator which contains all stable models.
     pub fn stable_with_prefilter<'a, 'c>(&'a mut self) -> impl Iterator<Item = Vec<Term>> + 'c
     where
         'a: 'c,
@@ -347,6 +349,137 @@ impl Adf {
             .map(|(int, _grd)| int)
     }
 
+    /// Computes the stable models.
+    /// Returns an iterator which contains all stable models.
+    /// This variant uses the computation of model and counter-model counts.
+    pub fn stable_count_optimisation<'a, 'c>(&'a mut self) -> impl Iterator<Item = Vec<Term>> + 'c
+    where
+        'a: 'c,
+    {
+        log::debug!("[Start] stable count optimisation");
+        let grounded = self.grounded();
+        self.two_val_model_counts(&grounded)
+            .into_iter()
+            .filter(|int| self.stability_check(int))
+    }
+
+    fn stability_check(&mut self, interpretation: &[Term]) -> bool {
+        let mut new_int = self.ac.clone();
+        for ac in new_int.iter_mut() {
+            *ac = interpretation
+                .iter()
+                .enumerate()
+                .fold(*ac, |acc, (var, term)| {
+                    if term.is_truth_value() && !term.is_true() {
+                        self.bdd.restrict(acc, Var(var), false)
+                    } else {
+                        acc
+                    }
+                });
+        }
+
+        let grd = self.grounded_internal(&new_int);
+        for (idx, grd) in grd.iter().enumerate() {
+            if !grd.compare_inf(&interpretation[idx]) {
+                return false;
+            }
+        }
+        true
+    }
+
+    fn two_val_model_counts(&mut self, interpr: &[Term]) -> Vec<Vec<Term>> {
+        log::trace!("two_val_model_counts({:?}) called ", interpr);
+        if let Some((idx, ac)) = interpr
+            .iter()
+            .enumerate()
+            .filter(|(_idx, val)| !val.is_truth_value())
+            .min_by(|(_idx_a, val_a), (_idx_b, val_b)| {
+                self.bdd
+                    .models(**val_a, true)
+                    .minimum()
+                    .cmp(&self.bdd.models(**val_b, true).minimum())
+            })
+        {
+            let mut result = Vec::new();
+            let check_models = !self.bdd.models(*ac, true).more_models();
+            log::trace!(
+                "Identified Var({}) with ac {:?} to be {}",
+                idx,
+                ac,
+                check_models
+            );
+            let _ = self // return value can be ignored, but must be catched
+                .bdd
+                .interpretations(*ac, check_models, Var(idx), &[], &[])
+                .iter()
+                .try_for_each(|(negative, positive)| {
+                    let mut new_int = interpr.to_vec();
+                    let res = negative
+                        .iter()
+                        .try_for_each(|var| {
+                            if new_int[var.value()].is_true() {
+                                return Err(());
+                            }
+                            new_int[var.value()] = Term::BOT;
+                            Ok(())
+                        })
+                        .and(positive.iter().try_for_each(|var| {
+                            if new_int[var.value()].is_truth_value()
+                                && !new_int[var.value()].is_true()
+                            {
+                                return Err(());
+                            }
+                            new_int[var.value()] = Term::TOP;
+                            Ok(())
+                        }));
+                    if res.is_ok() {
+                        new_int[idx] = if check_models { Term::TOP } else { Term::BOT };
+                        let upd_int = self.update_interpretation(&new_int);
+                        result.append(&mut self.two_val_model_counts(&upd_int));
+                    }
+                    res
+                });
+            // checked one alternative, we can now conclude that only the other option may work
+            log::trace!("checked one alternative, concluding the other value");
+            let new_int = interpr
+                .iter()
+                .map(|tree| self.bdd.restrict(*tree, Var(idx), !check_models))
+                .collect::<Vec<Term>>();
+            let mut upd_int = self.update_interpretation(&new_int);
+
+            // TODO: should be "must be true/false" instead of setting it to TOP/BOT and will need sanity checks at every iteration
+            log::trace!("\nnew_int {new_int:?}\nupd_int {upd_int:?}");
+            if new_int[idx].no_inf_decrease(&upd_int[idx]) {
+                upd_int[idx] = if check_models { Term::BOT } else { Term::TOP };
+                if new_int[idx].no_inf_decrease(&upd_int[idx]) {
+                    result.append(&mut self.two_val_model_counts(&upd_int));
+                }
+            }
+            result
+        } else {
+            // filter has created empty iterator
+            vec![interpr.to_vec()]
+        }
+    }
+
+    fn update_interpretation(&mut self, interpretation: &[Term]) -> Vec<Term> {
+        interpretation
+            .iter()
+            .map(|ac| {
+                interpretation
+                    .iter()
+                    .enumerate()
+                    .fold(*ac, |acc, (idx, val)| {
+                        if val.is_truth_value() {
+                            self.bdd.restrict(acc, Var(idx), val.is_true())
+                        } else {
+                            acc
+                        }
+                    })
+            })
+            .collect::<Vec<Term>>()
+    }
+
     /// Computes the complete models
     /// Returns an Iterator which contains all complete models
     pub fn complete<'a, 'c>(&'a mut self) -> impl Iterator<Item = Vec<Term>> + 'c
@@ -414,11 +547,11 @@ impl Adf {
 
                 let n_vdps = { |t| self.bdd.var_dependencies(t).len() };
 
-                let fc = match mcs.1 > 2 {
+                let fc = match mcs.models > 2 {
                     true => 2 * n_vdps(*t),
                     _ => 0,
                 };
-                let cfc = match mcs.0 > 2 {
+                let cfc = match mcs.cmodels > 2 {
                     true => 2 * n_vdps(*t),
                     _ => 0,
                 };
@@ -561,6 +694,50 @@ mod test {
         assert_eq!(adf.stable().next(), None);
     }
 
+    #[test]
+    fn stable_w_counts() {
+        let parser = AdfParser::default();
+        parser.parse()("s(a).s(b).s(c).s(d).ac(a,c(v)).ac(b,b).ac(c,and(a,b)).ac(d,neg(b)).\ns(e).ac(e,and(b,or(neg(b),c(f)))).s(f).\n\nac(f,xor(a,e)).")
+            .unwrap();
+        let mut adf = Adf::from_parser(&parser);
+
+        let mut stable = adf.stable_count_optimisation();
+        assert_eq!(
+            stable.next(),
+            Some(vec![
+                Term::TOP,
+                Term::BOT,
+                Term::BOT,
+                Term::TOP,
+                Term::BOT,
+                Term::TOP
+            ])
+        );
+        assert_eq!(stable.next(), None);
+        let parser = AdfParser::default();
+        parser.parse()("s(a).s(b).ac(a,neg(b)).ac(b,neg(a)).").unwrap();
+        let mut adf = Adf::from_parser(&parser);
+        let mut stable = adf.stable_count_optimisation();
+
+        assert_eq!(stable.next(), Some(vec![Term::BOT, Term::TOP]));
+        assert_eq!(stable.next(), Some(vec![Term::TOP, Term::BOT]));
+        assert_eq!(stable.next(), None);
+
+        let parser = AdfParser::default();
+        parser.parse()("s(a).s(b).ac(a,b).ac(b,a).").unwrap();
+        let mut adf = Adf::from_parser(&parser);
+
+        assert_eq!(
+            adf.stable_count_optimisation().collect::<Vec<_>>(),
+            vec![vec![Term::BOT, Term::BOT]]
+        );
+
+        let parser = AdfParser::default();
+        parser.parse()("s(a).s(b).ac(a,neg(a)).ac(b,a).").unwrap();
+        let mut adf = Adf::from_parser(&parser);
+        assert_eq!(adf.stable_count_optimisation().next(), None);
+    }
+
     #[test]
     fn complete() {
         let parser = AdfParser::default();
diff --git a/lib/src/datatypes/bdd.rs b/lib/src/datatypes/bdd.rs
index 2f1dedb..22cf282 100644
--- a/lib/src/datatypes/bdd.rs
+++ b/lib/src/datatypes/bdd.rs
@@ -72,6 +72,14 @@ impl Term {
         self.is_truth_value() == other.is_truth_value() && self.is_true() == other.is_true()
     }
 
+    /// Returns true if the information of *other* does not decrease and it is not inconsistent.
+    pub fn no_inf_decrease(&self, other: &Self) -> bool {
+        if self.compare_inf(other) {
+            return true;
+        }
+        !self.is_truth_value()
+    }
+
     /// Returns true, if the Term and the BDD have the same information-value
     pub fn cmp_information(&self, other: &biodivine_lib_bdd::Bdd) -> bool {
         self.is_truth_value() == other.is_truth_value() && self.is_true() == other.is_true()
@@ -178,7 +186,45 @@ impl BddNode {
 }
 
 /// Type alias for the pair of counter-models and models
-pub type ModelCounts = (usize, usize);
+#[derive(Debug, Clone, Copy, Eq, PartialEq, PartialOrd, Ord)]
+pub struct ModelCounts {
+    /// Contains the number of counter-models
+    pub cmodels: usize,
+    /// Contains the number of models
+    pub models: usize,
+}
+
+impl ModelCounts {
+    /// Represents the top-node model-counts
+    pub fn top() -> ModelCounts {
+        (0, 1).into()
+    }
+
+    /// Represents the bot-node model-counts
+    pub fn bot() -> ModelCounts {
+        (1, 0).into()
+    }
+
+    /// Returns the smaller size (models or counter-models).
+    pub fn minimum(&self) -> usize {
+        self.models.min(self.cmodels)
+    }
+
+    /// Returns true, if there are more models than counter-models.
+    /// If they are equal, the function returns true too.
+    pub fn more_models(&self) -> bool {
+        self.models >= self.minimum()
+    }
+}
+
+impl From<(usize, usize)> for ModelCounts {
+    fn from(tuple: (usize, usize)) -> Self {
+        ModelCounts {
+            cmodels: tuple.0,
+            models: tuple.1,
+        }
+    }
+}
 /// Type alias for the Modelcounts and the depth of a given Node in a BDD
 pub type CountNode = (ModelCounts, usize);
 /// Type alias for Facet counts, which contains number of facets and counter facets.
diff --git a/lib/src/obdd.rs b/lib/src/obdd.rs
index 0e2ac98..4cef5ab 100644
--- a/lib/src/obdd.rs
+++ b/lib/src/obdd.rs
@@ -51,11 +51,11 @@ impl Bdd {
             result
                 .count_cache
                 .borrow_mut()
-                .insert(Term::TOP, ((0, 1), 0));
+                .insert(Term::TOP, (ModelCounts::top(), 0));
             result
                 .count_cache
                 .borrow_mut()
-                .insert(Term::BOT, ((1, 0), 0));
+                .insert(Term::BOT, (ModelCounts::bot(), 0));
             result
         }
     }
@@ -102,6 +102,58 @@ impl Bdd {
         self.if_then_else(term_a, not_b, term_b)
     }
 
+    /// Computes the interpretations represented in the reduced BDD, which are either models or none.
+    /// *goal_var* is the variable to which the BDD is related to and it is ensured that the goal is consistent with the respective interpretation
+    /// *goal* is a boolean variable, which defines whether the models or inconsistent interpretations are of interest
+    pub fn interpretations(
+        &self,
+        tree: Term,
+        goal: bool,
+        goal_var: Var,
+        negative: &[Var],
+        positive: &[Var],
+    ) -> Vec<(Vec<Var>, Vec<Var>)> {
+        let mut result = Vec::new();
+        let node = self.nodes[tree.value()];
+        let var = node.var();
+        if tree.is_truth_value() {
+            return Vec::new();
+        }
+        // if the current var is the goal var, then only work with the hi-node if the goal is true
+        if (goal_var != var) || goal {
+            if node.hi().is_truth_value() {
+                if node.hi().is_true() == goal {
+                    result.push((negative.to_vec(), [positive, &[var]].concat()));
+                }
+            } else {
+                result.append(&mut self.interpretations(
+                    node.hi(),
+                    goal,
+                    goal_var,
+                    negative,
+                    &[positive, &[var]].concat(),
+                ));
+            }
+        }
+        // if the current var is the goal var, then only work with the lo-node if the goal is false
+        if (goal_var != var) || !goal {
+            if node.lo().is_truth_value() {
+                if node.lo().is_true() == goal {
+                    result.push(([negative, &[var]].concat(), positive.to_vec()));
+                }
+            } else {
+                result.append(&mut self.interpretations(
+                    node.lo(),
+                    goal,
+                    goal_var,
+                    &[negative, &[var]].concat(),
+                    positive,
+                ));
+            }
+        }
+        result
+    }
+
     pub fn restrict(&mut self, tree: Term, var: Var, val: bool) -> Term {
         let node = self.nodes[tree.0];
         #[cfg(feature = "variablelist")]
@@ -180,12 +232,20 @@ impl Bdd {
                     {
                         log::debug!("newterm: {} as {:?}", new_term, node);
                         let mut count_cache = self.count_cache.borrow_mut();
-                        let ((lo_cmodel, lo_model), lodepth) =
-                            *count_cache.get(&lo).expect("Cache corrupted");
-                        let ((hi_cmodel, hi_model), hidepth) =
-                            *count_cache.get(&hi).expect("Cache corrupted");
-                        log::debug!("lo (cm: {}, mo: {}, dp: {})", lo_cmodel, lo_model, lodepth);
-                        log::debug!("hi (cm: {}, mo: {}, dp: {})", hi_cmodel, hi_model, hidepth);
+                        let (lo_counts, lodepth) = *count_cache.get(&lo).expect("Cache corrupted");
+                        let (hi_counts, hidepth) = *count_cache.get(&hi).expect("Cache corrupted");
+                        log::debug!(
+                            "lo (cm: {}, mo: {}, dp: {})",
+                            lo_counts.cmodels,
+                            lo_counts.models,
+                            lodepth
+                        );
+                        log::debug!(
+                            "hi (cm: {}, mo: {}, dp: {})",
+                            hi_counts.cmodels,
+                            hi_counts.models,
+                            hidepth
+                        );
                         let (lo_exp, hi_exp) = if lodepth > hidepth {
                             (1, 2usize.pow((lodepth - hidepth) as u32))
                         } else {
@@ -196,9 +256,10 @@ impl Bdd {
                             new_term,
                             (
                                 (
-                                    lo_cmodel * lo_exp + hi_cmodel * hi_exp,
-                                    lo_model * lo_exp + hi_model * hi_exp,
-                                ),
+                                    lo_counts.cmodels * lo_exp + hi_counts.cmodels * hi_exp,
+                                    lo_counts.models * lo_exp + hi_counts.models * hi_exp,
+                                )
+                                    .into(),
                                 std::cmp::max(lodepth, hidepth) + 1,
                             ),
                         );
@@ -229,15 +290,15 @@ impl Bdd {
     /// Computes the number of counter-models, models, and variables for a given BDD-tree
     fn modelcount_naive(&self, term: Term) -> CountNode {
         if term == Term::TOP {
-            ((0, 1), 0)
+            (ModelCounts::top(), 0)
         } else if term == Term::BOT {
-            ((1, 0), 0)
+            (ModelCounts::bot(), 0)
         } else {
             let node = &self.nodes[term.0];
             let mut lo_exp = 0u32;
             let mut hi_exp = 0u32;
-            let ((lo_counter, lo_model), lodepth) = self.modelcount_naive(node.lo());
-            let ((hi_counter, hi_model), hidepth) = self.modelcount_naive(node.hi());
+            let (lo_counts, lodepth) = self.modelcount_naive(node.lo());
+            let (hi_counts, hidepth) = self.modelcount_naive(node.hi());
             if lodepth > hidepth {
                 hi_exp = (lodepth - hidepth) as u32;
             } else {
@@ -245,9 +306,10 @@ impl Bdd {
             }
             (
                 (
-                    lo_counter * 2usize.pow(lo_exp) + hi_counter * 2usize.pow(hi_exp),
-                    lo_model * 2usize.pow(lo_exp) + hi_model * 2usize.pow(hi_exp),
-                ),
+                    lo_counts.cmodels * 2usize.pow(lo_exp) + hi_counts.cmodels * 2usize.pow(hi_exp),
+                    lo_counts.models * 2usize.pow(lo_exp) + hi_counts.models * 2usize.pow(hi_exp),
+                )
+                    .into(),
                 std::cmp::max(lodepth, hidepth) + 1,
             )
         }
@@ -255,9 +317,9 @@ impl Bdd {
 
     fn modelcount_memoization(&self, term: Term) -> CountNode {
         if term == Term::TOP {
-            ((0, 1), 0)
+            (ModelCounts::top(), 0)
         } else if term == Term::BOT {
-            ((1, 0), 0)
+            (ModelCounts::bot(), 0)
         } else {
             if let Some(result) = self.count_cache.borrow().get(&term) {
                 return *result;
@@ -266,8 +328,8 @@ impl Bdd {
                 let node = &self.nodes[term.0];
                 let mut lo_exp = 0u32;
                 let mut hi_exp = 0u32;
-                let ((lo_counter, lo_model), lodepth) = self.modelcount_memoization(node.lo());
-                let ((hi_counter, hi_model), hidepth) = self.modelcount_memoization(node.hi());
+                let (lo_counts, lodepth) = self.modelcount_memoization(node.lo());
+                let (hi_counts, hidepth) = self.modelcount_memoization(node.hi());
                 if lodepth > hidepth {
                     hi_exp = (lodepth - hidepth) as u32;
                 } else {
@@ -275,9 +337,12 @@ impl Bdd {
                 }
                 (
                     (
-                        lo_counter * 2usize.pow(lo_exp) + hi_counter * 2usize.pow(hi_exp),
-                        lo_model * 2usize.pow(lo_exp) + hi_model * 2usize.pow(hi_exp),
-                    ),
+                        lo_counts.cmodels * 2usize.pow(lo_exp)
+                            + hi_counts.cmodels * 2usize.pow(hi_exp),
+                        lo_counts.models * 2usize.pow(lo_exp)
+                            + hi_counts.models * 2usize.pow(hi_exp),
+                    )
+                        .into(),
                     std::cmp::max(lodepth, hidepth) + 1,
                 )
             };
@@ -291,8 +356,12 @@ impl Bdd {
         self.generate_var_dependencies();
         #[cfg(feature = "adhoccounting")]
         {
-            self.count_cache.borrow_mut().insert(Term::TOP, ((0, 1), 0));
-            self.count_cache.borrow_mut().insert(Term::BOT, ((1, 0), 0));
+            self.count_cache
+                .borrow_mut()
+                .insert(Term::TOP, (ModelCounts::top(), 0));
+            self.count_cache
+                .borrow_mut()
+                .insert(Term::BOT, (ModelCounts::bot(), 0));
             for i in 0..self.nodes.len() {
                 log::debug!("fixing Term({})", i);
                 self.modelcount_memoization(Term(i));
@@ -462,7 +531,7 @@ mod test {
         let formula3 = bdd.xor(v1, v2);
         let formula4 = bdd.and(v3, formula2);
 
-        assert_eq!(bdd.models(v1, false), (1, 1));
+        assert_eq!(bdd.models(v1, false), (1, 1).into());
         let mut x = bdd.count_cache.get_mut().iter().collect::<Vec<_>>();
         x.sort();
         log::debug!("{:?}", formula1);
@@ -470,28 +539,28 @@ mod test {
             log::debug!("{:?}", x);
         }
         log::debug!("{:?}", x);
-        assert_eq!(bdd.models(formula1, false), (3, 1));
-        assert_eq!(bdd.models(formula2, false), (1, 3));
-        assert_eq!(bdd.models(formula3, false), (2, 2));
-        assert_eq!(bdd.models(formula4, false), (5, 3));
-        assert_eq!(bdd.models(Term::TOP, false), (0, 1));
-        assert_eq!(bdd.models(Term::BOT, false), (1, 0));
+        assert_eq!(bdd.models(formula1, false), (3, 1).into());
+        assert_eq!(bdd.models(formula2, false), (1, 3).into());
+        assert_eq!(bdd.models(formula3, false), (2, 2).into());
+        assert_eq!(bdd.models(formula4, false), (5, 3).into());
+        assert_eq!(bdd.models(Term::TOP, false), (0, 1).into());
+        assert_eq!(bdd.models(Term::BOT, false), (1, 0).into());
 
-        assert_eq!(bdd.models(v1, true), (1, 1));
-        assert_eq!(bdd.models(formula1, true), (3, 1));
-        assert_eq!(bdd.models(formula2, true), (1, 3));
-        assert_eq!(bdd.models(formula3, true), (2, 2));
-        assert_eq!(bdd.models(formula4, true), (5, 3));
-        assert_eq!(bdd.models(Term::TOP, true), (0, 1));
-        assert_eq!(bdd.models(Term::BOT, true), (1, 0));
+        assert_eq!(bdd.models(v1, true), (1, 1).into());
+        assert_eq!(bdd.models(formula1, true), (3, 1).into());
+        assert_eq!(bdd.models(formula2, true), (1, 3).into());
+        assert_eq!(bdd.models(formula3, true), (2, 2).into());
+        assert_eq!(bdd.models(formula4, true), (5, 3).into());
+        assert_eq!(bdd.models(Term::TOP, true), (0, 1).into());
+        assert_eq!(bdd.models(Term::BOT, true), (1, 0).into());
 
-        assert_eq!(bdd.modelcount_naive(v1), ((1, 1), 1));
-        assert_eq!(bdd.modelcount_naive(formula1), ((3, 1), 2));
-        assert_eq!(bdd.modelcount_naive(formula2), ((1, 3), 2));
-        assert_eq!(bdd.modelcount_naive(formula3), ((2, 2), 2));
-        assert_eq!(bdd.modelcount_naive(formula4), ((5, 3), 3));
-        assert_eq!(bdd.modelcount_naive(Term::TOP), ((0, 1), 0));
-        assert_eq!(bdd.modelcount_naive(Term::BOT), ((1, 0), 0));
+        assert_eq!(bdd.modelcount_naive(v1), ((1, 1).into(), 1));
+        assert_eq!(bdd.modelcount_naive(formula1), ((3, 1).into(), 2));
+        assert_eq!(bdd.modelcount_naive(formula2), ((1, 3).into(), 2));
+        assert_eq!(bdd.modelcount_naive(formula3), ((2, 2).into(), 2));
+        assert_eq!(bdd.modelcount_naive(formula4), ((5, 3).into(), 3));
+        assert_eq!(bdd.modelcount_naive(Term::TOP), ((0, 1).into(), 0));
+        assert_eq!(bdd.modelcount_naive(Term::BOT), ((1, 0).into(), 0));
 
         assert_eq!(
             bdd.modelcount_naive(formula4),
@@ -549,4 +618,58 @@ mod test {
                 assert!(left == right);
             });
     }
+
+    #[test]
+    fn interpretations() {
+        let mut bdd = Bdd::new();
+
+        let v1 = bdd.variable(Var(0));
+        let v2 = bdd.variable(Var(1));
+
+        let formula1 = bdd.and(v1, v2);
+        let formula2 = bdd.xor(v1, v2);
+
+        assert_eq!(
+            bdd.interpretations(formula1, true, Var(2), &[], &[]),
+            vec![(vec![], vec![Var(0), Var(1)])]
+        );
+
+        assert_eq!(
+            bdd.interpretations(formula1, true, Var(0), &[], &[]),
+            vec![(vec![], vec![Var(0), Var(1)])]
+        );
+
+        assert_eq!(
+            bdd.interpretations(formula1, false, Var(2), &[], &[]),
+            vec![(vec![Var(1)], vec![Var(0)]), (vec![Var(0)], vec![])]
+        );
+
+        assert_eq!(
+            bdd.interpretations(formula1, false, Var(0), &[], &[]),
+            vec![(vec![Var(0)], vec![])]
+        );
+
+        assert_eq!(
+            bdd.interpretations(formula2, false, Var(2), &[], &[]),
+            vec![
+                (vec![], vec![Var(0), Var(1)]),
+                (vec![Var(0), Var(1)], vec![])
+            ]
+        );
+
+        assert_eq!(
+            bdd.interpretations(formula2, true, Var(2), &[], &[]),
+            vec![(vec![Var(1)], vec![Var(0)]), (vec![Var(0)], vec![Var(1)])]
+        );
+
+        assert_eq!(
+            bdd.interpretations(formula2, true, Var(0), &[], &[]),
+            vec![(vec![Var(1)], vec![Var(0)])]
+        );
+
+        assert_eq!(
+            bdd.interpretations(Term::TOP, true, Var(0), &[], &[]),
+            vec![]
+        );
+    }
 }