module Math.ExpPairs.Ivic
( zetaOnS
, reverseZetaOnS
, mOnS
, reverseMOnS
, checkAbscissa
, findMinAbscissa
, mBigOnHalf
, reverseMBigOnHalf
, kolpakova2011
) where
import Data.Ratio
import Data.List (minimumBy)
import Data.Ord (comparing)
import Math.ExpPairs
zetaOnS :: Rational -> OptimizeResult
zetaOnS s
| s >= 1 = simulateOptimize 0
| s >= 1%2 = optimize
[K 1 + L 1 - M s :/: 2]
[L 1 >=. K 1 + M s]
| otherwise = optRes {optimalValue = r} where
optRes = zetaOnS (1-s)
r = Finite (1%2 - s) + optimalValue optRes
zetaOnHalf :: Rational
zetaOnHalf = 32%205
reverseZetaOnS :: Rational -> OptimizeResult
reverseZetaOnS mu
| mu >= 1%2 = simulateOptimize 0
| mu > zetaOnHalf = optimize [K 1 - L 1 + M 1 :/: 1] [M (1 + 2 * mu) >=. L 2]
| mu == zetaOnHalf = simulateOptimize (1 % 2)
| otherwise = optRes {optimalValue = negate $ optimalValue optRes} where
optRes = optimize [K 1 - L 1 :/: 1] [K 1 >=. M mu, L 2 >=. K 2 + 1]
lemma82_f :: Rational -> Rational
lemma82_f s
| s < 1%2 = undefined
| s<= 2%3 = 2/(3-4*s)
| s<=11%14 = 10/(7-8*s)
| s<=13%15 = 34/(15-16*s)
| s<=57%62 = 98/(31-32*s)
| otherwise = 5/(1-s)
mOnS :: Rational -> OptimizeResult
mOnS s
| s < 1%2 = simulateOptimize 0
| s < 5%8 = simulateOptimize $ 4/(3-4*s)
| s>= 1 = simulateOptimize' InfPlus
| otherwise = minimumBy (comparing optimalValue) [x1, x2, simulateOptimize (lemma82_f s * 2)] where
optRes = zetaOnS s
muS = toRational $ optimalValue optRes
alpha1 = (4-4*s)/(1+2*s)
beta1 = -12/(1+2*s)
x1 = optRes {optimalValue = Finite $ (1-alpha1)/muS - beta1}
t = scaleLF s (L 4 + 2) - 1 + K 2 - L 2
numer = t - scaleLF (4 * (1-s)) (K 1 + L 1) + scaleLF (4 * muS) (K 2 + L 2 + 1)
denom = scaleLF muS t
cons = if s >= 2%3 then [] else [scaleLF s (K 4 + L 8 + 2) >=. K 2 + L 6 + 1]
x2' = optimize [- numer :/: denom] cons
x2 = x2' {optimalValue = negate $ optimalValue x2'}
data Choice = Least | Median | Greatest
binarySearch :: (Rational -> Bool) -> Choice -> Rational -> Rational -> Rational -> Rational
binarySearch predicate choice precision = go
where
go a b
| b - a < precision = case choice of
Least -> a
Median -> c
Greatest -> b
| predicate c = go a c
| otherwise = go c b
where
c = (numerator a + numerator b) % (denominator a + denominator b)
mOnSTwoThird :: RationalInf
mOnSTwoThird = optimalValue $ mOnS $ 2 % 3
reverseMOnS :: Rational -> RationalInf -> Rational
reverseMOnS _ InfPlus = 1
reverseMOnS _ (Finite m)
| m <= 4 = 1 % 2
| m <= 8 = 3 % 4 - recip m
reverseMOnS prec m
| m < mOnSTwoThird = go (5 % 8) (2 % 3)
| otherwise = go (2 % 3) 1
where
go = binarySearch (\c -> optimalValue (mOnS c) > m) Greatest prec
kolpakova2011 :: Integer -> Double
kolpakova2011 k = 1 - 1/3 * 2**(2/3) * (4.45 * fromInteger k)**(-2/3)
checkAbscissa :: [(Rational, Rational)] -> Rational -> Bool
checkAbscissa xs s = sum rs < Finite 1 where
qs = map (\(n,m) -> optimalValue (mOnS (n*s)) / Finite m) xs
rs = map (\q -> 1/q) qs
findMinAbscissa :: Rational -> [(Rational, Rational)] -> Rational
findMinAbscissa prec xs = binarySearch (checkAbscissa xs) Greatest prec (1 % 2 / minimum (map fst xs)) 1
mBigOnHalf :: Rational -> OptimizeResult
mBigOnHalf a
| a < 4 = simulateOptimize 1
| a < 12 = simulateOptimize $ 1+(a-4)/8
| a > 16645467 / 972266 = simulateOptimize $ 1 + (a - 6) * zetaOnHalf
| otherwise = if Finite x >= optimalValue optRes
then simulateOptimize x
else optRes where
optRes = optimize [K 1 + L 1 :/: K 1]
[K (4 - a) + L 4 + 2 >=. 0]
x = 1 + 13*(a-6)/84
reverseMBigOnHalf :: Rational -> OptimizeResult
reverseMBigOnHalf m
| m <= 2 = simulateOptimize $ (m-1)*8 + 4
| otherwise = if Finite a <= optimalValue optRes
then simulateOptimize a
else optRes where
a = (m - 1) / zetaOnHalf + 6
optRes = optimize [K 4 + L 4 + 2 :/: K 1] [K (1 - m) + L 1 >=. 0]