```{-# LANGUAGE TypeSynonymInstances #-}
{-# LANGUAGE FlexibleInstances #-}

-- | Define various analytic functions for 'FDouble' and 'QDouble'.
module Algorithms.QLS.RealFunc where

import Data.List(mapAccumL)
import Quipper
import Algorithms.QLS.CircLiftingImport
import Algorithms.QLS.QSignedInt
import Algorithms.QLS.QDouble
import Algorithms.QLS.Utils

-- * Some analytic functions on 'FDouble'.

-- | Approximation of the sine function using Taylor series.
build_circuit
approx_sin :: FDouble -> FDouble
approx_sin x = let x2 = x * x in
let x3 = x2 * x in
let x4 = x2 * x2 in
let x5 = x4 * x in
let x7 = x2 * x5 in
let x9 = x2 * x7 in
let x11 = x2 * x9 in
x - (x3 / 6.0)
+ (x5 / 120.0)
- (x7 / 5040.0)
+ (x9 / 362880.0)
- (x11 / 39916800.0)

-- | Implementation of the sine function valid on the whole domain of
-- 'FDouble'.
build_circuit
local_sin :: FDouble -> FDouble
local_sin x = let n = fromIntegral \$ floor (x/(2.0*local_pi)) in
let y = x - 2.0*local_pi*n in
if (y < local_pi/2.0) then approx_sin y
else if (y > 3.0*local_pi/2.0) then approx_sin (y - 2.0*local_pi)
else approx_sin (local_pi - y)

-- | Implementation of the cosine function valid on the whole domain
-- of 'FDouble'.
build_circuit
local_cos :: FDouble -> FDouble
local_cos x = local_sin (x + local_pi/2.0)

-- listAngle :: [Double]
listAngle = snd \$ mapAccumL (\a x -> (a+1, x / (2 ** a))) 3.0 (replicate after_radix_length pi)

-- template_listAngle :: Circ [QDouble]
template_listAngle = mapM (\x -> qinit \$ fdouble x) listAngle

-- listCos :: [Double]
listCos = map cos listAngle

-- template_listCos :: Circ [QDouble]
template_listCos = mapM (\x -> qinit \$ fdouble x) listCos

-- listSin :: [Double]
listSin = map sin listAngle

-- template_listSin :: Circ [QDouble]
template_listSin = mapM (\x -> qinit \$ fdouble x) listCos

-- list_values :: [(FDouble,FDouble,FDouble)]
list_values = map (\(x,y,z) -> (fdouble x, fdouble y, fdouble z)) \$
zip3 listAngle listCos listSin

template_list_values = mapM (\(x,y,z) -> do
x' <- qinit \$ fdouble x
y' <- qinit \$ fdouble y
z' <- qinit \$ fdouble z
return (x',y',z')) \$ zip3 listAngle listCos listSin

-- | Auxiliary function for 'local_sqrt'.
build_circuit
approx_sqrt :: Int -> FDouble -> FDouble
approx_sqrt n x = case n of
0 -> x
n -> let s = approx_sqrt (paramPred n) x in (s + x/s)/2.0

-- | Approximation of the square root using iterative means.
build_circuit
local_sqrt :: FDouble -> FDouble
local_sqrt x = approx_sqrt paramTen x

-- | The function 'magnitude' defined for 'FDouble'.
build_circuit
local_mag :: FDouble -> FDouble -> FDouble
local_mag x y = local_sqrt (x * x + y * y)

-- | Apply the matrix
--
-- >  ( a b )
-- >  ( c d )
--
-- to the column vector (/x/,/y/).
build_circuit
rotate :: FDouble -> FDouble -> FDouble -> FDouble -> FDouble -> FDouble -> (FDouble,FDouble)
rotate a b c d x y = (a * x + b * y, c * x + d * y)
where
-- To help the GHC 8.0 typechecker
dummy = id_fdouble x

-- | Auxiliary function for 'approx_atan2'.
build_circuit
approx_atan2_aux :: FDouble -> FDouble -> (FDouble,FDouble, FDouble) -> (FDouble, FDouble, FDouble)
-> (FDouble,FDouble,FDouble)
approx_atan2_aux x y (angle, x', y') (r, cn, sn) =
let (a,(b,c)) = if (y' > y) then (angle - r, rotate cn sn (-sn) cn x' y')
else (angle + r, rotate cn (-sn) sn cn x' y')
in (a,b,c)
where
-- To help the GHC 8.0 typechecker
dummy_r = id_fdouble r

-- | Definition of 'atan2' using a CORDIC method.  Assume (/x/,/y/) is
-- in first quadrant and that /x/ > /y/.
build_circuit
approx_atan2 :: FDouble -> FDouble -> FDouble
approx_atan2 y x =
let list = list_values in
let (a,_,_) = foldl (approx_atan2_aux x y) (0.0, local_mag x y, 0.0) list in a

-- | Definition of 'atan2' using a CORDIC method. /x/ and /y/ can be any 'FDouble'.
build_circuit
local_atan2 :: FDouble -> FDouble -> FDouble
local_atan2 y' x' =
let (x,y,(pad,sign)) = if      (x' >= 0.0 && y' >= 0.0) then ( x',  y', (0.0, 1.0))
else if (x' >= 0.0 && y' <  0.0) then ( x', -y', (0.0, -1.0))
else if (x' <  0.0 && y' <  0.0) then (-x', -y', (-local_pi, 1.0))
else                                  (-x',  y', (local_pi,  -1.0))
in
let angle = if (x > y) then approx_atan2 y x
else            local_pi/2.0 - approx_atan2 x y
in sign * angle + pad

-- | The function 'mkPolar' defined for 'FDouble'.
build_circuit
local_mkPolar :: FDouble -> FDouble -> (FDouble,FDouble)
local_mkPolar p t = (p * local_cos t, p * local_sin t)

instance Floating FDouble where
pi    = fromRational \$ toRational pi
sin x = local_sin x
cos x = local_cos x

sinh x = (exp x - exp (-x)) / 2
cosh x = (exp x - exp (-x)) / 2

asinh x = error "asinh not defined for FDouble"
acosh x = error "acosh not defined for FDouble"
atanh x = error "atanh not defined for FDouble"

exp x = undefined
log x = undefined

asin x = error "asin not defined for FDouble"
acos x = error "acos not defined for FDouble"
atan x = atan2 x 1

instance RealFloat FDouble where
floatDigits x = xdouble_length x
floatRange _ = (0,0)
decodeFloat (XDouble k n) = (integer_of_fsint n, -k)
encodeFloat x k = XDouble (-k) (fromInteger x)
isNaN _ = False
isInfinite _ = False
isDenormalized _ = False
isNegativeZero _ = False
isIEEE _ = False
atan2 y x = local_atan2 y x

-- | A type class for quantum floating-point numbers.
class QFloating a where
-- | Quantum implementation of the sine function.
q_sin :: a -> Circ a
-- | Quantum implementation of the cosine function.
q_cos :: a -> Circ a

instance QFloating QDouble where
q_sin x = (unpack template_local_sin) x
q_cos x = (unpack template_local_cos) x

-- | Quantum implementation of 'atan2' on 'QDouble'.
q_atan2 :: QDouble -> QDouble -> Circ QDouble
q_atan2 = unpack template_local_atan2

-- | Quantum implementation of 'magnitude' on 'QDouble'.
q_magnitude :: (QDouble, QDouble) -> Circ QDouble
q_magnitude = Prelude.uncurry \$ unpack template_local_mag

-- | Quantum implementation of 'mkPolar' on 'QDouble'.
q_mkPolar :: QDouble -> QDouble -> Circ (QDouble,QDouble)
q_mkPolar = unpack template_local_mkPolar

-- | Quantum implementation of 'Re' on 'QDouble': a quantum complex is
-- a pair of two 'QDouble's.
q_Re :: (QDouble,QDouble) -> Circ QDouble
q_Re (x,y) = return x

-- | Quantum implementation of 'Im' on 'QDouble': a quantum complex is
-- a pair of two 'QDouble's.
q_Im :: (QDouble,QDouble) -> Circ QDouble
q_Im (x,y) = return y

my_test_fdouble = do
for 0 37 1 \$ \i -> do
let x = fromIntegral i
let a1 = fromRational \$ toRational (sin(x * pi/37))
let a2 = fromRational \$ toRational (cos(x * pi/37))
let z1 = local_atan2 a1 a2
let z2 = fromRational \$ toRational \$ atan2  (sin(x * pi/37)) (cos(x * pi/37))
putStrLn \$ show_fdouble \$ abs (z1 - z2)

-- * Template subroutines of analytic functions.

-- | Template version of 'sin'.
template_sin :: Circ (QDouble -> Circ QDouble)
template_sin = return \$ \x -> box "sin" q_sin x

-- | Template version of 'cos'.
template_cos :: Circ (QDouble -> Circ QDouble)
template_cos = return \$ \x -> box "cos" q_cos x

-- | Template version of 'atan2'.
template_atan2 :: Circ (QDouble -> Circ (QDouble -> Circ QDouble))
template_atan2 = return \$ \x -> return \$ \y -> box "atan" (uncurry q_atan2) (x,y)

-- * Template subroutines for dealing with quantum complex number encoded as pairs of 'QDouble'.

-- | Template version of 'mkPolar'.
template_mkPolar :: Circ (QDouble -> Circ (QDouble -> Circ (QDouble,QDouble)))
template_mkPolar = return \$ \x -> return \$ \y -> box "mkPolar" (uncurry q_mkPolar) (x,y)

-- | Template version of the constructor ':+' of 'Complex'.
template_symb_colon_symb_plus_ :: Circ (QDouble -> Circ (QDouble -> Circ (QDouble,QDouble)))
template_symb_colon_symb_plus_ = return \$ \x -> return \$ \y -> return (x,y)

-- | Template version of 'magnitude'.
template_magnitude :: Circ ((QDouble,QDouble) -> Circ QDouble)
template_magnitude = return \$ \p -> box "mag" q_magnitude p

-- | Template version of 'realPart'.
template_realPart :: Circ ((QDouble,QDouble) -> Circ QDouble)
template_realPart = return \$ \(x,y) -> return x

-- | Template version of 'imagPart'.
template_imagPart :: Circ ((QDouble,QDouble) -> Circ QDouble)
template_imagPart = return \$ \(x,y) -> return y

```