module Math.ExpPairs.Ivic (zetaOnS, reverseZetaOnS, mOnS) where
import Data.Ratio
import Data.List
import Data.Ord
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]
| 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 m = reverseMOnS' from to where
from = 1 % 2
to = 1 % 1
reverseMOnS' a b
| ba < 1%1000000 = a
| optimalValue (mOnS ((a+b)/2)) > m = reverseMOnS' a ((a+b)/2)
| otherwise = reverseMOnS' ((a+b)/2) 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
searchMinAbscissa :: [(Rational, Rational)] -> Rational
searchMinAbscissa xs = searchMinAbscissa' from to where
from = 1 % 2 / minimum (map fst xs)
to = 1 % 1
searchMinAbscissa' a b
| ba < 1%1000000 = a
| checkAbscissa xs ((a+b)/2) = searchMinAbscissa' a ((a+b)/2)
| otherwise = searchMinAbscissa' ((a+b)/2) b
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 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]
f l lambda = (llambda)*32/205
+ ((toRational . optimalValue . mBigOnHalf) (4*l/(4lambda)) 1) * (4lambda) /4
heckeZetaByHalf a = 1 xt where
d = toRational 12.571624917200547
ia 2 = 0
ia 3 = 1%4
ia a
| a<=12 = 32%205 * ((1 + 4 / d) * a 4) + (toRational (optimalValue (mBigOnHalf d)) 1) * a / d
| a<=15 = 32%205 * a + toRational (optimalValue (mBigOnHalf a)) 1
| otherwise = 32%205 * (2 * a 6)
xt = 1%2 / (1 + ia a)
bestLambda l = (\x -> (x, fromRational $ f l x)) (bestLambda' 0 (3999%1000)) where
bestLambda' a b
| ba < 1%1000000000 = a
| otherwise = if mx1 > mx2 then bestLambda' x1 b else bestLambda' a x2 where
x1 = (2*a+b)/3
x2 = (a+2*b)/3
ma = f l a
mb = f l b
mx1 = f l x1
mx2 = f l x2
difur a = a * m' m + (77%205) where
h = 1%(10^100)
m = toRational $ optimalValue $ mBigOnHalf a
mh = toRational $ optimalValue $ mBigOnHalf (a+h)
m' = (mhm) / h
solveD a b
| ba < 1%(10^6) = a
| otherwise = if difur c > 0 then solveD a c else solveD c b where
c = (a+b)/2