diff --git a/Statistics/Distribution/Poisson.hs b/Statistics/Distribution/Poisson.hs
--- a/Statistics/Distribution/Poisson.hs
+++ b/Statistics/Distribution/Poisson.hs
@@ -25,7 +25,7 @@
 import qualified Data.Vector.Unboxed as U
 import qualified Statistics.Distribution as D
 import Statistics.Constants (m_huge)
-import Statistics.Math (logGamma)
+import Statistics.Math (factorial, logGamma)
 
 newtype PoissonDistribution = PD {
       pdLambda :: Double
@@ -49,7 +49,13 @@
 {-# INLINE fromLambda #-}
 
 density :: PoissonDistribution -> Double -> Double
-density (PD l) x = exp (x * log l - l - logGamma x)
+density (PD l) x
+    | x < 0                   = 0
+    | l >= 100 && x >= l * 10 = 0
+    | l >= 3 && x >= l * 100  = 0
+    | x >= max 1 l * 200      = 0
+    | l < 20 && x <= 100      = exp (-l) * l ** x / factorial (floor x)
+    | otherwise               = x * log l - logGamma (x + 1) - l
 {-# INLINE density #-}
 
 cumulative :: PoissonDistribution -> Double -> Double
diff --git a/Statistics/Math.hs b/Statistics/Math.hs
--- a/Statistics/Math.hs
+++ b/Statistics/Math.hs
@@ -33,19 +33,18 @@
 import Statistics.Distribution.Normal (standard)
 import qualified Data.Vector.Unboxed as U
 
-data C = C {-# UNPACK #-} !Double {-# UNPACK #-} !Double {-# UNPACK #-} !Double
+data C = C {-# UNPACK #-} !Double {-# UNPACK #-} !Double
 
 -- | Evaluate a series of Chebyshev polynomials. Uses Clenshaw's
 -- algorithm.
 chebyshev :: Double             -- ^ Parameter of each function.
           -> U.Vector Double    -- ^ Coefficients of each polynomial
-          -- term, in increasing order.
+                                --   term, in increasing order.
           -> Double
-chebyshev x a = fini . U.foldl step (C 0 0 0) .
-                U.enumFromThenTo (U.length a - 1) (-1) $ 0
-    where step (C u v w) k = C (x2 * v - w + (a ! k)) u v
-          fini (C u _ w)   = (u - w) / 2
-          x2               = x * 2
+chebyshev x a = fini . U.foldl step (C 0 0) $ U.enumFromStepN (U.length a - 1) (-1) (U.length a - 1)
+    where step (C b1 b2) k = C ((a ! k) + x2 * b1 - b2) b1
+          fini (C b1 b2)   = (a ! 0) + x * b1 - b2
+          x2                 = x * 2
 
 -- | The binomial coefficient.
 --
diff --git a/Statistics/Sample.hs b/Statistics/Sample.hs
--- a/Statistics/Sample.hs
+++ b/Statistics/Sample.hs
@@ -78,7 +78,8 @@
     where
       fini (V a _) = a
       go (V m w) (x,xw) = V m' w'
-          where m' = m + xw * (x - m) / w'
+          where m' | w' == 0   = 0
+                   | otherwise = m + xw * (x - m) / w'
                 w' = w + xw
 {-# INLINE meanWeighted #-}
 
@@ -226,7 +227,7 @@
 {-# INLINE varianceUnbiased #-}
 
 -- | Standard deviation.  This is simply the square root of the
--- maximum likelihood estimate of the variance.
+-- unbiased estimate of the variance.
 stdDev :: Sample -> Double
 stdDev = sqrt . varianceUnbiased
 
diff --git a/statistics.cabal b/statistics.cabal
--- a/statistics.cabal
+++ b/statistics.cabal
@@ -1,5 +1,5 @@
 name:           statistics
-version:        0.5.1.0
+version:        0.5.1.1
 synopsis:       A library of statistical types, data, and functions
 description:
   This library provides a number of common functions and types useful