diff --git a/src/storm/modelchecker/csl/helper/SparseCtmcCslHelper.cpp b/src/storm/modelchecker/csl/helper/SparseCtmcCslHelper.cpp
index 41b768058..2d270c720 100644
--- a/src/storm/modelchecker/csl/helper/SparseCtmcCslHelper.cpp
+++ b/src/storm/modelchecker/csl/helper/SparseCtmcCslHelper.cpp
@@ -632,7 +632,9 @@ namespace storm {
                 }
                 
                 // Use Fox-Glynn to get the truncation points and the weights.
-                std::tuple<uint_fast64_t, uint_fast64_t, ValueType, std::vector<ValueType>> foxGlynnResult = storm::utility::numerical::getFoxGlynnCutoff(lambda, 1e+300, storm::settings::getModule<storm::settings::modules::GeneralSettings>().getPrecision() / 8.0);
+//                std::tuple<uint_fast64_t, uint_fast64_t, ValueType, std::vector<ValueType>> foxGlynnResult = storm::utility::numerical::getFoxGlynnCutoff(lambda, 1e+300, storm::settings::getModule<storm::settings::modules::GeneralSettings>().getPrecision() / 8.0);
+                
+                std::tuple<uint_fast64_t, uint_fast64_t, ValueType, std::vector<ValueType>> foxGlynnResult = storm::utility::numerical::foxGlynn(lambda, storm::settings::getModule<storm::settings::modules::GeneralSettings>().getPrecision() / 8.0);
                 STORM_LOG_DEBUG("Fox-Glynn cutoff points: left=" << std::get<0>(foxGlynnResult) << ", right=" << std::get<1>(foxGlynnResult));
                 
                 // Scale the weights so they add up to one.
diff --git a/src/storm/utility/numerical.h b/src/storm/utility/numerical.h
index bf9a533b4..78be4f6d6 100644
--- a/src/storm/utility/numerical.h
+++ b/src/storm/utility/numerical.h
@@ -2,10 +2,13 @@
 #include <tuple>
 #include <cmath>
 
+#include <boost/math/constants/constants.hpp>
+
 #include "storm/utility/macros.h"
 #include "storm/utility/constants.h"
 #include "storm/exceptions/InvalidArgumentException.h"
 #include "storm/exceptions/OutOfRangeException.h"
+#include "storm/exceptions/PrecisionExceededException.h"
 
 namespace storm {
     namespace utility {
@@ -125,7 +128,152 @@ namespace storm {
                     
                     return std::make_tuple(leftTruncationPoint, rightTruncationPoint, totalWeight, weights);
                 }
+ 
+            }
+
+            template<typename ValueType>
+            std::tuple<uint64_t, uint64_t, ValueType> foxGlynnFinder(ValueType lambda, ValueType underflow, ValueType overflow, ValueType accuracy) {
+                uint64_t m = static_cast<uint64_t>(lambda);
+                int64_t L = 0;
+                uint64_t R = 0;
+                uint64_t k = 4;
+                
+                ValueType pi = boost::math::constants::pi<ValueType>();
+
+                if (m < 25) {
+                    STORM_LOG_THROW(std::exp(-static_cast<double>(m)) >= underflow, storm::exceptions::PrecisionExceededException, "Underflow in Fox-Glynn.");
+                    L = 0;
+                } else {
+                    ValueType rhsb = (accuracy / 2) * std::sqrt(2 * pi) / ((1 + 1/lambda) * std::exp(1 / (8 * lambda)));
+                    k = 4;
+                    
+                    do {
+                        L = m - std::ceil(k * std::sqrt(lambda) + 0.5);
+                        if (L <= 0) {
+                            break;
+                        }
+                        if (std::exp(-(k*k) / 2) / k <= rhsb) {
+                            break;
+                        }
+                        ++k;
+                    } while (true);
+                }
+                
+                uint64_t mmax;
+                uint64_t maxr;
+                ValueType r2l;
+                ValueType rhsa;
+
+                if (m < 400) {
+                    mmax = 400;
+                    r2l = std::sqrt(2 * 400);
+                    rhsa = (accuracy / 2) * std::sqrt(2 * pi) / ((1 + 1.0/400) * std::sqrt(2.0) * std::exp(1.0 / 16));
+                    maxr = 400 + static_cast<uint64_t>(std::ceil((400 + 1) / 2.0));
+                } else {
+                    mmax = m;
+                    r2l = std::sqrt(2 * lambda);
+                    rhsa = (accuracy / 2) * std::sqrt(2 * pi) / ((1 + 1.0/lambda) * std::sqrt(2.0) * std::exp(1.0 / 16));
+                    maxr = m + static_cast<uint64_t>(std::ceil((lambda + 1) / 2.0));
+                }
+                
+                k = 4;
+                do {
+                    R = maxr + static_cast<uint64_t>(std::ceil(k * r2l + 0.5));
+                    STORM_LOG_THROW(R <= maxr, storm::exceptions::PrecisionExceededException, "Fox-Glynn: cannot bound right tail.");
+                    ValueType dkl = storm::utility::one<ValueType>();
+                    if (dkl * std::exp(-(k * k / 2.0)) / k <= rhsa) {
+                        break;
+                    }
+                    ++k;
+                } while (true);
+                
+                ValueType wm = std::pow<ValueType>(10, -10) * overflow * (R - L + 1);
+
+                if (m >= 25) {
+                    ValueType lcm = -1 - (1.0/12*25) - std::log(2 * pi) - 0.5 * std::log(static_cast<ValueType>(m));
+                    uint64_t i = m - L;
+                    ValueType llb;
+                    if (i <= static_cast<uint64_t>(L)) {
+                        llb = -static_cast<ValueType>(i)*(i + 1)*(((2 * i + 1) / (6 * lambda)) + 0.5) * (1 / lambda) + lcm;
+                    } else {
+                        llb = std::max(i * std::log(1 - (i / (m + 1.0))), -lambda);
+                    }
+                    STORM_LOG_THROW(llb >= std::log(underflow) - std::log(wm), storm::exceptions::PrecisionExceededException, "Fox-Glynn: Underflow at lower bound.");
+                    if (m >= 400) {
+                        llb = lcm - std::pow(R - m + 1, 2) / (2 * lambda);
+                        STORM_LOG_THROW(llb >= std::log(underflow) - std::log(wm), storm::exceptions::PrecisionExceededException, "Fox-Glynn: Underflow at upper bound.");
+                    }
+                }
+                
+                return std::make_tuple(static_cast<uint64_t>(L), R, wm);
+            }
+            
+            template<typename ValueType>
+            std::tuple<uint64_t, uint64_t, ValueType, std::vector<ValueType>> foxGlynnWeighter(ValueType lambda, ValueType accuracy) {
+                
+                STORM_LOG_THROW(lambda > storm::utility::zero<ValueType>(), storm::exceptions::InvalidArgumentException, "Fox-Glynn: Lambda must be positive.");
+                
+                ValueType underflow = std::numeric_limits<ValueType>::min();
+                ValueType overflow = std::numeric_limits<ValueType>::max();
+
+                // Initialize.
+                uint64_t m = static_cast<uint64_t>(lambda);
+                std::tuple<uint64_t, uint64_t, ValueType> finderResult = foxGlynnFinder(lambda, underflow, overflow, accuracy);
+                uint64_t L = std::get<0>(finderResult);
+                uint64_t R = std::get<1>(finderResult);
+                ValueType wm = std::get<2>(finderResult);
+                std::vector<ValueType> w(R - L + 1);
+                w.back() = wm;
+
+                // Down.
+                for (uint64_t j = m; j > L; --j) {
+                    w[j - L - 1] = (j / lambda) * w[j - L];
+                }
+
+                // Up.
+                if (lambda < 400) {
+                    // Special.
+                    STORM_LOG_THROW(R <= 600, storm::exceptions::PrecisionExceededException, "Potential underflow in Fox-Glynn.");
+                    for (uint64_t j = m; j < R; ++j) {
+                        auto q = lambda / (j + 1);
+                        if (w[j - L] > underflow / q) {
+                            w[j - L + 1] = q * w[j - L];
+                        } else {
+                            R = j;
+                            break;
+                        }
+                    }
+                } else {
+                    for (uint64_t j = m; j < R; ++j) {
+                        w[j - L + 1] = (lambda / (j + 1)) * w[j - L];
+                    }
+                }
+
+                // Compute W.
+                ValueType W = storm::utility::zero<ValueType>();
+                uint64_t s = L;
+                uint64_t t = R;
+                
+                while (s < t) {
+                    if (w[s - L] <= w[t - L]) {
+                        W += w[s - L];
+                        ++s;
+                    } else {
+                        W += w[t - L];
+                        --t;
+                    }
+                }
+                W += w[s - L];
+                
+                // Return.
+                return std::make_tuple(L, R, W, w);
             }
+            
+            template<typename ValueType>
+            std::tuple<uint_fast64_t, uint_fast64_t, ValueType, std::vector<ValueType>> foxGlynn(ValueType lambda, ValueType accuracy) {
+                return foxGlynnWeighter(lambda, accuracy);
+            }
+                
         }
     }
 }