diff --git a/src/sylvan_bdd.c b/src/sylvan_bdd.c
index 5be6d039..79d09f8c 100644
--- a/src/sylvan_bdd.c
+++ b/src/sylvan_bdd.c
@@ -131,6 +131,106 @@ TASK_IMPL_3(BDD, sylvan_and, BDD, a, BDD, b, BDDVAR, prev_level)
     return result;
 }
 
+/*
+    sylvan_disjoint could be implemented as "sylvan_and(a,b)==sylvan_false",
+    but this implementation avoids building new nodes and allows more short-circuitry.
+*/
+TASK_IMPL_3(char, sylvan_disjoint, BDD, a, BDD, b, BDDVAR, prev_level)
+{
+    /* Terminal cases */
+    if (a == sylvan_false || b == sylvan_false) return 1; 
+    if (a == sylvan_true || b == sylvan_true) return 0; /* since a,b != sylvan_false */
+    if (a == b) return 0; /* since a,b != sylvan_false */
+    if (a == BDD_TOGGLEMARK(b)) return 1;
+
+    sylvan_gc_test();
+
+    /* Count operation */
+    sylvan_stats_count(BDD_DISJOINT);
+
+    /* Improve for caching */
+    if (BDD_STRIPMARK(a) > BDD_STRIPMARK(b)) {
+        BDD t = b;
+        b = a;
+        a = t;
+    }
+
+    bddnode_t na = MTBDD_GETNODE(a);
+    bddnode_t nb = MTBDD_GETNODE(b);
+
+    BDDVAR va = bddnode_getvariable(na);
+    BDDVAR vb = bddnode_getvariable(nb);
+    BDDVAR level = va < vb ? va : vb;
+
+    int cachenow = granularity < 2 || prev_level == 0 ? 1 : prev_level / granularity != level / granularity;
+    if (cachenow) {
+        BDD result;
+        if (cache_get3(CACHE_BDD_DISJOINT, a, b, sylvan_false, &result)) {
+            sylvan_stats_count(BDD_DISJOINT_CACHED);
+            return (result==sylvan_false ? 0 : 1);
+        }
+    }
+
+    // Get cofactors
+    BDD aLow = a, aHigh = a;
+    BDD bLow = b, bHigh = b;
+    if (level == va) {
+        aLow = node_low(a, na);
+        aHigh = node_high(a, na);
+    }
+    if (level == vb) {
+        bLow = node_low(b, nb);
+        bHigh = node_high(b, nb);
+    }
+
+    int low=-1, high=-1, result;
+
+    // Try to obtain the subresults without recursion (short-circuiting)
+
+    if (aHigh == sylvan_false || bHigh == sylvan_false) {
+        high = 1;
+    } else if (aHigh == sylvan_true || bHigh == sylvan_true) {
+        high = 0; /* since none of them is sylvan_false */
+    } else if (aHigh == bHigh) {
+        high = 0; /* since none of them is sylvan_false */
+    } else if (aHigh == BDD_TOGGLEMARK(bHigh)) {
+        high = 1;
+    }
+
+    if (aLow == sylvan_false || bLow == sylvan_false) {
+        low = 1;
+    } else if (aLow == sylvan_true || bLow == sylvan_true) {
+        low = 0; /* since none of them is sylvan_false */
+    } else if (aLow == bLow) {
+        low = 0; /* since none of them is sylvan_false */
+    } else if (aLow == BDD_TOGGLEMARK(bLow)) {
+        low = 1;
+    }
+     
+    // Compute the result, if necessary, by parallel recursion
+
+    if (high==0 || low==0) {
+        result = 0;
+    }
+    else {
+        if (high==-1) SPAWN(sylvan_disjoint, aHigh, bHigh, level);
+        if (low ==-1) low = CALL(sylvan_disjoint, aLow, bLow, level);
+        if (high==-1) high = SYNC(sylvan_disjoint);
+        result = high && low;
+    }
+
+    // Store result in the cache and then return
+
+    if (cachenow) {
+        BDD to_cache = (result ? sylvan_true : sylvan_false);
+        if (cache_put3(CACHE_BDD_DISJOINT, a, b, sylvan_false, to_cache)) {
+            sylvan_stats_count(BDD_DISJOINT_CACHEDPUT);
+        }
+    }
+
+    return result;
+}
+
 TASK_IMPL_3(BDD, sylvan_xor, BDD, a, BDD, b, BDDVAR, prev_level)
 {
     /* Terminal cases */
diff --git a/src/sylvan_bdd.h b/src/sylvan_bdd.h
index 36bf8d0e..2f00fff5 100644
--- a/src/sylvan_bdd.h
+++ b/src/sylvan_bdd.h
@@ -78,6 +78,10 @@ TASK_DECL_3(BDD, sylvan_xor, BDD, BDD, BDDVAR);
 #define sylvan_diff(a,b) sylvan_and(a,sylvan_not(b))
 #define sylvan_less(a,b) sylvan_and(sylvan_not(a),b)
 
+TASK_DECL_3(char, sylvan_disjoint, BDD, BDD, BDDVAR);
+#define sylvan_disjoint(a,b) (RUN(sylvan_disjoint,a,b,0))
+#define sylvan_subset(a,b) (RUN(sylvan_disjoint,a,sylvan_not(b),0))
+
 /* Create a BDD representing just <var> or the negation of <var> */
 static inline BDD
 sylvan_nithvar(uint32_t var)
diff --git a/src/sylvan_int.h b/src/sylvan_int.h
index 84cc6ba9..a1609dbf 100644
--- a/src/sylvan_int.h
+++ b/src/sylvan_int.h
@@ -68,6 +68,7 @@ static const uint64_t CACHE_BDD_CLOSURE             = (13LL<<40);
 static const uint64_t CACHE_BDD_ISBDD               = (14LL<<40);
 static const uint64_t CACHE_BDD_SUPPORT             = (15LL<<40);
 static const uint64_t CACHE_BDD_PATHCOUNT           = (16LL<<40);
+static const uint64_t CACHE_BDD_DISJOINT            = (17LL<<40);
 
 // MDD operations
 static const uint64_t CACHE_MDD_RELPROD             = (20LL<<40);
diff --git a/src/sylvan_obj.cpp b/src/sylvan_obj.cpp
index 2acf0536..42d847d3 100644
--- a/src/sylvan_obj.cpp
+++ b/src/sylvan_obj.cpp
@@ -45,9 +45,7 @@ Bdd::operator=(const Bdd& right)
 bool
 Bdd::operator<=(const Bdd& other) const
 {
-    // TODO: better implementation, since we are not interested in the BDD result
-    BDD r = sylvan_ite(this->bdd, sylvan_not(other.bdd), sylvan_false);
-    return r == sylvan_false;
+    return sylvan_subset(this->bdd, other.bdd) == 1;
 }
 
 bool
@@ -219,12 +217,16 @@ Bdd::Xnor(const Bdd &g) const
     return sylvan_equiv(bdd, g.bdd);
 }
 
+bool
+Bdd::Disjoint(const Bdd &g) const
+{
+    return sylvan_disjoint(bdd, g.bdd) == 1;
+}
+
 bool
 Bdd::Leq(const Bdd &g) const
 {
-    // TODO: better implementation, since we are not interested in the BDD result
-    BDD r = sylvan_ite(bdd, sylvan_not(g.bdd), sylvan_false);
-    return r == sylvan_false;
+    return sylvan_subset(bdd, g.bdd) == 1;
 }
 
 Bdd
diff --git a/src/sylvan_obj.hpp b/src/sylvan_obj.hpp
index 6ddd7d08..4b1dc44e 100644
--- a/src/sylvan_obj.hpp
+++ b/src/sylvan_obj.hpp
@@ -186,6 +186,11 @@ class Bdd {
      */
     Bdd Xnor(const Bdd& g) const;
 
+    /**
+     * @brief Returns whether all f and g are disjoint, i.e. f * g is false
+     */
+    bool Disjoint(const Bdd& g) const;
+
     /**
      * @brief Returns whether all elements in f are also in g
      */
diff --git a/src/sylvan_stats.c b/src/sylvan_stats.c
index 31e7eb2b..456edf07 100644
--- a/src/sylvan_stats.c
+++ b/src/sylvan_stats.c
@@ -67,6 +67,7 @@ struct
     {2, BDD_SATCOUNT, "BDD satcount"},
     {2, BDD_PATHCOUNT, "BDD pathcount"},
     {2, BDD_ISBDD, "BDD isbdd"},
+    {2, BDD_DISJOINT, "BDD disjoint"},
 
     {2, MTBDD_APPLY, "MTBDD binary apply"},
     {2, MTBDD_UAPPLY, "MTBDD unary apply"},
diff --git a/src/sylvan_stats.h b/src/sylvan_stats.h
index 1715cde9..719afedf 100644
--- a/src/sylvan_stats.h
+++ b/src/sylvan_stats.h
@@ -53,6 +53,7 @@ typedef enum {
     OPCOUNTER(BDD_ISBDD),
     OPCOUNTER(BDD_SUPPORT),
     OPCOUNTER(BDD_PATHCOUNT),
+    OPCOUNTER(BDD_DISJOINT),
 
     /* MTBDD operations */
     OPCOUNTER(MTBDD_APPLY),
diff --git a/test/test_basic.c b/test/test_basic.c
index aafbda71..dc432c81 100644
--- a/test/test_basic.c
+++ b/test/test_basic.c
@@ -321,6 +321,39 @@ test_operators()
     return 0;
 }
 
+static int
+test_disjoint_subset()
+{
+    // We need to test: disjoint, subset
+    
+    int vars=3;
+    BDD v[vars];
+    for (int i=0; i<vars; i++) v[i] = sylvan_nithvar(i);
+
+    BDD test_input[] = {
+        sylvan_true, sylvan_false,
+        sylvan_false, sylvan_true,
+        v[0], v[1],
+        v[1], v[1],
+        v[0], sylvan_not(v[0]),
+        sylvan_and(v[0],v[1]), v[2],
+        sylvan_and(v[0],v[1]), sylvan_and(sylvan_not(v[0]),v[1]),
+        sylvan_and(v[0],v[1]), sylvan_or(sylvan_not(v[0]),v[2]),
+        sylvan_and(v[0],v[1]), sylvan_and(v[0],sylvan_not(v[1])),
+        sylvan_or(v[0],v[1]), sylvan_and(v[0],sylvan_not(v[1])),
+        sylvan_and(v[1],sylvan_or(v[0],v[2])), sylvan_or(v[1],sylvan_and(v[0],sylvan_not(v[2])))
+    };
+
+    for (int i=0; i<11; i++) {
+        BDD t1 = test_input[2*i];
+        BDD t2 = test_input[2*i+1];
+        test_assert(sylvan_disjoint(t1,t2) == (sylvan_and(t1,t2)==sylvan_false));
+        test_assert(sylvan_subset(t1,t2) == (sylvan_or(sylvan_not(t1),t2) == sylvan_true));
+    }
+
+    return 0;
+}
+
 int
 test_relprod()
 {
@@ -526,6 +559,8 @@ TASK_0(int, runtests)
     for (int j=0;j<10;j++) if (test_compose()) return 1;
     printf("Testing operators.\n");
     for (int j=0;j<10;j++) if (test_operators()) return 1;
+    printf("Testing disjoint and subset.\n");
+    for (int j=0;j<10;j++) if (test_disjoint_subset()) return 1;
 
     printf("Testing ldd.\n");
     if (test_ldd()) return 1;