diff --git a/Changes b/Changes
--- a/Changes
+++ b/Changes
@@ -1,3 +1,18 @@
+0.4.2.0:
+    This release supports GHC 7.6, 7.8 and 8.0.
+
+    Add new cabal flag check-bounds, which replaces all unsafe array functions with safe ones.
+
+    Add basic functions on Gaussian integers.
+    Add Moebius mu-function.
+
+    Forbid non-positive moduli in Math.NumberTheory.Moduli.
+
+    Fix out-of-bounds error in Math.NumberTheory.Primes.Heap, Math.NumberTheory.Primes.Sieve and Math.NumberTheory.MoebiusInversion.
+    Fix 32-bit build.
+    Fix binaryGCD on negative numbers.
+    Fix highestPower (various issues).
+
 0.4.1.0:
     Add integerLog10 variants at Bas van Dijk's request and expose
     Math.NumberTheory.Powers.Integer, with an added integerWordPower.
diff --git a/Math/NumberTheory/GCD.hs b/Math/NumberTheory/GCD.hs
--- a/Math/NumberTheory/GCD.hs
+++ b/Math/NumberTheory/GCD.hs
@@ -83,10 +83,13 @@
 binaryGCDImpl a 0 = abs a
 binaryGCDImpl 0 b = abs b
 binaryGCDImpl a b =
-    case shiftToOddCount a of
+    case shiftToOddCount a' of
       (!za, !oa) ->
-        case shiftToOddCount b of
+        case shiftToOddCount b' of
           (!zb, !ob) -> gcdOdd (abs oa) (abs ob) `shiftL` min za zb
+    where
+      a' = abs a
+      b' = abs b
 
 {-# SPECIALISE extendedGCD :: Int -> Int -> (Int, Int, Int),
                               Word -> Word -> (Word, Word, Word),
@@ -95,16 +98,21 @@
 -- | Calculate the greatest common divisor of two numbers and coefficients
 --   for the linear combination.
 --
---   Satisfies:
+--   For signed types satisfies:
 --
 -- > case extendedGCD a b of
 -- >   (d, u, v) -> u*a + v*b == d
--- >
--- > d == gcd a b
+-- >                && d == gcd a b
 --
---   and, for signed types,
+--   For unsigned and bounded types the property above holds, but since @u@ and @v@ must also be unsigned,
+--   the result may look weird. E. g., on 64-bit architecture
 --
--- >
+-- > extendedGCD (2 :: Word) (3 :: Word) == (1, 2^64-1, 1)
+--
+--   For unsigned and unbounded types (like 'Numeric.Natural.Natural') the result is undefined.
+--
+--   For signed types we also have
+--
 -- > abs u < abs b || abs b <= 1
 -- >
 -- > abs v < abs a || abs a <= 1
diff --git a/Math/NumberTheory/GaussianIntegers.hs b/Math/NumberTheory/GaussianIntegers.hs
new file mode 100644
--- /dev/null
+++ b/Math/NumberTheory/GaussianIntegers.hs
@@ -0,0 +1,238 @@
+-- |
+-- Module:      Math.NumberTheory.GaussianIntegers
+-- Copyright:   (c) 2016 Chris Fredrickson
+-- Licence:     MIT
+-- Maintainer:  Chris Fredrickson <chris.p.fredrickson@gmail.com>
+-- Stability:   Provisional
+-- Portability: Non-portable (GHC extensions)
+--
+-- This module exports functions for manipulating Gaussian integers, including
+-- computing their prime factorisations.
+--
+
+{-# LANGUAGE BangPatterns #-}
+module Math.NumberTheory.GaussianIntegers (
+    GaussianInteger((:+)),
+    ι,
+    real,
+    imag,
+    conjugate,
+    norm,
+    divModG,
+    divG,
+    modG,
+    quotRemG,
+    quotG,
+    remG,
+    (.^),
+    isPrime,
+    primes,
+    gcdG,
+    gcdG',
+    findPrime,
+    findPrime',
+    factorise,
+) where
+
+import qualified Math.NumberTheory.Moduli as Moduli
+import qualified Math.NumberTheory.Powers as Powers
+import qualified Math.NumberTheory.Primes.Factorisation as Factorisation
+import qualified Math.NumberTheory.Primes.Sieve as Sieve
+import qualified Math.NumberTheory.Primes.Testing as Testing
+
+infix 6 :+
+infixr 8 .^
+-- |A Gaussian integer is a+bi, where a and b are both integers.
+data GaussianInteger = (:+) { real :: !Integer, imag :: !Integer } deriving (Eq)
+
+-- |The imaginary unit, where
+--
+-- > ι .^ 2 == -1
+ι :: GaussianInteger
+ι = 0 :+ 1
+
+instance Show GaussianInteger where
+    show (a :+ b)
+        | b == 0     = show a
+        | a == 0     = s ++ b'
+        | otherwise  = show a ++ op ++ b'
+        where
+            b' = if abs b == 1 then "ι" else show (abs b) ++ "*ι"
+            op = if b > 0 then "+" else "-"
+            s  = if b > 0 then "" else "-"
+
+instance Num GaussianInteger where
+    (+) (a :+ b) (c :+ d) = (a + c) :+ (b + d)
+    (*) (a :+ b) (c :+ d) = (a * c - b * d) :+ (a * d + b * c)
+    abs z@(a :+ b)
+        | a == 0 && b == 0 =   z             -- origin
+        | a >  0 && b >= 0 =   z             -- first quadrant: (0, inf) x [0, inf)i
+        | a <= 0 && b >  0 =   b  :+ (-a)    -- second quadrant: (-inf, 0] x (0, inf)i
+        | a <  0 && b <= 0 = (-a) :+ (-b)    -- third quadrant: (-inf, 0) x (-inf, 0]i
+        | otherwise        = (-b) :+   a     -- fourth quadrant: [0, inf) x (-inf, 0)i
+    negate (a :+ b) = (-a) :+ (-b)
+    fromInteger n = n :+ 0
+    signum z@(a :+ b)
+        | a == 0 && b == 0 = z               -- hole at origin
+        | otherwise        = z `divG` abs z
+
+-- |Simultaneous 'quot' and 'rem'.
+quotRemG :: GaussianInteger -> GaussianInteger -> (GaussianInteger, GaussianInteger)
+quotRemG = divHelper quot
+
+-- |Gaussian integer division, truncating toward zero.
+quotG :: GaussianInteger -> GaussianInteger -> GaussianInteger
+n `quotG` d = q where (q,_) = quotRemG n d
+
+-- |Gaussian integer remainder, satisfying
+--
+-- > (x `quotG` y)*y + (x `remG` y) == x
+remG :: GaussianInteger -> GaussianInteger -> GaussianInteger
+n `remG`  d = r where (_,r) = quotRemG n d
+
+-- |Simultaneous 'div' and 'mod'.
+divModG :: GaussianInteger -> GaussianInteger -> (GaussianInteger, GaussianInteger)
+divModG = divHelper div
+
+-- |Gaussian integer division, truncating toward negative infinity.
+divG :: GaussianInteger -> GaussianInteger -> GaussianInteger
+n `divG` d = q where (q,_) = divModG n d
+
+-- |Gaussian integer remainder, satisfying
+--
+-- > (x `divG` y)*y + (x `modG` y) == x
+modG :: GaussianInteger -> GaussianInteger -> GaussianInteger
+n `modG` d = r where (_,r) = divModG n d
+
+divHelper :: (Integer -> Integer -> Integer) -> GaussianInteger -> GaussianInteger -> (GaussianInteger, GaussianInteger)
+divHelper divide g h =
+    let nr :+ ni = g * conjugate h
+        denom = norm h
+        q = divide nr denom :+ divide ni denom
+        p = h * q
+    in (q, g - p)
+
+-- |Conjugate a Gaussian integer.
+conjugate :: GaussianInteger -> GaussianInteger
+conjugate (r :+ i) = r :+ (-i)
+
+-- |The square of the magnitude of a Gaussian integer.
+norm :: GaussianInteger -> Integer
+norm (x :+ y) = x * x + y * y
+
+-- |Compute whether a given Gaussian integer is prime.
+isPrime :: GaussianInteger -> Bool
+isPrime g@(x :+ y)
+    | x == 0 && y /= 0 = abs y `mod` 4 == 3 && Testing.isPrime y
+    | y == 0 && x /= 0 = abs x `mod` 4 == 3 && Testing.isPrime x
+    | otherwise        = Testing.isPrime $ norm g
+
+-- |An infinite list of the Gaussian primes. Uses primes in Z to exhaustively
+-- generate all Gaussian primes, but not quite in order of ascending magnitude.
+primes :: [GaussianInteger]
+primes = [ g
+         | p <- Sieve.primes
+         , g <- if p `mod` 4 == 3
+                then [p :+ 0]
+                else
+                    if p == 2
+                    then [1 :+ 1]
+                    else let x :+ y = findPrime' p
+                         in [x :+ y, y :+ x]
+         ]
+
+-- |Compute the GCD of two Gaussian integers. Enforces the precondition that each
+-- integer must be in the first quadrant (or zero).
+gcdG :: GaussianInteger -> GaussianInteger -> GaussianInteger
+gcdG g h = gcdG' (abs g) (abs h)
+
+-- |Compute the GCD of two Gauss integers. Does not check the precondition.
+gcdG' :: GaussianInteger -> GaussianInteger -> GaussianInteger
+gcdG' g h
+    | h == 0    = g --done recursing
+    | otherwise = gcdG' h (abs (g `modG` h))
+
+-- |Find a Gaussian integer whose norm is the given prime number.
+-- Checks the precondition that p is prime and that p `mod` 4 /= 3.
+findPrime :: Integer -> GaussianInteger
+findPrime p
+    | p == 2 || (p `mod` 4 == 1 && Testing.isPrime p) = findPrime' p
+    | otherwise = error "p must be prime, and not congruent to 3 (mod 4)"
+
+-- |Find a Gaussian integer whose norm is the given prime number. Does not
+-- check the precondition.
+findPrime' :: Integer -> GaussianInteger
+findPrime' p =
+    let (Just c) = Moduli.sqrtModP (-1) p
+        k  = Powers.integerSquareRoot p
+        bs = [1 .. k]
+        asbs = map (\b' -> ((b' * c) `mod` p, b')) bs
+        (a, b) = head [ (a', b') | (a', b') <- asbs, a' <= k]
+    in a :+ b
+
+-- |Raise a Gaussian integer to a given power.
+(.^) :: (Integral a) => GaussianInteger -> a -> GaussianInteger
+a .^ e
+    | e < 0 && norm a == 1 =
+        case a of
+            1    :+ 0 -> 1
+            (-1) :+ 0 -> if even e then 1 else (-1)
+            0    :+ 1 -> (0 :+ (-1)) .^ (abs e `mod` 4)
+            _         -> (0 :+ 1) .^ (abs e `mod` 4)
+    | e < 0     = error "Cannot exponentiate non-unit Gaussian Int to negative power"
+    | a == 0    = 0
+    | e == 0    = 1
+    | even e    = s * s
+    | otherwise = a * a .^ (e - 1)
+    where
+    s = a .^ div e 2
+
+-- |Compute the prime factorization of a Gaussian integer. This is unique up to units (+/- 1, +/- i).
+factorise :: GaussianInteger -> [(GaussianInteger, Int)]
+factorise g
+    | g == 0    = error "0 has no prime factorisation"
+    | g == 1    = []
+    | otherwise =
+        let helper :: [(Integer, Int)] -> GaussianInteger -> [(GaussianInteger, Int)] -> [(GaussianInteger, Int)]
+            helper [] g' fs = (if g' == 1 then [] else [(g', 1)]) ++ fs    -- include the unit, if it isn't 1
+            helper ((!p, !e) : pt) g' fs
+                | p `mod` 4 == 3 =
+                    -- prime factors congruent to 3 mod 4 are simple.
+                    let pow = div e 2
+                        gp = fromInteger p
+                    in helper pt (g' `divG` (gp .^ pow)) ((gp, pow) : fs)
+                | otherwise      =
+                    -- general case: for every prime factor of the magnitude
+                    -- squared, find a Gaussian prime whose magnitude squared
+                    -- is that prime. Then find out how many times the original
+                    -- number is divisible by that Gaussian prime and its
+                    -- conjugate. The order that the prime factors are
+                    -- processed doesn't really matter, but it is reversed so
+                    -- that the Gaussian primes will be in order of increasing
+                    -- magnitude.
+                    let gp = findPrime' p
+                        (!gNext, !facs) = trialDivide g' [gp, abs $ conjugate gp] []
+                    in helper pt gNext (facs ++ fs)
+        in helper (reverse . Factorisation.factorise $ norm g) g []
+
+-- Divide a Gaussian integer by a set of (relatively prime) Gaussian integers,
+-- as many times as possible, and return the final quotient as well as a count
+-- of how many times each factor divided the original.
+trialDivide :: GaussianInteger -> [GaussianInteger] -> [(GaussianInteger, Int)] -> (GaussianInteger, [(GaussianInteger, Int)])
+trialDivide g [] fs = (g, fs)
+trialDivide g (pf : pft) fs
+    | g `modG` pf == 0 =
+        let (cnt, g') = countEvenDivisions g pf
+        in trialDivide g' pft ((pf, cnt) : fs)
+    | otherwise    = trialDivide g pft fs
+
+-- Divide a Gaussian integer by a possible factor, and return how many times
+-- the factor divided it evenly, as well as the result of dividing the original
+-- that many times.
+countEvenDivisions :: GaussianInteger -> GaussianInteger -> (Int, GaussianInteger)
+countEvenDivisions g pf = helper g 0
+    where
+    helper :: GaussianInteger -> Int -> (Int, GaussianInteger)
+    helper g' acc
+        | g' `modG` pf == 0 = helper (g' `divG` pf) (1 + acc)
+        | otherwise     = (acc, g')
diff --git a/Math/NumberTheory/Logarithms.hs b/Math/NumberTheory/Logarithms.hs
--- a/Math/NumberTheory/Logarithms.hs
+++ b/Math/NumberTheory/Logarithms.hs
@@ -36,10 +36,10 @@
 
 import Data.Bits
 import Data.Array.Unboxed
-import Data.Array.Base (unsafeAt)
 
 import Math.NumberTheory.Logarithms.Internal
 import Math.NumberTheory.Powers.Integer
+import Math.NumberTheory.Unsafe
 #if __GLASGOW_HASKELL__ < 707
 import Math.NumberTheory.Utils  (isTrue#)
 #endif
diff --git a/Math/NumberTheory/Moduli.hs b/Math/NumberTheory/Moduli.hs
--- a/Math/NumberTheory/Moduli.hs
+++ b/Math/NumberTheory/Moduli.hs
@@ -34,52 +34,50 @@
 
 #include "MachDeps.h"
 
-#if __GLASGOW_HASKELL__ < 709
+#if __GLASGOW_HASKELL__ < 709 || WORD_SIZE_IN_BITS == 32
 import Data.Word
 #endif
 import Data.Bits
 import Data.Array.Unboxed
-import Data.Array.Base (unsafeAt)
-import Data.Maybe (fromJust)
 import Data.List (nub)
 import Control.Monad (foldM, liftM2)
 
 import Math.NumberTheory.Utils (shiftToOddCount, splitOff)
 import Math.NumberTheory.GCD (extendedGCD)
 import Math.NumberTheory.Primes.Heap (sieveFrom)
+import Math.NumberTheory.Unsafe
+
 -- Guesstimated startup time for the Heap algorithm is lower than
 -- the cost to sieve an entire chunk.
 
--- | Invert a number relative to a modulus.
+-- | Invert a number relative to a positive modulus.
 --   If @number@ and @modulus@ are coprime, the result is
 --   @Just inverse@ where
 --
--- >    (number * inverse) `mod` (abs modulus) == 1
--- >    0 <= inverse < abs modulus
+-- >    (number * inverse) `mod` modulus == 1
+-- >    0 <= inverse < modulus
 --
---   unless @modulus == 0@ and @abs number == 1@, in which case the
---   result is @Just number@.
---   If @gcd number modulus > 1@, the result is @Nothing@.
+--   If @number `mod` modulus == 0@ or @gcd number modulus > 1@, the result is @Nothing@.
 invertMod :: Integer -> Integer -> Maybe Integer
-invertMod k 0 = if k == 1 || k == (-1) then Just k else Nothing
-invertMod k m = wrap $ go False 1 0 m' k'
+invertMod k m
+  | m <= 0 = error "Math.NumberTheory.Moduli.invertMod: non-positive modulus"
+  | otherwise = wrap $ go False 1 0 m k'
   where
-    m' = abs m
-    k' | r < 0     = r+m'
+    k' | r < 0     = r+m
        | otherwise = r
          where
-           r = k `rem` m'
-    wrap x = case (x*k') `rem` m' of
+           r = k `rem` m
+    wrap x = case (x*k') `rem` m of
                1 -> Just x
                _ -> Nothing
-    -- Calculate modular inverse of k' modulo m' by continued fraction expansion
-    -- of m'/k', say [a_0,a_1,...,a_s]. Let the convergents be p_j/q_j.
+    -- Calculate modular inverse of k' modulo m by continued fraction expansion
+    -- of m/k', say [a_0,a_1,...,a_s]. Let the convergents be p_j/q_j.
     -- Starting from j = -2, the arguments of go are
-    -- (p_j/q_j) > m'/k', p_{j+1}, p_j, and n, d with n/d = [a_{j+2},...,a_s].
-    -- Since m'/k' = p_s/q_s, and p_j*q_{j+1} - p_{j+1}*q_j = (-1)^(j+1), we have
-    -- p_{s-1}*k' - q_{s-1}*m' = (-1)^s * gcd m' k', so if the inverse exists,
+    -- (p_j/q_j) > m/k', p_{j+1}, p_j, and n, d with n/d = [a_{j+2},...,a_s].
+    -- Since m/k' = p_s/q_s, and p_j*q_{j+1} - p_{j+1}*q_j = (-1)^(j+1), we have
+    -- p_{s-1}*k' - q_{s-1}*m = (-1)^s * gcd m k', so if the inverse exists,
     -- it is either p_{s-1} or -p_{s-1}, depending on whether s is even or odd.
-    go !b _ po _ 0 = if b then po else (m'-po)
+    go !b _ po _ 0 = if b then po else (m-po)
     go b !pn po n d = case n `quotRem` d of
                         (q,r) -> go (not b) (q*pn+po) pn d r
 
@@ -161,7 +159,7 @@
 --
 -- > powerMod base exponent modulus
 --
---   calculates @(base ^ exponent) \`mod\` modulus@ by repeated squaring and reduction.
+--   calculates @(base ^ exponent) \`mod\` modulus@ by repeated squaring and reduction. Modulus must be positive.
 --   If @exponent < 0@ and @base@ is invertible modulo @modulus@, @(inverse ^ |exponent|) \`mod\` modulus@
 --   is calculated. This function does some input checking and sanitation before calling the unsafe worker.
 {-# RULES
@@ -176,18 +174,17 @@
   #-}
 powerModImpl :: (Integral a, Bits a) => Integer -> a -> Integer -> Integer
 powerModImpl base expo md
-  | md == 0     = base ^ expo
-  | md' == 1    = 0
+  | md <= 0     = error "Math.NumberTheory.Moduli.powerMod: non-positive modulus"
+  | md == 1     = 0
   | expo == 0   = 1
   | bse' == 1   = 1
-  | expo < 0    = case invertMod bse' md' of
-                    Just i  -> powerMod'Impl i (negate expo) md'
+  | expo < 0    = case invertMod bse' md of
+                    Just i  -> powerMod'Impl i (negate expo) md
                     Nothing -> error "Math.NumberTheory.Moduli.powerMod: Base isn't invertible with respect to modulus"
   | bse' == 0   = 0
-  | otherwise   = powerMod'Impl bse' expo md'
+  | otherwise   = powerMod'Impl bse' expo md
     where
-      md' = abs md
-      bse' = if base < 0 || md' <= base then base `mod` md' else base
+      bse' = if base < 0 || md <= base then base `mod` md else base
 
 -- | Modular power worker without input checking.
 --   Assumes all arguments strictly positive and modulus greater than 1.
@@ -213,21 +210,20 @@
 -- | Specialised version of 'powerMod' for 'Integer' exponents.
 --   Reduces the number of shifts of the exponent since shifting
 --   large 'Integer's is expensive. Call this function directly
---   if you don't want or can't rely on rewrite rules.
+--   if you don't want or can't rely on rewrite rules. Modulus must be positive.
 powerModInteger :: Integer -> Integer -> Integer -> Integer
 powerModInteger base ex mdl
-  | mdl == 0    = base ^ ex
-  | mdl' == 1   = 0
+  | mdl <= 0     = error "Math.NumberTheory.Moduli.powerModInteger: non-positive modulus"
+  | mdl == 1    = 0
   | ex == 0     = 1
-  | ex < 0      = case invertMod bse' mdl' of
-                    Just i  -> powerModInteger' i (negate ex) mdl'
+  | ex < 0      = case invertMod bse' mdl of
+                    Just i  -> powerModInteger' i (negate ex) mdl
                     Nothing -> error "Math.NumberTheory.Moduli.powerMod: Base isn't invertible with respect to modulus"
   | bse' == 0   = 0
   | bse' == 1   = 1
-  | otherwise   = powerModInteger' bse' ex mdl'
+  | otherwise   = powerModInteger' bse' ex mdl
     where
-      mdl' = abs mdl
-      bse' = if base < 0 || mdl' <= base then base `mod` mdl' else base
+      bse' = if base < 0 || mdl <= base then base `mod` mdl else base
 
 -- | Specialised worker without input checks. Makes the same assumptions
 --   as the general version 'powerMod''.
@@ -322,7 +318,7 @@
                         _      -> []
 
 -- | @sqrtModP' square prime@ finds a square root of @square@ modulo
---   prime. @prime@ /must/ be a (positive) prime, and @sqaure@ /must/ be a
+--   prime. @prime@ /must/ be a (positive) prime, and @square@ /must/ be a positive
 --   quadratic residue modulo @prime@, i.e. @'jacobi square prime == 1@.
 --   The precondition is /not/ checked.
 sqrtModP' :: Integer -> Integer -> Integer
@@ -333,7 +329,7 @@
 
 -- | @tonelliShanks square prime@ calculates a square root of @square@
 --   modulo @prime@, where @prime@ is a prime of the form @4*k + 1@ and
---   @square@ is a quadratic residue modulo @prime@, using the
+--   @square@ is a positive quadratic residue modulo @prime@, using the
 --   Tonelli-Shanks algorithm.
 --   No checks on the input are performed.
 tonelliShanks :: Integer -> Integer -> Integer
@@ -366,15 +362,15 @@
 sqrtModPP :: Integer -> (Integer,Int) -> Maybe Integer
 sqrtModPP n (2,e) = sqM2P n e
 sqrtModPP n (prime,expo) = case sqrtModP n prime of
-                             Just r -> Just $ fixup r
+                             Just r -> fixup r
                              _      -> Nothing
   where
     fixup r = let diff' = r*r-n
               in if diff' == 0
-                   then r
+                   then Just r
                    else case splitOff prime diff' of
-                          (e,q) | expo <= e -> r
-                                | otherwise -> hoist (fromJust $ invertMod (2*r) prime) r (q `mod` prime) (prime^e)
+                          (e,q) | expo <= e -> Just r
+                                | otherwise -> fmap (\inv -> hoist inv r (q `mod` prime) (prime^e)) (invertMod (2*r) prime)
                       --
     hoist inv root elim pp
         | diff' == 0    = root'
@@ -418,13 +414,17 @@
 -- | @sqrtModF n primePowers@ calculates a square root of @n@ modulo
 --   @product [p^k | (p,k) <- primePowers]@ if one exists and all primes
 --   are distinct.
+--   The list must be non-empty, @n@ must be coprime with all primes.
 sqrtModF :: Integer -> [(Integer,Int)] -> Maybe Integer
+sqrtModF _ []  = Nothing
 sqrtModF n pps = do roots <- mapM (sqrtModPP n) pps
                     chineseRemainder $ zip roots (map (uncurry (^)) pps)
 
 -- | @sqrtModFList n primePowers@ calculates all square roots of @n@ modulo
 --   @product [p^k | (p,k) <- primePowers]@ if all primes are distinct.
+--   The list must be non-empty, @n@ must be coprime with all primes.
 sqrtModFList :: Integer -> [(Integer,Int)] -> [Integer]
+sqrtModFList _ []  = []
 sqrtModFList n pps = map fst $ foldl1 (liftM2 comb) cs
   where
     ms :: [Integer]
@@ -439,6 +439,7 @@
 --   square roots of @n@ modulo @prime^expo@. The same restriction
 --   as in 'sqrtModPP' applies to the arguments.
 sqrtModPPList :: Integer -> (Integer,Int) -> [Integer]
+sqrtModPPList n (2,1) = [n `mod` 2]
 sqrtModPPList n (2,expo)
     = case sqM2P n expo of
         Just r -> let m = 1 `shiftL` (expo-1)
@@ -458,13 +459,14 @@
 -- > r ≡ r_k (mod m_k)
 -- >
 --
---   if all moduli are pairwise coprime. If not all moduli are
---   pairwise coprime, the result is @Nothing@ regardless of whether
+--   if all moduli are positive and pairwise coprime. Otherwise
+--   the result is @Nothing@ regardless of whether
 --   a solution exists.
 chineseRemainder :: [(Integer,Integer)] -> Maybe Integer
 chineseRemainder remainders = foldM addRem 0 remainders
   where
     !modulus = product (map snd remainders)
+    addRem acc (_,1) = Just acc
     addRem acc (r,m) = do
         let cf = modulus `quot` m
         inv <- invertMod cf m
@@ -489,21 +491,21 @@
                         Int -> Bool,
                         Word -> Bool
   #-}
-evenI :: (Integral a, Bits a) => a -> Bool
+evenI :: Integral a => a -> Bool
 evenI n = fromIntegral n .&. 1 == (0 :: Int)
 
 {-# SPECIALISE rem4 :: Integer -> Int,
                        Int -> Int,
                        Word -> Int
   #-}
-rem4 :: (Integral a, Bits a) => a -> Int
+rem4 :: Integral a => a -> Int
 rem4 n = fromIntegral n .&. 3
 
 {-# SPECIALISE rem8 :: Integer -> Int,
                        Int -> Int,
                        Word -> Int
   #-}
-rem8 :: (Integral a, Bits a) => a -> Int
+rem8 :: Integral a => a -> Int
 rem8 n = fromIntegral n .&. 7
 
 jac2 :: UArray Int Int
diff --git a/Math/NumberTheory/MoebiusInversion.hs b/Math/NumberTheory/MoebiusInversion.hs
--- a/Math/NumberTheory/MoebiusInversion.hs
+++ b/Math/NumberTheory/MoebiusInversion.hs
@@ -15,10 +15,11 @@
     ) where
 
 import Data.Array.ST
-import Data.Array.Base
+import Control.Monad
 import Control.Monad.ST
 
 import Math.NumberTheory.Powers.Squares
+import Math.NumberTheory.Unsafe
 
 -- | @totientSum n@ is, for @n > 0@, the sum of @[totient k | k <- [1 .. n]]@,
 --   computed via generalised Moebius inversion.
@@ -69,7 +70,7 @@
 --   Since the function arguments are used as array indices, the domain of
 --   @f@ is restricted to 'Int'.
 --
---   The value @f n@ is then computed by @generalInversion g n@). Note that when
+--   The value @f n@ is then computed by @generalInversion g n@. Note that when
 --   many values of @f@ are needed, there are far more efficient methods, this
 --   method is only appropriate to compute isolated values of @f@.
 generalInversion :: (Int -> Integer) -> Int -> Integer
@@ -90,7 +91,8 @@
         small <- newArray_ (0,mk0) :: ST s (STArray s Int Integer)
         unsafeWrite small 0 0
         unsafeWrite small 1 $! (fun 1)
-        unsafeWrite small 2 $! (fun 2 - fun 1)
+        when (mk0 >= 2) $
+            unsafeWrite small 2 $! (fun 2 - fun 1)
         let calcit switch change i
                 | mk0 < i   = return (switch,change)
                 | i == change = calcit (switch+1) (change + 4*switch+6) i
diff --git a/Math/NumberTheory/MoebiusInversion/Int.hs b/Math/NumberTheory/MoebiusInversion/Int.hs
--- a/Math/NumberTheory/MoebiusInversion/Int.hs
+++ b/Math/NumberTheory/MoebiusInversion/Int.hs
@@ -16,11 +16,11 @@
     ) where
 
 import Data.Array.ST
-import Data.Array.Base
+import Control.Monad
 import Control.Monad.ST
 
 import Math.NumberTheory.Powers.Squares
-
+import Math.NumberTheory.Unsafe
 
 -- | @totientSum n@ is, for @n > 0@, the sum of @[totient k | k <- [1 .. n]]@,
 --   computed via generalised Moebius inversion.
@@ -71,11 +71,12 @@
 --   That bears the risk of overflow, however, so be sure to use it only when it's
 --   safe.
 --
---   The value @f n@ is then computed by @generalInversion g n@). Note that when
+--   The value @f n@ is then computed by @generalInversion g n@. Note that when
 --   many values of @f@ are needed, there are far more efficient methods, this
 --   method is only appropriate to compute isolated values of @f@.
 generalInversion :: (Int -> Int) -> Int -> Int
 generalInversion fun n
+    | n < 1     = error "Moebius inversion only defined on positive domain"
     | n == 1    = fun 1
     | n == 2    = fun 2 - fun 1
     | n == 3    = fun 3 - 2*fun 1
@@ -91,7 +92,8 @@
         small <- newArray_ (0,mk0) :: ST s (STUArray s Int Int)
         unsafeWrite small 0 0
         unsafeWrite small 1 (fun 1)
-        unsafeWrite small 2 (fun 2 - fun 1)
+        when (mk0 >= 2) $
+            unsafeWrite small 2 (fun 2 - fun 1)
         let calcit switch change i
                 | mk0 < i   = return (switch,change)
                 | i == change = calcit (switch+1) (change + 4*switch+6) i
diff --git a/Math/NumberTheory/Powers/Cubes.hs b/Math/NumberTheory/Powers/Cubes.hs
--- a/Math/NumberTheory/Powers/Cubes.hs
+++ b/Math/NumberTheory/Powers/Cubes.hs
@@ -21,7 +21,6 @@
 #include "MachDeps.h"
 
 import Data.Array.Unboxed
-import Data.Array.Base
 import Data.Array.ST
 
 import Data.Bits
@@ -34,6 +33,7 @@
 import GHC.Integer.GMP.Internals
 
 import Math.NumberTheory.Logarithms.Internal (integerLog2#)
+import Math.NumberTheory.Unsafe
 #if __GLASGOW_HASKELL__ < 707
 import Math.NumberTheory.Utils (isTrue#)
 #endif
@@ -61,9 +61,9 @@
 --   that is, the largest integer @r@ such that @r^3 <= n@.
 --   The precondition @n >= 0@ is not checked.
 {-# RULES
-"integerCubeRoot'/Int"  integerCubeRoot' = cubeRootInt'
-"integerCubeRoot'/Word" integerCubeRoot' = cubeRootWord
-"integerCubeRoot'/Igr"  integerCubeRoot' = cubeRootIgr
+"integerCubeRoot'/Int"     integerCubeRoot' = cubeRootInt'
+"integerCubeRoot'/Word"    integerCubeRoot' = cubeRootWord
+"integerCubeRoot'/Integer" integerCubeRoot' = cubeRootIgr
   #-}
 {-# INLINE [1] integerCubeRoot' #-}
 integerCubeRoot' :: Integral a => a -> a
diff --git a/Math/NumberTheory/Powers/Fourth.hs b/Math/NumberTheory/Powers/Fourth.hs
--- a/Math/NumberTheory/Powers/Fourth.hs
+++ b/Math/NumberTheory/Powers/Fourth.hs
@@ -26,7 +26,6 @@
 
 import Data.Array.Unboxed
 import Data.Array.ST
-import Data.Array.Base (unsafeAt, unsafeWrite)
 
 import Data.Bits
 #if __GLASGOW_HASKELL__ < 705
@@ -34,6 +33,7 @@
 #endif
 
 import Math.NumberTheory.Logarithms.Internal (integerLog2#)
+import Math.NumberTheory.Unsafe
 #if __GLASGOW_HASKELL__ < 707
 import Math.NumberTheory.Utils (isTrue#)
 #endif
@@ -41,6 +41,10 @@
 -- | Calculate the integer fourth root of a nonnegative number,
 --   that is, the largest integer @r@ with @r^4 <= n@.
 --   Throws an error on negaitve input.
+{-# SPECIALISE integerFourthRoot :: Int -> Int,
+                                    Integer -> Integer,
+                                    Word -> Word
+  #-}
 integerFourthRoot :: Integral a => a -> a
 integerFourthRoot n
     | n < 0     = error "integerFourthRoot: negative argument"
@@ -50,9 +54,9 @@
 --   that is, the largest integer @r@ with @r^4 <= n@.
 --   The condition is /not/ checked.
 {-# RULES
-"integerFourthRoot'/Int"  integerFourthRoot' = biSqrtInt
-"integerFourthRoot'/Word" integerFourthRoot' = biSqrtWord
-"integerFourthRoot'/Igr"  integerFourthRoot' = biSqrtIgr
+"integerFourthRoot'/Int"     integerFourthRoot' = biSqrtInt
+"integerFourthRoot'/Word"    integerFourthRoot' = biSqrtWord
+"integerFourthRoot'/Integer" integerFourthRoot' = biSqrtIgr
   #-}
 {-# INLINE [1] integerFourthRoot' #-}
 integerFourthRoot' :: Integral a => a -> a
diff --git a/Math/NumberTheory/Powers/General.hs b/Math/NumberTheory/Powers/General.hs
--- a/Math/NumberTheory/Powers/General.hs
+++ b/Math/NumberTheory/Powers/General.hs
@@ -162,14 +162,17 @@
 --   of the prime exponents if some have been found, otherwise by trying
 --   prime exponents recursively.
 highestPower :: Integral a => a -> (a, Int)
-highestPower n
-  | abs n <= 1  = (n,3)
-  | n < 0       = case integerHighPower (toInteger $ negate n) of
+highestPower n'
+  | abs n <= 1  = (n', 3)
+  | n < 0       = case integerHighPower (negate n) of
                     (r,e) -> case shiftToOddCount e of
                                (k, o) -> (negate $ fromInteger (sqr k r), o)
-  | otherwise   = case integerHighPower (toInteger n) of
+  | otherwise   = case integerHighPower n of
                     (r,e) -> (fromInteger r, e)
     where
+      n :: Integer
+      n = toInteger n'
+
       sqr :: Int -> Integer -> Integer
       sqr 0 m = m
       sqr k m = sqr (k-1) (m*m)
@@ -199,7 +202,8 @@
 newtonK :: Integral a => a -> a -> a -> a
 newtonK k n a = go (step a)
   where
-    step m = ((k-1)*m + n `quot` (m^(k-1))) `quot` k
+    -- Beware integer overflow in m^(k-1)
+    step m = ((k-1)*m + fromInteger (toInteger n `quot` (toInteger m^(k-1)))) `quot` k
     go m
       | l < m     = go l
       | otherwise = m
@@ -321,13 +325,14 @@
                           | otherwise -> go (Set.fromDistinctAscList $ takeWhile (<= mx) $ take (k+1) (iterate (*2) 1)) o iops
   where
     mset k st = fst (Set.split (mx+1) (Set.mapMonotonic (*k) st))
+    -- unP p m = (k, m / p ^ k), where k is as large as possible such that p ^ k still divides m
     unP :: Int -> Int -> (Int,Int)
     unP p m = goP 0 m
       where
         goP :: Int -> Int -> (Int,Int)
-        goP !i j = case m `quotRem` p of
+        goP !i j = case j `quotRem` p of
                      (q,r) | r == 0 -> goP (i+1) q
-                           | otherwise -> (0,j)
+                           | otherwise -> (i,j)
     iops :: [Int]
     iops = 3:5:prs
     prs :: [Int]
@@ -343,5 +348,6 @@
       | otherwise =
         case unP p m of
           (0,_) -> go st m ps
-          (k,r) -> go (Set.unions (st:take k (iterate (mset p) st))) r ps
+          -- iterate f x = [x, f x, f (f x)...]
+          (k,r) -> go (Set.unions (take (k + 1) (iterate (mset p) st))) r ps
     go st m [] = go st m [m+1]
diff --git a/Math/NumberTheory/Powers/Squares.hs b/Math/NumberTheory/Powers/Squares.hs
--- a/Math/NumberTheory/Powers/Squares.hs
+++ b/Math/NumberTheory/Powers/Squares.hs
@@ -29,7 +29,6 @@
 
 import Data.Array.Unboxed
 import Data.Array.ST
-import Data.Array.Base (unsafeAt, unsafeWrite)
 
 import Data.Bits
 #if __GLASGOW_HASKELL__ < 705
@@ -37,6 +36,7 @@
 #endif
 
 import Math.NumberTheory.Logarithms.Internal (integerLog2#)
+import Math.NumberTheory.Unsafe
 #if __GLASGOW_HASKELL__ < 707
 import Math.NumberTheory.Utils (isTrue#)
 #endif
diff --git a/Math/NumberTheory/Primes/Counting.hs b/Math/NumberTheory/Primes/Counting.hs
--- a/Math/NumberTheory/Primes/Counting.hs
+++ b/Math/NumberTheory/Primes/Counting.hs
@@ -11,10 +11,14 @@
 module Math.NumberTheory.Primes.Counting
     ( -- * Exact functions
       primeCount
+    , primeCountMaxArg
     , nthPrime
+    , nthPrimeMaxArg
       -- * Approximations
     , approxPrimeCount
+    , approxPrimeCountOverestimateLimit
     , nthPrimeApprox
+    , nthPrimeApproxUnderestimateLimit
     ) where
 
 import Math.NumberTheory.Primes.Counting.Impl
diff --git a/Math/NumberTheory/Primes/Counting/Approximate.hs b/Math/NumberTheory/Primes/Counting/Approximate.hs
--- a/Math/NumberTheory/Primes/Counting/Approximate.hs
+++ b/Math/NumberTheory/Primes/Counting/Approximate.hs
@@ -12,24 +12,37 @@
 {-# OPTIONS_HADDOCK hide #-}
 module Math.NumberTheory.Primes.Counting.Approximate
     ( approxPrimeCount
+    , approxPrimeCountOverestimateLimit
     , nthPrimeApprox
+    , nthPrimeApproxUnderestimateLimit
     ) where
 
--- | @'approxPrimeCount' n@ gives (for @n > 0@) an
+-- For prime p = 3742914359 we have
+--   approxPrimeCount p = 178317879
+--         primeCount p = 178317880
+
+-- | Following property holds:
+--
+-- > approxPrimeCount n >= primeCount n || n >= approxPrimeCountOverestimateLimit
+approxPrimeCountOverestimateLimit :: Integral a => a
+approxPrimeCountOverestimateLimit = 3742914359
+
+-- | @'approxPrimeCount' n@ gives an
 --   approximation of the number of primes not exceeding
 --   @n@. The approximation is fairly good for @n@ large enough.
---   The number of primes should be slightly overestimated
---   (so it is suitable for allocation of storage) and is
---   never underestimated for @n <= 10^12@.
 approxPrimeCount :: Integral a => a -> a
-approxPrimeCount = truncate . appi . fromIntegral
+approxPrimeCount = truncate . max 0 . appi . fromIntegral
 
--- | @'nthPrimeApprox' n@ gives (for @n > 0@) an
+-- | Following property holds:
+--
+-- > nthPrimeApprox n <= nthPrime n || n >= nthPrimeApproxUnderestimateLimit
+nthPrimeApproxUnderestimateLimit :: Integer
+nthPrimeApproxUnderestimateLimit = 1000000000000
+
+-- | @'nthPrimeApprox' n@ gives an
 --   approximation to the n-th prime. The approximation
---   is fairly good for @n@ large enough. Dual to
---   @'approxPrimeCount'@, this estimate should err
---   on the low side (and does for @n < 10^12@).
-nthPrimeApprox :: Integral a => a -> a
+--   is fairly good for @n@ large enough.
+nthPrimeApprox :: Integer -> Integer
 nthPrimeApprox = max 1 . truncate . nthApp . fromIntegral . max 3
 
 -- Basically the approximation of the prime count by Li(x),
diff --git a/Math/NumberTheory/Primes/Counting/Impl.hs b/Math/NumberTheory/Primes/Counting/Impl.hs
--- a/Math/NumberTheory/Primes/Counting/Impl.hs
+++ b/Math/NumberTheory/Primes/Counting/Impl.hs
@@ -13,7 +13,12 @@
 {-# OPTIONS_GHC -fspec-constr-count=24 #-}
 #endif
 {-# OPTIONS_HADDOCK hide #-}
-module Math.NumberTheory.Primes.Counting.Impl (primeCount, nthPrime) where
+module Math.NumberTheory.Primes.Counting.Impl
+    ( primeCount
+    , primeCountMaxArg
+    , nthPrime
+    , nthPrimeMaxArg
+    ) where
 
 #include "MachDeps.h"
 
@@ -23,8 +28,8 @@
 import Math.NumberTheory.Powers.Squares
 import Math.NumberTheory.Powers.Cubes
 import Math.NumberTheory.Logarithms
+import Math.NumberTheory.Unsafe
 
-import Data.Array.Base
 import Data.Array.ST
 #if !MIN_VERSION_array(0,5,0)
     hiding (unsafeThaw)
@@ -39,10 +44,14 @@
 #define COUNT_T Int
 #endif
 
+-- | Maximal allowed argument of 'primeCount'. Currently 8e18.
+primeCountMaxArg :: Integer
+primeCountMaxArg = 8000000000000000000
+
 -- | @'primeCount' n == &#960;(n)@ is the number of (positive) primes not exceeding @n@.
 --
 --   For efficiency, the calculations are done on 64-bit signed integers, therefore @n@ must
---   not exceed @8 * 10^18@.
+--   not exceed 'primeCountMaxArg'.
 --
 --   Requires @/O/(n^0.5)@ space, the time complexity is roughly @/O/(n^0.7)@.
 --   For small bounds, @'primeCount' n@ simply counts the primes not exceeding @n@,
@@ -51,7 +60,7 @@
 --   <http://en.wikipedia.org/wiki/Prime_counting_function#Algorithms_for_evaluating_.CF.80.28x.29>.
 primeCount :: Integer -> Integer
 primeCount n
-    | n > 8000000000000000000   = error $ "primeCount: can't handle bound " ++ show n
+    | n > primeCountMaxArg = error $ "primeCount: can't handle bound " ++ show n
     | n < 2     = 0
     | n < 1000  = fromIntegral . length . takeWhile (<= n) . primeList . primeSieve $ max 242 n
     | n < 30000 = runST $ do
@@ -69,15 +78,19 @@
             !pdf = sieveCount ub cs sr
         in phn1 - pdf
 
+-- | Maximal allowed argument of 'nthPrime'. Currently 1.5e17.
+nthPrimeMaxArg :: Integer
+nthPrimeMaxArg = 150000000000000000
+
 -- | @'nthPrime' n@ calculates the @n@-th prime. Numbering of primes is
 --   @1@-based, so @'nthPrime' 1 == 2@.
 --
 --   Requires @/O/((n*log n)^0.5)@ space, the time complexity is roughly @/O/((n*log n)^0.7@.
---   The argument must be strictly positive, and must not exceed @1.5 * 10^17@.
+--   The argument must be strictly positive, and must not exceed 'nthPrimeMaxArg'.
 nthPrime :: Integer -> Integer
 nthPrime n
     | n < 1         = error "Prime indexing starts at 1"
-    | n > 150000000000000000    = error $ "nthPrime: can't handle index " ++ show n
+    | n > nthPrimeMaxArg = error $ "nthPrime: can't handle index " ++ show n
     | n < 200000    = nthPrimeCt n
     | ct0 < n       = tooLow n p0 (n-ct0) approxGap
     | otherwise     = tooHigh n p0 (ct0-n) approxGap
diff --git a/Math/NumberTheory/Primes/Factorisation.hs b/Math/NumberTheory/Primes/Factorisation.hs
--- a/Math/NumberTheory/Primes/Factorisation.hs
+++ b/Math/NumberTheory/Primes/Factorisation.hs
@@ -48,6 +48,10 @@
     , carmichaelSieve
     , sieveCarmichael
     , carmichaelFromCanonical
+      -- * Moebius function
+    , moebius
+    , μ
+    , moebiusFromCanonical
       -- * Divisors
     , divisors
     , tau
@@ -101,7 +105,7 @@
 φ = totient
 
 -- | Calculates the Carmichael function for a positive integer, that is,
---   the (smallest) exponent of the group of units in @&8484;/(n)@.
+--   the (smallest) exponent of the group of units in @&#8484;/(n)@.
 carmichael :: Integer -> Integer
 carmichael n
     | n < 1     = error "Carmichael function only defined for positive numbers"
@@ -111,6 +115,17 @@
 -- | Alias of 'carmichael' for people who prefer Greek letters.
 λ :: Integer -> Integer
 λ = carmichael
+
+-- | Calculates the Moebius function for a positive integer.
+moebius :: Integer -> Integer
+moebius n
+    | n < 1     = error "Carmichael function only defined for positive numbers"
+    | n == 1    = 1
+    | otherwise = moebiusFromCanonical (factorise' n)
+
+-- | Alias of 'moebius' for people who prefer Greek letters.
+μ :: Integer -> Integer
+μ = moebius
 
 -- | @'divisors' n@ is the set of all (positive) divisors of @n@.
 --   @'divisors' 0@ is an error because we can't create the set of all 'Integer's.
diff --git a/Math/NumberTheory/Primes/Factorisation/Montgomery.hs b/Math/NumberTheory/Primes/Factorisation/Montgomery.hs
--- a/Math/NumberTheory/Primes/Factorisation/Montgomery.hs
+++ b/Math/NumberTheory/Primes/Factorisation/Montgomery.hs
@@ -46,7 +46,6 @@
 #if __GLASGOW_HASKELL__ < 705
 import GHC.Word     -- Moved to GHC.Types
 #endif
-import Data.Array.Base
 
 import System.Random
 import Control.Monad.State.Strict
@@ -63,6 +62,7 @@
 import Math.NumberTheory.Primes.Sieve.Eratosthenes
 import Math.NumberTheory.Primes.Sieve.Indexing
 import Math.NumberTheory.Primes.Testing.Probabilistic
+import Math.NumberTheory.Unsafe
 import Math.NumberTheory.Utils
 
 -- | @'factorise' n@ produces the prime factorisation of @n@, including
diff --git a/Math/NumberTheory/Primes/Factorisation/Utils.hs b/Math/NumberTheory/Primes/Factorisation/Utils.hs
--- a/Math/NumberTheory/Primes/Factorisation/Utils.hs
+++ b/Math/NumberTheory/Primes/Factorisation/Utils.hs
@@ -13,6 +13,7 @@
     ( ppTotient
     , totientFromCanonical
     , carmichaelFromCanonical
+    , moebiusFromCanonical
     , divisorsFromCanonical
     , tauFromCanonical
     , divisorSumFromCanonical
@@ -49,6 +50,15 @@
     go !acc ((p,1):pps) = go (lcm acc (p-1)) pps
     go acc ((p,k):pps)  = go ((lcm acc (p-1))*integerPower p (k-1)) pps
     go acc []           = acc
+
+-- | Calculate the Moebius function from the canonical factorisation.
+moebiusFromCanonical :: [(a, Int)] -> Integer
+moebiusFromCanonical = go 1
+  where
+  go acc []            = acc
+  go acc ((_, 1) : xs) = go (negate acc) xs
+  go acc ((_, 0) : xs) = go acc xs          -- Should not really happen
+  go _   _             = 0                  -- Short circuit for powers > 1
 
 -- | The set of divisors, efficiently calculated from the canonical factorisation.
 divisorsFromCanonical :: [(Integer,Int)] -> Set Integer
diff --git a/Math/NumberTheory/Primes/Heap.hs b/Math/NumberTheory/Primes/Heap.hs
--- a/Math/NumberTheory/Primes/Heap.hs
+++ b/Math/NumberTheory/Primes/Heap.hs
@@ -22,7 +22,6 @@
 module Math.NumberTheory.Primes.Heap (primes, sieveFrom) where
 
 import Data.Array.Unboxed
-import Data.Array.Base (unsafeAt, unsafeRead, unsafeWrite)
 import Data.Array.ST
 import Control.Monad.ST
 import Data.List (foldl')
@@ -30,6 +29,8 @@
 import Data.Word
 #endif
 
+import Math.NumberTheory.Unsafe
+
 #ifndef SH_SIZE
 #define SH_SIZE 31
 #endif
@@ -331,12 +332,14 @@
     | r < m     = down (a-1)
     | otherwise = up a
       where
-        a = min 5758 ((192*r) `quot` 1001 - 1)
+        a = max 0 (min 5758 ((192*r) `quot` 1001 - 1))
         m = remainders `unsafeAt` a
         down k
+            | k < 0                         = 0
             | r < (remainders `unsafeAt` k) = down (k-1)
             | otherwise                     = k
         up k
+            | k+1 > 5759                        = 5759
             | r < (remainders `unsafeAt` (k+1)) = k
             | otherwise                         = up (k+1)
 
diff --git a/Math/NumberTheory/Primes/Sieve/Eratosthenes.hs b/Math/NumberTheory/Primes/Sieve/Eratosthenes.hs
--- a/Math/NumberTheory/Primes/Sieve/Eratosthenes.hs
+++ b/Math/NumberTheory/Primes/Sieve/Eratosthenes.hs
@@ -35,18 +35,18 @@
 #include "MachDeps.h"
 
 import Control.Monad.ST
-import Data.Array.Base
 import Data.Array.ST
 #if !MIN_VERSION_array(0,5,0)
                      hiding (unsafeFreeze, unsafeThaw, castSTUArray)
 #endif
 import Control.Monad (when)
 import Data.Bits
-#if __GLASGOW_HASKELL__ < 709
+#if __GLASGOW_HASKELL__ < 709 || WORD_SIZE_IN_BITS == 32
 import Data.Word
 #endif
 
 import Math.NumberTheory.Powers.Squares (integerSquareRoot)
+import Math.NumberTheory.Unsafe
 import Math.NumberTheory.Utils
 import Math.NumberTheory.Primes.Counting.Approximate
 import Math.NumberTheory.Primes.Sieve.Indexing
@@ -97,7 +97,7 @@
 -- | Compact store of primality flags.
 data PrimeSieve = PS !Integer {-# UNPACK #-} !(UArray Int Bool)
 
--- | Sieve primes up to (and including) a bound.
+-- | Sieve primes up to (and including) a bound (or 7, if bound is smaller).
 --   For small enough bounds, this is more efficient than
 --   using the segmented sieve.
 --
diff --git a/Math/NumberTheory/Primes/Sieve/Indexing.hs b/Math/NumberTheory/Primes/Sieve/Indexing.hs
--- a/Math/NumberTheory/Primes/Sieve/Indexing.hs
+++ b/Math/NumberTheory/Primes/Sieve/Indexing.hs
@@ -21,9 +21,10 @@
     ) where
 
 import Data.Array.Unboxed
-import Data.Array.Base (unsafeAt)
 import Data.Bits
 
+import Math.NumberTheory.Unsafe
+
 -- Auxiliary stuff, conversion between number and index,
 -- remainders modulo 30 and related things.
 
@@ -33,7 +34,9 @@
 --   #-}
 {-# INLINE idxPr #-}
 idxPr :: Integral a => a -> (Int,Int)
-idxPr n0 = (fromIntegral bytes0, rm3)
+idxPr n0
+    | n0 < 7    = (0, 0)
+    | otherwise = (fromIntegral bytes0, rm3)
   where
     n = if (fromIntegral n0 .&. 1 == (1 :: Int))
             then n0 else (n0-1)
diff --git a/Math/NumberTheory/Primes/Sieve/Misc.hs b/Math/NumberTheory/Primes/Sieve/Misc.hs
--- a/Math/NumberTheory/Primes/Sieve/Misc.hs
+++ b/Math/NumberTheory/Primes/Sieve/Misc.hs
@@ -31,7 +31,6 @@
     ) where
 
 import Control.Monad.ST
-import Data.Array.Base (unsafeRead, unsafeWrite, unsafeAt)
 import Data.Array.ST
 import Data.Array.Unboxed
 import Control.Monad (when)
@@ -44,6 +43,7 @@
 import Math.NumberTheory.Primes.Sieve.Indexing
 import Math.NumberTheory.Primes.Factorisation.Montgomery
 import Math.NumberTheory.Primes.Factorisation.Utils
+import Math.NumberTheory.Unsafe
 import Math.NumberTheory.Utils
 
 -- | A compact store of smallest prime factors.
@@ -149,7 +149,7 @@
                                                  | otherwise -> tdLoop j (integerSquareRoot' j) (ix+1)
           where
             p = toPrim ix
-            pix = unsafeAt sve ix
+            pix = unsafeAt sve $ fromIntegral p
     curve n = stdGenFactorisation (Just (bound*(bound+2))) (mkStdGen $ fromIntegral n `xor` 0xdecaf00d) Nothing n
 
 -- | @'totientSieve' n@ creates a store of the totients of the numbers not exceeding @n@.
diff --git a/Math/NumberTheory/Unsafe.hs b/Math/NumberTheory/Unsafe.hs
new file mode 100644
--- /dev/null
+++ b/Math/NumberTheory/Unsafe.hs
@@ -0,0 +1,73 @@
+-- |
+-- Module:      Math.NumberTheory.Unsafe
+-- Copyright:   (c) 2016 Andrew Lelechenko
+-- Licence:     MIT
+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>
+-- Stability:   Provisional
+--
+-- Layer to switch between safe and unsafe arrays.
+--
+
+{-# LANGUAGE CPP #-}
+
+module Math.NumberTheory.Unsafe
+  ( UArray
+  , bounds
+  , castSTUArray
+  , unsafeAt
+  , unsafeFreeze
+  , unsafeNewArray_
+  , unsafeRead
+  , unsafeThaw
+  , unsafeWrite
+  ) where
+
+#ifdef CheckBounds
+
+import Data.Array.Base
+  ( UArray
+  , castSTUArray
+  )
+import Data.Array.IArray
+  ( IArray
+  , bounds
+  , (!)
+  )
+import Data.Array.MArray
+#if !MIN_VERSION_array(0,5,0)
+  hiding (unsafeFreeze, unsafeThaw)
+#endif
+
+unsafeAt :: (IArray a e, Ix i) => a i e -> i -> e
+unsafeAt = (!)
+
+unsafeFreeze :: (Ix i, MArray a e m, IArray b e) => a i e -> m (b i e)
+unsafeFreeze = freeze
+
+unsafeNewArray_ :: (Ix i, MArray a e m) => (i, i) -> m (a i e)
+unsafeNewArray_ = newArray_
+
+unsafeRead :: (MArray a e m, Ix i) => a i e -> i -> m e
+unsafeRead = readArray
+
+unsafeThaw :: (Ix i, IArray a e, MArray b e m) => a i e -> m (b i e)
+unsafeThaw = thaw
+
+unsafeWrite :: (MArray a e m, Ix i) => a i e -> i -> e -> m ()
+unsafeWrite = writeArray
+
+#else
+
+import Data.Array.Base
+  ( UArray
+  , bounds
+  , castSTUArray
+  , unsafeAt
+  , unsafeFreeze
+  , unsafeNewArray_
+  , unsafeRead
+  , unsafeThaw
+  , unsafeWrite
+  )
+
+#endif
diff --git a/Math/NumberTheory/Utils.hs b/Math/NumberTheory/Utils.hs
--- a/Math/NumberTheory/Utils.hs
+++ b/Math/NumberTheory/Utils.hs
@@ -63,7 +63,7 @@
 "shiftToOddCount/Integer"   shiftToOddCount = shiftOCInteger
   #-}
 {-# INLINE [1] shiftToOddCount #-}
-shiftToOddCount :: (Integral a, Bits a) => a -> (Int, a)
+shiftToOddCount :: Integral a => a -> (Int, a)
 shiftToOddCount n = case shiftOCInteger (fromIntegral n) of
                       (z, o) -> (z, fromInteger o)
 
@@ -125,7 +125,7 @@
 "shiftToOdd/Integer"   shiftToOdd = shiftOInteger
   #-}
 {-# INLINE [1] shiftToOdd #-}
-shiftToOdd :: (Integral a, Bits a) => a -> a
+shiftToOdd :: Integral a => a -> a
 shiftToOdd n = fromInteger (shiftOInteger (fromIntegral n))
 
 -- | Specialised version for @'Int'@.
diff --git a/arithmoi.cabal b/arithmoi.cabal
--- a/arithmoi.cabal
+++ b/arithmoi.cabal
@@ -1,5 +1,5 @@
 name                : arithmoi
-version             : 0.4.1.3
+version             : 0.4.2.0
 cabal-version       : >= 1.10
 author              : Daniel Fischer
 copyright           : (c) 2011 Daniel Fischer
@@ -27,12 +27,12 @@
 
 category            : Math, Algorithms, Number Theory
 
-tested-with         : GHC == 7.4.2, GHC==7.6.3, GHC==7.8.3
+tested-with         : GHC==7.6.3, GHC==7.8.4, GHC==7.10.3, GHC==8.0.1
 
 extra-source-files  : Changes, TODO
 
-flag llvm
-    description         : Compile the library with the LLVM backend
+flag check-bounds
+    description         : Replace unsafe array operations with safe ones
     default             : False
     manual              : True
 
@@ -40,7 +40,7 @@
     default-language: Haskell2010
     build-depends       : base >= 4 && < 5
                           , array >= 0.3 && < 0.6
-                          , ghc-prim < 0.5
+                          , ghc-prim < 0.6
                           , integer-gmp < 1.1
                           , containers >= 0.3 && < 0.6
                           , random >= 1.0 && < 1.2
@@ -51,6 +51,7 @@
                           Math.NumberTheory.MoebiusInversion
                           Math.NumberTheory.MoebiusInversion.Int
                           Math.NumberTheory.Lucas
+                          Math.NumberTheory.GaussianIntegers
                           Math.NumberTheory.GCD
                           Math.NumberTheory.GCD.LowLevel
                           Math.NumberTheory.Powers
@@ -68,6 +69,7 @@
                           Math.NumberTheory.Primes.Testing.Certificates
                           Math.NumberTheory.Primes.Heap
     other-modules       : Math.NumberTheory.Utils
+                          Math.NumberTheory.Unsafe
                           Math.NumberTheory.Logarithms.Internal
                           Math.NumberTheory.Primes.Counting.Impl
                           Math.NumberTheory.Primes.Counting.Approximate
@@ -83,9 +85,9 @@
     other-extensions          : BangPatterns
 
     ghc-options         : -O2 -Wall
-    if flag(llvm)
-        ghc-options     : -fllvm
     ghc-prof-options    : -O2 -auto
+    if flag(check-bounds)
+        cpp-options     : -DCheckBounds
 
 source-repository head
   type:     git
@@ -105,9 +107,31 @@
   type:                 exitcode-stdio-1.0
   hs-source-dirs:       test-suite
   ghc-options:          -Wall
-  main-is:              Spec.hs
+  main-is:              Test.hs
   default-language: Haskell2010
-  build-depends:        base
-                        ,hspec
-                        ,arithmoi
-
+  build-depends:        base >= 4 && < 5
+                        , arithmoi >= 0.4 && < 0.5
+                        , tasty >= 0.10 && < 0.12
+                        , tasty-smallcheck >= 0.8 && < 0.9
+                        , tasty-quickcheck >= 0.8 && < 0.9
+                        , tasty-hunit >= 0.9 && < 0.10
+                        , QuickCheck >= 2.8 && < 2.9
+                        , smallcheck >= 1.1 && < 1.2
+  other-modules :   Math.NumberTheory.GaussianIntegersTests
+                  , Math.NumberTheory.GCDTests
+                  , Math.NumberTheory.GCD.LowLevelTests
+                  , Math.NumberTheory.LogarithmsTests
+                  , Math.NumberTheory.LucasTests
+                  , Math.NumberTheory.ModuliTests
+                  , Math.NumberTheory.Powers.CubesTests
+                  , Math.NumberTheory.MoebiusInversionTests
+                  , Math.NumberTheory.MoebiusInversion.IntTests
+                  , Math.NumberTheory.Powers.FourthTests
+                  , Math.NumberTheory.Powers.GeneralTests
+                  , Math.NumberTheory.Powers.IntegerTests
+                  , Math.NumberTheory.Powers.SquaresTests
+                  , Math.NumberTheory.PrimesTests
+                  , Math.NumberTheory.Primes.CountingTests
+                  , Math.NumberTheory.Primes.HeapTests
+                  , Math.NumberTheory.Primes.SieveTests
+                  , Math.NumberTheory.TestUtils
diff --git a/test-suite/Math/NumberTheory/GCD/LowLevelTests.hs b/test-suite/Math/NumberTheory/GCD/LowLevelTests.hs
new file mode 100644
--- /dev/null
+++ b/test-suite/Math/NumberTheory/GCD/LowLevelTests.hs
@@ -0,0 +1,74 @@
+-- |
+-- Module:      Math.NumberTheory.GCD.LowLevelTests
+-- Copyright:   (c) 2016 Andrew Lelechenko
+-- Licence:     MIT
+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>
+-- Stability:   Provisional
+--
+-- Tests for Math.NumberTheory.GCD.LowLevel
+--
+
+{-# LANGUAGE CPP       #-}
+{-# LANGUAGE MagicHash #-}
+
+{-# OPTIONS_GHC -fno-warn-type-defaults #-}
+
+module Math.NumberTheory.GCD.LowLevelTests
+  ( testSuite
+  ) where
+
+import Test.Tasty
+
+#if MIN_VERSION_base(4,8,0)
+#else
+import Data.Word
+#endif
+
+import GHC.Exts
+
+import Math.NumberTheory.GCD.LowLevel
+import Math.NumberTheory.TestUtils
+
+-- | Check that 'gcdInt' matches 'gcd'.
+gcdIntProperty :: Int -> Int -> Bool
+gcdIntProperty a b = gcdInt a b == gcd a b
+
+-- | Check that 'gcdWord' matches 'gcd'.
+gcdWordProperty :: Word -> Word -> Bool
+gcdWordProperty a b = gcdWord a b == gcd a b
+
+-- | Check that 'gcdInt#' matches 'gcd'.
+gcdIntProperty# :: Int -> Int -> Bool
+gcdIntProperty# a@(I# a') b@(I# b') = I# (gcdInt# a' b') == gcd a b
+
+-- | Check that 'gcdWord#' matches 'gcd'.
+gcdWordProperty# :: Word -> Word -> Bool
+gcdWordProperty# a@(W# a') b@(W# b') = W# (gcdWord# a' b') == gcd a b
+
+-- | Check that numbers are coprime iff their gcd equals to 1.
+coprimeIntProperty :: Int -> Int -> Bool
+coprimeIntProperty a b = coprimeInt a b == (gcd a b == 1)
+
+-- | Check that numbers are coprime iff their gcd equals to 1.
+coprimeWordProperty :: Word -> Word -> Bool
+coprimeWordProperty a b = coprimeWord a b == (gcd a b == 1)
+
+-- | Check that numbers are coprime iff their gcd equals to 1.
+coprimeIntProperty# :: Int -> Int -> Bool
+coprimeIntProperty# a@(I# a') b@(I# b') = coprimeInt# a' b' == (gcd a b == 1)
+
+-- | Check that numbers are coprime iff their gcd equals to 1.
+coprimeWordProperty# :: Word -> Word -> Bool
+coprimeWordProperty# a@(W# a') b@(W# b') = coprimeWord# a' b' == (gcd a b == 1)
+
+testSuite :: TestTree
+testSuite = testGroup "LowLevel"
+  [ testSmallAndQuick "gcdInt"       gcdIntProperty
+  , testSmallAndQuick "gcdWord"      gcdWordProperty
+  , testSmallAndQuick "gcdInt#"      gcdIntProperty#
+  , testSmallAndQuick "gcdWord#"     gcdWordProperty#
+  , testSmallAndQuick "coprimeInt"   coprimeIntProperty
+  , testSmallAndQuick "coprimeWord"  coprimeWordProperty
+  , testSmallAndQuick "coprimeInt#"  coprimeIntProperty#
+  , testSmallAndQuick "coprimeWord#" coprimeWordProperty#
+  ]
diff --git a/test-suite/Math/NumberTheory/GCDTests.hs b/test-suite/Math/NumberTheory/GCDTests.hs
new file mode 100644
--- /dev/null
+++ b/test-suite/Math/NumberTheory/GCDTests.hs
@@ -0,0 +1,50 @@
+-- |
+-- Module:      Math.NumberTheory.GCDTests
+-- Copyright:   (c) 2016 Andrew Lelechenko
+-- Licence:     MIT
+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>
+-- Stability:   Provisional
+--
+-- Tests for Math.NumberTheory.GCD
+--
+
+{-# LANGUAGE ScopedTypeVariables #-}
+
+{-# OPTIONS_GHC -fno-warn-type-defaults #-}
+
+module Math.NumberTheory.GCDTests
+  ( testSuite
+  ) where
+
+import Test.Tasty
+
+import Data.Bits
+
+import Math.NumberTheory.GCD
+import Math.NumberTheory.TestUtils
+
+-- | Check that 'binaryGCD' matches 'gcd'.
+binaryGCDProperty :: (Integral a, Bits a) => AnySign a -> AnySign a -> Bool
+binaryGCDProperty (AnySign a) (AnySign b) = binaryGCD a b == gcd a b
+
+-- | Check that 'extendedGCD' is consistent with documentation.
+extendedGCDProperty :: forall a. Integral a => AnySign a -> AnySign a -> Bool
+extendedGCDProperty (AnySign a) (AnySign b) =
+  u * a + v * b == d
+  && d == gcd a b
+  -- (-1) >= 0 is true for unsigned types
+  && (abs u < abs b || abs b <= 1 || (-1 :: a) >= 0)
+  && (abs v < abs a || abs a <= 1 || (-1 :: a) >= 0)
+  where
+    (d, u, v) = extendedGCD a b
+
+-- | Check that numbers are coprime iff their gcd equals to 1.
+coprimeProperty :: (Integral a, Bits a) => AnySign a -> AnySign a -> Bool
+coprimeProperty (AnySign a) (AnySign b) = coprime a b == (gcd a b == 1)
+
+testSuite :: TestTree
+testSuite = testGroup "GCD"
+  [ testSameIntegralProperty "binaryGCD"   binaryGCDProperty
+  , testSameIntegralProperty "extendedGCD" extendedGCDProperty
+  , testSameIntegralProperty "coprime"     coprimeProperty
+  ]
diff --git a/test-suite/Math/NumberTheory/GaussianIntegersTests.hs b/test-suite/Math/NumberTheory/GaussianIntegersTests.hs
new file mode 100644
--- /dev/null
+++ b/test-suite/Math/NumberTheory/GaussianIntegersTests.hs
@@ -0,0 +1,79 @@
+{-# OPTIONS_GHC -fno-warn-type-defaults #-}
+
+-- |
+-- Module:      Math.NumberTheory.GaussianIntegersTests
+-- Copyright:   (c) 2016 Chris Fredrickson
+-- Licence:     MIT
+-- Maintainer:  Chris Fredrickson <chris.p.fredrickson@gmail.com>
+-- Stability:   Provisional
+--
+-- Tests for Math.NumberTheory.GaussianIntegers
+--
+
+module Math.NumberTheory.GaussianIntegersTests
+  ( testSuite
+  ) where
+
+import Test.Tasty
+import Test.Tasty.HUnit
+
+import Math.NumberTheory.GaussianIntegers
+import Math.NumberTheory.TestUtils
+
+-- | Number is zero or is equal to the product of its factors.
+factoriseProperty :: Integer -> Integer -> Bool
+factoriseProperty x y
+  =  x == 0 && y == 0
+  || g == g'
+  where
+    g = x :+ y
+    factors = factorise g
+    g' = product $ map (uncurry (.^)) factors
+
+-- | Number is prime iff it is non-zero
+--   and has exactly one (non-unit) factor.
+isPrimeProperty :: Integer -> Integer -> Bool
+isPrimeProperty x y
+  =  x == 0 && y == 0
+  || isPrime g && n == 1
+  || not (isPrime g) && n /= 1
+  where
+    g = x :+ y
+    factors = factorise g
+    nonUnitFactors = filter (\(p, _) -> norm p /= 1) factors
+    -- Count factors taking into account multiplicity
+    n = sum $ map snd nonUnitFactors
+
+-- | The list of primes should include only primes.
+primesGeneratesPrimesProperty :: NonNegative Int -> Bool
+primesGeneratesPrimesProperty (NonNegative i) = isPrime (primes !! i)
+
+-- | signum and abs should satisfy: z == signum z * abs z
+signumAbsProperty :: Integer -> Integer -> Bool
+signumAbsProperty x y = z == signum z * abs z
+  where
+    z = x :+ y
+
+-- | abs maps a Gaussian integer to its associate in first quadrant.
+absProperty :: Integer -> Integer -> Bool
+absProperty x y = isOrigin || (inFirstQuadrant && isAssociate)
+  where
+    z = x :+ y
+    z'@(x' :+ y') = abs z
+    isOrigin = z' == 0 && z == 0
+    inFirstQuadrant = x' > 0 && y' >= 0     -- first quadrant includes the positive real axis, but not the origin or the positive imaginary axis
+    isAssociate = z' `elem` map (\e -> z * (0 :+ 1) .^ e) [0 .. 3]
+
+-- | a special case that tests rounding/truncating in GCD.
+gcdGSpecialCase1 :: Assertion
+gcdGSpecialCase1 = assertEqual "gcdG" 1 $ gcdG (12 :+ 23) (23 :+ 34)
+
+testSuite :: TestTree
+testSuite = testGroup "GaussianIntegers"
+  [ testSmallAndQuick "factorise"         factoriseProperty
+  , testSmallAndQuick "isPrime"           isPrimeProperty
+  , testSmallAndQuick "primes"            primesGeneratesPrimesProperty
+  , testSmallAndQuick "signumAbsProperty" signumAbsProperty
+  , testSmallAndQuick "absProperty"       absProperty
+  , testCase          "gcdG (12 :+ 23) (23 :+ 34)" gcdGSpecialCase1
+  ]
diff --git a/test-suite/Math/NumberTheory/LogarithmsTests.hs b/test-suite/Math/NumberTheory/LogarithmsTests.hs
new file mode 100644
--- /dev/null
+++ b/test-suite/Math/NumberTheory/LogarithmsTests.hs
@@ -0,0 +1,112 @@
+-- |
+-- Module:      Math.NumberTheory.LogarithmsTests
+-- Copyright:   (c) 2016 Andrew Lelechenko
+-- Licence:     MIT
+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>
+-- Stability:   Provisional
+--
+-- Tests for Math.NumberTheory.Logarithms
+--
+
+{-# LANGUAGE CPP       #-}
+
+{-# OPTIONS_GHC -fno-warn-type-defaults #-}
+
+module Math.NumberTheory.LogarithmsTests
+  ( testSuite
+  ) where
+
+import Test.Tasty
+
+#if MIN_VERSION_base(4,8,0)
+#else
+import Data.Word
+#endif
+
+import Math.NumberTheory.Logarithms
+import Math.NumberTheory.TestUtils
+
+-- | Check that 'integerLogBase' returns the largest integer @l@ such that @b@ ^ @l@ <= @n@ and @b@ ^ (@l@+1) > @n@.
+integerLogBaseProperty :: Positive Integer -> Positive Integer -> Bool
+integerLogBaseProperty (Positive b) (Positive n) = b == 1 || b ^ l <= n && b ^ (l + 1) > n
+  where
+    l = toInteger $ integerLogBase b n
+
+-- | Check that 'integerLog2' returns the largest integer @l@ such that 2 ^ @l@ <= @n@ and 2 ^ (@l@+1) > @n@.
+integerLog2Property :: Positive Integer -> Bool
+integerLog2Property (Positive n) = 2 ^ l <= n && 2 ^ (l + 1) > n
+  where
+    l = toInteger $ integerLog2 n
+
+-- | Check that 'integerLog10' returns the largest integer @l@ such that 10 ^ @l@ <= @n@ and 10 ^ (@l@+1) > @n@.
+integerLog10Property :: Positive Integer -> Bool
+integerLog10Property (Positive n) = 10 ^ l <= n && 10 ^ (l + 1) > n
+  where
+    l = toInteger $ integerLog10 n
+
+-- | Check that 'intLog2' returns the largest integer @l@ such that 2 ^ @l@ <= @n@ and 2 ^ (@l@+1) > @n@.
+intLog2Property :: Positive Int -> Bool
+intLog2Property (Positive n) = 2 ^ l <= n && (2 ^ (l + 1) > n || n > maxBound `div` 2)
+  where
+    l = intLog2 n
+
+-- | Check that 'wordLog2' returns the largest integer @l@ such that 2 ^ @l@ <= @n@ and 2 ^ (@l@+1) > @n@.
+wordLog2Property :: Positive Word -> Bool
+wordLog2Property (Positive n) = 2 ^ l <= n && (2 ^ (l + 1) > n || n > maxBound `div` 2)
+  where
+    l = wordLog2 n
+
+-- | Check that 'integerLogBase'' returns the largest integer @l@ such that @b@ ^ @l@ <= @n@ and @b@ ^ (@l@+1) > @n@.
+integerLogBase'Property :: Positive Integer -> Positive Integer -> Bool
+integerLogBase'Property (Positive b) (Positive n) = b == 1 || b ^ l <= n && b ^ (l + 1) > n
+  where
+    l = toInteger $ integerLogBase' b n
+
+-- | Check that 'integerLogBase'' returns the largest integer @l@ such that @b@ ^ @l@ <= @n@ and @b@ ^ (@l@+1) > @n@ for @b@ > 32 and @n@ >= @b@ ^ 2.
+integerLogBase'Property2 :: Positive Integer -> Positive Integer -> Bool
+integerLogBase'Property2 (Positive b') (Positive n') = b ^ l <= n && b ^ (l + 1) > n
+  where
+    b = b' + 32
+    n = n' + b ^ 2 - 1
+    l = toInteger $ integerLogBase' b n
+
+-- | Check that 'integerLog2'' returns the largest integer @l@ such that 2 ^ @l@ <= @n@ and 2 ^ (@l@+1) > @n@.
+integerLog2'Property :: Positive Integer -> Bool
+integerLog2'Property (Positive n) = 2 ^ l <= n && 2 ^ (l + 1) > n
+  where
+    l = toInteger $ integerLog2' n
+
+-- | Check that 'integerLog10'' returns the largest integer @l@ such that 10 ^ @l@ <= @n@ and 10 ^ (@l@+1) > @n@.
+integerLog10'Property :: Positive Integer -> Bool
+integerLog10'Property (Positive n) = 10 ^ l <= n && 10 ^ (l + 1) > n
+  where
+    l = toInteger $ integerLog10' n
+
+-- | Check that 'intLog2'' returns the largest integer @l@ such that 2 ^ @l@ <= @n@ and 2 ^ (@l@+1) > @n@.
+intLog2'Property :: Positive Int -> Bool
+intLog2'Property (Positive n) = 2 ^ l <= n && (2 ^ (l + 1) > n || n > maxBound `div` 2)
+  where
+    l = intLog2' n
+
+-- | Check that 'wordLog2'' returns the largest integer @l@ such that 2 ^ @l@ <= @n@ and 2 ^ (@l@+1) > @n@.
+wordLog2'Property :: Positive Word -> Bool
+wordLog2'Property (Positive n) = 2 ^ l <= n && (2 ^ (l + 1) > n || n > maxBound `div` 2)
+  where
+    l = wordLog2' n
+
+testSuite :: TestTree
+testSuite = testGroup "Logarithms"
+  [ testSmallAndQuick "integerLogBase"  integerLogBaseProperty
+  , testSmallAndQuick "integerLog2"     integerLog2Property
+  , testSmallAndQuick "integerLog10"    integerLog10Property
+  , testSmallAndQuick "intLog2"         intLog2Property
+  , testSmallAndQuick "wordLog2"        wordLog2Property
+
+  , testSmallAndQuick "integerLogBase'" integerLogBase'Property
+  , testSmallAndQuick "integerLogBase' with base > 32 and n >= base ^ 2"
+      integerLogBase'Property2
+  , testSmallAndQuick "integerLog2'"    integerLog2'Property
+  , testSmallAndQuick "integerLog10'"   integerLog10'Property
+  , testSmallAndQuick "intLog2'"        intLog2'Property
+  , testSmallAndQuick "wordLog2'"       wordLog2'Property
+  ]
diff --git a/test-suite/Math/NumberTheory/LucasTests.hs b/test-suite/Math/NumberTheory/LucasTests.hs
new file mode 100644
--- /dev/null
+++ b/test-suite/Math/NumberTheory/LucasTests.hs
@@ -0,0 +1,104 @@
+-- |
+-- Module:      Math.NumberTheory.LucasTests
+-- Copyright:   (c) 2016 Andrew Lelechenko
+-- Licence:     MIT
+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>
+-- Stability:   Provisional
+--
+-- Tests for Math.NumberTheory.Lucas
+--
+
+{-# LANGUAGE CPP       #-}
+
+{-# OPTIONS_GHC -fno-warn-type-defaults #-}
+
+module Math.NumberTheory.LucasTests
+  ( testSuite
+  ) where
+
+import Test.Tasty
+import Test.Tasty.HUnit
+
+import Math.NumberTheory.Lucas
+import Math.NumberTheory.TestUtils
+
+-- | Check that 'fibonacci' matches the definition of Fibonacci sequence.
+fibonacciProperty1 :: AnySign Int -> Bool
+fibonacciProperty1 (AnySign n) = fibonacci n + fibonacci (n + 1) == fibonacci (n +2)
+
+-- | Check that 'fibonacci' for negative indices is correctly defined.
+fibonacciProperty2 :: NonNegative Int -> Bool
+fibonacciProperty2 (NonNegative n) = fibonacci n == (if even n then negate else id) (fibonacci (- n))
+
+-- | Check that 'fibonacciPair' is a pair of consequent 'fibonacci'.
+fibonacciPairProperty :: AnySign Int -> Bool
+fibonacciPairProperty (AnySign n) = fibonacciPair n == (fibonacci n, fibonacci (n + 1))
+
+-- | Check that 'fibonacci 0' is 0.
+fibonacciSpecialCase0 :: Assertion
+fibonacciSpecialCase0 = assertEqual "fibonacci" (fibonacci 0) 0
+
+-- | Check that 'fibonacci 1' is 1.
+fibonacciSpecialCase1 :: Assertion
+fibonacciSpecialCase1 = assertEqual "fibonacci" (fibonacci 1) 1
+
+
+-- | Check that 'lucas' matches the definition of Lucas sequence.
+lucasProperty1 :: AnySign Int -> Bool
+lucasProperty1 (AnySign n) = lucas n + lucas (n + 1) == lucas (n +2)
+
+-- | Check that 'lucas' for negative indices is correctly defined.
+lucasProperty2 :: NonNegative Int -> Bool
+lucasProperty2 (NonNegative n) = lucas n == (if odd n then negate else id) (lucas (- n))
+
+-- | Check that 'lucasPair' is a pair of consequent 'lucas'.
+lucasPairProperty :: AnySign Int -> Bool
+lucasPairProperty (AnySign n) = lucasPair n == (lucas n, lucas (n + 1))
+
+-- | Check that 'lucas 0' is 2.
+lucasSpecialCase0 :: Assertion
+lucasSpecialCase0 = assertEqual "lucas" (lucas 0) 2
+
+-- | Check that 'lucas 1' is 1.
+lucasSpecialCase1 :: Assertion
+lucasSpecialCase1 = assertEqual "lucas" (lucas 1) 1
+
+-- | Check that 'generalLucas' matches its definition.
+generalLucasProperty1 :: AnySign Integer -> AnySign Integer -> NonNegative Int -> Bool
+generalLucasProperty1 (AnySign p) (AnySign q) (NonNegative n) = un1 == un1' && vn1 == vn1' && un2 == p * un1 - q * un && vn2 == p * vn1 - q * vn
+  where
+    (un, un1, vn, vn1) = generalLucas p q n
+    (un1', un2, vn1', vn2) = generalLucas p q (n + 1)
+
+-- | Check that 'generalLucas' 1 (-1) is 'fibonacciPair' plus 'lucasPair'.
+generalLucasProperty2 :: NonNegative Int -> Bool
+generalLucasProperty2 (NonNegative n) = (un, un1) == fibonacciPair n && (vn, vn1) == lucasPair n
+  where
+    (un, un1, vn, vn1) = generalLucas 1 (-1) n
+
+-- | Check that 'generalLucas' p _ 0 is (0, 1, 2, p).
+generalLucasProperty3 :: AnySign Integer -> AnySign Integer -> Bool
+generalLucasProperty3 (AnySign p) (AnySign q) = generalLucas p q 0 == (0, 1, 2, p)
+
+testSuite :: TestTree
+testSuite = testGroup "Lucas"
+  [ testGroup "fibonacci"
+    [ testSmallAndQuick "matches definition"  fibonacciProperty1
+    , testSmallAndQuick "negative indices"    fibonacciProperty2
+    , testSmallAndQuick "pair"                fibonacciPairProperty
+    , testCase          "fibonacci 0"         fibonacciSpecialCase0
+    , testCase          "fibonacci 1"         fibonacciSpecialCase1
+    ]
+  , testGroup "lucas"
+    [ testSmallAndQuick "matches definition"  lucasProperty1
+    , testSmallAndQuick "negative indices"    lucasProperty2
+    , testSmallAndQuick "pair"                lucasPairProperty
+    , testCase          "lucas 0"             lucasSpecialCase0
+    , testCase          "lucas 1"             lucasSpecialCase1
+    ]
+  , testGroup "generalLucas"
+    [ testSmallAndQuick "matches definition"  generalLucasProperty1
+    , testSmallAndQuick "generalLucas 1 (-1)" generalLucasProperty2
+    , testSmallAndQuick "generalLucas _ _ 0"  generalLucasProperty3
+    ]
+  ]
diff --git a/test-suite/Math/NumberTheory/ModuliTests.hs b/test-suite/Math/NumberTheory/ModuliTests.hs
new file mode 100644
--- /dev/null
+++ b/test-suite/Math/NumberTheory/ModuliTests.hs
@@ -0,0 +1,222 @@
+-- |
+-- Module:      Math.NumberTheory.ModuliTests
+-- Copyright:   (c) 2016 Andrew Lelechenko
+-- Licence:     MIT
+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>
+-- Stability:   Provisional
+--
+-- Tests for Math.NumberTheory.Moduli
+--
+
+{-# LANGUAGE CPP          #-}
+{-# LANGUAGE ViewPatterns #-}
+
+{-# OPTIONS_GHC -fno-warn-type-defaults #-}
+
+module Math.NumberTheory.ModuliTests
+  ( testSuite
+  ) where
+
+import Test.Tasty
+
+import Control.Arrow
+import Data.Bits
+import Data.List (tails, nub)
+import Data.Maybe
+
+import Math.NumberTheory.Moduli
+import Math.NumberTheory.TestUtils
+
+toOdd :: Num a => a -> a
+toOdd n = n * 2 + 1
+
+unwrapPP :: (Prime, Power Int) -> (Integer, Int)
+unwrapPP (Prime p, Power e) = (p, e)
+
+-- | Check that 'jacobi' matches 'jacobi''.
+jacobiProperty1 :: (Integral a, Bits a) => AnySign a -> NonNegative a -> Bool
+jacobiProperty1 (AnySign a) (NonNegative (toOdd -> n)) = n == 1 && j == 1 || n > 1 && j == j'
+  where
+    j = jacobi a n
+    j' = jacobi' a n
+
+-- https://en.wikipedia.org/wiki/Jacobi_symbol#Properties, item 2
+jacobiProperty2 :: (Integral a, Bits a) => AnySign a -> NonNegative a -> Bool
+jacobiProperty2 (AnySign a) (NonNegative (toOdd -> n)) = jacobi a n == jacobi (a + n) n
+
+-- https://en.wikipedia.org/wiki/Jacobi_symbol#Properties, item 3
+jacobiProperty3 :: (Integral a, Bits a) => AnySign a -> NonNegative a -> Bool
+jacobiProperty3 (AnySign a) (NonNegative (toOdd -> n)) = j == 0 && g /= 1 || abs j == 1 && g == 1
+  where
+    j = jacobi a n
+    g = gcd a n
+
+-- https://en.wikipedia.org/wiki/Jacobi_symbol#Properties, item 4
+jacobiProperty4 :: (Integral a, Bits a) => AnySign a -> AnySign a -> NonNegative a -> Bool
+jacobiProperty4 (AnySign a) (AnySign b) (NonNegative (toOdd -> n)) = jacobi (a * b) n == jacobi a n * jacobi b n
+
+jacobiProperty4_Integer :: AnySign Integer -> AnySign Integer -> NonNegative Integer -> Bool
+jacobiProperty4_Integer = jacobiProperty4
+
+-- https://en.wikipedia.org/wiki/Jacobi_symbol#Properties, item 5
+jacobiProperty5 :: (Integral a, Bits a) => AnySign a -> NonNegative a -> NonNegative a -> Bool
+jacobiProperty5 (AnySign a) (NonNegative (toOdd -> m)) (NonNegative (toOdd -> n)) = jacobi a (m * n) == jacobi a m * jacobi a n
+
+jacobiProperty5_Integer :: AnySign Integer -> NonNegative Integer -> NonNegative Integer -> Bool
+jacobiProperty5_Integer = jacobiProperty5
+
+-- https://en.wikipedia.org/wiki/Jacobi_symbol#Properties, item 6
+jacobiProperty6 :: (Integral a, Bits a) => NonNegative a -> NonNegative a -> Bool
+jacobiProperty6 (NonNegative (toOdd -> m)) (NonNegative (toOdd -> n)) = gcd m n /= 1 || jacobi m n * jacobi n m == (if m `mod` 4 == 1 || n `mod` 4 == 1 then 1 else -1)
+
+-- | Check that 'invertMod' inverts numbers modulo.
+invertModProperty :: AnySign Integer -> Positive Integer -> Bool
+invertModProperty (AnySign k) (Positive m) = case invertMod k m of
+  Nothing  -> k `mod` m == 0 || gcd k m > 1
+  Just inv -> gcd k m == 1
+      && k * inv `mod` m == 1 && 0 <= inv && inv < m
+
+-- | Check that the result of 'powerMod' is between 0 and modulo (non-inclusive).
+powerModProperty1 :: (Integral a, Bits a) => AnySign a -> AnySign Integer -> Positive Integer -> Bool
+powerModProperty1 (AnySign e) (AnySign b) (Positive m)
+  =  e < 0 && isNothing (invertMod b m)
+  || (0 <= pm && pm < m)
+  where
+    pm = powerMod b e m
+
+-- | Check that 'powerMod' is multiplicative by first argument.
+powerModProperty2 :: (Integral a, Bits a) => AnySign a -> AnySign Integer -> AnySign Integer -> Positive Integer -> Bool
+powerModProperty2 (AnySign e) (AnySign b1) (AnySign b2) (Positive m)
+  =  e < 0 && (isNothing (invertMod b1 m) || isNothing (invertMod b2 m))
+  || pm1 * pm2 `mod` m == pm12
+  where
+    pm1  = powerMod b1  e m
+    pm2  = powerMod b2  e m
+    pm12 = powerMod (b1 * b2) e m
+
+-- | Check that 'powerMod' is additive by second argument.
+powerModProperty3 :: (Integral a, Bits a) => AnySign a -> AnySign a -> AnySign Integer -> Positive Integer -> Bool
+powerModProperty3 (AnySign e1) (AnySign e2) (AnySign b) (Positive m)
+  =  (e1 < 0 || e2 < 0) && isNothing (invertMod b m)
+  || pm1 * pm2 `mod` m == pm12
+  where
+    pm1  = powerMod b e1 m
+    pm2  = powerMod b e2 m
+    pm12 = powerMod b (e1 + e2) m
+
+-- | Specialized to trigger 'powerModInteger'.
+powerModProperty1_Integer :: AnySign Integer -> AnySign Integer -> Positive Integer -> Bool
+powerModProperty1_Integer = powerModProperty1
+
+-- | Specialized to trigger 'powerModInteger'.
+powerModProperty2_Integer :: AnySign Integer -> AnySign Integer -> AnySign Integer -> Positive Integer -> Bool
+powerModProperty2_Integer = powerModProperty2
+
+-- | Specialized to trigger 'powerModInteger'.
+powerModProperty3_Integer :: AnySign Integer -> AnySign Integer -> AnySign Integer -> Positive Integer -> Bool
+powerModProperty3_Integer = powerModProperty3
+
+-- | Check that 'powerMod' matches 'powerMod''.
+powerMod'Property :: (Integral a, Bits a) => Positive a -> Positive Integer -> Positive Integer -> Bool
+powerMod'Property (Positive e) (Positive b) (Positive m) = m == 1 || powerMod' b e m == powerMod b e m
+
+-- | Specialized to trigger 'powerModInteger''.
+powerMod'Property_Integer :: Positive Integer -> Positive Integer -> Positive Integer -> Bool
+powerMod'Property_Integer = powerMod'Property
+
+-- | Check that 'chineseRemainder' is defined iff modulos are coprime.
+--   Also check that the result is a solution of input modular equations.
+chineseRemainderProperty :: [(Integer, Positive Integer)] -> Bool
+chineseRemainderProperty rms' = case chineseRemainder rms of
+  Nothing -> not areCoprime
+  Just n  -> areCoprime && map (n `mod`) ms == zipWith mod rs ms
+  where
+    rms = map (second getPositive) rms'
+    (rs, ms) = unzip rms
+    areCoprime = all (== 1) [ gcd m1 m2 | (m1 : m2s) <- tails ms, m2 <- m2s ]
+
+-- | Check that 'chineseRemainder' matches 'chineseRemainder2'.
+chineseRemainder2Property :: Integer -> Positive Integer -> Integer -> Positive Integer -> Bool
+chineseRemainder2Property r1 (Positive m1) r2 (Positive m2) = gcd m1 m2 /= 1
+  || Just (chineseRemainder2 (r1, m1) (r2, m2)) == chineseRemainder [(r1, m1), (r2, m2)]
+
+-- | Check that 'sqrtMod' is defined iff a quadratic residue exists.
+--   Also check that the result is a solution of input modular equation.
+sqrtModPProperty :: AnySign Integer -> Prime -> Bool
+sqrtModPProperty (AnySign n) (Prime p) = case sqrtModP n p of
+  Nothing -> jacobi n p == -1
+  Just rt -> (p == 2 || jacobi n p /= -1) && rt ^ 2 `mod` p == n `mod` p
+
+sqrtModPListProperty :: AnySign Integer -> Prime -> Bool
+sqrtModPListProperty (AnySign n) (Prime p) = all (\rt -> rt ^ 2 `mod` p == n `mod` p) (sqrtModPList n p)
+
+sqrtModP'Property :: Positive Integer -> Prime -> Bool
+sqrtModP'Property (Positive n) (Prime p) = (p /= 2 && jacobi n p /= 1) || rt ^ 2 `mod` p == n `mod` p
+  where
+    rt = sqrtModP' n p
+
+tonelliShanksProperty :: Positive Integer -> Prime -> Bool
+tonelliShanksProperty (Positive n) (Prime p) = p `mod` 4 /= 1 || jacobi n p /= 1 || rt ^ 2 `mod` p == n `mod` p
+  where
+    rt = tonelliShanks n p
+
+sqrtModPPProperty :: AnySign Integer -> (Prime, Power Int) -> Bool
+sqrtModPPProperty (AnySign n) (Prime p, Power e) = gcd n p > 1 || case sqrtModPP n (p, e) of
+  Nothing -> True
+  Just rt -> rt ^ 2 `mod` (p ^ e) == n `mod` (p ^ e)
+
+sqrtModPPListProperty :: AnySign Integer -> (Prime, Power Int) -> Bool
+sqrtModPPListProperty (AnySign n) (Prime p, Power e) = gcd n p > 1
+  || all (\rt -> rt ^ 2 `mod` (p ^ e) == n `mod` (p ^ e)) (sqrtModPPList n (p, e))
+
+sqrtModFProperty :: AnySign Integer -> [(Prime, Power Int)] -> Bool
+sqrtModFProperty (AnySign n) (map unwrapPP -> pes) = case sqrtModF n pes of
+  Nothing -> True
+  Just rt -> all (\(p, e) -> rt ^ 2 `mod` (p ^ e) == n `mod` (p ^ e)) pes
+
+sqrtModFListProperty :: AnySign Integer -> [(Prime, Power Int)] -> Bool
+sqrtModFListProperty (AnySign n) (map unwrapPP -> pes)
+  = nub ps /= ps || all
+    (\rt -> all (\(p, e) -> rt ^ 2 `mod` (p ^ e) == n `mod` (p ^ e)) pes)
+    (sqrtModFList n pes)
+  where
+    ps = map fst pes
+
+testSuite :: TestTree
+testSuite = testGroup "Moduli"
+  [ testGroup "jacobi"
+    [ testSameIntegralProperty "matches jacobi'"              jacobiProperty1
+    , testSameIntegralProperty "same modulo n"                jacobiProperty2
+    , testSameIntegralProperty "consistent with gcd"          jacobiProperty3
+    , testSmallAndQuick        "multiplicative 1"             jacobiProperty4_Integer
+    , testSmallAndQuick        "multiplicative 2"             jacobiProperty5_Integer
+    , testSameIntegralProperty "law of quadratic reciprocity" jacobiProperty6
+    ]
+  , testSmallAndQuick "invertMod" invertModProperty
+  , testGroup "powerMod"
+    [ testGroup "generic"
+      [ testIntegralProperty "bounded between 0 and m"  powerModProperty1
+      , testIntegralProperty "multiplicative by base"   powerModProperty2
+      , testSameIntegralProperty "additive by exponent" powerModProperty3
+      , testIntegralProperty "matches powerMod'"        powerMod'Property
+      ]
+    , testGroup "Integer"
+      [ testSmallAndQuick "bounded between 0 and m" powerModProperty1_Integer
+      , testSmallAndQuick "multiplicative by base"  powerModProperty2_Integer
+      , testSmallAndQuick "additive by exponent"    powerModProperty3_Integer
+      , testSmallAndQuick "matches powerMod'"       powerMod'Property_Integer
+      ]
+    ]
+    , testSmallAndQuick "chineseRemainder"  chineseRemainderProperty
+    , testSmallAndQuick "chineseRemainder2" chineseRemainder2Property
+    , testGroup "sqrtMod"
+      [ testSmallAndQuick "sqrtModP"      sqrtModPProperty
+      , testSmallAndQuick "sqrtModPList"  sqrtModPListProperty
+      , testSmallAndQuick "sqrtModP'"     sqrtModP'Property
+      , testSmallAndQuick "tonelliShanks" tonelliShanksProperty
+      , testSmallAndQuick "sqrtModPP"     sqrtModPPProperty
+      , testSmallAndQuick "sqrtModPPList" sqrtModPPListProperty
+      , testSmallAndQuick "sqrtModF"      sqrtModFProperty
+      , testSmallAndQuick "sqrtModFList"  sqrtModFListProperty
+      ]
+  ]
diff --git a/test-suite/Math/NumberTheory/MoebiusInversion/IntTests.hs b/test-suite/Math/NumberTheory/MoebiusInversion/IntTests.hs
new file mode 100644
--- /dev/null
+++ b/test-suite/Math/NumberTheory/MoebiusInversion/IntTests.hs
@@ -0,0 +1,45 @@
+-- |
+-- Module:      Math.NumberTheory.MoebiusInversion.IntTests
+-- Copyright:   (c) 2016 Andrew Lelechenko
+-- Licence:     MIT
+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>
+-- Stability:   Provisional
+--
+-- Tests for Math.NumberTheory.MoebiusInversion.Int
+--
+
+{-# OPTIONS_GHC -fno-warn-type-defaults #-}
+
+module Math.NumberTheory.MoebiusInversion.IntTests
+  ( testSuite
+  ) where
+
+import Test.Tasty
+import Test.Tasty.HUnit
+import Test.Tasty.QuickCheck as QC hiding (Positive)
+
+import Math.NumberTheory.MoebiusInversion.Int
+import Math.NumberTheory.Primes.Factorisation
+import Math.NumberTheory.TestUtils
+
+totientSumProperty :: Positive Int -> Bool
+totientSumProperty (Positive n) = toInteger (totientSum n) == sum (map totient [1 .. toInteger n])
+
+totientSumSpecialCase1 :: Assertion
+totientSumSpecialCase1 = assertEqual "totientSum" 4496 (totientSum 121)
+
+generalInversionProperty :: (Int -> Int) -> Positive Int -> Bool
+generalInversionProperty g (Positive n)
+  =  g n == sum [f (n `quot` k) | k <- [1 .. n]]
+  && f n == sum [fromInteger (moebius (toInteger k)) * g (n `quot` k) | k <- [1 .. n]]
+  where
+    f = generalInversion g
+
+testSuite :: TestTree
+testSuite = testGroup "Int"
+  [ testGroup "totientSum"
+    [ testSmallAndQuick "matches definitions" totientSumProperty
+    , testCase          "special case 1"      totientSumSpecialCase1
+    ]
+  , QC.testProperty "generalInversion" generalInversionProperty
+  ]
diff --git a/test-suite/Math/NumberTheory/MoebiusInversionTests.hs b/test-suite/Math/NumberTheory/MoebiusInversionTests.hs
new file mode 100644
--- /dev/null
+++ b/test-suite/Math/NumberTheory/MoebiusInversionTests.hs
@@ -0,0 +1,45 @@
+-- |
+-- Module:      Math.NumberTheory.MoebiusInversionTests
+-- Copyright:   (c) 2016 Andrew Lelechenko
+-- Licence:     MIT
+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>
+-- Stability:   Provisional
+--
+-- Tests for Math.NumberTheory.MoebiusInversion
+--
+
+{-# OPTIONS_GHC -fno-warn-type-defaults #-}
+
+module Math.NumberTheory.MoebiusInversionTests
+  ( testSuite
+  ) where
+
+import Test.Tasty
+import Test.Tasty.HUnit
+import Test.Tasty.QuickCheck as QC hiding (Positive)
+
+import Math.NumberTheory.MoebiusInversion
+import Math.NumberTheory.Primes.Factorisation
+import Math.NumberTheory.TestUtils
+
+totientSumProperty :: Positive Int -> Bool
+totientSumProperty (Positive n) = totientSum n == sum (map totient [1 .. toInteger n])
+
+totientSumSpecialCase1 :: Assertion
+totientSumSpecialCase1 = assertEqual "totientSum" 4496 (totientSum 121)
+
+generalInversionProperty :: (Int -> Integer) -> Positive Int -> Bool
+generalInversionProperty g (Positive n)
+  =  g n == sum [f (n `quot` k) | k <- [1 .. n]]
+  && f n == sum [moebius (toInteger k) * g (n `quot` k) | k <- [1 .. n]]
+  where
+    f = generalInversion g
+
+testSuite :: TestTree
+testSuite = testGroup "MoebiusInversion"
+  [ testGroup "totientSum"
+    [ testSmallAndQuick "matches definitions" totientSumProperty
+    , testCase          "special case 1"      totientSumSpecialCase1
+    ]
+  , QC.testProperty "generalInversion" generalInversionProperty
+  ]
diff --git a/test-suite/Math/NumberTheory/Powers/CubesTests.hs b/test-suite/Math/NumberTheory/Powers/CubesTests.hs
new file mode 100644
--- /dev/null
+++ b/test-suite/Math/NumberTheory/Powers/CubesTests.hs
@@ -0,0 +1,154 @@
+-- |
+-- Module:      Math.NumberTheory.Powers.CubesTests
+-- Copyright:   (c) 2016 Andrew Lelechenko
+-- Licence:     MIT
+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>
+-- Stability:   Provisional
+--
+-- Tests for Math.NumberTheory.Powers.Cubes
+--
+
+{-# LANGUAGE CPP #-}
+
+{-# OPTIONS_GHC -fno-warn-type-defaults #-}
+
+module Math.NumberTheory.Powers.CubesTests
+  ( testSuite
+  ) where
+
+import Test.Tasty
+import Test.Tasty.HUnit
+
+import Data.Maybe
+#if MIN_VERSION_base(4,8,0)
+#else
+import Data.Word
+#endif
+
+
+import Math.NumberTheory.Powers.Cubes
+import Math.NumberTheory.TestUtils
+
+#include "MachDeps.h"
+
+-- | Check that 'integerCubeRoot' returns the largest integer @m@ with @m^3 <= n@.
+--
+-- (m + 1) ^ 3 /= n && cond
+-- means
+-- (m + 1) ^ 3 > n
+-- but without overflow for bounded types
+integerCubeRootProperty :: Integral a => AnySign a -> Bool
+integerCubeRootProperty (AnySign n) = m ^ 3 <= n && (m + 1) ^ 3 /= n && cond
+  where
+    m = integerCubeRoot n
+    cond
+      | m == -1   = n == -1
+      | m < 0     = (m + 1) ^ 2 <= n `div` (m + 1)
+      | otherwise = (m + 1) ^ 2 >= n `div` (m + 1)
+
+-- | Specialized to trigger 'cubeRootInt''.
+integerCubeRootProperty_Int :: AnySign Int -> Bool
+integerCubeRootProperty_Int = integerCubeRootProperty
+
+-- | Specialized to trigger 'cubeRootWord'.
+integerCubeRootProperty_Word :: AnySign Word -> Bool
+integerCubeRootProperty_Word = integerCubeRootProperty
+
+-- | Specialized to trigger 'cubeRootIgr'.
+integerCubeRootProperty_Integer :: AnySign Integer -> Bool
+integerCubeRootProperty_Integer = integerCubeRootProperty
+
+-- | Check that 'integerCubeRoot' returns the largest integer @m@ with @m^3 <= n@, , where @n@ has form @k@^3-1.
+integerCubeRootProperty2 :: Integral a => AnySign a -> Bool
+integerCubeRootProperty2 (AnySign k) = m ^ 3 <= n && (m + 1) ^ 3 /= n && cond
+  where
+    n = k ^ 3 - 1
+    m = integerCubeRoot n
+    cond
+      | m == -1   = n == -1
+      | m < 0     = (m + 1) ^ 2 <= n `div` (m + 1)
+      | otherwise = (m + 1) ^ 2 >= n `div` (m + 1)
+
+-- | Specialized to trigger 'cubeRootInt''.
+integerCubeRootProperty2_Int :: AnySign Int -> Bool
+integerCubeRootProperty2_Int = integerCubeRootProperty2
+
+-- | Specialized to trigger 'cubeRootWord'.
+integerCubeRootProperty2_Word :: AnySign Word -> Bool
+integerCubeRootProperty2_Word = integerCubeRootProperty2
+
+#if WORD_SIZE_IN_BITS == 64
+
+-- | Check that 'integerCubeRoot' of 2^63-1 is 2^21-1, not 2^21.
+integerCubeRootSpecialCase1_Int :: Assertion
+integerCubeRootSpecialCase1_Int =
+  assertEqual "integerCubeRoot" (integerCubeRoot (maxBound :: Int)) (2 ^ 21 - 1)
+
+-- | Check that 'integerCubeRoot' of 2^63-1 is 2^21-1, not 2^21.
+integerCubeRootSpecialCase1_Word :: Assertion
+integerCubeRootSpecialCase1_Word =
+  assertEqual "integerCubeRoot" (integerCubeRoot (maxBound `div` 2 :: Word)) (2 ^ 21 - 1)
+
+-- | Check that 'integerCubeRoot' of 2^64-1 is 2642245.
+integerCubeRootSpecialCase2 :: Assertion
+integerCubeRootSpecialCase2 =
+  assertEqual "integerCubeRoot" (integerCubeRoot (maxBound :: Word)) 2642245
+
+#endif
+
+-- | Check that 'integerCubeRoot'' returns the largest integer @m@ with @m^3 <= n@.
+integerCubeRoot'Property :: Integral a => NonNegative a -> Bool
+integerCubeRoot'Property (NonNegative n) = m ^ 3 <= n && (m + 1) ^ 3 /= n && (m + 1) ^ 2 >= n `div` (m + 1)
+  where
+    m = integerCubeRoot' n
+
+-- | Check that the number 'isCube' iff its 'integerCubeRoot' is exact.
+isCubeProperty :: Integral a => AnySign a -> Bool
+isCubeProperty (AnySign n) = (n /= m ^ 3 && not t) || (n == m ^ 3 && t)
+  where
+    t = isCube n
+    m = integerCubeRoot n
+
+-- | Check that the number 'isCube'' iff its 'integerCubeRoot'' is exact.
+isCube'Property :: Integral a => NonNegative a -> Bool
+isCube'Property (NonNegative n) = (n /= m ^ 3 && not t) || (n == m ^ 3 && t)
+  where
+    t = isCube' n
+    m = integerCubeRoot' n
+
+-- | Check that 'exactCubeRoot' returns an exact integer cubic root
+-- and is consistent with 'isCube'.
+exactCubeRootProperty :: Integral a => AnySign a -> Bool
+exactCubeRootProperty (AnySign n) = case exactCubeRoot n of
+  Nothing -> not (isCube n)
+  Just m  -> isCube n && n == m ^ 3
+
+-- | Check that 'isPossibleCube' is consistent with 'exactCubeRoot'.
+isPossibleCubeProperty :: Integral a => NonNegative a -> Bool
+isPossibleCubeProperty (NonNegative n) = t || not t && isNothing m
+  where
+    t = isPossibleCube n
+    m = exactCubeRoot n
+
+testSuite :: TestTree
+testSuite = testGroup "Cubes"
+  [ testGroup "integerCubeRoot"
+    [ testIntegralProperty "generic"         integerCubeRootProperty
+    , testSmallAndQuick    "generic Int"     integerCubeRootProperty_Int
+    , testSmallAndQuick    "generic Word"    integerCubeRootProperty_Word
+    , testSmallAndQuick    "generic Integer" integerCubeRootProperty_Integer
+
+    , testIntegralProperty "almost cube"      integerCubeRootProperty2
+    , testSmallAndQuick    "almost cube Int"  integerCubeRootProperty2_Int
+    , testSmallAndQuick    "almost cube Word" integerCubeRootProperty2_Word
+
+    , testCase             "maxBound :: Int"      integerCubeRootSpecialCase1_Int
+    , testCase             "maxBound / 2 :: Word" integerCubeRootSpecialCase1_Word
+    , testCase             "maxBound :: Word"     integerCubeRootSpecialCase2
+    ]
+  , testIntegralProperty "integerCubeRoot'" integerCubeRoot'Property
+  , testIntegralProperty "isCube"           isCubeProperty
+  , testIntegralProperty "isCube'"          isCube'Property
+  , testIntegralProperty "exactCubeRoot"    exactCubeRootProperty
+  , testIntegralProperty "isPossibleCube"   isPossibleCubeProperty
+  ]
diff --git a/test-suite/Math/NumberTheory/Powers/FourthTests.hs b/test-suite/Math/NumberTheory/Powers/FourthTests.hs
new file mode 100644
--- /dev/null
+++ b/test-suite/Math/NumberTheory/Powers/FourthTests.hs
@@ -0,0 +1,145 @@
+-- |
+-- Module:      Math.NumberTheory.Powers.FourthTests
+-- Copyright:   (c) 2016 Andrew Lelechenko
+-- Licence:     MIT
+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>
+-- Stability:   Provisional
+--
+-- Tests for Math.NumberTheory.Powers.Fourth
+--
+
+{-# LANGUAGE CPP #-}
+
+{-# OPTIONS_GHC -fno-warn-type-defaults #-}
+
+module Math.NumberTheory.Powers.FourthTests
+  ( testSuite
+  ) where
+
+import Test.Tasty
+import Test.Tasty.HUnit
+
+import Data.Maybe
+#if MIN_VERSION_base(4,8,0)
+#else
+import Data.Word
+#endif
+
+import Math.NumberTheory.Powers.Fourth
+import Math.NumberTheory.TestUtils
+
+#include "MachDeps.h"
+
+-- | Check that 'integerFourthRoot' returns the largest integer @m@ with @m^4 <= n@.
+--
+-- (m + 1) ^ 4 /= n && (m + 1) ^ 3 >= n `div` (m + 1)
+-- means
+-- (m + 1) ^ 4 > n
+-- but without overflow for bounded types
+integerFourthRootProperty :: Integral a => NonNegative a -> Bool
+integerFourthRootProperty (NonNegative n) = m >= 0 && m ^ 4 <= n && (m + 1) ^ 4 /= n && (m + 1) ^ 3 >= n `div` (m + 1)
+  where
+    m = integerFourthRoot n
+
+-- | Specialized to trigger 'biSqrtInt'.
+integerFourthRootProperty_Int :: NonNegative Int -> Bool
+integerFourthRootProperty_Int = integerFourthRootProperty
+
+-- | Specialized to trigger 'biSqrtWord'.
+integerFourthRootProperty_Word :: NonNegative Word -> Bool
+integerFourthRootProperty_Word = integerFourthRootProperty
+
+-- | Specialized to trigger 'biSqrtIgr'.
+integerFourthRootProperty_Integer :: NonNegative Integer -> Bool
+integerFourthRootProperty_Integer = integerFourthRootProperty
+
+-- | Check that 'integerFourthRoot' returns the largest integer @m@ with @m^4 <= n@, , where @n@ has form @k@^4-1.
+integerFourthRootProperty2 :: Integral a => NonNegative a -> Bool
+integerFourthRootProperty2 (NonNegative k) = n < 0 || m >= 0 && m ^ 4 <= n && (m + 1) ^ 4 /= n && (m + 1) ^ 3 >= n `div` (m + 1)
+  where
+    n = k ^ 4 - 1
+    m = integerFourthRoot n
+
+-- | Specialized to trigger 'biSqrtInt.
+integerFourthRootProperty2_Int :: NonNegative Int -> Bool
+integerFourthRootProperty2_Int = integerFourthRootProperty2
+
+-- | Specialized to trigger 'biSqrtWord'.
+integerFourthRootProperty2_Word :: NonNegative Word -> Bool
+integerFourthRootProperty2_Word = integerFourthRootProperty2
+
+#if WORD_SIZE_IN_BITS == 64
+
+-- | Check that 'integerFourthRoot' of 2^60-1 is 2^15-1, not 2^15.
+integerFourthRootSpecialCase1_Int :: Assertion
+integerFourthRootSpecialCase1_Int =
+  assertEqual "integerFourthRoot" (integerFourthRoot (maxBound `div` 8 :: Int)) (2 ^ 15 - 1)
+
+-- | Check that 'integerFourthRoot' of 2^60-1 is 2^15-1, not 2^15.
+integerFourthRootSpecialCase1_Word :: Assertion
+integerFourthRootSpecialCase1_Word =
+  assertEqual "integerFourthRoot" (integerFourthRoot (maxBound `div` 16 :: Word)) (2 ^ 15 - 1)
+
+-- | Check that 'integerFourthRoot' of 2^64-1 is 2^16-1, not 2^16.
+integerFourthRootSpecialCase2 :: Assertion
+integerFourthRootSpecialCase2 =
+  assertEqual "integerFourthRoot" (integerFourthRoot (maxBound :: Word)) (2 ^ 16 - 1)
+
+#endif
+
+-- | Check that 'integerFourthRoot'' returns the largest integer @m@ with @m^4 <= n@.
+integerFourthRoot'Property :: Integral a => NonNegative a -> Bool
+integerFourthRoot'Property (NonNegative n) = m >= 0 && m ^ 4 <= n && (m + 1) ^ 4 /= n && (m + 1) ^ 3 >= n `div` (m + 1)
+  where
+    m = integerFourthRoot' n
+
+-- | Check that the number 'isFourthPower' iff its 'integerFourthRoot' is exact.
+isFourthPowerProperty :: Integral a => AnySign a -> Bool
+isFourthPowerProperty (AnySign n) = (n < 0 && not t) || (n /= m ^ 4 && not t) || (n == m ^ 4 && t)
+  where
+    t = isFourthPower n
+    m = integerFourthRoot n
+
+-- | Check that the number 'isFourthPower'' iff its 'integerFourthRoot'' is exact.
+isFourthPower'Property :: Integral a => NonNegative a -> Bool
+isFourthPower'Property (NonNegative n) = (n /= m ^ 4 && not t) || (n == m ^ 4 && t)
+  where
+    t = isFourthPower' n
+    m = integerFourthRoot' n
+
+-- | Check that 'exactFourthRoot' returns an exact integer root of fourth power
+-- and is consistent with 'isFourthPower'.
+exactFourthRootProperty :: Integral a => AnySign a -> Bool
+exactFourthRootProperty (AnySign n) = case exactFourthRoot n of
+  Nothing -> not (isFourthPower n)
+  Just m  -> isFourthPower n && n == m ^ 4
+
+-- | Check that 'isPossibleFourthPower' is consistent with 'exactFourthRoot'.
+isPossibleFourthPowerProperty :: Integral a => NonNegative a -> Bool
+isPossibleFourthPowerProperty (NonNegative n) = t || not t && isNothing m
+  where
+    t = isPossibleFourthPower n
+    m = exactFourthRoot n
+
+testSuite :: TestTree
+testSuite = testGroup "Fourth"
+  [ testGroup "integerFourthRoot"
+    [ testIntegralProperty "generic"         integerFourthRootProperty
+    , testSmallAndQuick    "generic Int"     integerFourthRootProperty_Int
+    , testSmallAndQuick    "generic Word"    integerFourthRootProperty_Word
+    , testSmallAndQuick    "generic Integer" integerFourthRootProperty_Integer
+
+    , testIntegralProperty "almost Fourth"      integerFourthRootProperty2
+    , testSmallAndQuick    "almost Fourth Int"  integerFourthRootProperty2_Int
+    , testSmallAndQuick    "almost Fourth Word" integerFourthRootProperty2_Word
+
+    , testCase             "maxBound / 8 :: Int"   integerFourthRootSpecialCase1_Int
+    , testCase             "maxBound / 16 :: Word" integerFourthRootSpecialCase1_Word
+    , testCase             "maxBound :: Word"      integerFourthRootSpecialCase2
+    ]
+  , testIntegralProperty "integerFourthRoot'"    integerFourthRoot'Property
+  , testIntegralProperty "isFourthPower"         isFourthPowerProperty
+  , testIntegralProperty "isFourthPower'"        isFourthPower'Property
+  , testIntegralProperty "exactFourthRoot"       exactFourthRootProperty
+  , testIntegralProperty "isPossibleFourthPower" isPossibleFourthPowerProperty
+  ]
diff --git a/test-suite/Math/NumberTheory/Powers/GeneralTests.hs b/test-suite/Math/NumberTheory/Powers/GeneralTests.hs
new file mode 100644
--- /dev/null
+++ b/test-suite/Math/NumberTheory/Powers/GeneralTests.hs
@@ -0,0 +1,128 @@
+-- |
+-- Module:      Math.NumberTheory.Powers.GeneralTests
+-- Copyright:   (c) 2016 Andrew Lelechenko
+-- Licence:     MIT
+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>
+-- Stability:   Provisional
+--
+-- Tests for Math.NumberTheory.Powers.General
+--
+
+{-# LANGUAGE CPP #-}
+
+{-# OPTIONS_GHC -fno-warn-type-defaults #-}
+
+module Math.NumberTheory.Powers.GeneralTests
+  ( testSuite
+  ) where
+
+import Test.Tasty
+import Test.Tasty.HUnit
+
+import Math.NumberTheory.Powers.General
+import Math.NumberTheory.TestUtils
+
+-- | Check that 'integerRoot' @pow@ returns the largest integer @m@ with @m^pow <= n@.
+integerRootProperty :: (Integral a, Integral b) => AnySign a -> Power b -> Bool
+integerRootProperty (AnySign n) (Power pow) = (even pow && n < 0)
+  || (toInteger root ^ pow <= toInteger n && toInteger n < toInteger (root + 1) ^ pow)
+    where
+      root = integerRoot pow n
+
+-- | Check that the number 'isKthPower' iff its 'integerRoot' is exact.
+isKthPowerProperty :: (Integral a, Integral b) => AnySign a -> Power b -> Bool
+isKthPowerProperty (AnySign n) (Power pow) = (even pow && n < 0 && not t) || (n /= root ^ pow && not t) || (n == root ^ pow && t)
+  where
+    t = isKthPower pow n
+    root = integerRoot pow n
+
+-- | Check that 'exactRoot' returns an exact integer root
+-- and is consistent with 'isKthPower'.
+exactRootProperty :: (Integral a, Integral b) => AnySign a -> Power b -> Bool
+exactRootProperty (AnySign n) (Power pow) = case exactRoot pow n of
+  Nothing   -> not (isKthPower pow n)
+  Just root -> isKthPower pow n && n == root ^ pow
+
+-- | Check that 'isPerfectPower' is consistent with 'highestPower'.
+isPerfectPowerProperty :: Integral a => AnySign a -> Bool
+isPerfectPowerProperty (AnySign n) = (k > 1 && t) || (k == 1 && not t)
+  where
+    t = isPerfectPower n
+    (_, k) = highestPower n
+
+-- | Check that the first component of 'highestPower' is square-free.
+highestPowerProperty :: Integral a => AnySign a -> Bool
+highestPowerProperty (AnySign n) = (n `elem` [-1, 0, 1] && k == 3) || (b ^ k == n && b' == b && k' == 1)
+  where
+    (b, k) = highestPower n
+    (b', k') = highestPower b
+
+-- | Check that 'largePFPower' is consistent with documentation.
+largePFPowerProperty :: Positive Integer -> Integer -> Bool
+largePFPowerProperty (Positive bd) n = bd == 1 || b == 0 || d' /= 0 || n <= b * d * d || any (\p -> gcd n p > 1) [2..bd] || b ^ k == n
+  where
+    (b, k) = largePFPower bd n
+    (d, d') = bd `quotRem` b
+
+highestPowerSpecialCases :: [Assertion]
+highestPowerSpecialCases =
+  -- Freezes before d44a13b.
+  [ a ( 1013582159576576
+      , 1013582159576576
+      , 1)
+  -- Freezes before d44a13b.
+  , a ( 1013582159576576 ^ 7
+      , 1013582159576576
+      , 7)
+
+  , a ( -2 ^ 63 :: Int
+      , -2 :: Int
+      , 63)
+
+  , a ( (2 ^ 63 - 1) ^ 21
+      , 2 ^ 63 - 1
+      , 21)
+
+  , a ( 576116746989720969230211509779286598589421531472851155101032940901763389787901933902294777750323196846498573545522289802689311975294763847414975335235584
+      , 576116746989720969230211509779286598589421531472851155101032940901763389787901933902294777750323196846498573545522289802689311975294763847414975335235584
+      , 1)
+
+  , a ( -340282366920938463500268095579187314689
+      , -340282366920938463500268095579187314689
+      , 1)
+
+  , a ( 268398749 :: Int
+      , 268398749 :: Int
+      , 1)
+
+  , a ( 118372752099 :: Int
+      , 118372752099 :: Int
+      , 1)
+
+  , a ( 1409777209 :: Int
+      , 37547 :: Int
+      , 2)
+
+  , a ( -6277101735386680764856636523970481806547819498980467802113
+      , -18446744073709551617
+      , 3)
+
+  , a ( -18446744073709551619 ^ 5
+      , -18446744073709551619
+      , 5)
+  ]
+  where
+    a (n, b, k) = assertEqual "highestPower" (b, k) (highestPower n)
+
+testSuite :: TestTree
+testSuite = testGroup "General"
+  [ testIntegral2Property "integerRoot"    integerRootProperty
+  , testIntegral2Property "isKthPower"     isKthPowerProperty
+  , testIntegral2Property "exactRoot"      exactRootProperty
+  , testIntegralProperty  "isPerfectPower" isPerfectPowerProperty
+  , testGroup "highestPower"
+    ( testIntegralProperty  "highestPower"   highestPowerProperty
+    : zipWith (\i a -> testCase ("special case " ++ show i) a) [1..] highestPowerSpecialCases
+    )
+  , testSmallAndQuick     "largePFPower"   largePFPowerProperty
+  ]
diff --git a/test-suite/Math/NumberTheory/Powers/IntegerTests.hs b/test-suite/Math/NumberTheory/Powers/IntegerTests.hs
new file mode 100644
--- /dev/null
+++ b/test-suite/Math/NumberTheory/Powers/IntegerTests.hs
@@ -0,0 +1,41 @@
+-- |
+-- Module:      Math.NumberTheory.Powers.IntegerTests
+-- Copyright:   (c) 2016 Andrew Lelechenko
+-- Licence:     MIT
+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>
+-- Stability:   Provisional
+--
+-- Tests for Math.NumberTheory.Powers.Integer
+--
+
+{-# LANGUAGE CPP #-}
+
+{-# OPTIONS_GHC -fno-warn-type-defaults #-}
+
+module Math.NumberTheory.Powers.IntegerTests
+  ( testSuite
+  ) where
+
+import Test.Tasty
+
+#if MIN_VERSION_base(4,8,0)
+#else
+import Data.Word
+#endif
+
+import Math.NumberTheory.Powers.Integer
+import Math.NumberTheory.TestUtils
+
+-- | Check that 'integerPower' == '^'.
+integerPowerProperty :: Integer -> Power Int -> Bool
+integerPowerProperty a (Power b) = integerPower a b == a ^ b
+
+-- | Check that 'integerWordPower' == '^'.
+integerWordPowerProperty :: Integer -> Power Word -> Bool
+integerWordPowerProperty a (Power b) = integerWordPower a b == a ^ b
+
+testSuite :: TestTree
+testSuite = testGroup "Integer"
+  [ testSmallAndQuick "integerPower"     integerPowerProperty
+  , testSmallAndQuick "integerWordPower" integerWordPowerProperty
+  ]
diff --git a/test-suite/Math/NumberTheory/Powers/SquaresTests.hs b/test-suite/Math/NumberTheory/Powers/SquaresTests.hs
new file mode 100644
--- /dev/null
+++ b/test-suite/Math/NumberTheory/Powers/SquaresTests.hs
@@ -0,0 +1,155 @@
+-- |
+-- Module:      Math.NumberTheory.Powers.SquaresTests
+-- Copyright:   (c) 2016 Andrew Lelechenko
+-- Licence:     MIT
+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>
+-- Stability:   Provisional
+--
+-- Tests for Math.NumberTheory.Powers.Squares
+--
+
+{-# LANGUAGE CPP #-}
+
+{-# OPTIONS_GHC -fno-warn-type-defaults #-}
+
+module Math.NumberTheory.Powers.SquaresTests
+  ( testSuite
+  ) where
+
+import Test.Tasty
+import Test.Tasty.HUnit
+
+import Data.Maybe
+#if MIN_VERSION_base(4,8,0)
+#else
+import Data.Word
+#endif
+
+import Math.NumberTheory.Powers.Squares
+import Math.NumberTheory.TestUtils
+
+#include "MachDeps.h"
+
+-- | Check that 'integerSquareRoot' returns the largest integer @m@ with @m*m <= n@.
+--
+-- (m + 1) ^ 2 /= n && m + 1 >= n `div` (m + 1)
+-- means
+-- (m + 1) ^ 2 > n
+-- but without overflow for bounded types
+integerSquareRootProperty :: Integral a => NonNegative a -> Bool
+integerSquareRootProperty (NonNegative n) = m >=0 && m * m <= n && (m + 1) ^ 2 /= n && m + 1 >= n `div` (m + 1)
+  where
+    m = integerSquareRoot n
+
+-- | Specialized to trigger 'isqrtInt''.
+integerSquareRootProperty_Int :: NonNegative Int -> Bool
+integerSquareRootProperty_Int = integerSquareRootProperty
+
+-- | Specialized to trigger 'isqrtWord'.
+integerSquareRootProperty_Word :: NonNegative Word -> Bool
+integerSquareRootProperty_Word = integerSquareRootProperty
+
+-- | Check that 'integerSquareRoot' returns the largest integer @m@ with @m*m <= n@, where @n@ has form @k@^2-1.
+integerSquareRootProperty2 :: Integral a => NonNegative a -> Bool
+integerSquareRootProperty2 (NonNegative k) = n < 0
+  || m >=0 && m * m <= n && (m + 1) ^ 2 /= n && m + 1 >= n `div` (m + 1)
+  where
+    n = k ^ 2 - 1
+    m = integerSquareRoot n
+
+-- | Specialized to trigger 'isqrtInt''.
+integerSquareRootProperty2_Int :: NonNegative Int -> Bool
+integerSquareRootProperty2_Int = integerSquareRootProperty2
+
+-- | Specialized to trigger 'isqrtWord'.
+integerSquareRootProperty2_Word :: NonNegative Word -> Bool
+integerSquareRootProperty2_Word = integerSquareRootProperty2
+
+#if WORD_SIZE_IN_BITS == 64
+
+-- | Check that 'integerSquareRoot' of 2^62-1 is 2^31-1, not 2^31.
+integerSquareRootSpecialCase1_Int :: Assertion
+integerSquareRootSpecialCase1_Int =
+  assertEqual "integerSquareRoot" (integerSquareRoot (maxBound `div` 2 :: Int)) (2 ^ 31 - 1)
+
+-- | Check that 'integerSquareRoot' of 2^62-1 is 2^31-1, not 2^31.
+integerSquareRootSpecialCase1_Word :: Assertion
+integerSquareRootSpecialCase1_Word =
+  assertEqual "integerSquareRoot" (integerSquareRoot (maxBound `div` 4 :: Word)) (2 ^ 31 - 1)
+
+-- | Check that 'integerSquareRoot' of 2^64-1 is 2^32-1, not 2^32.
+integerSquareRootSpecialCase2 :: Assertion
+integerSquareRootSpecialCase2 =
+  assertEqual "integerSquareRoot" (integerSquareRoot (maxBound :: Word)) (2 ^ 32 - 1)
+
+#endif
+
+-- | Check that 'integerSquareRoot'' returns the largest integer @r@ with @r*r <= n@.
+integerSquareRoot'Property :: Integral a => NonNegative a -> Bool
+integerSquareRoot'Property (NonNegative n) = m >=0 && m * m <= n && (m + 1) ^ 2 /= n && m + 1 >= n `div` (m + 1)
+  where
+    m = integerSquareRoot' n
+
+-- | Check that the number 'isSquare' iff its 'integerSquareRoot' is exact.
+isSquareProperty :: Integral a => AnySign a -> Bool
+isSquareProperty (AnySign n) = (n < 0 && not t) || (n /= m * m && not t) || (n == m * m && t)
+  where
+    t = isSquare n
+    m = integerSquareRoot n
+
+-- | Check that the number 'isSquare'' iff its 'integerSquareRoot'' is exact.
+isSquare'Property :: Integral a => NonNegative a -> Bool
+isSquare'Property (NonNegative n) = (n /= m * m && not t) || (n == m * m && t)
+  where
+    t = isSquare' n
+    m = integerSquareRoot' n
+
+-- | Check that 'exactSquareRoot' returns an exact integer square root
+-- and is consistent with 'isSquare'.
+exactSquareRootProperty :: Integral a => AnySign a -> Bool
+exactSquareRootProperty (AnySign n) = case exactSquareRoot n of
+  Nothing -> not (isSquare n)
+  Just m  -> isSquare n && n == m * m
+
+-- | Check that 'isPossibleSquare' is consistent with 'exactSquareRoot'
+-- and that 'isPossibleSquare2' is a refinement of 'isPossibleSquare'.
+isPossibleSquareProperty :: Integral a => NonNegative a -> Bool
+isPossibleSquareProperty (NonNegative n) = t || not t && not t2 && isNothing m
+  where
+    t = isPossibleSquare n
+    t2 = isPossibleSquare2 n
+    m = exactSquareRoot n
+
+-- | Check that 'isPossibleSquare2'' is consistent with 'exactSquareRoot'.
+isPossibleSquare2Property :: Integral a => NonNegative a -> Bool
+isPossibleSquare2Property (NonNegative n) = t || not t && isNothing m
+  where
+    t = isPossibleSquare2 n
+    m = exactSquareRoot n
+
+
+testSuite :: TestTree
+testSuite = testGroup "Squares"
+  [ testGroup "integerSquareRoor"
+    [ testIntegralProperty "generic"          integerSquareRootProperty
+    , testSmallAndQuick    "generic Int"      integerSquareRootProperty_Int
+    , testSmallAndQuick    "generic Word"     integerSquareRootProperty_Word
+
+    , testIntegralProperty "almost square"      integerSquareRootProperty2
+    , testSmallAndQuick    "almost square Int"  integerSquareRootProperty2_Int
+    , testSmallAndQuick    "almost square Word" integerSquareRootProperty2_Word
+
+#if WORD_SIZE_IN_BITS == 64
+    , testCase             "maxBound / 2 :: Int"  integerSquareRootSpecialCase1_Int
+    , testCase             "maxBound / 4 :: Word" integerSquareRootSpecialCase1_Word
+    , testCase             "maxBound :: Word"     integerSquareRootSpecialCase2
+#endif
+    ]
+
+  , testIntegralProperty "integerSquareRoot'" integerSquareRoot'Property
+  , testIntegralProperty "isSquare"           isSquareProperty
+  , testIntegralProperty "isSquare'"          isSquare'Property
+  , testIntegralProperty "exactSquareRoot"    exactSquareRootProperty
+  , testIntegralProperty "isPossibleSquare"   isPossibleSquareProperty
+  , testIntegralProperty "isPossibleSquare2"  isPossibleSquare2Property
+  ]
diff --git a/test-suite/Math/NumberTheory/Primes/CountingTests.hs b/test-suite/Math/NumberTheory/Primes/CountingTests.hs
new file mode 100644
--- /dev/null
+++ b/test-suite/Math/NumberTheory/Primes/CountingTests.hs
@@ -0,0 +1,148 @@
+-- |
+-- Module:      Math.NumberTheory.Primes.CountingTests
+-- Copyright:   (c) 2016 Andrew Lelechenko
+-- Licence:     MIT
+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>
+-- Stability:   Provisional
+--
+-- Tests for Math.NumberTheory.Primes.Counting
+--
+
+{-# OPTIONS_GHC -fno-warn-type-defaults #-}
+
+module Math.NumberTheory.Primes.CountingTests
+  ( testSuite
+  ) where
+
+import Test.Tasty
+import Test.Tasty.HUnit
+
+import Math.NumberTheory.Primes.Counting
+import Math.NumberTheory.Primes.Testing
+import Math.NumberTheory.TestUtils
+
+-- | https://en.wikipedia.org/wiki/Prime-counting_function#Table_of_.CF.80.28x.29.2C_x_.2F_ln_x.2C_and_li.28x.29
+table :: [(Integer, Integer)]
+table =
+  [ (10^1,   4)
+  , (10^2,   25)
+  , (10^3,   168)
+  , (10^4,   1229)
+  , (10^5,   9592)
+  , (10^6,   78498)
+  , (10^7,   664579)
+  , (10^8,   5761455)
+  , (10^9,   50847534)
+  , (10^10,  455052511)
+  , (10^11,  4118054813)
+  , (10^12,  37607912018)
+  , (10^13,  346065536839)
+  -- , (10^14,  3204941750802)
+  -- , (10^15,  29844570422669)
+  -- , (10^16,  279238341033925)
+  -- , (10^17,  2623557157654233)
+  -- , (10^18,  24739954287740860)
+  -- , (10^19,  234057667276344607)
+  -- , (10^20,  2220819602560918840)
+  ]
+
+-- | Check that values of 'primeCount' are non-negative.
+primeCountProperty1 :: Integer -> Bool
+primeCountProperty1 n = n > primeCountMaxArg
+  || n >  0 && primeCount n >= 0
+  || n <= 0 && primeCount n == 0
+
+-- | Check that 'primeCount' is monotonically increasing function.
+primeCountProperty2 :: Positive Integer -> Positive Integer -> Bool
+primeCountProperty2 (Positive n1) (Positive n2)
+  =  n1 > primeCountMaxArg
+  || n2 > primeCountMaxArg
+  || n1 <= n2 && p1 <= p2
+  || n1 >  n2 && p1 >= p2
+  where
+    p1 = primeCount n1
+    p2 = primeCount n2
+
+-- | Check that 'primeCount' is strictly increasing iff an argument is prime.
+primeCountProperty3 :: Positive Integer -> Bool
+primeCountProperty3 (Positive n)
+  =  isPrime n && primeCount (n - 1) + 1 == primeCount n
+  || not (isPrime n) && primeCount (n - 1) == primeCount n
+
+-- | Check tabulated values.
+primeCountSpecialCases :: [Assertion]
+primeCountSpecialCases = map a table
+  where
+  a (n, m) = assertEqual "primeCount" m (primeCount n)
+
+
+-- | Check that values of 'nthPrime' are positive.
+nthPrimeProperty1 :: Positive Integer -> Bool
+nthPrimeProperty1 (Positive n) = n > nthPrimeMaxArg
+  || nthPrime n > 0
+
+-- | Check that 'nthPrime' is monotonically increasing function.
+nthPrimeProperty2 :: Positive Integer -> Positive Integer -> Bool
+nthPrimeProperty2 (Positive n1) (Positive n2)
+  =  n1 > nthPrimeMaxArg
+  || n2 > nthPrimeMaxArg
+  || n1 <= n2 && p1 <= p2
+  || n1 >  n2 && p1 >= p2
+  where
+    p1 = nthPrime n1
+    p2 = nthPrime n2
+
+-- | Check that values of 'nthPrime' are prime.
+nthPrimeProperty3 :: Positive Integer -> Bool
+nthPrimeProperty3 (Positive n) = isPrime $ nthPrime n
+
+-- | Check tabulated values.
+nthPrimeSpecialCases :: [Assertion]
+nthPrimeSpecialCases = map a table
+  where
+  a (n, m) = assertBool "nthPrime" $ n > nthPrime m
+
+
+-- | Check that values of 'approxPrimeCount' are non-negative.
+approxPrimeCountProperty1 :: Integral a => AnySign a -> Bool
+approxPrimeCountProperty1 (AnySign a) = approxPrimeCount a >= 0
+
+-- | Check that 'approxPrimeCount' is consistent with 'approxPrimeCountOverestimateLimit'.
+approxPrimeCountProperty2 :: Integral a => Positive a -> Bool
+approxPrimeCountProperty2 (Positive a) = a >= approxPrimeCountOverestimateLimit
+  || toInteger (approxPrimeCount a) >= primeCount (toInteger a)
+
+
+-- | Check that values of 'nthPrimeApprox' are positive.
+nthPrimeApproxProperty1 :: AnySign Integer -> Bool
+nthPrimeApproxProperty1 (AnySign a) = nthPrimeApprox a > 0
+
+-- | Check that 'nthPrimeApprox' is consistent with 'nthPrimeApproxUnderestimateLimit'.
+nthPrimeApproxProperty2 :: Positive Integer -> Bool
+nthPrimeApproxProperty2 (Positive a) = a >= nthPrimeApproxUnderestimateLimit
+  || toInteger (nthPrimeApprox a) <= nthPrime (toInteger a)
+
+
+testSuite :: TestTree
+testSuite = testGroup "Counting"
+  [ testGroup "primeCount"
+    ( testSmallAndQuick "non-negative"        primeCountProperty1
+    : testSmallAndQuick "monotonic"           primeCountProperty2
+    : testSmallAndQuick "increases on primes" primeCountProperty3
+    : zipWith (\i a -> testCase ("special case " ++ show i) a) [1..] primeCountSpecialCases
+    )
+  , testGroup "nthPrime"
+    ( testSmallAndQuick "positive"  nthPrimeProperty1
+    : testSmallAndQuick "monotonic" nthPrimeProperty2
+    : testSmallAndQuick "is prime"  nthPrimeProperty3
+    : zipWith (\i a -> testCase ("special case " ++ show i) a) [1..] nthPrimeSpecialCases
+    )
+  , testGroup "approxPrimeCount"
+    [ testIntegralProperty "non-negative"             approxPrimeCountProperty1
+    , testIntegralProperty "overestimates primeCount" approxPrimeCountProperty2
+    ]
+  , testGroup "nthPrimeApprox"
+    [ testSmallAndQuick "positive"                nthPrimeApproxProperty1
+    , testSmallAndQuick "underestimates nthPrime" nthPrimeApproxProperty2
+    ]
+  ]
diff --git a/test-suite/Math/NumberTheory/Primes/HeapTests.hs b/test-suite/Math/NumberTheory/Primes/HeapTests.hs
new file mode 100644
--- /dev/null
+++ b/test-suite/Math/NumberTheory/Primes/HeapTests.hs
@@ -0,0 +1,67 @@
+-- |
+-- Module:      Math.NumberTheory.Primes.HeapTests
+-- Copyright:   (c) 2016 Andrew Lelechenko
+-- Licence:     MIT
+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>
+-- Stability:   Provisional
+--
+-- Tests for Math.NumberTheory.Primes.Heap
+--
+
+{-# LANGUAGE CPP #-}
+
+{-# OPTIONS_GHC -fno-warn-type-defaults #-}
+
+module Math.NumberTheory.Primes.HeapTests
+  ( testSuite
+  ) where
+
+import Prelude hiding (words)
+
+import Test.Tasty
+import Test.Tasty.HUnit
+
+#if MIN_VERSION_base(4,8,0)
+#else
+import Data.Word
+#endif
+
+import Math.NumberTheory.Primes.Heap
+import Math.NumberTheory.Primes.Testing
+import Math.NumberTheory.TestUtils
+
+-- | Check that 'primes' over different integral types matches with 'isPrime'.
+primesProperty1 :: Assertion
+primesProperty1 = do
+  assertEqual "ints  == integers" (trim ints)  (trim integers)
+  assertEqual "words == integers" (trim words) (trim integers)
+  assertEqual "naive == integers" (trim naive) (trim integers)
+  where
+    trim :: Integral a => [a] -> [Integer]
+    trim = map toInteger . take 100000
+
+    ints     = primes :: [Int]
+    words    = primes :: [Word]
+    integers = primes :: [Integer]
+    naive    = filter isPrime [1..] :: [Integer]
+
+-- | Check that 'sieveFrom' over different integral types matches with 'isPrime'.
+sieveFromProperty1 :: NonNegative Integer -> Bool
+sieveFromProperty1 (NonNegative lowBound)
+  =  trim ints  == trim integers
+  && trim words == trim integers
+  && trim naive == trim integers
+  where
+    trim :: Integral a => [a] -> [Integer]
+    trim = map toInteger . take 1000
+
+    ints     = sieveFrom (fromInteger lowBound) :: [Int]
+    words    = sieveFrom (fromInteger lowBound) :: [Word]
+    integers = sieveFrom lowBound               :: [Integer]
+    naive    = filter isPrime [lowBound..]      :: [Integer]
+
+testSuite :: TestTree
+testSuite = testGroup "Heap"
+  [ testCase          "primes"    primesProperty1
+  , testSmallAndQuick "sieveFrom" sieveFromProperty1
+  ]
diff --git a/test-suite/Math/NumberTheory/Primes/SieveTests.hs b/test-suite/Math/NumberTheory/Primes/SieveTests.hs
new file mode 100644
--- /dev/null
+++ b/test-suite/Math/NumberTheory/Primes/SieveTests.hs
@@ -0,0 +1,69 @@
+-- |
+-- Module:      Math.NumberTheory.Primes.SieveTests
+-- Copyright:   (c) 2016 Andrew Lelechenko
+-- Licence:     MIT
+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>
+-- Stability:   Provisional
+--
+-- Tests for Math.NumberTheory.Primes.Sieve
+--
+
+{-# LANGUAGE CPP #-}
+
+{-# OPTIONS_GHC -fno-warn-type-defaults #-}
+
+module Math.NumberTheory.Primes.SieveTests
+  ( testSuite
+  ) where
+
+import Prelude hiding (words)
+
+import Test.Tasty
+import Test.Tasty.HUnit
+
+import Math.NumberTheory.Primes.Sieve
+import qualified Math.NumberTheory.Primes.Heap as H
+import Math.NumberTheory.TestUtils
+
+-- | Check that both 'primes' produce the same.
+primesProperty1 :: Assertion
+primesProperty1 = do
+  assertEqual "Sieve == Heap" (trim primes) (trim H.primes)
+  where
+    trim = take 100000
+
+-- | Check that both 'sieveFrom' produce the same.
+sieveFromProperty1 :: AnySign Integer -> Bool
+sieveFromProperty1 (AnySign lowBound)
+  = trim (sieveFrom lowBound) == trim (H.sieveFrom lowBound)
+  where
+    trim = take 1000
+
+-- | Check that 'primeList' from 'primeSieve' matches truncated 'primes'.
+primeSieveProperty1 :: AnySign Integer -> Bool
+primeSieveProperty1 (AnySign highBound)
+  = primeList (primeSieve highBound) == takeWhile (<= (highBound `max` 7)) primes
+
+-- | Check that 'primeList' from 'psieveList' matches 'primes'.
+psieveListProperty1 :: Assertion
+psieveListProperty1 = do
+  assertEqual "primes == primeList . psieveList" (trim primes) (trim $ concatMap primeList psieveList)
+  where
+    trim = take 100000
+
+-- | Check that 'primeList' from 'psieveFrom' matches 'sieveFrom'.
+psieveFromProperty1 :: AnySign Integer -> Bool
+psieveFromProperty1 (AnySign lowBound)
+  = trim (sieveFrom lowBound) == trim (filter (>= lowBound) (concatMap primeList $ psieveFrom lowBound))
+  where
+    trim = take 1000
+
+
+testSuite :: TestTree
+testSuite = testGroup "Sieve"
+  [ testCase          "primes"     primesProperty1
+  , testSmallAndQuick "sieveFrom"  sieveFromProperty1
+  , testSmallAndQuick "primeSieve" primeSieveProperty1
+  , testCase          "psieveList" psieveListProperty1
+  , testSmallAndQuick "psieveFrom" psieveFromProperty1
+  ]
diff --git a/test-suite/Math/NumberTheory/PrimesTests.hs b/test-suite/Math/NumberTheory/PrimesTests.hs
new file mode 100644
--- /dev/null
+++ b/test-suite/Math/NumberTheory/PrimesTests.hs
@@ -0,0 +1,40 @@
+-- |
+-- Module:      Math.NumberTheory.PrimesTests
+-- Copyright:   (c) 2016 Andrew Lelechenko
+-- Licence:     MIT
+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>
+-- Stability:   Provisional
+--
+-- Tests for Math.NumberTheory.Primes
+--
+
+{-# OPTIONS_GHC -fno-warn-type-defaults #-}
+
+module Math.NumberTheory.PrimesTests
+  ( testSuite
+  ) where
+
+import Test.Tasty
+import Test.Tasty.HUnit
+
+import Math.NumberTheory.Primes
+import Math.NumberTheory.TestUtils
+
+primesSumWonk :: Int -> Int
+primesSumWonk upto = sum . takeWhile (< upto) . map fromInteger . primeList $ primeSieve (toInteger upto)
+
+primesSum :: Int -> Int
+primesSum upto = sum . takeWhile (< upto) . map fromInteger $ primes
+
+primesSumProperty :: NonNegative Int -> Bool
+primesSumProperty (NonNegative n) = primesSumWonk n == primesSum n
+
+
+sieveFactorSpecialCase1 :: Assertion
+sieveFactorSpecialCase1 = assertEqual "sieveFactor" [(29, 1), (73, 1)] $ sieveFactor (factorSieve 2048) (29*73)
+
+testSuite :: TestTree
+testSuite = testGroup "Primes"
+  [ testSmallAndQuick "primesSum"   primesSumProperty
+  , testCase          "sieveFactor" sieveFactorSpecialCase1
+  ]
diff --git a/test-suite/Math/NumberTheory/TestUtils.hs b/test-suite/Math/NumberTheory/TestUtils.hs
new file mode 100644
--- /dev/null
+++ b/test-suite/Math/NumberTheory/TestUtils.hs
@@ -0,0 +1,226 @@
+-- |
+-- Module:      Math.NumberTheory.TestUtils
+-- Copyright:   (c) 2016 Andrew Lelechenko
+-- Licence:     MIT
+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>
+-- Stability:   Provisional
+-- Portability: Non-portable (GHC extensions)
+--
+-- Utils to test Math.NumberTheory
+--
+
+{-# LANGUAGE ConstraintKinds            #-}
+{-# LANGUAGE CPP                        #-}
+{-# LANGUAGE DataKinds                  #-}
+{-# LANGUAGE DeriveFoldable             #-}
+{-# LANGUAGE DeriveFunctor              #-}
+{-# LANGUAGE DeriveTraversable          #-}
+{-# LANGUAGE FlexibleContexts           #-}
+{-# LANGUAGE FlexibleInstances          #-}
+{-# LANGUAGE GeneralizedNewtypeDeriving #-}
+{-# LANGUAGE KindSignatures             #-}
+{-# LANGUAGE MultiParamTypeClasses      #-}
+{-# LANGUAGE RankNTypes                 #-}
+{-# LANGUAGE ScopedTypeVariables        #-}
+{-# LANGUAGE TypeFamilies               #-}
+{-# LANGUAGE TypeOperators              #-}
+{-# LANGUAGE UndecidableInstances       #-}
+
+#if __GLASGOW_HASKELL__ >= 800
+{-# LANGUAGE UndecidableSuperClasses    #-}
+
+{-# OPTIONS_GHC -fconstraint-solver-iterations=0 #-}
+#endif
+
+{-# OPTIONS_GHC -fno-warn-orphans #-}
+{-# OPTIONS_GHC -fno-warn-type-defaults #-}
+
+module Math.NumberTheory.TestUtils
+  ( module Math.NumberTheory.TestUtils
+  , module Test.SmallCheck.Series
+  , Large(..)
+  ) where
+
+import Test.Tasty
+import Test.Tasty.SmallCheck as SC
+import Test.Tasty.QuickCheck as QC hiding (Positive, NonNegative, generate, getNonNegative)
+
+import Test.SmallCheck.Series (Positive(..), NonNegative(..), Serial(..), Series, generate)
+
+import Control.Applicative
+import Data.Bits
+#if MIN_VERSION_base(4,8,0)
+#else
+import Data.Foldable (Foldable)
+import Data.Traversable (Traversable)
+import Data.Word
+#endif
+import GHC.Exts
+
+import Math.NumberTheory.Primes
+
+newtype AnySign a = AnySign { getAnySign :: a }
+  deriving (Eq, Ord, Read, Show, Num, Enum, Bounded, Integral, Real, Functor, Foldable, Traversable, Arbitrary)
+
+instance (Monad m, Serial m a) => Serial m (AnySign a) where
+  series = AnySign <$> series
+
+instance (Num a, Ord a, Arbitrary a) => Arbitrary (Positive a) where
+  arbitrary = Positive <$> (arbitrary `suchThat` (> 0))
+  shrink (Positive x) = Positive <$> filter (> 0) (shrink x)
+
+instance (Num a, Ord a, Arbitrary a) => Arbitrary (NonNegative a) where
+  arbitrary = NonNegative <$> (arbitrary `suchThat` (>= 0))
+  shrink (NonNegative x) = NonNegative <$> filter (>= 0) (shrink x)
+
+instance (Num a, Bounded a) => Bounded (Positive a) where
+  minBound = Positive 1
+  maxBound = Positive (maxBound :: a)
+
+instance (Num a, Bounded a) => Bounded (NonNegative a) where
+  minBound = NonNegative 0
+  maxBound = NonNegative (maxBound :: a)
+
+newtype Huge a = Huge { getHuge :: a }
+  deriving (Eq, Ord, Enum, Bounded, Show, Num, Real, Integral)
+
+instance (Num a, Arbitrary a) => Arbitrary (Huge a) where
+  arbitrary = do
+    Positive l <- arbitrary
+    ds <- vector l
+    return $ Huge $ foldl1 (\acc n -> acc * 2^63 + n) ds
+
+newtype Power a = Power { getPower :: a }
+  deriving (Eq, Ord, Enum, Bounded, Show, Num, Real, Integral)
+
+instance (Monad m, Num a, Ord a, Serial m a) => Serial m (Power a) where
+  series = Power <$> series `suchThatSerial` (> 0)
+
+instance (Num a, Ord a, Integral a, Arbitrary a) => Arbitrary (Power a) where
+  arbitrary = Power <$> (getSmall <$> arbitrary) `suchThat` (> 0)
+  shrink (Power x) = Power <$> filter (> 0) (shrink x)
+
+newtype Prime = Prime { getPrime :: Integer }
+  deriving (Eq, Ord, Show)
+
+instance Arbitrary Prime where
+  arbitrary = Prime <$> arbitrary `suchThat` (\p -> p > 0 && isPrime p)
+
+instance Monad m => Serial m Prime where
+  series = Prime <$> series `suchThatSerial` (\p -> p > 0 && isPrime p)
+
+instance Monad m => Serial m Word where
+  series =
+    generate (\d -> if d >= 0 then pure 0 else empty) <|> nats
+    where
+      nats = generate $ \d -> if d > 0 then [1 .. fromInteger (toInteger d)] else empty
+
+suchThatSerial :: Series m a -> (a -> Bool) -> Series m a
+suchThatSerial s p = s >>= \x -> if p x then pure x else empty
+
+
+-- https://www.cs.ox.ac.uk/projects/utgp/school/andres.pdf, p. 21
+-- :k Compose = (k1 -> Constraint) -> (k2 -> k1) -> (k2 -> Constraint)
+class    (f (g x)) => (f `Compose` g) x
+instance (f (g x)) => (f `Compose` g) x
+
+type family ConcatMap (w :: * -> Constraint) (cs :: [*]) :: Constraint
+#if __GLASGOW_HASKELL__ >= 708
+  where
+    ConcatMap w '[] = ()
+    ConcatMap w (c ': cs) = (w c, ConcatMap w cs)
+#else
+type instance ConcatMap w '[] = ()
+type instance ConcatMap w (c ': cs) = (w c, ConcatMap w cs)
+#endif
+
+type family Matrix (as :: [* -> Constraint]) (w :: * -> *) (bs :: [*]) :: Constraint
+#if __GLASGOW_HASKELL__ >= 708
+  where
+    Matrix '[] w bs = ()
+    Matrix (a ': as) w bs = (ConcatMap (a `Compose` w) bs, Matrix as w bs)
+#else
+type instance Matrix '[] w bs = ()
+type instance Matrix (a ': as) w bs = (ConcatMap (a `Compose` w) bs, Matrix as w bs)
+#endif
+
+type TestableIntegral wrapper =
+  ( Matrix '[Arbitrary, Show, Serial IO] wrapper '[Int, Word, Integer]
+  , Matrix '[Bounded, Integral] wrapper '[Int, Word]
+  , Num (wrapper Integer)
+  )
+
+
+testIntegralProperty
+  :: forall wrapper bool. (TestableIntegral wrapper, SC.Testable IO bool, QC.Testable bool)
+  => String -> (forall a. (Integral a, Bits a) => wrapper a -> bool) -> TestTree
+testIntegralProperty name f = testGroup name
+  [ SC.testProperty "smallcheck Int"     (f :: wrapper Int     -> bool)
+  , SC.testProperty "smallcheck Word"    (f :: wrapper Word    -> bool)
+  , SC.testProperty "smallcheck Integer" (f :: wrapper Integer -> bool)
+  , QC.testProperty "quickcheck Int"     (f :: wrapper Int     -> bool)
+  , QC.testProperty "quickcheck Word"    (f :: wrapper Word    -> bool)
+  , QC.testProperty "quickcheck Integer" (f :: wrapper Integer -> bool)
+  , QC.testProperty "quickcheck Large Int"     ((f :: wrapper Int     -> bool) . getLarge)
+  , QC.testProperty "quickcheck Large Word"    ((f :: wrapper Word    -> bool) . getLarge)
+  , QC.testProperty "quickcheck Huge  Integer" ((f :: wrapper Integer -> bool) . getHuge)
+  ]
+
+testSameIntegralProperty
+  :: forall wrapper1 wrapper2 bool. (TestableIntegral wrapper1, TestableIntegral wrapper2, SC.Testable IO bool, QC.Testable bool)
+  => String -> (forall a. (Integral a, Bits a) => wrapper1 a -> wrapper2 a -> bool) -> TestTree
+testSameIntegralProperty name f = testGroup name
+  [ SC.testProperty "smallcheck Int"     (f :: wrapper1 Int     -> wrapper2 Int     -> bool)
+  , SC.testProperty "smallcheck Word"    (f :: wrapper1 Word    -> wrapper2 Word    -> bool)
+  , SC.testProperty "smallcheck Integer" (f :: wrapper1 Integer -> wrapper2 Integer -> bool)
+  , QC.testProperty "quickcheck Int"     (f :: wrapper1 Int     -> wrapper2 Int     -> bool)
+  , QC.testProperty "quickcheck Word"    (f :: wrapper1 Word    -> wrapper2 Word    -> bool)
+  , QC.testProperty "quickcheck Integer" (f :: wrapper1 Integer -> wrapper2 Integer -> bool)
+  , QC.testProperty "quickcheck Large Int"     (\(Large a) (Large b) -> (f :: wrapper1 Int     -> wrapper2 Int     -> bool) a b)
+  , QC.testProperty "quickcheck Large Word"    (\(Large a) (Large b) -> (f :: wrapper1 Word    -> wrapper2 Word    -> bool) a b)
+  , QC.testProperty "quickcheck Huge  Integer" (\(Huge  a) (Huge  b) -> (f :: wrapper1 Integer -> wrapper2 Integer -> bool) a b)
+  ]
+
+testIntegral2Property
+  :: forall wrapper1 wrapper2 bool. (TestableIntegral wrapper1, TestableIntegral wrapper2, SC.Testable IO bool, QC.Testable bool)
+  => String -> (forall a1 a2. (Integral a1, Integral a2, Bits a1, Bits a2) => wrapper1 a1 -> wrapper2 a2 -> bool) -> TestTree
+testIntegral2Property name f = testGroup name
+  [ SC.testProperty "smallcheck Int Int"         (f :: wrapper1 Int     -> wrapper2 Int     -> bool)
+  , SC.testProperty "smallcheck Int Word"        (f :: wrapper1 Int     -> wrapper2 Word    -> bool)
+  , SC.testProperty "smallcheck Int Integer"     (f :: wrapper1 Int     -> wrapper2 Integer -> bool)
+  , SC.testProperty "smallcheck Word Int"        (f :: wrapper1 Word    -> wrapper2 Int     -> bool)
+  , SC.testProperty "smallcheck Word Word"       (f :: wrapper1 Word    -> wrapper2 Word    -> bool)
+  , SC.testProperty "smallcheck Word Integer"    (f :: wrapper1 Word    -> wrapper2 Integer -> bool)
+  , SC.testProperty "smallcheck Integer Int"     (f :: wrapper1 Integer -> wrapper2 Int     -> bool)
+  , SC.testProperty "smallcheck Integer Word"    (f :: wrapper1 Integer -> wrapper2 Word    -> bool)
+  , SC.testProperty "smallcheck Integer Integer" (f :: wrapper1 Integer -> wrapper2 Integer -> bool)
+
+  , QC.testProperty "quickcheck Int Int"         (f :: wrapper1 Int     -> wrapper2 Int     -> bool)
+  , QC.testProperty "quickcheck Int Word"        (f :: wrapper1 Int     -> wrapper2 Word    -> bool)
+  , QC.testProperty "quickcheck Int Integer"     (f :: wrapper1 Int     -> wrapper2 Integer -> bool)
+  , QC.testProperty "quickcheck Word Int"        (f :: wrapper1 Word    -> wrapper2 Int     -> bool)
+  , QC.testProperty "quickcheck Word Word"       (f :: wrapper1 Word    -> wrapper2 Word    -> bool)
+  , QC.testProperty "quickcheck Word Integer"    (f :: wrapper1 Word    -> wrapper2 Integer -> bool)
+  , QC.testProperty "quickcheck Integer Int"     (f :: wrapper1 Integer -> wrapper2 Int     -> bool)
+  , QC.testProperty "quickcheck Integer Word"    (f :: wrapper1 Integer -> wrapper2 Word    -> bool)
+  , QC.testProperty "quickcheck Integer Integer" (f :: wrapper1 Integer -> wrapper2 Integer -> bool)
+
+  , QC.testProperty "quickcheck Large Int Int"         ((f :: wrapper1 Int     -> wrapper2 Int     -> bool) . getLarge)
+  , QC.testProperty "quickcheck Large Int Word"        ((f :: wrapper1 Int     -> wrapper2 Word    -> bool) . getLarge)
+  , QC.testProperty "quickcheck Large Int Integer"     ((f :: wrapper1 Int     -> wrapper2 Integer -> bool) . getLarge)
+  , QC.testProperty "quickcheck Large Word Int"        ((f :: wrapper1 Word    -> wrapper2 Int     -> bool) . getLarge)
+  , QC.testProperty "quickcheck Large Word Word"       ((f :: wrapper1 Word    -> wrapper2 Word    -> bool) . getLarge)
+  , QC.testProperty "quickcheck Large Word Integer"    ((f :: wrapper1 Word    -> wrapper2 Integer -> bool) . getLarge)
+  , QC.testProperty "quickcheck Huge  Integer Int"     ((f :: wrapper1 Integer -> wrapper2 Int     -> bool) . getHuge)
+  , QC.testProperty "quickcheck Huge  Integer Word"    ((f :: wrapper1 Integer -> wrapper2 Word    -> bool) . getHuge)
+  , QC.testProperty "quickcheck Huge  Integer Integer" ((f :: wrapper1 Integer -> wrapper2 Integer -> bool) . getHuge)
+  ]
+
+testSmallAndQuick
+  :: SC.Testable IO a
+  => QC.Testable a
+  => String -> a -> TestTree
+testSmallAndQuick name f = testGroup name
+  [ SC.testProperty "smallcheck" f
+  , QC.testProperty "quickcheck" f
+  ]
diff --git a/test-suite/Spec.hs b/test-suite/Spec.hs
deleted file mode 100644
--- a/test-suite/Spec.hs
+++ /dev/null
@@ -1,1 +0,0 @@
-{-# OPTIONS_GHC -F -pgmF hspec-discover #-}
diff --git a/test-suite/Test.hs b/test-suite/Test.hs
new file mode 100644
--- /dev/null
+++ b/test-suite/Test.hs
@@ -0,0 +1,66 @@
+import Test.Tasty
+
+import qualified Math.NumberTheory.GCDTests as GCD
+import qualified Math.NumberTheory.GCD.LowLevelTests as GCDLowLevel
+
+import qualified Math.NumberTheory.LogarithmsTests as Logarithms
+
+import qualified Math.NumberTheory.LucasTests as Lucas
+
+import qualified Math.NumberTheory.ModuliTests as Moduli
+
+import qualified Math.NumberTheory.MoebiusInversionTests as MoebiusInversion
+import qualified Math.NumberTheory.MoebiusInversion.IntTests as MoebiusInversionInt
+
+import qualified Math.NumberTheory.Powers.CubesTests as Cubes
+import qualified Math.NumberTheory.Powers.FourthTests as Fourth
+import qualified Math.NumberTheory.Powers.GeneralTests as General
+import qualified Math.NumberTheory.Powers.IntegerTests as Integer
+import qualified Math.NumberTheory.Powers.SquaresTests as Squares
+
+import qualified Math.NumberTheory.PrimesTests as Primes
+import qualified Math.NumberTheory.Primes.CountingTests as Counting
+import qualified Math.NumberTheory.Primes.HeapTests as Heap
+import qualified Math.NumberTheory.Primes.SieveTests as Sieve
+
+import qualified Math.NumberTheory.GaussianIntegersTests as Gaussian
+
+main :: IO ()
+main = defaultMain tests
+
+tests :: TestTree
+tests = testGroup "All"
+  [ testGroup "Powers"
+    [ Cubes.testSuite
+    , Fourth.testSuite
+    , General.testSuite
+    , Integer.testSuite
+    , Squares.testSuite
+    ]
+  , testGroup "GCD"
+    [ GCD.testSuite
+    , GCDLowLevel.testSuite
+    ]
+  , testGroup "Logarithms"
+    [ Logarithms.testSuite
+    ]
+  , testGroup "Lucas"
+    [ Lucas.testSuite
+    ]
+  , testGroup "Moduli"
+    [ Moduli.testSuite
+    ]
+  , testGroup "MoebiusInversion"
+    [ MoebiusInversion.testSuite
+    , MoebiusInversionInt.testSuite
+    ]
+  , testGroup "Primes"
+    [ Primes.testSuite
+    , Counting.testSuite
+    , Heap.testSuite
+    , Sieve.testSuite
+    ]
+  , testGroup "Gaussian"
+    [ Gaussian.testSuite
+    ]
+  ]