module Math.ExpPairs.Ivic
( zetaOnS
, reverseZetaOnS
, mOnS
, reverseMOnS
, checkAbscissa
, findMinAbscissa
, mBigOnHalf
, reverseMBigOnHalf
) 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
[RationalForm (LinearForm 1 1 (s)) 2]
[Constraint (LinearForm (1) 1 (s)) NonStrict]
| 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 [RationalForm (LinearForm 1 (1) 1) 1] [Constraint (LinearForm 0 (2) (1+2*mu)) NonStrict]
| mu == zetaOnHalf = simulateOptimize (1 % 2)
| otherwise = optRes {optimalValue = negate $ optimalValue optRes} where
optRes = optimize [RationalForm (LinearForm 1 (1) 0) 1] [Constraint (LinearForm 1 0 (mu)) NonStrict, Constraint (LinearForm (1) 1 (1%2)) NonStrict]
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}
numer = LinearForm
(4*s + (8*muS + 2))
(8*s + (8*muS + 6))
(2*s + (4*muS + 1))
denom = LinearForm
(2*muS)
(4*muS*s 2*muS)
(2*muS*s muS)
cons = if s >= 2%3 then [] else [Constraint
(LinearForm (4*s2) (8*s6) (2*s1)) NonStrict
]
x2' = optimize [RationalForm numer denom] cons
x2 = x2' {optimalValue = negate $ optimalValue x2'}
reverseMOnS :: Rational -> RationalInf -> Rational
reverseMOnS prec m = reverseMOnS' from to where
from = 1 % 2
to = 1
reverseMOnS' a b
| b a < prec = c
| optimalValue (mOnS c) > m = reverseMOnS' a c
| otherwise = reverseMOnS' c b
where
c = (numerator a + numerator b) % (denominator a + denominator b)
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 = searchMinAbscissa' from to where
from = 1 % 2 / minimum (map fst xs)
to = 1 % 1
searchMinAbscissa' a b
| b a < prec = b
| checkAbscissa xs c = searchMinAbscissa' a c
| otherwise = searchMinAbscissa' c b
where
c = (numerator a + numerator b) % (denominator a + denominator b)
mBigOnHalf :: Rational -> OptimizeResult
mBigOnHalf a
| a < 4 = simulateOptimize 1
| a < 12 = simulateOptimize $ 1+(a4)/8
| a > 41614060315296730740083860226662 % 2636743270445733804969041895717 = simulateOptimize $ 1 + 32*(a6)/205
| otherwise = if Finite x >= optimalValue optRes
then simulateOptimize x
else optRes where
optRes = optimize [RationalForm (LinearForm 1 1 0) (LinearForm 1 0 0)]
[Constraint (LinearForm (4a) 4 2) NonStrict]
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 = (m1)*205/32 + 6
optRes = optimize [RationalForm (LinearForm 4 4 2) (LinearForm 1 0 0)] [Constraint (LinearForm (1m) 1 0) NonStrict]