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