diff --git a/src/builder/DdPrismModelBuilder.cpp b/src/builder/DdPrismModelBuilder.cpp
index e2c12baa2..993fca443 100644
--- a/src/builder/DdPrismModelBuilder.cpp
+++ b/src/builder/DdPrismModelBuilder.cpp
@@ -1175,6 +1175,15 @@ namespace storm {
                 STORM_LOG_TRACE("Iteration " << iteration << " of reachability analysis.");
                 changed = false;
                 storm::dd::Bdd<Type> tmp = reachableStates.relationalProduct(transitionBdd, generationInfo.rowMetaVariables, generationInfo.columnMetaVariables);
+                
+                storm::dd::Bdd<Type> tmp2 = reachableStates.andExists(transitionBdd, generationInfo.rowMetaVariables);
+                tmp2 = tmp2.swapVariables(generationInfo.rowColumnMetaVariablePairs);
+                
+                tmp2.exportToDot("tmp2.dot");
+                tmp.exportToDot("tmp.dot");
+
+                assert(tmp == tmp2);
+                
                 storm::dd::Bdd<Type> newReachableStates = tmp && (!reachableStates);
                                 
                 // Check whether new states were indeed discovered.
@@ -1183,6 +1192,9 @@ namespace storm {
                 }
                 
                 reachableStates |= newReachableStates;
+                
+                std::cout << "iter: " << iteration << " with nnz: " << reachableStates.getNonZeroCount() << std::endl;
+                
                 ++iteration;
             } while (changed);
             
diff --git a/src/storage/dd/Bdd.cpp b/src/storage/dd/Bdd.cpp
index ad6fe77f2..87234b0e1 100644
--- a/src/storage/dd/Bdd.cpp
+++ b/src/storage/dd/Bdd.cpp
@@ -156,7 +156,15 @@ namespace storm {
                 }
             }
             
-            return Bdd<LibraryType>(this->getDdManager(), internalBdd.relationalProduct(relation, rowVariables), newMetaVariables);
+            std::vector<InternalBdd<LibraryType>> columnVariables;
+            for (auto const& metaVariable : columnMetaVariables) {
+                DdMetaVariable<LibraryType> const& variable = this->getDdManager()->getMetaVariable(metaVariable);
+                for (auto const& ddVariable : variable.getDdVariables()) {
+                    columnVariables.push_back(ddVariable);
+                }
+            }
+            
+            return Bdd<LibraryType>(this->getDdManager(), internalBdd.relationalProduct(relation, rowVariables, columnVariables), newMetaVariables);
         }
         
         template<DdType LibraryType>
@@ -165,6 +173,14 @@ namespace storm {
             std::set<storm::expressions::Variable> newMetaVariables;
             std::set_difference(tmpMetaVariables.begin(), tmpMetaVariables.end(), rowMetaVariables.begin(), rowMetaVariables.end(), std::inserter(newMetaVariables, newMetaVariables.begin()));
             
+            std::vector<InternalBdd<LibraryType>> rowVariables;
+            for (auto const& metaVariable : rowMetaVariables) {
+                DdMetaVariable<LibraryType> const& variable = this->getDdManager()->getMetaVariable(metaVariable);
+                for (auto const& ddVariable : variable.getDdVariables()) {
+                    rowVariables.push_back(ddVariable);
+                }
+            }
+            
             std::vector<InternalBdd<LibraryType>> columnVariables;
             for (auto const& metaVariable : columnMetaVariables) {
                 DdMetaVariable<LibraryType> const& variable = this->getDdManager()->getMetaVariable(metaVariable);
@@ -173,7 +189,7 @@ namespace storm {
                 }
             }
             
-            return Bdd<LibraryType>(this->getDdManager(), internalBdd.inverseRelationalProduct(relation, columnVariables), newMetaVariables);
+            return Bdd<LibraryType>(this->getDdManager(), internalBdd.inverseRelationalProduct(relation, rowVariables, columnVariables), newMetaVariables);
         }
         
         template<DdType LibraryType>
diff --git a/src/storage/dd/cudd/InternalCuddBdd.cpp b/src/storage/dd/cudd/InternalCuddBdd.cpp
index 36d99d5cf..a257c9640 100644
--- a/src/storage/dd/cudd/InternalCuddBdd.cpp
+++ b/src/storage/dd/cudd/InternalCuddBdd.cpp
@@ -29,39 +29,25 @@ namespace storm {
             return !(*this == other);
         }
         
-        InternalBdd<DdType::CUDD> InternalBdd<DdType::CUDD>::relationalProduct(InternalBdd<DdType::CUDD> const& relation, std::vector<InternalBdd<DdType::CUDD>> const& rowVariables) const {
+        InternalBdd<DdType::CUDD> InternalBdd<DdType::CUDD>::relationalProduct(InternalBdd<DdType::CUDD> const& relation, std::vector<InternalBdd<DdType::CUDD>> const& rowVariables, std::vector<InternalBdd<DdType::CUDD>> const& columnVariables) const {
             InternalBdd<DdType::CUDD> cube = ddManager->getBddOne();
             for (auto const& variable : rowVariables) {
                 cube &= variable;
             }
             
             InternalBdd<DdType::CUDD> result = this->andExists(relation, cube);
-            
-            // Create the corresponding column variable vector for the variable swap.
-            std::vector<InternalBdd<DdType::CUDD>> columnVariables;
-            for (auto const& variable : rowVariables) {
-                columnVariables.push_back(InternalBdd<DdType::CUDD>(ddManager, ddManager->getCuddManager().bddVar(variable.getIndex() + 1)));
-            }
             result = result.swapVariables(rowVariables, columnVariables);
-            
             return result;
         }
         
-        InternalBdd<DdType::CUDD> InternalBdd<DdType::CUDD>::inverseRelationalProduct(InternalBdd<DdType::CUDD> const& relation, std::vector<InternalBdd<DdType::CUDD>> const& columnVariables) const {
+        InternalBdd<DdType::CUDD> InternalBdd<DdType::CUDD>::inverseRelationalProduct(InternalBdd<DdType::CUDD> const& relation, std::vector<InternalBdd<DdType::CUDD>> const& rowVariables, std::vector<InternalBdd<DdType::CUDD>> const& columnVariables) const {
             InternalBdd<DdType::CUDD> cube = ddManager->getBddOne();
             for (auto const& variable : columnVariables) {
                 cube &= variable;
             }
             
             InternalBdd<DdType::CUDD> result = this->andExists(relation, cube);
-            
-            // Create the corresponding column variable vector for the variable swap.
-            std::vector<InternalBdd<DdType::CUDD>> rowVariables;
-            for (auto const& variable : rowVariables) {
-                rowVariables.push_back(InternalBdd<DdType::CUDD>(ddManager, ddManager->getCuddManager().bddVar(variable.getIndex() - 1)));
-            }
             result = result.swapVariables(rowVariables, columnVariables);
-            
             return result;
         }
         
diff --git a/src/storage/dd/cudd/InternalCuddBdd.h b/src/storage/dd/cudd/InternalCuddBdd.h
index 2c8951aba..4358e491f 100644
--- a/src/storage/dd/cudd/InternalCuddBdd.h
+++ b/src/storage/dd/cudd/InternalCuddBdd.h
@@ -81,18 +81,20 @@ namespace storm {
              *
              * @param relation The relation to use.
              * @param rowVariables The row variables of the relation represented as individual BDDs.
+             * @param columnVariables The column variables of the relation represented as individual BDDs.
              * @return The ralational product.
              */
-            InternalBdd<DdType::CUDD> relationalProduct(InternalBdd<DdType::CUDD> const& relation, std::vector<InternalBdd<DdType::CUDD>> const& rowVariables) const;
+            InternalBdd<DdType::CUDD> relationalProduct(InternalBdd<DdType::CUDD> const& relation, std::vector<InternalBdd<DdType::CUDD>> const& rowVariables, std::vector<InternalBdd<DdType::CUDD>> const& columnVariables) const;
 
             /*!
              * Computes the inverse relational product of the current BDD and the given BDD representing a relation.
              *
              * @param relation The relation to use.
+             * @param rowVariables The row variables of the relation represented as individual BDDs.
              * @param columnVariables The row variables of the relation represented as individual BDDs.
              * @return The ralational product.
              */
-            InternalBdd<DdType::CUDD> inverseRelationalProduct(InternalBdd<DdType::CUDD> const& relation, std::vector<InternalBdd<DdType::CUDD>> const& columnVariables) const;
+            InternalBdd<DdType::CUDD> inverseRelationalProduct(InternalBdd<DdType::CUDD> const& relation, std::vector<InternalBdd<DdType::CUDD>> const& rowVariables, std::vector<InternalBdd<DdType::CUDD>> const& columnVariables) const;
 
             /*!
              * Performs an if-then-else with the given operands, i.e. maps all valuations that are mapped to a non-zero
diff --git a/src/storage/dd/sylvan/InternalSylvanBdd.cpp b/src/storage/dd/sylvan/InternalSylvanBdd.cpp
index f08273001..c2e72a97a 100644
--- a/src/storage/dd/sylvan/InternalSylvanBdd.cpp
+++ b/src/storage/dd/sylvan/InternalSylvanBdd.cpp
@@ -30,16 +30,25 @@ namespace storm {
             return sylvanBdd != other.sylvanBdd;
         }
         
-        InternalBdd<DdType::Sylvan> InternalBdd<DdType::Sylvan>::relationalProduct(InternalBdd<DdType::Sylvan> const& relation, std::vector<InternalBdd<DdType::Sylvan>> const& rowVariables) const {
+        InternalBdd<DdType::Sylvan> InternalBdd<DdType::Sylvan>::relationalProduct(InternalBdd<DdType::Sylvan> const& relation, std::vector<InternalBdd<DdType::Sylvan>> const& rowVariables, std::vector<InternalBdd<DdType::Sylvan>> const& columnVariables) const {
             InternalBdd<DdType::Sylvan> cube = ddManager->getBddOne();
             for (auto const& variable : rowVariables) {
                 cube &= variable;
             }
+            for (auto const& variable : columnVariables) {
+                cube &= variable;
+            }
+            
+            this->exportToDot("set.dot", {});
+            relation.exportToDot("relation.dot", {});
+            cube.exportToDot("cube.dot", {});
             
-            return InternalBdd<DdType::Sylvan>(ddManager, this->sylvanBdd.RelNext(relation.sylvanBdd, cube.sylvanBdd));
+            InternalBdd<DdType::Sylvan> result(ddManager, this->sylvanBdd.RelNext(relation.sylvanBdd, cube.sylvanBdd));
+            result.exportToDot("result.dot", {});
+            return result;
         }
 
-        InternalBdd<DdType::Sylvan> InternalBdd<DdType::Sylvan>::inverseRelationalProduct(InternalBdd<DdType::Sylvan> const& relation, std::vector<InternalBdd<DdType::Sylvan>> const& columnVariables) const {
+        InternalBdd<DdType::Sylvan> InternalBdd<DdType::Sylvan>::inverseRelationalProduct(InternalBdd<DdType::Sylvan> const& relation, std::vector<InternalBdd<DdType::Sylvan>> const& rowVariables, std::vector<InternalBdd<DdType::Sylvan>> const& columnVariables) const {
             STORM_LOG_THROW(false, storm::exceptions::NotImplementedException, "Not yet implemented.");
         }
         
diff --git a/src/storage/dd/sylvan/InternalSylvanBdd.h b/src/storage/dd/sylvan/InternalSylvanBdd.h
index 092321510..5777f391d 100644
--- a/src/storage/dd/sylvan/InternalSylvanBdd.h
+++ b/src/storage/dd/sylvan/InternalSylvanBdd.h
@@ -70,18 +70,20 @@ namespace storm {
              *
              * @param relation The relation to use.
              * @param rowVariables The row variables of the relation represented as individual BDDs.
+             * @param columnVariables The row variables of the relation represented as individual BDDs.
              * @return The ralational product.
              */
-            InternalBdd<DdType::Sylvan> relationalProduct(InternalBdd<DdType::Sylvan> const& relation, std::vector<InternalBdd<DdType::Sylvan>> const& rowVariables) const;
+            InternalBdd<DdType::Sylvan> relationalProduct(InternalBdd<DdType::Sylvan> const& relation, std::vector<InternalBdd<DdType::Sylvan>> const& rowVariables, std::vector<InternalBdd<DdType::Sylvan>> const& columnVariables) const;
             
             /*!
              * Computes the inverse relational product of the current BDD and the given BDD representing a relation.
              *
              * @param relation The relation to use.
+             * @param rowVariables The row variables of the relation represented as individual BDDs.
              * @param columnVariables The row variables of the relation represented as individual BDDs.
              * @return The ralational product.
              */
-            InternalBdd<DdType::Sylvan> inverseRelationalProduct(InternalBdd<DdType::Sylvan> const& relation, std::vector<InternalBdd<DdType::Sylvan>> const& columnVariables) const;
+            InternalBdd<DdType::Sylvan> inverseRelationalProduct(InternalBdd<DdType::Sylvan> const& relation, std::vector<InternalBdd<DdType::Sylvan>> const& rowVariables, std::vector<InternalBdd<DdType::Sylvan>> const& columnVariables) const;
             
             /*!
              * Performs an if-then-else with the given operands, i.e. maps all valuations that are mapped to a non-zero
diff --git a/src/storage/dd/sylvan/InternalSylvanDdManager.cpp b/src/storage/dd/sylvan/InternalSylvanDdManager.cpp
index 9f6df3e70..c30da4863 100644
--- a/src/storage/dd/sylvan/InternalSylvanDdManager.cpp
+++ b/src/storage/dd/sylvan/InternalSylvanDdManager.cpp
@@ -9,9 +9,9 @@ namespace storm {
     namespace dd {
         uint_fast64_t InternalDdManager<DdType::Sylvan>::numberOfInstances = 0;
         
-        // We let the variables start at an odd offset, since some functions provided by sylvan assume that the source/row
-        // variables are at odd levels.
-        uint_fast64_t InternalDdManager<DdType::Sylvan>::nextFreeVariableIndex = 1;
+        // It is important that the variable pairs start at an even offset, because sylvan assumes this to be true for
+        // some operations.
+        uint_fast64_t InternalDdManager<DdType::Sylvan>::nextFreeVariableIndex = 0;
         
         InternalDdManager<DdType::Sylvan>::InternalDdManager() {
             if (numberOfInstances == 0) {