diff --git a/MathObj/FactoredRational.hs b/MathObj/FactoredRational.hs --- a/MathObj/FactoredRational.hs +++ b/MathObj/FactoredRational.hs @@ -2,11 +2,12 @@ -- allowing efficient representation, multiplication and division of -- large numbers, especially of the sort occurring in combinatorial -- computations. --- +-- -- The module also includes a method for generating factorials in --- factored form directly, and for computing Euler's totient and --- generating all divisors of factored integers. -module MathObj.FactoredRational +-- factored form directly, and for generating all divisors of +-- factored integers, or computing Euler's totient (phi) function +-- and the Möbius (mu) function. +module MathObj.FactoredRational ( -- * Type T @@ -14,7 +15,8 @@ , factorial , eulerPhi , divisors - + , mu + ) where import qualified Algebra.Additive as Additive @@ -30,11 +32,11 @@ import Data.Numbers.Primes -import PreludeBase +import PreludeBase import NumericPrelude -- Represent rational numbers by their prime factorizations. --- Perhaps this should use a sparse representation instead, using a Map from +-- Perhaps this should use a sparse representation instead, using a Map from -- primes to powers? Well, that should be easy enough to change later. -- | The type of factored rationals. @@ -54,7 +56,8 @@ -- XXX this ought to be improved. instance Show T where show FQZero = "0" - show (FQ True pows) = "(-1)" ++ showPows pows + show (FQ True pows) = "(-1)" ++ showPows pows + show (FQ False []) = "1" show (FQ False pows) = showPows pows showPows :: [Integer] -> String @@ -67,7 +70,7 @@ zero = FQZero FQZero + a = a a + FQZero = a - x + y = fromRational' (toRational x + toRational y) + x + y = fromRational' (toRational x + toRational y) negate FQZero = FQZero negate (FQ s e) = FQ (not s) e @@ -93,7 +96,7 @@ zipWithExt' as [] = zipWith f as (repeat db) zipWithExt' (a:as) (b:bs) = f a b : zipWithExt' as bs --- | A simple factoring method. +-- | A simple factoring method. -- -- We should probably just depend on another module with some -- dedicated, efficient factoring code written by someone really @@ -106,7 +109,7 @@ factor' 1 _ = [] factor' n (p:ps) = let (k,n') = logRem n p in k : factor' n' ps - + -- | @logRem n p@ computes (k,n'), where k is the highest power of p -- that divides n, and n' = n `div` p^k. logRem :: Integer -> Integer -> (Integer, Integer) @@ -114,7 +117,7 @@ where logRem' k n p | n `mod` p == 0 = logRem' (k+1) (n `div` p) p | otherwise = (k,n) -instance ZeroTestable.C T where +instance ZeroTestable.C T where isZero FQZero = True isZero _ = False @@ -152,7 +155,7 @@ then (undefined,a) else (a/b,0) -instance RealIntegral.C T +instance RealIntegral.C T -- default definition is fine instance ToInteger.C T where @@ -176,7 +179,7 @@ -- | @factorialFactors n p@ computes the power of prime p in the -- factorization of n!. factorialFactors :: Integer -> Integer -> Integer -factorialFactors n p = sum +factorialFactors n p = sum . takeWhile (>0) . map (n `div`) $ iterate (*p) p @@ -194,3 +197,13 @@ divisors :: T -> [T] divisors FQZero = [1] divisors (FQ b pows) = map (FQ b) $ mapM (enumFromTo 0) pows + +-- | Möbius (mu) function of a positive integer: mu(n) = 0 if one or +-- more repeated prime factors, 1 if n = 1, and (-1)^k if n is a +-- product of k distinct primes. +mu :: T -> Integer +mu FQZero = error "FactoredRational.mu: zero argument" +mu (FQ True _) = error "FactoredRational.mu: negative argument" +mu (FQ _ []) = 1 -- mu(1) = 1 +mu (FQ _ pows) | all (`elem` [0,1]) pows = (-1)^(sum pows) + | otherwise = 0 diff --git a/MathObj/Monomial.hs b/MathObj/Monomial.hs --- a/MathObj/Monomial.hs +++ b/MathObj/Monomial.hs @@ -6,6 +6,7 @@ T(..) -- * Creating monomials + , mkMonomial , constant , x @@ -51,11 +52,15 @@ , powers :: M.Map Integer Integer } +mkMonomial :: a -> [(Integer, Integer)] -> T a +mkMonomial a p = Cons a (M.fromList p) + instance (ZeroTestable.C a, Ring.C a, Eq a, Show a) => Show (T a) where - show (Cons a pows) | isZero a = "0" - | a == 1 = showVars pows - | a == (-1) = "-" ++ showVars pows - | otherwise = show a ++ " " ++ showVars pows + show (Cons a pows) | isZero a = "0" + | M.null pows = show a + | a == 1 = showVars pows + | a == (-1) = "-" ++ showVars pows + | otherwise = show a ++ " " ++ showVars pows showVars :: M.Map Integer Integer -> String showVars m = intercalate " " $ concatMap showVar (M.toList m) diff --git a/MathObj/MultiVarPolynomial.hs b/MathObj/MultiVarPolynomial.hs --- a/MathObj/MultiVarPolynomial.hs +++ b/MathObj/MultiVarPolynomial.hs @@ -15,6 +15,8 @@ , compose + , merge + ) where import qualified Algebra.Additive as Additive @@ -73,16 +75,20 @@ -- | Add two polynomials. We assume that they are already sorted, so -- that addition works on infinite polynomials. add :: (Ord a, Additive.C a) => [a] -> [a] -> [a] -add xs ys = merge (+) xs ys +add xs ys = merge True (+) xs ys --- | Merge two sorted lists, with a combining function for elements --- that are equal. -merge :: Ord a => (a -> a -> a) -> [a] -> [a] -> [a] -merge _ [] ys = ys -merge _ xs [] = xs -merge f xxs@(x:xs) yys@(y:ys) | x < y = x : merge f xs yys - | x > y = y : merge f xxs ys - | otherwise = (f x y) : merge f xs ys +-- | Merge two sorted lists, with a flag specifying whether to keep +-- singletons, and a combining function for elements that are equal. +merge :: Ord a => Bool -> (a -> a -> a) -> [a] -> [a] -> [a] +merge True _ [] ys = ys +merge False _ [] _ = [] +merge True _ xs [] = xs +merge False _ _ [] = [] +merge b f xxs@(x:xs) yys@(y:ys) | x < y = if' b (x:) id $ merge b f xs yys + | x > y = if' b (y:) id $ merge b f xxs ys + | otherwise = f x y : merge b f xs ys + where if' True x _ = x + if' False _ y = y instance (Additive.C a, ZeroTestable.C a) => Additive.C (T a) where zero = fromMonomials [] diff --git a/np-extras.cabal b/np-extras.cabal --- a/np-extras.cabal +++ b/np-extras.cabal @@ -1,5 +1,5 @@ name: np-extras -version: 0.1 +version: 0.2 license: BSD3 license-file: LICENSE build-type: Simple