{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE FlexibleInstances #-}
module Quantum.Synthesis.QuadraticEquation (
Quadratic (..)
) where
import Data.Number.FixedPrec
import Quantum.Synthesis.Ring
import Quantum.Synthesis.ToReal
class Quadratic t a where
quadratic :: t -> t -> t -> Maybe (a, a)
int_quadratic :: (Fractional t, Floor t, Ord t) => t -> t -> Maybe (Integer, Integer)
int_quadratic b c
| radix < 0 = Nothing
| otherwise = Just (t0, t1)
where
radix = b^2/4 - c
tm = -b / 2
rootradix' = intsqrt (floor_of radix)
t1' = floor_of tm + rootradix'
t1
| is_solution1 (t1'+2) = t1'+2
| is_solution1 (t1'+1) = t1'+1
| otherwise = t1'
t0' = ceiling_of tm - rootradix'
t0
| is_solution0 (t0'-2) = t0'-2
| is_solution0 (t0'-1) = t0'-1
| otherwise = t0'
is_solution1 x = f x' >= 0 && (f (x'-1) < 0 || x'-1 < tm) where
x' = fromInteger x
is_solution0 x = f x' >= 0 && (f (x'+1) < 0 || x'-1 > tm) where
x' = fromInteger x
f x = x^2 + b*x + c
quadratic_fixedprec :: (Fractional t, Floor t, Ord t, Precision e) => t -> t -> t -> Maybe (FixedPrec e, FixedPrec e)
quadratic_fixedprec a b c
| False = Just (r, r)
| otherwise = do
(x0, x1) <- int_quadratic b' c'
return (fromInteger x0 / prec, fromInteger x1 / prec)
where
r = 0
d = getprec r
prec = 10^d
prec' = 10^d
b' = prec' * b/a
c' = prec'^2 * c/a
q = int_quadratic b' c'
instance (Fractional t, Floor t, Ord t, Precision e) => Quadratic t (FixedPrec e) where
quadratic = quadratic_fixedprec
instance (ToReal t) => Quadratic t Double where
quadratic a' b' c'
| radix < 0 = Nothing
| b >= 0 = Just (t1, t2)
| otherwise = Just (t1', t2')
where
radix = b^2 - 4*a*c
s1 = -b - sqrt radix
s2 = -b + sqrt radix
t1 = s1 / (2*a)
t2 = (2*c) / s1
t1' = (2*c) / s2
t2' = s2 / (2*a)
a = to_real a'
b = to_real b'
c = to_real c'