diff --git a/src/builder/ExplicitPrismModelBuilder.cpp b/src/builder/ExplicitPrismModelBuilder.cpp
index 039048111..4a777c6f9 100644
--- a/src/builder/ExplicitPrismModelBuilder.cpp
+++ b/src/builder/ExplicitPrismModelBuilder.cpp
@@ -11,6 +11,8 @@
 
 #include "src/settings/modules/GeneralSettings.h"
 
+#include "src/generator/PrismNextStateGenerator.h"
+
 #include "src/utility/prism.h"
 #include "src/utility/macros.h"
 #include "src/utility/ConstantsComparator.h"
@@ -37,7 +39,7 @@ namespace storm {
                     stateRewardVector.resize(rowGroupCount);
                     optionalStateRewardVector = std::move(stateRewardVector);
                 }
-
+                
                 boost::optional<std::vector<ValueType>> optionalStateActionRewardVector;
                 if (hasStateActionRewards) {
                     stateActionRewardVector.resize(rowCount);
@@ -75,12 +77,12 @@ namespace storm {
         ExplicitPrismModelBuilder<ValueType, RewardModelType, StateType>::InternalStateInformation::InternalStateInformation(uint64_t bitsPerState) : stateStorage(bitsPerState, 10000000), bitsPerState(bitsPerState), numberOfStates() {
             // Intentionally left empty.
         }
-            
+        
         template <typename ValueType, typename RewardModelType, typename StateType>
         ExplicitPrismModelBuilder<ValueType, RewardModelType, StateType>::ModelComponents::ModelComponents() : transitionMatrix(), stateLabeling(), rewardModels(), choiceLabeling() {
             // Intentionally left empty.
         }
-
+        
         template <typename ValueType, typename RewardModelType, typename StateType>
         ExplicitPrismModelBuilder<ValueType, RewardModelType, StateType>::Options::Options() : buildCommandLabels(false), buildAllRewardModels(true), buildStateInformation(false), rewardModelsToBuild(), constantDefinitions(), buildAllLabels(true), labelsToBuild(), expressionLabels(), terminalStates(), negatedTerminalStates() {
             // Intentionally left empty.
@@ -135,7 +137,7 @@ namespace storm {
                 }
             }
         }
-
+        
         template <typename ValueType, typename RewardModelType, typename StateType>
         void ExplicitPrismModelBuilder<ValueType, RewardModelType, StateType>::Options::preserveFormula(storm::logic::Formula const& formula) {
             // If we already had terminal states, we need to erase them.
@@ -180,7 +182,7 @@ namespace storm {
         }
         
         template <typename ValueType, typename RewardModelType, typename StateType>
-        ExplicitPrismModelBuilder<ValueType, RewardModelType, StateType>::ExplicitPrismModelBuilder(storm::prism::Program const& program, Options const& options) : options(options), program(program) {
+        ExplicitPrismModelBuilder<ValueType, RewardModelType, StateType>::ExplicitPrismModelBuilder(storm::prism::Program const& program, Options const& options) : options(options), program(program), variableInformation(program), internalStateInformation(variableInformation.getTotalBitOffset()) {
             // Start by defining the undefined constants in the model.
             if (options.constantDefinitions) {
                 this->program = program.defineUndefinedConstants(options.constantDefinitions.get());
@@ -189,56 +191,48 @@ namespace storm {
             }
             
             // If the program still contains undefined constants and we are not in a parametric setting, assemble an appropriate error message.
-#ifdef STORM_HAVE_CARL
-            // If the program either has undefined constants or we are building a parametric model, but the parameters
-            // not only appear in the probabilities, we re
             if (!std::is_same<ValueType, storm::RationalFunction>::value && this->program.hasUndefinedConstants()) {
-#else
-                if (this->program->hasUndefinedConstants()) {
-#endif
-                    std::vector<std::reference_wrapper<storm::prism::Constant const>> undefinedConstants = this->program.getUndefinedConstants();
-                    std::stringstream stream;
-                    bool printComma = false;
-                    for (auto const& constant : undefinedConstants) {
-                        if (printComma) {
-                            stream << ", ";
-                        } else {
-                            printComma = true;
-                        }
-                        stream << constant.get().getName() << " (" << constant.get().getType() << ")";
+                std::vector<std::reference_wrapper<storm::prism::Constant const>> undefinedConstants = this->program.getUndefinedConstants();
+                std::stringstream stream;
+                bool printComma = false;
+                for (auto const& constant : undefinedConstants) {
+                    if (printComma) {
+                        stream << ", ";
+                    } else {
+                        printComma = true;
                     }
-                    stream << ".";
-                    STORM_LOG_THROW(false, storm::exceptions::InvalidArgumentException, "Program still contains these undefined constants: " + stream.str());
-#ifdef STORM_HAVE_CARL
-                } else if (std::is_same<ValueType, storm::RationalFunction>::value && !this->program.hasUndefinedConstantsOnlyInUpdateProbabilitiesAndRewards()) {
-                    STORM_LOG_THROW(false, storm::exceptions::InvalidArgumentException, "The program contains undefined constants that appear in some places other than update probabilities and reward value expressions, which is not admitted.");
-#endif
+                    stream << constant.get().getName() << " (" << constant.get().getType() << ")";
                 }
-                
-                // If the set of labels we are supposed to built is restricted, we need to remove the other labels from the program.
-                if (options.labelsToBuild) {
-                    if (!options.buildAllLabels) {
-                        this->program.filterLabels(options.labelsToBuild.get());
-                    }
+                stream << ".";
+                STORM_LOG_THROW(false, storm::exceptions::InvalidArgumentException, "Program still contains these undefined constants: " + stream.str());
+            } else if (std::is_same<ValueType, storm::RationalFunction>::value && !this->program.hasUndefinedConstantsOnlyInUpdateProbabilitiesAndRewards()) {
+                STORM_LOG_THROW(false, storm::exceptions::InvalidArgumentException, "The program contains undefined constants that appear in some places other than update probabilities and reward value expressions, which is not admitted.");
+            }
+            
+            // If the set of labels we are supposed to built is restricted, we need to remove the other labels from the program.
+            if (options.labelsToBuild) {
+                if (!options.buildAllLabels) {
+                    this->program.filterLabels(options.labelsToBuild.get());
                 }
+            }
+            
+            // If we need to build labels for expressions that may appear in some formula, we need to add appropriate
+            // labels to the program.
+            if (options.expressionLabels) {
+                std::map<storm::expressions::Variable, storm::expressions::Expression> constantsSubstitution = this->program.getConstantsSubstitution();
                 
-                // If we need to build labels for expressions that may appear in some formula, we need to add appropriate
-                // labels to the program.
-                if (options.expressionLabels) {
-                    std::map<storm::expressions::Variable, storm::expressions::Expression> constantsSubstitution = this->program.getConstantsSubstitution();
-                    
-                    for (auto const& expression : options.expressionLabels.get()) {
-                        std::stringstream stream;
-                        stream << expression.substitute(constantsSubstitution);
-                        std::string name = stream.str();
-                        if (!this->program.hasLabel(name)) {
-                            this->program.addLabel(name, expression);
-                        }
+                for (auto const& expression : options.expressionLabels.get()) {
+                    std::stringstream stream;
+                    stream << expression.substitute(constantsSubstitution);
+                    std::string name = stream.str();
+                    if (!this->program.hasLabel(name)) {
+                        this->program.addLabel(name, expression);
                     }
                 }
-                
-                // Now that the program is fixed, we we need to substitute all constants with their concrete value.
-                this->program = program.substituteConstants();
+            }
+            
+            // Now that the program is fixed, we we need to substitute all constants with their concrete value.
+            this->program = program.substituteConstants();
         }
         
         template <typename ValueType, typename RewardModelType, typename StateType>
@@ -255,15 +249,15 @@ namespace storm {
         template <typename ValueType, typename RewardModelType, typename StateType>
         std::shared_ptr<storm::models::sparse::Model<ValueType, RewardModelType>> ExplicitPrismModelBuilder<ValueType, RewardModelType, StateType>::translate() {
             STORM_LOG_DEBUG("Building representation of program:" << std::endl << *program << std::endl);
-                
+            
             // Select the appropriate reward models (after the constants have been substituted).
             std::vector<std::reference_wrapper<storm::prism::RewardModel const>> selectedRewardModels;
-                
+            
             // First, we make sure that all selected reward models actually exist.
             for (auto const& rewardModelName : options.rewardModelsToBuild) {
                 STORM_LOG_THROW(rewardModelName.empty() || program.hasRewardModel(rewardModelName), storm::exceptions::InvalidArgumentException, "Model does not possess a reward model with the name '" << rewardModelName << "'.");
             }
-                
+            
             for (auto const& rewardModel : program.getRewardModels()) {
                 if (options.buildAllRewardModels || options.rewardModelsToBuild.find(rewardModel.getName()) != options.rewardModelsToBuild.end()) {
                     selectedRewardModels.push_back(rewardModel);
@@ -274,7 +268,7 @@ namespace storm {
             if (selectedRewardModels.empty() && program.getNumberOfRewardModels() == 1 && options.rewardModelsToBuild.size() == 1 && *options.rewardModelsToBuild.begin() == "") {
                 selectedRewardModels.push_back(program.getRewardModel(0));
             }
-                
+            
             ModelComponents modelComponents = buildModelComponents(*program, selectedRewardModels, options);
             
             std::shared_ptr<storm::models::sparse::Model<ValueType, RewardModelType>> result;
@@ -295,10 +289,10 @@ namespace storm {
             
             return result;
         }
-                    
+        
         template <typename ValueType, typename RewardModelType, typename StateType>
         storm::expressions::SimpleValuation ExplicitPrismModelBuilder<ValueType, RewardModelType, StateType>::unpackStateIntoValuation(storm::storage::BitVector const& currentState) {
-            storm::expressions::SimpleValuation result(variableInformation.manager.getSharedPointer());
+            storm::expressions::SimpleValuation result(program.getManager().getSharedPointer());
             for (auto const& booleanVariable : variableInformation.booleanVariables) {
                 result.setBooleanValue(booleanVariable.variable, currentState.get(booleanVariable.bitOffset));
             }
@@ -322,17 +316,54 @@ namespace storm {
             
             return actualIndexBucketPair.first;
         }
-
+        
         template <typename ValueType, typename RewardModelType, typename StateType>
-        boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>> ExplicitPrismModelBuilder<ValueType, RewardModelType, StateType>::buildMatrices(std::vector<std::reference_wrapper<storm::prism::RewardModel const>> const& selectedRewardModels, InternalStateInformation& internalStateInformation, storm::storage::SparseMatrixBuilder<ValueType>& transitionMatrixBuilder, std::vector<RewardModelBuilder<typename RewardModelType::ValueType>>& rewardModelBuilders, boost::optional<storm::expressions::Expression> const& terminalExpression) {
+        boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>> ExplicitPrismModelBuilder<ValueType, RewardModelType, StateType>::buildMatrices(std::vector<std::reference_wrapper<storm::prism::RewardModel const>> const& selectedRewardModels, storm::storage::SparseMatrixBuilder<ValueType>& transitionMatrixBuilder, std::vector<RewardModelBuilder<typename RewardModelType::ValueType>>& rewardModelBuilders, boost::optional<storm::expressions::Expression> const& terminalExpression) {
             // Create choice labels, if requested,
             boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>> choiceLabels;
             if (options.buildCommandLabels) {
                 choiceLabels = std::vector<boost::container::flat_set<uint_fast64_t>>();
             }
+
+            // Create a generator that is able to expand states.
+            storm::generator::PrismNextStateGenerator<ValueType, StateType> generator(program, variableInformation, options.buildCommandLabels);
+            for (auto const& rewardModel : selectedRewardModels) {
+                generator.addRewardModel(rewardModel.get());
+            }
+
+            // Create a callback for the next-state generator to enable it to request the index of states.
+            std::function<StateType (CompressedState const&)> stateToIdCallback = std::bind1st(&ExplicitPrismModelBuilder<ValueType, RewardModelType, StateType>::getOrAddStateIndex, this);
             
-            // A comparator that can be used to check whether we actually have distributions.
-            storm::utility::ConstantsComparator<ValueType> comparator;
+            // Let the generator create all initial states.
+            generator.getInitialStates(stateToIdCallback);
+            
+            // Now explore the current state until there is no more reachable state.
+            uint_fast64_t currentRow = 0;
+
+            // Perform a DFS through the model.
+            while (!statesToExplore.empty()) {
+                // Get the first state in the queue.
+                std::pair<CompressedState, StateType> currentState = internalStateInformation.stateStorage.getBucketAndValue(statesToExplore.front());
+                statesToExplore.pop();
+                
+                STORM_LOG_TRACE("Exploring state with id " << currentState.second << ".");
+                
+                bool deadlockOnPurpose = false;
+                if (terminalExpression && evaluator.asBool(terminalExpression.get())) {
+                    // Nothing to do in this case.
+                    deadlockOnPurpose = true;
+                } else {
+                    // Retrieve all choices for the current state.
+                    allUnlabeledChoices = getUnlabeledTransitions(program, discreteTimeModel, internalStateInformation, variableInformation, currentState, commandLabels, evaluator, stateQueue, comparator);
+                    allLabeledChoices = getLabeledTransitions(program, discreteTimeModel, internalStateInformation, variableInformation, currentState, commandLabels, evaluator, stateQueue, comparator);
+                }
+
+
+            }
+
+            
+                
+                
             
             // Initialize a queue and insert the initial state.
             std::queue<storm::storage::BitVector> stateQueue;
@@ -420,7 +451,7 @@ namespace storm {
                         std::cout << unpackStateIntoValuation(currentState, variableInformation).toString(true) << std::endl;
                         STORM_LOG_THROW(false, storm::exceptions::WrongFormatException,
                                         "Error while creating sparse matrix from probabilistic program: found deadlock state. For fixing these, please provide the appropriate option.");
-
+                        
                     }
                 } else if (totalNumberOfChoices == 1) {
                     if (!deterministicModel) {
@@ -468,7 +499,7 @@ namespace storm {
                     // or compose them to one choice.
                     if (deterministicModel) {
                         Choice<ValueType> globalChoice;
-
+                        
                         // We need to prepare the entries of those vectors that are going to be used.
                         auto builderIt = rewardModelBuilders.begin();
                         for (auto rewardModelIt = selectedRewardModels.begin(), rewardModelIte = selectedRewardModels.end(); rewardModelIt != rewardModelIte; ++rewardModelIt, ++builderIt) {
@@ -571,7 +602,7 @@ namespace storm {
                     } else {
                         // If the model is nondeterministic, we add all choices individually.
                         transitionMatrixBuilder.newRowGroup(currentRow);
-
+                        
                         auto builderIt = rewardModelBuilders.begin();
                         for (auto rewardModelIt = selectedRewardModels.begin(), rewardModelIte = selectedRewardModels.end(); rewardModelIt != rewardModelIte; ++rewardModelIt, ++builderIt) {
                             if (rewardModelIt->get().hasStateRewards()) {
@@ -651,36 +682,7 @@ namespace storm {
         typename ExplicitPrismModelBuilder<ValueType, RewardModelType, StateType>::ModelComponents ExplicitPrismModelBuilder<ValueType, RewardModelType, StateType>::buildModelComponents(storm::prism::Program const& program, std::vector<std::reference_wrapper<storm::prism::RewardModel const>> const& selectedRewardModels, Options const& options) {
             ModelComponents modelComponents;
             
-            uint_fast64_t bitOffset = 0;
-            VariableInformation variableInformation(program.getManager());
-            for (auto const& booleanVariable : program.getGlobalBooleanVariables()) {
-                variableInformation.booleanVariables.emplace_back(booleanVariable.getExpressionVariable(), booleanVariable.getInitialValueExpression().evaluateAsBool(), bitOffset);
-                ++bitOffset;
-                variableInformation.booleanVariableToIndexMap[booleanVariable.getExpressionVariable()] = variableInformation.booleanVariables.size() - 1;
-            }
-            for (auto const& integerVariable : program.getGlobalIntegerVariables()) {
-                int_fast64_t lowerBound = integerVariable.getLowerBoundExpression().evaluateAsInt();
-                int_fast64_t upperBound = integerVariable.getUpperBoundExpression().evaluateAsInt();
-                uint_fast64_t bitwidth = static_cast<uint_fast64_t>(std::ceil(std::log2(upperBound - lowerBound + 1)));
-                variableInformation.integerVariables.emplace_back(integerVariable.getExpressionVariable(), integerVariable.getInitialValueExpression().evaluateAsInt(), lowerBound, upperBound, bitOffset, bitwidth);
-                bitOffset += bitwidth;
-                variableInformation.integerVariableToIndexMap[integerVariable.getExpressionVariable()] = variableInformation.integerVariables.size() - 1;
-            }
-            for (auto const& module : program.getModules()) {
-                for (auto const& booleanVariable : module.getBooleanVariables()) {
-                    variableInformation.booleanVariables.emplace_back(booleanVariable.getExpressionVariable(), booleanVariable.getInitialValueExpression().evaluateAsBool(), bitOffset);
-                    ++bitOffset;
-                    variableInformation.booleanVariableToIndexMap[booleanVariable.getExpressionVariable()] = variableInformation.booleanVariables.size() - 1;
-                }
-                for (auto const& integerVariable : module.getIntegerVariables()) {
-                    int_fast64_t lowerBound = integerVariable.getLowerBoundExpression().evaluateAsInt();
-                    int_fast64_t upperBound = integerVariable.getUpperBoundExpression().evaluateAsInt();
-                    uint_fast64_t bitwidth = static_cast<uint_fast64_t>(std::ceil(std::log2(upperBound - lowerBound + 1)));
-                    variableInformation.integerVariables.emplace_back(integerVariable.getExpressionVariable(), integerVariable.getInitialValueExpression().evaluateAsInt(), lowerBound, upperBound, bitOffset, bitwidth);
-                    bitOffset += bitwidth;
-                    variableInformation.integerVariableToIndexMap[integerVariable.getExpressionVariable()] = variableInformation.integerVariables.size() - 1;
-                }
-            }
+            VariableInformation variableInformation(program);
             
             // Create the structure for storing the reachable state space.
             uint64_t bitsPerState = ((bitOffset / 64) + 1) * 64;
diff --git a/src/builder/ExplicitPrismModelBuilder.h b/src/builder/ExplicitPrismModelBuilder.h
index 5e4f5b2c5..21c45e30b 100644
--- a/src/builder/ExplicitPrismModelBuilder.h
+++ b/src/builder/ExplicitPrismModelBuilder.h
@@ -242,7 +242,7 @@ namespace storm {
              * @return A tuple containing a vector with all rows at which the nondeterministic choices of each state begin
              * and a vector containing the labels associated with each choice.
              */
-            boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>> buildMatrices(std::vector<std::reference_wrapper<storm::prism::RewardModel const>> const& selectedRewardModels, InternalStateInformation& internalStateInformation, storm::storage::SparseMatrixBuilder<ValueType>& transitionMatrixBuilder, std::vector<RewardModelBuilder<typename RewardModelType::ValueType>>& rewardModelBuilders, boost::optional<storm::expressions::Expression> const& terminalExpression);
+            boost::optional<std::vector<boost::container::flat_set<uint_fast64_t>>> buildMatrices(std::vector<std::reference_wrapper<storm::prism::RewardModel const>> const& selectedRewardModels, storm::storage::SparseMatrixBuilder<ValueType>& transitionMatrixBuilder, std::vector<RewardModelBuilder<typename RewardModelType::ValueType>>& rewardModelBuilders, boost::optional<storm::expressions::Expression> const& terminalExpression);
             
             /*!
              * Explores the state space of the given program and returns the components of the model as a result.
@@ -280,8 +280,9 @@ namespace storm {
             // successful build.
             boost::optional<StateInformation> stateInformation;
             
-            // A queue of states that still need to be explored.
-            std::queue<storm::storage::BitVector> statesToExplore;
+            // A queue of states that still need to be explored. The indices in this queue are the bucket indices in the
+            // bit vector hash map holding the compressed states.
+            std::queue<std::size_t> statesToExplore;
 
         };
         
diff --git a/src/generator/NextStateGenerator.h b/src/generator/NextStateGenerator.h
index 8e8e4bbdb..4248ab90a 100644
--- a/src/generator/NextStateGenerator.h
+++ b/src/generator/NextStateGenerator.h
@@ -16,10 +16,10 @@ namespace storm {
         template<typename ValueType, typename StateType = uint32_t>
         class NextStateGenerator {
         public:
-            typedef StateType (*StateToIdCallback)(CompressedState const&);
-
-            virtual std::vector<StateType> getInitialStates(StateToIdCallback stateToIdCallback) = 0;
-            virtual StateBehavior<ValueType, StateType> expand(CompressedState const& state, StateToIdCallback stateToIdCallback) = 0;
+            typedef std::function<StateType (CompressedState const&)> StateToIdCallback;
+            
+            virtual std::vector<StateType> getInitialStates(StateToIdCallback const& stateToIdCallback) = 0;
+            virtual StateBehavior<ValueType, StateType> expand(CompressedState const& state, StateToIdCallback const& stateToIdCallback) = 0;
         };
     }
 }
diff --git a/src/generator/PrismNextStateGenerator.cpp b/src/generator/PrismNextStateGenerator.cpp
index a23ed29a6..602abb678 100644
--- a/src/generator/PrismNextStateGenerator.cpp
+++ b/src/generator/PrismNextStateGenerator.cpp
@@ -15,13 +15,108 @@ namespace storm {
         template<typename ValueType, typename StateType>
         void PrismNextStateGenerator<ValueType, StateType>::addRewardModel(storm::prism::RewardModel const& rewardModel) {
             selectedRewardModels.push_back(rewardModel);
+            hasStateActionRewards |= rewardModel.hasStateActionRewards();
         }
         
         template<typename ValueType, typename StateType>
-        StateBehavior<ValueType, StateType> PrismNextStateGenerator<ValueType, StateType>::expand(CompressedState const& state, typename NextStateGenerator<ValueType, StateType>::StateToIdCallback stateToIdCallback) {
-            // TODO
+        std::vector<StateType> PrismNextStateGenerator<ValueType, StateType>::getInitialStates(StateToIdCallback const& stateToIdCallback) {
+            // FIXME, TODO, whatever
         }
         
+        template<typename ValueType, typename StateType>
+        StateBehavior<ValueType, StateType> PrismNextStateGenerator<ValueType, StateType>::expand(CompressedState const& state, StateToIdCallback const& stateToIdCallback) {
+            // Start by unpacking the state into the evaluator so we can quickly evaluate expressions later.
+            unpackStateIntoEvaluator(state);
+            
+            // Prepare the result, in case we return early.
+            StateBehavior<ValueType, StateType> result;
+            
+            // First, construct the state rewards, as we may return early if there are no choices later and we already
+            // need the state rewards then.
+            for (auto const& rewardModel : selectedRewardModels) {
+                ValueType stateReward = storm::utility::zero<ValueType>();
+                if (rewardModel->hasStateRewards()) {
+                    for (auto const& stateReward : rewardModel->getStateRewards()) {
+                        if (evaluator.asBool(stateReward.getStatePredicateExpression())) {
+                            stateReward += ValueType(evaluator.asRational(stateReward.getRewardValueExpression()));
+                        }
+                    }
+                }
+                result.addStateReward(stateReward);
+            }
+            
+            // Get all choices for the state.
+            std::vector<Choice<ValueType>> allChoices = getUnlabeledChoices(state, stateToIdCallback);
+            std::vector<Choice<ValueType>> allLabeledChoices = getLabeledChoices(state, stateToIdCallback);
+            for (auto& choice : allLabeledChoices) {
+                allChoices.push_back(std::move(choice));
+            }
+            
+            std::size_t totalNumberOfChoices = allChoices.size();
+            
+            // If there is not a single choice, we return immediately, because the state has no behavior (other than
+            // the state reward).
+            if (totalNumberOfChoices == 0) {
+                return result;
+            }
+            
+            // If the model is a deterministic model, we need to fuse the choices into one.
+            if (program.isDeterministicModel() && totalNumberOfChoices > 1) {
+                Choice<ValueType> globalChoice;
+                
+                // For CTMCs, we need to keep track of the total exit rate to scale the action rewards later. For DTMCs
+                // this is equal to the number of choices, which is why we initialize it like this here.
+                ValueType totalExitRate = static_cast<ValueType>(totalNumberOfChoices);
+                
+                // Iterate over all choices and combine the probabilities/rates into one choice.
+                for (auto const& choice : allChoices) {
+                    for (auto const& stateProbabilityPair : choice) {
+                        if (program.isDiscreteTimeModel()) {
+                            globalChoice.getOrAddEntry(stateProbabilityPair.first) += stateProbabilityPair.second / totalNumberOfChoices;
+                            if (hasStateActionRewards) {
+                                totalExitRate += choice.getTotalMass();
+                            }
+                        } else {
+                            globalChoice.getOrAddEntry(stateProbabilityPair.first) += stateProbabilityPair.second;
+                        }
+                    }
+                    
+                    if (buildChoiceLabeling) {
+                        globalChoice.addChoiceLabels(choice.getChoiceLabels());
+                    }
+                }
+                
+                // Now construct the state-action reward for all selected reward models.
+                for (auto const& rewardModel : selectedRewardModels) {
+                    ValueType stateActionReward = storm::utility::zero<ValueType>();
+                    if (rewardModel->hasStateActionRewards()) {
+                        for (auto const& stateActionReward : rewardModel->getStateActionRewards()) {
+                            for (auto const& choice : allChoices) {
+                                if (stateActionReward.getActionIndex() == choice.getActionIndex() && evaluator.asBool(stateActionReward.getStatePredicateExpression())) {
+                                    stateActionReward += ValueType(evaluator.asRational(stateActionReward.getRewardValueExpression())) / totalExitRate;
+                                }
+                            }
+                            
+                        }
+                    }
+                    globalChoice.addChoiceReward(stateActionReward);
+                }
+                
+                // Move the newly fused choice in place.
+                allChoices.clear();
+                allChoices.push_back(std::move(globalChoice));
+            }
+            
+            // Move all remaining choices in place.
+            for (auto& choice : allChoices) {
+                result.addChoice(std::move(choice));
+            }
+            
+            return result;
+        }
+        
+        
+        
         template<typename ValueType, typename StateType>
         void PrismNextStateGenerator<ValueType, StateType>::unpackStateIntoEvaluator(storm::storage::BitVector const& state) {
             for (auto const& booleanVariable : variableInformation.booleanVariables) {
@@ -31,7 +126,7 @@ namespace storm {
                 evaluator.setIntegerValue(integerVariable.variable, state.getAsInt(integerVariable.bitOffset, integerVariable.bitWidth) + integerVariable.lowerBound);
             }
         }
-
+        
         template<typename ValueType, typename StateType>
         CompressedState PrismNextStateGenerator<ValueType, StateType>::applyUpdate(CompressedState const& state, storm::prism::Update const& update) {
             CompressedState newState(state);
@@ -110,7 +205,7 @@ namespace storm {
         }
         
         template<typename ValueType, typename StateType>
-        std::vector<Choice<ValueType>> PrismNextStateGenerator<ValueType, StateType>::getUnlabeledTransitions(CompressedState const& state, StateToIdCallback stateToIdCallback) {
+        std::vector<Choice<ValueType>> PrismNextStateGenerator<ValueType, StateType>::getUnlabeledChoices(CompressedState const& state, StateToIdCallback stateToIdCallback) {
             std::vector<Choice<ValueType>> result;
             
             // Iterate over all modules.
@@ -129,7 +224,7 @@ namespace storm {
                         continue;
                     }
                     
-                    result.push_back(Choice<ValueType>(0, buildChoiceLabeling));
+                    result.push_back(Choice<ValueType>(command.getActionIndex()));
                     Choice<ValueType>& choice = result.back();
                     
                     // Remember the command labels only if we were asked to.
@@ -152,6 +247,19 @@ namespace storm {
                         probabilitySum += probability;
                     }
                     
+                    // Create the state-action reward for the newly created choice.
+                    for (auto const& rewardModel : selectedRewardModels) {
+                        ValueType stateActionReward = storm::utility::zero<ValueType>();
+                        if (rewardModel->hasStateActionRewards()) {
+                            for (auto const& stateActionReward : rewardModel->getStateActionRewards()) {
+                                if (stateActionReward.getActionIndex() == choice.getActionIndex() && evaluator.asBool(stateActionReward.getStatePredicateExpression())) {
+                                    stateActionReward += ValueType(evaluator.asRational(stateActionReward.getRewardValueExpression())) * choice.getTotalMass();
+                                }
+                            }
+                        }
+                        choice.addChoiceReward(stateActionReward);
+                    }
+                    
                     // Check that the resulting distribution is in fact a distribution.
                     STORM_LOG_THROW(!program.isDiscreteTimeModel() || comparator.isOne(probabilitySum), storm::exceptions::WrongFormatException, "Probabilities do not sum to one for command '" << command << "' (actually sum to " << probabilitySum << ").");
                 }
@@ -161,7 +269,7 @@ namespace storm {
         }
         
         template<typename ValueType, typename StateType>
-        std::vector<Choice<ValueType>> PrismNextStateGenerator<ValueType, StateType>::getLabeledTransitions(CompressedState const& state, StateToIdCallback stateToIdCallback) {
+        std::vector<Choice<ValueType>> PrismNextStateGenerator<ValueType, StateType>::getLabeledChoices(CompressedState const& state, StateToIdCallback stateToIdCallback) {
             std::vector<Choice<ValueType>> result;
             
             for (uint_fast64_t actionIndex : program.getSynchronizingActionIndices()) {
@@ -208,7 +316,7 @@ namespace storm {
                         // At this point, we applied all commands of the current command combination and newTargetStates
                         // contains all target states and their respective probabilities. That means we are now ready to
                         // add the choice to the list of transitions.
-                        result.push_back(Choice<ValueType>(actionIndex, buildChoiceLabeling));
+                        result.push_back(Choice<ValueType>(actionIndex));
                         
                         // Now create the actual distribution.
                         Choice<ValueType>& choice = result.back();
@@ -221,6 +329,7 @@ namespace storm {
                             }
                         }
                         
+                        // Add the probabilities/rates to the newly created choice.
                         ValueType probabilitySum = storm::utility::zero<ValueType>();
                         for (auto const& stateProbabilityPair : *newTargetStates) {
                             StateType actualIndex = stateToIdCallback(stateProbabilityPair.first);
@@ -231,6 +340,19 @@ namespace storm {
                         // Check that the resulting distribution is in fact a distribution.
                         STORM_LOG_THROW(!program.isDiscreteTimeModel() || !comparator.isConstant(probabilitySum) || comparator.isOne(probabilitySum), storm::exceptions::WrongFormatException, "Sum of update probabilities do not some to one for some command (actually sum to " << probabilitySum << ").");
                         
+                        // Create the state-action reward for the newly created choice.
+                        for (auto const& rewardModel : selectedRewardModels) {
+                            ValueType stateActionReward = storm::utility::zero<ValueType>();
+                            if (rewardModel->hasStateActionRewards()) {
+                                for (auto const& stateActionReward : rewardModel->getStateActionRewards()) {
+                                    if (stateActionReward.getActionIndex() == choice.getActionIndex() && evaluator.asBool(stateActionReward.getStatePredicateExpression())) {
+                                        stateActionReward += ValueType(evaluator.asRational(stateActionReward.getRewardValueExpression())) * choice.getTotalMass();
+                                    }
+                                }
+                            }
+                            choice.addChoiceReward(stateActionReward);
+                        }
+                        
                         // Dispose of the temporary maps.
                         delete currentTargetStates;
                         delete newTargetStates;
diff --git a/src/generator/PrismNextStateGenerator.h b/src/generator/PrismNextStateGenerator.h
index d3660ce05..282a2f25f 100644
--- a/src/generator/PrismNextStateGenerator.h
+++ b/src/generator/PrismNextStateGenerator.h
@@ -22,8 +22,8 @@ namespace storm {
              */
             void addRewardModel(storm::prism::RewardModel const& rewardModel);
             
-            virtual std::vector<StateType> getInitialStates(StateToIdCallback stateToIdCallback) = 0;
-            virtual StateBehavior<ValueType, StateType> expand(CompressedState const& state, StateToIdCallback stateToIdCallback) override;
+            virtual std::vector<StateType> getInitialStates(StateToIdCallback const& stateToIdCallback) override;
+            virtual StateBehavior<ValueType, StateType> expand(CompressedState const& state, StateToIdCallback const& stateToIdCallback) override;
                         
         private:
             /*!
@@ -82,6 +82,9 @@ namespace storm {
             // The reward models that need to be considered.
             std::vector<std::reference_wrapper<storm::prism::RewardModel const>> selectedRewardModels;
             
+            // A flag that stores whether at least one of the selected reward models has state-action rewards.
+            bool hasStateActionRewards;
+            
             // A flag that stores whether or not to build the choice labeling.
             bool buildChoiceLabeling;
 
diff --git a/src/generator/VariableInformation.cpp b/src/generator/VariableInformation.cpp
index 7f7ed357e..bb20c5225 100644
--- a/src/generator/VariableInformation.cpp
+++ b/src/generator/VariableInformation.cpp
@@ -1,4 +1,4 @@
-#include "src/generator/prism/VariableInformation.h"
+#include "src/generator/VariableInformation.h"
 
 #include "src/utility/macros.h"
 #include "src/exceptions/InvalidArgumentException.h"
@@ -14,8 +14,39 @@ namespace storm {
             // Intentionally left empty.
         }
         
-        VariableInformation::VariableInformation(storm::expressions::ExpressionManager const& manager) : manager(manager) {
-            // Intentionally left empty.
+        VariableInformation::VariableInformation(storm::prism::Program const& program) : totalBitOffset(0) {
+            for (auto const& booleanVariable : program.getGlobalBooleanVariables()) {
+                booleanVariables.emplace_back(booleanVariable.getExpressionVariable(), booleanVariable.getInitialValueExpression().evaluateAsBool(), totalBitOffset);
+                ++totalBitOffset;
+                booleanVariableToIndexMap[booleanVariable.getExpressionVariable()] = booleanVariables.size() - 1;
+            }
+            for (auto const& integerVariable : program.getGlobalIntegerVariables()) {
+                int_fast64_t lowerBound = integerVariable.getLowerBoundExpression().evaluateAsInt();
+                int_fast64_t upperBound = integerVariable.getUpperBoundExpression().evaluateAsInt();
+                uint_fast64_t bitwidth = static_cast<uint_fast64_t>(std::ceil(std::log2(upperBound - lowerBound + 1)));
+                integerVariables.emplace_back(integerVariable.getExpressionVariable(), integerVariable.getInitialValueExpression().evaluateAsInt(), lowerBound, upperBound, totalBitOffset, bitwidth);
+                totalBitOffset += bitwidth;
+                integerVariableToIndexMap[integerVariable.getExpressionVariable()] = integerVariables.size() - 1;
+            }
+            for (auto const& module : program.getModules()) {
+                for (auto const& booleanVariable : module.getBooleanVariables()) {
+                    booleanVariables.emplace_back(booleanVariable.getExpressionVariable(), booleanVariable.getInitialValueExpression().evaluateAsBool(), totalBitOffset);
+                    ++totalBitOffset;
+                    booleanVariableToIndexMap[booleanVariable.getExpressionVariable()] = booleanVariables.size() - 1;
+                }
+                for (auto const& integerVariable : module.getIntegerVariables()) {
+                    int_fast64_t lowerBound = integerVariable.getLowerBoundExpression().evaluateAsInt();
+                    int_fast64_t upperBound = integerVariable.getUpperBoundExpression().evaluateAsInt();
+                    uint_fast64_t bitwidth = static_cast<uint_fast64_t>(std::ceil(std::log2(upperBound - lowerBound + 1)));
+                    integerVariables.emplace_back(integerVariable.getExpressionVariable(), integerVariable.getInitialValueExpression().evaluateAsInt(), lowerBound, upperBound, totalBitOffset, bitwidth);
+                    totalBitOffset += bitwidth;
+                    integerVariableToIndexMap[integerVariable.getExpressionVariable()] = integerVariables.size() - 1;
+                }
+            }
+        }
+        
+        uint_fast64_t VariableInformation::getTotalBitOffset() const {
+            return totalBitOffset;
         }
         
         uint_fast64_t VariableInformation::getBitOffset(storm::expressions::Variable const& variable) const {
diff --git a/src/generator/VariableInformation.h b/src/generator/VariableInformation.h
index c1f06e741..a06156336 100644
--- a/src/generator/VariableInformation.h
+++ b/src/generator/VariableInformation.h
@@ -5,6 +5,7 @@
 #include <boost/container/flat_map.hpp>
 
 #include "src/storage/expressions/Variable.h"
+#include "src/storage/prism/Program.h"
 
 namespace storm {
     namespace generator {
@@ -48,12 +49,16 @@ namespace storm {
         
         // A structure storing information about the used variables of the program.
         struct VariableInformation {
-            VariableInformation(storm::expressions::ExpressionManager const& manager);
+            VariableInformation(storm::prism::Program const& program);
+            uint_fast64_t getTotalBitOffset() const;
             
             // Provide methods to access the bit offset and width of variables in the compressed state.
             uint_fast64_t getBitOffset(storm::expressions::Variable const& variable) const;
             uint_fast64_t getBitWidth(storm::expressions::Variable const& variable) const;
             
+            // The total bit offset over all variables.
+            uint_fast64_t totalBitOffset;
+            
             // The known boolean variables.
             boost::container::flat_map<storm::expressions::Variable, uint_fast64_t> booleanVariableToIndexMap;
             std::vector<BooleanVariableInformation> booleanVariables;
@@ -61,8 +66,6 @@ namespace storm {
             // The known integer variables.
             boost::container::flat_map<storm::expressions::Variable, uint_fast64_t> integerVariableToIndexMap;
             std::vector<IntegerVariableInformation> integerVariables;
-            
-            storm::expressions::ExpressionManager const& manager;
         };
         
     }
diff --git a/src/storage/BitVector.cpp b/src/storage/BitVector.cpp
index cd5f46009..38b94d852 100644
--- a/src/storage/BitVector.cpp
+++ b/src/storage/BitVector.cpp
@@ -60,7 +60,7 @@ namespace storm {
             return currentIndex == other.currentIndex;
         }
 
-        BitVector::BitVector() : bitCount(0), bucketVector() {
+        BitVector::BitVector() : bitCount(0), buckets(nullptr) {
             // Intentionally left empty.
         }
 
@@ -73,10 +73,17 @@ namespace storm {
 
             // Initialize the storage with the required values.
             if (init) {
-                bucketVector = std::vector<uint64_t>(bucketCount, -1ll);
+                buckets = new uint64_t[bucketCount];
+                std::fill_n(buckets, bucketCount, -1ull);
                 truncateLastBucket();
             } else {
-                bucketVector = std::vector<uint64_t>(bucketCount, 0);
+                buckets = new uint64_t[bucketCount]();
+            }
+        }
+        
+        BitVector::~BitVector() {
+            if (buckets != nullptr) {
+                delete buckets;
             }
         }
 
@@ -89,23 +96,22 @@ namespace storm {
             // Intentionally left empty.
         }
 
-        BitVector::BitVector(uint_fast64_t bucketCount, uint_fast64_t bitCount) : bitCount(bitCount), bucketVector(bucketCount) {
+        BitVector::BitVector(uint_fast64_t bucketCount, uint_fast64_t bitCount) : bitCount(bitCount) {
             STORM_LOG_ASSERT((bucketCount << 6) == bitCount, "Bit count does not match number of buckets.");
+            buckets = new uint64_t[bucketCount]();
         }
 
-        BitVector::BitVector(BitVector const& other) : bitCount(other.bitCount), bucketVector(other.bucketVector) {
-            // Intentionally left empty.
+        BitVector::BitVector(BitVector const& other) : bitCount(other.bitCount) {
+            buckets = new uint64_t[other.bucketCount()];
+            std::copy_n(other.buckets, other.bucketCount(), buckets);
         }
 
-        //BitVector::BitVector(BitVector&& other) : bitCount(other.bitCount), bucketVector(std::move(other.bucketVector)) {
-        // Intentionally left empty.
-        //}
-
         BitVector& BitVector::operator=(BitVector const& other) {
             // Only perform the assignment if the source and target are not identical.
             if (this != &other) {
                 bitCount = other.bitCount;
-                bucketVector = std::vector<uint64_t>(other.bucketVector);
+                buckets = new uint64_t[other.bucketCount()];
+                std::copy_n(other.buckets, other.bucketCount(), buckets);
             }
             return *this;
         }
@@ -117,9 +123,9 @@ namespace storm {
                 return false;
             }
 
-            std::vector<uint64_t>::const_iterator first1 = this->bucketVector.begin();
-            std::vector<uint64_t>::const_iterator last1 = this->bucketVector.end();
-            std::vector<uint64_t>::const_iterator first2 = other.bucketVector.begin();
+            uint64_t* first1 = this->buckets;
+            uint64_t* last1 = this->buckets + this->bucketCount();
+            uint64_t* first2 = other.buckets;
 
             for (; first1 != last1; ++first1, ++first2) {
                 if (*first1 < *first2) {
@@ -135,7 +141,8 @@ namespace storm {
             // Only perform the assignment if the source and target are not identical.
             if (this != &other) {
                 bitCount = other.bitCount;
-                bucketVector = std::move(other.bucketVector);
+                this->buckets = other.buckets;
+                other.buckets = nullptr;
             }
 
             return *this;
@@ -146,14 +153,7 @@ namespace storm {
             if (this->bitCount != other.bitCount) return false;
 
             // If the lengths match, we compare the buckets one by one.
-            for (std::vector<uint64_t>::const_iterator it1 = bucketVector.begin(), it2 = other.bucketVector.begin(); it1 != bucketVector.end(); ++it1, ++it2) {
-                if (*it1 != *it2) {
-                    return false;
-                }
-            }
-
-            // All buckets were equal, so the bit vectors are equal.
-            return true;
+            return std::equal(this->buckets, this->buckets + this->bucketCount(), other.buckets);
         }
 
         bool BitVector::operator!=(BitVector const& other) const {
@@ -166,9 +166,9 @@ namespace storm {
 
             uint64_t mask = 1ull << (63 - (index & mod64mask));
             if (value) {
-                bucketVector[bucket] |= mask;
+                buckets[bucket] |= mask;
             } else {
-                bucketVector[bucket] &= ~mask;
+                buckets[bucket] &= ~mask;
             }
         }
 
@@ -182,7 +182,7 @@ namespace storm {
         bool BitVector::operator[](uint_fast64_t index) const {
             uint64_t bucket = index >> 6;
             uint64_t mask = 1ull << (63 - (index & mod64mask));
-            return (this->bucketVector[bucket] & mask) == mask;
+            return (this->buckets[bucket] & mask) == mask;
         }
 
         bool BitVector::get(uint_fast64_t index) const {
@@ -197,89 +197,76 @@ namespace storm {
                     ++newBucketCount;
                 }
 
-                if (newBucketCount > bucketVector.size()) {
+                if (newBucketCount > this->bucketCount()) {
+                    uint64_t* newBuckets = new uint64_t[newBucketCount];
+                    std::copy_n(buckets, this->bucketCount(), newBuckets);
                     if (init) {
-                        bucketVector.back() |= ((1ull << (64 - (bitCount & mod64mask))) - 1ull);
-                        bucketVector.resize(newBucketCount, -1ull);
+                        newBuckets[this->bucketCount() - 1] |= ((1ull << (64 - (bitCount & mod64mask))) - 1ull);
+                        std::fill_n(newBuckets, newBucketCount - this->bucketCount(), -1ull);
                     } else {
-                        bucketVector.resize(newBucketCount, 0);
+                        std::fill_n(newBuckets, newBucketCount - this->bucketCount(), 0);
                     }
+                    delete buckets;
+                    buckets = newBuckets;
                     bitCount = newLength;
                 } else {
                     // If the underlying storage does not need to grow, we have to insert the missing bits.
                     if (init) {
-                        bucketVector.back() |= ((1ull << (64 - (bitCount & mod64mask))) - 1ull);
+                        buckets[this->bucketCount() - 1] |= ((1ull << (64 - (bitCount & mod64mask))) - 1ull);
                     }
                     bitCount = newLength;
                 }
                 truncateLastBucket();
             } else {
-                bitCount = newLength;
                 uint_fast64_t newBucketCount = newLength >> 6;
                 if ((newLength & mod64mask) != 0) {
                     ++newBucketCount;
                 }
 
-                bucketVector.resize(newBucketCount);
+                // If the number of buckets needs to be reduced, we resize it now. Otherwise, we can just truncate the
+                // last bucket.
+                if (newBucketCount < this->bucketCount()) {
+                    uint64_t* newBuckets = new uint64_t[newBucketCount];
+                    std::copy_n(buckets, newBucketCount, newBuckets);
+                    delete buckets;
+                    buckets = newBuckets;
+                    bitCount = newLength;
+                }
+                bitCount = newLength;
                 truncateLastBucket();
             }
         }
 
         BitVector BitVector::operator&(BitVector const& other) const {
             STORM_LOG_ASSERT(bitCount == other.bitCount, "Length of the bit vectors does not match.");
-
             BitVector result(bitCount);
-            std::vector<uint64_t>::iterator it = result.bucketVector.begin();
-            for (std::vector<uint64_t>::const_iterator it1 = bucketVector.begin(), it2 = other.bucketVector.begin(); it != result.bucketVector.end(); ++it1, ++it2, ++it) {
-                *it = *it1 & *it2;
-            }
-
+            std::transform(this->buckets, this->buckets + this->bucketCount(), other.buckets, result.buckets, [] (uint64_t const& a, uint64_t const& b) { return a & b; });
             return result;
         }
 
         BitVector& BitVector::operator&=(BitVector const& other) {
             STORM_LOG_ASSERT(bitCount == other.bitCount, "Length of the bit vectors does not match.");
-
-            std::vector<uint64_t>::iterator it = bucketVector.begin();
-            for (std::vector<uint64_t>::const_iterator otherIt = other.bucketVector.begin(); it != bucketVector.end() && otherIt != other.bucketVector.end(); ++it, ++otherIt) {
-                *it &= *otherIt;
-            }
-
+            std::transform(this->buckets, this->buckets + this->bucketCount(), other.buckets, this->buckets, [] (uint64_t const& a, uint64_t const& b) { return a & b; });
             return *this;
         }
 
         BitVector BitVector::operator|(BitVector const& other) const {
             STORM_LOG_ASSERT(bitCount == other.bitCount, "Length of the bit vectors does not match.");
-
             BitVector result(bitCount);
-            std::vector<uint64_t>::iterator it = result.bucketVector.begin();
-            for (std::vector<uint64_t>::const_iterator it1 = bucketVector.begin(), it2 = other.bucketVector.begin(); it != result.bucketVector.end(); ++it1, ++it2, ++it) {
-                *it = *it1 | *it2;
-            }
-
+            std::transform(this->buckets, this->buckets + this->bucketCount(), other.buckets, result.buckets, [] (uint64_t const& a, uint64_t const& b) { return a | b; });
             return result;
         }
 
         BitVector& BitVector::operator|=(BitVector const& other) {
             STORM_LOG_ASSERT(bitCount == other.bitCount, "Length of the bit vectors does not match.");
-
-            std::vector<uint64_t>::iterator it = bucketVector.begin();
-            for (std::vector<uint64_t>::const_iterator otherIt = other.bucketVector.begin(); it != bucketVector.end() && otherIt != other.bucketVector.end(); ++it, ++otherIt) {
-                *it |= *otherIt;
-            }
-
+            std::transform(this->buckets, this->buckets + this->bucketCount(), other.buckets, this->buckets, [] (uint64_t const& a, uint64_t const& b) { return a | b; });
             return *this;
         }
 
         BitVector BitVector::operator^(BitVector const& other) const {
             STORM_LOG_ASSERT(bitCount == other.bitCount, "Length of the bit vectors does not match.");
-
             BitVector result(bitCount);
-            std::vector<uint64_t>::iterator it = result.bucketVector.begin();
-            for (std::vector<uint64_t>::const_iterator it1 = bucketVector.begin(), it2 = other.bucketVector.begin(); it != result.bucketVector.end(); ++it1, ++it2, ++it) {
-                *it = *it1 ^ *it2;
-            }
-
+            std::transform(this->buckets, this->buckets + this->bucketCount(), other.buckets, result.buckets, [] (uint64_t const& a, uint64_t const& b) { return a ^ b; });
             result.truncateLastBucket();
             return result;
         }
@@ -314,19 +301,13 @@ namespace storm {
 
         BitVector BitVector::operator~() const {
             BitVector result(this->bitCount);
-            std::vector<uint64_t>::iterator it = result.bucketVector.begin();
-            for (std::vector<uint64_t>::const_iterator it1 = bucketVector.begin(); it != result.bucketVector.end(); ++it1, ++it) {
-                *it = ~(*it1);
-            }
-
+            std::transform(this->buckets, this->buckets + this->bucketCount(), result.buckets, [] (uint64_t const& a) { return ~a; });
             result.truncateLastBucket();
             return result;
         }
 
         void BitVector::complement() {
-            for (auto& element : bucketVector) {
-                element = ~element;
-            }
+            std::transform(this->buckets, this->buckets + this->bucketCount(), this->buckets, [] (uint64_t const& a) { return ~a; });
             truncateLastBucket();
         }
 
@@ -334,11 +315,7 @@ namespace storm {
             STORM_LOG_ASSERT(bitCount == other.bitCount, "Length of the bit vectors does not match.");
 
             BitVector result(bitCount);
-            std::vector<uint64_t>::iterator it = result.bucketVector.begin();
-            for (std::vector<uint64_t>::const_iterator it1 = bucketVector.begin(), it2 = other.bucketVector.begin(); it != result.bucketVector.end(); ++it1, ++it2, ++it) {
-                *it = ~(*it1) | *it2;
-            }
-
+            std::transform(this->buckets, this->buckets + this->bucketCount(), other.buckets, result.buckets, [] (uint64_t const& a, uint64_t const& b) { return (~a | b); });
             result.truncateLastBucket();
             return result;
         }
@@ -346,7 +323,11 @@ namespace storm {
         bool BitVector::isSubsetOf(BitVector const& other) const {
             STORM_LOG_ASSERT(bitCount == other.bitCount, "Length of the bit vectors does not match.");
 
-            for (std::vector<uint64_t>::const_iterator it1 = bucketVector.begin(), it2 = other.bucketVector.begin(); it1 != bucketVector.end(); ++it1, ++it2) {
+            uint64_t const* it1 = buckets;
+            uint64_t const* ite1 = buckets + bucketCount();
+            uint64_t const* it2 = other.buckets;
+            
+            for (; it1 != ite1; ++it1, ++it2) {
                 if ((*it1 & *it2) != *it1) {
                     return false;
                 }
@@ -357,7 +338,11 @@ namespace storm {
         bool BitVector::isDisjointFrom(BitVector const& other) const {
             STORM_LOG_ASSERT(bitCount == other.bitCount, "Length of the bit vectors does not match.");
 
-            for (std::vector<uint64_t>::const_iterator it1 = bucketVector.begin(), it2 = other.bucketVector.begin(); it1 != bucketVector.end(); ++it1, ++it2) {
+            uint64_t const* it1 = buckets;
+            uint64_t const* ite1 = buckets + bucketCount();
+            uint64_t const* it2 = other.buckets;
+            
+            for (; it1 != ite1; ++it1, ++it2) {
                 if ((*it1 & *it2) != 0) {
                     return false;
                 }
@@ -372,9 +357,9 @@ namespace storm {
             // Compute the first bucket that needs to be checked and the number of buckets.
             uint64_t index = bitIndex >> 6;
 
-            std::vector<uint64_t>::const_iterator first1 = bucketVector.begin() + index;
-            std::vector<uint64_t>::const_iterator first2 = other.bucketVector.begin();
-            std::vector<uint64_t>::const_iterator last2 = other.bucketVector.end();
+            uint64_t const* first1 = buckets + index;
+            uint64_t const* first2 = other.buckets;
+            uint64_t const* last2 = other.buckets + other.bucketCount();
 
             for (; first2 != last2; ++first1, ++first2) {
                 if (*first1 != *first2) {
@@ -391,10 +376,10 @@ namespace storm {
             // Compute the first bucket that needs to be checked and the number of buckets.
             uint64_t index = bitIndex >> 6;
 
-            std::vector<uint64_t>::iterator first1 = bucketVector.begin() + index;
-            std::vector<uint64_t>::const_iterator first2 = other.bucketVector.begin();
-            std::vector<uint64_t>::const_iterator last2 = other.bucketVector.end();
-
+            uint64_t* first1 = buckets + index;
+            uint64_t const* first2 = other.buckets;
+            uint64_t const* last2 = other.buckets + other.bucketCount();
+            
             for (; first2 != last2; ++first1, ++first2) {
                 *first1 = *first2;
             }
@@ -406,8 +391,8 @@ namespace storm {
             STORM_LOG_ASSERT(index + numberOfBuckets <= this->bucketCount(), "Argument is out-of-range.");
 
             storm::storage::BitVector result(numberOfBuckets, numberOfBits);
-            std::copy(this->bucketVector.begin() + index, this->bucketVector.begin() + index + numberOfBuckets, result.bucketVector.begin());
-
+            std::copy(this->buckets + index, this->buckets + index + numberOfBuckets, result.buckets);
+            result.truncateLastBucket();
             return result;
         }
 
@@ -425,10 +410,10 @@ namespace storm {
             if (bitIndexInBucket + numberOfBits < 64) {
                 // If the value stops before the end of the bucket, we need to erase some lower bits.
                 mask &= ~((1ull << (64 - (bitIndexInBucket + numberOfBits))) - 1ull);
-                return (bucketVector[bucket] & mask) >> (64 - (bitIndexInBucket + numberOfBits));
+                return (buckets[bucket] & mask) >> (64 - (bitIndexInBucket + numberOfBits));
             } else if (bitIndexInBucket + numberOfBits > 64) {
                 // In this case, the integer "crosses" the bucket line.
-                uint64_t result = (bucketVector[bucket] & mask);
+                uint64_t result = (buckets[bucket] & mask);
                 ++bucket;
 
                 // Compute the remaining number of bits.
@@ -440,13 +425,13 @@ namespace storm {
                 // Strip away everything from the second bucket that is beyond the final index and add it to the
                 // intermediate result.
                 mask = ~((1ull << (64 - numberOfBits)) - 1ull);
-                uint64_t lowerBits = bucketVector[bucket] & mask;
+                uint64_t lowerBits = buckets[bucket] & mask;
                 result |= (lowerBits >> (64 - numberOfBits));
 
                 return result;
             } else {
                 // In this case, it suffices to take the current mask.
-                return bucketVector[bucket] & mask;
+                return buckets[bucket] & mask;
             }
         }
 
@@ -466,10 +451,10 @@ namespace storm {
             if (bitIndexInBucket + numberOfBits < 64) {
                 // If the value stops before the end of the bucket, we need to erase some lower bits.
                 mask &= ~((1ull << (64 - (bitIndexInBucket + numberOfBits))) - 1ull);
-                bucketVector[bucket] = (bucketVector[bucket] & ~mask) | (value << (64 - (bitIndexInBucket + numberOfBits)));
+                buckets[bucket] = (buckets[bucket] & ~mask) | (value << (64 - (bitIndexInBucket + numberOfBits)));
             } else if (bitIndexInBucket + numberOfBits > 64) {
                 // Write the part of the value that falls into the first bucket.
-                bucketVector[bucket] = (bucketVector[bucket] & ~mask) | (value >> (numberOfBits + (bitIndexInBucket - 64)));
+                buckets[bucket] = (buckets[bucket] & ~mask) | (value >> (numberOfBits + (bitIndexInBucket - 64)));
                 ++bucket;
 
                 // Compute the remaining number of bits.
@@ -480,45 +465,41 @@ namespace storm {
 
                 // Put the remaining bits in their place.
                 mask = ((1ull << (64 - numberOfBits)) - 1ull);
-                bucketVector[bucket] = (bucketVector[bucket] & mask) | value;
+                buckets[bucket] = (buckets[bucket] & mask) | value;
             } else {
-                bucketVector[bucket] = (bucketVector[bucket] & ~mask) | value;
+                buckets[bucket] = (buckets[bucket] & ~mask) | value;
             }
         }
 
         bool BitVector::empty() const {
-            for (auto& element : bucketVector) {
-                if (element != 0) {
-                    return false;
-                }
-            }
-            return true;
+            uint64_t* last = buckets + bucketCount();
+            uint64_t* it = std::find(buckets, last, 0);
+            return it != last;
         }
 
         bool BitVector::full() const {
             // Check that all buckets except the last one have all bits set.
-            for (uint_fast64_t index = 0; index < bucketVector.size() - 1; ++index) {
-                if (bucketVector[index] != -1ull) {
+            uint64_t* last = buckets + bucketCount() - 1;
+            for (uint64_t const* it = buckets; it != last; ++it) {
+                if (*it != -1ull) {
                     return false;
                 }
             }
 
             // Now check whether the relevant bits are set in the last bucket.
             uint64_t mask = ~((1ull << (64 - (bitCount & mod64mask))) - 1ull);
-            if ((bucketVector.back() & mask) != mask) {
+            if ((*last & mask) != mask) {
                 return false;
             }
             return true;
         }
 
         void BitVector::clear() {
-            for (auto& element : bucketVector) {
-                element = 0;
-            }
+            std::fill_n(buckets, this->bucketCount(), 0);
         }
 
         uint_fast64_t BitVector::getNumberOfSetBits() const {
-            return getNumberOfSetBitsBeforeIndex(bucketVector.size() << 6);
+            return getNumberOfSetBitsBeforeIndex(bitCount);
         }
 
         uint_fast64_t BitVector::getNumberOfSetBitsBeforeIndex(uint_fast64_t index) const {
@@ -529,14 +510,14 @@ namespace storm {
             for (uint_fast64_t i = 0; i < bucket; ++i) {
                 // Check if we are using g++ or clang++ and, if so, use the built-in function
 #if (defined (__GNUG__) || defined(__clang__))
-                result += __builtin_popcountll(bucketVector[i]);
+                result += __builtin_popcountll(buckets[i]);
 #elif defined WINDOWS
 #include <nmmintrin.h>
                 // If the target machine does not support SSE4, this will fail.
                 result += _mm_popcnt_u64(bucketVector[i]);
 #else
                 uint_fast32_t cnt;
-                uint_fast64_t bitset = bucketVector[i];
+                uint_fast64_t bitset = buckets[i];
                 for (cnt = 0; bitset; cnt++) {
                     bitset &= bitset - 1;
                 }
@@ -548,7 +529,7 @@ namespace storm {
             uint64_t tmp = index & mod64mask;
             if (tmp != 0) {
                 tmp = ~((1ll << (64 - (tmp & mod64mask))) - 1ll);
-                tmp &= bucketVector[bucket];
+                tmp &= buckets[bucket];
                 // Check if we are using g++ or clang++ and, if so, use the built-in function
 #if (defined (__GNUG__) || defined(__clang__))
                 result += __builtin_popcountll(tmp);
@@ -585,19 +566,19 @@ namespace storm {
         }
 
         std::size_t BitVector::getSizeInBytes() const {
-            return sizeof (*this) + sizeof (uint64_t) * bucketVector.size();
+            return sizeof (*this) + sizeof (uint64_t) * bucketCount();
         }
 
         BitVector::const_iterator BitVector::begin() const {
-            return const_iterator(bucketVector.data(), 0, bitCount);
+            return const_iterator(buckets, 0, bitCount);
         }
 
         BitVector::const_iterator BitVector::end() const {
-            return const_iterator(bucketVector.data(), bitCount, bitCount, false);
+            return const_iterator(buckets, bitCount, bitCount, false);
         }
 
         uint_fast64_t BitVector::getNextSetIndex(uint_fast64_t startingIndex) const {
-            return getNextSetIndex(bucketVector.data(), startingIndex, bitCount);
+            return getNextSetIndex(buckets, startingIndex, bitCount);
         }
 
         uint_fast64_t BitVector::getNextSetIndex(uint64_t const* dataPtr, uint_fast64_t startingIndex, uint_fast64_t endIndex) {
@@ -641,17 +622,17 @@ namespace storm {
 
         void BitVector::truncateLastBucket() {
             if ((bitCount & mod64mask) != 0) {
-                bucketVector.back() &= ~((1ll << (64 - (bitCount & mod64mask))) - 1ll);
+                buckets[bucketCount() - 1] &= ~((1ll << (64 - (bitCount & mod64mask))) - 1ll);
             }
         }
 
         size_t BitVector::bucketCount() const {
-            return bucketVector.size();
+            return bitCount >> 6;
         }
 
-        std::ostream& operator<<(std::ostream& out, BitVector const& bitVector) {
-            out << "bit vector(" << bitVector.getNumberOfSetBits() << "/" << bitVector.bitCount << ") [";
-            for (auto index : bitVector) {
+        std::ostream& operator<<(std::ostream& out, BitVector const& bitvector) {
+            out << "bit vector(" << bitvector.getNumberOfSetBits() << "/" << bitvector.bitCount << ") [";
+            for (auto index : bitvector) {
                 out << index << " ";
             }
             out << "]";
@@ -659,13 +640,13 @@ namespace storm {
             return out;
         }
 
-        std::size_t NonZeroBitVectorHash::operator()(storm::storage::BitVector const& bv) const {
-            STORM_LOG_ASSERT(bv.size() > 0, "Cannot hash bit vector of zero size.");
+        std::size_t NonZeroBitVectorHash::operator()(storm::storage::BitVector const& bitvector) const {
+            STORM_LOG_ASSERT(bitvector.size() > 0, "Cannot hash bit vector of zero size.");
             std::size_t result = 0;
 
-            for (uint_fast64_t index = 0; index < bv.bucketVector.size(); ++index) {
+            for (uint_fast64_t index = 0; index < bitvector.bucketCount(); ++index) {
                 result ^= result << 3;
-                result ^= result >> bv.getAsInt(index << 6, 5);
+                result ^= result >> bitvector.getAsInt(index << 6, 5);
             }
 
             // Erase the last bit and add one to definitely make this hash value non-zero.
diff --git a/src/storage/BitVector.h b/src/storage/BitVector.h
index f9c088ceb..a6ad89cc3 100644
--- a/src/storage/BitVector.h
+++ b/src/storage/BitVector.h
@@ -104,6 +104,11 @@ namespace storm {
              */
             BitVector();
             
+            /*
+             * Deconstructs a bit vector by deleting the underlying storage.
+             */
+            ~BitVector();
+            
             /*!
              * Constructs a bit vector which can hold the given number of bits and initializes all bits with the
              * provided truth value.
@@ -494,7 +499,7 @@ namespace storm {
             uint_fast64_t bitCount;
             
             // The underlying storage of 64-bit buckets for all bits of this bit vector.
-            std::vector<uint64_t> bucketVector;
+            uint64_t* buckets;
             
             // A bit mask that can be used to reduce a modulo 64 operation to a logical "and".
             static const uint_fast64_t mod64mask = (1 << 6) - 1;