{-# LANGUAGE NoImplicitPrelude #-}
{-
Complex Gaussian bell multiplied with a polynomial.
In order to make this free of @pi@ factors,
we have to choose @recip (sqrt pi)@
as unit for translations and modulations,
for linear factors and in the differentiation.
-}
module MathObj.Gaussian.Polynomial where
import qualified MathObj.Gaussian.Bell as Bell
import qualified MathObj.Polynomial as Poly
import qualified Number.Complex as Complex
import qualified Algebra.Differential as Differential
import qualified Algebra.Transcendental as Trans
import qualified Algebra.Field as Field
-- import qualified Algebra.Real as Real
import qualified Algebra.Ring as Ring
import qualified Algebra.Additive as Additive
import Algebra.Transcendental (pi, )
import Algebra.Ring ((*), )
-- import Algebra.Additive ((+))
import NumericPrelude
import PreludeBase hiding (reverse, )
-- import Prelude ()
data T a = Cons {bell :: Bell.T a, polynomial :: Poly.T (Complex.T a)}
deriving (Eq, Show)
{-# INLINE evaluate #-}
evaluate :: (Trans.C a) =>
T a -> a -> Complex.T a
evaluate f x =
Bell.evaluateSqRt (bell f) x *
Poly.evaluate (polynomial f) (Complex.fromReal $ sqrt pi * x)
{- ToDo: evaluating a complex polynomial for a real argument can be optimized -}
multiply :: (Ring.C a) =>
T a -> T a -> T a
multiply x y =
Cons
(Bell.multiply (bell x) (bell y))
(polynomial x * polynomial y)
convolve :: (Field.C a) =>
T a -> T a -> T a
convolve f g =
reverse $ fourier $ multiply (fourier f) (fourier g)
reverse :: Additive.C a => T a -> T a
reverse x =
Cons
(Bell.reverse $ bell x)
(Poly.reverse $ polynomial x)
{-
We use a Horner like scheme
in order to translate multiplications with @id@
to differentations on the Fourier side.
Quadratic runtime.
fourier (Cons bell (Poly.const a + Poly.shift x))
= fourier (Cons bell (Poly.const a)) + fourier (Cons bell (Poly.shift x))
= fourier (Cons bell (Poly.const a)) + differentiate (fourier (Cons bell x))
untested
-}
fourier :: (Field.C a) =>
T a -> T a
fourier x =
foldr
(\c p ->
let q = differentiate p
in q{polynomial =
Poly.const c +
fmap Complex.quarterLeft (polynomial q)})
(Cons (Bell.fourier $ bell x) zero) $
Poly.coeffs $ polynomial x
{-
Differentiate and divide by @sqrt pi@ in order to stay in a ring.
This way, we do not need to fiddle with pi factors.
-}
differentiate :: (Ring.C a) => T a -> T a
differentiate f =
f{polynomial =
Bell.exponentPolynomial (bell f) * polynomial f +
Differential.differentiate (polynomial f)}
{- No Ring instance for Gaussians
instance (Ring.C a) => Differential.C (T a) where
differentiate = differentiate
-}
{- laws
differentiate (f*g) =
(differentiate f) * g + f * (differentiate g)
-}