-- |
-- 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 = forall a. (a -> a) -> a -> Infinite a
Inf.iterate (\[a]
l -> forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall a. Semiring a => a -> a -> a
plus ([a]
l forall a. [a] -> [a] -> [a]
++ [forall a. Semiring a => a
zero]) (forall a. Semiring a => a
zero forall a. a -> [a] -> [a]
: [a]
l)) [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 = forall a. (a -> a) -> a -> Infinite a
Inf.iterate (forall a. Infinite a -> Infinite a
Inf.tail forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall b a. (b -> a -> b) -> b -> Infinite a -> Infinite b
Inf.scanl' forall a. Semiring a => a -> a -> a
plus forall a. Semiring a => a
zero) (forall a. a -> Infinite a
Inf.repeat 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 = forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl'
  (\a
x (a
k, a
nk1) -> forall a. HasCallStack => Maybe a -> a
fromJust forall a b. (a -> b) -> a -> b
$ (a
x forall a. Semiring a => a -> a -> a
`times` a
nk1) forall a. GcdDomain a => a -> a -> Maybe a
`divide` a
k)
  forall a. Semiring a => a
one
  (forall a b. [a] -> [b] -> [(a, b)]
zip [forall a. Semiring a => a
one..a
n] [a
n, forall a. Enum a => a -> a
pred a
n..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 = forall b a. (b -> a -> b) -> b -> Infinite a -> Infinite b
Inf.scanl'
  (\a
x a
k -> forall a. HasCallStack => Maybe a -> a
fromJust (a
x forall a. Semiring a => a -> a -> a
`times` (a
n forall a. Semiring a => a -> a -> a
`plus` a
k) forall a. GcdDomain a => a -> a -> Maybe a
`divide` a
k))
  forall a. Semiring a => a
one
  (forall a. Semiring a => a
one...)
{-# 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 forall a. Ord a => a -> a -> Bool
< Word
2
  = []
  | Bool
otherwise
  = forall a. (a -> Bool) -> [a] -> [a]
filter ((forall a. Eq a => a -> a -> Bool
/= Word
0) forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a b. (a, b) -> b
snd)
  forall a b. (a -> b) -> a -> b
$ forall a b. (a -> b) -> [a] -> [b]
map (\Prime Word
p -> (Prime Word
p, Word -> Word -> Word
mult (forall a. Prime a -> a
unPrime Prime Word
p) Word
n forall a. Num a => a -> a -> a
- Word -> Word -> Word
mult (forall a. Prime a -> a
unPrime Prime Word
p) (Word
n forall a. Num a => a -> a -> a
- Word
k) forall a. Num a => a -> a -> a
- Word -> Word -> Word
mult (forall a. Prime a -> a
unPrime Prime Word
p) Word
k))
    [forall a. Bounded a => a
minBound .. 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 forall a. Integral a => a -> a -> a
`quot` Word
p
        go :: Word -> Word -> Word
go !Word
acc !Word
x
          | Word
x forall a. Ord a => a -> a -> Bool
>= Word
p = let xp :: Word
xp = Word
x forall a. Integral a => a -> a -> a
`quot` Word
p in Word -> Word -> Word
go (Word
acc 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 = forall b a. (b -> a -> b) -> b -> Infinite a -> Infinite b
Inf.scanl 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 forall 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 forall a. Num a => a -> a -> a
+ a
n 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 = forall a. (a -> a) -> a -> Infinite a
Inf.iterate forall {a}. (Num a, Enum a) => [a] -> [a]
f [a
1]
  where
    f :: [a] -> [a]
f [a]
xs = a
0 forall 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 forall a. Num a => a -> a -> a
+ a
k 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 = forall a b c.
(a -> b -> c) -> Infinite a -> Infinite b -> Infinite c
Inf.zipWith forall {b}. Integral b => b -> b -> [b]
f (forall a. Infinite a -> Infinite a
Inf.tail forall a. (Num a, Enum a) => Infinite a
factorial) (a
1...)
  where
    f :: b -> b -> [b]
f b
nf b
n = forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl (\b
x b
k -> b
x forall a. Num a => a -> a -> a
* (b
n forall a. Num a => a -> a -> a
- b
k) forall a. Integral a => a -> a -> a
`div` (b
k forall a. Num a => a -> a -> a
* (b
k forall a. Num a => a -> a -> a
+ b
1))) b
nf [b
1..b
nforall 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 = forall b a. (b -> a -> b) -> b -> Infinite a -> Infinite b
Inf.scanl forall {a}. (Num a, Enum a) => [a] -> a -> [a]
f [] (a
1...)
  where
    f :: [a] -> a -> [a]
f [a]
xs a
n = a
1 forall 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 forall a. Num a => a -> a -> a
- a
k) forall a. Num a => a -> a -> a
* a
x forall a. Num a => a -> a -> a
+ (a
k forall a. Num a => a -> a -> a
+ a
1) 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 = forall b a. (b -> a -> b) -> b -> Infinite a -> Infinite b
Inf.scanl forall {a}. (Num a, Enum a) => [a] -> a -> [a]
f [] (a
1...)
  where
    f :: [a] -> a -> [a]
f [a]
xs a
n = a
1 forall 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 forall a. Num a => a -> a -> a
* a
n forall a. Num a => a -> a -> a
- a
k forall a. Num a => a -> a -> a
- a
1) forall a. Num a => a -> a -> a
* a
x forall a. Num a => a -> a -> a
+ (a
k forall a. Num a => a -> a -> a
+ a
1) 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 a.
Integral a =>
(forall b. [b] -> [b]) -> Infinite (Ratio a) -> Infinite (Ratio a)
helperForBEEP forall a. a -> a
id (forall a b. (a -> b) -> Infinite a -> Infinite b
Inf.map forall a. Fractional a => a -> a
recip (Ratio a
1...))
{-# 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
  = forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall a. Num a => a -> a -> a
(*) ((Ratio a
0:)
  forall a b. (a -> b) -> a -> b
$ forall b. [b] -> [b]
reverse
  forall a b. (a -> b) -> a -> b
$ forall a. Int -> Infinite a -> [a]
Inf.take (Int
p forall a. Num a => a -> a -> a
+ Int
1) forall a. Integral a => Infinite (Ratio a)
bernoulli)
  forall a b. (a -> b) -> a -> b
$ forall a b. (a -> b) -> [a] -> [b]
map (forall a. Integral a => a -> a -> Ratio a
% (forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
pforall a. Num a => a -> a -> a
+a
1))
  forall a b. (a -> b) -> a -> b
$ forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall a. Num a => a -> a -> a
(*) (forall a. (a -> a) -> a -> [a]
iterate forall a. Num a => a -> a
negate (if forall a. Integral a => a -> Bool
odd Int
p then a
1 else -a
1))
  forall a b. (a -> b) -> a -> b
$ forall a. Semiring a => Infinite [a]
binomial forall a. Infinite a -> Word -> a
Inf.!! (forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
pforall 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' = forall a. Infinite a -> Infinite a
Inf.tail forall a b. (a -> b) -> a -> b
$ forall a.
Integral a =>
(forall b. [b] -> [b]) -> Infinite (Ratio a) -> Infinite (Ratio a)
helperForBEEP (forall a. Int -> [a] -> [a]
drop Int
1) Infinite (Ratio a)
as
  where
    as :: Infinite (Ratio a)
    as :: Infinite (Ratio a)
as = 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 forall a. Num a => a -> a -> a
* a
ones) forall a. Integral a => a -> a -> Ratio a
% a
frac)
        (forall a. NonEmpty a -> Infinite a
Inf.cycle (a
1 forall a. a -> [a] -> NonEmpty a
:| [a
1, a
1, a
1, -a
1, -a
1, -a
1, -a
1]))
        (forall a. Infinite a -> Infinite a
dups (forall a. (a -> a) -> a -> Infinite a
Inf.iterate (a
2 *) a
1))
        (forall a. NonEmpty a -> Infinite a
Inf.cycle (a
1 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 = forall a b. (a -> b -> b) -> Infinite a -> b
Inf.foldr (\x
n Infinite x
list -> x
n forall a. a -> Infinite a -> Infinite a
:< x
n 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 = forall a b. (a -> b) -> Infinite a -> Infinite b
Inf.map forall a. Ratio a -> a
numerator 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 = forall a. Infinite a -> Infinite a
Inf.tail forall a b. (a -> b) -> a -> b
$ forall a.
Integral a =>
(forall b. [b] -> [b]) -> Infinite (Ratio a) -> Infinite (Ratio a)
helperForBEEP (forall a. Int -> [a] -> [a]
drop Int
1) (forall a b. (a -> b) -> Infinite a -> Infinite b
Inf.map forall a. Fractional a => a -> a
recip (forall a. (a -> a) -> a -> Infinite a
Inf.iterate (Ratio a
2 *) 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 forall a. a -> [a] -> [a]
: b -> a -> [a] -> [b]
go (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 = forall a b. (a -> b) -> Infinite a -> Infinite b
Inf.map ([Ratio a] -> Ratio a
f forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall b. [b] -> [b]
g) forall a. (Num a, Enum a) => Infinite [a]
stirling2
  where
    f :: [Ratio a] -> Ratio a
f = forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum forall b c a. (b -> c) -> (a -> b) -> a -> c
. 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 forall a. Num a => a -> a -> a
* Ratio a
fact forall a. Num a => a -> a -> a
* Ratio a
x forall a. Num a => a -> a -> a
* Ratio a
stir) (forall b. [b] -> [b]
cycle [Ratio a
1, -Ratio a
1]) (forall a. Infinite a -> [a]
Inf.toList forall a. (Num a, Enum a) => Infinite a
factorial) (forall a. Infinite a -> [a]
Inf.toList Infinite (Ratio a)
xs)
{-# SPECIALIZE helperForBEEP :: (forall b. [b] -> [b]) -> Infinite (Ratio Int) -> Infinite (Ratio Int) #-}
{-# SPECIALIZE helperForBEEP :: (forall b. [b] -> [b]) -> Infinite (Rational) -> Infinite (Rational)     #-}