module Math.ExpPairs.Kratzel
( TauabTheorem (..)
, tauab
, TauabcTheorem (..)
, tauabc
) where
import Control.Arrow
import Data.Ratio ((%))
import Data.Ord (comparing)
import Data.List (minimumBy)
import Text.PrettyPrint.Leijen
import Math.ExpPairs
data TauabTheorem
= Kr511a
| Kr511b
| Kr512a
| Kr512b
deriving (Eq, Ord, Enum, Bounded, Show)
instance Pretty TauabTheorem where
pretty = text . show
divideResult :: Real a => a -> (b, OptimizeResult) -> (b, OptimizeResult)
divideResult d = second (\o -> o {optimalValue = optimalValue o / Finite (toRational d)})
tauab :: Integer -> Integer -> (TauabTheorem, OptimizeResult)
tauab a' b'
| d /= 1 = divideResult d $ tauab (a'`div` d) (b' `div` d) where
d = gcd a' b'
tauab a' b' = minimumBy (comparing (optimalValue . snd)) [kr511a, kr511b, kr512a, kr512b] where
a = toRational a'
b = toRational b'
kr511a = (Kr511a, optimize
[RationalForm (LinearForm 2 2 (1)) (LinearForm 0 0 (a+b))]
[Constraint (LinearForm (2*b) (2*a) (a)) NonStrict])
kr511b = (Kr511b, optimize
[RationalForm (LinearForm 1 0 0) (LinearForm b (a) a)]
[Constraint (LinearForm (2*b) (2*a) a) Strict])
kr512a = (Kr512a, simulateOptimize r) where
r = if 11*a >= 8*b then 19/29/(a+b) else 1%1
kr512b = if 11*a >= 8*b then kr512a else (Kr512b, optimize
[
RationalForm (LinearForm (11) 8 (4)) (LinearForm (29*b) (29*a) (4*b20*a))
]
[
Constraint (LinearForm (2*b) (2*a) (a)) NonStrict,
Constraint (LinearForm (29) 0 4) Strict,
Constraint (LinearForm 29 29 (24)) Strict
])
data TauabcTheorem
= Kolesnik
| Kr61
| Kr62
| Kr63
| Kr64
| Kr65
| Kr66
| Tauab TauabTheorem
deriving (Eq, Ord, Show)
instance Pretty TauabcTheorem where
pretty (Tauab t) = pretty t
pretty t = pretty (show t)
tauabc :: Integer -> Integer -> Integer -> (TauabcTheorem, OptimizeResult)
tauabc a' b' c'
| d /= 1 = divideResult d $ tauabc (a'`div` d) (b' `div` d) (c' `div` d) where
d = gcd (gcd a' b') c'
tauabc 1 1 1 = (Kolesnik, simulateOptimize $ 43%96)
tauabc a' b' c' = minimumBy (comparing (optimalValue . snd)) [kr61, kr62, kr63, kr64, kr65, kr66] where
a = toRational a'
b = toRational b'
c = toRational c'
kr61
| c<a+b = (Kr61, simulateOptimize $ 2/(a+b+c))
| optimalValue optRes < Finite (recip c) = (Kr61, simulateOptimize $ 1/c)
| otherwise = (Tauab th, optRes)
where
(th, optRes) = tauab a' b'
kr62 = (Kr62, optimize
[RationalForm (LinearForm 2 2 0) (LinearForm 0 0 (a+b+c))]
[
Constraint (LinearForm (bc) a 0) NonStrict,
Constraint (LinearForm (2*c) (2*c) (a+b+c)) NonStrict
])
kr63 = (Kr63, optimize
[RationalForm (LinearForm 4 2 3) (LinearForm (2*(a+b+c)) 0 (3*(a+b+c)))]
[Constraint (LinearForm (2*(abc)) (2*a) (2*abc)) NonStrict])
kr64 = (Kr64, simulateOptimize r) where
r = recip (a+b+c) * minimum ((a+b+c):[24*(k1)%(3*2^k4) | k<-[1..maxk], (3*2^k2*k4)%1 * a >= 2 * (b+c), (3*2^k8)%1 * (a+b) >= (3*2^k4*k+4)%1 * c])
maxk = 4 `max` floor (logBase 2 (fromRational $ b+c) :: Double)
kr65 = (Kr65, simulateOptimize r) where
r = if 7*a>=2*(b+c) && 4*(a+b)>=5*c then 3%2/(a+b+c) else 1%1
kr66 = (Kr66, simulateOptimize r) where
r = if 18*a>=7*(b+c) && 2*(a+b)>=3*c then 25%17/(a+b+c) else 1%1