{-# 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) -}