diff --git a/Numeric/SpecFunctions/Internal.hs b/Numeric/SpecFunctions/Internal.hs
--- a/Numeric/SpecFunctions/Internal.hs
+++ b/Numeric/SpecFunctions/Internal.hs
@@ -352,7 +352,7 @@
 -- | Compute the natural logarithm of the beta function.
 --
 -- \[
--- B(a,b) = \int_0^1 t^{a-1}(1-t)^{1-b}\,dt = \frac{\Gamma{a}\Gamma{b}}{\Gamma{a+b}}
+-- B(a,b) = \int_0^1 t^{a-1}(1-t)^{1-b}\,dt = \frac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)}
 -- \]
 logBeta
   :: Double                     -- ^ /a/ > 0
@@ -490,12 +490,10 @@
   | a <  0 || a >  1 =
       modErr $ printf "invIncompleteBeta x must be in [0,1].  p=%g q=%g a=%g" p q a
   | a == 0 || a == 1 = a
-  | a > 0.5          = 1 - invIncompleteBetaWorker (logBeta p q) q p (1 - a)
-  | otherwise        =     invIncompleteBetaWorker (logBeta p q) p q  a
+  | otherwise        = invIncompleteBetaWorker (logBeta p q) p q  a
 
 
 invIncompleteBetaWorker :: Double -> Double -> Double -> Double -> Double
--- NOTE: p <= 0.5.
 invIncompleteBetaWorker beta a b p = loop (0::Int) (invIncBetaGuess beta a b p)
   where
     a1 = a - 1
@@ -562,14 +560,14 @@
     in x
   -- If both a and b larger or equal that 1 but not too big we use
   -- same approximation as above but calculate it a bit differently
-  | a+b <= 6 && a>=1 && b>=1 =
+  | a+b <= 6 && a>1 && b>1 =
     let x_infl = (a - 1) / (a + b - 2)
         p_infl = incompleteBeta a b x_infl
         x | p < p_infl = exp ((log(p * a) + beta) / a)
           | otherwise  = 1 - exp((log((1-p) * b) + beta) / b)
     in x
   -- For small a and not too big b we use approximation from boost.
-  | b < 5 && a < 1 =
+  | b < 5 && a <= 1 =
     let x | p**(1/a) < 0.5 = (p * a * exp beta) ** (1/a)
           | otherwise      = 1 - (1 - p ** (b * exp beta))**(1/b)
     in x
@@ -829,7 +827,9 @@
 -- | Calculate the error term of the Stirling approximation.  This is
 -- only defined for non-negative values.
 --
--- > stirlingError @n@ = @log(n!) - log(sqrt(2*pi*n)*(n/e)^n)
+-- \[
+-- \operatorname{stirlingError}(n) = \log(n!) - \log(\sqrt{2\pi n}\frac{n}{e}^n)
+-- \]
 stirlingError :: Double -> Double
 stirlingError n
   | n <= 15.0   = case properFraction (n+n) of
diff --git a/changelog.md b/changelog.md
--- a/changelog.md
+++ b/changelog.md
@@ -1,5 +1,10 @@
 Changes in 0.2.0.0
 
+  * Bug fixes and documentation tweaks
+
+
+Changes in 0.2.0.0
+
   * `logGamma` now uses Lancsoz approximation and same as `logGammaL`.  Old
      implementation of `logGamma` moved to `Numeric.SpecFunctions.Extra.logGammaAS245`.
 
diff --git a/math-functions.cabal b/math-functions.cabal
--- a/math-functions.cabal
+++ b/math-functions.cabal
@@ -1,5 +1,5 @@
 name:           math-functions
-version:        0.2.0.0
+version:        0.2.0.1
 cabal-version:  >= 1.10
 license:        BSD3
 license-file:   LICENSE
diff --git a/tests/Tests/SpecFunctions.hs b/tests/Tests/SpecFunctions.hs
--- a/tests/Tests/SpecFunctions.hs
+++ b/tests/Tests/SpecFunctions.hs
@@ -50,9 +50,8 @@
       $ and [ eq 1e-15 (logGammaL (fromIntegral n))
                        (logFactorial (n-1))
             | n <- [3..10000::Int]]
-    -- FIXME: Too low!
   , testAssertion "logGammaL is expected to be precise at 1e-10 level [fractional points]"
-      $ and [ eq 1e-10 (logGammaL x) lg | (x,lg) <- tableLogGamma ]
+      $ and [ eq (64*m_epsilon) (logGammaL x) lg | (x,lg) <- tableLogGamma ]
     -- FIXME: loss of precision when logBeta p q ≈ 0.
     --        Relative error doesn't work properly in this case.
   , testAssertion "logBeta is expected to be precise at 1e-6 level"
@@ -66,9 +65,8 @@
     -- Relative precision is lost when digamma(x) ≈ 0
   , testAssertion "digamma is expected to be precise at 1e-12"
       $ and [ eq 1e-12 r (digamma x) | (x,r) <- tableDigamma ]
-    -- FIXME: Why 1e-8? Is it due to poor precision of logBeta?
-  , testAssertion "incompleteBeta is expected to be precise at 1e-8 level"
-      $ and [ eq 1e-8 (incompleteBeta p q x) ib | (p,q,x,ib) <- tableIncompleteBeta ]
+  , testAssertion "incompleteBeta is expected to be precise at 32*m_epsilon level"
+      $ and [ eq (32 * m_epsilon) (incompleteBeta p q x) ib | (p,q,x,ib) <- tableIncompleteBeta ]
   , testAssertion "incompleteBeta with p > 3000 and q > 3000"
       $ and [ eq 1e-11 (incompleteBeta p q x) ib | (x,p,q,ib) <-
                  [ (0.495,  3001,  3001, 0.2192546757957825068677527085659175689142653854877723)