content: -- |
  -- Module:      Math.NumberTheory.Zeta.DirichletTests
  -- Copyright:   (c) 2018 Alexandre Rodrigues Baldé
  -- Licence:     MIT
  -- Maintainer:  Alexandre Rodrigues Baldé <alexandrer_b@outlook.com>
  --
  -- Tests for Math.NumberTheory.Zeta.Dirichlet
  --
  
  {-# OPTIONS_GHC -fno-warn-type-defaults #-}
  
  module Math.NumberTheory.Zeta.DirichletTests
    ( testSuite
    ) where
  
  import Data.ExactPi                (approximateValue)
  import Data.List.Infinite (Infinite(..))
  import qualified Data.List.Infinite as Inf
  
  import Test.Tasty
  import Test.Tasty.HUnit            (Assertion, assertEqual, testCase)
  
  import Math.NumberTheory.Zeta      (betas, betasOdd)
  import Math.NumberTheory.TestUtils
  
  epsilon :: Double
  epsilon = 1e-14
  
  betas' :: Infinite Double
  betas' = betas epsilon
  
  betasOddSpecialCase1 :: Assertion
  betasOddSpecialCase1
    = assertEqualUpToEps "beta(1) = pi/4" epsilon
      (approximateValue $ Inf.head betasOdd)
      (pi / 4)
  
  betasOddSpecialCase2 :: Assertion
  betasOddSpecialCase2
    = assertEqualUpToEps "beta(3) = pi^3/32" epsilon
      (approximateValue $ betasOdd Inf.!! 1)
      (pi^3 / 32)
  
  betasOddSpecialCase3 :: Assertion
  betasOddSpecialCase3
    = assertEqualUpToEps "beta(5) = 5*pi^5/1536" epsilon
      (approximateValue $ betasOdd Inf.!! 2)
      ((5 * pi^5) / 1536)
  
  betasOddProperty1 :: Positive Int -> Bool
  betasOddProperty1 (Positive m)
    =  betaM < 1
    || betaM < betaM1
    where
      betaM  = approximateValue (betasOdd Inf.!! fromIntegral m)
      betaM1 = approximateValue (betasOdd Inf.!! (fromIntegral m + 1))
  
  betasOddProperty2 :: Word -> Bool
  betasOddProperty2 m
    = abs (betaM - betaM') < epsilon
    where
      betaM  = approximateValue (betasOdd Inf.!! m)
      betaM' = betas' Inf.!! ((2 * m) + 1)
  
  betasSpecialCase1 :: Assertion
  betasSpecialCase1
    = assertEqual "beta(0) = 1/2"
      (Inf.head betas')
      (1 / 2)
  
  betasSpecialCase2 :: Assertion
  betasSpecialCase2
    = assertEqualUpToEps "beta(2) = 0.9159655" epsilon
      (betas' Inf.!! 2)
      0.9159655941772190150546035149323841107
  
  betasSpecialCase3 :: Assertion
  betasSpecialCase3
    = assertEqualUpToEps "beta(4) = 0.9889445" epsilon
      (betas' Inf.!! 4)
      0.9889445517411053361084226332283778213
  
  betasProperty1 :: Positive Int -> Bool
  betasProperty1 (Positive m)
    =  betaM <= betaM1
    && betaM1 <= 1
    where
      betaM  = betas' Inf.!! fromIntegral m
      betaM1 = betas' Inf.!! (fromIntegral m + 1)
  
  betasProperty2 :: NonNegative Int -> NonNegative Int -> Bool
  betasProperty2 (NonNegative e1) (NonNegative e2)
    = maximum (Inf.take 35 $ Inf.drop 2 $ Inf.zipWith ((abs .) . (-)) (betas eps1) (betas eps2)) <= eps1 + eps2
    where
      eps1, eps2 :: Double
      eps1 = max ((1.0 / 2) ^ e1) ((1.0 / 2) ^ 53)
      eps2 = max ((1.0 / 2) ^ e2) ((1.0 / 2) ^ 53)
  
  
  testSuite :: TestTree
  testSuite = testGroup "Beta"
    [ testGroup "betasOdd"
      [ testCase "beta(1)"                                        betasOddSpecialCase1
      , testCase "beta(3)"                                        betasOddSpecialCase2
      , testCase "beta(5)"                                        betasOddSpecialCase3
      , testSmallAndQuick "beta(2n-1) < beta(2n+1)"               betasOddProperty1
      , testSmallAndQuick "betasOdd matches betas"                betasOddProperty2
      ]
    , testGroup "betas"
      [ testCase "beta(0)"                                        betasSpecialCase1
      , testCase "beta(2)"                                        betasSpecialCase2
      , testCase "beta(4)"                                        betasSpecialCase3
      , testSmallAndQuick "beta(n) < beta(n+1)"                   betasProperty1
      , testSmallAndQuick "precision"                             betasProperty2
      ]
    ]
  