diff --git a/Numeric/SpecFunctions/Internal.hs b/Numeric/SpecFunctions/Internal.hs
--- a/Numeric/SpecFunctions/Internal.hs
+++ b/Numeric/SpecFunctions/Internal.hs
@@ -777,7 +777,7 @@
     in x
   -- For small a and not too big b we use approximation from boost.
   | b < 5 && a <= 1 =
-    let x | p**(1/a) < 0.5 = (p * a * exp beta) ** (1/a)
+    let x | p**(1/a) < 0.5 = (p ** (a * exp beta)) ** (1/a)
           | otherwise      = 1 - (1 - p ** (b * exp beta))**(1/b)
     in x
   -- When a>>b and both are large approximation from [Temme1992],
diff --git a/changelog.md b/changelog.md
--- a/changelog.md
+++ b/changelog.md
@@ -1,3 +1,9 @@
+## Changes in 0.3.4.2
+
+  * Fixed crash in `invIncompleteBeta` (#68) for some inputs initial approximation
+    was computed incorrectly.
+
+
 ## Changes in 0.3.4.1
 
   * Precision of `incompleteGamma` improved.
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.3.4.1
+version:        0.3.4.2
 cabal-version:  >= 1.10
 license:        BSD2
 license-file:   LICENSE
@@ -28,8 +28,9 @@
      || ==8.2.2
      || ==8.4.4
      || ==8.6.5
-     || ==8.8.3
-     || ==8.10.1
+     || ==8.8.4
+     || ==8.10.2
+     || ==9.0.1
   , GHCJS ==8.4
 
 extra-source-files:
diff --git a/tests/Tests/SpecFunctions.hs b/tests/Tests/SpecFunctions.hs
--- a/tests/Tests/SpecFunctions.hs
+++ b/tests/Tests/SpecFunctions.hs
@@ -26,14 +26,20 @@
 import Numeric.MathFunctions.Constants  (m_epsilon,m_tiny)
 
 erfTol,erfcTol,erfcLargeTol :: Int
-#if USE_SYSTEM_ERF && !defined(__GHCJS__)
-erfTol       = 1
-erfcTol      = 2
-erfcLargeTol = 2
-#else
+-- Pure haskell implementation is not very good
+#if !defined(USE_SYSTEM_ERF) || defined(__GHCJS__)
 erfTol       = 4
 erfcTol      = 4
 erfcLargeTol = 64
+-- Macos's erf is slightly worse that GNU one
+#elif defined(darwin_HOST_OS)
+erfTol       = 2
+erfcTol      = 2
+erfcLargeTol = 2
+#else
+erfTol       = 1
+erfcTol      = 2
+erfcLargeTol = 2
 #endif
 
 isGHCJS :: Bool
@@ -43,6 +49,14 @@
 isGHCJS = False
 #endif
 
+isWindows :: Bool
+#if defined(mingw32_HOST_OS)
+isWindows = True
+#else
+isWindows = False
+#endif
+
+
 tests :: TestTree
 tests = testGroup "Special functions"
   [ testGroup "erf"
@@ -88,7 +102,10 @@
     [ testCase "incompleteGamma table" $
         forTable "tests/tables/igamma.dat" $ \[a,x,exact] -> do
           let err | a < 10    = 16
-                  | a <= 101  = if isGHCJS then 64 else 32
+                  | a <= 101  = case () of
+                      _| isGHCJS   -> 64
+                       | isWindows -> 64
+                       | otherwise -> 32
                   | a == 201  = 200
                   | otherwise = 32
           checkTabularPure err (show (a,x)) exact (incompleteGamma a x)
@@ -137,6 +154,12 @@
     --
     , testProperty "0 <= I[B] <= 1" incompleteBetaInRange
     , testProperty "ibeta symmetry" incompleteBetaSymmetry
+    , testCase "Regression #68" $ do
+        let a = 1
+            b = 0.3
+            p = 0.3
+            x = invIncompleteBeta a b p
+        assertBool "Inversion OK" $ incompleteBeta a b x `ulpDistance` p < 4
     -- XXX FIXME DISABLED due to failures
     -- , testProperty "invIncompleteBeta  = B^-1" $ invIBetaIsInverse
     ]
@@ -177,7 +200,7 @@
     [ testCase "choose table" $
       forM_ [0 .. 1000] $ \n ->
         forM_ [0 .. n]  $ \k -> do
-          checkTabular 2048
+          checkTabular (if isWindows then 3072 else 2048)
             (show (n,k))
             (fromIntegral $ choose' n k)
             (choose (fromInteger n) (fromInteger k))
@@ -205,14 +228,18 @@
   roundtrip_invErfc_erfc,
   roundtrip_erf_invErf
   :: (Double,Double)
-#if USE_SYSTEM_ERF && !defined(__GHCJS__)
-roundtrip_erfc_invErfc = (2,2)
-roundtrip_invErfc_erfc = (2,2)
-roundtrip_erf_invErf   = (1,1)
-#else
+#if !defined(USE_SYSTEM_ERF) || defined(__GHCJS__)
 roundtrip_erfc_invErfc = (2,8)
 roundtrip_invErfc_erfc = (8,4)
 roundtrip_erf_invErf   = (128,128)
+#elif defined(darwin_HOST_OS)
+roundtrip_erfc_invErfc = (4,4)
+roundtrip_invErfc_erfc = (4,4)
+roundtrip_erf_invErf   = (2,2)
+#else
+roundtrip_erfc_invErfc = (2,2)
+roundtrip_invErfc_erfc = (2,2)
+roundtrip_erf_invErf   = (1,1)
 #endif
 
 -- id ≈ erfc . invErfc