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

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

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

import Data.ExactPi
import Data.List.Infinite (Infinite(..), (...), (....))
import qualified Data.List.Infinite as Inf
import Data.Ratio                     ((%))

import Math.NumberTheory.Recurrences  (bernoulli)
import Math.NumberTheory.Zeta.Hurwitz (zetaHurwitz)
import Math.NumberTheory.Zeta.Utils   (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 :: Infinite ExactPi
zetasEven :: Infinite ExactPi
zetasEven = (Integer -> Rational -> ExactPi)
-> Infinite Integer -> Infinite Rational -> Infinite ExactPi
forall a b c.
(a -> b -> c) -> Infinite a -> Infinite b -> Infinite c
Inf.zipWith Integer -> Rational -> ExactPi
Exact ((Integer
0, Integer
2)....) (Infinite Rational -> Infinite ExactPi)
-> Infinite Rational -> Infinite ExactPi
forall a b. (a -> b) -> a -> b
$ (Rational -> Rational -> Rational)
-> Infinite Rational -> Infinite Rational -> Infinite Rational
forall a b c.
(a -> b -> c) -> Infinite a -> Infinite b -> Infinite c
Inf.zipWith Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
(*) (Infinite Rational -> Infinite Rational
forall a. Infinite a -> Infinite a
skipOdds Infinite Rational
forall a. Integral a => Infinite (Ratio a)
bernoulli) Infinite Rational
cs
  where
    cs :: Infinite Rational
    cs :: Infinite Rational
cs = (- Integer
1 Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
% Integer
2) Rational -> Infinite Rational -> Infinite Rational
forall a. a -> Infinite a -> Infinite a
:< (Rational -> Integer -> Rational)
-> Infinite Rational -> Infinite Integer -> Infinite Rational
forall a b c.
(a -> b -> c) -> Infinite a -> Infinite b -> Infinite c
Inf.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))) Infinite 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 -> Infinite a
zetasOdd :: forall a. (Floating a, Ord a) => a -> Infinite a
zetasOdd a
eps = (a
1 a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
0) a -> Infinite a -> Infinite a
forall a. a -> Infinite a -> Infinite a
:< Infinite a -> Infinite a
forall a. Infinite a -> Infinite a
Inf.tail (Infinite a -> Infinite a
forall a. Infinite a -> Infinite a
skipEvens (Infinite a -> Infinite a) -> Infinite a -> Infinite a
forall a b. (a -> b) -> a -> b
$ a -> a -> Infinite a
forall a. (Floating a, Ord a) => a -> a -> Infinite 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 -> Infinite a
zetas :: forall a. (Floating a, Ord a) => a -> Infinite a
zetas a
eps = a
e a -> Infinite a -> Infinite a
forall a. a -> Infinite a -> Infinite a
:< a
o a -> Infinite a -> Infinite a
forall a. a -> Infinite a -> Infinite a
:< (a -> a -> a) -> Infinite a -> Infinite a
forall a. (a -> a -> a) -> Infinite a -> Infinite a
Inf.scanl1 a -> a -> a
forall {a}. (Ord a, Fractional a) => a -> a -> a
f (Infinite a -> Infinite a -> Infinite a
forall a. Infinite a -> Infinite a -> Infinite a
Inf.interleave Infinite a
es Infinite a
os)
  where
    a
e :< Infinite a
es = (ExactPi -> a) -> Infinite ExactPi -> Infinite a
forall a b. (a -> b) -> Infinite a -> Infinite b
Inf.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) Infinite ExactPi
zetasEven
    a
o :< Infinite a
os = a -> Infinite a
forall a. (Floating a, Ord a) => a -> Infinite 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))