diff --git a/ChangeLog b/ChangeLog
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,4 +1,25 @@
- Changes in 0.10.2.0
+Changes in 0.10.5.0
+
+  * Enthropy type class for distributions is added.
+
+  * Probability and probability density of distribution is given in
+    log domain too.
+
+Changes in 0.10.4.1
+
+  * Bugfix for GHC < 7.6
+
+Changes in 0.10.4.0
+
+  * Support for versions of GHC older than 7.2 is discontinued.
+
+  * All datatypes now support 'Data.Binary' and 'GHC.Generics'.
+
+Changes in 0.10.3.0
+
+  * Bug fixes
+
+Changes in 0.10.2.0
 
   * Bugs in DCT and IDCT are fixed.
 
diff --git a/Statistics/Distribution.hs b/Statistics/Distribution.hs
--- a/Statistics/Distribution.hs
+++ b/Statistics/Distribution.hs
@@ -21,6 +21,8 @@
     , Mean(..)
     , MaybeVariance(..)
     , Variance(..)
+    , MaybeEntropy(..)
+    , Entropy(..)
       -- ** Random number generation
     , ContGen(..)
     , DiscreteGen(..)
@@ -65,20 +67,36 @@
 class Distribution  d => DiscreteDistr d where
     -- | Probability of n-th outcome.
     probability :: d -> Int -> Double
+    probability d = exp . logProbability d
+    {-# INLINE probability #-}
 
+    -- | Logarithm of probability of n-th outcome
+    logProbability :: d -> Int -> Double
+    logProbability d = log . probability d
+    {-# INLINE logProbability #-}
 
--- | Continuous probability distributuion
+
+-- | Continuous probability distributuion.
+--
+--   Minimal complete definition is 'quantile' and either 'density' or
+--   'logDensity'.
 class Distribution d => ContDistr d where
     -- | Probability density function. Probability that random
     -- variable /X/ lies in the infinitesimal interval
     -- [/x/,/x+/&#948;/x/) equal to /density(x)/&#8901;&#948;/x/
     density :: d -> Double -> Double
+    density d = exp . logDensity d
+    {-# INLINE density #-}
 
     -- | Inverse of the cumulative distribution function. The value
     -- /x/ for which P(/X/&#8804;/x/) = /p/. If probability is outside
     -- of [0,1] range function should call 'error'
     quantile :: d -> Double -> Double
 
+    -- | Natural logarithm of density.
+    logDensity :: d -> Double -> Double
+    logDensity d = log . density d
+    {-# INLINE logDensity #-}
 
 
 -- | Type class for distributions with mean. 'maybeMean' should return
@@ -116,6 +134,22 @@
     stdDev   :: d -> Double
     stdDev = sqrt . variance
 
+-- | Type class for distributions with entropy, meaning Shannon entropy
+--   in the case of a discrete distribution, or differential entropy in the
+--   case of a continuous one.  'maybeEntropy' should return 'Nothing' if
+--   entropy is undefined for the chosen parameter values.
+class (Distribution d) => MaybeEntropy d where
+  -- | Returns the entropy of a distribution, in nats, if such is defined.
+  maybeEntropy :: d -> Maybe Double
+  
+-- | Type class for distributions with entropy, meaning Shannon
+--   entropy in the case of a discrete distribution, or differential
+--   entropy in the case of a continuous one.  If the distribution has
+--   well-defined entropy for all valid parameter values then it
+--   should be an instance of this type class.
+class (MaybeEntropy d) => Entropy d where
+  -- | Returns the entropy of a distribution, in nats.
+  entropy :: d -> Double
 
 -- | Generate discrete random variates which have given
 --   distribution.
diff --git a/Statistics/Distribution/Beta.hs b/Statistics/Distribution/Beta.hs
--- a/Statistics/Distribution/Beta.hs
+++ b/Statistics/Distribution/Beta.hs
@@ -23,7 +23,8 @@
 import Data.Binary (Binary)
 import Data.Data (Data, Typeable)
 import GHC.Generics (Generic)
-import Numeric.SpecFunctions           (incompleteBeta, invIncompleteBeta, logBeta)
+import Numeric.SpecFunctions (
+  incompleteBeta, invIncompleteBeta, logBeta, digamma)
 import Numeric.MathFunctions.Constants (m_NaN)
 import qualified Statistics.Distribution as D
 
@@ -81,6 +82,18 @@
 instance D.MaybeVariance BetaDistribution where
   maybeVariance = Just . D.variance
   {-# INLINE maybeVariance #-}
+
+instance D.Entropy BetaDistribution where
+  entropy (BD a b) =
+    logBeta a b 
+    - (a-1) * digamma a
+    - (b-1) * digamma b
+    + (a+b-2) * digamma (a+b)
+  {-# INLINE entropy #-}
+    
+instance D.MaybeEntropy BetaDistribution where
+  maybeEntropy = Just . D.entropy
+  {-# INLINE maybeEntropy #-}
 
 instance D.ContDistr BetaDistribution where
   density (BD a b) x
diff --git a/Statistics/Distribution/Binomial.hs b/Statistics/Distribution/Binomial.hs
--- a/Statistics/Distribution/Binomial.hs
+++ b/Statistics/Distribution/Binomial.hs
@@ -27,7 +27,9 @@
 import Data.Data (Data, Typeable)
 import GHC.Generics (Generic)
 import qualified Statistics.Distribution as D
-import Numeric.SpecFunctions (choose)
+import qualified Statistics.Distribution.Poisson.Internal as I
+import Numeric.SpecFunctions (choose,incompleteBeta)
+import Numeric.MathFunctions.Constants (m_epsilon)
 
 
 -- | The binomial distribution.
@@ -59,7 +61,15 @@
     maybeStdDev   = Just . D.stdDev
     maybeVariance = Just . D.variance
 
+instance D.Entropy BinomialDistribution where
+  entropy (BD n p)
+    | n == 0 = 0
+    | n <= 100 = directEntropy (BD n p)
+    | otherwise = I.poissonEntropy (fromIntegral n * p)
 
+instance D.MaybeEntropy BinomialDistribution where
+  maybeEntropy = Just . D.entropy
+
 -- This could be slow for big n
 probability :: BinomialDistribution -> Int -> Double
 probability (BD n p) k
@@ -70,15 +80,13 @@
 
 -- Summation from different sides required to reduce roundoff errors
 cumulative :: BinomialDistribution -> Double -> Double
-cumulative d@(BD n _) x
+cumulative (BD n p) x
   | isNaN x      = error "Statistics.Distribution.Binomial.cumulative: NaN input"
   | isInfinite x = if x > 0 then 1 else 0
   | k <  0       = 0
   | k >= n       = 1
-  | k <  m       = D.sumProbabilities d 0 k
-  | otherwise    = 1 - D.sumProbabilities d (k+1) n
+  | otherwise    = incompleteBeta (fromIntegral (n-k)) (fromIntegral (k+1)) (1 - p)
   where
-    m = floor (mean d)
     k = floor x
 {-# INLINE cumulative #-}
 
@@ -89,6 +97,13 @@
 variance :: BinomialDistribution -> Double
 variance (BD n p) = fromIntegral n * p * (1 - p)
 {-# INLINE variance #-}
+
+directEntropy :: BinomialDistribution -> Double
+directEntropy d@(BD n _) =   
+  negate . sum $
+  takeWhile (< negate m_epsilon) $
+  dropWhile (not . (< negate m_epsilon)) $
+  [ let x = probability d k in x * log x | k <- [0..n]]
 
 -- | Construct binomial distribution. Number of trials must be
 --   non-negative and probability must be in [0,1] range
diff --git a/Statistics/Distribution/CauchyLorentz.hs b/Statistics/Distribution/CauchyLorentz.hs
--- a/Statistics/Distribution/CauchyLorentz.hs
+++ b/Statistics/Distribution/CauchyLorentz.hs
@@ -69,3 +69,9 @@
 
 instance D.ContGen CauchyDistribution where
   genContVar = D.genContinous
+
+instance D.Entropy CauchyDistribution where
+  entropy (CD _ s) = log s + log (4*pi)
+    
+instance D.MaybeEntropy CauchyDistribution where
+  maybeEntropy = Just . D.entropy
diff --git a/Statistics/Distribution/ChiSquared.hs b/Statistics/Distribution/ChiSquared.hs
--- a/Statistics/Distribution/ChiSquared.hs
+++ b/Statistics/Distribution/ChiSquared.hs
@@ -21,7 +21,8 @@
 import Data.Binary (Binary)
 import Data.Data (Data, Typeable)
 import GHC.Generics (Generic)
-import Numeric.SpecFunctions (incompleteGamma,invIncompleteGamma,logGamma)
+import Numeric.SpecFunctions (
+  incompleteGamma,invIncompleteGamma,logGamma,digamma)
 
 import qualified Statistics.Distribution         as D
 import qualified System.Random.MWC.Distributions as MWC
@@ -68,6 +69,17 @@
 instance D.MaybeVariance ChiSquared where
     maybeStdDev   = Just . D.stdDev
     maybeVariance = Just . D.variance
+    
+instance D.Entropy ChiSquared where
+  entropy (ChiSquared ndf) =
+    let kHalf = 0.5 * fromIntegral ndf in
+    kHalf 
+    + log 2 
+    + logGamma kHalf
+    + (1-kHalf) * digamma kHalf
+
+instance D.MaybeEntropy ChiSquared where
+  maybeEntropy = Just . D.entropy
 
 instance D.ContGen ChiSquared where
     genContVar (ChiSquared n) = MWC.chiSquare n
diff --git a/Statistics/Distribution/Exponential.hs b/Statistics/Distribution/Exponential.hs
--- a/Statistics/Distribution/Exponential.hs
+++ b/Statistics/Distribution/Exponential.hs
@@ -26,11 +26,13 @@
 import Data.Binary (Binary)
 import Data.Data (Data, Typeable)
 import GHC.Generics (Generic)
+import Numeric.MathFunctions.Constants (m_neg_inf)
 import qualified Statistics.Distribution         as D
 import qualified Statistics.Sample               as S
 import qualified System.Random.MWC.Distributions as MWC
 import Statistics.Types (Sample)
 
+
 newtype ExponentialDistribution = ED {
       edLambda :: Double
     } deriving (Eq, Read, Show, Typeable, Data, Generic)
@@ -42,7 +44,12 @@
     complCumulative = complCumulative
 
 instance D.ContDistr ExponentialDistribution where
-    density  = density
+    density (ED l) x
+      | x < 0     = 0
+      | otherwise = l * exp (-l * x)
+    logDensity (ED l) x
+      | x < 0     = m_neg_inf
+      | otherwise = log l + (-l * x)
     quantile = quantile
 
 instance D.Mean ExponentialDistribution where
@@ -60,6 +67,12 @@
     maybeStdDev   = Just . D.stdDev
     maybeVariance = Just . D.variance
 
+instance D.Entropy ExponentialDistribution where
+  entropy (ED l) = 1 - log l
+
+instance D.MaybeEntropy ExponentialDistribution where
+  maybeEntropy = Just . D.entropy
+
 instance D.ContGen ExponentialDistribution where
   genContVar = MWC.exponential . edLambda
 
@@ -73,10 +86,6 @@
                          | otherwise = exp (-l * x)
 {-# INLINE complCumulative #-}
 
-density :: ExponentialDistribution -> Double -> Double
-density (ED l) x | x < 0     = 0
-                 | otherwise = l * exp (-l * x)
-{-# INLINE density #-}
 
 quantile :: ExponentialDistribution -> Double -> Double
 quantile (ED l) p
diff --git a/Statistics/Distribution/FDistribution.hs b/Statistics/Distribution/FDistribution.hs
--- a/Statistics/Distribution/FDistribution.hs
+++ b/Statistics/Distribution/FDistribution.hs
@@ -18,9 +18,11 @@
 
 import Data.Binary (Binary)
 import Data.Data (Data, Typeable)
+import Numeric.MathFunctions.Constants (m_neg_inf)
 import GHC.Generics (Generic)
 import qualified Statistics.Distribution as D
-import Numeric.SpecFunctions (logBeta, incompleteBeta, invIncompleteBeta)
+import Numeric.SpecFunctions (
+  logBeta, incompleteBeta, invIncompleteBeta, digamma)
 
 
 
@@ -47,19 +49,23 @@
   cumulative = cumulative
 
 instance D.ContDistr FDistribution where
-  density  = density
+  density d x
+    | x <= 0    = 0
+    | otherwise = exp $ logDensity d x
+  logDensity d x
+    | x <= 0    = m_neg_inf
+    | otherwise = logDensity d x
   quantile = quantile
 
 cumulative :: FDistribution -> Double -> Double
 cumulative (F n m _) x
   | x <= 0       = 0
   | isInfinite x = 1            -- Only matches +∞
-  | x > 0        = let y = n*x in incompleteBeta (0.5 * n) (0.5 * m) (y / (m + y))
+  | otherwise    = let y = n*x in incompleteBeta (0.5 * n) (0.5 * m) (y / (m + y))
 
-density :: FDistribution -> Double -> Double
-density (F n m fac) x
-  | x > 0     = exp $ fac + log x * (0.5 * n - 1) - log(m + n*x) * 0.5 * (n + m)
-  | otherwise = 0
+logDensity :: FDistribution -> Double -> Double
+logDensity (F n m fac) x
+  = fac + log x * (0.5 * n - 1) - log(m + n*x) * 0.5 * (n + m)
 
 quantile :: FDistribution -> Double -> Double
 quantile (F n m _) p
@@ -78,6 +84,19 @@
   maybeStdDev (F n m _)
     | m > 4     = Just $ 2 * sqr m * (m + n - 2) / (n * sqr (m - 2) * (m - 4))
     | otherwise = Nothing
+
+instance D.Entropy FDistribution where
+  entropy (F n m _) =
+    let nHalf = 0.5 * n
+        mHalf = 0.5 * m in
+    log (n/m) 
+    + logBeta nHalf mHalf
+    + (1 - nHalf) * digamma nHalf 
+    - (1 + mHalf) * digamma mHalf
+    + (nHalf + mHalf) * digamma (nHalf + mHalf)
+
+instance D.MaybeEntropy FDistribution where
+  maybeEntropy = Just . D.entropy
 
 instance D.ContGen FDistribution where
   genContVar = D.genContinous
diff --git a/Statistics/Distribution/Gamma.hs b/Statistics/Distribution/Gamma.hs
--- a/Statistics/Distribution/Gamma.hs
+++ b/Statistics/Distribution/Gamma.hs
@@ -28,8 +28,9 @@
 import Data.Binary (Binary)
 import Data.Data (Data, Typeable)
 import GHC.Generics (Generic)
-import Numeric.MathFunctions.Constants (m_pos_inf, m_NaN)
-import Numeric.SpecFunctions           (incompleteGamma, invIncompleteGamma)
+import Numeric.MathFunctions.Constants (m_pos_inf, m_NaN, m_neg_inf)
+import Numeric.SpecFunctions (
+  incompleteGamma, invIncompleteGamma, logGamma, digamma)
 import Statistics.Distribution.Poisson.Internal  as Poisson
 import qualified Statistics.Distribution         as D
 import qualified System.Random.MWC.Distributions as MWC
@@ -67,6 +68,9 @@
 
 instance D.ContDistr GammaDistribution where
     density    = density
+    logDensity (GD k theta) x
+      | x <= 0    = m_neg_inf
+      | otherwise = log x * (k - 1) - (x / theta) - logGamma k - log theta * k
     quantile   = quantile
 
 instance D.Variance GammaDistribution where
@@ -83,6 +87,16 @@
 instance D.MaybeVariance GammaDistribution where
     maybeStdDev   = Just . D.stdDev
     maybeVariance = Just . D.variance
+
+instance D.MaybeEntropy GammaDistribution where
+  maybeEntropy (GD a l)
+    | a > 0 && l > 0 = 
+      Just $
+      a 
+      + log l 
+      + logGamma a 
+      + (1-a) * digamma a
+    | otherwise = Nothing
 
 instance D.ContGen GammaDistribution where
     genContVar (GD a l) = MWC.gamma a l
diff --git a/Statistics/Distribution/Geometric.hs b/Statistics/Distribution/Geometric.hs
--- a/Statistics/Distribution/Geometric.hs
+++ b/Statistics/Distribution/Geometric.hs
@@ -8,29 +8,39 @@
 -- Stability   : experimental
 -- Portability : portable
 --
--- The Geometric distribution. This is the probability distribution of
--- the number of Bernoulli trials needed to get one success, supported
--- on the set [1,2..].
+-- The Geometric distribution. There are two variants of
+-- distribution. First is the probability distribution of the number
+-- of Bernoulli trials needed to get one success, supported on the set
+-- [1,2..] ('GeometricDistribution'). Sometimes it's referred to as
+-- the /shifted/ geometric distribution to distinguish from another
+-- one.
 --
--- This distribution is sometimes referred to as the /shifted/
--- geometric distribution, to distinguish it from a variant measuring
--- the number of failures before the first success, defined over the
--- set [0,1..].
-
+-- Second variant is probability distribution of the number of
+-- failures before first success, defined over the set [0,1..]
+-- ('GeometricDistribution0').
 module Statistics.Distribution.Geometric
     (
       GeometricDistribution
+    , GeometricDistribution0
     -- * Constructors
     , geometric
+    , geometric0
     -- ** Accessors
     , gdSuccess
+    , gdSuccess0
     ) where
 
-import Data.Binary (Binary)
-import Data.Data (Data, Typeable)
-import GHC.Generics (Generic)
-import qualified Statistics.Distribution as D
+import Control.Monad  (liftM)
+import Data.Binary    (Binary)
+import Data.Data      (Data, Typeable)
+import GHC.Generics   (Generic)
+import Numeric.MathFunctions.Constants(m_pos_inf,m_neg_inf)
+import qualified Statistics.Distribution         as D
+import qualified System.Random.MWC.Distributions as MWC
 
+----------------------------------------------------------------
+-- Distribution over [1..]
+
 newtype GeometricDistribution = GD {
       gdSuccess :: Double
     } deriving (Eq, Read, Show, Typeable, Data, Generic)
@@ -41,8 +51,14 @@
     cumulative = cumulative
 
 instance D.DiscreteDistr GeometricDistribution where
-    probability = probability
+    probability (GD s) n
+      | n < 1     = 0
+      | otherwise = s * (1-s) ** (fromIntegral n - 1)
+    logProbability (GD s) n
+       | n < 1     = m_neg_inf
+       | otherwise = log s + log (1-s) * (fromIntegral n - 1)
 
+
 instance D.Mean GeometricDistribution where
     mean (GD s) = 1 / s
     {-# INLINE mean #-}
@@ -58,7 +74,22 @@
     maybeStdDev   = Just . D.stdDev
     maybeVariance = Just . D.variance
 
+instance D.Entropy GeometricDistribution where
+  entropy (GD s)
+    | s == 0 = m_pos_inf
+    | s == 1 = 0
+    | otherwise = negate $ (s * log s + (1-s) * log (1-s)) / s
 
+instance D.MaybeEntropy GeometricDistribution where
+  maybeEntropy = Just . D.entropy
+
+instance D.DiscreteGen GeometricDistribution where
+  genDiscreteVar (GD s) g = MWC.geometric1 s g
+  {-# INLINE genDiscreteVar #-}
+instance D.ContGen GeometricDistribution where
+  genContVar d g = fromIntegral `liftM` D.genDiscreteVar d g
+  {-# INLINE genContVar #-}
+
 -- | Create geometric distribution.
 geometric :: Double                -- ^ Success rate
           -> GeometricDistribution
@@ -68,11 +99,6 @@
     error $ "Statistics.Distribution.Geometric.geometric: probability must be in [0,1] range. Got " ++ show x
 {-# INLINE geometric #-}
 
-probability :: GeometricDistribution -> Int -> Double
-probability (GD s) n | n < 1     = 0
-                     | otherwise = s * (1-s) ** (fromIntegral n - 1)
-{-# INLINE probability #-}
-
 cumulative :: GeometricDistribution -> Double -> Double
 cumulative (GD s) x
   | x < 1        = 0
@@ -80,3 +106,57 @@
   | isNaN      x = error "Statistics.Distribution.Geometric.cumulative: NaN input"
   | otherwise    = 1 - (1-s) ^ (floor x :: Int)
 {-# INLINE cumulative #-}
+
+
+----------------------------------------------------------------
+-- Distribution over [0..]
+
+newtype GeometricDistribution0 = GD0 {
+      gdSuccess0 :: Double
+    } deriving (Eq, Read, Show, Typeable, Data, Generic)
+
+instance Binary GeometricDistribution0
+
+instance D.Distribution GeometricDistribution0 where
+    cumulative (GD0 s) x = cumulative (GD s) (x + 1)
+
+instance D.DiscreteDistr GeometricDistribution0 where
+    probability    (GD0 s) n = D.probability    (GD s) (n + 1)
+    logProbability (GD0 s) n = D.logProbability (GD s) (n + 1)
+
+instance D.Mean GeometricDistribution0 where
+    mean (GD0 s) = 1 / s - 1
+    {-# INLINE mean #-}
+
+instance D.Variance GeometricDistribution0 where
+    variance (GD0 s) = D.variance (GD s)
+    {-# INLINE variance #-}
+
+instance D.MaybeMean GeometricDistribution0 where
+    maybeMean = Just . D.mean
+
+instance D.MaybeVariance GeometricDistribution0 where
+    maybeStdDev   = Just . D.stdDev
+    maybeVariance = Just . D.variance
+
+instance D.Entropy GeometricDistribution0 where
+  entropy (GD0 s) = D.entropy (GD s)
+
+instance D.MaybeEntropy GeometricDistribution0 where
+  maybeEntropy = Just . D.entropy
+
+instance D.DiscreteGen GeometricDistribution0 where
+  genDiscreteVar (GD0 s) g = MWC.geometric0 s g
+  {-# INLINE genDiscreteVar #-}
+instance D.ContGen GeometricDistribution0 where
+  genContVar d g = fromIntegral `liftM` D.genDiscreteVar d g
+  {-# INLINE genContVar #-}
+
+-- | Create geometric distribution.
+geometric0 :: Double                -- ^ Success rate
+           -> GeometricDistribution0
+geometric0 x
+  | x >= 0 && x <= 1 = GD0 x
+  | otherwise        =
+    error $ "Statistics.Distribution.Geometric.geometric: probability must be in [0,1] range. Got " ++ show x
+{-# INLINE geometric0 #-}
diff --git a/Statistics/Distribution/Hypergeometric.hs b/Statistics/Distribution/Hypergeometric.hs
--- a/Statistics/Distribution/Hypergeometric.hs
+++ b/Statistics/Distribution/Hypergeometric.hs
@@ -30,6 +30,7 @@
 import Data.Binary (Binary)
 import Data.Data (Data, Typeable)
 import GHC.Generics (Generic)
+import Numeric.MathFunctions.Constants (m_epsilon)
 import Numeric.SpecFunctions (choose)
 import qualified Statistics.Distribution as D
 
@@ -60,7 +61,11 @@
     maybeStdDev   = Just . D.stdDev
     maybeVariance = Just . D.variance
 
-
+instance D.Entropy HypergeometricDistribution where
+  entropy = directEntropy
+  
+instance D.MaybeEntropy HypergeometricDistribution where
+  maybeEntropy = Just . D.entropy
 
 variance :: HypergeometricDistribution -> Double
 variance (HD m l k) = (k' * ml) * (1 - ml) * (l' - k') / (l' - 1)
@@ -73,6 +78,14 @@
 mean :: HypergeometricDistribution -> Double
 mean (HD m l k) = fromIntegral k * fromIntegral m / fromIntegral l
 {-# INLINE mean #-}
+
+directEntropy :: HypergeometricDistribution -> Double
+directEntropy d@(HD m _ _) =
+    negate . sum $
+  takeWhile (< negate m_epsilon) $
+  dropWhile (not . (< negate m_epsilon)) $
+  [ let x = probability d n in x * log x | n <- [0..m]]
+
 
 hypergeometric :: Int               -- ^ /m/
                -> Int               -- ^ /l/
diff --git a/Statistics/Distribution/Normal.hs b/Statistics/Distribution/Normal.hs
--- a/Statistics/Distribution/Normal.hs
+++ b/Statistics/Distribution/Normal.hs
@@ -46,7 +46,7 @@
     complCumulative = complCumulative
 
 instance D.ContDistr NormalDistribution where
-    density    = density
+    logDensity = logDensity
     quantile   = quantile
 
 instance D.MaybeMean NormalDistribution where
@@ -62,6 +62,12 @@
 instance D.Variance NormalDistribution where
     stdDev = stdDev
 
+instance D.Entropy NormalDistribution where
+  entropy d = 0.5 * log (2 * pi * exp 1 * D.variance d)
+
+instance D.MaybeEntropy NormalDistribution where
+  maybeEntropy = Just . D.entropy
+
 instance D.ContGen NormalDistribution where
     genContVar d = MWC.normal (mean d) (stdDev d)
     {-# INLINE genContVar #-}
@@ -84,7 +90,7 @@
 normalDistr m sd
   | sd > 0    = ND { mean       = m
                    , stdDev     = sd
-                   , ndPdfDenom = m_sqrt_2_pi * sd
+                   , ndPdfDenom = log $ m_sqrt_2_pi * sd
                    , ndCdfDenom = m_sqrt_2 * sd
                    }
   | otherwise =
@@ -99,8 +105,8 @@
   where
     (m,v) = S.meanVariance xs
 
-density :: NormalDistribution -> Double -> Double
-density d x = exp (-xm * xm / (2 * sd * sd)) / ndPdfDenom d
+logDensity :: NormalDistribution -> Double -> Double
+logDensity d x = (-xm * xm / (2 * sd * sd)) - ndPdfDenom d
     where xm = x - mean d
           sd = stdDev d
 
diff --git a/Statistics/Distribution/Poisson.hs b/Statistics/Distribution/Poisson.hs
--- a/Statistics/Distribution/Poisson.hs
+++ b/Statistics/Distribution/Poisson.hs
@@ -29,8 +29,8 @@
 import GHC.Generics (Generic)
 import qualified Statistics.Distribution as D
 import qualified Statistics.Distribution.Poisson.Internal as I
-import Numeric.SpecFunctions (incompleteGamma)
-
+import Numeric.SpecFunctions (incompleteGamma,logFactorial)
+import Numeric.MathFunctions.Constants (m_neg_inf)
 
 
 newtype PoissonDistribution = PD {
@@ -49,6 +49,9 @@
 
 instance D.DiscreteDistr PoissonDistribution where
     probability (PD lambda) x = I.probability lambda (fromIntegral x)
+    logProbability (PD lambda) i
+      | i < 0     = m_neg_inf
+      | otherwise = log lambda * fromIntegral i - logFactorial i - lambda
     {-# INLINE probability #-}
 
 instance D.Variance PoissonDistribution where
@@ -65,7 +68,12 @@
 instance D.MaybeVariance PoissonDistribution where
     maybeStdDev   = Just . D.stdDev
 
+instance D.Entropy PoissonDistribution where
+  entropy (PD lambda) = I.poissonEntropy lambda
 
+instance D.MaybeEntropy PoissonDistribution where
+  maybeEntropy = Just . D.entropy
+
 -- | Create Poisson distribution.
 poisson :: Double -> PoissonDistribution
 poisson l
@@ -78,3 +86,6 @@
 --
 -- * Loader, C. (2000) Fast and Accurate Computation of Binomial
 --   Probabilities. <http://projects.scipy.org/scipy/raw-attachment/ticket/620/loader2000Fast.pdf>
+-- * Adell, J., Lekuona, A., and Yu, Y. (2010) Sharp Bounds on the
+--   Entropy of the Poisson Law and Related Quantities
+--   <http://arxiv.org/pdf/1001.2897.pdf>
diff --git a/Statistics/Distribution/Poisson/Internal.hs b/Statistics/Distribution/Poisson/Internal.hs
--- a/Statistics/Distribution/Poisson/Internal.hs
+++ b/Statistics/Distribution/Poisson/Internal.hs
@@ -11,11 +11,12 @@
 
 module Statistics.Distribution.Poisson.Internal
     (
-      probability
+      probability, poissonEntropy
     ) where
 
-import Numeric.MathFunctions.Constants (m_sqrt_2_pi, m_tiny)
-import Numeric.SpecFunctions           (logGamma, stirlingError)
+import Data.List(unfoldr)
+import Numeric.MathFunctions.Constants (m_sqrt_2_pi, m_tiny, m_epsilon)
+import Numeric.SpecFunctions (logGamma, stirlingError, choose, logFactorial)
 import Numeric.SpecFunctions.Extra     (bd0)
 
 -- | An unchecked, non-integer-valued version of Loader's saddle point
@@ -31,3 +32,147 @@
   | otherwise            = exp (-(stirlingError x) - bd0 x lambda) /
                            (m_sqrt_2_pi * sqrt x)
 {-# INLINE probability #-}
+
+-- | Compute entropy using Theorem 1 from "Sharp Bounds on the Entropy
+-- of the Poisson Law".  This function is unused because 'directEntorpy'
+-- is just as accurate and is faster by about a factor of 4.
+alyThm1 :: Double -> Double
+alyThm1 lambda =
+  sum (takeWhile (\x -> abs x >= m_epsilon * lll) alySeries) + lll
+  where lll = lambda * (1 - log lambda)
+        alySeries =
+          [ alyc k * exp (fromIntegral k * log lambda - logFactorial k)
+          | k <- [2..] ]
+
+alyc :: Int -> Double
+alyc k =
+  sum [ parity j * choose (k-1) j * log (fromIntegral j+1) | j <- [0..k-1] ]
+  where parity j
+          | even (k-j) = -1
+          | otherwise  = 1
+                         
+-- | Returns [x, x^2, x^3, x^4, ...]
+powers :: Double -> [Double]
+powers x = unfoldr (\y -> Just (y*x,y*x)) 1
+
+-- | Returns an upper bound according to theorem 2 of "Sharp Bounds on
+-- the Entropy of the Poisson Law"
+alyThm2Upper :: Double -> [Double] -> Double
+alyThm2Upper lambda coefficients =
+  1.4189385332046727 + 0.5 * log lambda +
+  zipCoefficients lambda coefficients
+
+-- | Returns the average of the upper and lower bounds accounding to
+-- theorem 2.
+alyThm2 :: Double -> [Double] -> [Double] -> Double
+alyThm2 lambda upper lower =
+  alyThm2Upper lambda upper + 0.5 * (zipCoefficients lambda lower)
+
+zipCoefficients :: Double -> [Double] -> Double 
+zipCoefficients lambda coefficients =
+  (sum $ map (uncurry (*)) (zip (powers $ recip lambda) coefficients))
+
+-- Mathematica code deriving the coefficients below:
+--
+-- poissonMoment[0, s_] := 1
+-- poissonMoment[1, s_] := 0
+-- poissonMoment[k_, s_] := 
+--   Sum[s * Binomial[k - 1, j] * poissonMoment[j, s], {j, 0, k - 2}]
+--
+-- upperSeries[m_]  :=
+--  Distribute[Integrate[
+--    Sum[(-1)^(j - 1) * 
+--      poissonMoment[j, \[Lambda]] / (j * (j - 1)* \[Lambda]^j),
+--     {j, 3, 2 m - 1}],
+--    \[Lambda]]]
+--
+-- lowerSeries[m_] :=
+--  Distribute[Integrate[
+--    poissonMoment[
+--      2 m + 2, \[Lambda]] / ((2 m + 
+--         1)*\[Lambda]^(2 m + 2)), \[Lambda]]]
+--
+-- upperBound[m_] := upperSeries[m] + (Log[2*Pi*\[Lambda]] + 1)/2 
+--
+-- lowerBound[m_] := upperBound[m] + lowerSeries[m]
+
+upperCoefficients4 :: [Double]
+upperCoefficients4 = [1/12, 1/24, -103/180, -13/40, -1/210]
+
+lowerCoefficients4 :: [Double]
+lowerCoefficients4 = [0,0,0, -105/4, -210, -2275/18, -167/21, -1/72]
+
+upperCoefficients6 :: [Double]
+upperCoefficients6 = [1/12, 1/24, 19/360, 9/80, -38827/2520,
+                      -74855/1008, -73061/2520, -827/720, -1/990]
+
+lowerCoefficients6 :: [Double]
+lowerCoefficients6 = [0,0,0,0,0, -3465/2, -45045, -466235/4, -531916/9,
+                      -56287/10, -629/11, -1/156]
+
+upperCoefficients8 :: [Double]
+upperCoefficients8 = [1/12, 1/24, 19/360, 9/80, 863/2520, 1375/1008,
+                      -3023561/2520, -15174047/720, -231835511/5940,
+                      -18927611/1320, -58315591/60060, -23641/3640,
+                      -1/2730]
+
+lowerCoefficients8 :: [Double]
+lowerCoefficients8 = [0,0,0,0,0,0,0, -2027025/8, -15315300, -105252147,
+                      -178127950, -343908565/4, -10929270, -3721149/14,
+                      -7709/15, -1/272]
+  
+upperCoefficients10 :: [Double]
+upperCoefficients10 = [1/12, 1/24, 19/360, 9,80, 863/2520, 1375/1008,
+                       33953/5040, 57281/1440, -2271071617/11880,
+                       -1483674219/176, -31714406276557/720720,
+                       -7531072742237/131040, -1405507544003/65520,
+                       -21001919627/10080, -1365808297/36720,
+                       -26059/544, -1/5814]
+                      
+lowerCoefficients10 :: [Double]
+lowerCoefficients10 = [0,0,0,0,0,0,0,0,0,-130945815/2, -7638505875,
+                       -438256243425/4, -435477637540, -3552526473925/6,
+                       -857611717105/3, -545654955967/12, -5794690528/3,
+                       -578334559/42, -699043/133, -1/420]
+                 
+upperCoefficients12 :: [Double]
+upperCoefficients12 = [1/12, 1/24, 19/360, 863/2520, 1375/1008,
+                       33953/5040, 57281/1440, 3250433/11880,
+                       378351/176, -37521922090657/720720,
+                       -612415657466657/131040, -3476857538815223/65520,
+                       -243882174660761/1440, -34160796727900637/183600,
+                       -39453820646687/544, -750984629069237/81396,
+                       -2934056300989/9576, -20394527513/12540,
+                       -3829559/9240, -1/10626]
+
+lowerCoefficients12 :: [Double]
+lowerCoefficients12 = [0,0,0,0,0,0,0,0,0,0,0,
+                       -105411381075/4, -5270569053750, -272908057767345/2,
+                       -1051953238104769, -24557168490009155/8,
+                       -3683261873403112, -5461918738302026/3,
+                       -347362037754732, -2205885452434521/100,
+                       -12237195698286/35, -16926981721/22,
+                       -6710881/155, -1/600]
+                      
+-- | Compute entropy directly from its definition. This is just as accurate
+-- as 'alyThm1' for lambda <= 1 and is faster, but is slow for large lambda,
+-- and produces some underestimation due to accumulation of floating point
+-- error.
+directEntropy :: Double -> Double
+directEntropy lambda =   
+  negate . sum $
+  takeWhile (< negate m_epsilon * lambda) $
+  dropWhile (not . (< negate m_epsilon * lambda)) $
+  [ let x = probability lambda k in x * log x | k <- [0..]]
+
+-- | Compute the entropy of a poisson distribution using the best available
+-- method.
+poissonEntropy :: Double -> Double
+poissonEntropy lambda
+  | lambda == 0 = 0
+  | lambda <= 10 = directEntropy lambda
+  | lambda <= 12 = alyThm2 lambda upperCoefficients4 lowerCoefficients4
+  | lambda <= 18 = alyThm2 lambda upperCoefficients6 lowerCoefficients6
+  | lambda <= 24 = alyThm2 lambda upperCoefficients8 lowerCoefficients8
+  | lambda <= 30 = alyThm2 lambda upperCoefficients10 lowerCoefficients10
+  | otherwise = alyThm2 lambda upperCoefficients12 lowerCoefficients12
diff --git a/Statistics/Distribution/StudentT.hs b/Statistics/Distribution/StudentT.hs
--- a/Statistics/Distribution/StudentT.hs
+++ b/Statistics/Distribution/StudentT.hs
@@ -21,7 +21,8 @@
 import GHC.Generics (Generic)
 import qualified Statistics.Distribution as D
 import Statistics.Distribution.Transform (LinearTransform (..))
-import Numeric.SpecFunctions (logBeta, incompleteBeta, invIncompleteBeta)
+import Numeric.SpecFunctions (
+  logBeta, incompleteBeta, invIncompleteBeta, digamma)
 
 -- | Student-T distribution
 newtype StudentT = StudentT { studentTndf :: Double }
@@ -39,7 +40,8 @@
   cumulative = cumulative
 
 instance D.ContDistr StudentT where
-  density  = density
+  density    d@(StudentT ndf) x = exp (logDensityUnscaled d x) / sqrt ndf
+  logDensity d@(StudentT ndf) x = logDensityUnscaled d x - log (sqrt ndf)
   quantile = quantile
 
 cumulative :: StudentT -> Double -> Double
@@ -49,9 +51,9 @@
   where
     ibeta = incompleteBeta (0.5 * ndf) 0.5 (ndf / (ndf + x*x))
 
-density :: StudentT -> Double -> Double
-density (StudentT ndf) x =
-    exp( log (ndf / (ndf + x*x)) * (0.5 * (1 + ndf)) - logBeta 0.5 (0.5 * ndf) ) / sqrt ndf
+logDensityUnscaled :: StudentT -> Double -> Double
+logDensityUnscaled (StudentT ndf) x =
+    log (ndf / (ndf + x*x)) * (0.5 * (1 + ndf)) - logBeta 0.5 (0.5 * ndf)
 
 quantile :: StudentT -> Double -> Double
 quantile (StudentT ndf) p
@@ -70,6 +72,15 @@
 instance D.MaybeVariance StudentT where
   maybeVariance (StudentT ndf) | ndf > 2   = Just $! ndf / (ndf - 2)
                                | otherwise = Nothing
+
+instance D.Entropy StudentT where
+  entropy (StudentT ndf) =
+    0.5 * (ndf+1) * (digamma ((1+ndf)/2) - digamma(ndf/2))
+    + log (sqrt ndf) 
+    + logBeta (ndf/2) 0.5
+
+instance D.MaybeEntropy StudentT where
+  maybeEntropy = Just . D.entropy
 
 instance D.ContGen StudentT where
   genContVar = D.genContinous
diff --git a/Statistics/Distribution/Transform.hs b/Statistics/Distribution/Transform.hs
--- a/Statistics/Distribution/Transform.hs
+++ b/Statistics/Distribution/Transform.hs
@@ -55,7 +55,8 @@
   cumulative (LinearTransform loc sc dist) x = D.cumulative dist $ (x-loc) / sc
 
 instance D.ContDistr d => D.ContDistr (LinearTransform d) where
-  density  (LinearTransform loc sc dist) x = D.density dist ((x-loc) / sc) / sc
+  density    (LinearTransform loc sc dist) x = D.density    dist ((x-loc) / sc) / sc
+  logDensity (LinearTransform loc sc dist) x = D.logDensity dist ((x-loc) / sc) - log sc
   quantile (LinearTransform loc sc dist) p = loc + sc * D.quantile dist p
 
 instance D.MaybeMean d => D.MaybeMean (LinearTransform d) where
@@ -71,6 +72,14 @@
 instance (D.Variance d) => D.Variance (LinearTransform d) where
   variance (LinearTransform _ sc dist) = sc * sc * D.variance dist
   stdDev   (LinearTransform _ sc dist) = sc * D.stdDev dist
+
+instance (D.MaybeEntropy d, D.DiscreteDistr d) 
+         => D.MaybeEntropy (LinearTransform d) where
+  maybeEntropy (LinearTransform _ _ dist) = D.maybeEntropy dist
+
+instance (D.Entropy d, D.DiscreteDistr d) 
+         => D.Entropy (LinearTransform d) where
+  entropy (LinearTransform _ _ dist) = D.entropy dist
 
 instance D.ContGen d => D.ContGen (LinearTransform d) where
   genContVar (LinearTransform loc sc d) g = do
diff --git a/Statistics/Distribution/Uniform.hs b/Statistics/Distribution/Uniform.hs
--- a/Statistics/Distribution/Uniform.hs
+++ b/Statistics/Distribution/Uniform.hs
@@ -73,5 +73,11 @@
 instance D.MaybeVariance UniformDistribution where
     maybeStdDev   = Just . D.stdDev
 
+instance D.Entropy UniformDistribution where
+  entropy (UniformDistribution a b) = log (b - a)
+
+instance D.MaybeEntropy UniformDistribution where
+  maybeEntropy = Just . D.entropy
+
 instance D.ContGen UniformDistribution where
     genContVar (UniformDistribution a b) gen = MWC.uniformR (a,b) gen
diff --git a/Statistics/Sample.hs b/Statistics/Sample.hs
--- a/Statistics/Sample.hs
+++ b/Statistics/Sample.hs
@@ -203,7 +203,7 @@
 data V = V {-# UNPACK #-} !Double {-# UNPACK #-} !Double
 
 robustSumVar :: (G.Vector v Double) => Double -> v Double -> Double
-robustSumVar m samp = G.sum . G.map (square . subtract m) $ samp
+robustSumVar m = G.sum . G.map (square . subtract m)
   where square x = x * x
 {-# INLINE robustSumVar #-}
 
diff --git a/statistics.cabal b/statistics.cabal
--- a/statistics.cabal
+++ b/statistics.cabal
@@ -1,5 +1,5 @@
 name:           statistics
-version:        0.10.4.1
+version:        0.10.5.0
 synopsis:       A library of statistical types, data, and functions
 description:
   This library provides a number of common functions and types useful
@@ -21,49 +21,6 @@
   .
   * Common statistical tests for significant differences between
     samples.
-  .
-  Changes in 0.10.4.0
-  .
-  * Support for versions of GHC older than 7.2 is discontinued.
-  .
-  * All datatypes now support 'Data.Binary' and 'GHC.Generics'.
-  .
-  Changes in 0.10.3.0
-  .
-  * Bug fixes
-  .
-  Changes in 0.10.2.0
-  .
-  * Bugs in DCT and IDCT are fixed.
-  .
-  * Accesors for uniform distribution are added.
-  .
-  * 'ContGen' instances for all continous distribtuions are added.
-  .
-  * Beta distribution is added.
-  .
-  * Constructor for improper gamma distribtuion is added.
-  .
-  * Binomial distribution allows zero trials.
-  .
-  * Poisson distribution now accept zero parameter.
-  .
-  * Integer overflow in caculation of Wilcoxon-T test is fixed.
-  .
-  * Bug in 'ContGen' instance for normal distribution is fixed.
-  .
-  Changes in 0.10.1.0
-  .
-  * Kolmogorov-Smirnov nonparametric test added.
-  .
-  * Pearson's chi squared test added.
-  .
-  * Type class for generating random variates for given distribution
-    is added.
-  .
-  * Modules 'Statistics.Math' and 'Statistics.Constants' are moved to
-    the @math-functions@ package. They are still available but marked
-    as deprecated.
 
 license:        BSD3
 license-file:   LICENSE
@@ -137,7 +94,7 @@
     deepseq >= 1.1.0.2,
     erf,
     monad-par         >= 0.3.4,
-    mwc-random        >= 0.11.0.0,
+    mwc-random        >= 0.13.0.0,
     math-functions    >= 0.1.2,
     primitive         >= 0.3,
     vector            >= 0.7.1,
diff --git a/tests/Tests/Distribution.hs b/tests/Tests/Distribution.hs
--- a/tests/Tests/Distribution.hs
+++ b/tests/Tests/Distribution.hs
@@ -59,6 +59,7 @@
 
   , discreteDistrTests (T :: T BinomialDistribution       )
   , discreteDistrTests (T :: T GeometricDistribution      )
+  , discreteDistrTests (T :: T GeometricDistribution0     )
   , discreteDistrTests (T :: T HypergeometricDistribution )
   , discreteDistrTests (T :: T PoissonDistribution        )
 
@@ -73,9 +74,10 @@
 contDistrTests :: (Param d, ContDistr d, QC.Arbitrary d, Typeable d, Show d) => T d -> Test
 contDistrTests t = testGroup ("Tests for: " ++ typeName t) $
   cdfTests t ++
-  [ testProperty "PDF sanity"              $ pdfSanityCheck   t
-  , testProperty "Quantile is CDF inverse" $ quantileIsInvCDF t
+  [ testProperty "PDF sanity"              $ pdfSanityCheck     t
+  , testProperty "Quantile is CDF inverse" $ quantileIsInvCDF   t
   , testProperty "quantile fails p<0||p>1" $ quantileShouldFail t
+  , testProperty "log density check"       $ logDensityCheck    t
   ]
 
 -- Tests for discrete distribution
@@ -85,6 +87,7 @@
   [ testProperty "Prob. sanity"         $ probSanityCheck       t
   , testProperty "CDF is sum of prob."  $ discreteCDFcorrect    t
   , testProperty "Discrete CDF is OK"   $ cdfDiscreteIsCorrect  t
+  , testProperty "log probabilty check" $ logProbabilityCheck   t
   ]
 
 -- Tests for distributions which have CDF
@@ -144,7 +147,7 @@
 -- CDF for discrete distribution uses <= for comparison
 cdfDiscreteIsCorrect :: (DiscreteDistr d) => T d -> d -> Property
 cdfDiscreteIsCorrect _ d
-  = printTestCase (unlines $ map show badN)
+  = printTestCase (unlines badN)
   $ null badN  
   where
     -- We are checking that:
@@ -154,7 +157,7 @@
     -- Apporixmate equality is tricky here. Scale is set by maximum
     -- value of CDF and probability. Case when all proabilities are
     -- zero should be trated specially.
-    badN = [ (i,p,p1,dp, (p1-p-dp) / max p1 dp)
+    badN = [ printf "N=%3i    p[i]=%g\tp[i+1]=%g\tdP=%g\trelerr=%g" i p p1 dp ((p1-p-dp) / max p1 dp)
            | i <- [0 .. 100]
            , let p      = cumulative d $ fromIntegral i - 1e-6
                  p1     = cumulative d $ fromIntegral i
@@ -164,6 +167,19 @@
            && relerr > 1e-14
            ]
 
+logDensityCheck :: (ContDistr d) => T d -> d -> Double -> Property
+logDensityCheck _ d x
+  = printTestCase (printf "density    = %g" p)
+  $ printTestCase (printf "logDensity = %g" logP)
+  $ printTestCase (printf "log p      = %g" (log p))
+  $ printTestCase (printf "eps        = %g" (abs (logP - log p) / max (abs (log p)) (abs logP)))
+  $ or [ p == 0     && logP == (-1/0)
+       , p < 1e-308 && logP < 609
+       , eq 1e-14 (log p) logP
+       ]
+  where
+    p    = density d x
+    logP = logDensity d x
 
 -- PDF is positive
 pdfSanityCheck :: (ContDistr d) => T d -> d -> Double -> Bool
@@ -214,7 +230,21 @@
     p1 = cumulative d (fromIntegral m + 0.5) - cumulative d (fromIntegral n - 0.5)
     p2 = sum $ map (probability d) [n .. m]
 
+logProbabilityCheck :: (DiscreteDistr d) => T d -> d -> Int -> Property
+logProbabilityCheck _ d x
+  = printTestCase (printf "probability    = %g" p)
+  $ printTestCase (printf "logProbability = %g" logP)
+  $ printTestCase (printf "log p          = %g" (log p))
+  $ printTestCase (printf "eps            = %g" (abs (logP - log p) / max (abs (log p)) (abs logP)))
+  $ or [ p == 0     && logP == (-1/0)
+       , p < 1e-308 && logP < 609
+       , eq 1e-14 (log p) logP
+       ]
+  where
+    p    = probability d x
+    logP = logProbability d x
 
+
     
 ----------------------------------------------------------------
 -- Arbitrary instances for ditributions
@@ -230,6 +260,8 @@
   arbitrary = betaDistr <$> QC.choose (1e-3,10) <*> QC.choose (1e-3,10)
 instance QC.Arbitrary GeometricDistribution where
   arbitrary = geometric <$> QC.choose (0,1)
+instance QC.Arbitrary GeometricDistribution0 where
+  arbitrary = geometric0 <$> QC.choose (0,1)
 instance QC.Arbitrary HypergeometricDistribution where
   arbitrary = do l <- QC.choose (1,20)
                  m <- QC.choose (0,l)