    -- |
    -- 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 PostfixOperators    #-}
    {-# LANGUAGE RankNTypes          #-}
    {-# 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.List.Infinite (Infinite(..), (...))
    import qualified Data.List.Infinite as Inf
    import Data.List.NonEmpty (NonEmpty(..))
    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 => Infinite [a]
    binomial :: forall a. Semiring a => Infinite [a]
    binomial = ([a] -> [a]) -> [a] -> Infinite [a]
    forall a. (a -> a) -> a -> Infinite a
    Inf.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 :: Infinite [Int]     #-}
    {-# SPECIALIZE binomial :: Infinite [Word]    #-}
    {-# SPECIALIZE binomial :: Infinite [Integer] #-}
    {-# SPECIALIZE binomial :: Infinite [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 => Infinite (Infinite a)
    binomialRotated :: forall a. Semiring a => Infinite (Infinite a)
    binomialRotated = (Infinite a -> Infinite a) -> Infinite a -> Infinite (Infinite a)
    forall a. (a -> a) -> a -> Infinite a
    Inf.iterate (Infinite a -> Infinite a
    forall a. Infinite a -> Infinite a
    Inf.tail (Infinite a -> Infinite a)
    -> (Infinite a -> Infinite a) -> Infinite a -> Infinite a
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . (a -> a -> a) -> a -> Infinite a -> Infinite a
    forall b a. (b -> a -> b) -> b -> Infinite a -> Infinite b
    Inf.scanl' a -> a -> a
    forall a. Semiring a => a -> a -> a
    plus a
    forall a. Semiring a => a
    zero) (a -> Infinite a
    forall a. a -> Infinite a
    Inf.repeat a
    forall a. Semiring a => a
    one)
    {-# SPECIALIZE binomialRotated :: Infinite (Infinite Int)     #-}
    {-# SPECIALIZE binomialRotated :: Infinite (Infinite Word)    #-}
    {-# SPECIALIZE binomialRotated :: Infinite (Infinite Integer) #-}
    {-# SPECIALIZE binomialRotated :: Infinite (Infinite 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 -> Infinite a
    binomialDiagonal :: forall a. (Enum a, GcdDomain a) => a -> Infinite a
    binomialDiagonal a
    n = (a -> a -> a) -> a -> Infinite a -> Infinite a
    forall b a. (b -> a -> b) -> b -> Infinite a -> Infinite b
    Inf.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
    onea -> Infinite a
    forall a. Enum a => a -> Infinite a
    ...)
    {-# SPECIALIZE binomialDiagonal :: Int     -> Infinite Int     #-}
    {-# SPECIALIZE binomialDiagonal :: Word    -> Infinite Word    #-}
    {-# SPECIALIZE binomialDiagonal :: Integer -> Infinite Integer #-}
    {-# SPECIALIZE binomialDiagonal :: Natural -> Infinite 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) => Infinite [a]
    stirling1 :: forall a. (Num a, Enum a) => Infinite [a]
    stirling1 = ([a] -> a -> [a]) -> [a] -> Infinite a -> Infinite [a]
    forall b a. (b -> a -> b) -> b -> Infinite a -> Infinite b
    Inf.scanl [a] -> a -> [a]
    forall {a}. (Num a, Enum a) => [a] -> a -> [a]
    f [a
    1] (a
    0a -> Infinite a
    forall a. Enum a => a -> Infinite a
    ...)
      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 :: Infinite [Int]     #-}
    {-# SPECIALIZE stirling1 :: Infinite [Word]    #-}
    {-# SPECIALIZE stirling1 :: Infinite [Integer] #-}
    {-# SPECIALIZE stirling1 :: Infinite [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) => Infinite [a]
    stirling2 :: forall a. (Num a, Enum a) => Infinite [a]
    stirling2 = ([a] -> [a]) -> [a] -> Infinite [a]
    forall a. (a -> a) -> a -> Infinite a
    Inf.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 :: Infinite [Int]     #-}
    {-# SPECIALIZE stirling2 :: Infinite [Word]    #-}
    {-# SPECIALIZE stirling2 :: Infinite [Integer] #-}
    {-# SPECIALIZE stirling2 :: Infinite [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 => Infinite [a]
    -- Implementation was derived from code by https://github.com/grandpascorpion
    lah :: forall a. Integral a => Infinite [a]
    lah = (a -> a -> [a]) -> Infinite a -> Infinite a -> Infinite [a]
    forall a b c.
    (a -> b -> c) -> Infinite a -> Infinite b -> Infinite c
    Inf.zipWith a -> a -> [a]
    forall {b}. Integral b => b -> b -> [b]
    f (Infinite a -> Infinite a
    forall a. Infinite a -> Infinite a
    Inf.tail Infinite a
    forall a. (Num a, Enum a) => Infinite a
    factorial) (a
    1a -> Infinite a
    forall a. Enum a => a -> Infinite a
    ...)
      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 :: Infinite [Int]     #-}
    {-# SPECIALIZE lah :: Infinite [Word]    #-}
    {-# SPECIALIZE lah :: Infinite [Integer] #-}
    {-# SPECIALIZE lah :: Infinite [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) => Infinite [a]
    eulerian1 :: forall a. (Num a, Enum a) => Infinite [a]
    eulerian1 = ([a] -> a -> [a]) -> [a] -> Infinite a -> Infinite [a]
    forall b a. (b -> a -> b) -> b -> Infinite a -> Infinite b
    Inf.scanl [a] -> a -> [a]
    forall {a}. (Num a, Enum a) => [a] -> a -> [a]
    f [] (a
    1a -> Infinite a
    forall a. Enum a => a -> Infinite a
    ...)
      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 :: Infinite [Int]     #-}
    {-# SPECIALIZE eulerian1 :: Infinite [Word]    #-}
    {-# SPECIALIZE eulerian1 :: Infinite [Integer] #-}
    {-# SPECIALIZE eulerian1 :: Infinite [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) => Infinite [a]
    eulerian2 :: forall a. (Num a, Enum a) => Infinite [a]
    eulerian2 = ([a] -> a -> [a]) -> [a] -> Infinite a -> Infinite [a]
    forall b a. (b -> a -> b) -> b -> Infinite a -> Infinite b
    Inf.scanl [a] -> a -> [a]
    forall {a}. (Num a, Enum a) => [a] -> a -> [a]
    f [] (a
    1a -> Infinite a
    forall a. Enum a => a -> Infinite a
    ...)
      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 :: Infinite [Int]     #-}
    {-# SPECIALIZE eulerian2 :: Infinite [Word]    #-}
    {-# SPECIALIZE eulerian2 :: Infinite [Integer] #-}
    {-# SPECIALIZE eulerian2 :: Infinite [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 => Infinite (Ratio a)
    bernoulli :: forall a. Integral a => Infinite (Ratio a)
    bernoulli = (forall b. [b] -> [b]) -> Infinite (Ratio a) -> Infinite (Ratio a)
    forall a.
    Integral a =>
    (forall b. [b] -> [b]) -> Infinite (Ratio a) -> Infinite (Ratio a)
    helperForBEEP [b] -> [b]
    forall a. a -> a
    forall b. [b] -> [b]
    id ((Ratio a -> Ratio a) -> Infinite (Ratio a) -> Infinite (Ratio a)
    forall a b. (a -> b) -> Infinite a -> Infinite b
    Inf.map Ratio a -> Ratio a
    forall a. Fractional a => a -> a
    recip (Ratio a
    1Ratio a -> Infinite (Ratio a)
    forall a. Enum a => a -> Infinite a
    ...))
    {-# SPECIALIZE bernoulli :: Infinite (Ratio Int) #-}
    {-# SPECIALIZE bernoulli :: Infinite (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 b. [b] -> [b]
    reverse
      ([Ratio a] -> [Ratio a]) -> [Ratio a] -> [Ratio a]
    forall a b. (a -> b) -> a -> b
    $ Int -> Infinite (Ratio a) -> [Ratio a]
    forall a. Int -> Infinite a -> [a]
    Inf.take (Int
    p Int -> Int -> Int
    forall a. Num a => a -> a -> a
    + Int
    1) Infinite (Ratio a)
    forall a. Integral a => Infinite (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
    $ Infinite [a]
    forall a. Semiring a => Infinite [a]
    binomial Infinite [a] -> Word -> [a]
    forall a. Infinite a -> Word -> a
    Inf.!! (Int -> Word
    forall a b. (Integral a, Num b) => a -> b
    fromIntegral (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 => Infinite (Ratio a)
    euler' :: forall a. Integral a => Infinite (Ratio a)
    euler' = Infinite (Ratio a) -> Infinite (Ratio a)
    forall a. Infinite a -> Infinite a
    Inf.tail (Infinite (Ratio a) -> Infinite (Ratio a))
    -> Infinite (Ratio a) -> Infinite (Ratio a)
    forall a b. (a -> b) -> a -> b
    $ (forall b. [b] -> [b]) -> Infinite (Ratio a) -> Infinite (Ratio a)
    forall a.
    Integral a =>
    (forall b. [b] -> [b]) -> Infinite (Ratio a) -> Infinite (Ratio a)
    helperForBEEP (Int -> [b] -> [b]
    forall a. Int -> [a] -> [a]
    drop Int
    1) Infinite (Ratio a)
    as
      where
        as :: Infinite (Ratio a)
        as :: Infinite (Ratio a)
    as = (a -> a -> a -> Ratio a)
    -> Infinite a -> Infinite a -> Infinite a -> Infinite (Ratio a)
    forall a b c d.
    (a -> b -> c -> d)
    -> Infinite a -> Infinite b -> Infinite c -> Infinite d
    Inf.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)
            (NonEmpty a -> Infinite a
    forall a. NonEmpty a -> Infinite a
    Inf.cycle (a
    1 a -> [a] -> NonEmpty a
    forall a. a -> [a] -> NonEmpty a
    :| [a
    1, a
    1, a
    1, -a
    1, -a
    1, -a
    1, -a
    1]))
            (Infinite a -> Infinite a
    forall a. Infinite a -> Infinite a
    dups ((a -> a) -> a -> Infinite a
    forall a. (a -> a) -> a -> Infinite a
    Inf.iterate (a
    2 a -> a -> a
    forall a. Num a => a -> a -> a
    *) a
    1))
            (NonEmpty a -> Infinite a
    forall a. NonEmpty a -> Infinite a
    Inf.cycle (a
    1 a -> [a] -> NonEmpty a
    forall a. a -> [a] -> NonEmpty a
    :| [a
    1, a
    1, a
    0]))

        dups :: forall x . Infinite x -> Infinite x
        dups :: forall a. Infinite a -> Infinite a
    dups = (x -> Infinite x -> Infinite x) -> Infinite x -> Infinite x
    forall a b. (a -> b -> b) -> Infinite a -> b
    Inf.foldr (\x
    n Infinite x
    list -> x
    n x -> Infinite x -> Infinite x
    forall a. a -> Infinite a -> Infinite a
    :< x
    n x -> Infinite x -> Infinite x
    forall a. a -> Infinite a -> Infinite a
    :< Infinite x
    list)
    {-# SPECIALIZE euler' :: Infinite (Ratio Int)     #-}
    {-# SPECIALIZE euler' :: Infinite (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 => Infinite a
    euler :: forall a. Integral a => Infinite a
    euler = (Ratio a -> a) -> Infinite (Ratio a) -> Infinite a
    forall a b. (a -> b) -> Infinite a -> Infinite b
    Inf.map Ratio a -> a
    forall a. Ratio a -> a
    numerator Infinite (Ratio a)
    forall a. Integral a => Infinite (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 => Infinite (Ratio a)
    eulerPolyAt1 :: forall a. Integral a => Infinite (Ratio a)
    eulerPolyAt1 = Infinite (Ratio a) -> Infinite (Ratio a)
    forall a. Infinite a -> Infinite a
    Inf.tail (Infinite (Ratio a) -> Infinite (Ratio a))
    -> Infinite (Ratio a) -> Infinite (Ratio a)
    forall a b. (a -> b) -> a -> b
    $ (forall b. [b] -> [b]) -> Infinite (Ratio a) -> Infinite (Ratio a)
    forall a.
    Integral a =>
    (forall b. [b] -> [b]) -> Infinite (Ratio a) -> Infinite (Ratio a)
    helperForBEEP (Int -> [b] -> [b]
    forall a. Int -> [a] -> [a]
    drop Int
    1) ((Ratio a -> Ratio a) -> Infinite (Ratio a) -> Infinite (Ratio a)
    forall a b. (a -> b) -> Infinite a -> Infinite b
    Inf.map Ratio a -> Ratio a
    forall a. Fractional a => a -> a
    recip ((Ratio a -> Ratio a) -> Ratio a -> Infinite (Ratio a)
    forall a. (a -> a) -> a -> Infinite a
    Inf.iterate (Ratio a
    2 Ratio a -> Ratio a -> Ratio a
    forall a. Num a => a -> a -> a
    *) Ratio a
    1))
    {-# SPECIALIZE eulerPolyAt1 :: Infinite (Ratio Int)     #-}
    {-# SPECIALIZE eulerPolyAt1 :: Infinite (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 => (forall b. [b] -> [b]) -> Infinite (Ratio a) -> Infinite (Ratio a)
    helperForBEEP :: forall a.
    Integral a =>
    (forall b. [b] -> [b]) -> Infinite (Ratio a) -> Infinite (Ratio a)
    helperForBEEP forall b. [b] -> [b]
    g Infinite (Ratio a)
    xs = ([Ratio a] -> Ratio a) -> Infinite [Ratio a] -> Infinite (Ratio a)
    forall a b. (a -> b) -> Infinite a -> Infinite b
    Inf.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]
    forall b. [b] -> [b]
    g) Infinite [Ratio a]
    forall a. (Num a, Enum a) => Infinite [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]) (Infinite (Ratio a) -> [Ratio a]
    forall a. Infinite a -> [a]
    Inf.toList Infinite (Ratio a)
    forall a. (Num a, Enum a) => Infinite a
    factorial) (Infinite (Ratio a) -> [Ratio a]
    forall a. Infinite a -> [a]
    Inf.toList Infinite (Ratio a)
    xs)
    {-# INLINABLE helperForBEEP #-}