exp-pairs 0.1.2.0 → 0.1.3.0
raw patch · 13 files changed
+581/−292 lines, 13 filesdep +deepseqdep +generic-derivingdep +wl-pprintdep ~basePVP: major bump suggested
API removals or changes: PVP suggests a major version bump
Dependencies added: deepseq, generic-deriving, wl-pprint
Dependency ranges changed: base
API changes (from Hackage documentation)
- Math.ExpPairs.Matrix3: instance Fractional t => Fractional (Matrix3 t)
- Math.ExpPairs.Matrix3: instance Num t => Num (Matrix3 t)
- Math.ExpPairs.Matrix3: prettyMatrix :: Show t => Matrix3 t -> String
- Math.ExpPairs.Process: A :: Process
- Math.ExpPairs.Process: BA :: Process
- Math.ExpPairs.Process: instance Enum Process
- Math.ExpPairs.Process: instance Eq Process
- Math.ExpPairs.Process: instance Memoizable Process
- Math.ExpPairs.Process: instance Ord Process
- Math.ExpPairs.Process: instance Read Process
- Math.ExpPairs.Process: instance Show Process
+ Math.ExpPairs.Ivic: checkAbscissa :: [(Rational, Rational)] -> Rational -> Bool
+ Math.ExpPairs.Ivic: findMinAbscissa :: Rational -> [(Rational, Rational)] -> Rational
+ Math.ExpPairs.Ivic: mBigOnHalf :: Rational -> OptimizeResult
+ Math.ExpPairs.Ivic: reverseMBigOnHalf :: Rational -> OptimizeResult
+ Math.ExpPairs.Ivic: reverseMOnS :: Rational -> RationalInf -> Rational
+ Math.ExpPairs.LinearForm: instance Bounded IneqType
+ Math.ExpPairs.LinearForm: instance Constructor C1_0Constraint
+ Math.ExpPairs.LinearForm: instance Constructor C1_0LinearForm
+ Math.ExpPairs.LinearForm: instance Constructor C1_0RationalForm
+ Math.ExpPairs.LinearForm: instance Datatype D1Constraint
+ Math.ExpPairs.LinearForm: instance Datatype D1LinearForm
+ Math.ExpPairs.LinearForm: instance Datatype D1RationalForm
+ Math.ExpPairs.LinearForm: instance Enum IneqType
+ Math.ExpPairs.LinearForm: instance Eq t => Eq (Constraint t)
+ Math.ExpPairs.LinearForm: instance Eq t => Eq (RationalForm t)
+ Math.ExpPairs.LinearForm: instance Foldable Constraint
+ Math.ExpPairs.LinearForm: instance Foldable LinearForm
+ Math.ExpPairs.LinearForm: instance Foldable RationalForm
+ Math.ExpPairs.LinearForm: instance Functor Constraint
+ Math.ExpPairs.LinearForm: instance Functor LinearForm
+ Math.ExpPairs.LinearForm: instance Functor RationalForm
+ Math.ExpPairs.LinearForm: instance Generic (Constraint t)
+ Math.ExpPairs.LinearForm: instance Generic (LinearForm t)
+ Math.ExpPairs.LinearForm: instance Generic (RationalForm t)
+ Math.ExpPairs.LinearForm: instance NFData t => NFData (Constraint t)
+ Math.ExpPairs.LinearForm: instance NFData t => NFData (LinearForm t)
+ Math.ExpPairs.LinearForm: instance NFData t => NFData (RationalForm t)
+ Math.ExpPairs.LinearForm: instance Ord IneqType
+ Math.ExpPairs.Matrix3: instance (Fractional t, Ord t) => Fractional (Matrix3 t)
+ Math.ExpPairs.Matrix3: instance (Num t, Ord t) => Num (Matrix3 t)
+ Math.ExpPairs.Matrix3: instance NFData t => NFData (Matrix3 t)
+ Math.ExpPairs.Matrix3: instance NFData t => NFData (Vector3 t)
+ Math.ExpPairs.Matrix3: ladermanMult :: Num t => Matrix3 t -> Matrix3 t -> Matrix3 t
+ Math.ExpPairs.Matrix3: makarovMult :: Num t => Matrix3 t -> Matrix3 t -> Matrix3 t
+ Math.ExpPairs.PrettyProcess: data PrettyProcess
+ Math.ExpPairs.PrettyProcess: instance Memoizable PrettyProcess
+ Math.ExpPairs.PrettyProcess: instance Pretty PrettyProcess
+ Math.ExpPairs.PrettyProcess: instance Show PrettyProcess
+ Math.ExpPairs.PrettyProcess: prettify :: [Process] -> PrettyProcess
+ Math.ExpPairs.PrettyProcess: uglify :: PrettyProcess -> [Process]
+ Math.ExpPairs.Process: Path :: !ProcessMatrix -> ![Process] -> Path
+ Math.ExpPairs.Process: instance Constructor C1_0Path
+ Math.ExpPairs.Process: instance Datatype D1Path
+ Math.ExpPairs.Process: instance Generic Path
+ Math.ExpPairs.ProcessMatrix: A :: Process
+ Math.ExpPairs.ProcessMatrix: BA :: Process
+ Math.ExpPairs.ProcessMatrix: aMatrix :: ProcessMatrix
+ Math.ExpPairs.ProcessMatrix: baMatrix :: ProcessMatrix
+ Math.ExpPairs.ProcessMatrix: data Process
+ Math.ExpPairs.ProcessMatrix: data ProcessMatrix
+ Math.ExpPairs.ProcessMatrix: evalMatrix :: Num t => ProcessMatrix -> (t, t, t) -> (t, t, t)
+ Math.ExpPairs.ProcessMatrix: instance Enum Process
+ Math.ExpPairs.ProcessMatrix: instance Eq Process
+ Math.ExpPairs.ProcessMatrix: instance Eq ProcessMatrix
+ Math.ExpPairs.ProcessMatrix: instance Memoizable Process
+ Math.ExpPairs.ProcessMatrix: instance Monoid ProcessMatrix
+ Math.ExpPairs.ProcessMatrix: instance Num ProcessMatrix
+ Math.ExpPairs.ProcessMatrix: instance Ord Process
+ Math.ExpPairs.ProcessMatrix: instance Read Process
+ Math.ExpPairs.ProcessMatrix: instance Show Process
+ Math.ExpPairs.ProcessMatrix: instance Show ProcessMatrix
- Math.ExpPairs: Finite :: (Ratio t) -> RatioInf t
+ Math.ExpPairs: Finite :: !(Ratio t) -> RatioInf t
- Math.ExpPairs.Matrix3: det :: Num t => Matrix3 t -> t
+ Math.ExpPairs.Matrix3: det :: (Num t, Ord t) => Matrix3 t -> t
- Math.ExpPairs.RatioInf: Finite :: (Ratio t) -> RatioInf t
+ Math.ExpPairs.RatioInf: Finite :: !(Ratio t) -> RatioInf t
Files
- Math/ExpPairs.hs +21/−9
- Math/ExpPairs/Ivic.hs +37/−73
- Math/ExpPairs/Kratzel.hs +10/−5
- Math/ExpPairs/LinearForm.hs +53/−34
- Math/ExpPairs/Matrix3.hs +170/−28
- Math/ExpPairs/MenzerNowak.hs +4/−2
- Math/ExpPairs/Pair.hs +8/−2
- Math/ExpPairs/PrettyProcess.hs +125/−0
- Math/ExpPairs/Process.hs +33/−99
- Math/ExpPairs/ProcessMatrix.hs +61/−0
- Math/ExpPairs/RatioInf.hs +48/−36
- exp-pairs.cabal +7/−2
- tests/Tests.hs +4/−2
Math/ExpPairs.hs view
@@ -14,12 +14,28 @@ A set of useful applications can be found in "Math.ExpPairs.Ivic", "Math.ExpPairs.Kratzel" and "Math.ExpPairs.MenzerNowak". -}-module Math.ExpPairs (optimize, OptimizeResult, optimalValue, optimalPair, optimalPath, simulateOptimize, simulateOptimize', LinearForm (..), RationalForm (..), IneqType (..), Constraint (..), InitPair, Path, RatioInf (..), RationalInf) where+module Math.ExpPairs+ ( optimize+ , OptimizeResult+ , optimalValue+ , optimalPair+ , optimalPath+ , simulateOptimize+ , simulateOptimize'+ , LinearForm (..)+ , RationalForm (..)+ , IneqType (..)+ , Constraint (..)+ , InitPair+ , Path+ , RatioInf (..)+ , RationalInf+ ) where -import Data.Ratio-import Data.Ord-import Data.List-import Data.Monoid+import Data.Ratio ((%), numerator, denominator)+import Data.Ord (comparing)+import Data.List (minimumBy)+import Data.Monoid (mempty, mappend) import Math.ExpPairs.LinearForm import Math.ExpPairs.Process@@ -33,10 +49,6 @@ m = lcm dq dr k = numerator q * (m `div` dq) l = numerator r * (m `div` dr)--proj2fracs :: (Integer, Integer, Integer) -> (Rational, Rational)-proj2fracs (k, l, m) = (k%m, l%m)- 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
Math/ExpPairs/Ivic.hs view
@@ -12,11 +12,20 @@ Mineola, New York: Dover Publications, 2003. -}-module Math.ExpPairs.Ivic (zetaOnS, reverseZetaOnS, mOnS) where+module Math.ExpPairs.Ivic+ ( zetaOnS+ , reverseZetaOnS+ , mOnS+ , reverseMOnS+ , checkAbscissa+ , findMinAbscissa+ , mBigOnHalf+ , reverseMBigOnHalf+ ) where -import Data.Ratio-import Data.List-import Data.Ord+import Data.Ratio ((%))+import Data.List (minimumBy)+import Data.Ord (comparing) import Math.ExpPairs @@ -61,7 +70,7 @@ -- and that alpha2 <= 1 for S >= 2/3 or S >= 5/8 and -- (4S-2)k + (8S-6)l + 2S-1 >=0 --- | Compute m(σ) such that ∫_1^T |ζ(σ+it)|^m(σ) dt ≪ T^(1+ε).+-- | Compute maximal m(σ) such that ∫_1^T |ζ(σ+it)|^m(σ) dt ≪ T^(1+ε). -- See equation (8.97) in Ivić2003. mOnS :: Rational -> OptimizeResult mOnS s@@ -97,48 +106,39 @@ x2' = optimize [RationalForm numer denom] cons x2 = x2' {optimalValue = negate $ optimalValue x2'} -reverseMOnS m = reverseMOnS' from to where+-- | 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 bigger value.+reverseMOnS :: Rational -> RationalInf -> Rational+reverseMOnS prec m = reverseMOnS' from to where from = 1 % 2 to = 1 % 1 reverseMOnS' a b- | b-a < 1%1000000 = a+ | b-a < prec = a | optimalValue (mOnS ((a+b)/2)) > m = reverseMOnS' a ((a+b)/2) | otherwise = reverseMOnS' ((a+b)/2) 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 σ. 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 -searchMinAbscissa :: [(Rational, Rational)] -> Rational-searchMinAbscissa xs = searchMinAbscissa' from to where+-- | 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+ε).+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 < 1%1000000 = a+ | b-a < prec = a | checkAbscissa xs ((a+b)/2) = searchMinAbscissa' a ((a+b)/2) | otherwise = searchMinAbscissa' ((a+b)/2) b --- % \begin{lemma}\label{l:pointwise-moments-on-1/2}--- % For $A\ge12$ let--- % $$--- % f(A) = 1+\inf \{ l/k \mid (4-A)k+4l+2 \ge 0 \}.--- % $$--- % Then--- % $$--- % R \ll T V^{-6} \log^8 T + T^{f(A)+\eps} V^{-A}--- % \ll T^{\max \{ 1+32(A-6)/205, f(A) \}} V^{-A}--- % $$--- % and thus--- % $$--- % M(A) \le \max \{ 1+32(A-6)/205, f(A) \}.--- % $$--- % \end{lemma}---- Constant--- is produced by--- optimize [RationalForm (LinearForm 4 4 2) (LinearForm 1 0 0)] [Constraint (LinearForm (-64) (-77) 64) Strict]-+-- | 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@@ -149,7 +149,14 @@ 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] +-- | 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 | otherwise = if Finite a <= optimalValue optRes@@ -158,47 +165,4 @@ a = (m-1)*205/32 + 6 optRes = optimize [RationalForm (LinearForm 4 4 2) (LinearForm 1 0 0)] [Constraint (LinearForm (1-m) 1 0) NonStrict] --f l lambda = (l-lambda)*32/205- + ((toRational . optimalValue . mBigOnHalf) (4*l/(4-lambda)) - 1) * (4-lambda) /4----bestLambda l = minimum $ map (f l) lambdas `zip` lambdas where--- lambdas = [0,1%100..4-1%100]--heckeZetaByHalf a = 1 - xt where- d = toRational 12.571624917200547- ia 2 = 0- ia 3 = 1%4- ia a- | a<=12 = 32%205 * ((1 + 4 / d) * a - 4) + (toRational (optimalValue (mBigOnHalf d)) - 1) * a / d- | a<=15 = 32%205 * a + toRational (optimalValue (mBigOnHalf a)) - 1- | otherwise = 32%205 * (2 * a - 6)- xt = 1%2 / (1 + ia a)---bestLambda l = (\x -> (x, fromRational $ f l x)) (bestLambda' 0 (3999%1000)) where- bestLambda' a b- | b-a < 1%1000000000 = a- | otherwise = if mx1 > mx2 then bestLambda' x1 b else bestLambda' a x2 where- x1 = (2*a+b)/3- x2 = (a+2*b)/3- ma = f l a- mb = f l b- mx1 = f l x1- mx2 = f l x2--difur a = a * m' - m + (77%205) where- h = 1%(10^100)- m = toRational $ optimalValue $ mBigOnHalf a- mh = toRational $ optimalValue $ mBigOnHalf (a+h)- m' = (mh-m) / h--solveD a b- | b-a < 1%(10^6) = a- | otherwise = if difur c > 0 then solveD a c else solveD c b where- c = (a+b)/2--{--D = 12.571624917200547--}
Math/ExpPairs/Kratzel.hs view
@@ -24,11 +24,16 @@ He also provided a set of theorems to estimate Θ(a, b) and Θ(a, b, c). -}-module Math.ExpPairs.Kratzel (TauabTheorem (..), tauab, TauabcTheorem (..), tauabc) where+module Math.ExpPairs.Kratzel+ ( TauabTheorem (..)+ , tauab+ , TauabcTheorem (..)+ , tauabc+ ) where -import Data.Ratio-import Data.Ord-import Data.List+import Data.Ratio ((%))+import Data.Ord (comparing)+import Data.List (minimumBy) import Math.ExpPairs @@ -117,7 +122,7 @@ [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))+ 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
Math/ExpPairs/LinearForm.hs view
@@ -1,3 +1,4 @@+{-# LANGUAGE DeriveFunctor, DeriveFoldable, DeriveGeneric #-} {-| Module : Math.ExpPairs.LinearForm Description : Linear forms, rational forms and constraints@@ -9,19 +10,35 @@ 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+module Math.ExpPairs.LinearForm+ ( LinearForm (..)+ , evalLF+ , substituteLF+ , RationalForm (..)+ , evalRF+ , IneqType (..)+ , Constraint (..)+ , checkConstraint+ ) where -import Data.List-import Data.Ratio-import Data.Monoid+import Control.DeepSeq+import Data.Foldable (Foldable (..), toList)+import Data.List (intercalate)+import Data.Ratio (numerator, denominator)+import Data.Monoid (Monoid, mempty, mappend)+import GHC.Generics (Generic (..))+ 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)+ deriving (Eq, Functor, Foldable, Generic) +instance NFData t => NFData (LinearForm t) where+ 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"@@ -34,23 +51,23 @@ instance Num t => Num (LinearForm t) where (LinearForm a b c) + (LinearForm d e f) = LinearForm (a+d) (b+e) (c+f)- (*) = undefined- negate (LinearForm a b c) = LinearForm (negate a) (negate b) (negate c)- abs = undefined- signum = undefined+ (*) = 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+ mempty = 0 mappend = (+) scaleLF :: (Num t, Eq t) => t -> LinearForm t -> LinearForm t-scaleLF 0 (LinearForm {}) = LinearForm 0 0 0-scaleLF s (LinearForm a b c) = LinearForm (a*s) (b*s) (c*s)+scaleLF 0 = const 0+scaleLF s = fmap (* s) -- |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+evalLF (k, l, m) (LinearForm a b c) = a * k + l * b + m * c -- |Substitute linear forms k, l and m into a given linear form -- a*k + b*l + c*m to obtain a new linear form.@@ -59,32 +76,33 @@ -- | Define a rational form of two variables, equal to the ratio of two 'LinearForm'. data RationalForm t = RationalForm (LinearForm t) (LinearForm t)- deriving (Show)+ deriving (Eq, Show, Functor, Foldable, Generic) +instance NFData t => NFData (RationalForm t) where+ rnf = rnf . toList+ instance Num t => Num (RationalForm t) where- (+) = undefined- (*) = undefined+ (+) = error "Addition of RationalForm is undefined"+ (*) = error "Multiplication of RationalForm is undefined" negate (RationalForm a b) = RationalForm (negate a) b- abs = undefined- signum = undefined+ 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 +mapTriple :: (a -> b) -> (a, a, a) -> (b, b, b)+mapTriple f (x, y, z) = (f x, f y, f z)+ -- |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- k = fromInteger k'- l = fromInteger l'- m = fromInteger m'- numer = toRational $ evalLF (k, l, m) num- denom = toRational $ evalLF (k, l, m) den--substituteRF :: (Eq t, Num t) => (LinearForm t, LinearForm t, LinearForm t) -> RationalForm t -> RationalForm t-substituteRF (k, l, m) (RationalForm num den) = RationalForm (substituteLF (k, l, m) num) (substituteLF (k, l, m) den)+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 -- |Constants to specify the strictness of 'Constraint'. data IneqType@@ -92,20 +110,21 @@ = Strict -- | Non-strict inequality (≥0). | NonStrict- deriving (Eq, Show)+ deriving (Eq, Ord, Show, Enum, Bounded) -- |A linear constraint of two variables. data Constraint t = Constraint (LinearForm t) IneqType- deriving (Show)+ deriving (Eq, Show, Functor, Foldable, Generic) +instance NFData t => NFData (Constraint t) where+ 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)+checkConstraint (k, l, m) (Constraint lf ineq) = if ineq==NonStrict then signum numer /= -1 else signum numer == 1 where- k = fromInteger k'- l = fromInteger l'- m = fromInteger m'- numer = evalLF (k, l, m) lf+ klm = mapTriple fromInteger (k, l, m)+ numer = evalLF klm lf
Math/ExpPairs/Matrix3.hs view
@@ -19,12 +19,15 @@ , det , multCol , normalize- , prettyMatrix+ , makarovMult+ , ladermanMult ) where import Prelude hiding (foldl1)-import Data.Foldable (Foldable (..), toList)-import GHC.Generics (Generic (..))+import Data.Foldable (Foldable (..), toList)+import GHC.Generics (Generic (..))+import Data.List (transpose)+import Control.DeepSeq -- |Three-component vector. data Vector3 t = Vector3 {@@ -34,6 +37,9 @@ } deriving (Eq, Show, Functor, Foldable, Generic) +instance NFData t => NFData (Vector3 t) where+ rnf = rnf . toList+ -- |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.,@@ -52,8 +58,11 @@ a32 :: !t, a33 :: !t }- deriving (Eq, Show, Functor, Foldable, Generic)+ deriving (Eq, Functor, Foldable, Generic) +instance NFData t => NFData (Matrix3 t) where+ rnf = rnf . toList+ diag :: Num t => t -> Matrix3 t diag n = Matrix3 { a11 = n,@@ -67,7 +76,7 @@ a33 = n } -instance Num t => Num (Matrix3 t) where+instance (Num t, Ord t) => Num (Matrix3 t) where a + b = Matrix3 { a11 = a11 a + a11 b, a12 = a12 a + a12 b,@@ -80,8 +89,18 @@ a33 = a33 a + a33 b } - -- intercalate ",\n" [ "a"++(show i)++(show j)++" = "++( intercalate " + " ["a"++(show i)++(show k)++" a * "++"a"++(show k)++(show j)++" b" | k<-[1..3]] ) | i<-[1..3], j<-[1..3]]- a * b = Matrix3 {+ (*) = usualMult++ negate = fmap negate++ abs = undefined++ signum = diag . signum . det++ 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,@@ -92,23 +111,153 @@ 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 #-} - negate = fmap negate+-- | 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.+--+-- 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.+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 - abs = undefined+ u33 = b21 - b11 - b23 - b31+ u34 = b22 - b12+ u35 = b22 - b32+ u36 = b33 - b31+ u37 = b13 - b23 - signum = diag . signum . det+ 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 - fromInteger = diag . fromInteger+ v33 = m12 + m14+ v34 = m16 + m14+ v35 = m4 + m6+ v36 = m7 + m6+ v37 = v33 + m15+ v38 = v35 + m5+ v39 = v34 + m18+ v40 = v36 + m9 --- |Computes the determinant of a matrix.-det :: Num t => Matrix3 t -> t+ 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.+-- 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.+--+-- 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.+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++ 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++ 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+{-# 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) -instance Fractional t => Fractional (Matrix3 t) where+instance (Fractional t, Ord t) => Fractional (Matrix3 t) where fromRational = diag . fromRational recip a@(Matrix3 {..}) = Matrix3 {@@ -123,7 +272,7 @@ a33 = (a11 * a22 - a12 * a21) / d } where d = det a --- |Convert a list of 9 elements into 'Matrix3'. Reverse conversion can be done using 'Foldable' instance.+-- |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,@@ -146,19 +295,12 @@ 0 -> a d -> fmap (`div` d) a --- |Print a matrix, separating rows with new lines and elements--- with spaces.-prettyMatrix :: Show t => Matrix3 t -> String-prettyMatrix Matrix3 {..} =- show a11 ++ ' ' :- show a12 ++ ' ' :- show a13 ++ '\n' :- show a21 ++ ' ' :- show a22 ++ ' ' :- show a23 ++ '\n' :- show a31 ++ ' ' :- show a32 ++ ' ' :- show a33+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 -- |Multiplicate a matrix by a vector (considered as a column). multCol :: Num t => Matrix3 t -> Vector3 t -> Vector3 t
Math/ExpPairs/MenzerNowak.hs view
@@ -17,9 +17,11 @@ Σ_{n ≤ x} τ_{a, b}(l_1, k_1; l_2, k_2; n) with an error term of order (x \/ k_1^a \/ k_2^b)^(Θ(a, b) + ε). They provided an expression for Θ(a, b) in terms of exponent pairs. -}-module Math.ExpPairs.MenzerNowak (menzerNowak) where+module Math.ExpPairs.MenzerNowak+ ( menzerNowak+ ) where -import Data.Ratio+import Data.Ratio ((%)) import Math.ExpPairs
Math/ExpPairs/Pair.hs view
@@ -12,9 +12,15 @@ Below /A/ and /B/ stands for van der Corput's processes. See "Math.ExpPairs.Process" for explanations. -}-module Math.ExpPairs.Pair (Triangle (..), InitPair' (..), InitPair, initPairs, initPairToValue) where+module Math.ExpPairs.Pair+ ( Triangle (..)+ , InitPair' (..)+ , InitPair+ , initPairs+ , initPairToValue+ ) where -import Data.Ratio+import Data.Ratio ((%)) -- |Vertices of the triangle of initial exponent pairs. data Triangle
+ Math/ExpPairs/PrettyProcess.hs view
@@ -0,0 +1,125 @@+{-|+Module : Math.ExpPairs.PrettyProcess+Description : Compact representation of process sequences+Copyright : (c) Andrew Lelechenko, 2015+License : GPL-3+Maintainer : andrew.lelechenko@gmail.com+Stability : experimental+Portability : TemplateHaskell++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 #-}+{-# LANGUAGE TemplateHaskell #-}+{-# OPTIONS_GHC -fno-warn-type-defaults #-}+module Math.ExpPairs.PrettyProcess+ ( prettify,+ uglify,+ PrettyProcess) where++import Data.List (minimumBy)+import Data.Ord (comparing)+import Data.Function.Memoize (memoize, deriveMemoizable)+import Text.PrettyPrint.Leijen++import Math.ExpPairs.ProcessMatrix++-- | Compact representation of the sequence of 'Process'.+data PrettyProcess+ = 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++-- | Width of the bracket.+bracketWidth :: Int+bracketWidth = 4++-- | Width of the subscript-sized character (e. g., power).+subscriptWidth :: Int+subscriptWidth = 4++-- | Width of the processes in typeset+processWidth :: Process -> Int+processWidth A = 10+processWidth BA = 20++-- | 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++-- | Convert 'PrettyProcess' to 'PrettyProcessWithWidth'.+annotateWithWidth :: PrettyProcess -> PrettyProcessWithWidth+annotateWithWidth p = PPWL p (printedWidth p)++-- | 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++-- | 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++-- | Find the most compact representation of the sequence of processes.+prettify :: [Process] -> PrettyProcess+prettify = ppwlProcess . prettifyP++-- | Find the most compact representation of the sequence of processes, keeping track of widthess.+prettifyP :: [Process] -> PrettyProcessWithWidth+prettifyP = memoize prettify' where++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)++ yss = xs'' : map f bcs++ bcs = takeWhile (not . null . snd) $ iterate bcf ([head xs], tail xs)++ 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++-- | 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
Math/ExpPairs/Process.hs view
@@ -9,31 +9,23 @@ 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 #-}-module Math.ExpPairs.Process (Process (..), Path (), aPath, baPath, evalPath, lengthPath) where+{-# LANGUAGE DeriveGeneric #-}+module Math.ExpPairs.Process+ ( Process ()+ , Path (Path)+ , aPath+ , baPath+ , evalPath+ , lengthPath+ ) where -import Data.Monoid-import Data.List-import Data.Ord-import Data.Function.Memoize+import GHC.Generics (Generic)+import Generics.Deriving.Monoid (Monoid, mempty, memptydefault, mappend, mappenddefault)+import Text.PrettyPrint.Leijen -import qualified Math.ExpPairs.Matrix3 as Mx --- | 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)--deriveMemoizable ''Process--type ProcessMatrix = Mx.Matrix3 Integer--process2matrix :: Process -> ProcessMatrix-process2matrix A = Mx.Matrix3 1 0 0 1 1 1 2 0 2-process2matrix BA = Mx.Matrix3 0 1 0 2 0 1 2 0 2+import Math.ExpPairs.ProcessMatrix+import Math.ExpPairs.PrettyProcess -- | Holds a list of 'Process' and a matrix of projective -- transformation, which they define. It also provides a fancy 'Show'@@ -41,22 +33,15 @@ -- -- > show (mconcat $ replicate 10 aPath) == "A^10" ---data Path = Path ProcessMatrix [Process]---- | Path consisting of a single process 'A'.-aPath :: Path-aPath = Path (process2matrix A) [ A]---- | Path consisting of a single process 'BA'.-baPath :: Path-baPath = Path (process2matrix BA) [BA]+data Path = Path !ProcessMatrix ![Process]+ deriving (Eq, Generic) instance Monoid Path where- mempty = Path 1 []- mappend (Path m1 l1) (Path m2 l2) = Path (Mx.normalize $ m1*m2) (l1++l2)+ mempty = memptydefault+ mappend = mappenddefault instance Show Path where- show (Path m l) = prettyProcesses l -- ++ "\n" ++ Mx.prettyMatrix m+ show (Path _ l) = show (pretty (prettify l)) -- ++ "\n" ++ Mx.prettyMatrix m instance Read Path where readsPrec _ zs = [reads' zs] where@@ -67,80 +52,29 @@ reads' ('B':xs) = (baPath, xs) reads' xs = (mempty, xs) -instance Eq Path where- (Path m1 _) == (Path m2 _) = Mx.normalize m1 == Mx.normalize m2- instance Ord Path where (Path _ q1) <= (Path _ q2) = cmp q1 q2 where- cmp (A:p1) (A:p2) = cmp p1 p2+ 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+ cmp (A:_) (BA:_) = True+ cmp (BA:_) (A:_) = False+ cmp [] _ = True+ cmp _ [] = False +-- | Path consisting of a single process 'A'.+aPath :: Path+aPath = Path aMatrix [ A]++-- | Path consisting of a single process 'BA'.+baPath :: Path+baPath = Path baMatrix [BA]+ -- |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 _) (a,b,c) = (a',b',c') where- m' = fmap fromInteger m- (Mx.Vector3 a' b' c') = Mx.multCol m' (Mx.Vector3 a b c)+evalPath (Path m _) = evalMatrix m -- | Count processes in the 'Path'. Note that 'BA' counts -- for one process, not two. lengthPath :: Path -> Int lengthPath (Path _ xs) = length xs--symbolWidth :: Int-symbolWidth = 10-bracketWidth :: Int-bracketWidth = 4-subscriptWidth :: Int-subscriptWidth = 4---- Пусть строка из n символов является повторением строки из l символов--- Какова длина такой записи?-len0 :: [Process] -> Int -> (Int, String)-len0 xs 1 = (lxs, pxs) where- lxs = length pxs * symbolWidth- pxs = concatMap show xs-len0 [A] n = (symbolWidth + subscriptWidth, show A ++ "^" ++ show n)-len0 [BA] n = (symbolWidth + bracketWidth*2 + subscriptWidth, "(" ++ show BA ++ ")^" ++ show n)-len0 xs n = (lxs + bracketWidth*2 + subscriptWidth, "(" ++ pxs ++ ")^" ++ show n) where- (lxs, pxs) = len2M xs--len0M :: [Process] -> Int -> (Int, String)-len0M = memoize len0---- Простейшая оптимизация: строка as, целиком состоящая из n повторений подстроки bs, может быть записана как bs^n-len1 :: [Process] -> (Int, String)-len1 as = if null inner- then len0M as 1- else len0M as 1 `min` minimumBy (comparing fst) inner where- l = length as- bs n = take n as- cs m xs = concat (replicate m xs)- inner = [len0M (bs n) (l`div`n) | n<-[1..l-1], l`mod`n==0, cs (l`div`n) (bs n) == as]--len1M :: [Process] -> (Int, String)-len1M = memoize len1---- Перебираем все способы разбить строку на две части и применить к каждой из них len1-len2 :: [Process] -> (Int, String)-len2 as = if null inner- then len1M as- else len1M as `min` minimumBy (comparing fst) inner where- l = length as- bs n = take n as- cs n = drop n as- add (x, xs) (y, ys) = (x+y, xs++ys)- inner = [ len2M (bs n) `add` len2M (cs n) | n<-[1..l-1] ]--len2M :: [Process] -> (Int, String)-len2M = memoize len2--prettyProcesses :: [Process] -> String-prettyProcesses = snd . len2M---
+ Math/ExpPairs/ProcessMatrix.hs view
@@ -0,0 +1,61 @@+{-|+Module : Math.ExpPairs.ProcessMatrix+Description : Monoidal wrapper for Matrix3+Copyright : (c) Andrew Lelechenko, 2014-2015+License : GPL-3+Maintainer : andrew.lelechenko@gmail.com+Stability : experimental+Portability : TemplateHaskell++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 #-}+module Math.ExpPairs.ProcessMatrix+ ( Process (..)+ , ProcessMatrix ()+ , aMatrix+ , baMatrix+ , evalMatrix+ ) where++import Data.Monoid (Monoid, mempty, mappend)+import Data.Function.Memoize (deriveMemoizable)++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)++deriveMemoizable ''Process++newtype ProcessMatrix = ProcessMatrix (Matrix3 Integer)+ deriving (Eq, Num, Show)++instance Monoid ProcessMatrix where+ 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+process2matrix BA = ProcessMatrix $ Matrix3 0 1 0 2 0 1 2 0 2++-- | Return process matrix for 'A'-process.+aMatrix :: ProcessMatrix+aMatrix = process2matrix A++-- | Return process matrix for 'BA'-process.+baMatrix :: ProcessMatrix+baMatrix = process2matrix BA++-- |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)+
Math/ExpPairs/RatioInf.hs view
@@ -1,3 +1,4 @@+{-# LANGUAGE BangPatterns #-} {-| Module : Math.ExpPairs.RatioInf Description : Rational numbers with infinities@@ -9,9 +10,12 @@ 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+module Math.ExpPairs.RatioInf+ ( RatioInf (..)+ , RationalInf+ ) where -import Data.Ratio+import Data.Ratio (Ratio) -- |Extends a rational type with positive and negative -- infinities.@@ -19,7 +23,7 @@ -- |Negative infinity = InfMinus -- |Finite value- | Finite (Ratio t)+ | Finite !(Ratio t) -- |Positive infinity | InfPlus deriving (Ord, Eq)@@ -33,7 +37,7 @@ show (Finite x) = show x show InfPlus = "+Inf" -instance (Integral t) => Num (RatioInf t) where+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@@ -42,7 +46,7 @@ _ + InfPlus = InfPlus (Finite a) + (Finite b) = Finite (a+b) - fromInteger n = Finite (fromInteger n)+ fromInteger = Finite . fromInteger signum InfMinus = Finite (-1) signum InfPlus = Finite 1@@ -56,48 +60,56 @@ negate InfPlus = InfMinus negate (Finite r) = Finite (negate r) - InfMinus * a- | signum a == Finite 0 = error "Cannot multiply infinity by zero"- | signum a == Finite 1 = InfMinus- | signum a == Finite (-1) = InfPlus- InfPlus * a- | signum a == Finite 0 = error "Cannot multiply infinity by zero"- | signum a == Finite 1 = InfPlus- | signum a == Finite (-1) = InfMinus- a * InfMinus- | signum a == Finite 0 = error "Cannot multiply infinity by zero"- | signum a == Finite 1 = InfMinus- | signum a == Finite (-1) = InfPlus- a * InfPlus- | signum a == Finite 0 = error "Cannot multiply infinity by zero"- | signum a == Finite 1 = InfPlus- | signum a == Finite (-1) = InfMinus- (Finite a) * (Finite b) = Finite (a*b)+ InfMinus * InfMinus = InfMinus+ InfMinus * InfPlus = InfMinus+ InfMinus * Finite a = case signum a of+ 1 -> InfMinus+ -1 -> InfPlus+ _ -> error "Cannot multiply infinity by zero" -instance (Integral t) => Fractional (RatioInf t) where+ 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 * InfPlus = case signum a of+ 1 -> InfPlus+ -1 -> InfMinus+ _ -> error "Cannot multiply infinity by zero"++ Finite a * Finite b = Finite (a * b)++instance Integral t => Fractional (RatioInf t) where fromRational = Finite . fromRational InfMinus / InfMinus = error "Cannot divide infinity by infinity" InfMinus / InfPlus = error "Cannot divide infinity by infinity"- InfMinus / (Finite a)- | signum a == 0 = error "Cannot divide infinity by zero"- | signum a == 1 = InfMinus- | signum a == -1 = InfPlus+ 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)- | signum a == 0 = error "Cannot divide infinity by zero"- | signum a == 1 = InfPlus- | signum a == -1 = InfMinus+ 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) -instance (Integral t) => Real (RatioInf t) where+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"
exp-pairs.cabal view
@@ -1,5 +1,5 @@ name: exp-pairs-version: 0.1.2.0+version: 0.1.3.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@@ -24,10 +24,15 @@ Math.ExpPairs.Matrix3, Math.ExpPairs.Pair, Math.ExpPairs.Process,+ Math.ExpPairs.PrettyProcess,+ Math.ExpPairs.ProcessMatrix, Math.ExpPairs.RatioInf build-depends: base >=4 && <5, memoize >=0.1,- ghc-prim+ ghc-prim,+ generic-deriving,+ wl-pprint >=1.2,+ deepseq >=1.3 default-language: Haskell2010 ghc-options: -Wall
tests/Tests.hs view
@@ -2,6 +2,7 @@ import qualified Matrix3 (testSuite) import qualified RatioInf (testSuite) import qualified Pair (testSuite)+import qualified PrettyProcess (testSuite) import qualified Ivic (testSuite) import qualified Kratzel (testSuite)@@ -14,10 +15,11 @@ tests :: TestTree tests = testGroup "Tests"- [ LinearForm.testSuite- , Matrix3.testSuite+ [ Matrix3.testSuite+ , LinearForm.testSuite , RatioInf.testSuite , Pair.testSuite+ , PrettyProcess.testSuite , Ivic.testSuite , Kratzel.testSuite , MenzerNowak.testSuite