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

{-# LANGUAGE BangPatterns, ScopedTypeVariables #-}
module Math.Combinat.Numbers.Sequences where

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

import Data.Array
import Data.Bits ( shiftL , shiftR , (.&.) )

import Math.Combinat.Helper 
import Math.Combinat.Sign

import Math.Combinat.Numbers.Primes ( primes , factorize , productOfFactors )

import qualified Data.Map.Strict as Map   -- used for factorialPrimeExponentsNaive

--------------------------------------------------------------------------------
-- * Factorial

-- | The factorial function (A000142).
factorial :: Integral a => a -> Integer
factorial :: forall a. Integral a => a -> Integer
factorial = forall a. Integral a => a -> Integer
factorialSplit

-- | Faster implementation of the factorial function
factorialSplit :: Integral a => a -> Integer
factorialSplit :: forall a. Integral a => a -> Integer
factorialSplit a
n = forall a. Integral a => a -> a -> Integer
productFromTo a
1 a
n

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

-- | \"Swing factorial\" algorithm
factorialSwing :: Integral a => a -> Integer
factorialSwing :: forall a. Integral a => a -> Integer
factorialSwing a
n = [(Integer, Int)] -> Integer
productOfFactors (Int -> [(Integer, Int)]
factorialPrimeExponents forall a b. (a -> b) -> a -> b
$ forall a b. (Integral a, Num b) => a -> b
fromIntegral a
n) where

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

-- | Prime factorization of the factorial (using the \"swing factorial\" algorithm)
factorialPrimeExponents :: Int -> [(Integer,Int)]
factorialPrimeExponents :: Int -> [(Integer, Int)]
factorialPrimeExponents Int
n = forall a. (a -> Bool) -> [a] -> [a]
filter forall {a} {a}. (Ord a, Num a) => (a, a) -> Bool
cond forall a b. (a -> b) -> a -> b
$ forall a b. [a] -> [b] -> [(a, b)]
zip [Integer]
primes (Int -> [Int]
factorialPrimeExponents_ Int
n) where
  cond :: (a, a) -> Bool
cond (a
_,!a
e) = a
e forall a. Ord a => a -> a -> Bool
> a
0

factorialPrimeExponentsNaive :: forall a. Integral a => a -> [(Integer,Int)]
factorialPrimeExponentsNaive :: forall a. Integral a => a -> [(Integer, Int)]
factorialPrimeExponentsNaive a
n = [(Integer, Int)]
result where
  fi :: a -> Integer
fi = forall a b. (Integral a, Num b) => a -> b
fromIntegral :: a -> Integer
  result :: [(Integer, Int)]
result = forall k a. Map k a -> [(k, a)]
Map.toList 
         forall a b. (a -> b) -> a -> b
$ forall (f :: * -> *) k a.
(Foldable f, Ord k) =>
(a -> a -> a) -> f (Map k a) -> Map k a
Map.unionsWith forall a. Num a => a -> a -> a
(+) 
         forall a b. (a -> b) -> a -> b
$ forall a b. (a -> b) -> [a] -> [b]
map forall k a. Ord k => [(k, a)] -> Map k a
Map.fromList 
         forall a b. (a -> b) -> a -> b
$ forall a b. (a -> b) -> [a] -> [b]
map Integer -> [(Integer, Int)]
factorize 
         forall a b. (a -> b) -> a -> b
$ forall a b. (a -> b) -> [a] -> [b]
map a -> Integer
fi [a
1..a
n] 

factorialPrimeExponents_ :: Int -> [Int]
factorialPrimeExponents_ :: Int -> [Int]
factorialPrimeExponents_ = Int -> [Int]
go where
  go :: Int -> [Int]
go  Int
0 = []
  go  Int
1 = []
  go  Int
2 = [Int
1]
  go !Int
n = [Int] -> [Int] -> [Int]
longAdd [Int]
half [Int]
swing where
    half :: [Int]
half  = forall a b. (a -> b) -> [a] -> [b]
map (forall a b c. (a -> b -> c) -> b -> a -> c
flip forall a. Bits a => a -> Int -> a
shiftL Int
1) forall a b. (a -> b) -> a -> b
$ Int -> [Int]
go (forall a. Bits a => a -> Int -> a
shiftR Int
n Int
1)
    swing :: [Int]
swing = Int -> [Int]
swingFactorialExponents_ Int
n

  longAdd :: [Int] -> [Int] -> [Int]
  longAdd :: [Int] -> [Int] -> [Int]
longAdd [Int]
xs [] = [Int]
xs
  longAdd [] [Int]
ys = [Int]
ys
  longAdd (!Int
x:[Int]
xs) (!Int
y:[Int]
ys) = (Int
xforall a. Num a => a -> a -> a
+Int
y) forall a. a -> [a] -> [a]
: [Int] -> [Int] -> [Int]
longAdd [Int]
xs [Int]
ys

-- | Prime factorizaiton of the \"swing factorial\" function)
swingFactorialExponents_ :: Int -> [Int]
swingFactorialExponents_ :: Int -> [Int]
swingFactorialExponents_ = Int -> [Int]
go where
  go :: Int -> [Int]
go Int
0 = []
  go Int
1 = []
  go Int
2 = [Int
1]
  go Int
n = Int
expo2 forall a. a -> [a] -> [a]
: forall a b. (a -> b) -> [a] -> [b]
map Integer -> Int
expo (forall a. [a] -> [a]
tail [Integer]
ps) where

    nn :: Integer
nn = forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n :: Integer

    ps :: [Integer]
    ps :: [Integer]
ps = forall a. (a -> Bool) -> [a] -> [a]
takeWhile (forall a. Ord a => a -> a -> Bool
<=Integer
nn) [Integer]
primes 

    expo2 :: Int
    expo2 :: Int
expo2 = Int -> Int -> Int
go Int
0 (forall a. Bits a => a -> Int -> a
shiftR Int
n Int
1) where
      go :: Int -> Int -> Int
      go :: Int -> Int -> Int
go !Int
acc !Int
r  
        | Int
r forall a. Ord a => a -> a -> Bool
< Int
1     = Int
acc
        | Bool
otherwise = Int -> Int -> Int
go Int
acc' Int
r' 
        where
          acc' :: Int
acc' = Int
acc forall a. Num a => a -> a -> a
+ (Int
r forall a. Bits a => a -> a -> a
.&. Int
1)
          r' :: Int
r'   = forall a. Bits a => a -> Int -> a
shiftR Int
r Int
1

    expo :: Integer -> Int
    expo :: Integer -> Int
expo Integer
pp = Int -> Int -> Int
go Int
0 (forall a. Integral a => a -> a -> a
div Int
n Int
p) where
      p :: Int
p = forall a. Num a => Integer -> a
fromInteger Integer
pp :: Int
      go :: Int -> Int -> Int
      go :: Int -> Int -> Int
go !Int
acc !Int
r  
        | Int
r forall a. Ord a => a -> a -> Bool
< Int
1     = Int
acc
        | Bool
otherwise = Int -> Int -> Int
go Int
acc' Int
r' 
        where
          acc' :: Int
acc' = Int
acc forall a. Num a => a -> a -> a
+ (Int
r forall a. Bits a => a -> a -> a
.&. Int
1)
          r' :: Int
r'   = forall a. Integral a => a -> a -> a
div Int
r Int
p

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

-- | The double factorial
doubleFactorial :: Integral a => a -> Integer
doubleFactorial :: forall a. Integral a => a -> Integer
doubleFactorial = forall a. Integral a => a -> Integer
doubleFactorialSplit

-- | Faster implementation of the double factorial function
doubleFactorialSplit :: Integral a => a -> Integer
doubleFactorialSplit :: forall a. Integral a => a -> Integer
doubleFactorialSplit a
n
  | a
n forall a. Ord a => a -> a -> Bool
<  a
0    = forall a. HasCallStack => [Char] -> a
error [Char]
"doubleFactorialSplit: input should be nonnegative"
  | a
n forall a. Eq a => a -> a -> Bool
== a
0    = Integer
1
  | forall a. Integral a => a -> Bool
odd a
n     = forall a. Integral a => a -> a -> Integer
productFromToStride2 a
2 a
n
  | Bool
otherwise = let halfn :: a
halfn = forall a. Integral a => a -> a -> a
div a
n a
2 
                in  forall a. Bits a => a -> Int -> a
shiftL (forall a. Integral a => a -> Integer
factorialSplit a
halfn) (forall a b. (Integral a, Num b) => a -> b
fromIntegral a
halfn)

-- | Naive implementation of the double factorial (A006882).
doubleFactorialNaive :: Integral a => a -> Integer
doubleFactorialNaive :: forall a. Integral a => a -> Integer
doubleFactorialNaive a
n
  | a
n forall a. Ord a => a -> a -> Bool
<  a
0    = forall a. HasCallStack => [Char] -> a
error [Char]
"doubleFactorialNaive: input should be nonnegative"
  | a
n forall a. Eq a => a -> a -> Bool
== a
0    = Integer
1
  | forall a. Integral a => a -> Bool
odd a
n     = forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [Integer
1,Integer
3..forall a b. (Integral a, Num b) => a -> b
fromIntegral a
n]
  | Bool
otherwise = forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [Integer
2,Integer
4..forall a b. (Integral a, Num b) => a -> b
fromIntegral a
n]

--------------------------------------------------------------------------------
-- * Binomial and multinomial

-- | Binomial numbers (A007318). Note: This is zero for @n<0@ or @k<0@; see also 'signedBinomial' below.
binomial :: Integral a => a -> a -> Integer
binomial :: forall a. Integral a => a -> a -> Integer
binomial = forall a. Integral a => a -> a -> Integer
binomialSplit

-- | Faster implementation of binomial
binomialSplit :: Integral a => a -> a -> Integer
binomialSplit :: forall a. Integral a => a -> a -> Integer
binomialSplit a
n a
k 
  | a
k forall a. Ord a => a -> a -> Bool
> a
n = Integer
0
  | a
k forall a. Ord a => a -> a -> Bool
< a
0 = Integer
0
  | a
k forall a. Ord a => a -> a -> Bool
> (a
n forall a. Integral a => a -> a -> a
`div` a
2) = forall a. Integral a => a -> a -> Integer
binomialSplit a
n (a
nforall a. Num a => a -> a -> a
-a
k)
  | Bool
otherwise = (forall a. Integral a => a -> a -> Integer
productFromTo (a
nforall a. Num a => a -> a -> a
-a
k) a
n) forall a. Integral a => a -> a -> a
`div` (forall a. Integral a => a -> a -> Integer
productFromTo a
1 a
k)

-- | A007318. Note: This is zero for @n<0@ or @k<0@; see also 'signedBinomial' below.
binomialNaive :: Integral a => a -> a -> Integer
binomialNaive :: forall a. Integral a => a -> a -> Integer
binomialNaive a
n a
k 
  | a
k forall a. Ord a => a -> a -> Bool
> a
n = Integer
0
  | a
k forall a. Ord a => a -> a -> Bool
< a
0 = Integer
0
  | a
k forall a. Ord a => a -> a -> Bool
> (a
n forall a. Integral a => a -> a -> a
`div` a
2) = forall a. Integral a => a -> a -> Integer
binomial a
n (a
nforall a. Num a => a -> a -> a
-a
k)
  | Bool
otherwise = (forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [Integer
n'forall a. Num a => a -> a -> a
-Integer
k'forall a. Num a => a -> a -> a
+Integer
1 .. Integer
n']) forall a. Integral a => a -> a -> a
`div` (forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [Integer
1..Integer
k'])
  where 
    k' :: Integer
k' = forall a b. (Integral a, Num b) => a -> b
fromIntegral a
k
    n' :: Integer
n' = 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 forall a. Ord a => a -> a -> Bool
>= Int
0     = forall a. Integral a => a -> a -> Integer
binomial Int
n Int
k
  | Int
k forall a. Ord a => a -> a -> Bool
>= Int
0     = forall a b. (Integral a, Num b) => a -> b -> b
negateIfOdd    Int
k  forall a b. (a -> b) -> a -> b
$ forall a. Integral a => a -> a -> Integer
binomial (Int
kforall a. Num a => a -> a -> a
-Int
nforall a. Num a => a -> a -> a
-Int
1)   Int
k  
  | Bool
otherwise  = forall a b. (Integral a, Num b) => a -> b -> b
negateIfOdd (Int
nforall a. Num a => a -> a -> a
+Int
k) forall a b. (a -> b) -> a -> b
$ forall a. Integral a => a -> a -> Integer
binomial (-Int
kforall a. Num a => a -> a -> a
-Int
1) (-Int
nforall 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 :: forall a. Integral a => a -> [Integer]
pascalRow a
n' = Integer -> Integer -> [Integer]
worker Integer
0 Integer
1 where
  n :: Integer
n = forall a b. (Integral a, Num b) => a -> b
fromIntegral a
n'
  worker :: Integer -> Integer -> [Integer]
worker Integer
j Integer
x
    | Integer
jforall a. Ord a => a -> a -> Bool
>Integer
n   = [] 
    | Bool
True  = let j' :: Integer
j'=Integer
jforall a. Num a => a -> a -> a
+Integer
1 in Integer
x forall a. a -> [a] -> [a]
: Integer -> Integer -> [Integer]
worker Integer
j' (forall a. Integral a => a -> a -> a
div (Integer
xforall a. Num a => a -> a -> a
*(Integer
nforall a. Num a => a -> a -> a
-Integer
j)) Integer
j') 

multinomial :: Integral a => [a] -> Integer
multinomial :: forall a. Integral a => [a] -> Integer
multinomial [a]
xs = forall a. Integral a => a -> a -> a
div
  (forall a. Integral a => a -> Integer
factorial (forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [a]
xs))
  (forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [ 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 :: forall a. Integral a => a -> Integer
catalan a
n 
  | a
n forall a. Ord a => a -> a -> Bool
< a
0     = Integer
0
  | Bool
otherwise = forall a. Integral a => a -> a -> Integer
binomial (a
nforall a. Num a => a -> a -> a
+a
n) a
n forall a. Integral a => a -> a -> a
`div` forall a b. (Integral a, Num b) => a -> b
fromIntegral (a
nforall 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 :: forall a. Integral a => a -> a -> Integer
catalanTriangle a
n a
k
  | a
k forall a. Ord a => a -> a -> Bool
> a
n     = Integer
0
  | a
k forall a. Ord a => a -> a -> Bool
< a
0     = Integer
0
  | Bool
otherwise = (forall a. Integral a => a -> a -> Integer
binomial (a
nforall a. Num a => a -> a -> a
+a
k) a
n forall a. Num a => a -> a -> a
* forall a b. (Integral a, Num b) => a -> b
fromIntegral (a
nforall a. Num a => a -> a -> a
-a
kforall a. Num a => a -> a -> a
+a
1)) forall a. Integral a => a -> a -> a
`div` forall a b. (Integral a, Num b) => a -> b
fromIntegral (a
nforall 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 :: forall a. Integral a => a -> Array Int Integer
signedStirling1stArray a
n
  | a
n forall a. Ord a => a -> a -> Bool
<  a
1    = forall a. HasCallStack => [Char] -> a
error [Char]
"stirling1stArray: n should be at least 1"
  | a
n forall a. Eq a => a -> a -> Bool
== a
1    = forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (Int
1,Int
1 ) [Integer
1]
  | Bool
otherwise = forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (Int
1,Int
n') [ Int -> Integer
lkp (Int
kforall a. Num a => a -> a -> a
-Int
1) forall a. Num a => a -> a -> a
- forall a b. (Integral a, Num b) => a -> b
fromIntegral (a
nforall a. Num a => a -> a -> a
-a
1) forall a. Num a => a -> a -> a
* Int -> Integer
lkp Int
k | Int
k<-[Int
1..Int
n'] ] 
  where
    prev :: Array Int Integer
prev = forall a. Integral a => a -> Array Int Integer
signedStirling1stArray (a
nforall a. Num a => a -> a -> a
-a
1)
    n' :: Int
n' = forall a b. (Integral a, Num b) => a -> b
fromIntegral a
n :: Int
    lkp :: Int -> Integer
lkp Int
j | Int
j forall a. Ord a => a -> a -> Bool
<  Int
1    = Integer
0
          | Int
j forall a. Ord a => a -> a -> Bool
>= Int
n'   = Integer
0
          | Bool
otherwise = Array Int Integer
prev 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 :: forall a. Integral a => a -> a -> Integer
signedStirling1st a
n a
k 
  | a
kforall a. Eq a => a -> a -> Bool
==a
0 Bool -> Bool -> Bool
&& a
nforall a. Eq a => a -> a -> Bool
==a
0 = Integer
1
  | a
k forall a. Ord a => a -> a -> Bool
< a
1        = Integer
0
  | a
k forall a. Ord a => a -> a -> Bool
> a
n        = Integer
0
  | Bool
otherwise    = forall a. Integral a => a -> Array Int Integer
signedStirling1stArray a
n forall i e. Ix i => Array i e -> i -> e
! (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 :: forall a. Integral a => a -> a -> Integer
unsignedStirling1st a
n a
k = forall a. Num a => a -> a
abs (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 :: forall a. Integral a => a -> a -> Integer
stirling2nd a
n a
k 
  | a
kforall a. Eq a => a -> a -> Bool
==a
0 Bool -> Bool -> Bool
&& a
nforall a. Eq a => a -> a -> Bool
==a
0 = Integer
1
  | a
k forall a. Ord a => a -> a -> Bool
< a
1        = Integer
0
  | a
k forall a. Ord a => a -> a -> Bool
> a
n        = Integer
0
  | Bool
otherwise = forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Integer]
xs forall a. Integral a => a -> a -> a
`div` forall a. Integral a => a -> Integer
factorial a
k where
      xs :: [Integer]
xs = [ forall a b. (Integral a, Num b) => a -> b -> b
negateIfOdd (a
kforall a. Num a => a -> a -> a
-a
i) forall a b. (a -> b) -> a -> b
$ forall a. Integral a => a -> a -> Integer
binomial a
k a
i forall a. Num a => a -> a -> a
* (forall a b. (Integral a, Num b) => a -> b
fromIntegral a
i)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 :: forall a. Integral a => a -> Rational
bernoulli a
n 
  | a
n forall a. Ord a => a -> a -> Bool
<  a
0    = forall a. HasCallStack => [Char] -> a
error [Char]
"bernoulli: n should be nonnegative"
  | a
n forall a. Eq a => a -> a -> Bool
== a
0    = Rational
1
  | a
n forall a. Eq a => a -> a -> Bool
== a
1    = -Rational
1forall a. Fractional a => a -> a -> a
/Rational
2
  | Bool
otherwise = 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 = forall a. Real a => a -> Rational
toRational (forall a b. (Integral a, Num b) => a -> b -> b
negateIfOdd (a
nforall a. Num a => a -> a -> a
+a
k) forall a b. (a -> b) -> a -> b
$ forall a. Integral a => a -> Integer
factorial a
k forall a. Num a => a -> a -> a
* forall a. Integral a => a -> a -> Integer
stirling2nd a
n a
k) 
        forall a. Fractional a => a -> a -> a
/ forall a. Real a => a -> Rational
toRational (a
kforall 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 :: forall a. Integral a => a -> Array Int Integer
bellNumbersArray a
nn = Array Int Integer
arr where
  arr :: Array Int Integer
arr = forall i e. Ix i => (i, i) -> [(i, e)] -> Array i e
array (Int
0::Int,Int
n) [(Int, Integer)]
kvs 
  n :: Int
n = forall a b. (Integral a, Num b) => a -> b
fromIntegral a
nn :: Int
  kvs :: [(Int, Integer)]
kvs = (Int
0,Integer
1) forall a. a -> [a] -> [a]
: [ (Int
k, Int -> Integer
f Int
k) | Int
k<-[Int
1..Int
n] ] 
  f :: Int -> Integer
f Int
n = forall a. Num a => [a] -> a
sum' [ forall a. Integral a => a -> a -> Integer
binomial (Int
nforall a. Num a => a -> a -> a
-Int
1) Int
k forall a. Num a => a -> a -> a
* Array Int Integer
arr forall i e. Ix i => Array i e -> i -> e
! Int
k | Int
k<-[Int
0..Int
nforall 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 :: forall a. Integral a => a -> Integer
bellNumber a
nn
  | Int
n forall a. Ord a => a -> a -> Bool
<  Int
0     = forall a. HasCallStack => [Char] -> a
error [Char]
"bellNumber: expecting a nonnegative index"
  | Int
n forall a. Eq a => a -> a -> Bool
== Int
0     = Integer
1
  | Bool
otherwise  = forall a. Num a => [a] -> a
sum' [ forall a. Integral a => a -> a -> Integer
stirling2nd Int
n Int
k | Int
k<-[Int
1..Int
n] ] 
  where
    n :: Int
n = forall a b. (Integral a, Num b) => a -> b
fromIntegral a
nn :: Int

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