    -- |
    -- Module:      Math.NumberTheory.Recurrences.Bilinear
    -- Copyright:   (c) 2016 Andrew Lelechenko
    -- Licence:     MIT
    -- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>
    --
    -- Bilinear recurrent sequences and Bernoulli numbers,
    -- roughly covering Ch. 5-6 of /Concrete Mathematics/
    -- by R. L. Graham, D. E. Knuth and O. Patashnik.
    --
    -- #memory# __Note on memory leaks and memoization.__
    -- Top-level definitions in this module are polymorphic, so the results of computations are not retained in memory.
    -- Make them monomorphic to take advantages of memoization. Compare
    --
    -- >>> binomial !! 1000 !! 1000 :: Integer -- (0.01 secs, 1,385,512 bytes)
    -- 1
    -- >>> binomial !! 1000 !! 1000 :: Integer -- (0.01 secs, 1,381,616 bytes)
    -- 1
    --
    -- against
    --
    -- >>> let binomial' = binomial :: [[Integer]]
    -- >>> binomial' !! 1000 !! 1000 :: Integer -- (0.01 secs, 1,381,696 bytes)
    -- 1
    -- >>> binomial' !! 1000 !! 1000 :: Integer -- (0.01 secs, 391,152 bytes)
    -- 1

    {-# LANGUAGE BangPatterns        #-}
    {-# LANGUAGE ScopedTypeVariables #-}

    module Math.NumberTheory.Recurrences.Bilinear
      ( -- * Pascal triangle
        binomial
      , binomialRotated
      , binomialLine
      , binomialDiagonal
      , binomialFactors
        -- * Other recurrences
      , stirling1
      , stirling2
      , lah
      , eulerian1
      , eulerian2
      , bernoulli
      , euler
      , eulerPolyAt1
      , faulhaberPoly
      ) where

    import Data.Euclidean (GcdDomain(..))
    import Data.List (scanl', zipWith4)
    import Data.Maybe
    import Data.Ratio
    import Data.Semiring (Semiring(..))
    import Numeric.Natural

    import Math.NumberTheory.Recurrences.Linear (factorial)
    import Math.NumberTheory.Primes

    -- | Infinite zero-based table of binomial coefficients (also known as Pascal triangle).
    --
    -- > binomial !! n !! k == n! / k! / (n - k)!
    --
    -- Note that 'binomial' !! n !! k is asymptotically slower
    -- than 'binomialLine' n !! k,
    -- but imposes only 'Semiring' constraint.
    --
    -- >>> take 6 binomial :: [[Int]]
    -- [[1],[1,1],[1,2,1],[1,3,3,1],[1,4,6,4,1],[1,5,10,10,5,1]]
    binomial :: Semiring a => [[a]]
    binomial :: forall a. Semiring a => [[a]]
    binomial = ([a] -> [a]) -> [a] -> [[a]]
    forall a. (a -> a) -> a -> [a]
    iterate (\[a]
    l -> (a -> a -> a) -> [a] -> [a] -> [a]
    forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
    zipWith a -> a -> a
    forall a. Semiring a => a -> a -> a
    plus ([a]
    l [a] -> [a] -> [a]
    forall a. [a] -> [a] -> [a]
    ++ [a
    forall a. Semiring a => a
    zero]) (a
    forall a. Semiring a => a
    zero a -> [a] -> [a]
    forall a. a -> [a] -> [a]
    : [a]
    l)) [a
    forall a. Semiring a => a
    one]
    {-# SPECIALIZE binomial :: [[Int]]     #-}
    {-# SPECIALIZE binomial :: [[Word]]    #-}
    {-# SPECIALIZE binomial :: [[Integer]] #-}
    {-# SPECIALIZE binomial :: [[Natural]] #-}

    -- | Pascal triangle, rotated by 45 degrees.
    --
    -- > binomialRotated !! n !! k == (n + k)! / n! / k! == binomial !! (n + k) !! k
    --
    -- Note that 'binomialRotated' !! n !! k is asymptotically slower
    -- than 'binomialDiagonal' n !! k,
    -- but imposes only 'Semiring' constraint.
    --
    -- >>> take 6 (map (take 6) binomialRotated) :: [[Int]]
    -- [[1,1,1,1,1,1],[1,2,3,4,5,6],[1,3,6,10,15,21],[1,4,10,20,35,56],[1,5,15,35,70,126],[1,6,21,56,126,252]]
    binomialRotated :: Semiring a => [[a]]
    binomialRotated :: forall a. Semiring a => [[a]]
    binomialRotated = ([a] -> [a]) -> [a] -> [[a]]
    forall a. (a -> a) -> a -> [a]
    iterate ([a] -> [a]
    forall a. HasCallStack => [a] -> [a]
    tail ([a] -> [a]) -> ([a] -> [a]) -> [a] -> [a]
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . (a -> a -> a) -> a -> [a] -> [a]
    forall b a. (b -> a -> b) -> b -> [a] -> [b]
    scanl' a -> a -> a
    forall a. Semiring a => a -> a -> a
    plus a
    forall a. Semiring a => a
    zero) (a -> [a]
    forall a. a -> [a]
    repeat a
    forall a. Semiring a => a
    one)
    {-# SPECIALIZE binomialRotated :: [[Int]]     #-}
    {-# SPECIALIZE binomialRotated :: [[Word]]    #-}
    {-# SPECIALIZE binomialRotated :: [[Integer]] #-}
    {-# SPECIALIZE binomialRotated :: [[Natural]] #-}

    -- | The n-th (zero-based) line of 'binomial'
    -- (and the n-th diagonal of 'binomialRotated').
    --
    -- >>> binomialLine 5
    -- [1,5,10,10,5,1]
    binomialLine :: (Enum a, GcdDomain a) => a -> [a]
    binomialLine :: forall a. (Enum a, GcdDomain a) => a -> [a]
    binomialLine a
    n = (a -> (a, a) -> a) -> a -> [(a, a)] -> [a]
    forall b a. (b -> a -> b) -> b -> [a] -> [b]
    scanl'
      (\a
    x (a
    k, a
    nk1) -> Maybe a -> a
    forall a. HasCallStack => Maybe a -> a
    fromJust (Maybe a -> a) -> Maybe a -> a
    forall a b. (a -> b) -> a -> b
    $ (a
    x a -> a -> a
    forall a. Semiring a => a -> a -> a
    `times` a
    nk1) a -> a -> Maybe a
    forall a. GcdDomain a => a -> a -> Maybe a
    `divide` a
    k)
      a
    forall a. Semiring a => a
    one
      ([a] -> [a] -> [(a, a)]
    forall a b. [a] -> [b] -> [(a, b)]
    zip [a
    forall a. Semiring a => a
    one..a
    n] [a
    n, a -> a
    forall a. Enum a => a -> a
    pred a
    n..a
    forall a. Semiring a => a
    one])
    {-# SPECIALIZE binomialLine :: Int     -> [Int]     #-}
    {-# SPECIALIZE binomialLine :: Word    -> [Word]    #-}
    {-# SPECIALIZE binomialLine :: Integer -> [Integer] #-}
    {-# SPECIALIZE binomialLine :: Natural -> [Natural] #-}

    -- | The n-th (zero-based) diagonal of 'binomial'
    -- (and the n-th line of 'binomialRotated').
    --
    -- >>> take 6 (binomialDiagonal 5)
    -- [1,6,21,56,126,252]
    binomialDiagonal :: (Enum a, GcdDomain a) => a -> [a]
    binomialDiagonal :: forall a. (Enum a, GcdDomain a) => a -> [a]
    binomialDiagonal a
    n = (a -> a -> a) -> a -> [a] -> [a]
    forall b a. (b -> a -> b) -> b -> [a] -> [b]
    scanl'
      (\a
    x a
    k -> Maybe a -> a
    forall a. HasCallStack => Maybe a -> a
    fromJust (a
    x a -> a -> a
    forall a. Semiring a => a -> a -> a
    `times` (a
    n a -> a -> a
    forall a. Semiring a => a -> a -> a
    `plus` a
    k) a -> a -> Maybe a
    forall a. GcdDomain a => a -> a -> Maybe a
    `divide` a
    k))
      a
    forall a. Semiring a => a
    one
      [a
    forall a. Semiring a => a
    one..]
    {-# SPECIALIZE binomialDiagonal :: Int     -> [Int]     #-}
    {-# SPECIALIZE binomialDiagonal :: Word    -> [Word]    #-}
    {-# SPECIALIZE binomialDiagonal :: Integer -> [Integer] #-}
    {-# SPECIALIZE binomialDiagonal :: Natural -> [Natural] #-}

    -- | Prime factors of a binomial coefficient.
    --
    -- > binomialFactors n k == factorise (binomial !! n !! k)
    --
    -- >>> binomialFactors 10 4
    -- [(Prime 2,1),(Prime 3,1),(Prime 5,1),(Prime 7,1)]
    binomialFactors :: Word -> Word -> [(Prime Word, Word)]
    binomialFactors :: Word -> Word -> [(Prime Word, Word)]
    binomialFactors Word
    n Word
    k
      | Word
    n Word -> Word -> Bool
    forall a. Ord a => a -> a -> Bool
    < Word
    2
      = []
      | Bool
    otherwise
      = ((Prime Word, Word) -> Bool)
    -> [(Prime Word, Word)] -> [(Prime Word, Word)]
    forall a. (a -> Bool) -> [a] -> [a]
    filter ((Word -> Word -> Bool
    forall a. Eq a => a -> a -> Bool
    /= Word
    0) (Word -> Bool)
    -> ((Prime Word, Word) -> Word) -> (Prime Word, Word) -> Bool
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . (Prime Word, Word) -> Word
    forall a b. (a, b) -> b
    snd)
      ([(Prime Word, Word)] -> [(Prime Word, Word)])
    -> [(Prime Word, Word)] -> [(Prime Word, Word)]
    forall a b. (a -> b) -> a -> b
    $ (Prime Word -> (Prime Word, Word))
    -> [Prime Word] -> [(Prime Word, Word)]
    forall a b. (a -> b) -> [a] -> [b]
    map (\Prime Word
    p -> (Prime Word
    p, Word -> Word -> Word
    mult (Prime Word -> Word
    forall a. Prime a -> a
    unPrime Prime Word
    p) Word
    n Word -> Word -> Word
    forall a. Num a => a -> a -> a
    - Word -> Word -> Word
    mult (Prime Word -> Word
    forall a. Prime a -> a
    unPrime Prime Word
    p) (Word
    n Word -> Word -> Word
    forall a. Num a => a -> a -> a
    - Word
    k) Word -> Word -> Word
    forall a. Num a => a -> a -> a
    - Word -> Word -> Word
    mult (Prime Word -> Word
    forall a. Prime a -> a
    unPrime Prime Word
    p) Word
    k))
        [Prime Word
    forall a. Bounded a => a
    minBound .. Word -> Prime Word
    forall a.
    (Bits a, Integral a, UniqueFactorisation a) =>
    a -> Prime a
    precPrime Word
    n]
      where
        mult :: Word -> Word -> Word
        mult :: Word -> Word -> Word
    mult Word
    p Word
    m = Word -> Word -> Word
    go Word
    mp Word
    mp
          where
            mp :: Word
    mp = Word
    m Word -> Word -> Word
    forall a. Integral a => a -> a -> a
    `quot` Word
    p
            go :: Word -> Word -> Word
    go !Word
    acc !Word
    x
              | Word
    x Word -> Word -> Bool
    forall a. Ord a => a -> a -> Bool
    >= Word
    p = let xp :: Word
    xp = Word
    x Word -> Word -> Word
    forall a. Integral a => a -> a -> a
    `quot` Word
    p in Word -> Word -> Word
    go (Word
    acc Word -> Word -> Word
    forall a. Num a => a -> a -> a
    + Word
    xp) Word
    xp
              | Bool
    otherwise = Word
    acc

    -- | Infinite zero-based table of <https://en.wikipedia.org/wiki/Stirling_numbers_of_the_first_kind Stirling numbers of the first kind>.
    --
    -- >>> take 5 (map (take 5) stirling1)
    -- [[1],[0,1],[0,1,1],[0,2,3,1],[0,6,11,6,1]]
    --
    -- Complexity: @stirling1 !! n !! k@ is O(n ln n) bits long, its computation
    -- takes O(k n^2 ln n) time and forces thunks @stirling1 !! i !! j@ for @0 <= i <= n@ and @max(0, k - n + i) <= j <= k@.
    --
    -- One could also consider 'Math.Combinat.Numbers.unsignedStirling1st' from <http://hackage.haskell.org/package/combinat combinat> package to compute stand-alone values.
    stirling1 :: (Num a, Enum a) => [[a]]
    stirling1 :: forall a. (Num a, Enum a) => [[a]]
    stirling1 = ([a] -> a -> [a]) -> [a] -> [a] -> [[a]]
    forall b a. (b -> a -> b) -> b -> [a] -> [b]
    scanl [a] -> a -> [a]
    forall {a}. (Num a, Enum a) => [a] -> a -> [a]
    f [a
    1] [a
    0..]
      where
        f :: [a] -> a -> [a]
    f [a]
    xs a
    n = a
    0 a -> [a] -> [a]
    forall a. a -> [a] -> [a]
    : (a -> a -> a -> a) -> a -> [a] -> a -> [a]
    forall b a. Enum b => (b -> a -> a -> b) -> b -> [a] -> a -> [b]
    zipIndexedListWithTail (\a
    _ a
    x a
    y -> a
    x a -> a -> a
    forall a. Num a => a -> a -> a
    + a
    n a -> a -> a
    forall a. Num a => a -> a -> a
    * a
    y) a
    1 [a]
    xs a
    0
    {-# SPECIALIZE stirling1 :: [[Int]]     #-}
    {-# SPECIALIZE stirling1 :: [[Word]]    #-}
    {-# SPECIALIZE stirling1 :: [[Integer]] #-}
    {-# SPECIALIZE stirling1 :: [[Natural]] #-}

    -- | Infinite zero-based table of <https://en.wikipedia.org/wiki/Stirling_numbers_of_the_second_kind Stirling numbers of the second kind>.
    --
    -- >>> take 5 (map (take 5) stirling2)
    -- [[1],[0,1],[0,1,1],[0,1,3,1],[0,1,7,6,1]]
    --
    -- Complexity: @stirling2 !! n !! k@ is O(n ln n) bits long, its computation
    -- takes O(k n^2 ln n) time and forces thunks @stirling2 !! i !! j@ for @0 <= i <= n@ and @max(0, k - n + i) <= j <= k@.
    --
    -- One could also consider 'Math.Combinat.Numbers.stirling2nd' from <http://hackage.haskell.org/package/combinat combinat> package to compute stand-alone values.
    stirling2 :: (Num a, Enum a) => [[a]]
    stirling2 :: forall a. (Num a, Enum a) => [[a]]
    stirling2 = ([a] -> [a]) -> [a] -> [[a]]
    forall a. (a -> a) -> a -> [a]
    iterate [a] -> [a]
    forall {a}. (Num a, Enum a) => [a] -> [a]
    f [a
    1]
      where
        f :: [a] -> [a]
    f [a]
    xs = a
    0 a -> [a] -> [a]
    forall a. a -> [a] -> [a]
    : (a -> a -> a -> a) -> a -> [a] -> a -> [a]
    forall b a. Enum b => (b -> a -> a -> b) -> b -> [a] -> a -> [b]
    zipIndexedListWithTail (\a
    k a
    x a
    y -> a
    x a -> a -> a
    forall a. Num a => a -> a -> a
    + a
    k a -> a -> a
    forall a. Num a => a -> a -> a
    * a
    y) a
    1 [a]
    xs a
    0
    {-# SPECIALIZE stirling2 :: [[Int]]     #-}
    {-# SPECIALIZE stirling2 :: [[Word]]    #-}
    {-# SPECIALIZE stirling2 :: [[Integer]] #-}
    {-# SPECIALIZE stirling2 :: [[Natural]] #-}

    -- | Infinite one-based table of <https://en.wikipedia.org/wiki/Lah_number Lah numbers>.
    -- @lah !! n !! k@ equals to lah(n + 1, k + 1).
    --
    -- >>> take 5 (map (take 5) lah)
    -- [[1],[2,1],[6,6,1],[24,36,12,1],[120,240,120,20,1]]
    --
    -- Complexity: @lah !! n !! k@ is O(n ln n) bits long, its computation
    -- takes O(k n ln n) time and forces thunks @lah !! n !! i@ for @0 <= i <= k@.
    lah :: Integral a => [[a]]
    -- Implementation was derived from code by https://github.com/grandpascorpion
    lah :: forall a. Integral a => [[a]]
    lah = (a -> a -> [a]) -> [a] -> [a] -> [[a]]
    forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
    zipWith a -> a -> [a]
    forall {b}. Integral b => b -> b -> [b]
    f ([a] -> [a]
    forall a. HasCallStack => [a] -> [a]
    tail [a]
    forall a. (Num a, Enum a) => [a]
    factorial) [a
    1..]
      where
        f :: b -> b -> [b]
    f b
    nf b
    n = (b -> b -> b) -> b -> [b] -> [b]
    forall b a. (b -> a -> b) -> b -> [a] -> [b]
    scanl (\b
    x b
    k -> b
    x b -> b -> b
    forall a. Num a => a -> a -> a
    * (b
    n b -> b -> b
    forall a. Num a => a -> a -> a
    - b
    k) b -> b -> b
    forall a. Integral a => a -> a -> a
    `div` (b
    k b -> b -> b
    forall a. Num a => a -> a -> a
    * (b
    k b -> b -> b
    forall a. Num a => a -> a -> a
    + b
    1))) b
    nf [b
    1..b
    nb -> b -> b
    forall a. Num a => a -> a -> a
    -b
    1]
    {-# SPECIALIZE lah :: [[Int]]     #-}
    {-# SPECIALIZE lah :: [[Word]]    #-}
    {-# SPECIALIZE lah :: [[Integer]] #-}
    {-# SPECIALIZE lah :: [[Natural]] #-}

    -- | Infinite zero-based table of <https://en.wikipedia.org/wiki/Eulerian_number Eulerian numbers of the first kind>.
    --
    -- >>> take 5 (map (take 5) eulerian1)
    -- [[],[1],[1,1],[1,4,1],[1,11,11,1]]
    --
    -- Complexity: @eulerian1 !! n !! k@ is O(n ln n) bits long, its computation
    -- takes O(k n^2 ln n) time and forces thunks @eulerian1 !! i !! j@ for @0 <= i <= n@ and @max(0, k - n + i) <= j <= k@.
    --
    eulerian1 :: (Num a, Enum a) => [[a]]
    eulerian1 :: forall a. (Num a, Enum a) => [[a]]
    eulerian1 = ([a] -> a -> [a]) -> [a] -> [a] -> [[a]]
    forall b a. (b -> a -> b) -> b -> [a] -> [b]
    scanl [a] -> a -> [a]
    forall {a}. (Num a, Enum a) => [a] -> a -> [a]
    f [] [a
    1..]
      where
        f :: [a] -> a -> [a]
    f [a]
    xs a
    n = a
    1 a -> [a] -> [a]
    forall a. a -> [a] -> [a]
    : (a -> a -> a -> a) -> a -> [a] -> a -> [a]
    forall b a. Enum b => (b -> a -> a -> b) -> b -> [a] -> a -> [b]
    zipIndexedListWithTail (\a
    k a
    x a
    y -> (a
    n a -> a -> a
    forall a. Num a => a -> a -> a
    - a
    k) a -> a -> a
    forall a. Num a => a -> a -> a
    * a
    x a -> a -> a
    forall a. Num a => a -> a -> a
    + (a
    k a -> a -> a
    forall a. Num a => a -> a -> a
    + a
    1) a -> a -> a
    forall a. Num a => a -> a -> a
    * a
    y) a
    1 [a]
    xs a
    0
    {-# SPECIALIZE eulerian1 :: [[Int]]     #-}
    {-# SPECIALIZE eulerian1 :: [[Word]]    #-}
    {-# SPECIALIZE eulerian1 :: [[Integer]] #-}
    {-# SPECIALIZE eulerian1 :: [[Natural]] #-}

    -- | Infinite zero-based table of <https://en.wikipedia.org/wiki/Eulerian_number#Eulerian_numbers_of_the_second_kind Eulerian numbers of the second kind>.
    --
    -- >>> take 5 (map (take 5) eulerian2)
    -- [[],[1],[1,2],[1,8,6],[1,22,58,24]]
    --
    -- Complexity: @eulerian2 !! n !! k@ is O(n ln n) bits long, its computation
    -- takes O(k n^2 ln n) time and forces thunks @eulerian2 !! i !! j@ for @0 <= i <= n@ and @max(0, k - n + i) <= j <= k@.
    --
    eulerian2 :: (Num a, Enum a) => [[a]]
    eulerian2 :: forall a. (Num a, Enum a) => [[a]]
    eulerian2 = ([a] -> a -> [a]) -> [a] -> [a] -> [[a]]
    forall b a. (b -> a -> b) -> b -> [a] -> [b]
    scanl [a] -> a -> [a]
    forall {a}. (Num a, Enum a) => [a] -> a -> [a]
    f [] [a
    1..]
      where
        f :: [a] -> a -> [a]
    f [a]
    xs a
    n = a
    1 a -> [a] -> [a]
    forall a. a -> [a] -> [a]
    : (a -> a -> a -> a) -> a -> [a] -> a -> [a]
    forall b a. Enum b => (b -> a -> a -> b) -> b -> [a] -> a -> [b]
    zipIndexedListWithTail (\a
    k a
    x a
    y -> (a
    2 a -> a -> a
    forall a. Num a => a -> a -> a
    * a
    n a -> a -> a
    forall a. Num a => a -> a -> a
    - a
    k a -> a -> a
    forall a. Num a => a -> a -> a
    - a
    1) a -> a -> a
    forall a. Num a => a -> a -> a
    * a
    x a -> a -> a
    forall a. Num a => a -> a -> a
    + (a
    k a -> a -> a
    forall a. Num a => a -> a -> a
    + a
    1) a -> a -> a
    forall a. Num a => a -> a -> a
    * a
    y) a
    1 [a]
    xs a
    0
    {-# SPECIALIZE eulerian2 :: [[Int]]     #-}
    {-# SPECIALIZE eulerian2 :: [[Word]]    #-}
    {-# SPECIALIZE eulerian2 :: [[Integer]] #-}
    {-# SPECIALIZE eulerian2 :: [[Natural]] #-}

    -- | Infinite zero-based sequence of <https://en.wikipedia.org/wiki/Bernoulli_number Bernoulli numbers>,
    -- computed via <https://en.wikipedia.org/wiki/Bernoulli_number#Connection_with_Stirling_numbers_of_the_second_kind connection>
    -- with 'stirling2'.
    --
    -- >>> take 5 bernoulli
    -- [1 % 1,(-1) % 2,1 % 6,0 % 1,(-1) % 30]
    --
    -- Complexity: @bernoulli !! n@ is O(n ln n) bits long, its computation
    -- takes O(n^3 ln n) time and forces thunks @stirling2 !! i !! j@ for @0 <= i <= n@ and @0 <= j <= i@.
    --
    -- One could also consider 'Math.Combinat.Numbers.bernoulli' from <http://hackage.haskell.org/package/combinat combinat> package to compute stand-alone values.
    bernoulli :: Integral a => [Ratio a]
    bernoulli :: forall a. Integral a => [Ratio a]
    bernoulli = ([Ratio a] -> [Ratio a]) -> [Ratio a] -> [Ratio a]
    forall a.
    Integral a =>
    ([Ratio a] -> [Ratio a]) -> [Ratio a] -> [Ratio a]
    helperForBEEP [Ratio a] -> [Ratio a]
    forall a. a -> a
    id ((Ratio a -> Ratio a) -> [Ratio a] -> [Ratio a]
    forall a b. (a -> b) -> [a] -> [b]
    map Ratio a -> Ratio a
    forall a. Fractional a => a -> a
    recip [Ratio a
    1..])
    {-# SPECIALIZE bernoulli :: [Ratio Int] #-}
    {-# SPECIALIZE bernoulli :: [Rational] #-}

    -- | <https://en.wikipedia.org/wiki/Faulhaber%27s_formula Faulhaber's formula>.
    --
    -- >>> sum (map (^ 10) [0..100])
    -- 959924142434241924250
    -- >>> sum $ zipWith (*) (faulhaberPoly 10) (iterate (* 100) 1)
    -- 959924142434241924250 % 1
    faulhaberPoly :: (GcdDomain a, Integral a) => Int -> [Ratio a]
    -- Implementation by https://github.com/CarlEdman
    faulhaberPoly :: forall a. (GcdDomain a, Integral a) => Int -> [Ratio a]
    faulhaberPoly Int
    p
      = (Ratio a -> Ratio a -> Ratio a)
    -> [Ratio a] -> [Ratio a] -> [Ratio a]
    forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
    zipWith Ratio a -> Ratio a -> Ratio a
    forall a. Num a => a -> a -> a
    (*) ((Ratio a
    0Ratio a -> [Ratio a] -> [Ratio a]
    forall a. a -> [a] -> [a]
    :)
      ([Ratio a] -> [Ratio a]) -> [Ratio a] -> [Ratio a]
    forall a b. (a -> b) -> a -> b
    $ [Ratio a] -> [Ratio a]
    forall a. [a] -> [a]
    reverse
      ([Ratio a] -> [Ratio a]) -> [Ratio a] -> [Ratio a]
    forall a b. (a -> b) -> a -> b
    $ Int -> [Ratio a] -> [Ratio a]
    forall a. Int -> [a] -> [a]
    take (Int
    p Int -> Int -> Int
    forall a. Num a => a -> a -> a
    + Int
    1) [Ratio a]
    forall a. Integral a => [Ratio a]
    bernoulli)
      ([Ratio a] -> [Ratio a]) -> [Ratio a] -> [Ratio a]
    forall a b. (a -> b) -> a -> b
    $ (a -> Ratio a) -> [a] -> [Ratio a]
    forall a b. (a -> b) -> [a] -> [b]
    map (a -> a -> Ratio a
    forall a. Integral a => a -> a -> Ratio a
    % (Int -> a
    forall a b. (Integral a, Num b) => a -> b
    fromIntegral Int
    pa -> a -> a
    forall a. Num a => a -> a -> a
    +a
    1))
      ([a] -> [Ratio a]) -> [a] -> [Ratio a]
    forall a b. (a -> b) -> a -> b
    $ (a -> a -> a) -> [a] -> [a] -> [a]
    forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
    zipWith a -> a -> a
    forall a. Num a => a -> a -> a
    (*) ((a -> a) -> a -> [a]
    forall a. (a -> a) -> a -> [a]
    iterate a -> a
    forall a. Num a => a -> a
    negate (if Int -> Bool
    forall a. Integral a => a -> Bool
    odd Int
    p then a
    1 else -a
    1))
      ([a] -> [a]) -> [a] -> [a]
    forall a b. (a -> b) -> a -> b
    $ [[a]]
    forall a. Semiring a => [[a]]
    binomial [[a]] -> Int -> [a]
    forall a. HasCallStack => [a] -> Int -> a
    !! (Int
    pInt -> Int -> Int
    forall a. Num a => a -> a -> a
    +Int
    1)

    -- | Infinite zero-based list of <https://en.wikipedia.org/wiki/Euler_number Euler numbers>.
    -- The algorithm used was derived from <http://www.emis.ams.org/journals/JIS/VOL4/CHEN/AlgBE2.pdf Algorithms for Bernoulli numbers and Euler numbers>
    -- by Kwang-Wu Chen, second formula of the Corollary in page 7.
    -- Sequence <https://oeis.org/A122045 A122045> in OEIS.
    --
    -- >>> take 10 euler' :: [Rational]
    -- [1 % 1,0 % 1,(-1) % 1,0 % 1,5 % 1,0 % 1,(-61) % 1,0 % 1,1385 % 1,0 % 1]
    euler' :: forall a . Integral a => [Ratio a]
    euler' :: forall a. Integral a => [Ratio a]
    euler' = [Ratio a] -> [Ratio a]
    forall a. HasCallStack => [a] -> [a]
    tail ([Ratio a] -> [Ratio a]) -> [Ratio a] -> [Ratio a]
    forall a b. (a -> b) -> a -> b
    $ ([Ratio a] -> [Ratio a]) -> [Ratio a] -> [Ratio a]
    forall a.
    Integral a =>
    ([Ratio a] -> [Ratio a]) -> [Ratio a] -> [Ratio a]
    helperForBEEP [Ratio a] -> [Ratio a]
    forall a. HasCallStack => [a] -> [a]
    tail [Ratio a]
    as
      where
        as :: [Ratio a]
        as :: [Ratio a]
    as = (a -> a -> a -> Ratio a) -> [a] -> [a] -> [a] -> [Ratio a]
    forall a b c d. (a -> b -> c -> d) -> [a] -> [b] -> [c] -> [d]
    zipWith3
            (\a
    sgn a
    frac a
    ones -> (a
    sgn a -> a -> a
    forall a. Num a => a -> a -> a
    * a
    ones) a -> a -> Ratio a
    forall a. Integral a => a -> a -> Ratio a
    % a
    frac)
            ([a] -> [a]
    forall a. HasCallStack => [a] -> [a]
    cycle [a
    1, a
    1, a
    1, a
    1, -a
    1, -a
    1, -a
    1, -a
    1])
            ([a] -> [a]
    forall a. [a] -> [a]
    dups ((a -> a) -> a -> [a]
    forall a. (a -> a) -> a -> [a]
    iterate (a
    2 a -> a -> a
    forall a. Num a => a -> a -> a
    *) a
    1))
            ([a] -> [a]
    forall a. HasCallStack => [a] -> [a]
    cycle [a
    1, a
    1, a
    1, a
    0])

        dups :: forall x . [x] -> [x]
        dups :: forall a. [a] -> [a]
    dups = (x -> [x] -> [x]) -> [x] -> [x] -> [x]
    forall a b. (a -> b -> b) -> b -> [a] -> b
    forall (t :: * -> *) a b.
    Foldable t =>
    (a -> b -> b) -> b -> t a -> b
    foldr (\x
    n [x]
    list -> x
    n x -> [x] -> [x]
    forall a. a -> [a] -> [a]
    : x
    n x -> [x] -> [x]
    forall a. a -> [a] -> [a]
    : [x]
    list) []
    {-# SPECIALIZE euler' :: [Ratio Int]     #-}
    {-# SPECIALIZE euler' :: [Rational]      #-}

    -- | The same sequence as @euler'@, but with type @[a]@ instead of @[Ratio a]@
    -- as the denominators in @euler'@ are always @1@.
    --
    -- >>> take 10 euler :: [Integer]
    -- [1,0,-1,0,5,0,-61,0,1385,0]
    euler :: forall a . Integral a => [a]
    euler :: forall a. Integral a => [a]
    euler = (Ratio a -> a) -> [Ratio a] -> [a]
    forall a b. (a -> b) -> [a] -> [b]
    map Ratio a -> a
    forall a. Ratio a -> a
    numerator [Ratio a]
    forall a. Integral a => [Ratio a]
    euler'

    -- | Infinite zero-based list of the @n@-th order Euler polynomials evaluated at @1@.
    -- The algorithm used was derived from <http://www.emis.ams.org/journals/JIS/VOL4/CHEN/AlgBE2.pdf Algorithms for Bernoulli numbers and Euler numbers>
    -- by Kwang-Wu Chen, third formula of the Corollary in page 7.
    -- Element-by-element division of sequences <https://oeis.org/A198631 A1986631>
    -- and <https://oeis.org/A006519 A006519> in OEIS.
    --
    -- >>> take 10 eulerPolyAt1 :: [Rational]
    -- [1 % 1,1 % 2,0 % 1,(-1) % 4,0 % 1,1 % 2,0 % 1,(-17) % 8,0 % 1,31 % 2]
    eulerPolyAt1 :: forall a . Integral a => [Ratio a]
    eulerPolyAt1 :: forall a. Integral a => [Ratio a]
    eulerPolyAt1 = [Ratio a] -> [Ratio a]
    forall a. HasCallStack => [a] -> [a]
    tail ([Ratio a] -> [Ratio a]) -> [Ratio a] -> [Ratio a]
    forall a b. (a -> b) -> a -> b
    $ ([Ratio a] -> [Ratio a]) -> [Ratio a] -> [Ratio a]
    forall a.
    Integral a =>
    ([Ratio a] -> [Ratio a]) -> [Ratio a] -> [Ratio a]
    helperForBEEP [Ratio a] -> [Ratio a]
    forall a. HasCallStack => [a] -> [a]
    tail ((Ratio a -> Ratio a) -> [Ratio a] -> [Ratio a]
    forall a b. (a -> b) -> [a] -> [b]
    map Ratio a -> Ratio a
    forall a. Fractional a => a -> a
    recip ((Ratio a -> Ratio a) -> Ratio a -> [Ratio a]
    forall a. (a -> a) -> a -> [a]
    iterate (Ratio a
    2 Ratio a -> Ratio a -> Ratio a
    forall a. Num a => a -> a -> a
    *) Ratio a
    1))
    {-# SPECIALIZE eulerPolyAt1 :: [Ratio Int]     #-}
    {-# SPECIALIZE eulerPolyAt1 :: [Rational]      #-}

    -------------------------------------------------------------------------------
    -- Utils

    -- zipIndexedListWithTail f n as a == zipWith3 f [n..] as (tail as ++ [a])
    -- but inlines much better and avoids checks for distinct sizes of lists.
    zipIndexedListWithTail :: Enum b => (b -> a -> a -> b) -> b -> [a] -> a -> [b]
    zipIndexedListWithTail :: forall b a. Enum b => (b -> a -> a -> b) -> b -> [a] -> a -> [b]
    zipIndexedListWithTail b -> a -> a -> b
    f b
    n [a]
    as a
    a = case [a]
    as of
      []       -> []
      (a
    x : [a]
    xs) -> b -> a -> [a] -> [b]
    go b
    n a
    x [a]
    xs
      where
        go :: b -> a -> [a] -> [b]
    go b
    m a
    y [a]
    ys = case [a]
    ys of
          []       -> let v :: b
    v = b -> a -> a -> b
    f b
    m a
    y a
    a in [b
    v]
          (a
    z : [a]
    zs) -> let v :: b
    v = b -> a -> a -> b
    f b
    m a
    y a
    z in (b
    v b -> [b] -> [b]
    forall a. a -> [a] -> [a]
    : b -> a -> [a] -> [b]
    go (b -> b
    forall a. Enum a => a -> a
    succ b
    m) a
    z [a]
    zs)
    {-# INLINE zipIndexedListWithTail #-}

    -- | Helper for common code in @bernoulli, euler, eulerPolyAt1. All three
    -- sequences rely on @stirling2@ and have the same general structure of
    -- zipping four lists together with multiplication, with one of those lists
    -- being the sublists in @stirling2@, and two of them being the factorial
    -- sequence and @cycle [1, -1]@. The remaining list is passed to
    -- @helperForBEEP@ as an argument.
    --
    -- Note: This function has a @([Ratio a] -> [Ratio a])@ argument because
    -- @bernoulli !! n@ will use, for all nonnegative @n@, every element in
    -- @stirling2 !! n@, while @euler, eulerPolyAt1@ only use
    -- @tail $ stirling2 !! n@. As such, this argument serves to pass @id@
    -- in the former case, and @tail@ in the latter.
    helperForBEEP :: Integral a => ([Ratio a] -> [Ratio a]) -> [Ratio a] -> [Ratio a]
    helperForBEEP :: forall a.
    Integral a =>
    ([Ratio a] -> [Ratio a]) -> [Ratio a] -> [Ratio a]
    helperForBEEP [Ratio a] -> [Ratio a]
    g [Ratio a]
    xs = ([Ratio a] -> Ratio a) -> [[Ratio a]] -> [Ratio a]
    forall a b. (a -> b) -> [a] -> [b]
    map ([Ratio a] -> Ratio a
    f ([Ratio a] -> Ratio a)
    -> ([Ratio a] -> [Ratio a]) -> [Ratio a] -> Ratio a
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . [Ratio a] -> [Ratio a]
    g) [[Ratio a]]
    forall a. (Num a, Enum a) => [[a]]
    stirling2
      where
        f :: [Ratio a] -> Ratio a
    f = [Ratio a] -> Ratio a
    forall a. Num a => [a] -> a
    forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
    sum ([Ratio a] -> Ratio a)
    -> ([Ratio a] -> [Ratio a]) -> [Ratio a] -> Ratio a
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . (Ratio a -> Ratio a -> Ratio a -> Ratio a -> Ratio a)
    -> [Ratio a] -> [Ratio a] -> [Ratio a] -> [Ratio a] -> [Ratio a]
    forall a b c d e.
    (a -> b -> c -> d -> e) -> [a] -> [b] -> [c] -> [d] -> [e]
    zipWith4 (\Ratio a
    sgn Ratio a
    fact Ratio a
    x Ratio a
    stir -> Ratio a
    sgn Ratio a -> Ratio a -> Ratio a
    forall a. Num a => a -> a -> a
    * Ratio a
    fact Ratio a -> Ratio a -> Ratio a
    forall a. Num a => a -> a -> a
    * Ratio a
    x Ratio a -> Ratio a -> Ratio a
    forall a. Num a => a -> a -> a
    * Ratio a
    stir) ([Ratio a] -> [Ratio a]
    forall a. HasCallStack => [a] -> [a]
    cycle [Ratio a
    1, -Ratio a
    1]) [Ratio a]
    forall a. (Num a, Enum a) => [a]
    factorial [Ratio a]
    xs
    {-# SPECIALIZE helperForBEEP :: ([Ratio Int] -> [Ratio Int]) -> [Ratio Int] -> [Ratio Int] #-}
    {-# SPECIALIZE helperForBEEP :: ([Rational] -> [Rational]) -> [Rational] -> [Rational]     #-}