diff --git a/CHANGELOG.md b/CHANGELOG.md
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,6 +1,11 @@
 Changes
 =======
 
+Version 0.2.1.1
+----------------
+
+Switch benchmarks from gauge to tasty-bench.
+
 Version 0.2.1.0
 ----------------
 
diff --git a/Math/ExpPairs.hs b/Math/ExpPairs.hs
--- a/Math/ExpPairs.hs
+++ b/Math/ExpPairs.hs
@@ -44,8 +44,8 @@
 import Data.Ord      (comparing)
 import Data.List     (minimumBy)
 import Data.Ratio
-import Data.Text.Prettyprint.Doc hiding ((<>))
-import qualified Data.Text.Prettyprint.Doc as PP
+import Prettyprinter hiding ((<>))
+import qualified Prettyprinter as PP
 import Text.Printf
 
 import Math.ExpPairs.LinearForm
diff --git a/Math/ExpPairs/Kratzel.hs b/Math/ExpPairs/Kratzel.hs
--- a/Math/ExpPairs/Kratzel.hs
+++ b/Math/ExpPairs/Kratzel.hs
@@ -38,11 +38,13 @@
 
 import Control.Arrow hiding ((<+>))
 import Data.Function
+import Data.List.NonEmpty (NonEmpty (..))
+import qualified Data.List.NonEmpty as NE
 import Data.Maybe
 import Data.Ratio
 import Data.Ord   (comparing)
 import Data.List  (minimumBy, sort, inits, tails)
-import Data.Text.Prettyprint.Doc
+import Prettyprinter
 
 import qualified Data.Map as M
 import qualified Data.Set as S
@@ -99,8 +101,12 @@
       d = a `gcd` b
 tauab a b = tauab' a' b'
   where
-    [a', b'] = sort [a, b]
+    (a', b') = sort2 (a, b)
 
+sort2 :: Ord a => (a, a) -> (a, a)
+sort2 (a, b)
+  | a <= b    = (a, b)
+  | otherwise = (b, a)
 
 -- |Special type to specify the theorem of Krätzel1988,
 -- which provided the best estimate of Θ(a, b, c)
@@ -170,8 +176,15 @@
       d = a `gcd` b `gcd` c
 tauabc a b c = tauabc' a' b' c'
   where
-    [a', b', c'] = sort [a, b, c]
+    (a', b', c') = sort3 (a, b, c)
 
+sort3 :: Ord a => (a, a, a) -> (a, a, a)
+sort3 (a, b', c')
+  | a <= b    = (a, b, c)
+  | a <= c    = (b, a, c)
+  | otherwise = (b, c, a)
+  where
+    (b, c) = sort2 (b', c')
 
 -- |Special type to specify the theorem of Krätzel1988,
 -- which provided the best estimate of Θ(a, b, c, d)
@@ -290,8 +303,16 @@
       d = a1 `gcd` a2 `gcd` a3 `gcd` a4
 tauabcd a1 a2 a3 a4 = tauabcd' a1' a2' a3' a4'
   where
-    [a1', a2', a3', a4'] = sort [a1, a2, a3, a4]
+    (a1', a2', a3', a4') = sort4 (a1, a2, a3, a4)
 
+sort4 :: Ord a => (a, a, a, a) -> (a, a, a, a)
+sort4 (a, b', c', d')
+  | a <= b    = (a, b, c, d)
+  | a <= c    = (b, a, c, d)
+  | a <= d    = (b, c, a, d)
+  | otherwise = (b, c, d, a)
+  where
+    (b, c, d) = sort3 (b', c', d')
 
 -- |Special type to specify the theorem of Krätzel1988,
 -- which provided the best estimate of Θ(a1, a2...)
@@ -350,22 +371,23 @@
     go [_] = Node NoTheorem (simulateOptimize 0)
     go [a, b]    = (\(t, o) -> Node (Ab t) o) $ tauab a b
     go [a, b, c] = (\(t, o) -> Node (Abc t) o) $ tauabc a b c
-    go as@[a,b,c,d] = (\(t, o) -> Node (Abcd t) o) (tauabcd a b c d) `min` go608 as
+    go [a,b,c,d] = (\(t, o) -> Node (Abcd t) o) (tauabcd a b c d) `min` go608 (a :| [b,c,d])
     go as@(a:_)
       | all (== a) as
       = Node Ivic $ simulateOptimize $ reverseMOnS 1e-6 (fromIntegral $ length as) / fi a
-    go as = go608 as
+    go (a : as) = go608 (a :| as)
 
+    go608 :: NonEmpty Integer -> TauAResult
     go608 as = minimum $ mapMaybe f [1 .. length as - 1]
       where
-        f q = if (alphaV `max` betaV) < 1 / fi (last as)
+        f q = if (alphaV `max` betaV) < 1 / fi (NE.last as)
           then Just $ Combination alpha beta ret
           else Nothing
           where
-            alpha = (M.!) cache $ take q as
+            alpha = (M.!) cache $ NE.take q as
             alphaV = extractValue alpha
-            beta  = (M.!) cache $ drop q as
+            beta  = (M.!) cache $ NE.drop q as
             betaV = extractValue beta
-            a0 = fi $ head as
-            aq = fi $ as !! q
+            a0 = fi $ NE.head as
+            aq = fi $ as NE.!! q
             ret = (1 - a0 * aq * alphaV * betaV) / (a0 + aq - a0 * aq * (alphaV + betaV))
diff --git a/Math/ExpPairs/LinearForm.hs b/Math/ExpPairs/LinearForm.hs
--- a/Math/ExpPairs/LinearForm.hs
+++ b/Math/ExpPairs/LinearForm.hs
@@ -29,7 +29,7 @@
 import Data.Maybe     (mapMaybe)
 import Data.Ratio     (numerator, denominator)
 import GHC.Generics   (Generic (..))
-import Data.Text.Prettyprint.Doc
+import Prettyprinter
 
 import Math.ExpPairs.RatioInf
 
diff --git a/Math/ExpPairs/Matrix3.hs b/Math/ExpPairs/Matrix3.hs
--- a/Math/ExpPairs/Matrix3.hs
+++ b/Math/ExpPairs/Matrix3.hs
@@ -30,7 +30,7 @@
 import Data.Foldable  (Foldable (..), toList)
 import Data.List      (transpose)
 import GHC.Generics   (Generic (..))
-import Data.Text.Prettyprint.Doc
+import Prettyprinter
 
 -- | Matrix of order 3.
 data Matrix3 t = Matrix3 {
diff --git a/Math/ExpPairs/Pair.hs b/Math/ExpPairs/Pair.hs
--- a/Math/ExpPairs/Pair.hs
+++ b/Math/ExpPairs/Pair.hs
@@ -28,7 +28,7 @@
 import Data.Maybe
 import Data.Ratio
 import GHC.Generics (Generic (..))
-import Data.Text.Prettyprint.Doc
+import Prettyprinter
 
 -- |Vertices of the triangle of initial exponent pairs.
 data Triangle
diff --git a/Math/ExpPairs/PrettyProcess.hs b/Math/ExpPairs/PrettyProcess.hs
--- a/Math/ExpPairs/PrettyProcess.hs
+++ b/Math/ExpPairs/PrettyProcess.hs
@@ -15,10 +15,10 @@
     uglify,
     PrettyProcess) where
 
-import Data.List                (minimumBy, inits, tails)
+import Data.List                (minimumBy, inits, tails, uncons)
 import Data.Monoid              (mempty)
 import Data.Ord                 (comparing)
-import Data.Text.Prettyprint.Doc
+import Prettyprinter
 
 import qualified Data.Map as M
 import qualified Data.Set as S
@@ -84,7 +84,8 @@
 asRepeat xs = pair where
   l = length xs
   candidates = [ (take d xs, l `div` d) | d <- divisors l ]
-  pair = head $ filter (\(ys, n) -> concat (replicate n ys) == xs) candidates
+  pair = maybe (xs, 1) fst $ uncons $
+    filter (\(ys, n) -> concat (replicate n ys) == xs) candidates
 
 -- | Find the most compact representation of the sequence of processes.
 prettify :: [Process] -> PrettyProcess
@@ -102,15 +103,18 @@
       []   -> annotateWithWidth (Simply [])
       [A]  -> annotateWithWidth (Simply [A])
       [BA] -> annotateWithWidth (Simply [BA])
-      xs   -> minimumBy (comparing ppwlWidth) yss where
+      xs@(xsHd : xsTl) -> minimumBy (comparing ppwlWidth) yss where
         xs'' = case asRepeat xs of
           (_, 1)   -> annotateWithWidth (Simply xs)
           (xs', n) -> annotateWithWidth (Repeat (ppwlProcess $ (M.!) cache xs') n)
 
+        yss :: [PrettyProcessWithWidth]
         yss = xs'' : map f bcs
 
-        bcs = takeWhile (not . null . snd) $ iterate bcf ([head xs], tail xs)
+        bcs :: [([Process], [Process])]
+        bcs = takeWhile (not . null . snd) $ iterate bcf ([xsHd], xsTl)
 
+        bcf :: ([a], [a]) -> ([a], [a])
         bcf (_, [])    = error "prettifyP: unexpected second argument of bcf"
         bcf (zs, y:ys) = (zs++[y], ys)
 
diff --git a/Math/ExpPairs/Process.hs b/Math/ExpPairs/Process.hs
--- a/Math/ExpPairs/Process.hs
+++ b/Math/ExpPairs/Process.hs
@@ -25,7 +25,7 @@
   ) where
 
 import GHC.Generics             (Generic)
-import Data.Text.Prettyprint.Doc hiding ((<>))
+import Prettyprinter hiding ((<>))
 
 import Math.ExpPairs.ProcessMatrix
 import Math.ExpPairs.PrettyProcess
diff --git a/Math/ExpPairs/ProcessMatrix.hs b/Math/ExpPairs/ProcessMatrix.hs
--- a/Math/ExpPairs/ProcessMatrix.hs
+++ b/Math/ExpPairs/ProcessMatrix.hs
@@ -25,7 +25,7 @@
   ) where
 
 import GHC.Generics          (Generic (..))
-import Data.Text.Prettyprint.Doc
+import Prettyprinter
 
 import Math.ExpPairs.Matrix3
 
diff --git a/Math/ExpPairs/RatioInf.hs b/Math/ExpPairs/RatioInf.hs
--- a/Math/ExpPairs/RatioInf.hs
+++ b/Math/ExpPairs/RatioInf.hs
@@ -15,7 +15,7 @@
   ) where
 
 import Data.Ratio (Ratio, numerator, denominator)
-import Data.Text.Prettyprint.Doc
+import Prettyprinter
 
 -- | Extend 'Ratio' @t@ with \( \pm \infty  \) positive and negative
 -- infinities.
@@ -69,24 +69,24 @@
   InfMinus * Finite a = case signum a of
     1  -> InfMinus
     -1 -> InfPlus
-    _  -> error "Cannot multiply infinity by zero"
+    _  -> error "Cannot multiply -infinity by zero"
 
   InfPlus * InfMinus = InfMinus
   InfPlus * InfPlus  = InfPlus
   InfPlus * Finite a = case signum a of
     1  -> InfPlus
     -1 -> InfMinus
-    _  -> error "Cannot multiply infinity by zero"
+    _  -> error "Cannot multiply +infinity by zero"
 
   Finite a * InfMinus = case signum a of
     1  -> InfMinus
     -1 -> InfPlus
-    _  -> error "Cannot multiply infinity by zero"
+    _  -> error "Cannot multiply -infinity by zero"
 
   Finite a * InfPlus = case signum a of
     1  -> InfPlus
     -1 -> InfMinus
-    _  -> error "Cannot multiply infinity by zero"
+    _  -> error "Cannot multiply +infinity by zero"
 
   Finite a * Finite b = Finite (a * b)
 
diff --git a/README.md b/README.md
new file mode 100644
--- /dev/null
+++ b/README.md
@@ -0,0 +1,9 @@
+# exp-pairs [![Hackage](http://img.shields.io/hackage/v/exp-pairs.svg)](https://hackage.haskell.org/package/exp-pairs) [![Stackage LTS](http://stackage.org/package/exp-pairs/badge/lts)](http://stackage.org/lts/package/exp-pairs) [![Stackage Nightly](http://stackage.org/package/exp-pairs/badge/nightly)](http://stackage.org/nightly/package/exp-pairs)
+
+Package implements an algorithm to minimize the maximum of a list of rational objective functions over the set of exponent pairs. See full description in
+> A. V. Lelechenko, Linear programming over exponent pairs.
+> Acta Univ. Sapientiae, Inform. 5, No. 2, 271-287 (2013).
+> http://www.acta.sapientia.ro/acta-info/C5-2/info52-7.pdf
+
+A set of useful applications is also included,
+covering estimates of asymmetric divisor functions and moments of Riemann zeta-function.
diff --git a/Setup.hs b/Setup.hs
deleted file mode 100644
--- a/Setup.hs
+++ /dev/null
@@ -1,2 +0,0 @@
-import Distribution.Simple
-main = defaultMain
diff --git a/auxiliary/BenchMatrix.hs b/auxiliary/BenchMatrix.hs
--- a/auxiliary/BenchMatrix.hs
+++ b/auxiliary/BenchMatrix.hs
@@ -1,7 +1,7 @@
 module Main where
 
 import qualified Math.ExpPairs.Matrix3 as M3
-import Gauge.Main
+import Test.Tasty.Bench
 
 testm3 :: Int -> M3.Matrix3 Integer
 testm3 k = M3.fromList $ map (100*10^k `div`) [100..108]
diff --git a/auxiliary/OptimizeSum.hs b/auxiliary/OptimizeSum.hs
--- a/auxiliary/OptimizeSum.hs
+++ b/auxiliary/OptimizeSum.hs
@@ -13,14 +13,15 @@
 import Data.String
 import Text.RawString.QQ     (r)
 
-import Data.List (sortOn, subsequences, mapAccumL)
+import Data.List (sortOn, subsequences)
 import Data.Monoid
 import Data.Foldable
 import Data.Maybe
 import Data.Ord
 import Data.Function
+import Data.Traversable
 
-import Data.Text.Prettyprint.Doc hiding ((<>), group)
+import Prettyprinter hiding ((<>), group)
 
 import Prelude hiding (foldl, foldl1, maximum, mapM_, minimum, sum, concat)
 
diff --git a/exp-pairs.cabal b/exp-pairs.cabal
--- a/exp-pairs.cabal
+++ b/exp-pairs.cabal
@@ -1,5 +1,5 @@
 name:                exp-pairs
-version:             0.2.1.0
+version:             0.2.1.1
 synopsis:            Linear programming over exponent pairs
 description:         Package implements an algorithm to minimize rational objective function over the set of exponent pairs
 homepage:            https://github.com/Bodigrim/exp-pairs
@@ -9,14 +9,15 @@
 maintainer:          andrew.lelechenko@gmail.com
 category:            Math
 build-type:          Simple
-extra-source-files:  CHANGELOG.md, tests/*.txt
-cabal-version:       >=1.10
+extra-source-files:  tests/*.txt
+extra-doc-files:     README.md, CHANGELOG.md
+cabal-version:       1.18
 tested-with:
-  GHC ==8.4.4 GHC ==8.6.5 GHC ==8.8.3 GHC ==8.10.1
+  GHC ==8.4.4 GHC ==8.6.5 GHC ==8.8.4 GHC ==8.10.7 GHC ==9.0.2 GHC ==9.2.8 GHC ==9.4.8 GHC ==9.6.6 GHC ==9.8.4 GHC ==9.10.1 GHC ==9.12.1
 
 source-repository head
   type:     git
-  location: git://github.com/Bodigrim/exp-pairs.git
+  location: https://github.com/Bodigrim/exp-pairs.git
 
 library
   exposed-modules:     Math.ExpPairs,
@@ -31,14 +32,12 @@
                        Math.ExpPairs.ProcessMatrix,
                        Math.ExpPairs.RatioInf
   build-depends:       base >=4.11 && <5,
-                       ghc-prim,
-                       prettyprinter,
-                       deepseq >=1.3,
-                       containers
-  if impl(ghc < 8.0)
-    build-depends:     semigroups >= 0.8
+                       ghc-prim < 1,
+                       prettyprinter >= 1.7 && < 1.8,
+                       deepseq >=1.3 && <1.6,
+                       containers < 0.8
   default-language:    Haskell2010
-  ghc-options:         -Wall -fno-warn-type-defaults
+  ghc-options:         -Wall -fno-warn-type-defaults -Wcompat
 
 test-suite tests
   type:                exitcode-stdio-1.0
@@ -62,26 +61,26 @@
                        QuickCheck >=2.9,
                        smallcheck >=1.2,
                        exp-pairs,
-                       matrix >=0.1,
+                       matrix >=0.3.4,
                        random
   hs-source-dirs:      tests
   default-language:    Haskell2010
-  ghc-options:         -Wall -fno-warn-type-defaults
+  ghc-options:         -Wall -fno-warn-type-defaults -Wcompat -Wno-incomplete-uni-patterns
 
 benchmark matrix-bench
   build-depends:
     base,
     exp-pairs,
-    gauge
+    tasty-bench
   type: exitcode-stdio-1.0
   main-is: BenchMatrix.hs
   default-language: Haskell2010
   hs-source-dirs: auxiliary
-  ghc-options: -Wall -fno-warn-type-defaults
+  ghc-options: -Wall -fno-warn-type-defaults -Wcompat
 
 benchmark optimize-sum
   build-depends:
-    base,
+    base < 5,
     bimap,
     containers,
     prettyprinter,
@@ -90,4 +89,4 @@
   main-is: OptimizeSum.hs
   default-language: Haskell2010
   hs-source-dirs: auxiliary
-  ghc-options: -Wall -fno-warn-type-defaults
+  ghc-options: -Wall -fno-warn-type-defaults -Wcompat
diff --git a/tests/Etalon.hs b/tests/Etalon.hs
--- a/tests/Etalon.hs
+++ b/tests/Etalon.hs
@@ -2,7 +2,7 @@
 
 import System.Random
 import Data.Ord
-import Data.List
+import Data.List (sortBy)
 import Test.Tasty.HUnit
 
 unsort :: (RandomGen g) => g -> [x] -> [x]
@@ -22,4 +22,6 @@
   tests <- fetchRandomLines n filename
   let results = map f tests
   let fails = filter (not . snd) (zip tests results)
-  assertBool ("failed at " ++ show (fst $ head fails)) (null fails)
+  case fails of
+    [] -> pure ()
+    failure : _ -> assertFailure $ "failed at " ++ show (fst failure)
diff --git a/tests/Matrix3.hs b/tests/Matrix3.hs
--- a/tests/Matrix3.hs
+++ b/tests/Matrix3.hs
@@ -1,14 +1,19 @@
+{-# LANGUAGE CPP #-}
+
 module Matrix3 where
 
+import qualified Math.ExpPairs.Matrix3 as M3
+#ifdef MIN_VERSION_matrix
 import Data.Foldable
 import qualified Data.Matrix as M
-import qualified Math.ExpPairs.Matrix3 as M3
+#endif
 
 import Test.Tasty
 import Test.Tasty.QuickCheck as QC
 
 import Instances ()
 
+#ifdef MIN_VERSION_matrix
 toM :: M3.Matrix3 a -> M.Matrix a
 toM = M.fromList 3 3 . toList
 
@@ -20,47 +25,50 @@
   m'  = toM $ m1 `op1` m2
   m'' = toM m1 `op2` toM m2
 
-testMakarov :: M3.Matrix3 Integer -> M3.Matrix3 Integer -> Bool
-testMakarov m1 m2 = m1 * m2 == m1 `M3.makarovMult` m2
-
-testLaderman :: M3.Matrix3 Integer -> M3.Matrix3 Integer -> Bool
-testLaderman m1 m2 = m1 * m2 == m1 `M3.ladermanMult` m2
-
 testDet1 :: M3.Matrix3 Integer -> Bool
 testDet1 m = M3.det m == M.detLaplace (toM m)
 
 testDet2 :: M3.Matrix3 Rational -> Bool
 testDet2 m = M3.det m == M.detLU (toM m)
 
+testConv :: M3.Matrix3 Integer -> Bool
+testConv m = (toM3 . toM) m == m
+
+testMultCol :: M3.Matrix3 Integer -> (Integer, Integer, Integer) -> Bool
+testMultCol m v@(v1, v2, v3) = a==a' && b==b' && c==c' where
+  (a, b, c) = M3.multCol m v
+  [a', b', c'] = M.toList $ toM m * M.fromList 3 1 [v1, v2, v3]
+#endif
+
+testMakarov :: M3.Matrix3 Integer -> M3.Matrix3 Integer -> Bool
+testMakarov m1 m2 = m1 * m2 == m1 `M3.makarovMult` m2
+
+testLaderman :: M3.Matrix3 Integer -> M3.Matrix3 Integer -> Bool
+testLaderman m1 m2 = m1 * m2 == m1 `M3.ladermanMult` m2
+
 testRecip :: M3.Matrix3 Rational -> Bool
 testRecip m = M3.det m==0 || m/=m' && m==m'' && M3.det m * M3.det m' == 1 where
   m' = recip m
   m'' = recip m'
 
-testConv :: M3.Matrix3 Integer -> Bool
-testConv m = (toM3 . toM) m == m
-
 testNormalize :: Integer -> M3.Matrix3 Integer -> Bool
 testNormalize a m = (M3.normalize m' == m') && (a==0 || a>0 && m'==m'' || a<0 && m'==negate m'') where
   m' = M3.normalize m
   m'' = M3.normalize (m * fromInteger a)
 
-testMultCol :: M3.Matrix3 Integer -> (Integer, Integer, Integer) -> Bool
-testMultCol m v@(v1, v2, v3) = a==a' && b==b' && c==c' where
-  (a, b, c) = M3.multCol m v
-  [a', b', c'] = M.toList $ toM m * M.fromList 3 1 [v1, v2, v3]
-
 testSuite :: TestTree
 testSuite = testGroup "Matrix3"
-  [ QC.testProperty "plus"      $ testOp (+) (+)
+  [ QC.testProperty "makarov"     testMakarov
+  , QC.testProperty "laderman"    testLaderman
+  , QC.testProperty "recip"       testRecip
+  , QC.testProperty "normalize"   testNormalize
+#ifdef MIN_VERSION_matrix
+  , QC.testProperty "plus"      $ testOp (+) (+)
   , QC.testProperty "minus"     $ testOp (-) (-)
   , QC.testProperty "mult"      $ testOp (*) (*)
-  , QC.testProperty "makarov"     testMakarov
-  , QC.testProperty "laderman"    testLaderman
   , QC.testProperty "det1"        testDet1
-  , QC.testProperty "conversion"  testConv
   , QC.testProperty "det2"        testDet2
-  , QC.testProperty "recip"       testRecip
-  , QC.testProperty "normalize"   testNormalize
+  , QC.testProperty "conversion"  testConv
   , QC.testProperty "mult column" testMultCol
+#endif
   ]