    -- |
    -- Module:      Math.NumberTheory.Zeta.Riemann
    -- Copyright:   (c) 2016 Andrew Lelechenko
    -- Licence:     MIT
    -- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>
    --
    -- Riemann zeta-function.

    {-# LANGUAGE ScopedTypeVariables #-}

    module Math.NumberTheory.Zeta.Riemann
      ( zetas
      , zetasEven
      , zetasOdd
      ) where

    import Data.ExactPi
    import Data.Ratio                     ((%))

    import Math.NumberTheory.Recurrences  (bernoulli)
    import Math.NumberTheory.Zeta.Hurwitz (zetaHurwitz)
    import Math.NumberTheory.Zeta.Utils   (intertwine, skipEvens, skipOdds)

    -- | Infinite sequence of exact values of Riemann zeta-function at even arguments, starting with @ζ(0)@.
    -- Note that due to numerical errors conversion to 'Double' may return values below 1:
    --
    -- >>> approximateValue (zetasEven !! 25) :: Double
    -- 0.9999999999999996
    --
    -- Use your favorite type for long-precision arithmetic. For instance, 'Data.Number.Fixed.Fixed' works fine:
    --
    -- >>> import Data.Number.Fixed
    -- >>> approximateValue (zetasEven !! 25) :: Fixed Prec50
    -- 1.00000000000000088817842111574532859293035196051773
    --
    zetasEven :: [ExactPi]
    zetasEven :: [ExactPi]
    zetasEven = (Integer -> Rational -> ExactPi)
    -> [Integer] -> [Rational] -> [ExactPi]
    forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
    zipWith Integer -> Rational -> ExactPi
    Exact [Integer
    0, Integer
    2 ..] ([Rational] -> [ExactPi]) -> [Rational] -> [ExactPi]
    forall a b. (a -> b) -> a -> b
    $ (Rational -> Rational -> Rational)
    -> [Rational] -> [Rational] -> [Rational]
    forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
    zipWith Rational -> Rational -> Rational
    forall a. Num a => a -> a -> a
    (*) ([Rational] -> [Rational]
    forall a. [a] -> [a]
    skipOdds [Rational]
    forall a. Integral a => [Ratio a]
    bernoulli) [Rational]
    cs
      where
        cs :: [Rational]
    cs = (- Integer
    1 Integer -> Integer -> Rational
    forall a. Integral a => a -> a -> Ratio a
    % Integer
    2) Rational -> [Rational] -> [Rational]
    forall a. a -> [a] -> [a]
    : (Rational -> Integer -> Rational)
    -> [Rational] -> [Integer] -> [Rational]
    forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
    zipWith (\Rational
    i Integer
    f -> Rational
    i Rational -> Rational -> Rational
    forall a. Num a => a -> a -> a
    * (-Rational
    4) Rational -> Rational -> Rational
    forall a. Fractional a => a -> a -> a
    / Integer -> Rational
    forall a. Num a => Integer -> a
    fromInteger (Integer
    2 Integer -> Integer -> Integer
    forall a. Num a => a -> a -> a
    * Integer
    f Integer -> Integer -> Integer
    forall a. Num a => a -> a -> a
    * (Integer
    2 Integer -> Integer -> Integer
    forall a. Num a => a -> a -> a
    * Integer
    f Integer -> Integer -> Integer
    forall a. Num a => a -> a -> a
    - Integer
    1))) [Rational]
    cs [Integer
    1..]

    -- | Infinite sequence of approximate values of Riemann zeta-function
    -- at odd arguments, starting with @ζ(1)@.
    zetasOdd :: forall a. (Floating a, Ord a) => a -> [a]
    zetasOdd :: forall a. (Floating a, Ord a) => a -> [a]
    zetasOdd a
    eps = (a
    1 a -> a -> a
    forall a. Fractional a => a -> a -> a
    / a
    0) a -> [a] -> [a]
    forall a. a -> [a] -> [a]
    : [a] -> [a]
    forall a. HasCallStack => [a] -> [a]
    tail ([a] -> [a]
    forall a. [a] -> [a]
    skipEvens ([a] -> [a]) -> [a] -> [a]
    forall a b. (a -> b) -> a -> b
    $ a -> a -> [a]
    forall a. (Floating a, Ord a) => a -> a -> [a]
    zetaHurwitz a
    eps a
    1)

    -- | Infinite sequence of approximate (up to given precision)
    -- values of Riemann zeta-function at integer arguments, starting with @ζ(0)@.
    --
    -- >>> take 5 (zetas 1e-14) :: [Double]
    -- [-0.5,Infinity,1.6449340668482264,1.2020569031595942,1.0823232337111381]
    --
    -- Beware to force evaluation of @zetas !! 1@ if the type @a@ does not support infinite values
    -- (for instance, 'Data.Number.Fixed.Fixed').
    --
    zetas :: (Floating a, Ord a) => a -> [a]
    zetas :: forall a. (Floating a, Ord a) => a -> [a]
    zetas a
    eps = a
    e a -> [a] -> [a]
    forall a. a -> [a] -> [a]
    : a
    o a -> [a] -> [a]
    forall a. a -> [a] -> [a]
    : (a -> a -> a) -> [a] -> [a]
    forall a. (a -> a -> a) -> [a] -> [a]
    scanl1 a -> a -> a
    forall {a}. (Ord a, Fractional a) => a -> a -> a
    f ([a] -> [a] -> [a]
    forall a. [a] -> [a] -> [a]
    intertwine [a]
    es [a]
    os)
      where
        a
    e : [a]
    es = (ExactPi -> a) -> [ExactPi] -> [a]
    forall a b. (a -> b) -> [a] -> [b]
    map ((a -> a -> Bool) -> [Rational] -> a
    forall a. Fractional a => (a -> a -> Bool) -> [Rational] -> a
    getRationalLimit (\a
    a a
    b -> a -> a
    forall a. Num a => a -> a
    abs (a
    a a -> a -> a
    forall a. Num a => a -> a -> a
    - a
    b) a -> a -> Bool
    forall a. Ord a => a -> a -> Bool
    < a
    eps) ([Rational] -> a) -> (ExactPi -> [Rational]) -> ExactPi -> a
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . ExactPi -> [Rational]
    rationalApproximations) [ExactPi]
    zetasEven
        a
    o : [a]
    os = a -> [a]
    forall a. (Floating a, Ord a) => a -> [a]
    zetasOdd a
    eps

        -- Cap-and-floor to improve numerical stability:
        -- 0 < zeta(n + 1) - 1 < (zeta(n) - 1) / 2
        f :: a -> a -> a
    f a
    x a
    y = a
    1 a -> a -> a
    forall a. Ord a => a -> a -> a
    `max` (a
    y a -> a -> a
    forall a. Ord a => a -> a -> a
    `min` (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
    1) a -> a -> a
    forall a. Fractional a => a -> a -> a
    / a
    2))