    -- |
    -- Module:      Math.NumberTheory.Zeta.Hurwitz
    -- Copyright:   (c) 2018 Alexandre Rodrigues Baldé
    -- Licence:     MIT
    -- Maintainer:  Alexandre Rodrigues Baldé <alexandrer_b@outlook.com>
    --
    -- Hurwitz zeta function.

    {-# LANGUAGE PostfixOperators    #-}
    {-# LANGUAGE ScopedTypeVariables #-}

    module Math.NumberTheory.Zeta.Hurwitz
      ( zetaHurwitz
      ) where

    import Data.List.Infinite (Infinite(..), (....))
    import qualified Data.List.Infinite as Inf

    import Math.NumberTheory.Recurrences (bernoulli, factorial)
    import Math.NumberTheory.Zeta.Utils  (skipEvens, skipOdds)

    -- | Values of Hurwitz zeta function evaluated at @ζ(s, a)@ for @s ∈ [0, 1 ..]@.
    --
    -- The algorithm used was based on the Euler-Maclaurin formula and was derived
    -- from <http://fredrikj.net/thesis/thesis.pdf Fast and Rigorous Computation of Special Functions to High Precision>
    -- by F. Johansson, chapter 4.8, formula 4.8.5.
    -- The error for each value in this recurrence is given in formula 4.8.9 as an
    --  indefinite integral, and in formula 4.8.12 as a closed form formula.
    --
    -- It is the __user's responsibility__ to provide an appropriate precision for
    -- the type chosen.
    --
    -- For instance, when using @Double@s, it does not make sense
    -- to provide a number @ε < 1e-53@ as the desired precision. For @Float@s,
    -- providing an @ε < 1e-24@ also does not make sense.
    -- Example of how to call the function:
    --
    -- >>> zetaHurwitz 1e-15 0.25 !! 5
    -- 1024.3489745265808
    zetaHurwitz :: forall a . (Floating a, Ord a) => a -> a -> Infinite a
    zetaHurwitz :: forall a. (Floating a, Ord a) => a -> a -> Infinite a
    zetaHurwitz a
    eps a
    a = (a -> a -> a -> a)
    -> Infinite a -> Infinite a -> Infinite a -> Infinite a
    forall a b c d.
    (a -> b -> c -> d)
    -> Infinite a -> Infinite b -> Infinite c -> Infinite d
    Inf.zipWith3 (\a
    s a
    i a
    t -> a
    s a -> a -> a
    forall a. Num a => a -> a -> a
    + a
    i a -> a -> a
    forall a. Num a => a -> a -> a
    + a
    t) Infinite a
    ss Infinite a
    is Infinite a
    ts
      where
        -- When given @1e-14@ as the @eps@ argument, this'll be
        -- @div (33 * (length . takeWhile (>= 1) . iterate (/ 10) . recip) 1e-14) 10 == div (33 * 14) 10@
        -- @div (33 * 14) 10 == 46.
        -- meaning @N,M@ in formula 4.8.5 will be @46@.
        -- Multiplying by 33 and dividing by 10 is because asking for @14@ digits
        -- of decimal precision equals asking for @(log 10 / log 2) * 14 ~ 3.3 * 14 ~ 46@
        -- bits of precision.
        digitsOfPrecision :: Integer
        digitsOfPrecision :: Integer
    digitsOfPrecision =
           let magnitude :: Integer
    magnitude = Int -> Integer
    forall a. Integral a => a -> Integer
    toInteger (Int -> Integer) -> (a -> Int) -> a -> Integer
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . [a] -> Int
    forall a. [a] -> Int
    forall (t :: * -> *) a. Foldable t => t a -> Int
    length ([a] -> Int) -> (a -> [a]) -> a -> Int
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . (a -> Bool) -> [a] -> [a]
    forall a. (a -> Bool) -> [a] -> [a]
    takeWhile (a -> a -> Bool
    forall a. Ord a => a -> a -> Bool
    >= a
    1) ([a] -> [a]) -> (a -> [a]) -> a -> [a]
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . (a -> a) -> a -> [a]
    forall a. (a -> a) -> a -> [a]
    iterate (a -> a -> a
    forall a. Fractional a => a -> a -> a
    / a
    10) (a -> [a]) -> (a -> a) -> a -> [a]
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . a -> a
    forall a. Fractional a => a -> a
    recip (a -> Integer) -> a -> Integer
    forall a b. (a -> b) -> a -> b
    $ a
    eps
           in  Integer -> Integer -> Integer
    forall a. Integral a => a -> a -> a
    div (Integer
    magnitude Integer -> Integer -> Integer
    forall a. Num a => a -> a -> a
    * Integer
    33) Integer
    10

        -- @a + n@
        aPlusN :: a
        aPlusN :: a
    aPlusN = a
    a a -> a -> a
    forall a. Num a => a -> a -> a
    + Integer -> a
    forall a. Num a => Integer -> a
    fromInteger Integer
    digitsOfPrecision

        -- @[(a + n)^s | s <- [0, 1, 2 ..]]@
        powsOfAPlusN :: Infinite a
        powsOfAPlusN :: Infinite a
    powsOfAPlusN = (a -> a) -> a -> Infinite a
    forall a. (a -> a) -> a -> Infinite a
    Inf.iterate (a
    aPlusN *) a
    1

        -- [                   [      1      ] |                   ]
        -- | \sum_{k=0}^\(n-1) | ----------- | | s <- [0, 1, 2 ..] |
        -- [                   [ (a + k) ^ s ] |                   ]
        -- @S@ value in 4.8.5 formula.
        ss :: Infinite a
        ss :: Infinite a
    ss = let numbers :: [a]
    numbers = (Integer -> a) -> [Integer] -> [a]
    forall a b. (a -> b) -> [a] -> [b]
    map ((a
    a +) (a -> a) -> (Integer -> a) -> Integer -> a
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Integer -> a
    forall a. Num a => Integer -> a
    fromInteger) [Integer
    0..Integer
    digitsOfPrecisionInteger -> Integer -> Integer
    forall a. Num a => a -> a -> a
    -Integer
    1]
                 denoms :: Infinite [a]
    denoms  = Int -> a -> [a]
    forall a. Int -> a -> [a]
    replicate (Integer -> Int
    forall a. Num a => Integer -> a
    fromInteger Integer
    digitsOfPrecision) a
    1 [a] -> Infinite [a] -> Infinite [a]
    forall a. a -> Infinite a -> Infinite a
    :<
                           ([a] -> [a]) -> [a] -> Infinite [a]
    forall a. (a -> a) -> a -> Infinite a
    Inf.iterate ((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]
    numbers) [a]
    numbers
             in ([a] -> a) -> Infinite [a] -> Infinite a
    forall a b. (a -> b) -> Infinite a -> Infinite b
    Inf.map ([a] -> a
    forall a. Num a => [a] -> a
    forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
    sum ([a] -> a) -> ([a] -> [a]) -> [a] -> a
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . (a -> a) -> [a] -> [a]
    forall a b. (a -> b) -> [a] -> [b]
    map a -> a
    forall a. Fractional a => a -> a
    recip) Infinite [a]
    denoms

        -- [ (a + n) ^ (1 - s)            a + n         |                   ]
        -- | ----------------- = ---------------------- | s <- [0, 1, 2 ..] |
        -- [       s - 1          (a + n) ^ s * (s - 1) |                   ]
        -- @I@ value in 4.8.5 formula.
        is :: Infinite a
        is :: Infinite a
    is = let denoms :: Infinite a
    denoms = (a -> Integer -> a) -> Infinite a -> Infinite Integer -> Infinite a
    forall a b c.
    (a -> b -> c) -> Infinite a -> Infinite b -> Infinite c
    Inf.zipWith
                          (\a
    powOfA Integer
    int -> a
    powOfA a -> a -> a
    forall a. Num a => a -> a -> a
    * Integer -> a
    forall a. Num a => Integer -> a
    fromInteger Integer
    int)
                          Infinite a
    powsOfAPlusN
                          ((-Integer
    1, Integer
    0)....)
             in (a -> a) -> Infinite a -> Infinite a
    forall a b. (a -> b) -> Infinite a -> Infinite b
    Inf.map (a
    aPlusN /) Infinite a
    denoms

        -- [      1      |             ]
        -- [ ----------- | s <- [0 ..] ]
        -- [ (a + n) ^ s |             ]
        constants2 :: Infinite a
        constants2 :: Infinite a
    constants2 = (a -> a) -> Infinite a -> Infinite a
    forall a b. (a -> b) -> Infinite a -> Infinite b
    Inf.map a -> a
    forall a. Fractional a => a -> a
    recip Infinite a
    powsOfAPlusN

        -- [ [(s)_(2*k - 1) | k <- [1 ..]], s <- [0 ..]], i.e. odd indices of
        -- infinite rising factorial sequences, each sequence starting at a
        -- positive integer.
        pochhammers :: Infinite (Infinite Integer)
        pochhammers :: Infinite (Infinite Integer)
    pochhammers = let -- [ [(s)_k | k <- [1 ..]], s <- [1 ..]]
                          pochhs :: Infinite (Infinite Integer)
                          pochhs :: Infinite (Infinite Integer)
    pochhs = (Infinite Integer -> Infinite Integer)
    -> Infinite Integer -> Infinite (Infinite Integer)
    forall a. (a -> a) -> a -> Infinite a
    Inf.iterate (\(Integer
    x :< Infinite Integer
    xs) -> (Integer -> Integer) -> Infinite Integer -> Infinite Integer
    forall a b. (a -> b) -> Infinite a -> Infinite b
    Inf.map (Integer -> Integer -> Integer
    forall a. Integral a => a -> a -> a
    `div` Integer
    x) Infinite Integer
    xs) (Infinite Integer -> Infinite Integer
    forall a. Infinite a -> Infinite a
    Inf.tail Infinite Integer
    forall a. (Num a, Enum a) => Infinite a
    factorial)
                      in -- When @s@ is @0@, the infinite sequence of rising
                         -- factorials starting at @s@ is @[0,0,0,0..]@.
                         Integer -> Infinite Integer
    forall a. a -> Infinite a
    Inf.repeat Integer
    0 Infinite Integer
    -> Infinite (Infinite Integer) -> Infinite (Infinite Integer)
    forall a. a -> Infinite a -> Infinite a
    :< (Infinite Integer -> Infinite Integer)
    -> Infinite (Infinite Integer) -> Infinite (Infinite Integer)
    forall a b. (a -> b) -> Infinite a -> Infinite b
    Inf.map Infinite Integer -> Infinite Integer
    forall a. Infinite a -> Infinite a
    skipOdds Infinite (Infinite Integer)
    pochhs

        -- [            B_2k           |             ]
        -- | ------------------------- | k <- [1 ..] |
        -- [ (2k)! (a + n) ^ (2*k - 1) |             ]
        second :: [a]
        second :: [a]
    second =
            Int -> Infinite a -> [a]
    forall a. Int -> Infinite a -> [a]
    Inf.take (Integer -> Int
    forall a. Num a => Integer -> a
    fromInteger Integer
    digitsOfPrecision) (Infinite a -> [a]) -> Infinite a -> [a]
    forall a b. (a -> b) -> a -> b
    $
            (Rational -> Integer -> a -> a)
    -> Infinite Rational
    -> Infinite Integer
    -> Infinite a
    -> Infinite a
    forall a b c d.
    (a -> b -> c -> d)
    -> Infinite a -> Infinite b -> Infinite c -> Infinite d
    Inf.zipWith3
            (\Rational
    bern Integer
    evenFac a
    denom -> Rational -> a
    forall a. Fractional a => Rational -> a
    fromRational Rational
    bern a -> a -> a
    forall a. Fractional a => a -> a -> a
    / (a
    denom a -> a -> a
    forall a. Num a => a -> a -> a
    * Integer -> a
    forall a. Num a => Integer -> a
    fromInteger Integer
    evenFac))
            (Infinite Rational -> Infinite Rational
    forall a. Infinite a -> Infinite a
    Inf.tail (Infinite Rational -> Infinite Rational)
    -> Infinite Rational -> Infinite Rational
    forall a b. (a -> b) -> a -> b
    $ Infinite Rational -> Infinite Rational
    forall a. Infinite a -> Infinite a
    skipOdds Infinite Rational
    forall a. Integral a => Infinite (Ratio a)
    bernoulli)
            (Infinite Integer -> Infinite Integer
    forall a. Infinite a -> Infinite a
    Inf.tail (Infinite Integer -> Infinite Integer)
    -> Infinite Integer -> Infinite Integer
    forall a b. (a -> b) -> a -> b
    $ Infinite Integer -> Infinite Integer
    forall a. Infinite a -> Infinite a
    skipOdds Infinite Integer
    forall a. (Num a, Enum a) => Infinite a
    factorial)
            -- Recall that @powsOfAPlusN = [(a + n) ^ s | s <- [0 ..]]@, so this
            -- is @[(a + n) ^ (2 * s - 1) | s <- [1 ..]]@
            (Infinite a -> Infinite a
    forall a. Infinite a -> Infinite a
    skipEvens Infinite a
    powsOfAPlusN)

        fracs :: Infinite a
        fracs :: Infinite a
    fracs = (Infinite Integer -> a)
    -> Infinite (Infinite Integer) -> Infinite a
    forall a b. (a -> b) -> Infinite a -> Infinite b
    Inf.map
                ([a] -> a
    forall a. Num a => [a] -> a
    forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
    sum ([a] -> a) -> (Infinite Integer -> [a]) -> Infinite Integer -> a
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . (a -> Integer -> a) -> [a] -> [Integer] -> [a]
    forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
    zipWith (\a
    s Integer
    p -> a
    s a -> a -> a
    forall a. Num a => a -> a -> a
    * Integer -> a
    forall a. Num a => Integer -> a
    fromInteger Integer
    p) [a]
    second ([Integer] -> [a])
    -> (Infinite Integer -> [Integer]) -> Infinite Integer -> [a]
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Infinite Integer -> [Integer]
    forall a. Infinite a -> [a]
    Inf.toList)
                Infinite (Infinite Integer)
    pochhammers

        -- Infinite list of @T@ values in 4.8.5 formula, for every @s@ in
        -- @[0, 1, 2 ..]@.
        ts :: Infinite a
        ts :: Infinite a
    ts = (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
    constant2 a
    frac -> a
    constant2 a -> a -> a
    forall a. Num a => a -> a -> a
    * (a
    0.5 a -> a -> a
    forall a. Num a => a -> a -> a
    + a
    frac))
             Infinite a
    constants2
             Infinite a
    fracs