module Domain.Math.Data.PrimeFactors
( PrimeFactors
, splitPower, greatestPower, allPowers
) where
import Data.Maybe
import Domain.Math.Data.Primes
import qualified Data.IntMap as IM
data PrimeFactors = PF Integer Factors
type Factors = IM.IntMap Int
toFactors :: Integer -> Factors
toFactors a
| a == 0 = IM.singleton 0 1
| otherwise = rec $ primeFactors $ abs $ fromInteger a
where
rec [] = IM.empty
rec (x:xs) = IM.insert x (length ys + 1) (rec zs)
where
(ys, zs) = break (/= x) xs
fromFactors :: Factors -> Integer
fromFactors = product . map f . IM.toList
where f (a, i) = toInteger a ^ toInteger i
instance Show PrimeFactors where
show (PF a m) = show a ++ " (factors = " ++ show (IM.toList m) ++ ")"
instance Eq PrimeFactors where
PF a _ == PF b _ = a==b
instance Ord PrimeFactors where
PF a _ `compare` PF b _ = a `compare` b
instance Num PrimeFactors where
PF a m1 + PF b m2
| a==0 = PF b m2
| b==0 = PF a m1
| otherwise = fromInteger (a+b)
PF a m1 * PF b m2
| a==0 || b==0 = 0
| otherwise = PF (a*b) (IM.unionWith (+) m1 m2)
negate (PF a m) = PF (negate a) m
abs (PF a m) = PF (abs a) m
signum (PF a _) = fromInteger (signum a)
fromInteger n = PF n (toFactors n)
instance Enum PrimeFactors where
toEnum = fromIntegral
fromEnum = fromIntegral . toInteger
instance Real PrimeFactors where
toRational = toRational . toInteger
instance Integral PrimeFactors where
toInteger (PF a _) = a
quotRem = quotRemPF
greatestPower :: Integer -> Maybe (Integer, Integer)
greatestPower n = f 2 1
where
f b e | n == b ^ e = Just (b, e)
| b > n = Nothing
| b ^ e > n = f (b + 1) 1
| otherwise = f b (e + 1)
allPowers :: Integer -> [(Integer, Integer)]
allPowers n = do
(b, e) <- maybeToList $ greatestPower n
let f i = let (a, r) = e `divMod` i
in if a > 1 && r == 0 then Just (b^i, a) else Nothing
mapMaybe f [1..e]
splitPower :: Int -> PrimeFactors -> (PrimeFactors, PrimeFactors)
splitPower i (PF a m) = (PF b p1, PF c p2)
where
pairs = IM.map (`quotRem` i) m
p1 = IM.filter (>0) (fmap fst pairs)
p2 = IM.filter (>0) (fmap snd pairs)
b = fromFactors p1
c = a `div` (b^i)
quotRemPF :: PrimeFactors -> PrimeFactors -> (PrimeFactors, PrimeFactors)
quotRemPF (PF a m1) (PF b m2)
| b==0 = error "PrimeFactors: division by zero"
| a==0 = (0,0)
| otherwise = sign $
case (IM.null up, IM.null dn) of
(True, True) -> (1, 0)
(False, True) -> (PF (fromFactors up) up, 0)
(True, False) -> (0, PF a m1)
_ -> (fromInteger qn, fromInteger rn)
where
(up, dn) = IM.partition (>0) $ IM.filter (/=0) $ IM.unionWith (+) m1 (IM.map negate m2)
(qn, rn) = fromFactors up `quotRem` fromFactors (IM.map negate dn)
sign (q, r) = ( fromInteger (signum a*signum b) * q
, fromInteger (signum a) * r
)