-- | Some auxilary functions used internally

{-# LANGUAGE CPP, BangPatterns, TypeSynonymInstances, FlexibleInstances, DeriveFunctor #-}
module Math.Algebra.Polynomial.Misc where

--------------------------------------------------------------------------------

import Data.List
import Data.Monoid
import Data.Ratio
import Data.Array

-- Semigroup became a superclass of Monoid
#if MIN_VERSION_base(4,11,0)     
import Data.Foldable
import Data.Semigroup
#endif

import Control.Monad

import qualified Data.Map.Strict as Map ; import Data.Map (Map)

--------------------------------------------------------------------------------

{-
-- * Pairs

data Pair a 
  = Pair a a 
  deriving (Eq,Ord,Show,Functor)
-}

--------------------------------------------------------------------------------

equating :: Eq b => (a -> b) -> a -> a -> Bool  
equating :: (a -> b) -> a -> a -> Bool
equating a -> b
f a
x a
y = (a -> b
f a
x b -> b -> Bool
forall a. Eq a => a -> a -> Bool
== a -> b
f a
y)

--------------------------------------------------------------------------------
-- * Lists

unique :: Ord a => [a] -> [a]
unique :: [a] -> [a]
unique = ([a] -> a) -> [[a]] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map [a] -> a
forall a. [a] -> a
head ([[a]] -> [a]) -> ([a] -> [[a]]) -> [a] -> [a]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [a] -> [[a]]
forall a. Eq a => [a] -> [[a]]
group ([a] -> [[a]]) -> ([a] -> [a]) -> [a] -> [[a]]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [a] -> [a]
forall a. Ord a => [a] -> [a]
sort

-- | Synonym for histogram
count :: Ord b => [b] -> Map b Integer
count :: [b] -> Map b Integer
count = [b] -> Map b Integer
forall b. Ord b => [b] -> Map b Integer
histogram

histogram :: Ord b => [b] -> Map b Integer
histogram :: [b] -> Map b Integer
histogram [b]
xs = (Map b Integer -> b -> Map b Integer)
-> Map b Integer -> [b] -> Map b Integer
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl' Map b Integer -> b -> Map b Integer
forall k a. (Ord k, Num a) => Map k a -> k -> Map k a
f Map b Integer
forall k a. Map k a
Map.empty [b]
xs where
  f :: Map k a -> k -> Map k a
f Map k a
old k
x = (a -> a -> a) -> k -> a -> Map k a -> Map k a
forall k a. Ord k => (a -> a -> a) -> k -> a -> Map k a -> Map k a
Map.insertWith a -> a -> a
forall a. Num a => a -> a -> a
(+) k
x a
1 Map k a
old

{-# SPECIALIZE sum' :: [Int]     -> Int     #-}
{-# SPECIALIZE sum' :: [Integer] -> Integer #-}
sum' :: Num a => [a] -> a
sum' :: [a] -> a
sum' = (a -> a -> a) -> a -> [a] -> a
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl' a -> a -> a
forall a. Num a => a -> a -> a
(+) a
0

longZipWith :: (a -> c) -> (b -> c) -> (a -> b -> c) -> [a] -> [b] -> [c]
longZipWith :: (a -> c) -> (b -> c) -> (a -> b -> c) -> [a] -> [b] -> [c]
longZipWith a -> c
f b -> c
g a -> b -> c
h = [a] -> [b] -> [c]
go where
  go :: [a] -> [b] -> [c]
go (a
x:[a]
xs) (b
y:[b]
ys) = a -> b -> c
h a
x b
y c -> [c] -> [c]
forall a. a -> [a] -> [a]
: [a] -> [b] -> [c]
go [a]
xs [b]
ys
  go [a]
xs     []     = (a -> c) -> [a] -> [c]
forall a b. (a -> b) -> [a] -> [b]
map a -> c
f [a]
xs
  go []     [b]
ys     = (b -> c) -> [b] -> [c]
forall a b. (a -> b) -> [a] -> [b]
map b -> c
g [b]
ys

longReplaceListElem :: a -> Int -> a -> [a] -> [a]
longReplaceListElem :: a -> Int -> a -> [a] -> [a]
longReplaceListElem a
x0 Int
j a
y [a]
xs = Int -> [a] -> [a]
forall a. (Eq a, Num a) => a -> [a] -> [a]
go Int
j [a]
xs  where
  go :: a -> [a] -> [a]
go  a
0 (a
x:[a]
xs) = a
y  a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [a]
xs
  go !a
i (a
x:[a]
xs) = a
x  a -> [a] -> [a]
forall a. a -> [a] -> [a]
: a -> [a] -> [a]
go (a
ia -> a -> a
forall a. Num a => a -> a -> a
-a
1) [a]
xs
  go  a
0 []     = a
y  a -> [a] -> [a]
forall a. a -> [a] -> [a]
: []
  go !a
i []     = a
x0 a -> [a] -> [a]
forall a. a -> [a] -> [a]
: a -> [a] -> [a]
go (a
ia -> a -> a
forall a. Num a => a -> a -> a
-a
1) []

--------------------------------------------------------------------------------
-- * Maps
  
deleteLookup :: Ord a => a -> Map a b -> (Maybe b, Map a b)
deleteLookup :: a -> Map a b -> (Maybe b, Map a b)
deleteLookup a
k Map a b
table = (a -> Map a b -> Maybe b
forall k a. Ord k => k -> Map k a -> Maybe a
Map.lookup a
k Map a b
table, a -> Map a b -> Map a b
forall k a. Ord k => k -> Map k a -> Map k a
Map.delete a
k Map a b
table)  

unsafeDeleteLookup :: Ord a => a -> Map a b -> (b, Map a b)
unsafeDeleteLookup :: a -> Map a b -> (b, Map a b)
unsafeDeleteLookup a
k Map a b
table = (Maybe b -> b
forall p. Maybe p -> p
fromJust (a -> Map a b -> Maybe b
forall k a. Ord k => k -> Map k a -> Maybe a
Map.lookup a
k Map a b
table), a -> Map a b -> Map a b
forall k a. Ord k => k -> Map k a -> Map k a
Map.delete a
k Map a b
table) where
  fromJust :: Maybe p -> p
fromJust Maybe p
mb = case Maybe p
mb of
    Just p
y  -> p
y
    Maybe p
Nothing -> [Char] -> p
forall a. HasCallStack => [Char] -> a
error [Char]
"unsafeDeleteLookup: key not found"

-- | Example usage: @insertMap (:[]) (:) ...@
insertMap :: Ord k => (b -> a) -> (b -> a -> a) -> k -> b -> Map k a -> Map k a
insertMap :: (b -> a) -> (b -> a -> a) -> k -> b -> Map k a -> Map k a
insertMap b -> a
f b -> a -> a
g k
k b
y = (Maybe a -> Maybe a) -> k -> Map k a -> Map k a
forall k a.
Ord k =>
(Maybe a -> Maybe a) -> k -> Map k a -> Map k a
Map.alter Maybe a -> Maybe a
h k
k where
  h :: Maybe a -> Maybe a
h Maybe a
mb = case Maybe a
mb of
    Maybe a
Nothing -> a -> Maybe a
forall a. a -> Maybe a
Just (b -> a
f b
y)
    Just a
x  -> a -> Maybe a
forall a. a -> Maybe a
Just (b -> a -> a
g b
y a
x)    

-- | Example usage: @buildMap (:[]) (:) ...@
buildMap :: Ord k => (b -> a) -> (b -> a -> a) -> [(k,b)] -> Map k a
buildMap :: (b -> a) -> (b -> a -> a) -> [(k, b)] -> Map k a
buildMap b -> a
f b -> a -> a
g [(k, b)]
xs = (Map k a -> (k, b) -> Map k a) -> Map k a -> [(k, b)] -> Map k a
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl' Map k a -> (k, b) -> Map k a
forall k. Ord k => Map k a -> (k, b) -> Map k a
worker Map k a
forall k a. Map k a
Map.empty [(k, b)]
xs where
  worker :: Map k a -> (k, b) -> Map k a
worker !Map k a
old (k
k,b
y) = (Maybe a -> Maybe a) -> k -> Map k a -> Map k a
forall k a.
Ord k =>
(Maybe a -> Maybe a) -> k -> Map k a -> Map k a
Map.alter Maybe a -> Maybe a
h k
k Map k a
old where
    h :: Maybe a -> Maybe a
h Maybe a
mb = case Maybe a
mb of
      Maybe a
Nothing -> a -> Maybe a
forall a. a -> Maybe a
Just (b -> a
f b
y)
      Just a
x  -> a -> Maybe a
forall a. a -> Maybe a
Just (b -> a -> a
g b
y a
x)    

--------------------------------------------------------------------------------
-- * Signs

data Sign
  = Plus                            -- hmm, this way @Plus < Minus@, not sure about that
  | Minus
  deriving (Sign -> Sign -> Bool
(Sign -> Sign -> Bool) -> (Sign -> Sign -> Bool) -> Eq Sign
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: Sign -> Sign -> Bool
$c/= :: Sign -> Sign -> Bool
== :: Sign -> Sign -> Bool
$c== :: Sign -> Sign -> Bool
Eq,Eq Sign
Eq Sign
-> (Sign -> Sign -> Ordering)
-> (Sign -> Sign -> Bool)
-> (Sign -> Sign -> Bool)
-> (Sign -> Sign -> Bool)
-> (Sign -> Sign -> Bool)
-> (Sign -> Sign -> Sign)
-> (Sign -> Sign -> Sign)
-> Ord Sign
Sign -> Sign -> Bool
Sign -> Sign -> Ordering
Sign -> Sign -> Sign
forall a.
Eq a
-> (a -> a -> Ordering)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> a)
-> (a -> a -> a)
-> Ord a
min :: Sign -> Sign -> Sign
$cmin :: Sign -> Sign -> Sign
max :: Sign -> Sign -> Sign
$cmax :: Sign -> Sign -> Sign
>= :: Sign -> Sign -> Bool
$c>= :: Sign -> Sign -> Bool
> :: Sign -> Sign -> Bool
$c> :: Sign -> Sign -> Bool
<= :: Sign -> Sign -> Bool
$c<= :: Sign -> Sign -> Bool
< :: Sign -> Sign -> Bool
$c< :: Sign -> Sign -> Bool
compare :: Sign -> Sign -> Ordering
$ccompare :: Sign -> Sign -> Ordering
$cp1Ord :: Eq Sign
Ord,Int -> Sign -> ShowS
[Sign] -> ShowS
Sign -> [Char]
(Int -> Sign -> ShowS)
-> (Sign -> [Char]) -> ([Sign] -> ShowS) -> Show Sign
forall a.
(Int -> a -> ShowS) -> (a -> [Char]) -> ([a] -> ShowS) -> Show a
showList :: [Sign] -> ShowS
$cshowList :: [Sign] -> ShowS
show :: Sign -> [Char]
$cshow :: Sign -> [Char]
showsPrec :: Int -> Sign -> ShowS
$cshowsPrec :: Int -> Sign -> ShowS
Show,ReadPrec [Sign]
ReadPrec Sign
Int -> ReadS Sign
ReadS [Sign]
(Int -> ReadS Sign)
-> ReadS [Sign] -> ReadPrec Sign -> ReadPrec [Sign] -> Read Sign
forall a.
(Int -> ReadS a)
-> ReadS [a] -> ReadPrec a -> ReadPrec [a] -> Read a
readListPrec :: ReadPrec [Sign]
$creadListPrec :: ReadPrec [Sign]
readPrec :: ReadPrec Sign
$creadPrec :: ReadPrec Sign
readList :: ReadS [Sign]
$creadList :: ReadS [Sign]
readsPrec :: Int -> ReadS Sign
$creadsPrec :: Int -> ReadS Sign
Read)

oppositeSign :: Sign -> Sign
oppositeSign :: Sign -> Sign
oppositeSign Sign
s = case Sign
s of
  Sign
Plus  -> Sign
Minus
  Sign
Minus -> Sign
Plus

mulSign :: Sign -> Sign -> Sign
mulSign :: Sign -> Sign -> Sign
mulSign Sign
s1 Sign
s2 = case Sign
s1 of
  Sign
Plus  -> Sign
s2
  Sign
Minus -> Sign -> Sign
oppositeSign Sign
s2

productOfSigns :: [Sign] -> Sign
productOfSigns :: [Sign] -> Sign
productOfSigns = Sign -> [Sign] -> Sign
go Sign
Plus where
  go :: Sign -> [Sign] -> Sign
go !Sign
acc []     = Sign
acc
  go !Sign
acc (Sign
x:[Sign]
xs) = case Sign
x of
    Sign
Plus  -> Sign -> [Sign] -> Sign
go Sign
acc [Sign]
xs
    Sign
Minus -> Sign -> [Sign] -> Sign
go (Sign -> Sign
oppositeSign Sign
acc) [Sign]
xs

-- | Negate the second argument if the first is odd
negateIfOdd :: (Integral a, Num b) => a -> b -> b
negateIfOdd :: a -> b -> b
negateIfOdd a
k b
y = if a -> Bool
forall a. Integral a => a -> Bool
even a
k then b
y else b -> b
forall a. Num a => a -> a
negate b
y

--------------------------------------------------------------------------------

-- Semigroup became a superclass of Monoid
#if MIN_VERSION_base(4,11,0)        

instance Semigroup Sign where
  <> :: Sign -> Sign -> Sign
(<>)    = Sign -> Sign -> Sign
mulSign
  sconcat :: NonEmpty Sign -> Sign
sconcat = (Sign -> Sign -> Sign) -> NonEmpty Sign -> Sign
forall (t :: * -> *) a. Foldable t => (a -> a -> a) -> t a -> a
foldl1 Sign -> Sign -> Sign
mulSign

instance Monoid Sign where
  mempty :: Sign
mempty  = Sign
Plus
  mconcat :: [Sign] -> Sign
mconcat = [Sign] -> Sign
productOfSigns

#else

instance Monoid Sign where
  mempty  = Plus
  mappend = mulSign
  mconcat = productOfSigns

#endif

--------------------------------------------------------------------------------

{-# SPECIALIZE signValue :: Sign -> Int     #-}
{-# SPECIALIZE signValue :: Sign -> Integer #-}

-- | @+1@ or @-1@
signValue :: Num a => Sign -> a
signValue :: Sign -> a
signValue Sign
s = case Sign
s of
  Sign
Plus  ->  a
1
  Sign
Minus -> -a
1

{-# SPECIALIZE signed :: Sign -> Int     -> Int     #-}
{-# SPECIALIZE signed :: Sign -> Integer -> Integer #-}

-- | Negate the second argument if the first is 'Minus'
signed :: Num a => Sign -> a -> a
signed :: Sign -> a -> a
signed Sign
s a
y = case Sign
s of
  Sign
Plus  -> a
y
  Sign
Minus -> a -> a
forall a. Num a => a -> a
negate a
y

class IsSigned a where
  signOf :: a -> Maybe Sign

signOfNum :: (Ord a, Num a) => a -> Maybe Sign 
signOfNum :: a -> Maybe Sign
signOfNum a
x = case a -> a -> Ordering
forall a. Ord a => a -> a -> Ordering
compare a
x a
0 of
  Ordering
LT -> Sign -> Maybe Sign
forall a. a -> Maybe a
Just Sign
Minus
  Ordering
GT -> Sign -> Maybe Sign
forall a. a -> Maybe a
Just Sign
Plus
  Ordering
EQ -> Maybe Sign
forall a. Maybe a
Nothing

instance IsSigned Int      where signOf :: Int -> Maybe Sign
signOf = Int -> Maybe Sign
forall a. (Ord a, Num a) => a -> Maybe Sign
signOfNum
instance IsSigned Integer  where signOf :: Integer -> Maybe Sign
signOf = Integer -> Maybe Sign
forall a. (Ord a, Num a) => a -> Maybe Sign
signOfNum
instance IsSigned Rational where signOf :: Rational -> Maybe Sign
signOf = Rational -> Maybe Sign
forall a. (Ord a, Num a) => a -> Maybe Sign
signOfNum

--------------------------------------------------------------------------------
-- * Numbers

fromRat :: Rational -> Integer
fromRat :: Rational -> Integer
fromRat Rational
r = case Rational -> Integer
forall a. Ratio a -> a
denominator Rational
r of
  Integer
1 -> Rational -> Integer
forall a. Ratio a -> a
numerator Rational
r
  Integer
_ -> [Char] -> Integer
forall a. HasCallStack => [Char] -> a
error [Char]
"fromRat: not an integer"    

safeDiv :: Integer -> Integer -> Integer
safeDiv :: Integer -> Integer -> Integer
safeDiv Integer
a Integer
b = case Integer -> Integer -> (Integer, Integer)
forall a. Integral a => a -> a -> (a, a)
divMod Integer
a Integer
b of
  (Integer
q,Integer
0) -> Integer
q
  (Integer
q,Integer
r) -> [Char] -> Integer
forall a. HasCallStack => [Char] -> a
error ([Char] -> Integer) -> [Char] -> Integer
forall a b. (a -> b) -> a -> b
$ [Char]
"saveDiv: " [Char] -> ShowS
forall a. [a] -> [a] -> [a]
++ Integer -> [Char]
forall a. Show a => a -> [Char]
show Integer
a [Char] -> ShowS
forall a. [a] -> [a] -> [a]
++ [Char]
" = " [Char] -> ShowS
forall a. [a] -> [a] -> [a]
++ Integer -> [Char]
forall a. Show a => a -> [Char]
show Integer
b [Char] -> ShowS
forall a. [a] -> [a] -> [a]
++ [Char]
" * " [Char] -> ShowS
forall a. [a] -> [a] -> [a]
++ Integer -> [Char]
forall a. Show a => a -> [Char]
show Integer
q [Char] -> ShowS
forall a. [a] -> [a] -> [a]
++ [Char]
" + " [Char] -> ShowS
forall a. [a] -> [a] -> [a]
++ Integer -> [Char]
forall a. Show a => a -> [Char]
show Integer
r

--------------------------------------------------------------------------------
-- * Basic number theory

-- | A000142.
factorial :: Integral a => a -> Integer
factorial :: a -> Integer
factorial a
n
  | a
n a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<  a
0    = [Char] -> Integer
forall a. HasCallStack => [Char] -> a
error [Char]
"factorial: input should be nonnegative"
  | a
n a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0    = Integer
1
  | Bool
otherwise = [Integer] -> Integer
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [Integer
1..a -> Integer
forall a b. (Integral a, Num b) => a -> b
fromIntegral a
n]

-- | A007318. Note: This is zero for @n<0@ or @k<0@; see also 'signedBinomial' below.
binomial :: Integral a => a -> a -> Integer
binomial :: a -> a -> Integer
binomial a
n a
k 
  | a
k a -> a -> Bool
forall a. Ord a => a -> a -> Bool
> a
n = Integer
0
  | a
k a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< a
0 = Integer
0
  | a
k a -> a -> Bool
forall a. Ord a => a -> a -> Bool
> (a
n a -> a -> a
forall a. Integral a => a -> a -> a
`div` a
2) = a -> a -> Integer
forall a. Integral a => a -> a -> Integer
binomial a
n (a
na -> a -> a
forall a. Num a => a -> a -> a
-a
k)
  | Bool
otherwise = ([Integer] -> Integer
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [Integer
n'Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
-Integer
k'Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+Integer
1 .. Integer
n']) Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`div` ([Integer] -> Integer
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [Integer
1..Integer
k'])
  where 
    k' :: Integer
k' = a -> Integer
forall a b. (Integral a, Num b) => a -> b
fromIntegral a
k
    n' :: Integer
n' = a -> Integer
forall a b. (Integral a, Num b) => a -> b
fromIntegral a
n

moebiusMu :: Num c => Int -> c
moebiusMu :: Int -> c
moebiusMu Int
n 
  | (Int -> Bool) -> [Int] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any (Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>Int
1) [Int]
expos       =  c
0
  | Int -> Bool
forall a. Integral a => a -> Bool
even ([Integer] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Integer]
primes) =  c
1
  | Bool
otherwise            = -c
1
  where
    factors :: [(Integer, Int)]
factors = [Integer] -> [(Integer, Int)]
groupIntegerFactors ([Integer] -> [(Integer, Int)]) -> [Integer] -> [(Integer, Int)]
forall a b. (a -> b) -> a -> b
$ Integer -> [Integer]
integerFactorsTrialDivision (Int -> Integer
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n)
    ([Integer]
primes,[Int]
expos) = [(Integer, Int)] -> ([Integer], [Int])
forall a b. [(a, b)] -> ([a], [b])
unzip [(Integer, Int)]
factors

divisors :: Int -> [Int]
divisors :: Int -> [Int]
divisors Int
n = [ [Int] -> Int
forall b. Integral b => [b] -> Int
f [Int]
tup | [Int]
tup <- [Int] -> [[Int]]
tuples' [Int]
expos ] where
  grps :: [(Integer, Int)]
grps = [Integer] -> [(Integer, Int)]
groupIntegerFactors ([Integer] -> [(Integer, Int)]) -> [Integer] -> [(Integer, Int)]
forall a b. (a -> b) -> a -> b
$ Integer -> [Integer]
integerFactorsTrialDivision (Integer -> [Integer]) -> Integer -> [Integer]
forall a b. (a -> b) -> a -> b
$ Int -> Integer
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n
  ([Integer]
primes,[Int]
expos) = [(Integer, Int)] -> ([Integer], [Int])
forall a b. [(a, b)] -> ([a], [b])
unzip [(Integer, Int)]
grps
  int_ps :: [Int]
int_ps = (Integer -> Int) -> [Integer] -> [Int]
forall a b. (a -> b) -> [a] -> [b]
map Integer -> Int
forall a. Num a => Integer -> a
fromInteger [Integer]
primes :: [Int]
  f :: [b] -> Int
f [b]
es = (Int -> Int -> Int) -> Int -> [Int] -> Int
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl' Int -> Int -> Int
forall a. Num a => a -> a -> a
(*) Int
1 ([Int] -> Int) -> [Int] -> Int
forall a b. (a -> b) -> a -> b
$ (Int -> b -> Int) -> [Int] -> [b] -> [Int]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Int -> b -> Int
forall a b. (Num a, Integral b) => a -> b -> a
(^) [Int]
int_ps [b]
es

-- | Square-free divisors together with their Mobius mu value
squareFreeDivisors :: Int -> [(Int,Sign)]
squareFreeDivisors :: Int -> [(Int, Sign)]
squareFreeDivisors Int
n = ([Int] -> (Int, Sign)) -> [[Int]] -> [(Int, Sign)]
forall a b. (a -> b) -> [a] -> [b]
map [Int] -> (Int, Sign)
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> (a, Sign)
f ([Int] -> [[Int]]
forall a. [a] -> [[a]]
sublists [Int]
int_ps) where
  grps :: [(Integer, Int)]
grps = [Integer] -> [(Integer, Int)]
groupIntegerFactors ([Integer] -> [(Integer, Int)]) -> [Integer] -> [(Integer, Int)]
forall a b. (a -> b) -> a -> b
$ Integer -> [Integer]
integerFactorsTrialDivision (Integer -> [Integer]) -> Integer -> [Integer]
forall a b. (a -> b) -> a -> b
$ Int -> Integer
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n
  primes :: [Integer]
primes = ((Integer, Int) -> Integer) -> [(Integer, Int)] -> [Integer]
forall a b. (a -> b) -> [a] -> [b]
map (Integer, Int) -> Integer
forall a b. (a, b) -> a
fst [(Integer, Int)]
grps
  int_ps :: [Int]
int_ps = (Integer -> Int) -> [Integer] -> [Int]
forall a b. (a -> b) -> [a] -> [b]
map Integer -> Int
forall a. Num a => Integer -> a
fromInteger [Integer]
primes :: [Int]
  f :: t a -> (a, Sign)
f t a
ps = ( (a -> a -> a) -> a -> t a -> a
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl' a -> a -> a
forall a. Num a => a -> a -> a
(*) a
1 t a
ps , if Int -> Bool
forall a. Integral a => a -> Bool
even (t a -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length t a
ps) then Sign
Plus else Sign
Minus)

-- | List of primes, using tree merge with wheel. Code by Will Ness.
primes :: [Integer]
primes :: [Integer]
primes = Integer
2Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
:Integer
3Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
:Integer
5Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
:Integer
7Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
: Integer -> [Integer] -> [Integer] -> [Integer]
forall a. (Eq a, Num a) => a -> [a] -> [a] -> [a]
gaps Integer
11 [Integer]
wheel ([[Integer]] -> [Integer]
forall a. Ord a => [[a]] -> [a]
fold3t ([[Integer]] -> [Integer]) -> [[Integer]] -> [Integer]
forall a b. (a -> b) -> a -> b
$ Integer -> [Integer] -> [Integer] -> [[Integer]]
forall a. (Eq a, Num a) => a -> [a] -> [a] -> [[a]]
roll Integer
11 [Integer]
wheel [Integer]
primes') where                                                             

  primes' :: [Integer]
primes' = Integer
11Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
: Integer -> [Integer] -> [Integer] -> [Integer]
forall a. (Eq a, Num a) => a -> [a] -> [a] -> [a]
gaps Integer
13 ([Integer] -> [Integer]
forall a. [a] -> [a]
tail [Integer]
wheel) ([[Integer]] -> [Integer]
forall a. Ord a => [[a]] -> [a]
fold3t ([[Integer]] -> [Integer]) -> [[Integer]] -> [Integer]
forall a b. (a -> b) -> a -> b
$ Integer -> [Integer] -> [Integer] -> [[Integer]]
forall a. (Eq a, Num a) => a -> [a] -> [a] -> [[a]]
roll Integer
11 [Integer]
wheel [Integer]
primes')
  fold3t :: [[a]] -> [a]
fold3t ((a
x:[a]
xs): ~([a]
ys:[a]
zs:[[a]]
t)) 
    = a
x a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [a] -> [a] -> [a]
forall a. Ord a => [a] -> [a] -> [a]
union [a]
xs ([a] -> [a] -> [a]
forall a. Ord a => [a] -> [a] -> [a]
union [a]
ys [a]
zs) [a] -> [a] -> [a]
forall a. Ord a => [a] -> [a] -> [a]
`union` [[a]] -> [a]
fold3t ([[a]] -> [[a]]
forall a. Ord a => [[a]] -> [[a]]
pairs [[a]]
t)            
  pairs :: [[a]] -> [[a]]
pairs ((a
x:[a]
xs):[a]
ys:[[a]]
t) = (a
x a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [a] -> [a] -> [a]
forall a. Ord a => [a] -> [a] -> [a]
union [a]
xs [a]
ys) [a] -> [[a]] -> [[a]]
forall a. a -> [a] -> [a]
: [[a]] -> [[a]]
pairs [[a]]
t 
  wheel :: [Integer]
wheel = Integer
2Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
:Integer
4Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
:Integer
2Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
:Integer
4Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
:Integer
6Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
:Integer
2Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
:Integer
6Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
:Integer
4Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
:Integer
2Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
:Integer
4Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
:Integer
6Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
:Integer
6Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
:Integer
2Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
:Integer
6Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
:Integer
4Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
:Integer
2Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
:Integer
6Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
:Integer
4Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
:Integer
6Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
:Integer
8Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
:Integer
4Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
:Integer
2Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
:Integer
4Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
:Integer
2Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
:  
          Integer
4Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
:Integer
8Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
:Integer
6Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
:Integer
4Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
:Integer
6Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
:Integer
2Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
:Integer
4Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
:Integer
6Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
:Integer
2Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
:Integer
6Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
:Integer
6Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
:Integer
4Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
:Integer
2Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
:Integer
4Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
:Integer
6Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
:Integer
2Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
:Integer
6Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
:Integer
4Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
:Integer
2Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
:Integer
4Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
:Integer
2Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
:Integer
10Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
:Integer
2Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
:Integer
10Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
:[Integer]
wheel 
  gaps :: a -> [a] -> [a] -> [a]
gaps a
k ws :: [a]
ws@(a
w:[a]
t) cs :: [a]
cs@ ~(a
c:[a]
u) 
    | a
ka -> a -> Bool
forall a. Eq a => a -> a -> Bool
==a
c  = a -> [a] -> [a] -> [a]
gaps (a
ka -> a -> a
forall a. Num a => a -> a -> a
+a
w) [a]
t [a]
u              
    | Bool
True  = a
k a -> [a] -> [a]
forall a. a -> [a] -> [a]
: a -> [a] -> [a] -> [a]
gaps (a
ka -> a -> a
forall a. Num a => a -> a -> a
+a
w) [a]
t [a]
cs  
  roll :: a -> [a] -> [a] -> [[a]]
roll a
k ws :: [a]
ws@(a
w:[a]
t) ps :: [a]
ps@ ~(a
p:[a]
u) 
    | a
ka -> a -> Bool
forall a. Eq a => a -> a -> Bool
==a
p  = (a -> a -> a) -> a -> [a] -> [a]
forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl (\a
c a
d->a
ca -> a -> a
forall a. Num a => a -> a -> a
+a
pa -> a -> a
forall a. Num a => a -> a -> a
*a
d) (a
pa -> a -> a
forall a. Num a => a -> a -> a
*a
p) [a]
ws [a] -> [[a]] -> [[a]]
forall a. a -> [a] -> [a]
: a -> [a] -> [a] -> [[a]]
roll (a
ka -> a -> a
forall a. Num a => a -> a -> a
+a
w) [a]
t [a]
u              
    | Bool
True  = a -> [a] -> [a] -> [[a]]
roll (a
ka -> a -> a
forall a. Num a => a -> a -> a
+a
w) [a]
t [a]
ps   

  minus :: [a] -> [a] -> [a]
minus xxs :: [a]
xxs@(a
x:[a]
xs) yys :: [a]
yys@(a
y:[a]
ys) = case a -> a -> Ordering
forall a. Ord a => a -> a -> Ordering
compare a
x a
y of 
    Ordering
LT -> a
x a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [a] -> [a] -> [a]
minus [a]
xs  [a]
yys
    Ordering
EQ ->     [a] -> [a] -> [a]
minus [a]
xs  [a]
ys 
    Ordering
GT ->     [a] -> [a] -> [a]
minus [a]
xxs [a]
ys
  minus [a]
xs [] = [a]
xs
  minus [] [a]
_  = []
  
  union :: [a] -> [a] -> [a]
union xxs :: [a]
xxs@(a
x:[a]
xs) yys :: [a]
yys@(a
y:[a]
ys) = case a -> a -> Ordering
forall a. Ord a => a -> a -> Ordering
compare a
x a
y of 
    Ordering
LT -> a
x a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [a] -> [a] -> [a]
union [a]
xs  [a]
yys
    Ordering
EQ -> a
x a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [a] -> [a] -> [a]
union [a]
xs  [a]
ys 
    Ordering
GT -> a
y a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [a] -> [a] -> [a]
union [a]
xxs [a]
ys
  union [a]
xs [] = [a]
xs
  union [] [a]
ys =[a]
ys

--------------------------------------------------------------------------------
-- Prime factorization

-- | Groups integer factors. Example: from [2,2,2,3,3,5] we produce [(2,3),(3,2),(5,1)]  
groupIntegerFactors :: [Integer] -> [(Integer,Int)]
groupIntegerFactors :: [Integer] -> [(Integer, Int)]
groupIntegerFactors = ([Integer] -> (Integer, Int)) -> [[Integer]] -> [(Integer, Int)]
forall a b. (a -> b) -> [a] -> [b]
map [Integer] -> (Integer, Int)
forall a. [a] -> (a, Int)
f ([[Integer]] -> [(Integer, Int)])
-> ([Integer] -> [[Integer]]) -> [Integer] -> [(Integer, Int)]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Integer] -> [[Integer]]
forall a. Eq a => [a] -> [[a]]
group ([Integer] -> [[Integer]])
-> ([Integer] -> [Integer]) -> [Integer] -> [[Integer]]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Integer] -> [Integer]
forall a. Ord a => [a] -> [a]
sort where
  f :: [a] -> (a, Int)
f [a]
xs = ([a] -> a
forall a. [a] -> a
head [a]
xs, [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
xs)

-- | The naive trial division algorithm.
integerFactorsTrialDivision :: Integer -> [Integer]
integerFactorsTrialDivision :: Integer -> [Integer]
integerFactorsTrialDivision Integer
n 
  | Integer
nInteger -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<Integer
1 = [Char] -> [Integer]
forall a. HasCallStack => [Char] -> a
error [Char]
"integerFactorsTrialDivision: n should be at least 1"
  | Bool
otherwise = [Integer] -> Integer -> [Integer]
forall a. Integral a => [a] -> a -> [a]
go [Integer]
primes Integer
n 
  where
    go :: [a] -> a -> [a]
go [a]
_  a
1 = []
    go [a]
rs a
k = [a] -> a -> [a]
sub [a]
ps a
k where
      sub :: [a] -> a -> [a]
sub [] a
k = [a
k]
      sub qqs :: [a]
qqs@(a
q:[a]
qs) a
k = case a -> a -> a
forall a. Integral a => a -> a -> a
mod a
k a
q of
        a
0 -> a
q a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [a] -> a -> [a]
go [a]
qqs (a -> a -> a
forall a. Integral a => a -> a -> a
div a
k a
q)
        a
_ -> [a] -> a -> [a]
sub [a]
qs a
k
      ps :: [a]
ps = (a -> Bool) -> [a] -> [a]
forall a. (a -> Bool) -> [a] -> [a]
takeWhile (\a
p -> a
pa -> a -> a
forall a. Num a => a -> a -> a
*a
p a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<= a
k) [a]
rs  

--------------------------------------------------------------------------------
-- * Basic combinatorics

tuples' :: [Int] -> [[Int]]
tuples' :: [Int] -> [[Int]]
tuples' [] = [[]]
tuples' (Int
s:[Int]
ss) = [ Int
xInt -> [Int] -> [Int]
forall a. a -> [a] -> [a]
:[Int]
xs | Int
x <- [Int
0..Int
s] , [Int]
xs <- [Int] -> [[Int]]
tuples' [Int]
ss ] 

-- | All sublists of a list.
sublists :: [a] -> [[a]]
sublists :: [a] -> [[a]]
sublists [] = [[]]
sublists (a
x:[a]
xs) = [a] -> [[a]]
forall a. [a] -> [[a]]
sublists [a]
xs [[a]] -> [[a]] -> [[a]]
forall a. [a] -> [a] -> [a]
++ ([a] -> [a]) -> [[a]] -> [[a]]
forall a b. (a -> b) -> [a] -> [b]
map (a
xa -> [a] -> [a]
forall a. a -> [a] -> [a]
:) ([a] -> [[a]]
forall a. [a] -> [[a]]
sublists [a]
xs)  

--------------------------------------------------------------------------------
-- * Integer-indexed cache

intCache :: ((Int -> a) -> (Int -> a)) -> (Int -> a)
intCache :: ((Int -> a) -> Int -> a) -> Int -> a
intCache (Int -> a) -> Int -> a
compute = Int -> a
cached where
  cached :: Int -> a
cached Int
n = Int -> ITable a -> a
forall a. Int -> ITable a -> a
lkpITable Int
n ITable a
table
  table :: ITable a
table    = [a] -> ITable a
forall a. [a] -> ITable a
mkITable ((Int -> a) -> [Int] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map ((Int -> a) -> Int -> a
compute Int -> a
cached) [Int
0..])
  
newtype ITable a = ITable [Array Int a] 

mkITable :: [a] -> ITable a
mkITable :: [a] -> ITable a
mkITable = [Array Int a] -> ITable a
forall a. [Array Int a] -> ITable a
ITable ([Array Int a] -> ITable a)
-> ([a] -> [Array Int a]) -> [a] -> ITable a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> [a] -> [Array Int a]
forall a. Int -> [a] -> [Array Int a]
go Int
1 where
  go :: Int -> [a] -> [Array Int a]
go !Int
siz [a]
list = Array Int a
arr Array Int a -> [Array Int a] -> [Array Int a]
forall a. a -> [a] -> [a]
: Int -> [a] -> [Array Int a]
go (Int
2Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
siz) [a]
rest where
    ([a]
this,[a]
rest) = Int -> [a] -> ([a], [a])
forall a. Int -> [a] -> ([a], [a])
splitAt Int
siz [a]
list
    arr :: Array Int a
arr = (Int, Int) -> [a] -> Array Int a
forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (Int
0,Int
sizInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1) [a]
this

lkpITable :: Int -> ITable a -> a        
lkpITable :: Int -> ITable a -> a
lkpITable Int
idx (ITable [Array Int a]
list) = Int -> Int -> [Array Int a] -> a
forall t p. (Ix t, Num t) => t -> t -> [Array t p] -> p
go Int
1 Int
idx [Array Int a]
list where
  go :: t -> t -> [Array t p] -> p
go !t
siz !t
idx (Array t p
this:[Array t p]
rest) = if t
idx t -> t -> Bool
forall a. Ord a => a -> a -> Bool
< t
siz
    then (Array t p
this Array t p -> t -> p
forall i e. Ix i => Array i e -> i -> e
! t
idx)
    else t -> t -> [Array t p] -> p
go (t
2t -> t -> t
forall a. Num a => a -> a -> a
*t
siz) (t
idxt -> t -> t
forall a. Num a => a -> a -> a
-t
siz) [Array t p]
rest

--------------------------------------------------------------------------------
-- * Stirling numbers

-- | Rows of (signed) Stirling numbers of the first kind. OEIS:A008275.
-- Coefficients of the polinomial @(x-1)*(x-2)*...*(x-n+1)@.
-- This function uses the recursion formula.
signedStirling1stArray :: Integral a => a -> Array Int Integer
signedStirling1stArray :: a -> Array Int Integer
signedStirling1stArray a
n
  | a
n a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<  a
1    = [Char] -> Array Int Integer
forall a. HasCallStack => [Char] -> a
error [Char]
"stirling1stArray: n should be at least 1"
  | a
n a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
1    = (Int, Int) -> [Integer] -> Array Int Integer
forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (Int
1,Int
1 ) [Integer
1]
  | Bool
otherwise = (Int, Int) -> [Integer] -> Array Int Integer
forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (Int
1,Int
n') [ Int -> Integer
lkp (Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1) Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- a -> Integer
forall a b. (Integral a, Num b) => a -> b
fromIntegral (a
na -> a -> a
forall a. Num a => a -> a -> a
-a
1) Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Int -> Integer
lkp Int
k | Int
k<-[Int
1..Int
n'] ] 
  where
    prev :: Array Int Integer
prev = a -> Array Int Integer
forall a. Integral a => a -> Array Int Integer
signedStirling1stArray (a
na -> a -> a
forall a. Num a => a -> a -> a
-a
1)
    n' :: Int
n' = a -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral a
n :: Int
    lkp :: Int -> Integer
lkp Int
j | Int
j Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<  Int
1    = Integer
0
          | Int
j Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
n'   = Integer
0
          | Bool
otherwise = Array Int Integer
prev Array Int Integer -> Int -> Integer
forall i e. Ix i => Array i e -> i -> e
! Int
j 

-- | Stirling numbers of the second kind. OEIS:A008277.
-- This function uses an explicit formula.
-- 
-- Argument order: @stirling2nd n k@
--
stirling2nd :: Integral a => a -> a -> Integer
stirling2nd :: a -> a -> Integer
stirling2nd a
n a
k 
  | a
ka -> a -> Bool
forall a. Eq a => a -> a -> Bool
==a
0 Bool -> Bool -> Bool
&& a
na -> a -> Bool
forall a. Eq a => a -> a -> Bool
==a
0 = Integer
1
  | a
k a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< a
1        = Integer
0
  | a
k a -> a -> Bool
forall a. Ord a => a -> a -> Bool
> a
n        = Integer
0
  | Bool
otherwise = [Integer] -> Integer
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Integer]
xs Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`div` a -> Integer
forall a. Integral a => a -> Integer
factorial a
k where
      xs :: [Integer]
xs = [ a -> Integer -> Integer
forall a b. (Integral a, Num b) => a -> b -> b
negateIfOdd (a
ka -> a -> a
forall a. Num a => a -> a -> a
-a
i) (Integer -> Integer) -> Integer -> Integer
forall a b. (a -> b) -> a -> b
$ a -> a -> Integer
forall a. Integral a => a -> a -> Integer
binomial a
k a
i Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* (a -> Integer
forall a b. (Integral a, Num b) => a -> b
fromIntegral a
i)Integer -> a -> Integer
forall a b. (Num a, Integral b) => a -> b -> a
^a
n | a
i<-[a
0..a
k] ]

--------------------------------------------------------------------------------