module Euler.Util.Prime
(
factorizations
, factorizationsAsc
, primes
, isPrime
, primeFactors
, smallestPrimeFactor
, largestPrimeFactor
, divisors
, countDivisors
, divisorsOfPrimeProduct
, properDivisors
, properDivisorsOfPrimeProduct
) where
import Euler.Prelude
import Euler.Util.Arithmetic (divides)
import Euler.Util.List (untilNothing)
import qualified Euler.Util.FrontierSearch as FS
import qualified Euler.Util.Inf as Inf
import Data.List (length, product, takeWhile)
import qualified Data.Foldable as Foldable
import qualified Data.List as List
import qualified Data.Map.Strict as Map'
import qualified Data.MultiSet as MultiSet
import qualified Math.Combinatorics.Multiset as MS
factorizations :: Integral a => a -> Map a [a]
factorizations n = toMap $ [] : f []
where
f tail = (fold . untilNothing . fmap g) primes
where
g p = let fs = p : tail
in if product fs > n then Nothing else Just $ fs : f fs
toMap = Map'.fromList . fmap (liftA2 (,) Foldable.product sort)
factorizationsAsc :: Integral a => [[a]]
factorizationsAsc = (\(_, fs, _) -> MultiSet.toAscList fs) <$> FS.searchNodes FS.Conf
{ FS.start = [(1, MultiSet.empty, Inf.fromList primes)]
, FS.next = \(n, factors, higherPrimes) ->
let h:hs = Inf.toList higherPrimes in
(if MultiSet.size factors <= 1
then [(h, MultiSet.singleton h, Inf.fromList hs)]
else []) <>
((\f -> (n*f, MultiSet.insert f factors, higherPrimes)) <$> takeWhile (<h) primes)
, FS.nodeValue = \(n, _, _) -> n
}
primeFactors :: (Integral a, Integral b) => a -> [b]
primeFactors n =
sort $ primeFactors' [] n
where
primeFactors' :: (Integral x, Integral y) => [x] -> y -> [x]
primeFactors' f' 1 = f'
primeFactors' f' n' = primeFactors' f'' n''
where
f'' = d : f'
n'' = fromIntegral n' `div` fromIntegral d
d = smallestPrimeFactor n'
smallestPrimeFactor :: (Integral a, Integral b) => a -> b
smallestPrimeFactor n =
List.head $ List.filter (`divides` fromIntegral n) primes
largestPrimeFactor :: (Integral a, Integral b) => a -> b
largestPrimeFactor n
| isPrime n = fromIntegral n
| otherwise = largestPrimeFactor $ n `div` smallestPrimeFactor n
countDivisors :: (Integral a, Integral b) => a -> b
countDivisors =
Foldable.product . fmap (fromIntegral . (+1) . List.length) . group . primeFactors
divisors :: Integral a => a -> [a]
divisorsOfPrimeProduct :: Integral a => [a] -> [a]
properDivisors :: Integral a => a -> [a]
properDivisorsOfPrimeProduct :: Integral a => [a] -> [a]
divisorsOfPrimeProduct xs = fmap product $ maxKSubMultisets (length xs) xs
properDivisorsOfPrimeProduct xs = fmap product $ maxKSubMultisets (length xs - 1) xs
divisors = divisorsOfPrimeProduct . primeFactors
properDivisors = properDivisorsOfPrimeProduct . primeFactors
maxKSubMultisets :: Ord a => MS.Count -> [a] -> [[a]]
maxKSubMultisets maxK xs =
let set = MS.fromList xs
in do k <- [0 .. maxK]
subset <- MS.kSubsets k set
return $ MS.toList subset