1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
|
-- |
-- 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 = Inf.zipWith Exact ((1, 3)....) $ Inf.zipWith4
(\sgn denom eul twos -> sgn * (eul % (twos * denom)))
(Inf.cycle (1 :| [-1]))
(skipOdds factorial)
(skipOdds euler)
(Inf.iterate (4 *) 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 eps = (1 / 2) :< hurwitz
where
hurwitz :: Infinite a
hurwitz =
Inf.zipWith3 (\quarter threeQuarters four ->
(quarter - threeQuarters) / four)
(Inf.tail . skipOdds $ zetaHurwitz eps 0.25)
(Inf.tail . skipOdds $ zetaHurwitz eps 0.75)
(Inf.iterate (16 *) 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 eps = e :< o :< Inf.scanl1 f (Inf.interleave es os)
where
e :< es = betasEven eps
o :< os = Inf.map (getRationalLimit (\a b -> abs (a - b) < eps) . rationalApproximations) 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 x y = 1 `min` (y `max` (1 + (x - 1) / 2))
|
