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 (1s)
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/(34*s)
| s<=11%14 = 10/(78*s)
| s<=13%15 = 34/(1516*s)
| s<=57%62 = 98/(3132*s)
| otherwise = 5/(1s)
mOnS :: Rational -> OptimizeResult
mOnS s
| s < 1%2 = simulateOptimize 0
| s < 5%8 = simulateOptimize $ 4/(34*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 = (44*s)/(1+2*s)
beta1 = 12/(1+2*s)
x1 = optRes {optimalValue = Finite $ (1alpha1)/muS beta1}
t = scaleLF s (L 4 + 2) 1 + K 2 L 2
numer = t scaleLF (4 * (1s)) (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+(a4)/8
| a > 41614060315296730740083860226662 % 2636743270445733804969041895717 = 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 + 32*(a6)/205
reverseMBigOnHalf :: Rational -> OptimizeResult
reverseMBigOnHalf m
| m <= 2 = simulateOptimize $ (m1)*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]