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.1.1.0
+version:        0.1.1.1
 cabal-version:  >= 1.8
 license:        BSD3
 license-file:   LICENSE
@@ -16,6 +16,8 @@
   useful in statistical and numerical computing.
 extra-source-files:
   README.markdown
+  tests/*.hs
+  tests/Tests/*.hs
   tests/Tests/SpecFunctions/gen.py
 
 library
diff --git a/tests/Tests/Chebyshev.hs b/tests/Tests/Chebyshev.hs
new file mode 100644
--- /dev/null
+++ b/tests/Tests/Chebyshev.hs
@@ -0,0 +1,44 @@
+module Tests.Chebyshev (
+  tests
+  ) where
+
+import Data.Vector.Unboxed                  (fromList)
+import Test.Framework
+import Test.Framework.Providers.QuickCheck2
+import Test.QuickCheck                      (Arbitrary(..))
+
+import Tests.Helpers
+import Numeric.Polynomial.Chebyshev
+
+
+tests :: Test
+tests = testGroup "Chebyshev polynomials"
+  [ testProperty "Chebyshev 0" $ \a0 (Ch x) ->
+      (ch0 x * a0) ≈ (chebyshev x $ fromList [a0])
+  , testProperty "Chebyshev 1" $ \a0 a1 (Ch x) ->
+      (a0*ch0 x + a1*ch1 x) ≈  (chebyshev x $ fromList [a0,a1])
+  , testProperty "Chebyshev 2" $ \a0 a1 a2 (Ch x) ->
+       (a0*ch0 x + a1*ch1 x + a2*ch2 x) ≈ (chebyshev x $ fromList [a0,a1,a2])
+  , testProperty "Chebyshev 3" $ \a0 a1 a2 a3 (Ch x) ->
+       (a0*ch0 x + a1*ch1 x + a2*ch2 x + a3*ch3 x) ≈ (chebyshev x $ fromList [a0,a1,a2,a3])
+  , testProperty "Chebyshev 4" $ \a0 a1 a2 a3 a4 (Ch x) ->
+       (a0*ch0 x + a1*ch1 x + a2*ch2 x + a3*ch3 x + a4*ch4 x) ≈ (chebyshev x $ fromList [a0,a1,a2,a3,a4])
+  ]
+  where (≈) = eq 1e-12
+
+
+-- Chebyshev polynomials of low order
+ch0,ch1,ch2,ch3,ch4 :: Double -> Double
+ch0 _ = 1
+ch1 x = x
+ch2 x = 2*x^2 - 1
+ch3 x = 4*x^3 - 3*x
+ch4 x = 8*x^4 - 8*x^2 + 1
+
+
+-- Double in the [-1 .. 1] range
+newtype Ch = Ch Double
+             deriving Show
+instance Arbitrary Ch  where
+  arbitrary = do x <- arbitrary
+                 return $ Ch $ 2 * (snd . properFraction) x - 1
diff --git a/tests/Tests/Helpers.hs b/tests/Tests/Helpers.hs
new file mode 100644
--- /dev/null
+++ b/tests/Tests/Helpers.hs
@@ -0,0 +1,94 @@
+-- | Helpers for testing
+module Tests.Helpers (
+    -- * helpers
+    T(..)
+  , typeName
+  , eq
+  , eqC
+    -- * Generic QC tests
+  , monotonicallyIncreases
+  , monotonicallyIncreasesIEEE
+    -- * HUnit helpers
+  , testAssertion
+  , testEquality
+  ) where
+
+import Data.Complex
+import Data.Typeable
+
+import qualified Numeric.IEEE    as IEEE
+
+import qualified Test.HUnit      as HU
+import Test.Framework
+import Test.Framework.Providers.HUnit
+
+
+
+
+----------------------------------------------------------------
+-- Helpers
+----------------------------------------------------------------
+
+-- | Phantom typed value used to select right instance in QC tests
+data T a = T
+
+-- | String representation of type name
+typeName :: Typeable a => T a -> String
+typeName = show . typeOf . typeParam
+  where
+    typeParam :: T a -> a
+    typeParam _ = undefined
+
+-- | Approximate equality for 'Double'. Doesn't work well for numbers
+--   which are almost zero.
+eq :: Double                    -- ^ Relative error
+   -> Double -> Double -> Bool
+eq eps a b 
+  | a == 0 && b == 0 = True
+  | otherwise        = abs (a - b) <= eps * max (abs a) (abs b)
+
+-- | Approximate equality for 'Complex Double'
+eqC :: Double                   -- ^ Relative error
+    -> Complex Double
+    -> Complex Double
+    -> Bool
+eqC eps a@(ar :+ ai) b@(br :+ bi)
+  | a == 0 && b == 0 = True
+  | otherwise        = abs (ar - br) <= eps * d
+                    && abs (ai - bi) <= eps * d
+  where
+    d = max (realPart $ abs a) (realPart $ abs b)
+
+
+
+----------------------------------------------------------------
+-- Generic QC
+----------------------------------------------------------------
+
+-- Check that function is nondecreasing
+monotonicallyIncreases :: (Ord a, Ord b) => (a -> b) -> a -> a -> Bool
+monotonicallyIncreases f x1 x2 = f (min x1 x2) <= f (max x1 x2)
+
+-- Check that function is nondecreasing taking rounding errors into
+-- account.
+--
+-- In fact funstion is allowed to decrease less than one ulp in order
+-- to guard againist problems with excess precision. On x86 FPU works
+-- with 80-bit numbers but doubles are 64-bit so rounding happens
+-- whenever values are moved from registers to memory
+monotonicallyIncreasesIEEE :: (Ord a, IEEE.IEEE b)  => (a -> b) -> a -> a -> Bool
+monotonicallyIncreasesIEEE f x1 x2 =
+  y1 <= y2 || (y1 - y2) < y2 * IEEE.epsilon
+  where
+    y1 = f (min x1 x2)
+    y2 = f (max x1 x2)
+
+----------------------------------------------------------------
+-- HUnit helpers
+----------------------------------------------------------------
+
+testAssertion :: String -> Bool -> Test
+testAssertion str cont = testCase str $ HU.assertBool str cont
+
+testEquality :: (Show a, Eq a) => String -> a -> a -> Test
+testEquality msg a b = testCase msg $ HU.assertEqual msg a b
diff --git a/tests/Tests/SpecFunctions.hs b/tests/Tests/SpecFunctions.hs
new file mode 100644
--- /dev/null
+++ b/tests/Tests/SpecFunctions.hs
@@ -0,0 +1,127 @@
+{-# LANGUAGE ViewPatterns #-}
+-- | Tests for Statistics.Math
+module Tests.SpecFunctions (
+  tests
+  ) where
+
+import qualified Data.Vector as V
+import           Data.Vector   ((!))
+
+import Test.QuickCheck  hiding (choose)
+import Test.Framework
+import Test.Framework.Providers.QuickCheck2
+
+import Tests.Helpers
+import Tests.SpecFunctions.Tables
+import Numeric.SpecFunctions
+
+
+tests :: Test
+tests = testGroup "Special functions"
+  [ testProperty "Γ(x+1) = x·Γ(x) logGamma"  $ gammaReccurence logGamma  3e-8
+  , testProperty "Γ(x+1) = x·Γ(x) logGammaL" $ gammaReccurence logGammaL 2e-13
+  , testProperty "γ(1,x) = 1 - exp(-x)"      $ incompleteGammaAt1Check
+  , testProperty "γ - increases"             $
+      \s x y -> s > 0 && x > 0 && y > 0 ==> monotonicallyIncreases (incompleteGamma s) x y
+  , testProperty "invIncompleteGamma = γ^-1" $ invIGammaIsInverse
+  , testProperty "invIncompleteBeta  = B^-1" $ invIBetaIsInverse
+    -- Unit tests
+  , testAssertion "Factorial is expected to be precise at 1e-15 level"
+      $ and [ eq 1e-15 (factorial (fromIntegral n))
+                       (fromIntegral (factorial' n))
+            |n <- [0..170]]
+  , testAssertion "Log factorial is expected to be precise at 1e-15 level"
+      $ and [ eq 1e-15 (logFactorial (fromIntegral n))
+                       (log $ fromIntegral $ factorial' n)
+            | n <- [2..170]]
+  , testAssertion "logGamma is expected to be precise at 1e-9 level [integer points]"
+      $ and [ eq 1e-9 (logGamma (fromIntegral n))
+                      (logFactorial (n-1))
+            | n <- [3..10000]]
+  , testAssertion "logGamma is expected to be precise at 1e-9 level [fractional points]"
+      $ and [ eq 1e-9 (logGamma x) lg | (x,lg) <- tableLogGamma ]
+  , testAssertion "logGammaL is expected to be precise at 1e-15 level"
+      $ and [ eq 1e-15 (logGammaL (fromIntegral n))
+                       (logFactorial (n-1))
+            | n <- [3..10000]]
+  , testAssertion "logGammaL is expected to be precise at 1e-9 level [fractional points]"
+      $ and [ eq 1e-10 (logGammaL x) lg | (x,lg) <- tableLogGamma ]
+  , testAssertion "logBeta is expected to be precise at 1e-6 level"
+      $ and [ eq 1e-6 (logBeta p q)
+                      (logGammaL p + logGammaL q - logGammaL (p+q))
+            | p <- [0.1,0.2 .. 0.9] ++ [2 .. 20]
+            , q <- [0.1,0.2 .. 0.9] ++ [2 .. 20]]
+  -- 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 "choose is expected to precise at 1e-12 level"
+      $ and [ eq 1e-12 (choose (fromIntegral n) (fromIntegral k)) (fromIntegral $ choose' n k)
+            | n <- [0..300], k <- [0..n]]
+  ]
+
+----------------------------------------------------------------
+-- QC tests
+----------------------------------------------------------------
+
+-- Γ(x+1) = x·Γ(x)
+gammaReccurence :: (Double -> Double) -> Double -> Double -> Property
+gammaReccurence logG ε x =
+  (x > 0 && x < 100)  ==>  (abs (g2 - g1 - log x) < ε)
+    where
+      g1 = logG x
+      g2 = logG (x+1)
+
+
+-- γ(1,x) = 1 - exp(-x)
+-- Since Γ(1) = 1 normalization doesn't make any difference
+incompleteGammaAt1Check :: Double -> Property
+incompleteGammaAt1Check x =
+  x > 0 ==> (incompleteGamma 1 x + exp(-x)) ≈ 1
+  where
+    (≈) = eq 1e-13
+
+-- invIncompleteGamma is inverse of incompleteGamma
+invIGammaIsInverse :: Double -> Double -> Property
+invIGammaIsInverse (abs -> a) (abs . snd . properFraction -> p) =
+  a > 0 && p > 0 && p < 1  ==> ( printTestCase ("x  = " ++ show x )
+                               $ printTestCase ("p' = " ++ show p')
+                               $ printTestCase ("Δp = " ++ show (p - p'))
+                               $ abs (p - p') <= 1e-12
+                               )
+  where
+    x  = invIncompleteGamma a p
+    p' = incompleteGamma    a x
+
+-- invIncompleteBeta is inverse of incompleteBeta
+invIBetaIsInverse :: Double -> Double -> Double -> Property
+invIBetaIsInverse (abs -> p) (abs -> q) (abs . snd . properFraction -> x) =
+  p > 0 && q > 0  ==> ( printTestCase ("p   = " ++ show p )
+                      $ printTestCase ("q   = " ++ show q )
+                      $ printTestCase ("x   = " ++ show x )
+                      $ printTestCase ("x'  = " ++ show x')
+                      $ printTestCase ("a   = " ++ show a)  
+                      $ printTestCase ("err = " ++ (show $ abs $ (x - x') / x))
+                      $ abs (x - x') <= 1e-12
+                      )
+  where
+    x' = incompleteBeta    p q a
+    a  = invIncompleteBeta p q x
+  
+
+
+----------------------------------------------------------------
+-- Unit tests
+----------------------------------------------------------------
+
+-- Lookup table for fact factorial calculation. It has fixed size
+-- which is bad but it's OK for this particular case
+factorial_table :: V.Vector Integer
+factorial_table = V.generate 2000 (\n -> product [1..fromIntegral n])
+
+-- Exact implementation of factorial
+factorial' :: Integer -> Integer
+factorial' n = factorial_table ! fromIntegral n
+
+-- Exact albeit slow implementation of choose
+choose' :: Integer -> Integer -> Integer
+choose' n k = factorial' n `div` (factorial' k * factorial' (n-k))
diff --git a/tests/tests.hs b/tests/tests.hs
new file mode 100644
--- /dev/null
+++ b/tests/tests.hs
@@ -0,0 +1,8 @@
+import Test.Framework       (defaultMain)
+import qualified Tests.SpecFunctions
+import qualified Tests.Chebyshev
+
+main :: IO ()
+main = defaultMain [ Tests.SpecFunctions.tests
+                   , Tests.Chebyshev.tests
+                   ]