-- | A representation of rational numbers as lists of prime powers, -- 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 generating all divisors of -- factored integers, or computing Euler's totient (phi) function -- and the Möbius (mu) function. module MathObj.FactoredRational ( -- * Type T -- * Utilities , factorial , eulerPhi , divisors , mu ) where import qualified Algebra.Additive as Additive import qualified Algebra.Ring as Ring import qualified Algebra.Field as Field import qualified Algebra.IntegralDomain as Integral import qualified Algebra.ZeroTestable as ZeroTestable import qualified Algebra.ToRational as ToRational import qualified Algebra.RealIntegral as RealIntegral import qualified Algebra.ToInteger as ToInteger import Data.Numbers.Primes import qualified Algebra.Absolute as Real import NumericPrelude -- Represent rational numbers by their prime factorizations. -- 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. -- -- Instances are provided for Eq, Ord, Additive, Ring, ZeroTestable, -- Real, ToRational, Integral, RealIntegral, ToInteger, and Field. -- -- Note that currently, addition is performed on factored rationals -- by converting them to normal rationals, performing the addition, -- and factoring. This could probably be made more efficient by -- finding a common denominator, pulling out common factors from the -- numerators, and performing the addition and factoring only on the -- relatively prime parts. data T = FQZero -- ^ zero | FQ Bool [Integer] -- ^ prime exponents with sign bit, True = negative. -- XXX this ought to be improved. instance Show T where show FQZero = "0" show (FQ True pows) = "(-1)" ++ showPows pows show (FQ False []) = "1" show (FQ False pows) = showPows pows showPows :: [Integer] -> String showPows pows = concat $ zipWith showPow primes pows where showPow :: Integer -> Integer -> String showPow _ 0 = "" showPow p 1 = "(" ++ show p ++ ")" showPow p n = "(" ++ show p ++ "^" ++ show n ++ ")" instance Additive.C T where zero = FQZero FQZero + a = a a + FQZero = a x + y = fromRational' (toRational x + toRational y) negate FQZero = FQZero negate (FQ s e) = FQ (not s) e instance Ring.C T where FQZero * _ = FQZero _ * FQZero = FQZero (FQ s1 e1) * (FQ s2 e2) = FQ (s1 /= s2) (zipWithExt 0 0 (+) e1 e2) fromInteger 0 = FQZero fromInteger n | n < 0 = FQ True (factor (negate n)) | otherwise = FQ False (factor n) _ ^ 0 = one FQZero ^ _ = FQZero (FQ s e) ^ n = FQ s (map (*n) e) -- | Zip two lists together with a combining function, using default -- values to extend the lists if one is shorter than the other. zipWithExt :: a -> b -> (a -> b -> c) -> [a] -> [b] -> [c] zipWithExt da db f = zipWithExt' where zipWithExt' [] bs = zipWith f (repeat da) bs zipWithExt' as [] = zipWith f as (repeat db) zipWithExt' (a:as) (b:bs) = f a b : zipWithExt' as bs -- | A simple factoring method. -- -- We should probably just depend on another module with some -- dedicated, efficient factoring code written by someone really -- smart, but this simple method works OK for now. -- -- Precondition: argument is positive. factor :: Integer -> [Integer] factor m = factor' m primes where factor' 1 _ = [] factor' n (p:ps) = let (k,n') = logRem n p in k : factor' n' ps factor' _ [] = error "Oh noes, Euclid was wrong!" -- | @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) logRem = logRem' 0 where logRem' k n p | n `mod` p == 0 = logRem' (k+1) (n `div` p) p | otherwise = (k,n) instance ZeroTestable.C T where isZero FQZero = True isZero _ = False instance Eq T where FQZero == FQZero = True FQZero == (FQ _ _) = False (FQ _ _) == FQZero = False (FQ s1 e1) == (FQ s2 e2) = s1 == s2 && dropZeros e1 == dropZeros e2 where dropZeros = reverse . dropWhile (==0) . reverse instance Ord T where compare FQZero FQZero = EQ compare FQZero (FQ False _) = LT compare FQZero (FQ True _) = GT compare (FQ False _) FQZero = GT compare (FQ True _) FQZero = LT compare (FQ False _) (FQ True _) = GT compare (FQ True _) (FQ False _) = LT compare fq1 fq2 = compare (toRational fq1) (toRational fq2) instance Real.C T where abs FQZero = FQZero abs (FQ _ e) = FQ False e signum = Real.signumOrd instance ToRational.C T where toRational FQZero = 0 toRational (FQ s e) = (if s then negate else id) . product . zipWith (^-) (map (%1) primes) $ e instance Integral.C T where divMod a b = if isZero b then (undefined,a) else (a/b,zero) instance RealIntegral.C T -- default definition is fine instance ToInteger.C T where toInteger FQZero = 0 toInteger (FQ s e) | any (<0) e = error "non-integer in FactoredRational.toInteger" | otherwise = (if s then negate else id) . product . zipWith (^) primes $ e instance Field.C T where recip FQZero = error "division by zero" recip (FQ s e) = FQ s (map negate e) -- | Efficiently compute n! directly as a factored rational. factorial :: Integer -> T factorial 0 = one factorial 1 = one factorial n = FQ False (takeWhile (>0) . map (factorialFactors n) $ primes) -- | @factorialFactors n p@ computes the power of prime p in the -- factorization of n!. factorialFactors :: Integer -> Integer -> Integer factorialFactors n p = sum . takeWhile (>0) . map (n `div`) $ iterate (*p) p -- | Compute Euler's totient function (@eulerPhi n@ is the number of -- integers less than and relatively prime to n). Only makes sense -- for (nonnegative) integers. eulerPhi :: T -> Integer eulerPhi FQZero = 1 eulerPhi (FQ _ pows) = product $ zipWith phiP primes pows where phiP _ 0 = 1 phiP p a = p^(a-1) * (p-1) -- | List of the divisors of n. Only makes sense for integers. divisors :: T -> [T] divisors FQZero = [one] 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