-- | A few important number sequences. 
--  
-- See the \"On-Line Encyclopedia of Integer Sequences\",
-- <https://oeis.org> .

module Math.Combinat.Numbers where

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

import Data.Array

import Math.Combinat.Helper ( sum' )
import Math.Combinat.Sign

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

-- | 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]

-- | A006882.
doubleFactorial :: Integral a => a -> Integer
doubleFactorial :: a -> Integer
doubleFactorial 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]
"doubleFactorial: input should be nonnegative"
  | a
n a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0    = Integer
1
  | a -> Bool
forall a. Integral a => a -> Bool
odd a
n     = [Integer] -> Integer
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [Integer
1,Integer
3..a -> Integer
forall a b. (Integral a, Num b) => a -> b
fromIntegral a
n]
  | Bool
otherwise = [Integer] -> Integer
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [Integer
2,Integer
4..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

-- | The extension of the binomial function to negative inputs. This should satisfy the following properties:
--
-- > for n,k >=0 : signedBinomial n k == binomial n k
-- > for any n,k : signedBinomial n k == signedBinomial n (n-k) 
-- > for k >= 0  : signedBinomial (-n) k == (-1)^k * signedBinomial (n+k-1) k
--
-- Note: This is compatible with Mathematica's @Binomial@ function.
--
signedBinomial :: Int -> Int -> Integer
signedBinomial :: Int -> Int -> Integer
signedBinomial Int
n Int
k
  | Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
0     = Int -> Int -> Integer
forall a. Integral a => a -> a -> Integer
binomial Int
n Int
k
  | Int
k Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
0     = Int -> Integer -> Integer
forall a b. (Integral a, Num b) => a -> b -> b
negateIfOdd    Int
k  (Integer -> Integer) -> Integer -> Integer
forall a b. (a -> b) -> a -> b
$ Int -> Int -> Integer
forall a. Integral a => a -> a -> Integer
binomial (Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1)   Int
k  
  | Bool
otherwise  = Int -> Integer -> Integer
forall a b. (Integral a, Num b) => a -> b -> b
negateIfOdd (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
k) (Integer -> Integer) -> Integer -> Integer
forall a b. (a -> b) -> a -> b
$ Int -> Int -> Integer
forall a. Integral a => a -> a -> Integer
binomial (-Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1) (-Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1)

{-
test_signed_0 = [ signedBinomial ( n) k == signedBinomial ( n) ( n-k)                | n<-[-30..40] , k<-[-30..40] ]
test_signed_1 = [ signedBinomial (-n) k == signedBinomial (-n) (-n-k)                | n<-[-30..40] , k<-[-30..40] ]
test_signed_2 = [ signedBinomial (-n) k == negateIfOdd k $ signedBinomial (n+k-1) k  | n<-[-30..40] , k<-[0..30] ]
-}

-- | A given row of the Pascal triangle; equivalent to a sequence of binomial 
-- numbers, but much more efficient. You can also left-fold over it.
--
-- > pascalRow n == [ binomial n k | k<-[0..n] ]
pascalRow :: Integral a => a -> [Integer]
pascalRow :: a -> [Integer]
pascalRow a
n' = Integer -> Integer -> [Integer]
worker Integer
0 Integer
1 where
  n :: Integer
n = a -> Integer
forall a b. (Integral a, Num b) => a -> b
fromIntegral a
n'
  worker :: Integer -> Integer -> [Integer]
worker Integer
j Integer
x
    | Integer
jInteger -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
>Integer
n   = [] 
    | Bool
True  = let j' :: Integer
j'=Integer
jInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
+Integer
1 in Integer
x Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
: Integer -> Integer -> [Integer]
worker Integer
j' (Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
div (Integer
xInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*(Integer
nInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
-Integer
j)) Integer
j') 

multinomial :: Integral a => [a] -> Integer
multinomial :: [a] -> Integer
multinomial [a]
xs = Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
div
  (a -> Integer
forall a. Integral a => a -> Integer
factorial ([a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [a]
xs))
  ([Integer] -> Integer
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [ a -> Integer
forall a. Integral a => a -> Integer
factorial a
x | a
x<-[a]
xs ])  
  
--------------------------------------------------------------------------------
-- * Catalan numbers

-- | Catalan numbers. OEIS:A000108.
catalan :: Integral a => a -> Integer
catalan :: a -> Integer
catalan a
n 
  | a
n a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< a
0     = Integer
0
  | Bool
otherwise = a -> a -> Integer
forall a. Integral a => a -> a -> Integer
binomial (a
na -> a -> a
forall a. Num a => a -> a -> a
+a
n) a
n Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`div` 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)

-- | Catalan's triangle. OEIS:A009766.
-- Note:
--
-- > catalanTriangle n n == catalan n
-- > catalanTriangle n k == countStandardYoungTableaux (toPartition [n,k])
--
catalanTriangle :: Integral a => a -> a -> Integer
catalanTriangle :: a -> a -> Integer
catalanTriangle 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
  | Bool
otherwise = (a -> a -> Integer
forall a. Integral a => a -> a -> Integer
binomial (a
na -> a -> a
forall a. Num a => a -> a -> a
+a
k) a
n 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
ka -> a -> a
forall a. Num a => a -> a -> a
+a
1)) Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`div` 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)

--------------------------------------------------------------------------------
-- * 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 
        
-- | (Signed) Stirling numbers of the first kind. OEIS:A008275.
-- This function uses 'signedStirling1stArray', so it shouldn't be used
-- to compute /many/ Stirling numbers.
--
-- Argument order: @signedStirling1st n k@
--
signedStirling1st :: Integral a => a -> a -> Integer
signedStirling1st :: a -> a -> Integer
signedStirling1st 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    = a -> Array Int Integer
forall a. Integral a => a -> Array Int Integer
signedStirling1stArray a
n Array Int Integer -> Int -> Integer
forall i e. Ix i => Array i e -> i -> e
! (a -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral a
k)

-- | (Unsigned) Stirling numbers of the first kind. See 'signedStirling1st'.
unsignedStirling1st :: Integral a => a -> a -> Integer
unsignedStirling1st :: a -> a -> Integer
unsignedStirling1st a
n a
k = Integer -> Integer
forall a. Num a => a -> a
abs (a -> a -> Integer
forall a. Integral a => a -> a -> Integer
signedStirling1st a
n a
k)

-- | 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] ]

--------------------------------------------------------------------------------
-- * Bernoulli numbers

-- | Bernoulli numbers. @bernoulli 1 == -1%2@ and @bernoulli k == 0@ for
-- k>2 and /odd/. This function uses the formula involving Stirling numbers
-- of the second kind. Numerators: A027641, denominators: A027642.
bernoulli :: Integral a => a -> Rational
bernoulli :: a -> Rational
bernoulli a
n 
  | a
n a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<  a
0    = [Char] -> Rational
forall a. HasCallStack => [Char] -> a
error [Char]
"bernoulli: n should be nonnegative"
  | a
n a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0    = Rational
1
  | a
n a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
1    = -Rational
1Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/Rational
2
  | Bool
otherwise = [Rational] -> Rational
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [ a -> Rational
f a
k | a
k<-[a
1..a
n] ] 
  where
    f :: a -> Rational
f a
k = Integer -> Rational
forall a. Real a => a -> Rational
toRational (a -> Integer -> Integer
forall a b. (Integral a, Num b) => a -> b -> b
negateIfOdd (a
na -> a -> a
forall a. Num a => a -> a -> a
+a
k) (Integer -> Integer) -> Integer -> Integer
forall a b. (a -> b) -> a -> b
$ a -> Integer
forall a. Integral a => a -> Integer
factorial a
k Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* a -> a -> Integer
forall a. Integral a => a -> a -> Integer
stirling2nd a
n a
k) 
        Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/ a -> Rational
forall a. Real a => a -> Rational
toRational (a
ka -> a -> a
forall a. Num a => a -> a -> a
+a
1)

--------------------------------------------------------------------------------
-- * Bell numbers

-- | Bell numbers (Sloane's A000110) from B(0) up to B(n). B(0)=B(1)=1, B(2)=2, etc. 
--
-- The Bell numbers count the number of /set partitions/ of a set of size @n@
-- 
-- See <http://en.wikipedia.org/wiki/Bell_number>
--
bellNumbersArray :: Integral a => a -> Array Int Integer
bellNumbersArray :: a -> Array Int Integer
bellNumbersArray a
nn = Array Int Integer
arr where
  arr :: Array Int Integer
arr = (Int, Int) -> [(Int, Integer)] -> Array Int Integer
forall i e. Ix i => (i, i) -> [(i, e)] -> Array i e
array (Int
0::Int,Int
n) [(Int, Integer)]
kvs 
  n :: Int
n = a -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral a
nn :: Int
  kvs :: [(Int, Integer)]
kvs = (Int
0,Integer
1) (Int, Integer) -> [(Int, Integer)] -> [(Int, Integer)]
forall a. a -> [a] -> [a]
: [ (Int
k, Int -> Integer
f Int
k) | Int
k<-[Int
1..Int
n] ] 
  f :: Int -> Integer
f Int
n = [Integer] -> Integer
forall a. Num a => [a] -> a
sum' [ Int -> Int -> Integer
forall a. Integral a => a -> a -> Integer
binomial (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1) Int
k Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Array Int Integer
arr Array Int Integer -> Int -> Integer
forall i e. Ix i => Array i e -> i -> e
! Int
k | Int
k<-[Int
0..Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1] ]

-- | The n-th Bell number B(n), using the Stirling numbers of the second kind.
-- This may be slower than using 'bellNumbersArray'.
bellNumber :: Integral a => a -> Integer
bellNumber :: a -> Integer
bellNumber a
nn
  | Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<  Int
0     = [Char] -> Integer
forall a. HasCallStack => [Char] -> a
error [Char]
"bellNumber: expecting a nonnegative index"
  | Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0     = Integer
1
  | Bool
otherwise  = [Integer] -> Integer
forall a. Num a => [a] -> a
sum' [ Int -> Int -> Integer
forall a. Integral a => a -> a -> Integer
stirling2nd Int
n Int
k | Int
k<-[Int
1..Int
n] ] 
  where
    n :: Int
n = a -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral a
nn :: Int

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