exp-pairs 0.2.0.0 → 0.2.1.0

raw patch · 20 files changed

+726/−279 lines, 20 filesdep +bimapdep +gaugedep +raw-strings-qqdep ~basedep ~smallcheckPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

Dependencies added: bimap, gauge, raw-strings-qq

Dependency ranges changed: base, smallcheck

API changes (from Hackage documentation)

Files

• CHANGELOG.md +6/−0
• Math/ExpPairs.hs +10/−8
• Math/ExpPairs/Ivic.hs +45/−29
• Math/ExpPairs/Kratzel.hs +39/−16
• Math/ExpPairs/LinearForm.hs +10/−10
• Math/ExpPairs/Matrix3.hs +85/−75
• Math/ExpPairs/MenzerNowak.hs +3/−4
• Math/ExpPairs/Pair.hs +4/−6
• Math/ExpPairs/PrettyProcess.hs +9/−12
• Math/ExpPairs/Process.hs +10/−10
• Math/ExpPairs/ProcessMatrix.hs +18/−16
• Math/ExpPairs/RatioInf.hs +7/−13
• auxiliary/BenchMatrix.hs +17/−0
• auxiliary/OptimizeSum.hs +289/−0
• exp-pairs.cabal +30/−4
• tests/Etalon.hs +1/−1
• tests/Instances.hs +86/−52
• tests/Ivic.hs +28/−20
• tests/Kratzel.hs +26/−1
• tests/Matrix3.hs +3/−2

@@ -1,6 +1,12 @@ Changes ======= +Version 0.2.1.0+----------------++Improvements for 4D asymmetric divisor problems.+Refine matrix manipulations.+ Version 0.2.0.0 ---------------- 

@@ -1,12 +1,11 @@ {-| Module      : Math.ExpPairs-Description : Linear programming over exponent pairs-Copyright   : (c) Andrew Lelechenko, 2014-2015+Copyright   : (c) Andrew Lelechenko, 2014-2020 License     : GPL-3 Maintainer  : andrew.lelechenko@gmail.com-Stability   : experimental-Portability : POSIX +Linear programming over exponent pairs+ Package implements an algorithm to minimize the maximum of a list of rational objective functions over the set of exponent pairs. See full description in A. V. Lelechenko, Linear programming over exponent pairs. Acta Univ. Sapientiae, Inform. 5, No. 2, 271-287 (2013). <http://www.acta.sapientia.ro/acta-info/C5-2/info52-7.pdf>@@ -15,7 +14,6 @@ "Math.ExpPairs.Ivic", "Math.ExpPairs.Kratzel" and "Math.ExpPairs.MenzerNowak". -} -{-# LANGUAGE CPP             #-} {-# LANGUAGE PatternSynonyms #-}  module Math.ExpPairs@@ -45,7 +43,6 @@ import Data.Function (on) import Data.Ord      (comparing) import Data.List     (minimumBy)-import Data.Monoid import Data.Ratio import Data.Text.Prettyprint.Doc hiding ((<>)) import qualified Data.Text.Prettyprint.Doc as PP@@ -57,10 +54,15 @@ import Math.ExpPairs.RatioInf  -- | For a given @c@ returns linear form @c * k@+pattern K :: (Eq a, Num a) => a -> LinearForm a pattern K n = LinearForm n 0 0+ -- | For a given @c@ returns linear form @c * l@+pattern L :: (Eq a, Num a) => a -> LinearForm a pattern L n = LinearForm 0 n 0+ -- | For a given @c@ returns linear form @c * m@+pattern M :: (Eq a, Num a) => a -> LinearForm a pattern M n = LinearForm 0 0 n  -- | Build a constraint, which states that the value of the first linear form is greater than the value of the second one.@@ -85,7 +87,7 @@  evalFunctional :: [InitPair] -> [InitPair] -> [RationalForm Rational] -> [Constraint Rational] -> Path -> (RationalInf, InitPair) evalFunctional corners interiors rfs cons path = case rs of-  [] -> (InfPlus, error "evalFunctional: cannot find any exponential pair, which satisfies constraints")+  [] -> (InfPlus, error $ "evalFunctional: cannot find any exponential pair, which satisfies constraints " ++ show (map pretty cons))   _  -> minimumBy (comparing fst) rs   where     applyPath  = map (evalPath path . initPairToProjValue &&& id)@@ -168,5 +170,5 @@       pathba  = path <> baPath       branchB@(OptimizeResult r2' _ _) = optimize' rfs cons (OptimizeResult r1 ip1 pathba) -    consBuilder rr (num :/: den) = (substituteLF (num, den, 1) (L (toRational rr) - K 1)) >. 0+    consBuilder rr (num :/: den) = substituteLF (num, den, 1) (L (toRational rr) - K 1) >. 0 

@@ -1,15 +1,16 @@ {-| Module      : Math.ExpPairs.Ivic-Description : Riemann zeta-function-Copyright   : (c) Andrew Lelechenko, 2014-2015+Copyright   : (c) Andrew Lelechenko, 2014-2020 License     : GPL-3 Maintainer  : andrew.lelechenko@gmail.com-Stability   : experimental-Portability : POSIX -Provides functions to compute estimates Riemann zeta-function-ζ in a critical strip, given in  /Ivić A./ `The Riemann zeta-function: Theory and applications',-Mineola, New York: Dover Publications, 2003.+Estimates of the Riemann zeta-function+in a critical strip, according to+Ivić A. "The Riemann zeta-function: theory and applications",+Mineola, New York: Dover Publications, 2003,+and Lelechenko A. V. "Dirichlet divisor problem on Gaussian integers" in+Proceedings of the 6th international conference on analytic number theory+and spatial tesselations, Kyiv, 2018, vol. 1, p. 76-86.  -} module Math.ExpPairs.Ivic@@ -30,8 +31,9 @@  import Math.ExpPairs --- | Compute µ(σ) such that |ζ(σ+it)| ≪ |t|^µ(σ) .--- See equation (7.57) in Ivić2003.+-- | Compute \( \mu(\sigma) \) such that+-- \( \zeta(\sigma+it) \ll |t|^{\mu(\sigma)} \).+-- See equation (7.57) in Ivić, 2003. zetaOnS :: Rational -> OptimizeResult zetaOnS s   | s >= 1  = simulateOptimize 0@@ -63,8 +65,9 @@   | 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. Further justification will be published elsewhere.+-- | Compute maximal \( m(\sigma) \) such that+--  \( \int_1^T | \zeta(\sigma+it) |^{m(\sigma)} dt \ll T^{1+\varepsilon} \).+-- See equation (8.97) in Ivić, 2003. mOnS :: Rational -> OptimizeResult mOnS s   | s < 1%2 = simulateOptimize 0@@ -85,7 +88,7 @@     numer = t - scaleLF (4 * (1-s)) (K 1 + L 1) + scaleLF (4 * muS) (K 2 + L 2 + 1)     denom = scaleLF muS t -    cons = if s >= 2%3 then [] else [scaleLF s (K 4 + L 8 + 2) >=. K 2 + L 6 + 1]+    cons = [ scaleLF s (K 4 + L 8 + 2) >=. K 2 + L 6 + 1 | s < 2 % 3 ]      x2' = optimize [- numer :/: denom] cons     x2 = x2' {optimalValue = negate $ optimalValue x2'}@@ -108,13 +111,15 @@ mOnSTwoThird :: RationalInf mOnSTwoThird = optimalValue $ mOnS $ 2 % 3 --- | 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.--- Sometimes, when 'mOnS' gets especially lucky exponent pair, 'reverseMOnS' can miss--- real σ and returns significantly bigger value.+-- | Try to reverse 'mOnS': for a given precision and \( m \) compute \( \sigma \).+-- Implemented as a binary search, so its performance is very poor.+-- Since 'mOnS' is not monotonic, the result is not guaranteed to be neither+-- minimal nor maximal possible, but usually is close enough. ----- For integer m>=4 this function corresponds to the multidimensional Dirichlet problem--- and returns σ from error term O(x^{σ+ε}). See Ch. 13 in Ivić2003.+-- For integer \( m \ge 4 \) this function corresponds+-- to the multidimensional Dirichlet problem+-- and returns \( \sigma \) from error term \( O(x^{\sigma+\varepsilon}) \).+-- See Ch. 13 in Ivić, 2003. reverseMOnS :: Rational -> RationalInf -> Rational reverseMOnS _ InfPlus = 1 reverseMOnS _ (Finite m)@@ -126,23 +131,32 @@   where     go = binarySearch (\c -> optimalValue (mOnS c) > m) Greatest prec --- | An estimate of the symmetric multidimensional divisor function from Kolpakova, 2011.+-- | An estimate of the symmetric multidimensional divisor function from+-- Kolpakova O. V.,+-- "New estimates of the remainder in an asymptotic formula+-- in the multidimensional Dirichlet divisor problem", Mathematical Notes,+-- vol. 89, p. 504-518, 2011. kolpakova2011 :: Integer -> Double kolpakova2011 k = 1 - 1/3 * 2**(2/3) * (4.45 * fromInteger k)**(-2/3) --- | 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 σ.+-- | Check whether+-- \( \int_1^T \prod_i |\zeta(n_i\sigma+it|^{m_i} dt \ll T^{1+\varepsilon} \)+-- for a given list of pairs \( [(n_1, m_1), ...] \) and fixed \( \sigma \). 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 recip qs --- | Find for a given precision and list of pairs [(n_1, m_1), ...] the minimal σ--- such that ∫_1^T   Π_i|ζ(n_i*σ+it)|^m_i dt ≪ T^(1+ε).+-- | Find for a given precision and list of pairs \( [(n_1, m_1), ...] \)+-- the minimal \( \sigma \)+-- such that+-- \( \int_1^T \prod_i |\zeta(n_i\sigma+it|^{m_i} dt \ll T^{1+\varepsilon} \). findMinAbscissa :: Rational -> [(Rational, Rational)] -> Rational findMinAbscissa prec xs = binarySearch (checkAbscissa xs) Greatest prec (1 % 2 / minimum (map fst xs)) 1 --- | 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.+-- | For a given \( A \) compute minimal \( M(A) \) such that+-- \( \int_1^T |\zeta(1/2+it)|^A \ll T^{M(A)+\varepsilon} \)+-- See Ch. 8 in Ivić, 2003 and Th. 1 in Lelechenko, 2018. mBigOnHalf :: Rational -> OptimizeResult mBigOnHalf a   | a < 4     = simulateOptimize 1@@ -158,9 +172,11 @@ -- is produced by -- optimize [K 4 + L 4 + 2 :/: K 1] [26 >. K 26 + L 32] --- | Try to reverse 'mBigOnHalf': for a given M(A) find maximal possible A.--- Sometimes, when 'mBigOnHalf' gets especially lucky exponent pair, 'reverseMBigOnHalf' can miss--- real A and returns lower value.+-- | Try to reverse 'mBigOnHalf':+-- for a given \( M(A) \) find maximal possible \( A \).+-- Sometimes, when 'mBigOnHalf' gets especially lucky exponent pair,+-- 'reverseMBigOnHalf' can miss+-- real \( A \) and returns lower value. reverseMBigOnHalf :: Rational -> OptimizeResult reverseMBigOnHalf m   | m <= 2 = simulateOptimize $ (m-1)*8 + 4

@@ -1,12 +1,11 @@ {-| Module      : Math.ExpPairs.Kratzel-Description : Asymmetric divisor problem-Copyright   : (c) Andrew Lelechenko, 2014-2015+Copyright   : (c) Andrew Lelechenko, 2014-2020 License     : GPL-3 Maintainer  : andrew.lelechenko@gmail.com-Stability   : experimental-Portability : POSIX +Asymmetric divisor problem+ Let τ_{a, b}(n) denote the number of integer (v, w) with v^a w^b = n. @@ -70,9 +69,12 @@ divideResult :: Real a => a -> (b, OptimizeResult) -> (b, OptimizeResult) divideResult d = second (\o -> o {optimalValue = optimalValue o / Finite (toRational d)}) -tauab' :: Rational -> Rational -> (TauabTheorem, OptimizeResult)-tauab' a b = minimumBy (comparing snd) [kr511a, kr511b, kr512a, kr512b]+tauab' :: Integer -> Integer -> (TauabTheorem, OptimizeResult)+tauab' a' b' = minimumBy (comparing snd) [kr511a, kr511b, kr512a, kr512b]   where+    a = toRational a'+    b = toRational b'+     kr511a = (Kr511a, optimize       [K 2 + L 2 - 1 :/: M (a+b)]       [L (2 * a) >=. K (2 * b) + M a])@@ -97,7 +99,7 @@       d = a `gcd` b tauab a b = tauab' a' b'   where-    [a', b'] = sort $ map toRational [a, b]+    [a', b'] = sort [a, b]   -- |Special type to specify the theorem of Krätzel1988,@@ -129,16 +131,20 @@   pretty (Tauab t) = pretty t   pretty t         = pretty (show t) -tauabc' :: Rational -> Rational -> Rational -> (TauabcTheorem, OptimizeResult)-tauabc' a b c = minimumBy (comparing snd) [kr61, kr62, kr63, kr64, kr65, kr66]+tauabc' :: Integer -> Integer -> Integer -> (TauabcTheorem, OptimizeResult)+tauabc' a' b' c' = minimumBy (comparing snd) [kr61, kr62, kr63, kr64, kr65, kr66]   where+    a = toRational a'+    b = toRational b'+    c = toRational c'+     abc = a + b + c     kr61       | c<a+b = (Kr61, simulateOptimize $ 2/abc)       | optimalValue optRes < Finite (recip c) = (Kr61, simulateOptimize $ 1/c)       | otherwise = (Tauab th, optRes)       where-        (th, optRes) = tauab' a b+        (th, optRes) = tauab a' b'     kr62 = (Kr62, optimize       [K 2 + L 2 :/: M (a + b + c)]       [ L a >=. K (b + c)@@ -164,17 +170,24 @@       d = a `gcd` b `gcd` c tauabc a b c = tauabc' a' b' c'   where-    [a', b', c'] = sort $ map toRational [a, b, c]+    [a', b', c'] = sort [a, b, c]   -- |Special type to specify the theorem of Krätzel1988, -- which provided the best estimate of Θ(a, b, c, d) data TauabcdTheorem+  -- | Heath-Brown, 1978   = HeathBrown   | Tauabc TauabcTheorem+  -- | Theorem 6.11   | Kr611+  -- | Krätzel, Estimates in the general divisor problem,+  -- Abh. Math. Sem. Univ. Hamburg 62 (1992), 191-206,+  -- Theorem 2 for p = 4   | Kr1992_2+  -- | Ibidem, Theorem 3 for p = 4 under condition 3.1   | Kr1992_31+  -- | Ibidem, Theorem 3 for p = 4 under condition 3.2   | Kr1992_32   | Kr2010_1a   | Kr2010_1b@@ -187,18 +200,28 @@   pretty (Tauabc t) = pretty t   pretty t          = pretty (show t) -tauabcd' :: Rational -> Rational -> Rational -> Rational -> (TauabcdTheorem, OptimizeResult)-tauabcd' a1 a2 a3 a4 = minimumBy (comparing snd) [kr611, kr1992_2, kr1992_31, kr1992_32, kr2010_1a, kr2010_1b, kr2010_2, kr2010_3]+tauabcd' :: Integer -> Integer -> Integer -> Integer -> (TauabcdTheorem, OptimizeResult)+tauabcd' a1' a2' a3' a4' = minimumBy (comparing snd) [fallback, kr611, kr1992_2, kr1992_31, kr1992_32, kr2010_1a, kr2010_1b, kr2010_2, kr2010_3]   where+    a1 = toRational a1'+    a2 = toRational a2'+    a3 = toRational a3'+    a4 = toRational a4'+     a12 = a1 + a2     a123 = a1 + a2 + a3     a1234 = a1 + a2 + a3 + a4 +    (th3, optRes3) = tauabc a1' a2' a3'++    fallback = (Tauabc th3, optRes3 { optimalValue = optVal })+      where+        optVal = optimalValue optRes3 `max` Finite (1 % a4')+     kr611       | optimalValue optRes3 < Finite (recip a4) = (Kr611, optimize [form] cons)       | otherwise = (Tauabc th3, optRes3)       where-        (th3, optRes3) = tauabc' a1 a2 a3         form = K 2 + L 2 + 1 :/: M (a1 + a2 + a3 + a4)                              -- (6.46)         cons =           [ scaleLF a1 (L 2 - 1) >. K (2 * a4)                                      -- (6.41)@@ -229,7 +252,7 @@      kr1992_32 = (Kr1992_32, simulateOptimize $ if cond then val else 1)       where-        k = 32 % 205+        k = 13 % 84 -- Bourgain 2017         l = k + 1 % 2         val = (k + l + 2) / ((k + l) * a1 + a1234)         cond = (k + l - 2) * a4 < (k + l) * a1 + a1234@@ -267,7 +290,7 @@       d = a1 `gcd` a2 `gcd` a3 `gcd` a4 tauabcd a1 a2 a3 a4 = tauabcd' a1' a2' a3' a4'   where-    [a1', a2', a3', a4'] = sort $ map toRational [a1, a2, a3, a4]+    [a1', a2', a3', a4'] = sort [a1, a2, a3, a4]   -- |Special type to specify the theorem of Krätzel1988,

@@ -1,15 +1,17 @@-{-# LANGUAGE DeriveFunctor, DeriveFoldable, DeriveGeneric, CPP #-} {-| Module      : Math.ExpPairs.LinearForm-Description : Linear forms, rational forms and constraints-Copyright   : (c) Andrew Lelechenko, 2014-2015+Copyright   : (c) Andrew Lelechenko, 2014-2020 License     : GPL-3 Maintainer  : andrew.lelechenko@gmail.com-Stability   : experimental-Portability : POSIX +Linear forms, rational forms and constraints+ 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. -}++{-# LANGUAGE DeriveTraversable #-}+{-# LANGUAGE DeriveGeneric     #-}+ module Math.ExpPairs.LinearForm   ( LinearForm (..)   , scaleLF@@ -25,9 +27,7 @@ import Control.DeepSeq import Data.Foldable  (Foldable (..), toList) import Data.Maybe     (mapMaybe)-import Data.Monoid    (Monoid, mempty, mappend) import Data.Ratio     (numerator, denominator)-import Data.Semigroup (Semigroup, (<>)) import GHC.Generics   (Generic (..)) import Data.Text.Prettyprint.Doc @@ -37,7 +37,7 @@ -- 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, Show, Functor, Foldable, Generic)+  deriving (Eq, Show, Functor, Foldable, Traversable, Generic)  instance NFData t => NFData (LinearForm t) where   rnf = rnf . toList@@ -82,7 +82,7 @@  -- | Define a rational form of two variables, equal to the ratio of two 'LinearForm'. data RationalForm t = (LinearForm t) :/: (LinearForm t)-  deriving (Eq, Show, Functor, Foldable, Generic)+  deriving (Eq, Show, Functor, Foldable, Traversable, Generic) infix 5 :/:  instance (Num t, Eq t, Pretty t) => Pretty (RationalForm t) where@@ -129,7 +129,7 @@  -- |A linear constraint of two variables. data Constraint t = Constraint !(LinearForm t) !IneqType-  deriving (Eq, Show, Functor, Foldable, Generic)+  deriving (Eq, Show, Functor, Foldable, Traversable, Generic)  instance (Num t, Eq t, Pretty t) => Pretty (Constraint t) where   pretty (Constraint lf ineq) = pretty lf <+> pretty ineq <+> pretty "0"

@@ -1,25 +1,22 @@ {-| Module      : Math.ExpPairs.Matrix3-Description : Implements matrices of order 3-Copyright   : (c) Andrew Lelechenko, 2014-2015+Copyright   : (c) Andrew Lelechenko, 2014-2020 License     : GPL-3 Maintainer  : andrew.lelechenko@gmail.com-Stability   : experimental-Portability : POSIX -Provides types and functions for matrices and vectors of order 3.-Can be used instead of "Data.Matrix" to reduce overhead and simplify code.+Matrices of order 3+and efficient multiplication algorithms.+ -} -{-# LANGUAGE DeriveFunctor   #-}-{-# LANGUAGE DeriveFoldable  #-}-{-# LANGUAGE DeriveGeneric   #-}-{-# LANGUAGE RecordWildCards #-}+{-# LANGUAGE DeriveGeneric     #-}+{-# LANGUAGE DeriveTraversable #-}+{-# LANGUAGE RecordWildCards   #-}+{-# LANGUAGE Safe              #-}  module Math.ExpPairs.Matrix3   ( Matrix3 (..)   , fromList-  , toList   , det   , multCol   , normalize@@ -35,13 +32,7 @@ import GHC.Generics   (Generic (..)) import Data.Text.Prettyprint.Doc --- |Matrix of order 3. Instances of 'Num' and 'Fractional'--- are given in terms of the multiplicative group of matrices,--- not the additive one. E. g.,------ > toList 1 == [1,0,0,0,1,0,0,0,1]--- > toList 1 /= [1,1,1,1,1,1,1,1,1]---+-- | Matrix of order 3. data Matrix3 t = Matrix3 {   a11 :: !t,   a12 :: !t,@@ -53,7 +44,7 @@   a32 :: !t,   a33 :: !t   }-  deriving (Eq, Show, Functor, Foldable, Generic)+  deriving (Eq, Show, Functor, Foldable, Traversable, Generic)  instance NFData t => NFData (Matrix3 t) where   rnf = rnf . toList@@ -71,7 +62,8 @@   a33 = n   } -instance (Num t, Ord t) => Num (Matrix3 t) where+-- | 'fromInteger' returns a diagonal matrix+instance Num t => Num (Matrix3 t) where   a + b = Matrix3 {     a11 = a11 a + a11 b,     a12 = a12 a + a12 b,@@ -84,14 +76,10 @@     a33 = a33 a + a33 b     } -  (*) = usualMult--  negate = fmap negate--  abs = error "abs of Matrix3 is undefined"--  signum = diag . signum . det-+  (*)         = usualMult+  negate      = fmap negate+  abs         = id+  signum      = id   fromInteger = diag . fromInteger  usualMult :: Num t => Matrix3 t -> Matrix3 t -> Matrix3 t@@ -109,49 +97,59 @@ {-# SPECIALIZE usualMult :: Matrix3 Int -> Matrix3 Int -> Matrix3 Int #-} {-# SPECIALIZE usualMult :: Matrix3 Integer -> Matrix3 Integer -> Matrix3 Integer #-} --- | Multiplicate matrices by 23 multiplications and 68 additions.--- It becomes faster than usual multiplication (which requires 27 multiplications and 18 additions),--- when matrix's elements are large (several hundred digits) integers.+-- | Multiplicate matrices. Requires 23 multiplications and 62 additions.+-- It becomes faster than vanilla multiplication '(*)',+-- which requires 27 multiplications and 18 additions,+-- when matrix's elements are large (> 700 digits) integers. -- -- An algorithm follows--- /J. Laderman./ A noncommutative algorithm for multiplying 3 × 3 matrices using 23 multiplications. Bull. Amer. Math. Soc., 82:126–128, 1976.------ We were able to reduce the number of additions from 98 to 68 by sofisticated choice of intermediate variables.+-- /J. Laderman./+-- A noncommutative algorithm for multiplying \( 3 \times 3 \)+-- matrices using 23 multiplications.+-- Bull. Amer. Math. Soc., 82:126–128, 1976.+-- Our contribution is reducing the number of additions+-- from 98 to 62 by well-thought 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+    t33 = t37 + a22+    t34 = t38 + a32+    t35 = t39 - a32+    t36 = a22 - t40+    t37 = a23 - a13+    t38 = a31 - a11+    t39 = a13 - a33+    t40 = a11 - a21 -    u33 = b21 - b11 - b23 - b31-    u34 = b22 - b12-    u35 = b22 - b32-    u36 = b33 - b31-    u37 = b13 - b23+    u33 = b23 - u37+    u34 = u38 - b22+    u35 = b11 - u39+    u36 = b22 - u40+    u37 = b33 - b31+    u38 = b32 - b31+    u39 = b13 - b23+    u40 = b12 - b11 -    m1 = (t35 + t33 - a21) * b22-    m2 = (a11 - a21) * u34-    m3 = a22 * (u33 + b33 - u34)-    m4 = (a21 - t37) * (b11 + u34)-    m5 = (a22 + a21) * (b12 - b11)+    m1 = (t35 + a12 - t36) * b22+    m2 = t40 * (b22 - b12)+    m3 = a22 * (b21 - u36 - u33)+    m4 = t36 * u36+    m5 = (a22 + a21) * u40     m6 = a11 * b11-    m7 = (t36 + a32) * (b11 - u37)-    m8 = t36 * u37+    m7 = t34 * u35+    m8 = t38 * u39     m9 = (a31 + a32) * (b13 - b11)-    m10 = (t33 - a31 + t34) * b23-    m11 = a32 * (u33 + b13 - u35)-    m12 = (a32 - t35) * (b31 + u35)-    m13 = t35 * u35+    m10 = (a12 - t33 - t34) * b23+    m11 = a32 * (u34 + b21 - u35)+    m12 = t35 * u34+    m13 = t39 * (b22 - b32)     m14 = a13 * b31-    m15 = (a33 + a32) * (b32 - b31)-    m16 = (a22 - t34) * (b23 - u36)-    m17 = t34 * (b23 - b33)-    m18 = (a23 + a22) * u36+    m15 = (a33 + a32) * u38+    m16 = t33 * u33+    m17 = t37 * (b33 - b23)+    m18 = (a23 + a22) * u37     m19 = a12 * b21     m20 = a23 * b32     m21 = a21 * b13@@ -180,13 +178,16 @@  -- | Multiplicate matrices under assumption that multiplication of elements is commutative. -- Requires 22 multiplications and 66 additions.--- It becomes faster than usual multiplication (which requires 27 multiplications and 18 additions),--- when matrix's elements are large (several hundred digits) integers.+-- It becomes faster than vanilla multiplication '(*)',+-- which requires 27 multiplications and 18 additions,+-- when matrix's elements are large (> 700 digits) integers. -- -- An algorithm follows--- /O. M. Makarov./ An algorithm for multiplication of 3 × 3 matrices. Zh. Vychisl. Mat. i Mat. Fiz., 26(2):293–294, 320, 1986.------ We were able to reduce the number of additions from 105 to 66 by sofisticated choice of intermediate variables.+-- /O. M. Makarov./+-- An algorithm for multiplication of \( 3 \times 3 \) matrices.+-- Zh. Vychisl. Mat. i Mat. Fiz., 26(2):293–294, 320, 1986.+-- Our contribution is reducing the number of additions+-- from 105 to 66 by well-thought choice of intermediate variables. makarovMult :: Num t => Matrix3 t -> Matrix3 t -> Matrix3 t makarovMult   (Matrix3 k1 b1 c1 k2 b2 c2 k3 b3 c3)@@ -245,17 +246,18 @@     c33 = v32 + m4 + m21 {-# SPECIALIZE makarovMult :: Matrix3 Integer -> Matrix3 Integer -> Matrix3 Integer #-} --- |Compute the determinant of a matrix.+-- | Compute the determinant of a matrix. det :: Num t => Matrix3 t -> t det Matrix3 {..} =   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' returns a diagonal matrix+instance Fractional t => Fractional (Matrix3 t) where   fromRational = diag . fromRational -  recip a@(Matrix3 {..}) = Matrix3 {+  recip a@Matrix3{..} = Matrix3 {     a11 =  (a22 * a33 - a32 * a23) / d,     a12 = -(a21 * a33 - a23 * a31) / d,     a13 =  (a21 * a32 - a22 * a31) / d,@@ -267,7 +269,8 @@     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".+-- | Convert a list of 9 elements into 'Matrix3'.+-- Reverse conversion can be done by 'Data.Foldable.toList'. fromList :: [t] -> Matrix3 t fromList [a11, a12, a13, a21, a22, a23, a31, a32, a33] = Matrix3 {   a11 = a11,@@ -280,23 +283,30 @@   a32 = a32,   a33 = a33   }-fromList _ = error "The list must contain exactly 9 elements"+fromList _ = error "fromList: input must contain exactly 9 elements" --- |Divide all elements of the matrix by their greatest common+-- | Divide all elements of the matrix by their greatest common -- divisor. This is useful for matrices of projective -- 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+normalize a@Matrix3{..}+  | d <= 1    = a+  | otherwise = fmap (`quot` d) a+  where+    d = go a11 [a12, a13, a21, a22, a23, a31, a32, a33]+    go 1 _ = 1+    go 2 xs = if all even xs then 2 else 1+    go acc [] = acc+    go acc (x : xs) = go (gcd acc x) xs+{-# SPECIALIZE normalize :: Matrix3 Integer -> Matrix3 Integer #-}  instance Pretty t => Pretty (Matrix3 t) where   pretty = vsep . map hsep . pad . fmap pretty where-    pad (Matrix3 {..}) = map (zipWith fill ls) table 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).+-- | Multiplicate a matrix by a column vector. multCol :: Num t => Matrix3 t -> (t, t, t) -> (t, t, t) multCol Matrix3 {..} (a1, a2, a3) = (   a11 * a1 + a12 * a2 + a13 * a3,

@@ -1,11 +1,10 @@ {-| Module      : Math.ExpPairs.MenzerNowak-Description : Asymmetric divisor problem with congruence conditions-Copyright   : (c) Andrew Lelechenko, 2014-2015+Copyright   : (c) Andrew Lelechenko, 2014-2020 License     : GPL-3 Maintainer  : andrew.lelechenko@gmail.com-Stability   : experimental-Portability : POSIX++Asymmetric divisor problem with congruence conditions  Let τ_{a, b}(l_1, k_1; l_2, k_2; n) denote the number of integer (v, w) with v^a w^b = n, v ≡ l_1 (mod k_1), w ≡ l_2 (mod k_2).

@@ -1,20 +1,18 @@ {-| Module      : Math.ExpPairs.RatioInf-Description : Initial exponent pairs-Copyright   : (c) Andrew Lelechenko, 2014-2015+Copyright   : (c) Andrew Lelechenko, 2014-2020 License     : GPL-3 Maintainer  : andrew.lelechenko@gmail.com-Stability   : experimental-Portability : POSIX +Initial exponent pairs.+ Provides a set of initial exponent pairs, consisting-of two points (0, 1), (1\/2, 1\/2) and a triangle with vertices in (1\/6, 2\/3), (2\/13, 35\/52) and (32\/205, 269\/410). The triangle is represented as a list of nodes of a net, covering the triangle.+of two points (0, 1), (1\/2, 1\/2) and a triangle with vertices in (1\/6, 2\/3), (2\/13, 35\/52) and (13\/84, 55\/84). The triangle is represented as a list of nodes of a net, covering the triangle.  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 #-} 

@@ -1,16 +1,11 @@ {-| Module      : Math.ExpPairs.PrettyProcess-Description : Compact representation of process sequences-Copyright   : (c) Andrew Lelechenko, 2015+Copyright   : (c) Andrew Lelechenko, 2015-2020 License     : GPL-3 Maintainer  : andrew.lelechenko@gmail.com-Stability   : experimental-Portability : POSIX  Transforms sequences of 'Process' into most compact (by the means of typesetting) representation using brackets and powers. E. g., AAAABABABA -> A^4(BA)^3.--This module uses memoization extensively. -} {-# LANGUAGE LambdaCase      #-} {-# OPTIONS_GHC -fno-warn-type-defaults #-}@@ -30,7 +25,8 @@  import Math.ExpPairs.ProcessMatrix --- | Compact representation of the sequence of 'Process'.+-- | Compact representation of the sequence of 'Process',+-- using brackets and powers. data PrettyProcess   = Simply [Process]   | Repeat PrettyProcess Int@@ -48,20 +44,21 @@     Repeat xs n  -> parens (pretty xs) <> pretty "^" <> pretty n     Sequence a b -> pretty a <+> pretty b --- | Width of the bracket.+-- | Width of the bracket symbol. bracketWidth :: Int bracketWidth = 4 --- | Width of the subscript-sized character (e. g., power).+-- | Width of the subscript-sized digit (e. g., power). subscriptWidth :: Int subscriptWidth = 4 --- | Width of the processes in typeset+-- | Width of the process symbols. processWidth :: Process -> Int processWidth  A = 10 processWidth BA = 20 --- | Compute the width of the 'PrettyProcess' according to 'bracketWidth', 'subscriptWidth' and 'printedWidth''.+-- | Compute the width of the 'PrettyProcess'+-- according to 'bracketWidth', 'subscriptWidth' and 'printedWidth''. printedWidth :: PrettyProcess -> Int printedWidth = \case   Simply xs             -> sum (map processWidth xs)@@ -75,7 +72,7 @@ annotateWithWidth :: PrettyProcess -> PrettyProcessWithWidth annotateWithWidth p = PPWL p (printedWidth p) --- | Return non-trivial divisors of an argument.+-- | List divisors. divisors :: Int -> [Int] divisors n = ds1 ++ reverse ds2 where   (ds1, ds2) = unzip [ (a, n `div` a) | a <- [1 .. sqrtint n], n `mod` a == 0 ]

@@ -1,16 +1,18 @@ {-| Module      : Math.ExpPairs.Process-Description : Processes of van der Corput-Copyright   : (c) Andrew Lelechenko, 2014-2015+Copyright   : (c) Andrew Lelechenko, 2014-2020 License     : GPL-3 Maintainer  : andrew.lelechenko@gmail.com-Stability   : experimental-Portability : POSIX -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.+Sequences of \( A \)- and \( B \)-processes+of van der Corput's method of exponential sums.+A good reference 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 CPP             #-} {-# LANGUAGE DeriveGeneric   #-}  module Math.ExpPairs.Process@@ -23,8 +25,6 @@   ) where  import GHC.Generics             (Generic)-import Data.Monoid              (Monoid, mempty, mappend)-import Data.Semigroup           (Semigroup, (<>)) import Data.Text.Prettyprint.Doc hiding ((<>))  import Math.ExpPairs.ProcessMatrix@@ -72,12 +72,12 @@ baPath :: Path baPath = Path baMatrix [BA] --- |Apply a projective transformation, defined by 'Path',+-- | Apply a projective transformation, defined by 'Path', -- to a given point in two-dimensional projective space. evalPath :: (Num t) => Path -> (t, t, t) -> (t, t, t) evalPath (Path m _) = evalMatrix m  -- | Count processes in the 'Path'. Note that 'BA' counts--- for one process, not two.+-- for a single process. lengthPath :: Path -> Int lengthPath (Path _ xs) = length xs

@@ -1,17 +1,18 @@ {-| Module      : Math.ExpPairs.ProcessMatrix-Description : Monoidal wrapper for Matrix3-Copyright   : (c) Andrew Lelechenko, 2014-2015+Copyright   : (c) Andrew Lelechenko, 2014-2020 License     : GPL-3 Maintainer  : andrew.lelechenko@gmail.com-Stability   : experimental-Portability : POSIX -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.+Sequences of \( A \)- and \( B \)-processes+of van der Corput's method of exponential sums.+A good reference 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 BangPatterns               #-}-{-# LANGUAGE CPP                        #-} {-# LANGUAGE DeriveGeneric              #-} {-# LANGUAGE GeneralizedNewtypeDeriving #-} @@ -23,25 +24,26 @@   , evalMatrix   ) where -import Data.Monoid           (Monoid, mempty, mappend)-import Data.Semigroup        (Semigroup, (<>)) import GHC.Generics          (Generic (..)) import Data.Text.Prettyprint.Doc  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'.+-- | Since \( B \)-process is+-- <https://en.wikipedia.org/wiki/Involution_(mathematics) involutive>,+-- a sequence of \( A \)- and \( B \)-processes can be rewritten as a sequence+-- of 'A' and 'BA'. data Process-  -- | /A/-process+  -- | \( A \)-process   = A-  -- | /BA/-process+  -- | \( BA \)-process   | BA   deriving (Eq, Show, Read, Ord, Enum, Generic)  instance Pretty Process where   pretty = pretty . show --- | Sequence of processes, represented as a matrix 3x3.+-- | Sequence of processes, represented as a matrix \( 3 \times 3 \). newtype ProcessMatrix = ProcessMatrix (Matrix3 Integer)   deriving (Eq, Num, Show, Pretty) @@ -56,15 +58,15 @@ process2matrix  A = ProcessMatrix $ Matrix3 1 0 0 1 1 1  2 0 2 process2matrix BA = ProcessMatrix $ Matrix3 0 1 0 2 0 1  2 0 2 --- | Return process matrix for 'A'-process.+-- | Return process matrix for \( A \)-process. aMatrix :: ProcessMatrix aMatrix = process2matrix A --- | Return process matrix for 'BA'-process.+-- | Return process matrix for \( BA \)-process. baMatrix :: ProcessMatrix baMatrix = process2matrix BA --- |Apply a projective transformation, defined by 'Path',+-- | 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) = multCol (fmap fromInteger m)

@@ -1,16 +1,13 @@ {-| Module      : Math.ExpPairs.RatioInf-Description : Rational numbers with infinities-Copyright   : (c) Andrew Lelechenko, 2014-2015+Copyright   : (c) Andrew Lelechenko, 2014-2020 License     : GPL-3 Maintainer  : andrew.lelechenko@gmail.com-Stability   : experimental-Portability : POSIX -Provides types and necessary instances for rational numbers, extended with infinite values. Just use 'RationalInf' instead of 'Rational' from "Data.Ratio".+Rational numbers extended with infinities. -} -{-# OPTIONS_GHC -fno-warn-orphans #-}+{-# LANGUAGE Safe #-}  module Math.ExpPairs.RatioInf   ( RatioInf (..)@@ -20,15 +17,12 @@ import Data.Ratio (Ratio, numerator, denominator) import Data.Text.Prettyprint.Doc --- |Extends a rational type with positive and negative+-- | Extend 'Ratio' @t@ with \( \pm \infty  \) positive and negative -- infinities. data RatioInf t-  -- |Negative infinity-  = InfMinus-  -- |Finite value-  | Finite !(Ratio t)-  -- |Positive infinity-  | InfPlus+  = InfMinus          -- ^ \( - \infty  \)+  | Finite !(Ratio t) -- ^ Finite value+  | InfPlus           -- ^ \( + \infty  \)   deriving (Eq, Ord, Show)  -- |Arbitrary-precision rational numbers with positive and negative

@@ -0,0 +1,17 @@+module Main where++import qualified Math.ExpPairs.Matrix3 as M3+import Gauge.Main++testm3 :: Int -> M3.Matrix3 Integer+testm3 k = M3.fromList $ map (100*10^k `div`) [100..108]++compareMults :: Int -> Benchmark+compareMults k = bgroup (show k)+  [ bench "vanillaMult"  $ nf (\x -> x * x)                 (testm3 k)+  , bench "makarovMult"  $ nf (\x -> x `M3.makarovMult` x)  (testm3 k)+  , bench "ladermanMult" $ nf (\x -> x `M3.ladermanMult` x) (testm3 k)+  ]++main :: IO ()+main = defaultMain $ map compareMults [400,450..800]

@@ -0,0 +1,289 @@+{-# LANGUAGE OverloadedStrings, QuasiQuotes #-}+module Main where++import Data.Set (Set)+import qualified Data.Set as S+import Data.IntSet (IntSet)+import qualified Data.IntSet as IS+import Data.IntMap.Strict (IntMap)+import qualified Data.IntMap.Strict as IM+import Data.Bimap (Bimap)+import qualified Data.Bimap as B++import Data.String+import Text.RawString.QQ     (r)++import Data.List (sortOn, subsequences, mapAccumL)+import Data.Monoid+import Data.Foldable+import Data.Maybe+import Data.Ord+import Data.Function++import Data.Text.Prettyprint.Doc hiding ((<>), group)++import Prelude hiding (foldl, foldl1, maximum, mapM_, minimum, sum, concat)++data Sign = Minus | Plus+  deriving (Eq, Ord, Show)++data Form = Form String (Set (Sign, String))+  deriving (Eq, Ord, Show)++type Factor = Int+type Expr   = IntSet+data Forms = Forms (Bimap String Int) (IntMap Expr)+  deriving (Eq, Show)++instance IsString Forms where+  fromString = foldl addForm (Forms B.empty mempty) . map parseForm . filter (not . null) . lines++parseForm :: String -> Form+parseForm xs = case words xs of+  (lhs : "=" : ys) -> Form lhs (f ys) where+    f [] = mempty+    f ("+":z:zs) = S.singleton (Plus, z) <> f zs+    f ("-":z:zs) = S.singleton (Minus, z) <> f zs+    f (z:zs) = S.singleton (Plus, z) <> f zs+  as -> error (show as)++addForm :: Forms -> Form -> Forms+addForm (Forms ids exprs) (Form lhs rhs) =+  Forms (B.insert lhs lhsId ids') (IM.insert lhsId rhsExpr exprs) where+    (ids', rhsExpr) = foldl appendVar (ids, mempty) rhs+    lhsId = nextId ids'++appendVar :: (Bimap String Int, Expr) -> (Sign, String) -> (Bimap String Int, Expr)+appendVar (ids, expr) (sign, str) = (ids', expr') where+  maybeId = B.lookup str ids+  (ids', newId) = case maybeId of+    Nothing -> (B.insert str t ids, t) where+      t = nextId ids+    Just i  -> (ids, i)+  expr' = IS.insert (if sign==Plus then newId else negate newId) expr++nextId :: (Ord a, Num b, Enum b) => Bimap a b -> b+nextId bm+  | B.null bm = 1+  | otherwise = succ . fst . B.findMaxR $ bm++instance Pretty Forms where+  pretty (Forms ids exprs) = vsep (IM.elems $ IM.mapWithKey prettyExpr exprs) where+    toStr i = fromJust $ B.lookupR i ids+    prettyExpr lhsId expr = pretty (toStr lhsId) <+> equals <+> hsep (map prettyVar (IS.toDescList expr))+    prettyVar n = pretty (if n < 0 then '-' else '+') <+> pretty (toStr (abs n))+++flipSign :: Expr -> Expr+flipSign = IS.map negate++nonTrivial :: [a] -> Bool+nonTrivial []  = False+nonTrivial [_] = False+nonTrivial _   = True++-- | Input must be in ascending order+absExpr :: (Num a, Ord a) => [a] -> [a]+absExpr [] = []+absExpr [x] = [x]+absExpr xxs@(x : xs) = if negate x < last xs then xxs else map negate xxs++nub :: Ord a => [a] -> [a]+nub = toList . S.fromList++suspicious :: Forms -> [IntSet]+suspicious fs@(Forms _ exprs) = take 8 $ sortOn (substitutionWeight fs) subexprs+  where+    subexprs :: [Expr]+    subexprs = toList $ fold $ snd $ mapAccumL doExpr mempty (toList exprs)++    doExpr :: [IntSet] -> Expr -> ([IntSet], Set Expr)+    doExpr acc expr = (expr : acc+      , S.fromList+      $ map (IS.fromList . absExpr)+      $ filter nonTrivial+      $ concatMap (subsequences . IS.toList)+      $ filter ((> 1) . IS.size)+      $ map (IS.intersection expr) acc ++ map (IS.intersection (flipSign expr)) acc+      )++substitutionWeight:: Forms -> Expr -> Int+substitutionWeight (Forms _ exprs) newExpr = size * count where+  newExprNeg = flipSign newExpr+  size = IS.size newExpr - 1+  count = 1 - IM.size (IM.filter (\expr -> newExpr    `IS.isSubsetOf` expr || newExprNeg `IS.isSubsetOf` expr) exprs)++substitute :: Char -> Forms -> Expr -> Forms+substitute ch (Forms ids exprs) newExpr = Forms ids' exprs' where+  lhsId  = nextId ids+  lhs    = ch : show lhsId+  ids'   = B.insert lhs lhsId ids+  exprs' = IM.insert lhsId newExpr (IM.map eliminate exprs)++  newExprNeg = flipSign newExpr+  eliminate expr+    | newExpr    `IS.isSubsetOf` expr = IS.insert lhsId (expr `IS.difference` newExpr)+    | newExprNeg `IS.isSubsetOf` expr = IS.insert (-lhsId) (expr `IS.difference` newExprNeg)+    | otherwise                       = expr++weight :: Forms -> Int+weight (Forms _ exprs) = getSum $ foldMap (Sum . pred . IS.size) exprs++optimize :: Char -> Forms -> Forms+optimize ch fs = case suspicious fs of+  []   -> fs+  susp -> minimumBy (comparing weight) . map (optimize ch . substitute ch fs) $ susp++main :: IO ()+main = do+  let before =+        [ ('t', makarov1)+        , ('u', makarov2)+        , ('v', makarov3)+        , ('t', laderman1)+        , ('u', laderman2)+        , ('v', laderman3)+        ]+  let weightBefore = map (weight . snd) before+  let after = map (uncurry optimize) before+  let weightAfter = map weight after+  mapM_ (print . pretty) after+  putStrLn $ show weightBefore ++ " = " ++ show (sum weightBefore)+  putStrLn $ show weightAfter  ++ " = " ++ show (sum weightAfter)++makarov1 :: Forms+makarov1 = [r|+l1  = a3 + c1 - c2+l2  = a2 + b1 + b2+l3  = a2 + b1 + b3+l4  = a3 - c2 - c3+l5  = a1 - c1 + c2+l6  = a1 + b1 + b2+l7  = a1 + b1 + b3 + c2 + c3+l8  = a2+l9  = a3+l10 = b1+l11 = c2+l12 = c1 - c2+l13 = b1 + b2+l14 = a2 + b1+l15 = b2+l16 = a3 - c2+l17 = c2+l18 = b3 - c2 - c3+l19 = c1 + c3 - b1 - b3+l20 = b1 + b3+l21 = c2 + c3+l22 = c2 + c3 - b1 - b3+|]+++makarov2 :: Forms+makarov2 = [r|+r1 =  k1 + k7 - k8 + k9+r2 =  k2 - k4 + k5 - k6+r3 =  k3 - k4 + k5 - k6+r4 =  k3 - k7 + k8 - k9+r5 =  k1+r6 =  k2+r7 =  k3+r8 =  k1 + k4 - k5 + k6+r9 =  k2 + k7 - k8 + k9+r10 =  k4+r11 =  k7+r12 =  k1 + k7+r13 =  k4 - k2+r14 =  k4 - k5 + k6+r15 =  k6+r16 =  k7 - k8 + k9+r17 =  k8+r18 =  k6+r19 =  k8+r20 =  k4 - k3 + k6 + k8+r21 =  k3 + k6 - k7 + k8+r22 =  k6 + k8+|]++makarov3 :: Forms+makarov3 = [r|+c11 = m5 + m10 + m11 + m12+c12 = m8 + m10 - m14 + m17 - m18 + m19 - m22+c13 = m1 - m11 - m12 - m16 + m17 - m18 + m19 - m22+c21 = m6 - m10 + m11 + m13+c22 = m2 - m10 + m13 + m14 + m15 + m17+r23 = m9 - m11 + m15 - m16 + m17+c31 = m7 - m10 - m11 + m20 - m21 + m22+c32 = m3 - m10 + m14 - m17 + m18 + m20 + m22+c33 = m4 + m11 + m16 - m17 + m18 + m21+|]++laderman1 :: Forms+laderman1 = [r|+l1 = a11 + a12 + a13 - a21 - a22 - a32 - a33+l2 = a11 - a21+l3 = a22+l4 = - a11 + a21 + a22+l5 = a21 + a22+l6 = a11+l7 = - a11 + a31 + a32+l8 = - a11 + a31+l9 = a31 + a32+l10 = a11 + a12 + a13 - a22 - a23 - a31 - a32+l11 = a32+l12 = - a13 + a32 + a33+l13 = a13 - a33+l14 = a13+l15 = a32 + a33+l16 = - a13 + a22 + a23+l17 = a13 - a23+l18 = a22 + a23+l19 = a12+l20 = a23+l21 = a21+l22 = a31+l23 = a33+|]++laderman2 :: Forms+laderman2 = [r|+r1 = b22+r2 = - b12 + b22+r3 = - b11 + b12 + b21 - b22 - b23 - b31 + b33+r4 = b11 - b12 + b22+r5 = - b11 + b12+r6 = b11+r7 = b11 - b13 + b23+r8 = b13 - b23+r9 = - b11 + b13+r10 = b23+r11 = - b11 + b13 + b21 - b22 - b23 - b31 + b32+r12 = b22 + b31 - b32+r13 = b22 - b32+r14 = b31+r15 = - b31 + b32+r16 = b23 + b31 - b33+r17 = b23 - b33+r18 = - b31 + b33+r19 = b21+r20 = b32+r21 = b13+r22 = b12+r23 = b33+|]++laderman3 :: Forms+laderman3 = [r|+c11 = m6 + m14 + m19+c12 = m1 + m4 + m5 + m6 + m12 + m14 + m15+c13 = m6 + m7 + m9 + m10 + m14 + m16 + m18+c21 = m2 + m3 + m4 + m6 + m14 + m16 + m17+c22 = m2 + m4 + m5 + m6 + m20+c23 = m14 + m16 + m17 + m18 + m21+c31 = m6 + m7 + m8 + m11 + m12 + m13 + m14+c32 = m12 + m13 + m14 + m15 + m22+c33 = m6 + m7 + m8 + m9 + m23+|]++++

@@ -1,5 +1,5 @@ name:                exp-pairs-version:             0.2.0.0+version:             0.2.1.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@@ -11,6 +11,8 @@ build-type:          Simple extra-source-files:  CHANGELOG.md, tests/*.txt cabal-version:       >=1.10+tested-with:+  GHC ==8.4.4 GHC ==8.6.5 GHC ==8.8.3 GHC ==8.10.1  source-repository head   type:     git@@ -28,7 +30,7 @@                        Math.ExpPairs.PrettyProcess,                        Math.ExpPairs.ProcessMatrix,                        Math.ExpPairs.RatioInf-  build-depends:       base >=4 && <5,+  build-depends:       base >=4.11 && <5,                        ghc-prim,                        prettyprinter,                        deepseq >=1.3,@@ -52,16 +54,40 @@                        Process,                        PrettyProcess,                        RatioInf-  build-depends:       base >=4 && <5,+  build-depends:       base,                        tasty >=0.7,                        tasty-quickcheck,                        tasty-smallcheck,                        tasty-hunit,                        QuickCheck >=2.9,-                       smallcheck >=0.2.1,+                       smallcheck >=1.2,                        exp-pairs,                        matrix >=0.1,                        random   hs-source-dirs:      tests   default-language:    Haskell2010   ghc-options:         -Wall -fno-warn-type-defaults++benchmark matrix-bench+  build-depends:+    base,+    exp-pairs,+    gauge+  type: exitcode-stdio-1.0+  main-is: BenchMatrix.hs+  default-language: Haskell2010+  hs-source-dirs: auxiliary+  ghc-options: -Wall -fno-warn-type-defaults++benchmark optimize-sum+  build-depends:+    base,+    bimap,+    containers,+    prettyprinter,+    raw-strings-qq+  type: exitcode-stdio-1.0+  main-is: OptimizeSum.hs+  default-language: Haskell2010+  hs-source-dirs: auxiliary+  ghc-options: -Wall -fno-warn-type-defaults

@@ -13,7 +13,7 @@ fetchRandomLines n filename = do   etalon <- readFile filename   gen <- newStdGen-  let items = lines etalon -- (take n . unsort gen . lines) etalon+  let items = (take n . unsort gen . lines) etalon   let tests = map (map read . words) items   return tests 

@@ -1,14 +1,17 @@+{-# LANGUAGE DeriveGeneric         #-}+{-# LANGUAGE FlexibleInstances     #-}+{-# LANGUAGE MultiParamTypeClasses #-}+ {-# OPTIONS_GHC -fno-warn-orphans #-}-{-# LANGUAGE FlexibleInstances, MultiParamTypeClasses, DeriveGeneric, CPP #-}+ module Instances (Ratio01 (..), Positive (..), Sorted(..)) where  import Test.QuickCheck (Arbitrary(..), Gen, genericShrink, suchThat, vectorOf) import Test.SmallCheck.Series import Control.Applicative import Control.Monad-#if __GLASGOW_HASKELL__ < 710-import Data.Foldable-#endif+import Data.List (sort)+import Data.Ratio import GHC.Generics          (Generic (..))  import Math.ExpPairs.LinearForm@@ -54,20 +57,24 @@ instance Monad m => Serial m Process where   series = cons0 A \/ cons0 BA +-- | Does not have a 'Num' instance! 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))+instance (Ord t, Integral t, Arbitrary t) => Arbitrary (Ratio01 (Ratio t)) where+  arbitrary = do+    denom <- arbitrary `suchThat` (> 0)+    numer <- arbitrary+    pure $ Ratio01 $ (numer `mod` (denom + 1)) % denom   shrink (Ratio01 y) = Ratio01 <$> filter (\x -> 0 <= x && x <= 1) (shrink y) -instance (Ord t, Fractional t, Serial m t) => Serial m (Ratio01 t) where+instance (Ord t, Integral t, Serial m t) => Serial m (Ratio01 (Ratio 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+instance (Ord t, Integral t, Arbitrary t) => Arbitrary (InitPair' (Ratio t)) where   arbitrary = f <$> liftM2 (,) arbitrary arbitrary where     f :: (Ord t, Fractional t) => (Ratio01 t, Ratio01 t) -> InitPair' t     f (Ratio01 x, Ratio01 y)@@ -78,7 +85,7 @@         y' = y*(1-x)   shrink = genericShrink -instance (Ord t, Fractional t, Serial m t) => Serial m (InitPair' t) where+instance (Ord t, Integral t, Serial m t) => Serial m (InitPair' (Ratio t)) where   series = cons0 Corput01 \/ cons0 Corput12 \/ mseries     where       mseries = do@@ -97,31 +104,6 @@ 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 (,,,,,)- newtype Sorted t = Sorted t   deriving (Show, Generic) @@ -131,30 +113,82 @@ 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 (Num t, Ord t, Arbitrary t) => Arbitrary (Sorted (t, t, t)) where+  arbitrary = do+    a  <- arbitrary+    ab <- arbitrary `suchThat` (>= 0)+    bc <- arbitrary `suchThat` (>= 0)+    let b = a + ab; c = b + bc+    pure $ Sorted (a, b, c) -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 (Num t, Ord t, Serial m t) => Serial m (Sorted (t, t, t)) where+  series = do+    a  <- series+    ab <- series `suchThatSerial` (>= 0)+    bc <- series `suchThatSerial` (>= 0)+    let b = a + ab; c = b + bc+    pure $ Sorted (a, b, c) -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 (Num t, Ord t, Arbitrary t) => Arbitrary (Sorted (t, t, t, t)) where+  arbitrary = do+    a  <- arbitrary+    ab <- arbitrary `suchThat` (>= 0)+    bc <- arbitrary `suchThat` (>= 0)+    cd <- arbitrary `suchThat` (>= 0)+    let b = a + ab; c = b + bc; d = c + cd+    pure $ Sorted (a, b, c, d)+  shrink (Sorted (aa, bb, cc, dd))+    = map ((\[a, b, c, d] -> Sorted (a, b, c, d)) . sort)+    $ filter ((== 4) . length)+    $ shrink [aa, bb, cc, dd] -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 (Num t, Ord t, Serial m t) => Serial m (Sorted (t, t, t, t)) where+  series = do+    a  <- series+    ab <- series `suchThatSerial` (>= 0)+    bc <- series `suchThatSerial` (>= 0)+    cd <- series `suchThatSerial` (>= 0)+    let b = a + ab; c = b + bc; d = c + cd+    pure $ Sorted (a, b, c, d) -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 (Num t, Ord t, Arbitrary t) => Arbitrary (Sorted (t, t, t, t, t, t)) where+  arbitrary = do+    a  <- arbitrary+    ab <- arbitrary `suchThat` (>= 0)+    bc <- arbitrary `suchThat` (>= 0)+    cd <- arbitrary `suchThat` (>= 0)+    de <- arbitrary `suchThat` (>= 0)+    ef <- arbitrary `suchThat` (>= 0)+    let b = a + ab; c = b + bc; d = c + cd; e = d + de; f = e + ef+    pure $ Sorted (a, b, c, d, 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))+instance (Num t, Ord t, Serial m t) => Serial m (Sorted (t, t, t, t, t, t)) where+  series = do+    a  <- series+    ab <- series `suchThatSerial` (>= 0)+    bc <- series `suchThatSerial` (>= 0)+    cd <- series `suchThatSerial` (>= 0)+    de <- series `suchThatSerial` (>= 0)+    ef <- series `suchThatSerial` (>= 0)+    let b = a + ab; c = b + bc; d = c + cd; e = d + de; f = e + ef+    pure $ Sorted (a, b, c, d, e, f) +instance (Num t, Ord t, Arbitrary t) => Arbitrary (Sorted (t, t, t, t, t, t, t, t)) where+  arbitrary = do+    a  <- arbitrary+    ab <- arbitrary `suchThat` (>= 0)+    bc <- arbitrary `suchThat` (>= 0)+    cd <- arbitrary `suchThat` (>= 0)+    de <- arbitrary `suchThat` (>= 0)+    ef <- arbitrary `suchThat` (>= 0)+    fg <- arbitrary `suchThat` (>= 0)+    gh <- arbitrary `suchThat` (>= 0)+    let b = a + ab; c = b + bc; d = c + cd; e = d + de; f = e + ef; g = f + fg; h = g + gh+    pure $ Sorted (a, b, c, d, e, f, g, h)+  shrink (Sorted (aa, bb, cc, dd, ee, ff, gg, hh))+    = map ((\[a, b, c, d, e, f, g, h] -> Sorted (a, b, c, d, e, f, g, h)) . sort)+    $ filter ((== 8) . length)+    $ shrink [aa, bb, cc, dd, ee, ff, gg, hh]  instance Arbitrary Path where   arbitrary = foldMap (\x -> if x then aPath else baPath) <$> (arbitrary :: Gen [Bool])

@@ -59,20 +59,24 @@   a = fromMinus3To3 a'  testZetaReverse1 :: Ratio01 Rational -> Bool-testZetaReverse1 (Ratio01 s') = if t <= s + 2.1e-2 && s <= t + 2e-3 then True else trace (show $ fromRational $ s-t) False where-  s = fromHalfToOne s'-  zs = zetaOnS s-  t = toRational $ optimalValue $ reverseZetaOnS $ toRational $ optimalValue zs+testZetaReverse1 (Ratio01 s') = t <= s + 2.1e-2 && s <= t + 2e-3 ||+  trace (show $ fromRational $ s-t) False+  where+    s = fromHalfToOne s'+    zs = zetaOnS s+    t = toRational $ optimalValue $ reverseZetaOnS $ toRational $ optimalValue zs  testZetaReverse2 :: Ratio01 Rational -> Bool-testZetaReverse2 (Ratio01 s') = if t <= s + 1e-10 && s <= t + 4e-3 then True else trace (show $ fromRational $ s-t) False where-  s = s' * 32 / 205-  zs = reverseZetaOnS s-  t = toRational $ optimalValue $ zetaOnS $ toRational $ optimalValue zs+testZetaReverse2 (Ratio01 s') = t <= s + 1e-10 && s <= t + 4e-3 ||+  trace (show $ fromRational $ s-t) False+  where+    s = s' * 32 / 205+    zs = reverseZetaOnS s+    t = toRational $ optimalValue $ zetaOnS $ toRational $ optimalValue zs  testMOnSReverse1 :: Ratio01 Rational -> Bool-testMOnSReverse1 (Ratio01 s') =-  if t <= s + 4e-2 && s <= t + 1e-3 then True else trace (show $ fromRational $ s-t) False+testMOnSReverse1 (Ratio01 s') = t <= s + 4e-2 && s <= t + 1.4e-3 ||+  trace (show $ fromRational $ s-t) False   where     s = fromHalfToOne s'     zs = mOnS s@@ -82,7 +86,8 @@ testMOnSReverse2 (Ratio01 s')   =  s' == 0   || t' == InfPlus || t' == InfMinus-  || if recip t <= recip s + 1e-3 && recip s <= recip t + 1e-3 then True else trace (show $ fromRational $ recip s - recip t) False+  || recip t <= recip s + 1e-3 && recip s <= recip t + 1e-3+  || trace (show $ fromRational $ recip s - recip t) False     where       s  = 4 * recip s'       zs = reverseMOnS 1e-3 (Finite s)@@ -90,17 +95,20 @@       t  = toRational t'  testMBigOnHalfReverse1 :: Positive Rational -> Bool-testMBigOnHalfReverse1 (Positive s') = if recip t <= recip s + 2e-3 && recip s <= recip t + 1e-10 then True else trace (show $ fromRational $ recip s - recip t) False where-  s = s' + 4-  zs = mBigOnHalf s-  t = toRational $ optimalValue $ reverseMBigOnHalf $ toRational $ optimalValue zs+testMBigOnHalfReverse1 (Positive s') = recip t <= recip s + 2e-3 && recip s <= recip t + 1e-10 ||+  trace (show $ fromRational $ recip s - recip t) False+  where+    s = s' + 4+    zs = mBigOnHalf s+    t = toRational $ optimalValue $ reverseMBigOnHalf $ toRational $ optimalValue zs  testMBigOnHalfReverse2 :: Positive Rational -> Bool-testMBigOnHalfReverse2 (Positive s') = if recip t <= recip s + 2e-3 && recip s <= recip t + 1e-10 then True else trace (show $ fromRational $ recip s - recip t) False where-  s = s' + 1-  zs = reverseMBigOnHalf s-  t = toRational $ optimalValue $ mBigOnHalf $ toRational $ optimalValue zs-+testMBigOnHalfReverse2 (Positive s') = recip t <= recip s + 2e-3 && recip s <= recip t + 1e-10 ||+  trace (show $ fromRational $ recip s - recip t) False+  where+    s = s' + 1+    zs = reverseMBigOnHalf s+    t = toRational $ optimalValue $ mBigOnHalf $ toRational $ optimalValue zs  etalonZetaOnS :: Integer -> Integer -> Integer -> Integer -> Bool etalonZetaOnS a b c d = Finite (c%d) >= optimalValue (zetaOnS $ a%b)

@@ -37,12 +37,29 @@  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))+  = 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)) +testAbcdMonotonic :: Sorted (Positive Integer, Positive Integer, Positive Integer, Positive Integer, Positive Integer, Positive Integer, Positive Integer, Positive Integer) -> Bool+testAbcdMonotonic (Sorted (Positive a, Positive e, Positive b, Positive f, Positive c, Positive g, Positive d, Positive h))+  = (a == e && b == f && c == g && d == h) || theoremAbcd `elem` [HeathBrown, Kr1992_32] || zabcd >= zefgh+    where+      (theoremAbcd, resultAbcd) = tauabcd a b c d+      zabcd = optimalValue resultAbcd+      zefgh = optimalValue $ snd $ tauabcd e f g h++testAbcdCompareLow :: Sorted (Positive Integer, Positive Integer, Positive Integer, Positive Integer) -> Bool+testAbcdCompareLow (Sorted (Positive a, Positive b, Positive c, Positive d))+  = optimalValue (snd $ tauabcd a b c d) >= Finite (1 % (a + b + c + d))++-- | Kratzel1988, Eq. (6.31)+testAbcdCompareHigh :: Sorted (Positive Integer, Positive Integer, Positive Integer, Positive Integer) -> Bool+testAbcdCompareHigh (Sorted (Positive a, Positive b, Positive c, Positive d))+  = d >= a + b + c || optimalValue (snd $ tauabcd a b c d) <= Finite ((a + b + c) % (a * (a + b + c + d)))+ etalonTauab :: Integer -> Integer -> Integer -> Integer -> Bool etalonTauab a b c d = Finite (c % d) >= (optimalValue . snd) (tauab a b) @@ -55,6 +72,13 @@     (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 "tauabcd compare with 1/(a+b+c+d)" testAbcdCompareLow+  , QC.testProperty "tauabcd compare with 1/(a+b+c+d)" testAbcdCompareLow+  , SC.testProperty "tauabcd compare with (6.31)" testAbcdCompareHigh+  , QC.testProperty "tauabcd compare with (6.31)" testAbcdCompareHigh+  , QC.testProperty "tauabcd monotonic" testAbcdMonotonic+   , 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@@ -62,6 +86,7 @@   , 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

@@ -1,5 +1,6 @@ module Matrix3 where +import Data.Foldable import qualified Data.Matrix as M import qualified Math.ExpPairs.Matrix3 as M3 @@ -9,10 +10,10 @@ import Instances ()  toM :: M3.Matrix3 a -> M.Matrix a-toM = M.fromList 3 3 . M3.toList+toM = M.fromList 3 3 . toList  toM3 :: M.Matrix a -> M3.Matrix3 a-toM3 = M3.fromList . M.toList+toM3 = M3.fromList . 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