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

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

module Math.NumberTheory.Zeta.Dirichlet
  ( betas
  , betasEven
  , betasOdd
  ) where

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

import Math.NumberTheory.Recurrences  (euler, factorial)
import Math.NumberTheory.Zeta.Hurwitz (zetaHurwitz)
import Math.NumberTheory.Zeta.Utils   (skipOdds)

-- | Infinite sequence of exact values of Dirichlet beta-function at odd arguments, starting with @β(1)@.
--
-- >>> import Data.ExactPi
-- >>> approximateValue (betasOdd !! 25) :: Double
-- 0.9999999999999987
-- >>> import Data.Number.Fixed
-- >>> approximateValue (betasOdd !! 25) :: Fixed Prec50
-- 0.99999999999999999999999960726927497384196726751694
betasOdd :: Infinite ExactPi
betasOdd :: Infinite ExactPi
betasOdd = (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
1, Integer
3)....) (Infinite Rational -> Infinite ExactPi)
-> Infinite Rational -> Infinite ExactPi
forall a b. (a -> b) -> a -> b
$ (Rational -> Integer -> Integer -> Integer -> Rational)
-> Infinite Rational
-> Infinite Integer
-> Infinite Integer
-> Infinite Integer
-> Infinite Rational
forall a b c d e.
(a -> b -> c -> d -> e)
-> Infinite a
-> Infinite b
-> Infinite c
-> Infinite d
-> Infinite e
Inf.zipWith4
                                     (\Rational
sgn Integer
denom Integer
eul Integer
twos -> Rational
sgn Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* (Integer
eul Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
% (Integer
twos Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Integer
denom)))
                                     (NonEmpty Rational -> Infinite Rational
forall a. NonEmpty a -> Infinite a
Inf.cycle (Rational
1 Rational -> [Rational] -> NonEmpty Rational
forall a. a -> [a] -> NonEmpty a
:| [-Rational
1]))
                                     (Infinite Integer -> Infinite Integer
forall a. Infinite a -> Infinite a
skipOdds Infinite Integer
forall a. (Num a, Enum a) => Infinite a
factorial)
                                     (Infinite Integer -> Infinite Integer
forall a. Infinite a -> Infinite a
skipOdds Infinite Integer
forall a. Integral a => Infinite a
euler)
                                     ((Integer -> Integer) -> Integer -> Infinite Integer
forall a. (a -> a) -> a -> Infinite a
Inf.iterate (Integer
4 *) Integer
4)

-- | Infinite sequence of approximate values of the Dirichlet @β@ function at
-- positive even integer arguments, starting with @β(0)@.
betasEven :: forall a. (Floating a, Ord a) => a -> Infinite a
betasEven :: forall a. (Floating a, Ord a) => a -> Infinite a
betasEven a
eps = (a
1 a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
2) a -> Infinite a -> Infinite a
forall a. a -> Infinite a -> Infinite a
:< Infinite a
hurwitz
  where
    hurwitz :: Infinite a
    hurwitz :: Infinite a
hurwitz =
        (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
quarter a
threeQuarters a
four ->
            (a
quarter a -> a -> a
forall a. Num a => a -> a -> a
- a
threeQuarters) a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
four)
        (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
. Infinite a -> Infinite a
forall a. Infinite a -> Infinite a
skipOdds (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
0.25)
        (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
. Infinite a -> Infinite a
forall a. Infinite a -> Infinite a
skipOdds (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
0.75)
        ((a -> a) -> a -> Infinite a
forall a. (a -> a) -> a -> Infinite a
Inf.iterate (a
16 *) a
16)

-- | Infinite sequence of approximate (up to given precision)
-- values of Dirichlet beta-function at integer arguments, starting with @β(0)@.
--
-- >>> take 5 (betas 1e-14) :: [Double]
-- [0.5,0.7853981633974483,0.9159655941772189,0.9689461462593694,0.9889445517411051]
betas :: (Floating a, Ord a) => a -> Infinite a
betas :: forall a. (Floating a, Ord a) => a -> Infinite a
betas 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 = a -> Infinite a
forall a. (Floating a, Ord a) => a -> Infinite a
betasEven a
eps
    a
o :< Infinite a
os = (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
betasOdd

    -- Cap-and-floor to improve numerical stability:
    -- 1 > beta(n + 1) - 1 > (beta(n) - 1) / 2
    -- A similar method is used in @Math.NumberTheory.Zeta.Riemann.zetas@.
    f :: a -> a -> a
f a
x a
y = a
1 a -> a -> a
forall a. Ord a => a -> a -> a
`min` (a
y a -> a -> a
forall a. Ord a => a -> a -> a
`max` (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))