diff --git a/CHANGELOG.md b/CHANGELOG.md
new file mode 100644
--- /dev/null
+++ b/CHANGELOG.md
@@ -0,0 +1,15 @@
+Changes
+=======
+
+Version 0.1.4.0
+----------------
+
+Improve overall performance.
+Use Stern-Brocot tree for binary searches in Math.ExpPairs.Ivic.
+
+Version 0.1.3.0
+----------------
+
+New functions in Math.ExpPairs.Ivic: reverseMOnS, checkAbscissa, findMinAbscissa, mBigOnHalf, reverseMBigOnHalf.
+Fast mastrix multiplication via Makarov and Laderman algorithms.
+Rewrite from the scratch pretty printer of processes.
diff --git a/Math/ExpPairs.hs b/Math/ExpPairs.hs
--- a/Math/ExpPairs.hs
+++ b/Math/ExpPairs.hs
@@ -14,78 +14,82 @@
 A set of useful applications can be found in
 "Math.ExpPairs.Ivic", "Math.ExpPairs.Kratzel" and "Math.ExpPairs.MenzerNowak".
 -}
+{-# LANGUAGE CPP  #-}
+
 module Math.ExpPairs
-	( optimize
-	, OptimizeResult
-	, optimalValue
-	, optimalPair
-	, optimalPath
-	, simulateOptimize
-	, simulateOptimize'
-	, LinearForm (..)
-	, RationalForm (..)
-	, IneqType (..)
-	, Constraint (..)
-	, InitPair
-	, Path
-	, RatioInf (..)
-	, RationalInf
-	) where
+  ( optimize
+  , OptimizeResult
+  , optimalValue
+  , optimalPair
+  , optimalPath
+  , simulateOptimize
+  , simulateOptimize'
+  , LinearForm (..)
+  , RationalForm (..)
+  , IneqType (..)
+  , Constraint (..)
+  , InitPair
+  , Path
+  , RatioInf (..)
+  , RationalInf
+  ) where
 
-import Data.Ratio  ((%), numerator, denominator)
-import Data.Ord    (comparing)
-import Data.List   (minimumBy)
-import Data.Monoid (mempty, mappend)
+import Control.Arrow hiding ((<+>))
+import Data.Function (on)
+import Data.Ord      (comparing)
+import Data.List     (minimumBy)
+import Data.Monoid
+import Text.PrettyPrint.Leijen hiding ((<$>), (<>))
+import qualified Text.PrettyPrint.Leijen as PP
+import Text.Printf
 
 import Math.ExpPairs.LinearForm
 import Math.ExpPairs.Process
 import Math.ExpPairs.Pair
 import Math.ExpPairs.RatioInf
 
-fracs2proj :: (Rational, Rational) -> (Integer, Integer, Integer)
-fracs2proj (q, r) = (k, l, m) where
-	dq = denominator q
-	dr = denominator r
-	m = lcm dq dr
-	k = numerator q * (m `div` dq)
-	l = numerator r * (m `div` dr)
-
 evalFunctional :: [InitPair] -> [InitPair] -> [RationalForm Rational] -> [Constraint Rational] -> Path -> (RationalInf, InitPair)
-evalFunctional corners interiors rfs cons path = if null rs then (InfPlus, undefined) else minimumBy (comparing fst) rs where
-	applyPath ips = map (evalPath path . fracs2proj . initPairToValue) ips `zip` ips
-	corners'   = applyPath corners
-	interiors' = applyPath interiors
+evalFunctional corners interiors rfs cons path = case rs of
+  [] -> (InfPlus, undefined)
+  _  -> minimumBy (comparing fst) rs
+  where
+    applyPath  = map (evalPath path . initPairToProjValue &&& id)
+    corners'   = applyPath corners
+    interiors' = applyPath interiors
 
-	predicate (p, _) = all (checkConstraint p) cons
-	qs = if all predicate corners' then corners' else filter predicate interiors'
+    predicate (p, _) = all (checkConstraint p) cons
+    qs = if all predicate corners'
+          then corners'
+          else filter predicate interiors'
 
-	rs = map (\(p, ip) -> (maximum $ map (evalRF p) rfs, ip)) qs
+    rs = map (first $ \p -> maximum (map (evalRF p) rfs)) qs
 
 checkMConstraints :: Path -> [Constraint Rational] -> Bool
-checkMConstraints path = all (\con -> any (\p -> checkConstraint (evalPath path p) con ) triangleT) where
-	triangleT = map fracs2proj [ (0%1,1%1), (0%1,1%2), (1%2,1%2)]
+checkMConstraints path = all (\con -> any (\p -> checkConstraint (evalPath path p) con) triangleT) where
+  triangleT = [(0, 1, 1), (0, 1, 2), (1, 1, 2)]
 
 -- |Container for the result of optimization.
 data OptimizeResult = OptimizeResult {
-	-- | The minimal value of objective function.
-	optimalValue :: RationalInf,
-	-- | The initial exponent pair, on which minimal value was achieved.
-	optimalPair  :: InitPair,
-	-- | The sequence of processes, after which minimal value was
-	-- achieved.
-	optimalPath  :: Path
-	}
+  -- | The minimal value of objective function.
+  optimalValue :: RationalInf,
+  -- | The initial exponent pair, on which minimal value was achieved.
+  optimalPair  :: InitPair,
+  -- | The sequence of processes, after which minimal value was
+  -- achieved.
+  optimalPath  :: Path
+  }
+  deriving (Show)
 
-instance Show OptimizeResult where
-	show (OptimizeResult r' ip p) = show' r' ++ "\n" ++ show ip ++ "\t" ++ show p where
-		show' (Finite r) = show (fromRational r :: Double) ++ " = " ++ show r
-		show' r = show r
+instance Pretty OptimizeResult where
+  pretty (OptimizeResult r' ip p) = pretty' r' PP.<$> pretty ip </> pretty p where
+    pretty' r@(Finite rr) = text (printf "%.6f" (fromRational rr :: Double)) <+> equals <+> pretty r
+    pretty' r = pretty r
 
 instance Eq OptimizeResult where
-	a==b = optimalValue a == optimalValue b
+  (==) = (==) `on` optimalValue
 
 instance Ord OptimizeResult where
-	compare a b = compare (optimalValue a) (optimalValue b)
+  compare = compare `on` optimalValue
 
 -- |Wrap 'Rational' into 'OptimizeResult'.
 simulateOptimize :: Rational -> OptimizeResult
@@ -100,28 +104,28 @@
 -- all constraints and minimizes the maximum of all rational forms.
 optimize :: [RationalForm Rational] -> [Constraint Rational] -> OptimizeResult
 optimize rfs cons = optimize' rfs cons (OptimizeResult r0 ip0 mempty) where
-	(r0, ip0) = evalFunctional [Corput01, Corput12] [Corput01, Corput12] rfs cons mempty
+  (r0, ip0) = evalFunctional [Corput01, Corput12] [Corput01, Corput12] rfs cons mempty
 
 optimize' :: [RationalForm Rational] -> [Constraint Rational] -> OptimizeResult -> OptimizeResult
 optimize' rfs cons ret@(OptimizeResult r _ path)
-	| lengthPath path > 100 = ret
-	| otherwise = retBA where
-		ret0@(OptimizeResult r0 ip0 _) = if r0' < r then OptimizeResult r0' ip0' path else ret where
-			(r0', ip0') = evalFunctional corners interiors rfs cons path
-			corners = [Mix 1 0, Mix 0 1, Mix 0 0]
-			interiors = initPairs
+  | lengthPath path > 100 = ret
+  | otherwise = retBA where
+    ret0@(OptimizeResult r0 ip0 _) = if r0' < r then OptimizeResult r0' ip0' path else ret where
+      (r0', ip0') = evalFunctional corners interiors rfs cons path
+      corners = [Mix 1 0, Mix 0 1, Mix 0 0]
+      interiors = initPairs
 
-		cons0 = if r0==InfPlus then cons else cons ++ map (consBuilder r0) rfs
+    cons0 = if r0==InfPlus then cons else cons ++ map (consBuilder r0) rfs
 
-		retA@(OptimizeResult r1 ip1 _) = if checkMConstraints patha cons0 && r1' < r0 then branchA else ret0 where
-			patha  = path `mappend` aPath
-			branchA@(OptimizeResult r1' _ _) = optimize' rfs cons (OptimizeResult r0 ip0 patha)
+    retA@(OptimizeResult r1 ip1 _) = if checkMConstraints patha cons0 && r1' < r0 then branchA else ret0 where
+      patha  = path <> aPath
+      branchA@(OptimizeResult r1' _ _) = optimize' rfs cons (OptimizeResult r0 ip0 patha)
 
-		cons1 = if r1==r0	then cons0 else cons ++ map (consBuilder r1) rfs
+    cons1 = if r1==r0  then cons0 else cons ++ map (consBuilder r1) rfs
 
-		retBA = if checkMConstraints pathba cons1 && r2' < r1 then branchB else retA where
-			pathba  = path `mappend` baPath
-			branchB@(OptimizeResult r2' _ _) = optimize' rfs cons (OptimizeResult r1 ip1 pathba)
+    retBA = if checkMConstraints pathba cons1 && r2' < r1 then branchB else retA where
+      pathba  = path <> baPath
+      branchB@(OptimizeResult r2' _ _) = optimize' rfs cons (OptimizeResult r1 ip1 pathba)
 
-		consBuilder rr (RationalForm num den) = Constraint (substituteLF (num, den, 1) (LinearForm (-1) (toRational rr) 0)) Strict
+    consBuilder rr (RationalForm num den) = Constraint (substituteLF (num, den, 1) (LinearForm (-1) (toRational rr) 0)) Strict
 
diff --git a/Math/ExpPairs/Ivic.hs b/Math/ExpPairs/Ivic.hs
--- a/Math/ExpPairs/Ivic.hs
+++ b/Math/ExpPairs/Ivic.hs
@@ -13,17 +13,17 @@
 
 -}
 module Math.ExpPairs.Ivic
-	( zetaOnS
-	, reverseZetaOnS
-	, mOnS
-	, reverseMOnS
-	, checkAbscissa
-	, findMinAbscissa
-	, mBigOnHalf
-	, reverseMBigOnHalf
-	) where
+  ( zetaOnS
+  , reverseZetaOnS
+  , mOnS
+  , reverseMOnS
+  , checkAbscissa
+  , findMinAbscissa
+  , mBigOnHalf
+  , reverseMBigOnHalf
+  ) where
 
-import Data.Ratio ((%))
+import Data.Ratio
 import Data.List  (minimumBy)
 import Data.Ord   (comparing)
 
@@ -33,13 +33,13 @@
 -- See equation (7.57) in Ivić2003.
 zetaOnS :: Rational -> OptimizeResult
 zetaOnS s
-	| s >= 1  = simulateOptimize 0
-	| s >= 1%2 = optimize
-		[RationalForm (LinearForm 1 1 (-s)) 2]
-		[Constraint (LinearForm (-1) 1 (-s)) NonStrict]
-	| otherwise = optRes {optimalValue = r} where
-		optRes = zetaOnS (1-s)
-		r = Finite (1%2 - s) + optimalValue optRes
+  | s >= 1  = simulateOptimize 0
+  | s >= 1%2 = optimize
+    [RationalForm (LinearForm 1 1 (-s)) 2]
+    [Constraint (LinearForm (-1) 1 (-s)) NonStrict]
+  | otherwise = optRes {optimalValue = r} where
+    optRes = zetaOnS (1-s)
+    r = Finite (1%2 - s) + optimalValue optRes
 
 zetaOnHalf :: Rational
 zetaOnHalf = 32%205
@@ -47,64 +47,56 @@
 -- | An attempt to reverse 'zetaOnS'.
 reverseZetaOnS :: Rational -> OptimizeResult
 reverseZetaOnS mu
-	| mu >= 1%2   = simulateOptimize 0
-	| mu > zetaOnHalf = optimize [RationalForm (LinearForm 1 (-1) 1) 1] [Constraint (LinearForm 0 (-2) (1+2*mu)) NonStrict]
-	| otherwise = optRes {optimalValue = negate $ optimalValue optRes} where
-	optRes = optimize [RationalForm (LinearForm 1 (-1) 0) 1] [Constraint (LinearForm 1 0 (-mu)) NonStrict, Constraint (LinearForm (-1) 1 (-1%2)) NonStrict]
+  | mu >= 1%2   = simulateOptimize 0
+  | mu > zetaOnHalf = optimize [RationalForm (LinearForm 1 (-1) 1) 1] [Constraint (LinearForm 0 (-2) (1+2*mu)) NonStrict]
+  | mu == zetaOnHalf = simulateOptimize (1 % 2)
+  | otherwise = optRes {optimalValue = negate $ optimalValue optRes} where
+  optRes = optimize [RationalForm (LinearForm 1 (-1) 0) 1] [Constraint (LinearForm 1 0 (-mu)) NonStrict, Constraint (LinearForm (-1) 1 (-1%2)) NonStrict]
 
 lemma82_f :: Rational -> Rational
 lemma82_f s
-	| s < 1%2   = undefined
-	| s<= 2%3   =  2/(3-4*s)
-	| s<=11%14  = 10/(7-8*s)
-	| s<=13%15  = 34/(15-16*s)
-	| s<=57%62  = 98/(31-32*s)
-	| otherwise =	 5/(1-s)
-
--- Ivic, (8.97)
--- R << T V^{-2f(sigma)} + T^alpha1 V^beta1 + T^alpha2 V^beta2
---
--- If a<1 then T^a V^b << T V^{b+(a-1)/muS}
---
--- (8.97) implies that alpha1 <= 1 for S >= 1/2
--- and that alpha2 <= 1 for S >= 2/3 or S >= 5/8 and
---          (4S-2)k + (8S-6)l + 2S-1 >=0
+  | s < 1%2   = undefined
+  | s<= 2%3   =  2/(3-4*s)
+  | s<=11%14  = 10/(7-8*s)
+  | s<=13%15  = 34/(15-16*s)
+  | s<=57%62  = 98/(31-32*s)
+  | otherwise =   5/(1-s)
 
 -- | Compute maximal m(σ) such that ∫_1^T |ζ(σ+it)|^m(σ) dt ≪ T^(1+ε).
--- See equation (8.97) in Ivić2003.
+-- See equation (8.97) in Ivić2003. Further justification will be published elsewhere.
 mOnS :: Rational -> OptimizeResult
 mOnS s
-	| s < 1%2 = simulateOptimize 0
-	| s < 5%8 = simulateOptimize $ 4/(3-4*s)
-	| s>= 1   = simulateOptimize' InfPlus
-	| otherwise = minimumBy (comparing optimalValue) [x1, x2, simulateOptimize (lemma82_f s * 2)] where
+  | s < 1%2 = simulateOptimize 0
+  | s < 5%8 = simulateOptimize $ 4/(3-4*s)
+  | s>= 1   = simulateOptimize' InfPlus
+  | otherwise = minimumBy (comparing optimalValue) [x1, x2, simulateOptimize (lemma82_f s * 2)] where
 
-		optRes = zetaOnS s
-		muS    = toRational $ optimalValue optRes
-		alpha1 = (4-4*s)/(1+2*s)
-		beta1  = -12/(1+2*s)
-		x1 = optRes {optimalValue = Finite $ (1-alpha1)/muS - beta1}
+    optRes = zetaOnS s
+    muS    = toRational $ optimalValue optRes
+    alpha1 = (4-4*s)/(1+2*s)
+    beta1  = -12/(1+2*s)
+    x1 = optRes {optimalValue = Finite $ (1-alpha1)/muS - beta1}
 
-		--alpha2 = 4*(1-s)*(k+l)/((2*m+4*l)*s-m+2*k-2*l)
-		--beta2  = -4*(m+2*k+2*l)/((2*m+4*l)*s-m+2*k-2*l)
-		--ratio = (1-alpha2)/muS - beta2
-		--numer = numerator ratio
-		--denom = denominator ratio
-		numer = LinearForm
-			(-4*s + (-8*muS + 2))
-			(-8*s + (-8*muS + 6))
-			(-2*s + (-4*muS + 1))
-		denom = LinearForm
-			(2*muS)
-			(4*muS*s - 2*muS)
-			(2*muS*s - muS)
+    --alpha2 = 4*(1-s)*(k+l)/((2*m+4*l)*s-m+2*k-2*l)
+    --beta2  = -4*(m+2*k+2*l)/((2*m+4*l)*s-m+2*k-2*l)
+    --ratio = (1-alpha2)/muS - beta2
+    --numer = numerator ratio
+    --denom = denominator ratio
+    numer = LinearForm
+      (-4*s + (-8*muS + 2))
+      (-8*s + (-8*muS + 6))
+      (-2*s + (-4*muS + 1))
+    denom = LinearForm
+      (2*muS)
+      (4*muS*s - 2*muS)
+      (2*muS*s - muS)
 
-		cons = if s >= 2%3 then [] else [Constraint
-			(LinearForm (4*s-2) (8*s-6) (2*s-1)) NonStrict
-			]
+    cons = if s >= 2%3 then [] else [Constraint
+      (LinearForm (4*s-2) (8*s-6) (2*s-1)) NonStrict
+      ]
 
-		x2' = optimize [RationalForm numer denom] cons
-		x2 = x2' {optimalValue = negate $ optimalValue x2'}
+    x2' = optimize [RationalForm numer denom] cons
+    x2 = x2' {optimalValue = negate $ optimalValue x2'}
 
 -- | Try to reverse 'mOnS': for a given precision and m compute minimal possible σ.
 -- Implementation is usual try-and-divide search, so performance is very poor.
@@ -112,43 +104,47 @@
 -- real σ and returns bigger value.
 reverseMOnS :: Rational -> RationalInf -> Rational
 reverseMOnS prec m = reverseMOnS' from to where
-	from = 1 % 2
-	to   = 1 % 1
-	reverseMOnS' a b
-		| b-a < prec = a
-		| optimalValue (mOnS ((a+b)/2)) > m = reverseMOnS' a ((a+b)/2)
-		| otherwise = reverseMOnS' ((a+b)/2) b
+  from = 1 % 2
+  to   = 1
+  reverseMOnS' a b
+    | b - a < prec = c
+    | optimalValue (mOnS c) > m = reverseMOnS' a c
+    | otherwise = reverseMOnS' c b
+    where
+      c = (numerator a + numerator b) % (denominator a + denominator b)
 
--- | Check whether ∫_1^T 	Π |ζ(n_i*σ+it)|^m_i dt ≪ T^(1+ε) for a given list of pairs [(n_1, m_1), ...] and fixed σ.
+-- | Check whether ∫_1^T   Π_i |ζ(n_i*σ+it)|^m_i dt ≪ T^(1+ε) for a given list of pairs [(n_1, m_1), ...] and fixed σ.
 checkAbscissa :: [(Rational, Rational)] -> Rational -> Bool
 checkAbscissa xs s = sum rs < Finite 1 where
-	qs = map (\(n,m) -> optimalValue (mOnS (n*s)) / Finite m) xs
-	rs = map (\q -> 1/q) qs
+  qs = map (\(n,m) -> optimalValue (mOnS (n*s)) / Finite m) xs
+  rs = map (\q -> 1/q) qs
 
 -- | Find for a given precision and list of pairs [(n_1, m_1), ...] the minimal σ
--- such that ∫_1^T 	Π |ζ(n_i*σ+it)|^m_i dt ≪ T^(1+ε).
+-- such that ∫_1^T   Π_i|ζ(n_i*σ+it)|^m_i dt ≪ T^(1+ε).
 findMinAbscissa :: Rational -> [(Rational, Rational)] -> Rational
 findMinAbscissa prec xs = searchMinAbscissa' from to where
-	from = 1 % 2 / minimum (map fst xs)
-	to   = 1 % 1
-	searchMinAbscissa' a b
-		| b-a < prec = a
-		| checkAbscissa xs ((a+b)/2) = searchMinAbscissa' a ((a+b)/2)
-		| otherwise = searchMinAbscissa' ((a+b)/2) b
+  from = 1 % 2 / minimum (map fst xs)
+  to   = 1 % 1
+  searchMinAbscissa' a b
+    | b - a < prec = b
+    | checkAbscissa xs c = searchMinAbscissa' a c
+    | otherwise = searchMinAbscissa' c b
+    where
+      c = (numerator a + numerator b) % (denominator a + denominator b)
 
 -- | Compute minimal M(A) such that ∫_1^T |ζ(1/2+it)|^A dt ≪ T^(M(A)+ε).
 -- See Ch. 8 in Ivić2003. Further justification will be published elsewhere.
 mBigOnHalf :: Rational -> OptimizeResult
 mBigOnHalf a
-	| a < 4     = simulateOptimize 1
-	| a < 12    = simulateOptimize $ 1+(a-4)/8
-	| a > 41614060315296730740083860226662 % 2636743270445733804969041895717 = simulateOptimize $ 1 + 32*(a-6)/205
-	| otherwise = if Finite x >= optimalValue optRes
-		then simulateOptimize x
-		else optRes where
-			optRes = optimize [RationalForm (LinearForm 1 1 0) (LinearForm 1 0 0)]
-				[Constraint (LinearForm (4-a) 4 2) NonStrict]
-			x = 1 + 32*(a-6)/205
+  | a < 4     = simulateOptimize 1
+  | a < 12    = simulateOptimize $ 1+(a-4)/8
+  | a > 41614060315296730740083860226662 % 2636743270445733804969041895717 = simulateOptimize $ 1 + 32*(a-6)/205
+  | otherwise = if Finite x >= optimalValue optRes
+    then simulateOptimize x
+    else optRes where
+      optRes = optimize [RationalForm (LinearForm 1 1 0) (LinearForm 1 0 0)]
+        [Constraint (LinearForm (4-a) 4 2) NonStrict]
+      x = 1 + 32*(a-6)/205
 -- Constant 41614060315296730740083860226662 % 2636743270445733804969041895717
 -- is produced by
 -- optimize [RationalForm (LinearForm 4 4 2) (LinearForm 1 0 0)] [Constraint (LinearForm (-64) (-77) 64) Strict]
@@ -158,11 +154,11 @@
 -- real A and returns lower value.
 reverseMBigOnHalf :: Rational -> OptimizeResult
 reverseMBigOnHalf m
-	| m <= 2 = simulateOptimize $ (m-1)*8 + 4
-	| otherwise = if Finite a <= optimalValue optRes
-		then simulateOptimize a
-		else optRes where
-		a = (m-1)*205/32 + 6
-		optRes = optimize [RationalForm (LinearForm 4 4 2) (LinearForm 1 0 0)] [Constraint (LinearForm (1-m) 1 0) NonStrict]
+  | m <= 2 = simulateOptimize $ (m-1)*8 + 4
+  | otherwise = if Finite a <= optimalValue optRes
+    then simulateOptimize a
+    else optRes where
+    a = (m-1)*205/32 + 6
+    optRes = optimize [RationalForm (LinearForm 4 4 2) (LinearForm 1 0 0)] [Constraint (LinearForm (1-m) 1 0) NonStrict]
 
 
diff --git a/Math/ExpPairs/Kratzel.hs b/Math/ExpPairs/Kratzel.hs
--- a/Math/ExpPairs/Kratzel.hs
+++ b/Math/ExpPairs/Kratzel.hs
@@ -14,9 +14,9 @@
 (v, w, z) with v^a w^b z^c = n.
 
 Krätzel
-	(/Krätzel E./
-	`Lattice points'.
-	Dordrecht: Kluwer, 1988)
+  (/Krätzel E./
+  `Lattice points'.
+  Dordrecht: Kluwer, 1988)
 proved asymptotic formulas for
 Σ_{n ≤ x} τ_{a, b}(n) with an error term of order x^(Θ(a, b) + ε)
 and for
@@ -25,105 +25,123 @@
 
 -}
 module Math.ExpPairs.Kratzel
-	( TauabTheorem (..)
-	, tauab
-	, TauabcTheorem (..)
-	, tauabc
-	) where
+  ( TauabTheorem (..)
+  , tauab
+  , TauabcTheorem (..)
+  , tauabc
+  ) where
 
+import Control.Arrow
 import Data.Ratio ((%))
 import Data.Ord   (comparing)
 import Data.List  (minimumBy)
+import Text.PrettyPrint.Leijen
 
 import Math.ExpPairs
 
 -- |Special type to specify the theorem of Krätzel1988,
 -- which provided the best estimate of Θ(a, b)
 data TauabTheorem
-	-- | Theorem 5.11, case a)
-	= Kr511a
-	-- | Theorem 5.11, case b)
-	| Kr511b
-	-- | Theorem 5.12, case a)
-	| Kr512a
-	-- | Theorem 5.12, case b)
-	| Kr512b
-	deriving (Show)
+  -- | Theorem 5.11, case a)
+  = Kr511a
+  -- | Theorem 5.11, case b)
+  | Kr511b
+  -- | Theorem 5.12, case a)
+  | Kr512a
+  -- | Theorem 5.12, case b)
+  | Kr512b
+  deriving (Eq, Ord, Enum, Bounded, Show)
 
+instance Pretty TauabTheorem where
+  pretty = text . show
+
+divideResult :: Real a => a -> (b, OptimizeResult) -> (b, OptimizeResult)
+divideResult d = second (\o -> o {optimalValue = optimalValue o / Finite (toRational d)})
+
 -- |Compute Θ(a, b) for given a and b.
 tauab :: Integer -> Integer -> (TauabTheorem, OptimizeResult)
+tauab a' b'
+  | d /= 1 = divideResult d $ tauab (a'`div` d) (b' `div` d) where
+      d = gcd a' b'
 tauab a' b' = minimumBy (comparing (optimalValue . snd)) [kr511a, kr511b, kr512a, kr512b] where
-	a = a'%1
-	b = b'%1
-	kr511a = (Kr511a, optimize
-		[RationalForm (LinearForm 2 2 (-1)) (LinearForm 0 0 (a+b))]
-		[Constraint (LinearForm (-2*b) (2*a) (-a)) NonStrict])
-	kr511b = (Kr511b, optimize
-		[RationalForm (LinearForm 1 0 0) (LinearForm b (-a) a)]
-		[Constraint (LinearForm (2*b) (-2*a) a) Strict])
-	kr512a = (Kr512a, simulateOptimize r) where
-		r = if 11*a >= 8*b then 19/29/(a+b) else 1%1
-	kr512b = if 11*a >= 8*b then kr512a else (Kr512b, optimize
-		[
-			RationalForm (LinearForm (-11) 8 (-4)) (LinearForm (-29*b) (29*a) (4*b-20*a))
-		]
-		[
-			Constraint (LinearForm (-2*b) (2*a) (-a)) NonStrict,
-			Constraint (LinearForm (-29) 0 4) Strict,
-			Constraint (LinearForm 29 29 (-24)) Strict
-		])
+  a = toRational a'
+  b = toRational b'
+  kr511a = (Kr511a, optimize
+    [RationalForm (LinearForm 2 2 (-1)) (LinearForm 0 0 (a+b))]
+    [Constraint (LinearForm (-2*b) (2*a) (-a)) NonStrict])
+  kr511b = (Kr511b, optimize
+    [RationalForm (LinearForm 1 0 0) (LinearForm b (-a) a)]
+    [Constraint (LinearForm (2*b) (-2*a) a) Strict])
+  kr512a = (Kr512a, simulateOptimize r) where
+    r = if 11*a >= 8*b then 19/29/(a+b) else 1%1
+  kr512b = if 11*a >= 8*b then kr512a else (Kr512b, optimize
+    [
+      RationalForm (LinearForm (-11) 8 (-4)) (LinearForm (-29*b) (29*a) (4*b-20*a))
+    ]
+    [
+      Constraint (LinearForm (-2*b) (2*a) (-a)) NonStrict,
+      Constraint (LinearForm (-29) 0 4) Strict,
+      Constraint (LinearForm 29 29 (-24)) Strict
+    ])
 
 -- |Special type to specify the theorem of Krätzel1988,
 -- which provided the best estimate of Θ(a, b, c)
 data TauabcTheorem
-	-- | Kolesnik
-	-- (/Kolesnik G./ `On the estimation of multiple exponential sums'
-	-- \/\/ Recent progress in analytic number theory,
-	-- London: Academic Press, 1981, Vol. 1, P. 231–246)
-	-- proved that  Θ(1, 1, 1) = 43 \/96.
-	= Kolesnik
-	-- | Theorem 6.1
-	| Kr61
-	-- | Theorem 6.2
-	| Kr62
-	-- | Theorem 6.3
-	| Kr63
-	-- | Theorem 6.4
-	| Kr64
-	-- | Theorem 6.5
-	| Kr65
-	-- | Theorem 6.6
-	| Kr66
-	-- | In certain cases Θ(a, b, c) = Θ(a, b).
-	| Tauab TauabTheorem
-	deriving (Show)
+  -- | Kolesnik
+  -- (/Kolesnik G./ `On the estimation of multiple exponential sums'
+  -- \/\/ Recent progress in analytic number theory,
+  -- London: Academic Press, 1981, Vol. 1, P. 231–246)
+  -- proved that  Θ(1, 1, 1) = 43 \/96.
+  = Kolesnik
+  -- | Theorem 6.1
+  | Kr61
+  -- | Theorem 6.2
+  | Kr62
+  -- | Theorem 6.3
+  | Kr63
+  -- | Theorem 6.4
+  | Kr64
+  -- | Theorem 6.5
+  | Kr65
+  -- | Theorem 6.6
+  | Kr66
+  -- | In certain cases Θ(a, b, c) = Θ(a, b).
+  | Tauab TauabTheorem
+  deriving (Eq, Ord, Show)
 
+instance Pretty TauabcTheorem where
+  pretty (Tauab t) = pretty t
+  pretty t         = pretty (show t)
+
 -- |Compute Θ(a, b, c) for given a, b and c.
 tauabc :: Integer -> Integer -> Integer -> (TauabcTheorem, OptimizeResult)
+tauabc a' b' c'
+  | d /= 1 = divideResult d $ tauabc (a'`div` d) (b' `div` d) (c' `div` d) where
+      d = gcd (gcd a' b') c'
 tauabc 1 1 1 = (Kolesnik, simulateOptimize $ 43%96)
 tauabc a' b' c' = minimumBy (comparing (optimalValue . snd)) [kr61, kr62, kr63, kr64, kr65, kr66] where
-	a = a'%1
-	b = b'%1
-	c = c'%1
-	kr61
-		| c<a+b = (Kr61, simulateOptimize $ 2/(a+b+c))
-		| optimalValue optRes < Finite (recip c) = (Kr61, simulateOptimize $ 1/c)
-		| otherwise = (Tauab th, optRes)
-		where
-			(th, optRes) = tauab a' b'
-	kr62 = (Kr62, optimize
-		[RationalForm (LinearForm 2 2 0) (LinearForm 0 0 (a+b+c))]
-		[
-			Constraint (LinearForm (-b-c) a 0) NonStrict,
-			Constraint (LinearForm (-2*c) (-2*c) (a+b+c)) NonStrict
-		])
-	kr63 = (Kr63, optimize
-		[RationalForm (LinearForm 4 2 3) (LinearForm (2*(a+b+c)) 0 (3*(a+b+c)))]
-		[Constraint (LinearForm (2*(a-b-c)) (2*a) (2*a-b-c)) NonStrict])
-	kr64 = (Kr64, simulateOptimize r) where
-		r = recip (a+b+c) * minimum ((a+b+c):[2-4*(k-1)%(3*2^k-4) | k<-[1..maxk], (3*2^k-2*k-4)%1 * a >= 2 * (b+c), (3*2^k-8)%1 * (a+b) >= (3*2^k-4*k+4)%1 * c])
-		maxk = 4 `max` floor (logBase 2 (fromRational $ b+c) :: Double)
-	kr65 = (Kr65, simulateOptimize r) where
-		r = if 7*a>=2*(b+c) && 4*(a+b)>=5*c then 3%2/(a+b+c) else 1%1
-	kr66 = (Kr66, simulateOptimize r) where
-		r = if 18*a>=7*(b+c) && 2*(a+b)>=3*c then 25%17/(a+b+c) else 1%1
+  a = toRational a'
+  b = toRational b'
+  c = toRational c'
+  kr61
+    | c<a+b = (Kr61, simulateOptimize $ 2/(a+b+c))
+    | optimalValue optRes < Finite (recip c) = (Kr61, simulateOptimize $ 1/c)
+    | otherwise = (Tauab th, optRes)
+    where
+      (th, optRes) = tauab a' b'
+  kr62 = (Kr62, optimize
+    [RationalForm (LinearForm 2 2 0) (LinearForm 0 0 (a+b+c))]
+    [
+      Constraint (LinearForm (-b-c) a 0) NonStrict,
+      Constraint (LinearForm (-2*c) (-2*c) (a+b+c)) NonStrict
+    ])
+  kr63 = (Kr63, optimize
+    [RationalForm (LinearForm 4 2 3) (LinearForm (2*(a+b+c)) 0 (3*(a+b+c)))]
+    [Constraint (LinearForm (2*(a-b-c)) (2*a) (2*a-b-c)) NonStrict])
+  kr64 = (Kr64, simulateOptimize r) where
+    r = recip (a+b+c) * minimum ((a+b+c):[2-4*(k-1)%(3*2^k-4) | k<-[1..maxk], (3*2^k-2*k-4)%1 * a >= 2 * (b+c), (3*2^k-8)%1 * (a+b) >= (3*2^k-4*k+4)%1 * c])
+    maxk = 4 `max` floor (logBase 2 (fromRational $ b+c) :: Double)
+  kr65 = (Kr65, simulateOptimize r) where
+    r = if 7*a>=2*(b+c) && 4*(a+b)>=5*c then 3%2/(a+b+c) else 1%1
+  kr66 = (Kr66, simulateOptimize r) where
+    r = if 18*a>=7*(b+c) && 2*(a+b)>=3*c then 25%17/(a+b+c) else 1%1
diff --git a/Math/ExpPairs/LinearForm.hs b/Math/ExpPairs/LinearForm.hs
--- a/Math/ExpPairs/LinearForm.hs
+++ b/Math/ExpPairs/LinearForm.hs
@@ -1,4 +1,4 @@
-{-# LANGUAGE DeriveFunctor, DeriveFoldable, DeriveGeneric #-}
+{-# LANGUAGE DeriveFunctor, DeriveFoldable, DeriveGeneric, CPP #-}
 {-|
 Module      : Math.ExpPairs.LinearForm
 Description : Linear forms, rational forms and constraints
@@ -11,55 +11,56 @@
 Provides types for rational forms (to hold objective functions in "Math.ExpPairs") and linear contraints (to hold constraints of optimization). Both of them are built atop of projective linear forms.
 -}
 module Math.ExpPairs.LinearForm
-	( LinearForm (..)
-	, evalLF
-	, substituteLF
-	, RationalForm (..)
-	, evalRF
-	, IneqType (..)
-	, Constraint (..)
-	, checkConstraint
-	) where
+  ( LinearForm (..)
+  , evalLF
+  , substituteLF
+  , RationalForm (..)
+  , evalRF
+  , IneqType (..)
+  , Constraint (..)
+  , checkConstraint
+  ) where
 
 import Control.DeepSeq
 import Data.Foldable  (Foldable (..), toList)
-import Data.List      (intercalate)
-import Data.Ratio     (numerator, denominator)
+import Data.Maybe     (mapMaybe)
+#if __GLASGOW_HASKELL__ < 710
 import Data.Monoid    (Monoid, mempty, mappend)
+#endif
+import Data.Ratio     (numerator, denominator)
 import GHC.Generics   (Generic (..))
+import Text.PrettyPrint.Leijen
 
 import Math.ExpPairs.RatioInf
 
 -- |Define an affine linear form of two variables: a*k + b*l + c*m.
 -- First argument of 'LinearForm' stands for a, second for b
 -- and third for c. Linear forms form a monoid by addition.
-data LinearForm t = LinearForm t t t
-	deriving (Eq, Functor, Foldable, Generic)
+data LinearForm t = LinearForm !t !t !t
+  deriving (Eq, Show, Functor, Foldable, Generic)
 
 instance NFData t => NFData (LinearForm t) where
-	rnf = rnf . toList
+  rnf = rnf . toList
 
-instance (Num t, Eq t, Show t) => Show (LinearForm t) where
-	show (LinearForm a b c) = if (a==0) && (b==0) && (c==0)
-		then "0"
-		else "(" ++ intercalate " + " (filter (/=[]) $
-			[if a/= 0 then show a ++ "k" else []] ++
-			[if b/= 0 then show b ++ "l" else []] ++
-			[if c/= 0 then show c ++ "m" else []] ) ++ ")" -- where
-			-- show' :: Rational -> String
-			-- show' z = if denominator z==1 then show (numerator z) else show z
+instance (Num t, Eq t, Pretty t) => Pretty (LinearForm t) where
+  pretty (LinearForm 0 0 0) = char '0'
+  pretty (LinearForm a b c) = cat $ punctuate plus $ mapMaybe f [(a, 'k'), (b, 'l'), (c, 'm')] where
+    plus = space <> char '+' <> space
+    f (0, _) = Nothing
+    f (1, t) = Just (char t)
+    f (r, t) = Just (pretty r <+> char '*' <+> char t)
 
 instance Num t => Num (LinearForm t) where
-	(LinearForm a b c) + (LinearForm d e f) = LinearForm (a+d) (b+e) (c+f)
-	(*) = error "Multiplication of LinearForm is undefined"
-	negate = fmap negate
-	abs = error "Absolute value of LinearForm is undefined"
-	signum = error "Signum of LinearForm is undefined"
-	fromInteger n = LinearForm 0 0 (fromInteger n)
+  (LinearForm a b c) + (LinearForm d e f) = LinearForm (a+d) (b+e) (c+f)
+  (*)    = error "Multiplication of LinearForm is undefined"
+  negate = fmap negate
+  abs    = error "Absolute value of LinearForm is undefined"
+  signum = error "Signum of LinearForm is undefined"
+  fromInteger n = LinearForm 0 0 (fromInteger n)
 
 instance Num t => Monoid (LinearForm t) where
-	mempty  = 0
-	mappend = (+)
+  mempty  = 0
+  mappend = (+)
 
 scaleLF :: (Num t, Eq t) => t -> LinearForm t -> LinearForm t
 scaleLF 0 = const 0
@@ -68,6 +69,7 @@
 -- |Evaluate a linear form a*k + b*l + c*m for given k, l and m.
 evalLF :: Num t => (t, t, t) -> LinearForm t -> t
 evalLF (k, l, m) (LinearForm a b c) = a * k + l * b + m * c
+{-# INLINE evalLF #-}
 
 -- |Substitute linear forms k, l and m into a given linear form
 -- a*k + b*l + c*m to obtain a new linear form.
@@ -76,55 +78,67 @@
 
 -- | Define a rational form of two variables, equal to the ratio of two 'LinearForm'.
 data RationalForm t = RationalForm (LinearForm t) (LinearForm t)
-	deriving (Eq, Show, Functor, Foldable, Generic)
+  deriving (Eq, Show, Functor, Foldable, Generic)
 
+instance (Num t, Eq t, Pretty t) => Pretty (RationalForm t) where
+  pretty (RationalForm l1 l2) = parens (pretty l1) </> parens (pretty l2)
+
 instance NFData t => NFData (RationalForm t) where
-	rnf = rnf . toList
+  rnf = rnf . toList
 
 instance Num t => Num (RationalForm t) where
-	(+) = error "Addition of RationalForm is undefined"
-	(*) = error "Multiplication of RationalForm is undefined"
-	negate (RationalForm a b) = RationalForm (negate a) b
-	abs = error "Absolute value of RationalForm is undefined"
-	signum = error "Signum of RationalForm is undefined"
-	fromInteger n = RationalForm (fromInteger n) 1
+  (+) = error "Addition of RationalForm is undefined"
+  (*) = error "Multiplication of RationalForm is undefined"
+  negate (RationalForm a b) = RationalForm (negate a) b
+  abs = error "Absolute value of RationalForm is undefined"
+  signum = error "Signum of RationalForm is undefined"
+  fromInteger n = RationalForm (fromInteger n) 1
 
 instance Num t => Fractional (RationalForm t) where
-	fromRational r = RationalForm (fromInteger $ numerator r) (fromInteger $ denominator r)
-	recip (RationalForm a b) = RationalForm b a
+  fromRational r = RationalForm (fromInteger $ numerator r) (fromInteger $ denominator r)
+  recip (RationalForm a b) = RationalForm b a
 
 mapTriple :: (a -> b) -> (a, a, a) -> (b, b, b)
 mapTriple f (x, y, z) = (f x, f y, f z)
+{-# INLINE mapTriple #-}
 
 -- |Evaluate a rational form (a*k + b*l + c*m) \/ (a'*k + b'*l + c'*m)
 -- for given k, l and m.
 evalRF :: (Real t, Num t) => (Integer, Integer, Integer) -> RationalForm t -> RationalInf
 evalRF (k, l, m) (RationalForm num den) = if denom==0 then InfPlus else Finite (numer / denom) where
-	klm = mapTriple fromInteger (k, l, m)
-	numer = toRational $ evalLF klm num
-	denom = toRational $ evalLF klm den
+  klm = mapTriple fromInteger (k, l, m)
+  numer = toRational $ evalLF klm num
+  denom = toRational $ evalLF klm den
 
 -- |Constants to specify the strictness of 'Constraint'.
 data IneqType
-	-- | Strict inequality (>0).
-	= Strict
-	-- | Non-strict inequality (≥0).
-	| NonStrict
-	deriving (Eq, Ord, Show, Enum, Bounded)
+  -- | Strict inequality (>0).
+  = Strict
+  -- | Non-strict inequality (≥0).
+  | NonStrict
+  deriving (Eq, Ord, Show, Enum, Bounded, Generic)
 
+instance Pretty IneqType where
+  pretty Strict    = text ">"
+  pretty NonStrict = text ">="
+
 -- |A linear constraint of two variables.
-data Constraint t = Constraint (LinearForm t) IneqType
-	deriving (Eq, Show, Functor, Foldable, Generic)
+data Constraint t = Constraint !(LinearForm t) !IneqType
+  deriving (Eq, Show, Functor, Foldable, Generic)
 
+instance (Num t, Eq t, Pretty t) => Pretty (Constraint t) where
+  pretty (Constraint lf ineq) = pretty lf <+> pretty ineq <+> int 0
+
 instance NFData t => NFData (Constraint t) where
-	rnf (Constraint l i) = i `seq` rnf l
+  rnf (Constraint l i) = i `seq` rnf l
 
 -- |Evaluate a rational form of constraint and compare
 -- its value with 0. Strictness depends on the given 'IneqType'.
-checkConstraint :: (Num t, Eq t) => (Integer, Integer, Integer) -> Constraint t -> Bool
-checkConstraint (k, l, m) (Constraint lf ineq)
-	= if ineq==NonStrict
-		then signum numer /= -1
-		else signum numer == 1 where
-			klm = mapTriple fromInteger (k, l, m)
-			numer = evalLF klm lf
+checkConstraint :: (Num t, Ord t) => (Integer, Integer, Integer) -> Constraint t -> Bool
+checkConstraint (k, l, m) (Constraint lf ineq) = case ineq of
+  NonStrict -> numer >= 0
+  Strict    -> numer >  0
+  where
+    klm   = mapTriple fromInteger (k, l, m)
+    numer = evalLF klm lf
+{-# SPECIALIZE checkConstraint :: (Integer, Integer, Integer) -> Constraint Rational -> Bool #-}
diff --git a/Math/ExpPairs/Matrix3.hs b/Math/ExpPairs/Matrix3.hs
--- a/Math/ExpPairs/Matrix3.hs
+++ b/Math/ExpPairs/Matrix3.hs
@@ -1,4 +1,4 @@
-{-# LANGUAGE BangPatterns, RecordWildCards, DeriveFunctor, DeriveFoldable, DeriveGeneric #-}
+{-# LANGUAGE RecordWildCards, DeriveFunctor, DeriveFoldable, DeriveGeneric #-}
 {-|
 Module      : Math.ExpPairs.Matrix3
 Description : Implements matrices of order 3
@@ -12,33 +12,23 @@
 Can be used instead of "Data.Matrix" to reduce overhead and simplify code.
 -}
 module Math.ExpPairs.Matrix3
-	( Matrix3 (..)
-	, Vector3 (..)
-	, fromList
-	, toList
-	, det
-	, multCol
-	, normalize
-	, makarovMult
-	, ladermanMult
-	) where
+  ( Matrix3 (..)
+  , fromList
+  , toList
+  , det
+  , multCol
+  , normalize
+  , makarovMult
+  , ladermanMult
+  ) where
 
 import Prelude hiding (foldl1)
+
+import Control.DeepSeq
 import Data.Foldable  (Foldable (..), toList)
-import GHC.Generics   (Generic (..))
 import Data.List      (transpose)
-import Control.DeepSeq
-
--- |Three-component vector.
-data Vector3 t = Vector3 {
-	a1 :: !t,
-	a2 :: !t,
-	a3 :: !t
-	}
-	deriving (Eq, Show, Functor, Foldable, Generic)
-
-instance NFData t => NFData (Vector3 t) where
-	rnf = rnf . toList
+import GHC.Generics   (Generic (..))
+import Text.PrettyPrint.Leijen
 
 -- |Matrix of order 3. Instances of 'Num' and 'Fractional'
 -- are given in terms of the multiplicative group of matrices,
@@ -48,69 +38,69 @@
 -- > toList 1 /= [1,1,1,1,1,1,1,1,1]
 --
 data Matrix3 t = Matrix3 {
-	a11 :: !t,
-	a12 :: !t,
-	a13 :: !t,
-	a21 :: !t,
-	a22 :: !t,
-	a23 :: !t,
-	a31 :: !t,
-	a32 :: !t,
-	a33 :: !t
-	}
-	deriving (Eq, Functor, Foldable, Generic)
+  a11 :: !t,
+  a12 :: !t,
+  a13 :: !t,
+  a21 :: !t,
+  a22 :: !t,
+  a23 :: !t,
+  a31 :: !t,
+  a32 :: !t,
+  a33 :: !t
+  }
+  deriving (Eq, Show, Functor, Foldable, Generic)
 
 instance NFData t => NFData (Matrix3 t) where
-	rnf = rnf . toList
+  rnf = rnf . toList
 
 diag :: Num t => t -> Matrix3 t
 diag n = Matrix3 {
-	a11 = n,
-	a12 = 0,
-	a13 = 0,
-	a21 = 0,
-	a22 = n,
-	a23 = 0,
-	a31 = 0,
-	a32 = 0,
-	a33 = n
-	}
+  a11 = n,
+  a12 = 0,
+  a13 = 0,
+  a21 = 0,
+  a22 = n,
+  a23 = 0,
+  a31 = 0,
+  a32 = 0,
+  a33 = n
+  }
 
 instance (Num t, Ord t) => Num (Matrix3 t) where
-	a + b = Matrix3 {
-		a11 = a11 a + a11 b,
-		a12 = a12 a + a12 b,
-		a13 = a13 a + a13 b,
-		a21 = a21 a + a21 b,
-		a22 = a22 a + a22 b,
-		a23 = a23 a + a23 b,
-		a31 = a31 a + a31 b,
-		a32 = a32 a + a32 b,
-		a33 = a33 a + a33 b
-		}
+  a + b = Matrix3 {
+    a11 = a11 a + a11 b,
+    a12 = a12 a + a12 b,
+    a13 = a13 a + a13 b,
+    a21 = a21 a + a21 b,
+    a22 = a22 a + a22 b,
+    a23 = a23 a + a23 b,
+    a31 = a31 a + a31 b,
+    a32 = a32 a + a32 b,
+    a33 = a33 a + a33 b
+    }
 
-	(*) = usualMult
+  (*) = usualMult
 
-	negate = fmap negate
+  negate = fmap negate
 
-	abs = undefined
+  abs = undefined
 
-	signum = diag . signum . det
+  signum = diag . signum . det
 
-	fromInteger = diag . fromInteger
+  fromInteger = diag . fromInteger
 
 usualMult :: Num t => Matrix3 t -> Matrix3 t -> Matrix3 t
 usualMult a b = Matrix3 {
-		a11 = a11 a * a11 b + a12 a * a21 b + a13 a * a31 b,
-		a12 = a11 a * a12 b + a12 a * a22 b + a13 a * a32 b,
-		a13 = a11 a * a13 b + a12 a * a23 b + a13 a * a33 b,
-		a21 = a21 a * a11 b + a22 a * a21 b + a23 a * a31 b,
-		a22 = a21 a * a12 b + a22 a * a22 b + a23 a * a32 b,
-		a23 = a21 a * a13 b + a22 a * a23 b + a23 a * a33 b,
-		a31 = a31 a * a11 b + a32 a * a21 b + a33 a * a31 b,
-		a32 = a31 a * a12 b + a32 a * a22 b + a33 a * a32 b,
-		a33 = a31 a * a13 b + a32 a * a23 b + a33 a * a33 b
-		}
+    a11 = a11 a * a11 b + a12 a * a21 b + a13 a * a31 b,
+    a12 = a11 a * a12 b + a12 a * a22 b + a13 a * a32 b,
+    a13 = a11 a * a13 b + a12 a * a23 b + a13 a * a33 b,
+    a21 = a21 a * a11 b + a22 a * a21 b + a23 a * a31 b,
+    a22 = a21 a * a12 b + a22 a * a22 b + a23 a * a32 b,
+    a23 = a21 a * a13 b + a22 a * a23 b + a23 a * a33 b,
+    a31 = a31 a * a11 b + a32 a * a21 b + a33 a * a31 b,
+    a32 = a31 a * a12 b + a32 a * a22 b + a33 a * a32 b,
+    a33 = a31 a * a13 b + a32 a * a23 b + a33 a * a33 b
+    }
 {-# SPECIALIZE usualMult :: Matrix3 Int -> Matrix3 Int -> Matrix3 Int #-}
 {-# SPECIALIZE usualMult :: Matrix3 Integer -> Matrix3 Integer -> Matrix3 Integer #-}
 
@@ -124,63 +114,63 @@
 -- We were able to reduce the number of additions from 98 to 68 by sofisticated choice of intermediate variables.
 ladermanMult :: Num t => Matrix3 t -> Matrix3 t -> Matrix3 t
 ladermanMult
-	(Matrix3 a11 a12 a13 a21 a22 a23 a31 a32 a33)
-	(Matrix3 b11 b12 b13 b21 b22 b23 b31 b32 b33)
-	= Matrix3 c11 c12 c13 c21 c22 c23 c31 c32 c33 where
-		t33 = t37 + a12 - a32
-		t34 = a13 - a23
-		t35 = a13 - a33
-		t36 = a31 - a11
-		t37 = a11 - a22
+  (Matrix3 a11 a12 a13 a21 a22 a23 a31 a32 a33)
+  (Matrix3 b11 b12 b13 b21 b22 b23 b31 b32 b33)
+  = Matrix3 c11 c12 c13 c21 c22 c23 c31 c32 c33 where
+    t33 = t37 + a12 - a32
+    t34 = a13 - a23
+    t35 = a13 - a33
+    t36 = a31 - a11
+    t37 = a11 - a22
 
-		u33 = b21 - b11 - b23 - b31
-		u34 = b22 - b12
-		u35 = b22 - b32
-		u36 = b33 - b31
-		u37 = b13 - b23
+    u33 = b21 - b11 - b23 - b31
+    u34 = b22 - b12
+    u35 = b22 - b32
+    u36 = b33 - b31
+    u37 = b13 - b23
 
-		m1 = (t35 + t33 - a21) * b22
-		m2 = (a11 - a21) * u34
-		m3 = a22 * (u33 + b33 - u34)
-		m4 = (a21 - t37) * (b11 + u34)
-		m5 = (a22 + a21) * (b12 - b11)
-		m6 = a11 * b11
-		m7 = (t36 + a32) * (b11 - u37)
-		m8 = t36 * u37
-		m9 = (a31 + a32) * (b13 - b11)
-		m10 = (t33 - a31 + t34) * b23
-		m11 = a32 * (u33 + b13 - u35)
-		m12 = (a32 - t35) * (b31 + u35)
-		m13 = t35 * u35
-		m14 = a13 * b31
-		m15 = (a33 + a32) * (b32 - b31)
-		m16 = (a22 - t34) * (b23 - u36)
-		m17 = t34 * (b23 - b33)
-		m18 = (a23 + a22) * u36
-		m19 = a12 * b21
-		m20 = a23 * b32
-		m21 = a21 * b13
-		m22 = a31 * b12
-		m23 = a33 * b33
+    m1 = (t35 + t33 - a21) * b22
+    m2 = (a11 - a21) * u34
+    m3 = a22 * (u33 + b33 - u34)
+    m4 = (a21 - t37) * (b11 + u34)
+    m5 = (a22 + a21) * (b12 - b11)
+    m6 = a11 * b11
+    m7 = (t36 + a32) * (b11 - u37)
+    m8 = t36 * u37
+    m9 = (a31 + a32) * (b13 - b11)
+    m10 = (t33 - a31 + t34) * b23
+    m11 = a32 * (u33 + b13 - u35)
+    m12 = (a32 - t35) * (b31 + u35)
+    m13 = t35 * u35
+    m14 = a13 * b31
+    m15 = (a33 + a32) * (b32 - b31)
+    m16 = (a22 - t34) * (b23 - u36)
+    m17 = t34 * (b23 - b33)
+    m18 = (a23 + a22) * u36
+    m19 = a12 * b21
+    m20 = a23 * b32
+    m21 = a21 * b13
+    m22 = a31 * b12
+    m23 = a33 * b33
 
-		v33 = m12 + m14
-		v34 = m16 + m14
-		v35 = m4 + m6
-		v36 = m7 + m6
-		v37 = v33 + m15
-		v38 = v35 + m5
-		v39 = v34 + m18
-		v40 = v36 + m9
+    v33 = m12 + m14
+    v34 = m16 + m14
+    v35 = m4 + m6
+    v36 = m7 + m6
+    v37 = v33 + m15
+    v38 = v35 + m5
+    v39 = v34 + m18
+    v40 = v36 + m9
 
-		c11 = m6 + m19 + m14
-		c12 = v38 + v37 + m1
-		c13 = v40 + v39 + m10
-		c21 = v35 + v34 + m3 + m2 + m17
-		c22 = v38 + m20 + m2
-		c23 = v39 + m21 + m17
-		c31 = v36 + v33 + m8 + m13 + m11
-		c32 = v37 + m22 + m13
-		c33 = v40 + m23 + m8
+    c11 = m6 + m19 + m14
+    c12 = v38 + v37 + m1
+    c13 = v40 + v39 + m10
+    c21 = v35 + v34 + m3 + m2 + m17
+    c22 = v38 + m20 + m2
+    c23 = v39 + m21 + m17
+    c31 = v36 + v33 + m8 + m13 + m11
+    c32 = v37 + m22 + m13
+    c33 = v40 + m23 + m8
 {-# SPECIALIZE ladermanMult :: Matrix3 Integer -> Matrix3 Integer -> Matrix3 Integer #-}
 
 -- | Multiplicate matrices under assumption that multiplication of elements is commutative.
@@ -194,97 +184,97 @@
 -- We were able to reduce the number of additions from 105 to 66 by sofisticated choice of intermediate variables.
 makarovMult :: Num t => Matrix3 t -> Matrix3 t -> Matrix3 t
 makarovMult
-	(Matrix3 k1 b1 c1 k2 b2 c2 k3 b3 c3)
-	(Matrix3 a1 a2 a3 k4 k5 k6 k7 k8 k9)
-	= Matrix3 c11 c12 c13 c21 c22 c23 c31 c32 c33 where
-		t32 = c3 + c2
-		t33 = b3 + b1
-		t34 = c1 - c2
-		t35 = b2 + b1
+  (Matrix3 k1 b1 c1 k2 b2 c2 k3 b3 c3)
+  (Matrix3 a1 a2 a3 k4 k5 k6 k7 k8 k9)
+  = Matrix3 c11 c12 c13 c21 c22 c23 c31 c32 c33 where
+    t32 = c3 + c2
+    t33 = b3 + b1
+    t34 = c1 - c2
+    t35 = b2 + b1
 
-		u32 = k4 + k6 - k5
-		u33 = k9 + k7 - k8
-		u34 = k6 + k8
+    u32 = k4 + k6 - k5
+    u33 = k9 + k7 - k8
+    u34 = k6 + k8
 
-		m1 = (t34 + a3) * (u33 + k1)
-		m2 = (t35 + a2) * (k2 - u32)
-		m3 = (t33 + a2) * (k3 - u32)
-		m4 = (a3 - t32) * (k3 - u33)
-		m5 = (a1 - t34) * k1
-		m6 = (t35 + a1) * k2
-		m7 = (t33 + t32 + a1) * k3
-		m8 = a2 * (k1 + u32)
-		m9 = a3 * (u33 + k2)
-		m10 = b1 * k4
-		m11 = c2 * k7
-		m12 = t34 * (k7 + k1)
-		m13 = t35 * (k4 - k2)
-		m14 = (b1 + a2) * u32
-		m15 = b2 * k6
-		m16 = (a3 - c2) * u33
-		m17 = c2 * k8
-		m18 = (b3 - t32) * k6
-		m19 = (c3 + c1 - t33) * k8
-		m20 = t33 * (u34 + k4 - k3)
-		m21 = t32 * (u34 + k3 - k7)
-		m22 = (t32 - t33) * u34
+    m1 = (t34 + a3) * (u33 + k1)
+    m2 = (t35 + a2) * (k2 - u32)
+    m3 = (t33 + a2) * (k3 - u32)
+    m4 = (a3 - t32) * (k3 - u33)
+    m5 = (a1 - t34) * k1
+    m6 = (t35 + a1) * k2
+    m7 = (t33 + t32 + a1) * k3
+    m8 = a2 * (k1 + u32)
+    m9 = a3 * (u33 + k2)
+    m10 = b1 * k4
+    m11 = c2 * k7
+    m12 = t34 * (k7 + k1)
+    m13 = t35 * (k4 - k2)
+    m14 = (b1 + a2) * u32
+    m15 = b2 * k6
+    m16 = (a3 - c2) * u33
+    m17 = c2 * k8
+    m18 = (b3 - t32) * k6
+    m19 = (c3 + c1 - t33) * k8
+    m20 = t33 * (u34 + k4 - k3)
+    m21 = t32 * (u34 + k3 - k7)
+    m22 = (t32 - t33) * u34
 
-		v32 = v38 - v35
-		v33 = v35 - v36
-		v34 = m19 - m22
-		v35 = m17 - m18
-		v36 = m14 - m10
-		v37 = m11 + m10
-		v38 = m16 + m11
-		v39 = m20 + m22
-		v40 = m15 + m17
+    v32 = v38 - v35
+    v33 = v35 - v36
+    v34 = m19 - m22
+    v35 = m17 - m18
+    v36 = m14 - m10
+    v37 = m11 + m10
+    v38 = m16 + m11
+    v39 = m20 + m22
+    v40 = m15 + m17
 
-		c11 = v37 + m5 + m12
-		c12 = v34 + v33 + m8
-		c13 = v34 + m1 - m12 - v32
-		c21 = m6 + m13 + m11 - m10
-		c22 = v40 + v36 + m2 + m13
-		c23 = v40 + m9 - v38
-		c31 = v39 + m7 - m21 - v37
-		c32 = v39 + m3 - v33
-		c33 = v32 + m4 + m21
+    c11 = v37 + m5 + m12
+    c12 = v34 + v33 + m8
+    c13 = v34 + m1 - m12 - v32
+    c21 = m6 + m13 + m11 - m10
+    c22 = v40 + v36 + m2 + m13
+    c23 = v40 + m9 - v38
+    c31 = v39 + m7 - m21 - v37
+    c32 = v39 + m3 - v33
+    c33 = v32 + m4 + m21
 {-# SPECIALIZE makarovMult :: Matrix3 Integer -> Matrix3 Integer -> Matrix3 Integer #-}
 
 -- |Compute the determinant of a matrix.
 det :: (Num t, Ord t) => Matrix3 t -> t
 det Matrix3 {..} =
-	a11 * (a22 * a33 - a32 * a23)
-	- a12 * (a21 * a33 - a23 * a31)
-	+ a13 * (a21 * a32 - a22 * a31)
+  a11 * (a22 * a33 - a32 * a23)
+  - a12 * (a21 * a33 - a23 * a31)
+  + a13 * (a21 * a32 - a22 * a31)
 
 instance (Fractional t, Ord t) => Fractional (Matrix3 t) where
-	fromRational = diag . fromRational
+  fromRational = diag . fromRational
 
-	recip a@(Matrix3 {..}) = Matrix3 {
-		a11 =  (a22 * a33 - a32 * a23) / d,
-		a12 = -(a21 * a33 - a23 * a31) / d,
-		a13 =  (a21 * a32 - a22 * a31) / d,
-		a21 = -(a12 * a33 - a13 * a32) / d,
-		a22 =  (a11 * a33 - a13 * a31) / d,
-		a23 = -(a11 * a32 - a12 * a31) / d,
-		a31 =  (a12 * a23 - a13 * a22) / d,
-		a32 = -(a11 * a23 - a13 * a21) / d,
-		a33 =  (a11 * a22 - a12 * a21) / d
-		} where d = det a
+  recip a@(Matrix3 {..}) = Matrix3 {
+    a11 =  (a22 * a33 - a32 * a23) / d,
+    a12 = -(a21 * a33 - a23 * a31) / d,
+    a13 =  (a21 * a32 - a22 * a31) / d,
+    a21 = -(a12 * a33 - a13 * a32) / d,
+    a22 =  (a11 * a33 - a13 * a31) / d,
+    a23 = -(a11 * a32 - a12 * a31) / d,
+    a31 =  (a12 * a23 - a13 * a22) / d,
+    a32 = -(a11 * a23 - a13 * a21) / d,
+    a33 =  (a11 * a22 - a12 * a21) / d
+    } where d = det a
 
 -- |Convert a list of 9 elements into 'Matrix3'. Reverse conversion can be done by 'toList' from "Data.Foldable".
 fromList :: [t] -> Matrix3 t
 fromList [a11, a12, a13, a21, a22, a23, a31, a32, a33] = Matrix3 {
-	a11 = a11,
-	a12 = a12,
-	a13 = a13,
-	a21 = a21,
-	a22 = a22,
-	a23 = a23,
-	a31 = a31,
-	a32 = a32,
-	a33 = a33
-	}
+  a11 = a11,
+  a12 = a12,
+  a13 = a13,
+  a21 = a21,
+  a22 = a22,
+  a23 = a23,
+  a31 = a31,
+  a32 = a32,
+  a33 = a33
+  }
 fromList _ = error "The list must contain exactly 9 elements"
 
 -- |Divide all elements of the matrix by their greatest common
@@ -292,20 +282,20 @@
 -- transformations to reduce the magnitude of computations.
 normalize :: Integral t => Matrix3 t -> Matrix3 t
 normalize a = case foldl1 gcd a of
-	0 -> a
-	d -> fmap (`div` d) a
+  0 -> a
+  d -> fmap (`div` d) a
 
-instance Show t => Show (Matrix3 t) where
-	show = unlines . map unwords . pad . fmap show where
-		pad (Matrix3 {..}) = map (zipWith padCell ls) table where
-			table = [[a11, a12, a13], [a21, a22, a23], [a31, a32, a33]]
-			ls = map (maximum . map length) (transpose table)
-			padCell l xs = replicate (l - length xs) ' ' ++ xs
+instance Pretty t => Pretty (Matrix3 t) where
+  pretty = vsep . map hsep . pad . fmap pretty where
+    pad (Matrix3 {..}) = map (zipWith fill ls) table where
+      table = [[a11, a12, a13], [a21, a22, a23], [a31, a32, a33]]
+      ls = map (maximum . map (length . show)) (transpose table)
 
 -- |Multiplicate a matrix by a vector (considered as a column).
-multCol :: Num t => Matrix3 t -> Vector3 t -> Vector3 t
-multCol Matrix3 {..} Vector3 {..} = Vector3 {
-	a1 = a11 * a1 + a12 * a2 + a13 * a3,
-	a2 = a21 * a1 + a22 * a2 + a23 * a3,
-	a3 = a31 * a1 + a32 * a2 + a33 * a3
-	}
+multCol :: Num t => Matrix3 t -> (t, t, t) -> (t, t, t)
+multCol Matrix3 {..} (a1, a2, a3) = (
+  a11 * a1 + a12 * a2 + a13 * a3,
+  a21 * a1 + a22 * a2 + a23 * a3,
+  a31 * a1 + a32 * a2 + a33 * a3
+  )
+{-# INLINE multCol #-}
diff --git a/Math/ExpPairs/MenzerNowak.hs b/Math/ExpPairs/MenzerNowak.hs
--- a/Math/ExpPairs/MenzerNowak.hs
+++ b/Math/ExpPairs/MenzerNowak.hs
@@ -18,8 +18,8 @@
 
 -}
 module Math.ExpPairs.MenzerNowak
-	( menzerNowak
-	) where
+  ( menzerNowak
+  ) where
 
 import Data.Ratio    ((%))
 
@@ -28,10 +28,10 @@
 -- |Compute Θ(a, b) for given a and b.
 menzerNowak :: Integer -> Integer -> OptimizeResult
 menzerNowak a' b' = optimize
-	[
-		RationalForm (LinearForm 1 1 0) (LinearForm (a+b) 0 (a+b)),
-		RationalForm (LinearForm 1 0 0) (LinearForm (a+b) (-a) a)
-	]
-	[] where
-		a = a'%1
-		b = b'%1
+  [
+    RationalForm (LinearForm 1 1 0) (LinearForm (a+b) 0 (a+b)),
+    RationalForm (LinearForm 1 0 0) (LinearForm (a+b) (-a) a)
+  ]
+  [] where
+    a = a'%1
+    b = b'%1
diff --git a/Math/ExpPairs/Pair.hs b/Math/ExpPairs/Pair.hs
--- a/Math/ExpPairs/Pair.hs
+++ b/Math/ExpPairs/Pair.hs
@@ -12,63 +12,75 @@
 
 Below /A/ and /B/ stands for van der Corput's processes. See "Math.ExpPairs.Process" for explanations.
 -}
+{-# LANGUAGE DeriveGeneric        #-}
+{-# LANGUAGE FlexibleInstances    #-}
+{-# LANGUAGE TypeSynonymInstances #-}
+
+{-# OPTIONS_GHC -fno-warn-orphans #-}
+
 module Math.ExpPairs.Pair
-	( Triangle (..)
-	, InitPair' (..)
-	, InitPair
-	, initPairs
-	, initPairToValue
-	) where
+  ( Triangle (..)
+  , InitPair' (..)
+  , InitPair
+  , initPairs
+  , initPairToValue
+  , initPairToProjValue
+  ) where
 
-import Data.Ratio ((%))
+import Data.Maybe
+import Data.Ratio
+import GHC.Generics (Generic (..))
+import Text.PrettyPrint.Leijen
 
 -- |Vertices of the triangle of initial exponent pairs.
 data Triangle
-	-- |Usual van der Corput exponent pair
-	-- (1\/6, 2\/3) = /AB/(0, 1).
-	= Corput16
-	-- |An exponent pair (2\/13, 35\/52) from /Huxley M. N./
-	-- `Exponential sums and the Riemann zeta function'
-	-- \/\/ Proceedings of the International Number
-	-- Theory Conference held at Universite Laval in 1987, Walter de Gruyter, 1989, P. 417-423.
-	| HuxW87b1
-	-- | An exponent pair (32\/205, 269\/410) from /Huxley M. N./
-	-- `Exponential sums and the Riemann zeta function V' \/\/
+  -- |Usual van der Corput exponent pair
+  -- (1\/6, 2\/3) = /AB/(0, 1).
+  = Corput16
+  -- |An exponent pair (2\/13, 35\/52) from /Huxley M. N./
+  -- `Exponential sums and the Riemann zeta function'
+  -- \/\/ Proceedings of the International Number
+  -- Theory Conference held at Universite Laval in 1987, Walter de Gruyter, 1989, P. 417-423.
+  | HuxW87b1
+  -- | An exponent pair (32\/205, 269\/410) from /Huxley M. N./
+  -- `Exponential sums and the Riemann zeta function V' \/\/
   -- Proc. Lond. Math. Soc., 2005, Vol. 90, no. 1., P. 1--41.
-	| Hux05
-	deriving (Show, Bounded, Enum, Eq, Ord)
+  | Hux05
+  deriving (Show, Bounded, Enum, Eq, Ord, Generic)
 
+instance Pretty Triangle where
+  pretty = text . show
+
 -- |Type to hold an initial exponent pair.
 data InitPair' t
-	-- |Usual van der Corput exponent pair
-	-- (0, 1).
-	= Corput01
-	-- |Usual van der Corput exponent pair
-	-- (1\/2, 1\/2) = /B/(0, 1).
-	| Corput12
-	-- |Point from the interior of 'Triangle'.
-	-- Exactly
-	-- 'Mix' a b = a * 'Corput16' + b * 'HuxW87b1' + (1-a-b) * 'Hux05'
-	| Mix t t
-	deriving (Eq)
+  -- |Usual van der Corput exponent pair
+  -- (0, 1).
+  = Corput01
+  -- |Usual van der Corput exponent pair
+  -- (1\/2, 1\/2) = /B/(0, 1).
+  | Corput12
+  -- |Point from the interior of 'Triangle'.
+  -- Exactly
+  -- 'Mix' a b = a * 'Corput16' + b * 'HuxW87b1' + (1-a-b) * 'Hux05'
+  | Mix !t !t
+  deriving (Eq, Show, Generic)
 
 -- |Exponent pair built from rational fractions of
 -- 'Corput16', 'HuxW87b1' and 'Hux05'
 type InitPair = InitPair' Rational
 
-instance (Show t, Num t, Eq t) => Show (InitPair' t) where
-	show Corput01 = "(0, 1)"
-	show Corput12 = "(1/2, 1/2)"
-	show (Mix r1 r2) =
-		s1 ++ (if s1/="" && (s2/=""||s3/="") then " + " else "")
-		++ s2 ++ (if s2/="" && s3/="" then " + " else "") ++ s3
-		where
-			r3 = 1 - r1 - r2
-			f r t = if r==0 then "" else (if r==1 then "" else show r ++ " * ") ++ show t
-			s1 = f r1 Corput16
-			s2 = f r2 HuxW87b1
-			s3 = f r3 Hux05
+instance Pretty Rational where
+  pretty = rational
 
+instance (Pretty t, Num t, Eq t) => Pretty (InitPair' t) where
+  pretty Corput01 = parens (rational 0     <> comma <+> rational 1)
+  pretty Corput12 = parens (rational (1%2) <> comma <+> rational (1%2))
+  pretty (Mix r1 r2) = cat $ punctuate plus $ mapMaybe f [(r1, Corput16), (r2, HuxW87b1), (1 - r1 - r2, Hux05)] where
+    plus = space <> char '+' <> space
+    f (0, _) = Nothing
+    f (1, t) = Just (pretty t)
+    f (r, t) = Just (pretty r <+> char '*' <+> pretty t)
+
 sect :: Integer
 sect = 30
 
@@ -76,18 +88,35 @@
 -- 'Corput01', 'Corput12' and 496 = sum [1..31] 'Mix'-points,
 -- which forms a uniform net over 'Triangle'.
 initPairs :: [InitPair]
-initPairs = Corput01 : Corput12 : [Mix (r1%sect) (r2%sect) | r1<-[0..sect], r2<-[0..sect-r1]]
+initPairs = Corput01 : Corput12 : [Mix (r1 % sect) (r2 % sect) | r1 <- [0 .. sect], r2 <- [0 .. sect - r1]]
 
 -- |Convert initial exponent pair from its symbolic representation
 -- as 'InitPair' to pair of rationals.
 initPairToValue :: InitPair -> (Rational, Rational)
+initPairToValue (Mix r1 r2) = (x, y) where
+  r3 = 1 - r1 - r2
+  (x1, y1) = (1%6, 2%3)
+  (x2, y2) = ( 2 %  13,  35 %  52)
+  (x3, y3) = (32 % 205, 269 % 410)
+  x = x1*r1 + x2*r2 + x3*r3
+  y = y1*r1 + y2*r2 + y3*r3
+--initPairToValue (Mix r1 r2) = (13 % 1230 * r1 - 6 % 2665 * r2 + 32 % 205, 13 % 1230 * r1 + 181 % 10660 * r2 + 269 % 410)
 initPairToValue Corput01 = (0, 1)
 initPairToValue Corput12 = (1%2, 1%2)
-initPairToValue (Mix r1 r2) = (x, y) where
-	r3 = 1 - r1 - r2
-	(x1, y1) = (1%6, 2%3)
-	(x2, y2) = ( 2 %  13,  35 %  52)
-	(x3, y3) = (32 % 205, 269 % 410)
-	x = x1*r1 + x2*r2 + x3*r3
-	y = y1*r1 + y2*r2 + y3*r3
 
+-- | Same as 'initPairToValue', but immediately convert from Q^2 to PN^3.
+initPairToProjValue :: InitPair -> (Integer, Integer, Integer)
+initPairToProjValue (Mix r1 r2) = (k `div` d , l `div` d, m `div` d)
+  where
+    dr1 = denominator r1
+    dr2 = denominator r2
+    m = 31980 * dr1 * dr2
+    k = 338 * numerator r1 * dr2 -  72 * numerator r2 * dr1 +  4992 * dr1 * dr2
+    l = 338 * numerator r1 * dr2 + 543 * numerator r2 * dr1 + 20982 * dr1 * dr2
+
+    d = k `gcd` l `gcd` m
+
+initPairToProjValue Corput01 = (0, 1, 1)
+initPairToProjValue Corput12 = (1, 1, 2)
+
+{-# INLINABLE initPairToProjValue #-}
diff --git a/Math/ExpPairs/PrettyProcess.hs b/Math/ExpPairs/PrettyProcess.hs
--- a/Math/ExpPairs/PrettyProcess.hs
+++ b/Math/ExpPairs/PrettyProcess.hs
@@ -16,9 +16,9 @@
 {-# LANGUAGE TemplateHaskell #-}
 {-# OPTIONS_GHC -fno-warn-type-defaults #-}
 module Math.ExpPairs.PrettyProcess
-	( prettify,
-		uglify,
-		PrettyProcess) where
+  ( prettify,
+    uglify,
+    PrettyProcess) where
 
 import Data.List                (minimumBy)
 import Data.Ord                 (comparing)
@@ -29,23 +29,23 @@
 
 -- | Compact representation of the sequence of 'Process'.
 data PrettyProcess
-	= Simply [Process]
-	| Repeat PrettyProcess Int
-	| Sequence PrettyProcess PrettyProcess
-	deriving (Show)
+  = Simply [Process]
+  | Repeat PrettyProcess Int
+  | Sequence PrettyProcess PrettyProcess
+  deriving (Show)
 
 data PrettyProcessWithWidth = PPWL { ppwlProcess :: PrettyProcess, ppwlWidth :: Int }
 
 deriveMemoizable ''PrettyProcess
 
 instance Pretty PrettyProcess where
-	pretty = \case
-		Simply xs    -> hsep (map (text . show) xs)
-		Repeat _  0  -> empty
-		Repeat xs 1  -> pretty xs
-		Repeat (Simply [A]) n -> text (show A) <> char '^' <> int n
-		Repeat xs n  -> parens (pretty xs) <> char '^' <> int n
-		Sequence a b -> pretty a <+> pretty b
+  pretty = \case
+    Simply xs    -> hsep (map (text . show) xs)
+    Repeat _  0  -> empty
+    Repeat xs 1  -> pretty xs
+    Repeat (Simply [A]) n -> text (show A) <> char '^' <> int n
+    Repeat xs n  -> parens (pretty xs) <> char '^' <> int n
+    Sequence a b -> pretty a <+> pretty b
 
 -- | Width of the bracket.
 bracketWidth :: Int
@@ -63,12 +63,12 @@
 -- | Compute the width of the 'PrettyProcess' according to 'bracketWidth', 'subscriptWidth' and 'printedWidth''.
 printedWidth :: PrettyProcess -> Int
 printedWidth = \case
-	Simply xs             -> sum (map processWidth xs)
-	Repeat _ 0            -> 0
-	Repeat xs 1           -> printedWidth xs
-	Repeat (Simply [A]) _ -> processWidth A + subscriptWidth
-	Repeat xs _           -> printedWidth xs + bracketWidth * 2 + subscriptWidth
-	Sequence a b          -> printedWidth a + printedWidth b
+  Simply xs             -> sum (map processWidth xs)
+  Repeat _ 0            -> 0
+  Repeat xs 1           -> printedWidth xs
+  Repeat (Simply [A]) _ -> processWidth A + subscriptWidth
+  Repeat xs _           -> printedWidth xs + bracketWidth * 2 + subscriptWidth
+  Sequence a b          -> printedWidth a + printedWidth b
 
 -- | Convert 'PrettyProcess' to 'PrettyProcessWithWidth'.
 annotateWithWidth :: PrettyProcess -> PrettyProcessWithWidth
@@ -77,16 +77,16 @@
 -- | Return non-trivial divisors of an argument.
 divisors :: Int -> [Int]
 divisors n = ds1 ++ reverse ds2 where
-	(ds1, ds2) = unzip [ (a, n `div` a) | a <- [1 .. sqrtint n], n `mod` a == 0 ]
-	sqrtint = round . sqrt . fromIntegral
+  (ds1, ds2) = unzip [ (a, n `div` a) | a <- [1 .. sqrtint n], n `mod` a == 0 ]
+  sqrtint = round . sqrt . fromIntegral
 
 -- | Try to represent list as a replication of list.
 asRepeat :: [Process] -> ([Process], Int)
 asRepeat [] = ([], 0)
 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
+  l = length xs
+  candidates = [ (take d xs, l `div` d) | d <- divisors l ]
+  pair = head $ filter (\(ys, n) -> concat (replicate n ys) == xs) candidates
 
 -- | Find the most compact representation of the sequence of processes.
 prettify :: [Process] -> PrettyProcess
@@ -98,28 +98,28 @@
 
 prettify' :: [Process] -> PrettyProcessWithWidth
 prettify' = \case
-	[]   -> annotateWithWidth (Simply [])
-	[A]  -> annotateWithWidth (Simply [A])
-	[BA] -> annotateWithWidth (Simply [BA])
-	xs   -> minimumBy (comparing ppwlWidth) yss where
-		xs'' = case asRepeat xs of
-			(_, 1)   -> annotateWithWidth (Simply xs)
-			(xs', n) -> annotateWithWidth (Repeat (prettify xs') n)
+  []   -> annotateWithWidth (Simply [])
+  [A]  -> annotateWithWidth (Simply [A])
+  [BA] -> annotateWithWidth (Simply [BA])
+  xs   -> minimumBy (comparing ppwlWidth) yss where
+    xs'' = case asRepeat xs of
+      (_, 1)   -> annotateWithWidth (Simply xs)
+      (xs', n) -> annotateWithWidth (Repeat (prettify xs') n)
 
-		yss = xs'' : map f bcs
+    yss = xs'' : map f bcs
 
-		bcs = takeWhile (not . null . snd) $ iterate bcf ([head xs], tail xs)
+    bcs = takeWhile (not . null . snd) $ iterate bcf ([head xs], tail xs)
 
-		bcf (_, [])    = undefined
-		bcf (zs, y:ys) = (zs++[y], ys)
+    bcf (_, [])    = undefined
+    bcf (zs, y:ys) = (zs++[y], ys)
 
-		f (bs, cs) = PPWL (Sequence bsP csP) (bsW + csW) where
-			PPWL bsP bsW = prettifyP bs
-			PPWL csP csW = prettifyP cs
+    f (bs, cs) = PPWL (Sequence bsP csP) (bsW + csW) where
+      PPWL bsP bsW = prettifyP bs
+      PPWL csP csW = prettifyP cs
 
 -- | Unfold back 'PrettyProcess' into the sequence of 'Process'.
 uglify :: PrettyProcess -> [Process]
 uglify = \case
-	Simply xs      -> xs
-	Repeat xs n    -> concat . replicate n . uglify $ xs
-	Sequence xs ys -> uglify xs ++ uglify ys
+  Simply xs      -> xs
+  Repeat xs n    -> concat . replicate n . uglify $ xs
+  Sequence xs ys -> uglify xs ++ uglify ys
diff --git a/Math/ExpPairs/Process.hs b/Math/ExpPairs/Process.hs
--- a/Math/ExpPairs/Process.hs
+++ b/Math/ExpPairs/Process.hs
@@ -9,20 +9,19 @@
 
 Provides types for sequences of /A/- and /B/-processes of van der Corput. A good account on this topic can be found in /Graham S. W.,  Kolesnik G. A./ Van Der Corput's Method of Exponential Sums, Cambridge University Press, 1991, especially Ch. 5.
 -}
-{-# LANGUAGE DeriveGeneric #-}
+{-# LANGUAGE DeriveGeneric, CPP #-}
 module Math.ExpPairs.Process
-	( Process ()
-	, Path (Path)
-	, aPath
-	, baPath
-	, evalPath
-	, lengthPath
-	) where
+  ( Process ()
+  , Path (Path)
+  , aPath
+  , baPath
+  , evalPath
+  , lengthPath
+  ) where
 
 import GHC.Generics             (Generic)
-import Generics.Deriving.Monoid (Monoid, mempty, memptydefault, mappend, mappenddefault)
-import Text.PrettyPrint.Leijen
-
+import Data.Monoid
+import Text.PrettyPrint.Leijen hiding ((<>))
 
 import Math.ExpPairs.ProcessMatrix
 import Math.ExpPairs.PrettyProcess
@@ -34,32 +33,32 @@
 -- > show (mconcat $ replicate 10 aPath) == "A^10"
 --
 data Path = Path !ProcessMatrix ![Process]
-	deriving (Eq, Generic)
+  deriving (Eq, Show, Generic)
 
 instance Monoid Path where
-	mempty  = memptydefault
-	mappend = mappenddefault
+  mempty  = Path mempty mempty
+  mappend (Path m1 p1) (Path m2 p2) = Path (m1 <> m2) (p1 <> p2)
 
-instance Show Path where
-	show (Path _ l) = show (pretty (prettify l)) -- ++ "\n" ++ Mx.prettyMatrix m
+instance Pretty Path where
+  pretty (Path _ l) = pretty (prettify l)
 
 instance Read Path where
-	readsPrec _ zs = [reads' zs] where
-		reads' ('A':xs) = (aPath `mappend` path, ys) where
-			(path, ys) = reads' xs
-		reads' ('B':'A':xs) = (baPath `mappend` path, ys) where
-			(path, ys) = reads' xs
-		reads' ('B':xs) = (baPath, xs)
-		reads' xs = (mempty, xs)
+  readsPrec _ zs = [reads' zs] where
+    reads' ('A':xs) = (aPath <> path, ys) where
+      (path, ys) = reads' xs
+    reads' ('B':'A':xs) = (baPath <> path, ys) where
+      (path, ys) = reads' xs
+    reads' ('B':xs) = (baPath, xs)
+    reads' xs = (mempty, xs)
 
 instance Ord Path where
-	(Path _ q1) <= (Path _ q2) = cmp q1 q2 where
-		cmp (A:p1)  (A:p2)  = cmp p1 p2
-		cmp (BA:p1) (BA:p2) = cmp p2 p1
-		cmp (A:_)   (BA:_)  = True
-		cmp (BA:_)  (A:_)   = False
-		cmp []      _       = True
-		cmp _       []      = False
+  (Path _ q1) <= (Path _ q2) = cmp q1 q2 where
+    cmp (A:p1)  (A:p2)  = cmp p1 p2
+    cmp (BA:p1) (BA:p2) = cmp p2 p1
+    cmp (A:_)   (BA:_)  = True
+    cmp (BA:_)  (A:_)   = False
+    cmp []      _       = True
+    cmp _       []      = False
 
 -- | Path consisting of a single process 'A'.
 aPath :: Path
diff --git a/Math/ExpPairs/ProcessMatrix.hs b/Math/ExpPairs/ProcessMatrix.hs
--- a/Math/ExpPairs/ProcessMatrix.hs
+++ b/Math/ExpPairs/ProcessMatrix.hs
@@ -9,36 +9,43 @@
 
 Provides types for sequences of /A/- and /B/-processes of van der Corput. A good account on this topic can be found in /Graham S. W.,  Kolesnik G. A./ Van Der Corput's Method of Exponential Sums, Cambridge University Press, 1991, especially Ch. 5.
 -}
-{-# LANGUAGE TemplateHaskell, BangPatterns, GeneralizedNewtypeDeriving  #-}
+{-# LANGUAGE TemplateHaskell, BangPatterns, GeneralizedNewtypeDeriving, CPP, DeriveGeneric #-}
 module Math.ExpPairs.ProcessMatrix
-	( Process (..)
-	, ProcessMatrix ()
-	, aMatrix
-	, baMatrix
-	, evalMatrix
-	) where
+  ( Process (..)
+  , ProcessMatrix ()
+  , aMatrix
+  , baMatrix
+  , evalMatrix
+  ) where
 
+#if __GLASGOW_HASKELL__ < 710
 import Data.Monoid           (Monoid, mempty, mappend)
+#endif
 import Data.Function.Memoize (deriveMemoizable)
+import GHC.Generics          (Generic (..))
+import Text.PrettyPrint.Leijen
 
 import Math.ExpPairs.Matrix3
 
 -- | Since B^2 = id, B 'Corput16' = 'Corput16', B 'Hux05' = 'Hux05' and B 'HuxW87b1' = ???, the sequence of /A/- and /B/-processes, applied to 'initPairs' can be rewritten as a sequence of 'A' and 'BA'.
 data Process
-	-- | /A/-process
-	= A
-	-- | /BA/-process
-	| BA
-	deriving (Eq, Show, Read, Ord, Enum)
+  -- | /A/-process
+  = A
+  -- | /BA/-process
+  | BA
+  deriving (Eq, Show, Read, Ord, Enum, Generic)
 
+instance Pretty Process where
+  pretty = text . show
+
 deriveMemoizable ''Process
 
 newtype ProcessMatrix = ProcessMatrix (Matrix3 Integer)
-	deriving (Eq, Num, Show)
+  deriving (Eq, Num, Show, Pretty)
 
 instance Monoid ProcessMatrix where
-	mempty = 1
-	mappend (ProcessMatrix a) (ProcessMatrix b) = ProcessMatrix $ normalize $ a * b
+  mempty = 1
+  mappend (ProcessMatrix a) (ProcessMatrix b) = ProcessMatrix $ normalize $ a * b
 
 process2matrix :: Process -> ProcessMatrix
 process2matrix  A = ProcessMatrix $ Matrix3 1 0 0 1 1 1  2 0 2
@@ -55,7 +62,7 @@
 -- |Apply a projective transformation, defined by 'Path',
 -- to a given point in two-dimensional projective space.
 evalMatrix :: Num t => ProcessMatrix -> (t, t, t) -> (t, t, t)
-evalMatrix (ProcessMatrix m) (a,b,c) = (a',b',c') where
-	m' = fmap fromInteger m
-	(Vector3 a' b' c') = multCol m' (Vector3 a b c)
+evalMatrix (ProcessMatrix m) = multCol (fmap fromInteger m)
+{-# INLINABLE evalMatrix #-}
+{-# SPECIALIZE evalMatrix :: ProcessMatrix -> (Integer, Integer, Integer) -> (Integer, Integer, Integer) #-}
 
diff --git a/Math/ExpPairs/RatioInf.hs b/Math/ExpPairs/RatioInf.hs
--- a/Math/ExpPairs/RatioInf.hs
+++ b/Math/ExpPairs/RatioInf.hs
@@ -1,4 +1,3 @@
-{-# LANGUAGE BangPatterns #-}
 {-|
 Module      : Math.ExpPairs.RatioInf
 Description : Rational numbers with infinities
@@ -11,105 +10,119 @@
 Provides types and necessary instances for rational numbers, extended with infinite values. Just use 'RationalInf' instead of 'Rational' from "Data.Ratio".
 -}
 module Math.ExpPairs.RatioInf
-	( RatioInf (..)
-	, RationalInf
-	) where
+  ( RatioInf (..)
+  , RationalInf
+  ) where
 
-import Data.Ratio (Ratio)
+import Data.Ratio (Ratio, numerator, denominator)
+import Text.PrettyPrint.Leijen
 
 -- |Extends a rational type with positive and negative
 -- infinities.
 data RatioInf t
-	-- |Negative infinity
-	= InfMinus
-	-- |Finite value
-	| Finite !(Ratio t)
-	-- |Positive infinity
-	| InfPlus
-	deriving (Ord, Eq)
+  -- |Negative infinity
+  = InfMinus
+  -- |Finite value
+  | Finite !(Ratio t)
+  -- |Positive infinity
+  | InfPlus
+  deriving (Eq, Ord, Show)
 
 -- |Arbitrary-precision rational numbers with positive and negative
 -- infinities.
 type RationalInf = RatioInf Integer
 
-instance (Integral t, Show t) => Show (RatioInf t) where
-	show InfMinus   = "-Inf"
-	show (Finite x) = show x
-	show InfPlus    = "+Inf"
+instance (Integral t, Pretty t) => Pretty (RatioInf t) where
+  pretty InfMinus   = text "-Inf"
+  pretty (Finite x)
+    | denominator x == 1 = pretty (numerator x)
+    | otherwise          = pretty (numerator x) <+> char '/' <+> pretty (denominator x)
+  pretty InfPlus    = text "+Inf"
 
 instance Integral t => Num (RatioInf t) where
-	InfMinus + InfPlus = error "Cannot add up negative and positive infinities"
-	InfPlus + InfMinus = error "Cannot add up negative and positive infinities"
-	InfMinus + _ = InfMinus
-	InfPlus + _  = InfPlus
-	_ + InfMinus = InfMinus
-	_ + InfPlus  = InfPlus
-	(Finite a) + (Finite b) = Finite (a+b)
+  InfMinus + InfPlus = error "Cannot add up negative and positive infinities"
+  InfPlus + InfMinus = error "Cannot add up negative and positive infinities"
+  InfMinus + _ = InfMinus
+  InfPlus + _  = InfPlus
+  _ + InfMinus = InfMinus
+  _ + InfPlus  = InfPlus
+  (Finite a) + (Finite b) = Finite (a+b)
+  {-# SPECIALIZE (+) :: RationalInf -> RationalInf -> RationalInf #-}
 
-	fromInteger = Finite . fromInteger
+  fromInteger = Finite . fromInteger
+  {-# SPECIALIZE fromInteger :: Integer -> RationalInf #-}
 
-	signum InfMinus   = Finite (-1)
-	signum InfPlus    = Finite 1
-	signum (Finite r) = Finite (signum r)
+  signum InfMinus   = Finite (-1)
+  signum InfPlus    = Finite 1
+  signum (Finite r) = Finite (signum r)
+  {-# SPECIALIZE signum :: RationalInf -> RationalInf #-}
 
-	abs InfMinus   = InfPlus
-	abs InfPlus    = InfPlus
-	abs (Finite r) = Finite (abs r)
+  abs InfMinus   = InfPlus
+  abs InfPlus    = InfPlus
+  abs (Finite r) = Finite (abs r)
+  {-# SPECIALIZE abs :: RationalInf -> RationalInf #-}
 
-	negate InfMinus   = InfPlus
-	negate InfPlus    = InfMinus
-	negate (Finite r) = Finite (negate r)
+  negate InfMinus   = InfPlus
+  negate InfPlus    = InfMinus
+  negate (Finite r) = Finite (negate r)
+  {-# SPECIALIZE negate :: RationalInf -> RationalInf #-}
 
-	InfMinus * InfMinus = InfMinus
-	InfMinus * InfPlus  = InfMinus
-	InfMinus * Finite a = case signum a of
-		1  -> InfMinus
-		-1 -> InfPlus
-		_  -> error "Cannot multiply infinity by zero"
+  InfMinus * InfMinus = InfMinus
+  InfMinus * InfPlus  = InfMinus
+  InfMinus * Finite a = case signum a of
+    1  -> InfMinus
+    -1 -> InfPlus
+    _  -> 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"
+  InfPlus * InfMinus = InfMinus
+  InfPlus * InfPlus  = InfPlus
+  InfPlus * Finite a = case signum a of
+    1  -> InfPlus
+    -1 -> InfMinus
+    _  -> error "Cannot multiply infinity by zero"
 
-	Finite a * InfMinus = case signum a of
-		1  -> InfMinus
-		-1 -> InfPlus
-		_  -> error "Cannot multiply infinity by zero"
+  Finite a * InfMinus = case signum a of
+    1  -> InfMinus
+    -1 -> InfPlus
+    _  -> error "Cannot multiply infinity by zero"
 
-	Finite a * InfPlus = case signum a of
-		1  -> InfPlus
-		-1 -> InfMinus
-		_  -> error "Cannot multiply infinity by zero"
+  Finite a * InfPlus = case signum a of
+    1  -> InfPlus
+    -1 -> InfMinus
+    _  -> error "Cannot multiply infinity by zero"
 
-	Finite a * Finite b = Finite (a * b)
+  Finite a * Finite b = Finite (a * b)
 
+  {-# SPECIALIZE (*) :: RationalInf -> RationalInf -> RationalInf #-}
+
 instance Integral t => Fractional (RatioInf t) where
-	fromRational = Finite . fromRational
+  fromRational = Finite . fromRational
+  {-# SPECIALIZE fromRational :: Rational -> RationalInf #-}
 
-	InfMinus / InfMinus = error "Cannot divide infinity by infinity"
-	InfMinus / InfPlus  = error "Cannot divide infinity by infinity"
-	InfMinus / Finite a = case signum a of
-		1  -> InfMinus
-		-1 -> InfPlus
-		_  -> error "Cannot divide infinity by zero"
+  InfMinus / InfMinus = error "Cannot divide infinity by infinity"
+  InfMinus / InfPlus  = error "Cannot divide infinity by infinity"
+  InfMinus / Finite a = case signum a of
+    1  -> InfMinus
+    -1 -> InfPlus
+    _  -> error "Cannot divide infinity by zero"
 
-	InfPlus  / InfMinus = error "Cannot divide infinity by infinity"
-	InfPlus  / InfPlus  = error "Cannot divide infinity by infinity"
-	InfPlus / Finite a  = case signum a of
-		1  -> InfPlus
-		-1 -> InfMinus
-		_  -> error "Cannot divide infinity by zero"
+  InfPlus  / InfMinus = error "Cannot divide infinity by infinity"
+  InfPlus  / InfPlus  = error "Cannot divide infinity by infinity"
+  InfPlus / Finite a  = case signum a of
+    1  -> InfPlus
+    -1 -> InfMinus
+    _  -> error "Cannot divide infinity by zero"
 
-	Finite _ / InfPlus  = Finite 0
-	Finite _ / InfMinus = Finite 0
+  Finite _ / InfPlus  = Finite 0
+  Finite _ / InfMinus = Finite 0
 
-	Finite _ / Finite 0 = error "Cannot divide finite value by zero"
-	Finite a / Finite b = Finite (a/b)
+  Finite _ / Finite 0 = error "Cannot divide finite value by zero"
+  Finite a / Finite b = Finite (a/b)
 
+  {-# SPECIALIZE (/) :: RationalInf -> RationalInf -> RationalInf #-}
+
 instance Integral t => Real (RatioInf t) where
-	toRational (Finite r) = toRational r
-	toRational InfPlus    = error "Cannot map infinity into Rational"
-	toRational InfMinus   = error "Cannot map infinity into Rational"
+  toRational (Finite r) = toRational r
+  toRational InfPlus    = error "Cannot map infinity into Rational"
+  toRational InfMinus   = error "Cannot map infinity into Rational"
+  {-# SPECIALIZE toRational :: RationalInf -> Rational #-}
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.1.3.0
+version:             0.1.4.0
 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,6 +9,7 @@
 maintainer:          andrew.lelechenko@gmail.com
 category:            Math
 build-type:          Simple
+extra-source-files:  CHANGELOG.md
 cabal-version:       >=1.10
 
 source-repository head
@@ -30,15 +31,24 @@
   build-depends:       base >=4 && <5,
                        memoize >=0.1,
                        ghc-prim,
-                       generic-deriving,
                        wl-pprint >=1.2,
                        deepseq >=1.3
   default-language:    Haskell2010
-  ghc-options:         -Wall
+  ghc-options:         -Wall -fno-warn-type-defaults
 
 test-suite tests
   type:                exitcode-stdio-1.0
   main-is:             Tests.hs
+  other-modules:       Etalon,
+                       Instances,
+                       Ivic,
+                       Kratzel,
+                       LinearForm,
+                       Matrix3,
+                       MenzerNowak,
+                       Pair,
+                       PrettyProcess,
+                       RatioInf
   build-depends:       base >=4 && <5,
                        tasty >=0.7,
                        tasty-quickcheck,
@@ -52,4 +62,4 @@
                        random
   hs-source-dirs:      tests
   default-language:    Haskell2010
-  ghc-options:         -Wall
+  ghc-options:         -Wall -fno-warn-type-defaults
diff --git a/tests/Etalon.hs b/tests/Etalon.hs
new file mode 100644
--- /dev/null
+++ b/tests/Etalon.hs
@@ -0,0 +1,25 @@
+module Etalon (testEtalon) where
+
+import System.Random
+import Data.Ord
+import Data.List
+import Test.Tasty.HUnit
+
+unsort :: (RandomGen g) => g -> [x] -> [x]
+unsort g es = map snd . sortBy (comparing fst) $ zip rs es
+  where rs = randoms g :: [Integer]
+
+fetchRandomLines :: Read b => Int -> FilePath -> IO [[b]]
+fetchRandomLines n filename = do
+  etalon <- readFile filename
+  gen <- newStdGen
+  let items = (take n . unsort gen . lines) etalon
+  let tests = map (map read . words) items
+  return tests
+
+testEtalon :: Int -> ([Integer] -> Bool) -> String -> Assertion
+testEtalon n f filename = do
+  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)
diff --git a/tests/Instances.hs b/tests/Instances.hs
new file mode 100644
--- /dev/null
+++ b/tests/Instances.hs
@@ -0,0 +1,165 @@
+{-# OPTIONS_GHC -fno-warn-orphans #-}
+{-# LANGUAGE FlexibleInstances, MultiParamTypeClasses, DeriveGeneric, CPP #-}
+module Instances (Ratio01 (..), Positive (..), Sorted(..)) where
+
+import Test.QuickCheck hiding (Positive)
+import Test.SmallCheck.Series
+import Control.Applicative
+import Control.Monad
+import GHC.Generics          (Generic (..))
+
+import Math.ExpPairs.LinearForm
+import Math.ExpPairs.ProcessMatrix
+import Math.ExpPairs.Pair (InitPair' (..))
+import Math.ExpPairs.Matrix3 as M3 (Matrix3, fromList)
+
+instance Arbitrary a => Arbitrary (LinearForm a) where
+  arbitrary = LinearForm <$> arbitrary <*> arbitrary <*> arbitrary
+  shrink = genericShrink
+
+instance (Monad m, Serial m a) => Serial m (LinearForm a) where
+  series = cons3 LinearForm
+
+instance Arbitrary a => Arbitrary (RationalForm a) where
+  arbitrary = RationalForm <$> arbitrary <*> arbitrary
+  shrink = genericShrink
+
+instance (Monad m, Serial m a) => Serial m (RationalForm a) where
+  series = cons2 RationalForm
+
+instance Arbitrary a => Arbitrary (Constraint a) where
+  arbitrary = Constraint <$> arbitrary <*> arbitrary
+  shrink = genericShrink
+
+instance (Monad m, Serial m a) => Serial m (Constraint a) where
+  series = cons2 Constraint
+
+instance Arbitrary IneqType where
+  arbitrary = f <$> arbitrary where
+    f x = if x then Strict else NonStrict
+  shrink = genericShrink
+
+instance Monad m => Serial m IneqType where
+  series = cons0 Strict \/ cons0 NonStrict
+
+instance Arbitrary Process where
+  arbitrary = f <$> arbitrary where
+    f x = if x then A else BA
+  shrink = genericShrink
+
+instance Monad m => Serial m Process where
+  series = cons0 A \/ cons0 BA
+
+newtype Ratio01 t = Ratio01 t
+  deriving (Eq, Ord, Generic)
+
+instance (Ord t, Fractional t, Arbitrary t) => Arbitrary (Ratio01 t) where
+  arbitrary = Ratio01 <$> (arbitrary `suchThat` (\x -> 0 <= x && x <= 1))
+  shrink = genericShrink
+
+instance (Ord t, Fractional t, Serial m t) => Serial m (Ratio01 t) where
+  series = Ratio01 <$> (series `suchThatSerial` (\x -> 0 <= x && x <= 1))
+
+instance Show t => Show (Ratio01 t) where
+  showsPrec n (Ratio01 x) = showsPrec n x
+
+instance (Ord t, Fractional t, Arbitrary t) => Arbitrary (InitPair' t) where
+  arbitrary = f <$> liftM2 (,) arbitrary arbitrary where
+    f :: (Num t, Ord t, Fractional t) => (Ratio01 t, Ratio01 t) -> InitPair' t
+    f (Ratio01 x, Ratio01 y)
+      | 100*x<5   = Corput01
+      | 100*x<10  = Corput12
+      | otherwise = Mix x' y' where
+        x' = x*10/9
+        y' = y*(1-x)
+  shrink = genericShrink
+
+instance (Ord t, Fractional t, Serial m t) => Serial m (InitPair' t) where
+  series = cons0 Corput01 \/ cons0 Corput12 \/ mseries
+    where
+      mseries = do
+        (Ratio01 x) <- series
+        (Ratio01 y) <- series
+        return $ Mix x (y * (1-x))
+
+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 (Arbitrary a) => Arbitrary (M3.Matrix3 a) where
+  arbitrary = M3.fromList <$> vectorOf 9 arbitrary
+  shrink = genericShrink
+
+suchThatSerial :: Series m a -> (a -> Bool) -> Series m a
+suchThatSerial s p = s >>= \x -> if p x then pure x else empty
+
+cons5 :: (Serial m a, Serial m b, Serial m c, Serial m d, Serial m e) =>
+         (a->b->c->d->e->f) -> Series m f
+cons5 f = decDepth $
+  f <$> series
+    <~> series
+    <~> series
+    <~> series
+    <~> series
+
+instance (Serial m a, Serial m b, Serial m c, Serial m d, Serial m e) => Serial m (a,b,c,d,e) where
+  series = cons5 (,,,,)
+
+cons6 :: (Serial m a, Serial m b, Serial m c, Serial m d, Serial m e, Serial m f) =>
+         (a->b->c->d->e->f->g) -> Series m g
+cons6 f = decDepth $
+  f <$> series
+    <~> series
+    <~> series
+    <~> series
+    <~> series
+    <~> series
+
+instance (Serial m a, Serial m b, Serial m c, Serial m d, Serial m e, Serial m f) => Serial m (a,b,c,d,e,f) where
+  series = cons6 (,,,,,)
+
+liftM6  :: (Monad m) => (a1 -> a2 -> a3 -> a4 -> a5 -> a6 -> r) -> m a1 -> m a2 -> m a3 -> m a4 -> m a5 -> m a6 -> m r
+liftM6 f m1 m2 m3 m4 m5 m6 = do { x1 <- m1; x2 <- m2; x3 <- m3; x4 <- m4; x5 <- m5; x6 <- m6; return (f x1 x2 x3 x4 x5 x6) }
+
+instance (Arbitrary a, Arbitrary b, Arbitrary c, Arbitrary d, Arbitrary e, Arbitrary f)
+      => Arbitrary (a,b,c,d,e,f)
+ where
+  arbitrary = liftM6 (,,,,,) arbitrary arbitrary arbitrary arbitrary arbitrary arbitrary
+
+  shrink (u, v, w, x, y, z) =
+    [ (u', v', w', x', y', z')
+    | (u', (v', (w', (x', (y', z'))))) <- shrink (u, (v, (w, (x, (y, z))))) ]
+
+newtype Sorted t = Sorted t
+  deriving (Show, Generic)
+
+instance (Ord t, Arbitrary t) => Arbitrary (Sorted (t, t)) where
+  arbitrary = Sorted <$> (arbitrary `suchThat` uncurry (<=))
+
+instance (Ord t, Serial m t) => Serial m (Sorted (t, t)) where
+  series = Sorted <$> (series `suchThatSerial` uncurry (<=))
+
+instance (Ord t, Arbitrary t) => Arbitrary (Sorted (t, t, t)) where
+  arbitrary = Sorted <$> (arbitrary `suchThat` (\(a, b, c) -> a <= b && b <= c))
+
+instance (Ord t, Serial m t) => Serial m (Sorted (t, t, t)) where
+  series = Sorted <$> (series `suchThatSerial` (\(a, b, c) -> a <= b && b <= c))
+
+instance (Ord t, Arbitrary t) => Arbitrary (Sorted (t, t, t, t)) where
+  arbitrary = Sorted <$> (arbitrary `suchThat` (\(a, b, c, d) -> a <= b && b <= c && c <= d))
+
+instance (Ord t, Serial m t) => Serial m (Sorted (t, t, t, t)) where
+  series = Sorted <$> (series `suchThatSerial` (\(a, b, c, d) -> a <= b && b <= c && c <= d))
+
+instance (Ord t, Arbitrary t) => Arbitrary (Sorted (t, t, t, t, t)) where
+  arbitrary = Sorted <$> (arbitrary `suchThat` (\(a, b, c, d, e) -> a <= b && b <= c && c <= d && d <= e))
+
+instance (Ord t, Serial m t) => Serial m (Sorted (t, t, t, t, t)) where
+  series = Sorted <$> (series `suchThatSerial` (\(a, b, c, d, e) -> a <= b && b <= c && c <= d && d <= e))
+
+instance (Ord t, Arbitrary t) => Arbitrary (Sorted (t, t, t, t, t, t)) where
+  arbitrary = Sorted <$> (arbitrary `suchThat` (\(a, b, c, d, e, f) -> a <= b && b <= c && c <= d && d <= e && e <= f))
+
+instance (Ord t, Serial m t) => Serial m (Sorted (t, t, t, t, t, t)) where
+  series = Sorted <$> (series `suchThatSerial` (\(a, b, c, d, e, f) -> a <= b && b <= c && c <= d && d <= e && e <= f))
+
diff --git a/tests/Ivic.hs b/tests/Ivic.hs
new file mode 100644
--- /dev/null
+++ b/tests/Ivic.hs
@@ -0,0 +1,119 @@
+module Ivic where
+
+import Data.Ratio
+import Math.ExpPairs
+import Math.ExpPairs.Ivic
+
+import Test.Tasty
+import Test.Tasty.SmallCheck as SC
+import Test.Tasty.QuickCheck as QC
+import Test.Tasty.HUnit
+
+import Instances
+import Etalon (testEtalon)
+
+fromMinus3To3 :: Rational -> Rational
+fromMinus3To3 n = (n - 1 % 2) * 6
+
+fromHalfToOne :: Rational -> Rational
+fromHalfToOne n = n / 2 + 1 % 2
+
+testZetaOnS1 :: Sorted (Ratio01 Rational, Ratio01 Rational) -> Bool
+testZetaOnS1 (Sorted (Ratio01 a', Ratio01 b')) = a == b || za >= zb where
+  [ a,  b] = map fromMinus3To3 [a', b']
+  [za, zb] = map (optimalValue . zetaOnS) [a, b]
+
+-- May fail due to the granularity of 'sect'.
+testZetaOnS2 :: Sorted (Ratio01 Rational, Ratio01 Rational) -> Bool
+testZetaOnS2 (Sorted (Ratio01 a, Ratio01 b)) = a == b || za > zb where
+  [za, zb] = map (optimalValue . zetaOnS) [a, b]
+
+testZetaOnSsym :: Ratio01 Rational -> Bool
+testZetaOnSsym (Ratio01 a') = (toRational . abs) (za - za') == abs (a - 1 % 2) where
+  a   = fromMinus3To3 a'
+  za  = optimalValue $ zetaOnS a
+  za' = optimalValue $ zetaOnS (1 - a)
+
+testZetaOnSZero :: Ratio01 Rational -> Bool
+testZetaOnSZero (Ratio01 a') = a < 1 || optimalValue (zetaOnS a) == 0 where
+  a = fromMinus3To3 a'
+
+testMOnS1 :: Sorted (Ratio01 Rational, Ratio01 Rational) -> Bool
+testMOnS1 (Sorted (Ratio01 a', Ratio01 b')) = a == b || za <= zb where
+  [ a,  b] = map fromMinus3To3 [a', b']
+  [za, zb] = map (optimalValue . mOnS) [a, b]
+
+testMOnS2 :: Sorted (Ratio01 Rational, Ratio01 Rational) -> Bool
+testMOnS2 (Sorted (Ratio01 a', Ratio01 b')) = a == b || za < zb where
+  [ a,  b] = map fromHalfToOne [a', b']
+  [za, zb] = map (optimalValue . mOnS) [a, b]
+
+testMOnSZero :: Ratio01 Rational -> Bool
+testMOnSZero (Ratio01 a') = a >= 1%2 || (optimalValue . mOnS) a == 0 where
+  a = fromMinus3To3 a'
+
+testMOnSInf :: Ratio01 Rational -> Bool
+testMOnSInf (Ratio01 a') = a < 1 || (optimalValue . mOnS) a == InfPlus where
+  a = fromMinus3To3 a'
+
+testZetaReverse :: Ratio01 Rational -> Bool
+testZetaReverse (Ratio01 s') = abs (s - t) <= 5 % 1000 where
+  s = s' / 2
+  zs = zetaOnS s
+  t = toRational $ optimalValue $ reverseZetaOnS $ toRational $ optimalValue zs
+
+-- Convexity tests - they fail and it is OK
+testZetaConvex :: Sorted (Ratio01 Rational, Ratio01 Rational, Ratio01 Rational) -> Bool
+testZetaConvex (Sorted (Ratio01 a, Ratio01 b, Ratio01 c)) = a == b || b == c || zb <= k * Finite b + l where
+  [za, zb, zc] = map (optimalValue . zetaOnS) [a, b, c]
+  k = (za - zc) / Finite (a - c)
+  l = za - k * Finite a
+
+-- Ivic, Th. 8.1, p. 205
+testMConvex :: Sorted (Ratio01 Rational, Ratio01 Rational, Ratio01 Rational) -> Bool
+testMConvex (Sorted (Ratio01 a', Ratio01 b', Ratio01 c')) = a==b || b==c || za==InfPlus || zc==InfPlus
+  || zb>= za*zc*Finite(c-a)/(zc*Finite(c-b) + za*Finite(b-a)) where
+    [a,b,c] = map fromHalfToOne [a', b', c']
+    [za, zb, zc] = map (optimalValue . mOnS) [a,b,c] :: [RationalInf]
+
+etalonZetaOnS :: Integer -> Integer -> Integer -> Integer -> Bool
+etalonZetaOnS a b c d = Finite (c%d) >= optimalValue (zetaOnS $ a%b)
+
+etalonMOnS :: Integer -> Integer -> Integer -> Integer -> Bool
+etalonMOnS a b c d = Finite (c%d) <= (optimalValue . mOnS) (a%b)
+
+testSuite :: TestTree
+testSuite = testGroup "Ivic"
+  [ testCase "etalon zetaOnS"
+    (testEtalon 100 (\(a:b:c:d:_) -> etalonZetaOnS a b c d) "tests/etalon-zetaOnS.txt")
+  , testCase "etalon mOnS"
+    (testEtalon 100 (\(a:b:c:d:_) -> etalonMOnS a b c d) "tests/etalon-mOnS.txt")
+  , adjustOption (\(SC.SmallCheckDepth n) -> SC.SmallCheckDepth (n `div` 2)) $
+      SC.testProperty "zetaOnS monotonic" testZetaOnS1
+  , QC.testProperty "zetaOnS monotonic" testZetaOnS1
+  , adjustOption (\(SC.SmallCheckDepth n) -> SC.SmallCheckDepth (n `div` 2)) $
+      SC.testProperty "zetaOnS strict monotonic" testZetaOnS2
+  , QC.testProperty "zetaOnS strict monotonic" testZetaOnS2
+  , adjustOption (\(SC.SmallCheckDepth n) -> SC.SmallCheckDepth (n `div` 2)) $
+      SC.testProperty "mOnS monotonic" testMOnS1
+  , QC.testProperty "mOnS monotonic" testMOnS1
+  -- , adjustOption (\(SC.SmallCheckDepth n) -> SC.SmallCheckDepth (n `div` 2)) $
+  --     SC.testProperty "mOnS strict monotonic" testMOnS2
+  -- , QC.testProperty "mOnS strict monotonic" testMOnS2
+  , SC.testProperty "zetaOnS reverse" testZetaReverse
+  , QC.testProperty "zetaOnS reverse" testZetaReverse
+  , SC.testProperty "zetaOnS symmetry" testZetaOnSsym
+  , QC.testProperty "zetaOnS symmetry" testZetaOnSsym
+  , SC.testProperty "zetaOnS above s=1" testZetaOnSZero
+  , QC.testProperty "zetaOnS above s=1" testZetaOnSZero
+  , SC.testProperty "mOnS below s=1/2" testMOnSZero
+  , QC.testProperty "mOnS below s=1/2" testMOnSZero
+  , SC.testProperty "mOnS above s=1" testMOnSInf
+  , QC.testProperty "mOnS above s=1" testMOnSInf
+  -- , SC.testProperty "mOnS convex" testMConvex
+  -- , QC.testProperty "mOnS convex" testMConvex
+  -- , SC.testProperty "zetaOnS convex" testZetaConvex
+  -- , QC.testProperty "zetaOnS convex" testZetaConvex
+  ]
+
+
diff --git a/tests/Kratzel.hs b/tests/Kratzel.hs
new file mode 100644
--- /dev/null
+++ b/tests/Kratzel.hs
@@ -0,0 +1,72 @@
+module Kratzel where
+
+import Data.Ratio
+import Math.ExpPairs
+import Math.ExpPairs.Kratzel
+
+import Test.Tasty
+import Test.Tasty.SmallCheck as SC
+import Test.Tasty.QuickCheck as QC hiding (Positive)
+import Test.Tasty.HUnit
+
+import Instances
+import Etalon (testEtalon)
+
+testAbMonotonic :: Sorted (Positive Integer, Positive Integer, Positive Integer, Positive Integer) -> Bool
+testAbMonotonic (Sorted (Positive a, Positive c, Positive b, Positive d))
+  = (a == c && b == d) || zab > zcd
+    where
+      zab = optimalValue $ snd $ tauab a b
+      zcd = optimalValue $ snd $ tauab c d
+
+testAbCompareLow :: Sorted (Positive Integer, Positive Integer) -> Bool
+testAbCompareLow (Sorted (Positive a, Positive b))
+  = optimalValue (snd $ tauab a b) >= Finite (1 % (2 * a + 2 * b))
+
+testAbCompareHigh :: Sorted (Positive Integer, Positive Integer) -> Bool
+testAbCompareHigh (Sorted (Positive a, Positive b))
+  = optimalValue (snd $ tauab a b) < Finite (1 % (a + b))
+
+testAbcMonotonic :: Sorted (Positive Integer, Positive Integer, Positive Integer, Positive Integer, Positive Integer, Positive Integer) -> Bool
+testAbcMonotonic (Sorted (Positive a, Positive d, Positive b, Positive e, Positive c, Positive f))
+  = (a == d && b == e && c == f) || theoremAbc `elem` [Kolesnik, Kr64] || zabc >= zdef
+    where
+      (theoremAbc, resultAbc) = tauabc a b c
+      zabc = optimalValue resultAbc
+      zdef = optimalValue $ snd $ tauabc d e f
+
+testAbcCompareLow :: Sorted (Positive Integer, Positive Integer, Positive Integer) -> Bool
+testAbcCompareLow (Sorted (Positive a, Positive b, Positive c))
+  = c >= a + b || optimalValue (snd $ tauabc a b c) >= Finite (1 % (a + b + c))
+
+testAbcCompareHigh :: Sorted (Positive Integer, Positive Integer, Positive Integer) -> Bool
+testAbcCompareHigh (Sorted (Positive a, Positive b, Positive c))
+  = c >= a + b || optimalValue (snd $ tauabc a b c) < Finite (2 % (a + b + c))
+
+etalonTauab :: Integer -> Integer -> Integer -> Integer -> Bool
+etalonTauab a b c d = Finite (c % d) >= (optimalValue . snd) (tauab a b)
+
+etalonTauabc :: Integer -> Integer -> Integer -> Integer -> Integer -> Bool
+etalonTauabc a b c d e = Finite (d % e) >= (optimalValue . snd) (tauabc a b c)
+
+testSuite :: TestTree
+testSuite = testGroup "Kratzel"
+  [ testCase "etalon tauab"
+    (testEtalon 100 (\[x1, x2, x3, x4] -> etalonTauab x1 x2 x3 x4) "tests/etalon-tauab.txt")
+  , testCase "etalon tauabc"
+    (testEtalon 100 (\[x1, x2, x3, x4, x5] -> etalonTauabc x1 x2 x3 x4 x5) "tests/etalon-tauabc.txt")
+  , SC.testProperty "tauabc compare with 1/(a+b+c)" testAbcCompareLow
+  , QC.testProperty "tauabc compare with 1/(a+b+c)" testAbcCompareLow
+  , SC.testProperty "tauabc compare with 2/(a+b+c)" testAbcCompareHigh
+  , QC.testProperty "tauabc compare with 2/(a+b+c)" testAbcCompareHigh
+  , adjustOption (\(SC.SmallCheckDepth n) -> SC.SmallCheckDepth (n `div` 3)) $
+      SC.testProperty "tauabc monotonic" testAbcMonotonic
+  , QC.testProperty "tauabc monotonic" testAbcMonotonic
+  , SC.testProperty "tauab compare with 1/2(a+b)" testAbCompareLow
+  , QC.testProperty "tauab compare with 1/2(a+b)" testAbCompareLow
+  , SC.testProperty "tauab compare with 1/(a+b)" testAbCompareHigh
+  , QC.testProperty "tauab compare with 1/(a+b)" testAbCompareHigh
+  , adjustOption (\(SC.SmallCheckDepth n) -> SC.SmallCheckDepth (n `div` 2)) $
+      SC.testProperty "tauab monotonic" testAbMonotonic
+  , QC.testProperty "tauab monotonic" testAbMonotonic
+  ]
diff --git a/tests/LinearForm.hs b/tests/LinearForm.hs
new file mode 100644
--- /dev/null
+++ b/tests/LinearForm.hs
@@ -0,0 +1,86 @@
+{-# OPTIONS_GHC -fno-warn-type-defaults #-}
+
+module LinearForm where
+
+import Data.Ratio
+
+import Math.ExpPairs.LinearForm
+import Math.ExpPairs.RatioInf
+
+import Test.Tasty
+import Test.Tasty.SmallCheck as SC
+import Test.Tasty.QuickCheck as QC
+
+import Instances ()
+
+extractCoeffs :: Num t => LinearForm t -> (t, t, t)
+extractCoeffs lf =
+  ( evalLF (1, 0, 0) lf
+  , evalLF (0, 1, 0) lf
+  , evalLF (0, 0, 1) lf
+  )
+
+testPlus :: Rational -> Rational -> Rational -> Rational -> Rational -> Rational -> Bool
+testPlus a b c d e f = a+d==ad && b+e==be && c+f==cf where
+  l1 = LinearForm a b c
+  l2 = LinearForm d e f
+  (ad, be, cf) = extractCoeffs (l1 + l2)
+
+testMinus :: Rational -> Rational -> Rational -> Rational -> Rational -> Rational -> Bool
+testMinus a b c d e f = a-d==ad && b-e==be && c-f==cf where
+  l1 = LinearForm a b c
+  l2 = LinearForm d e f
+  (ad, be, cf) = extractCoeffs (l1 - l2)
+
+testFromInteger :: Integer -> Bool
+testFromInteger a = evalLF (0, 0, 1) (fromInteger a) == a
+
+testSubstitute1 :: LinearForm Rational -> Bool
+testSubstitute1 a
+  =  substituteLF (a, 0, 0) (LinearForm 1 0 0) == a
+  && substituteLF (0, a, 0) (LinearForm 0 1 0) == a
+  && substituteLF (0, 0, a) (LinearForm 0 0 1) == a
+
+testSubstitute2 :: LinearForm Rational -> LinearForm Rational
+                -> LinearForm Rational -> LinearForm Rational
+                -> LinearForm Rational -> LinearForm Rational
+                -> LinearForm Rational -> Bool
+testSubstitute2 a1 a2 b1 b2 c1 c2 lf
+  =  substituteLF (a1 + a2, b1 + b2, c1 + c2) lf
+  == substituteLF (a1, b1, c1) lf + substituteLF (a2, b2, c2) lf
+
+testNegateRF :: RationalForm Rational -> Integer -> Integer -> Integer -> Bool
+testNegateRF rf k l m = case evalRF (k, l, m) rf of
+  x@Finite{} -> x == negate (evalRF (k, l, m) (negate rf))
+  _          -> True
+
+testNegateVarsRF :: RationalForm Rational -> Integer -> Integer -> Integer -> Bool
+testNegateVarsRF rf k l m =
+  evalRF (k, l, m) rf == evalRF (-k, -l, -m) rf
+
+testFromIntegerRF :: Integer -> Bool
+testFromIntegerRF a = evalRF (0, 0, 1) (fromInteger a) == Finite (a % 1)
+
+testCheckConstraint :: Integer -> Integer -> Integer -> Constraint Rational -> Bool
+testCheckConstraint k l m c@(Constraint lf ineq)
+  =  (ineq==Strict    && isZero || x || y)
+  && (ineq==NonStrict && isZero || not (x && y))
+  where
+    x = checkConstraint (k, l, m) c
+    y = checkConstraint (k, l, m) (Constraint (negate lf) ineq)
+    isZero = evalLF (fromInteger k, fromInteger l, fromInteger m) lf == 0
+
+testSuite :: TestTree
+testSuite = testGroup "LinearForm"
+  [ QC.testProperty "plus" testPlus
+  , QC.testProperty "minus" testMinus
+  , SC.testProperty "from integer LF" testFromInteger
+  , QC.testProperty "from integer LF" testFromInteger
+  , QC.testProperty "substitute component" testSubstitute1
+  , QC.testProperty "substitution is linear" testSubstitute2
+  , QC.testProperty "negate RF" testNegateRF
+  , QC.testProperty "negate vars RF" testNegateVarsRF
+  , SC.testProperty "from integer RF" testFromIntegerRF
+  , QC.testProperty "from integer RF" testFromIntegerRF
+  , QC.testProperty "constraint" testCheckConstraint
+  ]
diff --git a/tests/Matrix3.hs b/tests/Matrix3.hs
new file mode 100644
--- /dev/null
+++ b/tests/Matrix3.hs
@@ -0,0 +1,65 @@
+module Matrix3 where
+
+import qualified Data.Matrix as M
+import qualified Math.ExpPairs.Matrix3 as M3
+
+import Test.Tasty
+import Test.Tasty.QuickCheck as QC
+
+import Instances ()
+
+toM :: M3.Matrix3 a -> M.Matrix a
+toM = M.fromList 3 3 . M3.toList
+
+toM3 :: M.Matrix a -> M3.Matrix3 a
+toM3 = M3.fromList . M.toList
+
+testOp :: (M3.Matrix3 Integer -> M3.Matrix3 Integer -> M3.Matrix3 Integer) -> (M.Matrix Integer -> M.Matrix Integer -> M.Matrix Integer) -> M3.Matrix3 Integer -> M3.Matrix3 Integer -> Bool
+testOp op1 op2 m1 m2 = m'==m'' where
+  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)
+
+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 "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 "mult column" testMultCol
+  ]
diff --git a/tests/MenzerNowak.hs b/tests/MenzerNowak.hs
new file mode 100644
--- /dev/null
+++ b/tests/MenzerNowak.hs
@@ -0,0 +1,40 @@
+module MenzerNowak where
+
+import Data.Ratio
+import Math.ExpPairs
+import Math.ExpPairs.MenzerNowak
+import Math.ExpPairs.Kratzel
+
+import Test.Tasty
+import Test.Tasty.SmallCheck as SC
+import Test.Tasty.QuickCheck as QC hiding (Positive)
+
+import Instances
+
+testMonotonic :: Sorted (Positive Integer, Positive Integer, Positive Integer, Positive Integer) -> Bool
+testMonotonic (Sorted (Positive a, Positive c, Positive b, Positive d))
+  =  (a == c && b == d) || zab > zcd
+    where
+      zab = optimalValue $ menzerNowak a b
+      zcd = optimalValue $ menzerNowak c d
+
+testCompareLow :: Sorted (Positive Integer, Positive Integer) -> Bool
+testCompareLow (Sorted (Positive a, Positive b))
+  = optimalValue (snd $ tauab a b) <= optimalValue (menzerNowak a b) + Finite eps
+    where
+      eps = 1 % (10 ^ (30::Integer))
+
+testCompareHigh :: Sorted (Positive Integer, Positive Integer) -> Bool
+testCompareHigh (Sorted (Positive a, Positive b))
+  = optimalValue (menzerNowak a b) < 1
+
+testSuite :: TestTree
+testSuite = testGroup "MenzerNowak"
+  [ SC.testProperty "compare with tauab" testCompareLow
+  , QC.testProperty "compare with tauab" testCompareLow
+  , SC.testProperty "compare with 1" testCompareHigh
+  , QC.testProperty "compare with 1" testCompareHigh
+  , adjustOption (\(SC.SmallCheckDepth n) -> SC.SmallCheckDepth (n `div` 2)) $
+      SC.testProperty "monotonic" testMonotonic
+  , QC.testProperty "monotonic" testMonotonic
+  ]
diff --git a/tests/Pair.hs b/tests/Pair.hs
new file mode 100644
--- /dev/null
+++ b/tests/Pair.hs
@@ -0,0 +1,33 @@
+module Pair where
+
+import Math.ExpPairs.Pair (InitPair, initPairToValue, initPairToProjValue)
+
+import Data.Ratio
+import Test.Tasty
+import Test.Tasty.SmallCheck as SC
+import Test.Tasty.QuickCheck as QC
+
+import Instances ()
+
+testBounds :: InitPair -> Bool
+testBounds ip = k>=0 && k<=1%2 && l>=1%2 && l<=1 where
+  (k, l) = initPairToValue ip
+
+fracs2proj :: (Rational, Rational) -> (Integer, Integer, Integer)
+fracs2proj (q, r) = (k, l, m) where
+  dq = denominator q
+  dr = denominator r
+  m = lcm dq dr
+  k = numerator q * (m `div` dq)
+  l = numerator r * (m `div` dr)
+
+testProjective :: InitPair -> Bool
+testProjective ip = initPairToProjValue ip == fracs2proj (initPairToValue ip)
+
+testSuite :: TestTree
+testSuite = testGroup "Pair"
+  [ SC.testProperty "bounds" testBounds
+  , QC.testProperty "bounds" testBounds
+  , SC.testProperty "projective" testProjective
+  , QC.testProperty "projective" testProjective
+  ]
diff --git a/tests/PrettyProcess.hs b/tests/PrettyProcess.hs
new file mode 100644
--- /dev/null
+++ b/tests/PrettyProcess.hs
@@ -0,0 +1,20 @@
+module PrettyProcess where
+
+import Math.ExpPairs.ProcessMatrix
+import Math.ExpPairs.PrettyProcess
+
+import Test.Tasty
+import Test.Tasty.SmallCheck as SC
+import Test.Tasty.QuickCheck as QC hiding (Positive)
+
+import Instances ()
+
+testUglifyPrettify :: [Process] -> Bool
+testUglifyPrettify xs = uglify (prettify xs) == xs
+
+testSuite :: TestTree
+testSuite = testGroup "PrettyProcess"
+  [ adjustOption (\(SC.SmallCheckDepth n) -> SC.SmallCheckDepth (n `min` 13)) $
+      SC.testProperty "uglify . prettify == id" testUglifyPrettify
+  , QC.testProperty "uglify . prettify == id" testUglifyPrettify
+  ]
diff --git a/tests/RatioInf.hs b/tests/RatioInf.hs
new file mode 100644
--- /dev/null
+++ b/tests/RatioInf.hs
@@ -0,0 +1,51 @@
+module RatioInf where
+
+import Math.ExpPairs.RatioInf (RatioInf (..), RationalInf)
+
+import Test.Tasty
+import Test.Tasty.SmallCheck as SC
+
+testPlus :: Rational -> Rational -> Bool
+testPlus a b = Finite (a+b) == Finite a + Finite b
+
+testMinus :: Rational -> Rational -> Bool
+testMinus a b = Finite (a-b) == Finite a - Finite b
+
+testMultiply :: Rational -> Rational -> Bool
+testMultiply a b = Finite (a*b) == Finite a * Finite b
+
+testDivide :: Rational -> Rational -> Bool
+testDivide a b = b==0 || (Finite (a/b) == Finite a / Finite b)
+
+testInfPlus :: RationalInf -> Rational -> Bool
+testInfPlus a b =  a + Finite b == a
+
+testInfMinus :: RationalInf -> Rational -> Bool
+testInfMinus a b = a - Finite b == a
+
+testInfMultiply :: RationalInf -> Rational -> Bool
+testInfMultiply a b = b==0 || a * Finite b * Finite b == a
+
+testInfDivide :: RationalInf -> Rational -> Bool
+testInfDivide a b =  b==0 || a / Finite b / Finite b == a
+
+testConversion :: Rational -> Bool
+testConversion a = toRational (Finite a) == a
+
+testSuite :: TestTree
+testSuite = testGroup "RatioInf"
+  [ SC.testProperty "plus"                testPlus
+  , SC.testProperty "minus"               testMinus
+  , SC.testProperty "multiply"            testMultiply
+  , SC.testProperty "divide"              testDivide
+  , SC.testProperty "infplus plus"      $ testInfPlus InfPlus
+  , SC.testProperty "infplus minus"     $ testInfPlus InfMinus
+  , SC.testProperty "infminus plus"     $ testInfMinus InfPlus
+  , SC.testProperty "infminus minus"    $ testInfMinus InfMinus
+  , SC.testProperty "infmultiply plus"  $ testInfMultiply InfPlus
+  , SC.testProperty "infmultiply minus" $ testInfMultiply InfMinus
+  , SC.testProperty "infdivide plus"    $ testInfDivide InfPlus
+  , SC.testProperty "infdivide minus"   $ testInfDivide InfMinus
+  , SC.testProperty "conversion"          testConversion
+  ]
+
diff --git a/tests/Tests.hs b/tests/Tests.hs
--- a/tests/Tests.hs
+++ b/tests/Tests.hs
@@ -15,15 +15,12 @@
 
 tests :: TestTree
 tests = testGroup "Tests"
-	[ Matrix3.testSuite
-	, LinearForm.testSuite
-	, RatioInf.testSuite
-	, Pair.testSuite
-	, PrettyProcess.testSuite
-	, Ivic.testSuite
-	, Kratzel.testSuite
-	, MenzerNowak.testSuite
-	]
-
-
-
+  [ Matrix3.testSuite
+  , LinearForm.testSuite
+  , RatioInf.testSuite
+  , Pair.testSuite
+  , PrettyProcess.testSuite
+  , Ivic.testSuite
+  , Kratzel.testSuite
+  , MenzerNowak.testSuite
+  ]