-- | 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 :: a -> Integer
factorial = a -> Integer
forall a. Integral a => a -> Integer
factorialSplit

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

-- | Naive implementation of factorial
factorialNaive :: Integral a => a -> Integer
factorialNaive :: a -> Integer
factorialNaive 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]
"factorialNaive: 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]

-- | \"Swing factorial\" algorithm
factorialSwing :: Integral a => a -> Integer
factorialSwing :: a -> Integer
factorialSwing a
n = [(Integer, Int)] -> Integer
productOfFactors (Int -> [(Integer, Int)]
factorialPrimeExponents (Int -> [(Integer, Int)]) -> Int -> [(Integer, Int)]
forall a b. (a -> b) -> a -> b
$ a -> Int
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 = ((Integer, Int) -> Bool) -> [(Integer, Int)] -> [(Integer, Int)]
forall a. (a -> Bool) -> [a] -> [a]
filter (Integer, Int) -> Bool
forall a a. (Ord a, Num a) => (a, a) -> Bool
cond ([(Integer, Int)] -> [(Integer, Int)])
-> [(Integer, Int)] -> [(Integer, Int)]
forall a b. (a -> b) -> a -> b
$ [Integer] -> [Int] -> [(Integer, Int)]
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 a -> a -> Bool
forall a. Ord a => a -> a -> Bool
> a
0

factorialPrimeExponentsNaive :: forall a. Integral a => a -> [(Integer,Int)]
factorialPrimeExponentsNaive :: a -> [(Integer, Int)]
factorialPrimeExponentsNaive a
n = [(Integer, Int)]
result where
  fi :: a -> Integer
fi = a -> Integer
forall a b. (Integral a, Num b) => a -> b
fromIntegral :: a -> Integer
  result :: [(Integer, Int)]
result = Map Integer Int -> [(Integer, Int)]
forall k a. Map k a -> [(k, a)]
Map.toList 
         (Map Integer Int -> [(Integer, Int)])
-> Map Integer Int -> [(Integer, Int)]
forall a b. (a -> b) -> a -> b
$ (Int -> Int -> Int) -> [Map Integer Int] -> Map Integer Int
forall (f :: * -> *) k a.
(Foldable f, Ord k) =>
(a -> a -> a) -> f (Map k a) -> Map k a
Map.unionsWith Int -> Int -> Int
forall a. Num a => a -> a -> a
(+) 
         ([Map Integer Int] -> Map Integer Int)
-> [Map Integer Int] -> Map Integer Int
forall a b. (a -> b) -> a -> b
$ ([(Integer, Int)] -> Map Integer Int)
-> [[(Integer, Int)]] -> [Map Integer Int]
forall a b. (a -> b) -> [a] -> [b]
map [(Integer, Int)] -> Map Integer Int
forall k a. Ord k => [(k, a)] -> Map k a
Map.fromList 
         ([[(Integer, Int)]] -> [Map Integer Int])
-> [[(Integer, Int)]] -> [Map Integer Int]
forall a b. (a -> b) -> a -> b
$ (Integer -> [(Integer, Int)]) -> [Integer] -> [[(Integer, Int)]]
forall a b. (a -> b) -> [a] -> [b]
map Integer -> [(Integer, Int)]
factorize 
         ([Integer] -> [[(Integer, Int)]])
-> [Integer] -> [[(Integer, Int)]]
forall a b. (a -> b) -> a -> b
$ (a -> Integer) -> [a] -> [Integer]
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  = (Int -> Int) -> [Int] -> [Int]
forall a b. (a -> b) -> [a] -> [b]
map ((Int -> Int -> Int) -> Int -> Int -> Int
forall a b c. (a -> b -> c) -> b -> a -> c
flip Int -> Int -> Int
forall a. Bits a => a -> Int -> a
shiftL Int
1) ([Int] -> [Int]) -> [Int] -> [Int]
forall a b. (a -> b) -> a -> b
$ Int -> [Int]
go (Int -> Int -> Int
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
xInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
y) Int -> [Int] -> [Int]
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 Int -> [Int] -> [Int]
forall a. a -> [a] -> [a]
: (Integer -> Int) -> [Integer] -> [Int]
forall a b. (a -> b) -> [a] -> [b]
map Integer -> Int
expo ([Integer] -> [Integer]
forall a. [a] -> [a]
tail [Integer]
ps) where

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

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

    expo2 :: Int
    expo2 :: Int
expo2 = Int -> Int -> Int
go Int
0 (Int -> Int -> Int
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 Int -> Int -> Bool
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 Int -> Int -> Int
forall a. Num a => a -> a -> a
+ (Int
r Int -> Int -> Int
forall a. Bits a => a -> a -> a
.&. Int
1)
          r' :: Int
r'   = Int -> Int -> Int
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 (Int -> Int -> Int
forall a. Integral a => a -> a -> a
div Int
n Int
p) where
      p :: Int
p = Integer -> Int
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 Int -> Int -> Bool
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 Int -> Int -> Int
forall a. Num a => a -> a -> a
+ (Int
r Int -> Int -> Int
forall a. Bits a => a -> a -> a
.&. Int
1)
          r' :: Int
r'   = Int -> Int -> Int
forall a. Integral a => a -> a -> a
div Int
r Int
p

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

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

-- | Faster implementation of the double factorial function
doubleFactorialSplit :: Integral a => a -> Integer
doubleFactorialSplit :: a -> Integer
doubleFactorialSplit 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]
"doubleFactorialSplit: 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     = a -> a -> Integer
forall a. Integral a => a -> a -> Integer
productFromToStride2 a
2 a
n
  | Bool
otherwise = let halfn :: a
halfn = a -> a -> a
forall a. Integral a => a -> a -> a
div a
n a
2 
                in  Integer -> Int -> Integer
forall a. Bits a => a -> Int -> a
shiftL (a -> Integer
forall a. Integral a => a -> Integer
factorialSplit a
halfn) (a -> Int
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 :: a -> Integer
doubleFactorialNaive 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]
"doubleFactorialNaive: 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]

--------------------------------------------------------------------------------
-- * 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 :: a -> a -> Integer
binomial = a -> a -> Integer
forall a. Integral a => a -> a -> Integer
binomialSplit

-- | Faster implementation of binomial
binomialSplit :: Integral a => a -> a -> Integer
binomialSplit :: a -> a -> Integer
binomialSplit 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
binomialSplit a
n (a
na -> a -> a
forall a. Num a => a -> a -> a
-a
k)
  | Bool
otherwise = (a -> a -> Integer
forall a. Integral a => a -> a -> Integer
productFromTo (a
na -> a -> a
forall a. Num a => a -> a -> a
-a
k) a
n) Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`div` (a -> a -> Integer
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 :: a -> a -> Integer
binomialNaive 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

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