packages feed

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 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