-- |
-- 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
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 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
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 :: 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
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 :: 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
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 :: 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
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 :: 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
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
  = (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
0:)
  ([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
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
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)
{-# SPECIALIZE helperForBEEP :: (forall b. [b] -> [b]) -> Infinite (Ratio Int) -> Infinite (Ratio Int) #-}
{-# SPECIALIZE helperForBEEP :: (forall b. [b] -> [b]) -> Infinite (Rational) -> Infinite (Rational)     #-}