-- |
-- 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