diff --git a/Statistics/Distribution/Binomial.hs b/Statistics/Distribution/Binomial.hs
--- a/Statistics/Distribution/Binomial.hs
+++ b/Statistics/Distribution/Binomial.hs
@@ -27,8 +27,10 @@
 import Data.Array.Vector
 import Data.Int (Int64)
 import Data.Typeable (Typeable)
+import Statistics.Constants (m_epsilon)
 import qualified Statistics.Distribution as D
-import Statistics.Math (choose)
+import Statistics.Distribution.Normal (standard)
+import Statistics.Math (choose, logFactorial)
 
 -- | The binomial distribution.
 data BinomialDistribution = BD {
@@ -50,17 +52,69 @@
     mean = mean
 
 density :: BinomialDistribution -> Double -> Double
-density (BD n p) x =
-    fromIntegral (n `choose` floor x) * p ** x * (1-p) ** (fromIntegral n-x)
+density (BD n p) x
+    | not (isIntegral x) = integralError "density"
+    | n == 0             = 1
+    | x < 0 || x > n'    = 0
+    | n <= 50 || x < 2   = sign * p'' ** x' * (n `choose` fx) * q'' ** nx'
+    | otherwise          = sign * exp (x' * log p'' + nx' * log q'' + lf)
+  where sign = oddX * oddNX
+        (x',p',q') | x > n' / 2 = (n'-x, q, p)
+                   | otherwise  = (x,    p, q)
+        oddX | p' < 0 && odd fx     = -1
+             | otherwise            = 1
+        oddNX | q' < 0 && odd nx    = -1
+              | otherwise           = 1
+        p'' = abs p'
+        q'' = abs q'
+        q   = 1 - p
+        nx  = n - fx
+        nx' = fromIntegral nx
+        fx  = floor x'
+        n'  = fromIntegral n
+        lf  = logFactorial n - logFactorial nx - logFactorial fx
 
 cumulative :: BinomialDistribution -> Double -> Double
-cumulative d =
-    sumU . mapU (density d . fromIntegral) . enumFromToU (0::Int) . floor
+cumulative d x
+  | isIntegral x = sumU . mapU (density d . fromIntegral) . enumFromToU (0::Int) . floor $ x
+  | otherwise    = integralError "cumulative"
 
+isIntegral :: Double -> Bool
+isIntegral x = x == floorf x
+
+floorf :: Double -> Double
+floorf = fromIntegral . (floor :: Double -> Int64)
+
 quantile :: BinomialDistribution -> Double -> Double
-quantile d@(BD n _p) p = fromIntegral . r64 $ D.findRoot d p (n'/2) 0 n'
-    where n'  = fromIntegral n
-          r64 = round :: Double -> Int64
+quantile dist@(BD n p) prob
+    | isNaN prob = prob
+    | p == 1     = n'
+    | n' < 1e5   = fst (search 1 y0 z0)
+    | otherwise  = let dy = floorf (n' / 1000)
+                   in  narrow dy (search dy y0 z0)
+  where q  = 1 - p
+        n' = fromIntegral n
+        y0 = n' `min` floorf (µ + σ * (d + γ * (d * d - 1) / 6) + 0.5)
+          where µ  = n' * p
+                σ  = sqrt (n' * p * q)
+                d = D.quantile standard prob
+                γ  = (q - p) / σ
+        z0 = cumulative dist y0
+        search dy y1 z1 | z0 >= prob' = left y1 z1
+                        | otherwise   = right y1
+          where
+            prob' = prob * (1 - 64 * m_epsilon)
+            left y oldZ | y == 0 || z < prob' = (y, oldZ)
+                        | otherwise           = left (max 0 y') z
+                where z  = cumulative dist y'
+                      y' = y - dy
+            right y | y' >= n' || z >= prob' = (y', z)
+                    | otherwise              = right y'
+                where z  = cumulative dist y'
+                      y' = y + dy
+        narrow dy (y,z) | dy <= 1 || dy' <= n'/1e15 = y
+                        | otherwise                 = narrow dy' (search dy y z)
+            where dy' = floorf (dy / 100)
 
 mean :: BinomialDistribution -> Double
 mean (BD n p) = fromIntegral n * p
@@ -78,3 +132,7 @@
     assert (p > 0 && p < 1) $
     BD n p
 {-# INLINE binomial #-}
+
+integralError :: String -> a
+integralError f = error ("Statistics.Distribution.Binomial." ++ f ++
+                         ": non-integer-valued input")
diff --git a/Statistics/Distribution/Hypergeometric.hs b/Statistics/Distribution/Hypergeometric.hs
--- a/Statistics/Distribution/Hypergeometric.hs
+++ b/Statistics/Distribution/Hypergeometric.hs
@@ -80,7 +80,7 @@
     | r > maxVal = 1/0
     | otherwise  = exp r
   where
-    a <> b = fromIntegral (a `choose` b)
+    a <> b = a `choose` b
     r = f m + f (l-m) - f l - f xi - f (k-xi) + f k -
         f (m-xi) - f (l-m-k+xi) + f (l-k)
     f = logFactorial
diff --git a/Statistics/Math.hs b/Statistics/Math.hs
--- a/Statistics/Math.hs
+++ b/Statistics/Math.hs
@@ -49,15 +49,17 @@
 -- | The binomial coefficient.
 --
 -- > 7 `choose` 3 == 35
-choose :: Int -> Int -> Int
+choose :: Int -> Int -> Double
 n `choose` k
-    | k > n = 0
-    | otherwise = ceiling . foldlU go 1 . enumFromToU 1 $ k'
+    | k > n     = 0
+    | k < 30    = foldlU go 1 . enumFromToU 1 $ k'
+    | otherwise = exp $ lg (n+1) - lg (k+1) - lg (n-k+1)
     where go a i = a * (nk + j) / j
               where j = fromIntegral i :: Double
           k' | k > n `div` 2 = n - k
              | otherwise     = k
           nk = fromIntegral (n - k')
+          lg = logGamma . fromIntegral
 {-# INLINE choose #-}
 
 data F = F {-# UNPACK #-} !Word64 {-# UNPACK #-} !Word64
diff --git a/Statistics/Sample/Powers.hs b/Statistics/Sample/Powers.hs
--- a/Statistics/Sample/Powers.hs
+++ b/Statistics/Sample/Powers.hs
@@ -107,7 +107,7 @@
     | k == 0    = 1
     | otherwise = (/n) . sumU . mapU go . indexedU . takeU (k+1) $ pa
   where
-    go (i :*: e) = fromIntegral (k `choose` i) * ((-m) ^ (k-i)) * e
+    go (i :*: e) = (k `choose` i) * ((-m) ^ (k-i)) * e
     n = indexU pa 0
     m = mean p
 {-# INLINE centralMoment #-}
diff --git a/statistics.cabal b/statistics.cabal
--- a/statistics.cabal
+++ b/statistics.cabal
@@ -1,5 +1,5 @@
 name:           statistics
-version:        0.3.6
+version:        0.4.0
 synopsis:       A library of statistical types, data, and functions
 description:
   This library provides a number of common functions and types useful
@@ -65,6 +65,6 @@
   -- gather extensive profiling data for now
   ghc-prof-options: -auto-all
 
-  ghc-options: -Wall -funbox-strict-fields -O2
+  ghc-options: -Wall -funbox-strict-fields
   if impl(ghc >= 6.8)
     ghc-options: -fwarn-tabs