-- |
-- 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 :: [[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 :: [[a]]
binomialRotated = ([a] -> [a]) -> [a] -> [[a]]
forall a. (a -> a) -> a -> [a]
iterate ([a] -> [a]
forall a. [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 :: 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 :: 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 :: [[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 :: [[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 :: [[a]]
lah = (a -> a -> [a]) -> [a] -> [a] -> [[a]]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> a -> [a]
forall a. Integral a => a -> a -> [a]
f ([a] -> [a]
forall a. [a] -> [a]
tail [a]
forall a. (Num a, Enum a) => [a]
factorial) [a
1..]
  where
    f :: a -> a -> [a]
f a
nf a
n = (a -> a -> a) -> a -> [a] -> [a]
forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl (\a
x a
k -> 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
k) a -> a -> a
forall a. Integral a => a -> a -> a
`div` (a
k 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
nf [a
1..a
na -> a -> a
forall a. Num a => a -> a -> a
-a
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 :: [[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 :: [[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 :: [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 :: 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. [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' :: [Ratio a]
euler' = [Ratio a] -> [Ratio a]
forall a. [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. [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. [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. [a] -> [a]
cycle [a
1, a
1, a
1, a
0])

    dups :: forall x . [x] -> [x]
    dups :: [x] -> [x]
dups = (x -> [x] -> [x]) -> [x] -> [x] -> [x]
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 :: [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 :: [Ratio a]
eulerPolyAt1 = [Ratio a] -> [Ratio a]
forall a. [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. [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 :: (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 :: ([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 (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. [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]     #-}