PrimeFactors.hs

-----------------------------------------------------------------------------
-- Copyright 2015, Open Universiteit Nederland. This file is distributed
-- under the terms of the GNU General Public License. For more information,
-- see the file "LICENSE.txt", which is included in the distribution.
-----------------------------------------------------------------------------
-- |
-- Maintainer  :  bastiaan.heeren@ou.nl
-- Stability   :  provisional
-- Portability :  portable (depends on ghc)
--
-----------------------------------------------------------------------------
--  $Id: PrimeFactors.hs 7527 2015-04-08 07:58:06Z bastiaan $

module Domain.Math.Data.PrimeFactors
   ( PrimeFactors
   , splitPower, greatestPower, allPowers
   ) where

import Data.Maybe
import Domain.Math.Data.Primes
import qualified Data.IntMap as IM

-------------------------------------------------------------
-- Representation

-- Invariants:
-- * Keys in map are prime numbers only (exception: representation of 0)
-- * Elements in map are positive (non-zero)
-- * Zero is represented by [(0,1)] (since 0^1 equals 0)
-- * The number can be negative, in which case we use the factors of
--   its absolute value
data PrimeFactors = PF Integer Factors

type Factors = IM.IntMap Int

-------------------------------------------------------------
-- Conversion to and from factors

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

-------------------------------------------------------------
-- Type class instances

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 -- prevent recomputing prime factors
      | 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

-------------------------------------------------------------
-- Utility functions

-- brute force, ugly
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)

-- -- n == a^x with (a,x) == greatestPower n
-- prop_greatestPower n = traceShow n $
--    maybe True (\(a,x) -> fromIntegral a ^ fromIntegral x == n) $ greatestPower n

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]

-- prop_allPowers n = traceShow n $
--   and (map (\(a,x) -> fromIntegral a ^ fromIntegral x == n) (allPowers n))

-- splitPower i a = (b,c)
--  => b^i * c = a
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
                 )