diff --git a/numeric-prelude.cabal b/numeric-prelude.cabal
--- a/numeric-prelude.cabal
+++ b/numeric-prelude.cabal
@@ -1,5 +1,5 @@
 Name:           numeric-prelude
-Version:        0.1
+Version:        0.1.1
 License:        GPL
 License-File:   LICENSE
 Author:         Dylan Thurston <dpt@math.harvard.edu>, Henning Thielemann <numericprelude@henning-thielemann.de>, Mikael Johansson
@@ -131,9 +131,12 @@
   default:     False
 
 Library
-  Build-Depends: parsec >=1 && <3, HUnit >=1 && <2, QuickCheck >=1 && <2
-  Build-Depends: non-negative >=0.0.2 && <0.1
-  Build-Depends: utility-ht >=0.0.4 && <0.1
+  Build-Depends:
+    parsec >=1 && <3,
+    QuickCheck >=1 && <2,
+    storable-record >=0.0.1 && <0.1,
+    non-negative >=0.0.2 && <0.1,
+    utility-ht >=0.0.4 && <0.1
   If flag(splitBase)
     Build-Depends:
       base >= 2 && <5,
@@ -222,10 +225,12 @@
     Number.Physical
     Number.Physical.Read
     Number.Physical.Show
+    NumericPrelude.Elementwise
     NumericPrelude
     PreludeBase
   Other-modules:
     NumericPrelude.List
+    Algebra.AffineSpace
     MathObj.Gaussian.Variance
     MathObj.Gaussian.Bell
     MathObj.Gaussian.Polynomial
@@ -247,5 +252,19 @@
     Test.MathObj.Gaussian.Variance
     Test.MathObj.Gaussian.Bell
   Main-Is: Test/Run.hs
-  If !flag(buildTests)
-    Buildable:         False
+  If flag(buildTests)
+    Build-Depends: HUnit >=1 && <2
+  Else
+    Buildable: False
+
+Executable test-gaussian
+  Hs-Source-Dirs: src, test
+  Main-Is: Gaussian.hs
+  Other-Modules:
+    MathObj.Gaussian.Example
+  If flag(buildTests)
+    Build-Depends:
+      gnuplot >=0.3 && <0.4,
+      HTam >=0.0.2 && <0.1
+  Else
+    Buildable: False
diff --git a/src/Algebra/Additive.hs b/src/Algebra/Additive.hs
--- a/src/Algebra/Additive.hs
+++ b/src/Algebra/Additive.hs
@@ -1,15 +1,19 @@
 {-# LANGUAGE NoImplicitPrelude #-}
 module Algebra.Additive (
-    {- * Class -}
+    -- * Class
     C,
     zero,
     (+), (-),
     negate, subtract,
 
-    {- * Complex functions -}
+    -- * Complex functions
     sum, sum1,
 
-    {- * Instances for atomic types -}
+    -- * Instance definition helpers
+    elementAdd, elementSub, elementNeg,
+    (<*>.+), (<*>.-), (<*>.-$),
+
+    -- * Instances for atomic types
     propAssociative,
     propCommutative,
     propIdentity,
@@ -21,6 +25,10 @@
 import Data.Int  (Int,  Int8,  Int16,  Int32,  Int64,  )
 import Data.Word (Word, Word8, Word16, Word32, Word64, )
 
+import qualified NumericPrelude.Elementwise as Elem
+import Control.Applicative (Applicative(pure, (<*>)), )
+import Data.Tuple.HT (fst3, snd3, thd3, )
+
 import qualified Data.Ratio as Ratio98
 import qualified Prelude as P
 import Prelude(Int, Integer, Float, Double, fromInteger, )
@@ -86,9 +94,63 @@
 
 
 
+{- |
+Instead of baking the add operation into the element function,
+we could use higher rank types
+and pass a generic @uncurry (+)@ to the run function.
+We do not do so in order to stay Haskell 98
+at least for parts of NumericPrelude.
+-}
+{-# INLINE elementAdd #-}
+elementAdd ::
+   (C x) =>
+   (v -> x) -> Elem.T (v,v) x
+elementAdd f =
+   Elem.element (\(x,y) -> f x + f y)
 
-{-* Instances for atomic types -}
+{-# INLINE elementSub #-}
+elementSub ::
+   (C x) =>
+   (v -> x) -> Elem.T (v,v) x
+elementSub f =
+   Elem.element (\(x,y) -> f x - f y)
 
+{-# INLINE elementNeg #-}
+elementNeg ::
+   (C x) =>
+   (v -> x) -> Elem.T v x
+elementNeg f =
+   Elem.element (negate . f)
+
+
+-- like <*>
+infixl 4 <*>.+, <*>.-, <*>.-$
+
+{- |
+> addPair :: (Additive.C a, Additive.C b) => (a,b) -> (a,b) -> (a,b)
+> addPair = Elem.run2 $ Elem.with (,) <*>.+  fst <*>.+  snd
+-}
+(<*>.+) ::
+   (C x) =>
+   Elem.T (v,v) (x -> a) -> (v -> x) -> Elem.T (v,v) a
+(<*>.+) f acc =
+   f <*> elementAdd acc
+
+(<*>.-) ::
+   (C x) =>
+   Elem.T (v,v) (x -> a) -> (v -> x) -> Elem.T (v,v) a
+(<*>.-) f acc =
+   f <*> elementSub acc
+
+(<*>.-$) ::
+   (C x) =>
+   Elem.T v (x -> a) -> (v -> x) -> Elem.T v a
+(<*>.-$) f acc =
+   f <*> elementNeg acc
+
+
+-- * Instances for atomic types
+
 instance C Integer where
    {-# INLINE zero #-}
    {-# INLINE negate #-}
@@ -224,27 +286,27 @@
 
 
 
-{-* Instances for composed types -}
+-- * Instances for composed types
 
 instance (C v0, C v1) => C (v0, v1) where
    {-# INLINE zero #-}
    {-# INLINE negate #-}
    {-# INLINE (+) #-}
    {-# INLINE (-) #-}
-   zero                   = (zero, zero)
-   (+)    (x0,x1) (y0,y1) = ((+) x0 y0, (+) x1 y1)
-   (-)    (x0,x1) (y0,y1) = ((-) x0 y0, (-) x1 y1)
-   negate (x0,x1)         = (negate x0, negate x1)
+   zero   = (,) zero zero
+   (+)    = Elem.run2 $ pure (,) <*>.+  fst <*>.+  snd
+   (-)    = Elem.run2 $ pure (,) <*>.-  fst <*>.-  snd
+   negate = Elem.run  $ pure (,) <*>.-$ fst <*>.-$ snd
 
 instance (C v0, C v1, C v2) => C (v0, v1, v2) where
    {-# INLINE zero #-}
    {-# INLINE negate #-}
    {-# INLINE (+) #-}
    {-# INLINE (-) #-}
-   zero                         = (zero, zero, zero)
-   (+)    (x0,x1,x2) (y0,y1,y2) = ((+) x0 y0, (+) x1 y1, (+) x2 y2)
-   (-)    (x0,x1,x2) (y0,y1,y2) = ((-) x0 y0, (-) x1 y1, (-) x2 y2)
-   negate (x0,x1,x2)            = (negate x0, negate x1, negate x2)
+   zero   = (,,) zero zero zero
+   (+)    = Elem.run2 $ pure (,,) <*>.+  fst3 <*>.+  snd3 <*>.+  thd3
+   (-)    = Elem.run2 $ pure (,,) <*>.-  fst3 <*>.-  snd3 <*>.-  thd3
+   negate = Elem.run  $ pure (,,) <*>.-$ fst3 <*>.-$ snd3 <*>.-$ thd3
 
 
 instance (C v) => C [v] where
@@ -268,7 +330,7 @@
    (-)    f g x = (-) (f x) (g x)
    negate f   x = negate (f x)
 
-{- * Properties -}
+-- * Properties
 
 propAssociative :: (Eq a, C a) => a -> a -> a -> Bool
 propCommutative :: (Eq a, C a) => a -> a -> Bool
diff --git a/src/Algebra/AffineSpace.hs b/src/Algebra/AffineSpace.hs
new file mode 100644
--- /dev/null
+++ b/src/Algebra/AffineSpace.hs
@@ -0,0 +1,247 @@
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE FlexibleInstances #-}
+{- |
+This module is not yet exported
+since its interface is not mature.
+There are two approaches for representing affine spaces:
+
+[1] Two sets: A set of points and a set of vectors.
+    Examples: Absolute potential and voltage,
+    absolute temperature and temperature difference.
+    Operations are
+      add :: vector -> point -> point
+      sub :: point -> point -> vector
+
+[2] One set for points, no vectors.
+    Examples: Interpolation
+    Operation:
+      combine :: [(coefficient, vector)] -> vector
+    Where it must be asserted,
+    that the coefficients sum up to 1.
+
+The second one is the one we follow here.
+It is more similar to Module and VectorSpace.
+-}
+module Algebra.AffineSpace where
+
+import qualified Algebra.PrincipalIdealDomain as PID
+import qualified Algebra.Additive as Additive
+import qualified Algebra.Module as Module
+import qualified Number.Ratio as Ratio
+
+import qualified Number.Complex as Complex
+
+import Control.Applicative (Applicative(pure, (<*>)), )
+
+import NumericPrelude hiding (zero, )
+import PreludeBase
+import Prelude ()
+
+{- |
+The type class is for representing affine spaces via affine combinations.
+However, we didn't find a way to both ensure the property
+that the combination coefficients sum up to 1,
+and keep it efficient.
+
+We propose this class instead of a combination of Additive and Module
+for interpolation for those types,
+where scaling and addition alone makes no sense.
+Such types are e.g. internal filter parameters in signal processing:
+For these types interpolation makes definitely sense,
+but addition and scaling make not.
+
+That is, both classes are isomorphic
+(you can define one in terms of the other),
+but instances of this class are more easily defined,
+and using an AffineSpace constraint instead of Module in a type signature
+is important for documentation purposes.
+AffineSpace should be superclass of Module.
+(But then you may ask, why not adding another superclass for Convex spaces.
+This class would provide a linear combination operation,
+where the coefficients sum up to one
+and all of them are non-negative.)
+
+We may add a safety layer that ensures
+that the coefficients sum up to 1,
+using start points on the simplex
+and functions to move on the simplex.
+Start points have components that sum up to 1, e.g.
+
+> (1, 0, 0, 0)
+> (0, 1, 0, 0)
+> (0, 0, 1, 0)
+> (0, 0, 0, 1)
+> (1/4, 1/4, 1/4, 1/4)
+
+Then you may move along the simplex in the directions
+
+> (1,  -1, 0,  0)
+> (0,   1, 0, -1)
+> (-1, -1, 3, -1)
+
+which are characterized by components that sum up to 0.
+
+For example linear combination is defined by
+
+> lerp k (a,b) = (1-k)*>a + k*>b
+
+that is the coefficients are (1-k) and k.
+The pair (1-k, k) can be constructed
+by starting at pair (1,0)
+and moving k units in direction (-1,1).
+
+> (1-k, k) = (1,0) + k*(-1,1)
+
+It is however a challenge to manage the coefficient tuples
+in a type safe and efficient way.
+For small numbers of interpolation nodes
+(constant, linear, cubic interpolation)
+a type level list would appropriate,
+but what to do for large tuples
+like for Whittaker interpolation?
+
+
+As side note:
+In principle it would be sufficient
+to provide an affine combination of two points,
+since all affine combinations of more points
+can be decomposed into such simple combinations.
+
+> lerp a x y = (1-a)*>x + a*>y
+
+E.g. @a*>x + b*>y + c*>z@ with @a+b+c=1@
+can be written as @lerp c (lerp (b/(1-c)) x y) z@.
+More generally you can use
+
+> lerpnorm a b x y
+>    = lerp (b/(a+b)) x y
+>    = (a/(a+b))*>x + (b/(a+b))*>y
+
+for writing
+
+> a*>x + b*>y + c*>z ==
+>    lerpnorm (a+b) c (lerpnorm a b x y) z
+
+or
+
+> a*>x + b*>y + c*>z + d*>w ==
+>    lerpnorm (a+b+c) d (lerpnorm (a+b) c (lerpnorm a b x y) z) w
+
+with @a+b+c+d=1@.
+
+The downside is, that lerpnorm requires division, that is, a field,
+whereas the computation of the coefficients
+sometimes only requires ring operations.
+-}
+class Zero v => C a v where
+   multiplyAccumulate :: (a,v) -> v -> v
+
+class Zero v where
+   zero :: v
+
+
+instance Zero Float where
+   {-# INLINE zero #-}
+   zero = Additive.zero
+
+instance Zero Double where
+   {-# INLINE zero #-}
+   zero = Additive.zero
+
+instance (Zero a) => Zero (Complex.T a) where
+   {-# INLINE zero #-}
+   zero = zero Complex.+: zero
+
+instance (PID.C a) => Zero (Ratio.T a) where
+   {-# INLINE zero #-}
+   zero = Additive.zero
+
+
+instance C Float Float where
+   {-# INLINE multiplyAccumulate #-}
+   multiplyAccumulate (a,x) y = a*x+y
+
+instance C Double Double where
+   {-# INLINE multiplyAccumulate #-}
+   multiplyAccumulate (a,x) y = a*x+y
+
+instance (C a v) => C a (Complex.T v) where
+   {-# INLINE multiplyAccumulate #-}
+   multiplyAccumulate =
+      makeMac2 (Complex.+:) Complex.real Complex.imag
+
+instance (PID.C a) => C (Ratio.T a) (Ratio.T a) where
+   {-# INLINE multiplyAccumulate #-}
+   multiplyAccumulate (a,x) y = a*x+y
+
+
+infixl 6 *.+
+
+{- |
+Infix variant of 'multiplyAccumulate'.
+-}
+{-# INLINE (*.+) #-}
+(*.+) :: C a v => v -> (a,v) -> v
+(*.+) = flip multiplyAccumulate
+
+
+-- * convenience functions for defining multiplyAccumulate
+
+{-# INLINE multiplyAccumulateModule #-}
+multiplyAccumulateModule ::
+   Module.C a v =>
+   (a,v) -> v -> v
+multiplyAccumulateModule (a,x) y =
+   a *> x + y
+
+
+{- |
+A special reader monad.
+-}
+newtype MAC a v x = MAC {runMac :: (a,v) -> v -> x}
+
+{-# INLINE element #-}
+element ::
+   (C a x) =>
+   (v -> x) -> MAC a v x
+element f =
+   MAC (\(a,x) y -> multiplyAccumulate (a, f x) (f y))
+
+instance Functor (MAC a v) where
+   {-# INLINE fmap #-}
+   fmap f (MAC x) =
+      MAC $ \av v -> f $ x av v
+
+instance Applicative (MAC a v) where
+   {-# INLINE pure #-}
+   {-# INLINE (<*>) #-}
+   pure x = MAC $ \ _av _v -> x
+   MAC f <*> MAC x =
+      MAC $ \av v -> f av v $ x av v
+
+{-# INLINE makeMac #-}
+makeMac ::
+   (C a x) =>
+   (x -> v) ->
+   (v -> x) ->
+   (a,v) -> v -> v
+makeMac cons x =
+   runMac $ pure cons <*> element x
+
+{-# INLINE makeMac2 #-}
+makeMac2 ::
+   (C a x, C a y) =>
+   (x -> y -> v) ->
+   (v -> x) -> (v -> y) ->
+   (a,v) -> v -> v
+makeMac2 cons x y =
+   runMac $ pure cons <*> element x <*> element y
+
+{-# INLINE makeMac3 #-}
+makeMac3 ::
+   (C a x, C a y, C a z) =>
+   (x -> y -> z -> v) ->
+   (v -> x) -> (v -> y) -> (v -> z) ->
+   (a,v) -> v -> v
+makeMac3 cons x y z =
+   runMac $ pure cons <*> element x <*> element y <*> element z
diff --git a/src/Algebra/Module.hs b/src/Algebra/Module.hs
--- a/src/Algebra/Module.hs
+++ b/src/Algebra/Module.hs
@@ -23,23 +23,31 @@
 import qualified Algebra.Laws as Laws
 
 import Algebra.Ring     ((*), fromInteger, )
-import Algebra.Additive ((+), zero, )
+import Algebra.Additive ((+), zero, sum, )
 
+import qualified NumericPrelude.Elementwise as Elem
+import Control.Applicative (Applicative(pure, (<*>)), )
+
 import Data.Function.HT (powerAssociative, )
-import Data.List (map, zipWith, foldl, )
+import Data.List (map, zipWith, )
+import Data.Tuple.HT (fst3, snd3, thd3, )
+import Data.Tuple (fst, snd, )
 
-import Prelude((.), Eq, Bool, Int, Integer, Float, Double)
+import Prelude((.), Eq, Bool, Int, Integer, Float, Double, ($), )
 -- import qualified Prelude as P
 
 
 -- Is this right?
 infixr 7 *>
 
-{- Functional dependency can't be used
-   since the instance (Algebra.Module.C a a)
-   would conflict with all others.
-   class Algebra.Module.C b a | b -> a where -}
+{-
+Functional dependency can't be used
+since @Complex.T a@ is a module
+with respect to both @a@ and @Complex.T a@.
 
+class Algebra.Module.C a v | v -> a where
+-}
+
 {-|
 A Module over a ring satisfies:
 
@@ -47,10 +55,19 @@
 >   (a * b) *> c === a *> (b *> c)
 >   (a + b) *> c === a *> c + b *> c
 -}
-class (Additive.C b, Ring.C a) => C a b where
+class (Ring.C a, Additive.C v) => C a v where
     -- | scale a vector by a scalar
-    (*>) :: a -> b -> b
+    (*>) :: a -> v -> v
 
+
+(<*>.*>) ::
+   (C a x) =>
+   Elem.T (a,v) (x -> c) -> (v -> x) -> Elem.T (a,v) c
+(<*>.*>) f acc =
+   f <*> Elem.element (\(a,v) -> a *> acc v)
+
+
+
 {-* Instances for atomic types -}
 
 instance C Float Float where
@@ -83,17 +100,19 @@
 
 instance (C a b0, C a b1) => C a (b0, b1) where
    {-# INLINE (*>) #-}
-   s *> (x0,x1)   = (s *> x0, s *> x1)
+   (*>) = Elem.run2 $ pure (,) <*>.*> fst <*>.*> snd
+   -- s *> (x0,x1)   = (s *> x0, s *> x1)
 
 instance (C a b0, C a b1, C a b2) => C a (b0, b1, b2) where
    {-# INLINE (*>) #-}
-   s *> (x0,x1,x2) = (s *> x0, s *> x1, s *> x2)
+   (*>) = Elem.run2 $ pure (,,) <*>.*> fst3 <*>.*> snd3 <*>.*> thd3
+   -- s *> (x0,x1,x2) = (s *> x0, s *> x1, s *> x2)
 
-instance (C a b) => C a [b] where
+instance (C a v) => C a [v] where
    {-# INLINE (*>) #-}
    (*>) = map . (*>)
 
-instance (C a b) => C a (c -> b) where
+instance (C a v) => C a (c -> v) where
    {-# INLINE (*>) #-}
    (*>) s f = (*>) s . f
 
@@ -106,8 +125,8 @@
 ToDo:
 Should it use 'NumericPrelude.List.zipWithMatch' ?
 -}
-linearComb :: C a b => [a] -> [b] -> b
-linearComb c = foldl (+) zero . zipWith (*>) c
+linearComb :: C a v => [a] -> [v] -> v
+linearComb c = sum . zipWith (*>) c
 
 {-|
 This function can be used to define any
@@ -116,18 +135,18 @@
 Better move to "Algebra.Additive"?
 -}
 {-# INLINE integerMultiply #-}
-integerMultiply :: (ToInteger.C a, Additive.C b) => a -> b -> b
-integerMultiply a b =
-   powerAssociative (+) zero b (ToInteger.toInteger a)
+integerMultiply :: (ToInteger.C a, Additive.C v) => a -> v -> v
+integerMultiply a v =
+   powerAssociative (+) zero v (ToInteger.toInteger a)
 
 
 {- * Properties -}
 
-propCascade :: (Eq b, C a b) => b -> a -> a -> Bool
+propCascade :: (Eq v, C a v) => v -> a -> a -> Bool
 propCascade  =  Laws.leftCascade (*) (*>)
 
-propRightDistributive :: (Eq b, C a b) => a -> b -> b -> Bool
+propRightDistributive :: (Eq v, C a v) => a -> v -> v -> Bool
 propRightDistributive  =  Laws.rightDistributive (*>) (+)
 
-propLeftDistributive :: (Eq b, C a b) => b -> a -> a -> Bool
+propLeftDistributive :: (Eq v, C a v) => v -> a -> a -> Bool
 propLeftDistributive x  =  Laws.homomorphism (*>x) (+) (+)
diff --git a/src/MathObj/Gaussian/Bell.hs b/src/MathObj/Gaussian/Bell.hs
--- a/src/MathObj/Gaussian/Bell.hs
+++ b/src/MathObj/Gaussian/Bell.hs
@@ -20,40 +20,50 @@
 import Algebra.Additive ((+), )
 
 import Test.QuickCheck (Arbitrary, arbitrary, coarbitrary, )
-import Control.Monad (liftM3, )
+import Control.Monad (liftM4, )
 
 -- import Prelude (($))
 import NumericPrelude
 import PreludeBase hiding (reverse, )
 
 
-data T a = Cons {c0, c1 :: Complex.T a, c2 :: a}
+data T a = Cons {amp :: a, c0, c1 :: Complex.T a, c2 :: a}
    deriving (Eq, Show)
 
 instance (Real.C a, Arbitrary a) => Arbitrary (T a) where
    arbitrary =
-      liftM3
-         (\a b c -> Cons a b (1 + abs c))
-         arbitrary arbitrary arbitrary
+      liftM4
+         (\k a b c -> Cons k a b (1 + abs c))
+         arbitrary arbitrary arbitrary arbitrary
    coarbitrary = undefined
 
 
-constant :: Additive.C a => T a
-constant = Cons zero zero zero
+constant :: Ring.C a => T a
+constant = Cons one zero zero zero
 
+{- |
+eigenfunction of 'fourier'
+-}
+unit :: Ring.C a => T a
+unit = Cons one zero zero one
+
 {-# INLINE evaluate #-}
 evaluate :: (Trans.C a) =>
    T a -> a -> Complex.T a
 evaluate f x =
-   Complex.exp $ Complex.scale (-pi) $
-   c0 f + Complex.scale x (c1 f) + Complex.fromReal (c2 f * x^2)
+   Complex.scale
+     (sqrt (amp f))
+     (Complex.exp $ Complex.scale (-pi) $
+      c0 f + Complex.scale x (c1 f) + Complex.fromReal (c2 f * x^2))
 
 evaluateSqRt :: (Trans.C a) =>
    T a -> a -> Complex.T a
 evaluateSqRt f x0 =
-   let x = sqrt pi * x0
-   in  Complex.exp $ negate $
-       c0 f + Complex.scale x (c1 f) + Complex.fromReal (c2 f * x^2)
+   Complex.scale
+     (sqrt (amp f))
+     (let x = sqrt pi * x0
+      in  Complex.exp $ negate $
+          c0 f + Complex.scale x (c1 f) + Complex.fromReal (c2 f * x^2))
 
 exponentPolynomial :: (Additive.C a) =>
    T a -> Poly.T (Complex.T a)
@@ -61,10 +71,17 @@
    Poly.fromCoeffs [c0 f, c1 f, Complex.fromReal (c2 f)]
 
 
-multiply :: (Additive.C a) =>
+variance :: (Trans.C a) =>
+   T a -> a
+variance f =
+   recip $ c2 f * 2*pi
+
+multiply :: (Ring.C a) =>
    T a -> T a -> T a
 multiply f g =
-   Cons (c0 f + c0 g) (c1 f + c1 g) (c2 f + c2 g)
+   Cons
+      (amp f * amp g)
+      (c0 f + c0 g) (c1 f + c1 g) (c2 f + c2 g)
 
 
 {-
@@ -89,6 +106,7 @@
          = -(f1 - g1)^2/(4*(f2 + g2))
        -}
    in  Cons
+          ((amp f * amp g) / (c2 f + c2 g))
           (c0 f + c0 g
               - Complex.scale (recip (4*s)) ((c1 f - c1 g)^2))
           (Complex.scale (c2 g / s) (c1 f) +
@@ -110,6 +128,7 @@
        g1 = translateComplex gd g0
    in  translateComplex (negate $ fd + gd) $
        Cons
+          ((amp f0 * amp g0) / (c2 f0 + c2 g0))
           (c0 f1 + c0 g1) zero
           (recip $ recip (c2 f1) + recip (c2 g1))
 
@@ -126,6 +145,7 @@
        c = c2 f
        rc = recip c
    in  Cons
+          (amp f / c2 f)
           (Complex.scale (rc/4) (-b^2) + a)
           (Complex.scale rc $ Complex.quarterRight b)
           rc
@@ -134,7 +154,7 @@
    T a -> T a
 fourierByTranslation f =
    translateComplex (Complex.scale (1/2) $ Complex.quarterLeft $ c1 f) $
-   Cons (c0 f) zero (recip $ c2 f)
+   Cons (amp f / c2 f) (c0 f) zero (recip $ c2 f)
 
 {-
 a + b*x + c*x^2
@@ -201,6 +221,7 @@
        b = c1 f
        c = c2 f
    in  Cons
+          (amp f)
           (Complex.fromReal (c*d^2) - Complex.scale d b + a)
           (Complex.fromReal (-2*c*d) + b)
           c
@@ -211,6 +232,7 @@
        b = c1 f
        c = c2 f
    in  Cons
+          (amp f)
           (Complex.scale c (d^2) - b*d + a)
           (Complex.scale (-2*c) d + b)
           c
@@ -218,6 +240,7 @@
 modulate :: Ring.C a => a -> T a -> T a
 modulate d f =
    Cons
+      (amp f)
       (c0 f)
       (c1 f + (zero +: 2*d))
       (c2 f)
@@ -225,6 +248,7 @@
 turn :: Ring.C a => a -> T a -> T a
 turn d f =
    Cons
+      (amp f)
       (c0 f + (zero +: 2*d))
       (c1 f)
       (c2 f)
@@ -237,6 +261,7 @@
 dilate :: Field.C a => a -> T a -> T a
 dilate k f =
    Cons
+      (amp f)
       (c0 f)
       (Complex.scale (recip k) $ c1 f)
       (c2 f / k^2)
@@ -244,9 +269,18 @@
 shrink :: Ring.C a => a -> T a -> T a
 shrink k f =
    Cons
+      (amp f)
       (c0 f)
       (Complex.scale k $ c1 f)
       (k^2 * c2 f)
+
+amplify :: (Ring.C a) => a -> T a -> T a
+amplify k f =
+   Cons
+      (k^2 * amp f)
+      (c0 f)
+      (c1 f)
+      (c2 f)
 
 
 {- laws
diff --git a/src/MathObj/Gaussian/Example.hs b/src/MathObj/Gaussian/Example.hs
new file mode 100644
--- /dev/null
+++ b/src/MathObj/Gaussian/Example.hs
@@ -0,0 +1,227 @@
+{-# LANGUAGE NoImplicitPrelude #-}
+{-
+Reciprocal of variance of a Gaussian bell curve.
+We describe the curve only in terms of its variance
+thus we represent a bell curve at the coordinate origin
+neglecting its amplitude.
+
+We could also define the amplitude as @root 4 c@,
+thus preserving L2 norm being one,
+but then @dilate@ and @shrink@ also include an amplification.
+
+We could do some projective geometry in the exponent
+in order to also have zero variance,
+which corresponds to the dirac impulse.
+-}
+module MathObj.Gaussian.Example where
+
+import qualified MathObj.Gaussian.Polynomial as PolyBell
+import qualified MathObj.Gaussian.Bell as Bell
+import qualified MathObj.Gaussian.Variance as Var
+
+import qualified MathObj.Polynomial as Poly
+
+import qualified Algebra.Transcendental as Trans
+import qualified Algebra.Algebraic      as Algebraic
+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 qualified Number.Complex as Complex
+
+import Algebra.Transcendental (pi, )
+import Algebra.Algebraic (root, )
+import Algebra.Ring ((*), (^), )
+
+import Number.Complex ((+:), )
+
+import qualified Numerics.Function as Func
+import qualified Numerics.Fourier as Fourier
+import qualified Numerics.Integration as Integ
+import qualified Numerics.Differentiation as Diff
+
+import qualified Graphics.Gnuplot.Simple as GP
+
+import Control.Applicative (liftA2, )
+
+-- import System.Exit (ExitCode, )
+
+-- import Prelude (($))
+import NumericPrelude
+import PreludeBase
+import qualified Prelude as P
+
+
+curve0 :: Var.T Double
+curve0 = curve0a
+
+curve0a :: Var.T Double
+curve0a = Var.Cons 1.4 3.3
+
+curve0b :: Var.T Double
+curve0b = Var.Cons 2.2 1.7
+
+variance0 :: (Double, Double)
+variance0 =
+   (Var.variance curve0,
+    (Integ.rectangular 1000 (-2,2) $ liftA2 (*) (^2) (Var.evaluate curve0)) /
+    (Integ.rectangular 1000 (-2,2) $ Var.evaluate curve0))
+
+norm10 :: (Double, Double)
+norm10 =
+   (Integ.rectangular 1000 (-2,2) $ Var.evaluate curve0,
+    Var.norm1 curve0)
+
+norm20 :: (Double, Double)
+norm20 =
+   (sqrt $ Integ.rectangular 1000 (-2,2) $ (^2) . Var.evaluate curve0,
+    Var.norm2 curve0)
+
+norm30 :: (Double, Double)
+norm30 =
+   (root 3 $ Integ.rectangular 1000 (-2,2) $ (^3) . Var.evaluate curve0,
+    Var.normP 3 curve0)
+
+fourier0 :: IO ()
+fourier0 =
+   GP.plotFuncs []
+      (GP.linearScale 100 (-2,2))
+      [Var.evaluate $ Var.fourier curve0,
+       Fourier.analysisTransformOneReal 100 (-2,2) $ Var.evaluate curve0]
+
+multiply0 :: IO ()
+multiply0 =
+   GP.plotFuncs []
+      (GP.linearScale 100 (-1,1))
+      [Var.evaluate $ Var.multiply curve0a curve0b,
+       liftA2 (*) (Var.evaluate curve0a) (Var.evaluate curve0b)]
+
+convolve0 :: IO ()
+convolve0 =
+   GP.plotFuncs []
+      (GP.linearScale 100 (-2,2))
+      [Var.evaluate $ Var.convolve curve0a curve0b,
+       Integ.convolve 1000 (-3,3) (Var.evaluate curve0a) (Var.evaluate curve0b)]
+
+
+curve1 :: Bell.T Double
+curve1 = curve1a
+
+curve1a :: Bell.T Double
+curve1a = Bell.Cons 1.4 (0.1+:0.3) ((-0.2)+:1.4) 2.3
+
+curve1b :: Bell.T Double
+curve1b = Bell.Cons 2.2 ((-0.3)+:2.1) (0.2+:(-0.4)) 1.7
+
+variance1 :: (Double, Double)
+variance1 =
+   (Bell.variance curve1,
+    (Integ.rectangular 1000 (-2,2) $
+        liftA2 (*) (^2)
+           (Complex.magnitudeSqr .
+            Func.translateRight
+               (Complex.real (Bell.c1 curve1) / (2 * Bell.c2 curve1))
+               (Bell.evaluate curve1))) /
+    (Integ.rectangular 1000 (-2,2) $ Complex.magnitude . Bell.evaluate curve1))
+
+{- the norm depends on too much things
+norm0vs1 :: (Double, Double)
+norm0vs1 =
+   ((Integ.rectangular 1000 (-5,5) $ Var.evaluate curve0)
+         * exp (- Complex.real (Bell.c0 curve1)),
+    Integ.rectangular 1000 (-5,5) $ Complex.magnitude . Bell.evaluate curve1)
+-}
+
+fourier1 :: IO ()
+fourier1 =
+   GP.plotFuncs []
+      (GP.linearScale 100 (-5,5))
+      [Complex.real . (Bell.evaluate $ Bell.fourier curve1),
+       fourierAnalysisReal 100 (-2,2) $ Bell.evaluate curve1]
+
+
+curve2 :: PolyBell.T Double
+curve2 =
+   PolyBell.Cons
+--      Bell.unit
+--      (Bell.Cons 1.4 (0.1+:0.3) 0 1.2)
+--      (Bell.Cons 1.4 (0.1+:0.3) ((-0.2)+:1.4) 1)
+      curve1
+--      (Poly.fromCoeffs [one])
+--      (Poly.fromCoeffs [zero,one])
+--      (Poly.fromCoeffs [zero,zero,one])
+--      (Poly.fromCoeffs [0,Complex.imaginaryUnit])
+      (Poly.fromCoeffs [1.4+:(-0.1),0.8+:(0.1),(-1.1)+:0.3])
+
+differentiate2 :: IO ()
+differentiate2 =
+   GP.plotFuncs []
+      (GP.linearScale 100 (-2,2))
+      [Complex.real . (PolyBell.evaluateSqRt $ PolyBell.differentiate curve2),
+       ((/ sqrt pi) . ) $ Diff.diff (1e-5) $ Complex.real . PolyBell.evaluateSqRt curve2]
+
+fourier2 :: IO ()
+fourier2 =
+   GP.plotFuncs []
+      (GP.linearScale 100 (-5,5))
+      [Complex.real . (PolyBell.evaluateSqRt $ PolyBell.fourier curve2),
+       fourierAnalysisReal 100 (-2,2) $ PolyBell.evaluateSqRt curve2]
+
+
+
+fourierAnalysisReal ::
+   (P.Floating a) =>
+   Integer -> (a, a) -> (a -> Complex.T a) -> a -> a
+fourierAnalysisReal n rng f =
+   liftA2 (P.-)
+      (Fourier.analysisTransformOneReal n rng (Complex.real . f))
+      (Fourier.analysisTransformOneImag n rng (Complex.imag . f))
+
+
+{- |
+Try to approximate @\x -> exp (-x^2) * x@
+by a difference of translated Gaussian bells.
+
+exp(-x^2) * x
+  ==  exp(-(a+b*x+c*x^2)) - exp(-(a-b*x+c*x^2))
+  ==  exp(-(a+c*x^2)) * (exp(-b*x) - exp(b*x))
+  ==  exp(-(a+c*x^2)) * 2*sinh (b*x)
+
+It holds
+  lim (\b x -> sinh (b*x) / b)  =  id
+-}
+diffApprox :: IO ()
+diffApprox =
+   let amp = (2*b)^- (-2)
+       a = 0
+       {-
+       amp = 1
+       a = log (2 * abs b)
+       -}
+       b = -0.1
+       c = 1
+       ac = Complex.fromReal a
+       bc = Complex.fromReal b
+   in  GP.plotFuncs []
+          (GP.linearScale 100 (-2,2::Double))
+          [Complex.real .
+           (PolyBell.evaluateSqRt $
+              PolyBell.Cons Bell.unit (Poly.fromCoeffs [zero,one])),
+           Complex.real .
+           liftA2 (-)
+             (PolyBell.evaluateSqRt $
+                PolyBell.Cons (Bell.Cons amp ac bc c) (Poly.fromCoeffs [one]))
+             (PolyBell.evaluateSqRt $
+                PolyBell.Cons (Bell.Cons amp ac (-bc) c) (Poly.fromCoeffs [one]))]
+
+
+polyApprox :: IO ()
+polyApprox =
+   GP.plotFuncs []
+      (GP.linearScale 100 (-2,2::Double))
+      [Complex.real .
+         PolyBell.evaluateSqRt curve2,
+       Complex.real . sum .
+         mapM (\(amp,b) -> \x -> amp * Bell.evaluateSqRt b x)
+         (PolyBell.approximateByBells 0.1 curve2)]
diff --git a/src/MathObj/Gaussian/Polynomial.hs b/src/MathObj/Gaussian/Polynomial.hs
--- a/src/MathObj/Gaussian/Polynomial.hs
+++ b/src/MathObj/Gaussian/Polynomial.hs
@@ -7,17 +7,30 @@
 as unit for translations and modulations,
 for linear factors and in the differentiation.
 -}
+{-
+ToDo:
+
+* In order to avoid the weird @sqrt pi@ factor,
+  use a polynomial expression in @pi@.
+
+* sum of multiple bells using Data.Map from exponent polynomial to coefficient polynomial
+  use of Algebra object.
+
+* Projective geometry in order to support Dirac impulse.
+-}
 module MathObj.Gaussian.Polynomial where
 
 import qualified MathObj.Gaussian.Bell as Bell
 
+import qualified MathObj.LaurentPolynomial as LPoly
 import qualified MathObj.Polynomial as Poly
 import qualified Number.Complex     as Complex
 
+import qualified Algebra.ZeroTestable   as ZeroTestable
 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.Real           as Real
 import qualified Algebra.Ring           as Ring
 import qualified Algebra.Additive       as Additive
 
@@ -25,65 +38,165 @@
 import Algebra.Ring ((*), )
 -- import Algebra.Additive ((+))
 
+import qualified Data.Record.HT as Rec
+import qualified Data.List as List
+import Data.Function.HT (nest, )
+import Data.Eq.HT (equating, )
+import Data.List.HT (mapAdjacent, )
+import Data.Tuple.HT (forcePair, )
+
+import Test.QuickCheck (Arbitrary, arbitrary, coarbitrary, )
+import Control.Monad (liftM2, )
+
 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)
+   deriving (Show)
 
+instance Real.C a => Eq (T a) where
+   (==) = equal
 
-{-# INLINE evaluate #-}
-evaluate :: (Trans.C a) =>
+{-
+The derived Eq is not correct.
+We have to combine the amplitude of the bell with the polynomial,
+respecting signs and the square root of the bell amplitude.
+-}
+equal :: Real.C a => T a -> T a -> Bool
+equal x y =
+   let bx = bell x
+       by = bell y
+       csign c =
+          Complex.real c > 0 ||
+          (Complex.real c == 0 && Complex.imag c > 0)
+       scaleSqr b =
+          map (\c -> (Complex.scale (Bell.amp b) (c^2), csign c)) .
+          Poly.coeffs . polynomial
+   in  Rec.equal
+          (equating Bell.c0 :
+           equating Bell.c1 :
+           equating Bell.c2 :
+           [])
+          bx by
+       &&
+       scaleSqr by x == scaleSqr bx y
+
+
+instance (Real.C a, Arbitrary a) => Arbitrary (T a) where
+   arbitrary =
+--      liftM2 Cons arbitrary arbitrary
+      liftM2 Cons
+         arbitrary
+         -- we have to restrict the number of polynomial coefficients,
+         -- since with the quadratic time algorithms like fourier and convolve,
+         -- in connection with Rational slow down tests too much.
+         (fmap (Poly.fromCoeffs . take 5 . Poly.coeffs) arbitrary)
+   coarbitrary = undefined
+
+
+
+{-# INLINE evaluateSqRt #-}
+evaluateSqRt :: (Trans.C a) =>
    T a -> a -> Complex.T a
-evaluate f x =
+evaluateSqRt 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 -}
 
 
+constant :: (Ring.C a) => T a
+constant =
+   Cons Bell.constant (Poly.const one)
+
+scale :: (Ring.C a) => a -> T a -> T a
+scale x f =
+   f{polynomial = fmap (Complex.scale x) $ polynomial f}
+
+scaleComplex :: (Ring.C a) => Complex.T a -> T a -> T a
+scaleComplex x f =
+   f{polynomial = fmap (x*) $ polynomial f}
+
+
+eigenfunction :: (Field.C a) => Int -> T a
+eigenfunction =
+   eigenfunctionDifferential
+
+eigenfunction0 :: (Ring.C a) => T a
+eigenfunction0 =
+   Cons Bell.unit (Poly.fromCoeffs [one])
+
+eigenfunction1 :: (Ring.C a) => T a
+eigenfunction1 =
+   Cons Bell.unit (Poly.fromCoeffs [zero, one])
+
+eigenfunction2 :: (Field.C a) => T a
+eigenfunction2 =
+   Cons Bell.unit (Poly.fromCoeffs [-(1/4), zero, one])
+
+eigenfunction3 :: (Field.C a) => T a
+eigenfunction3 =
+   Cons Bell.unit (Poly.fromCoeffs [zero, -(3/4), zero, one])
+
+
+eigenfunctionDifferential :: (Field.C a) => Int -> T a
+eigenfunctionDifferential n =
+   (\f -> f{bell = Bell.unit}) $
+   nest n (scale (-1/4) . differentiate) $
+   Cons (Bell.Cons one zero zero 2) one
+
+eigenfunctionIterative :: (Field.C a, Real.C a) => Int -> T a
+eigenfunctionIterative n =
+   fst . head . dropWhile (uncurry (/=)) . mapAdjacent (,) $
+   eigenfunctionIteration $
+   Cons
+      Bell.unit
+      (Poly.fromCoeffs $ replicate n zero ++ [one])
+
+eigenfunctionIteration :: (Field.C a) => T a -> [T a]
+eigenfunctionIteration =
+   iterate (\x ->
+      let y = fourier x
+          px = polynomial x
+          py = polynomial y
+          c = last (Poly.coeffs px) / last (Poly.coeffs py)
+      in  y{polynomial = fmap (0.5*) (px + fmap (c*) py)})
+
+
 multiply :: (Ring.C a) =>
    T a -> T a -> T a
-multiply x y =
+multiply f g =
    Cons
-      (Bell.multiply (bell x) (bell y))
-      (polynomial x * polynomial y)
+      (Bell.multiply (bell f) (bell g))
+      (polynomial f * polynomial g)
 
 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 (Cons bell (Poly.const a + Poly.shift f))
+  = fourier (Cons bell (Poly.const a)) + fourier (Cons bell (Poly.shift f))
+  = fourier (Cons bell (Poly.const a)) + C * differentiate (fourier (Cons bell f))
 -}
 fourier :: (Field.C a) =>
    T a -> T a
-fourier x =
+fourier f =
    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
+                   fmap (Complex.scale (1/2) . Complex.quarterLeft) (polynomial q)})
+      (Cons (Bell.fourier $ bell f) zero) $
+   Poly.coeffs $ polynomial f
 
 {-
 Differentiate and divide by @sqrt pi@ in order to stay in a ring.
@@ -92,8 +205,174 @@
 differentiate :: (Ring.C a) => T a -> T a
 differentiate f =
    f{polynomial =
-        Bell.exponentPolynomial (bell f) * polynomial f +
-        Differential.differentiate (polynomial f)}
+        Differential.differentiate (polynomial f)
+        - Differential.differentiate (Bell.exponentPolynomial (bell f))
+           * polynomial f}
+
+{-
+snd $ integrate $ differentiate (Cons Bell.unit (Poly.fromCoeffs [7,7,7,7]) :: T Double)
+-}
+integrate ::
+   (Field.C a, ZeroTestable.C a) =>
+   T a -> (Complex.T a, T a)
+integrate f =
+   let fs = Poly.coeffs $ polynomial f
+       (ys,~[r]) =
+          Poly.divModRev
+             {-
+             We need the shortening convention of 'zipWith'
+             in order to limit the result list,
+             we cannot use list instance for (-).
+             -}
+             (zipWith (-)
+                (0 : 0 : diffRev ys)
+                (List.reverse fs))
+             (List.reverse $ Poly.coeffs $
+              Differential.differentiate $
+              Bell.exponentPolynomial $ bell f)
+   in  forcePair $
+       if null fs
+         then (zero, f)
+         else (r, f{polynomial = Poly.fromCoeffs $ List.reverse ys})
+
+diffRev :: Ring.C a => [a] -> [a]
+diffRev xs =
+   zipWith (*) xs
+      (drop 1 (iterate (subtract 1) (fromIntegral $ length xs)))
+
+translate :: Ring.C a => a -> T a -> T a
+translate d =
+   translateComplex (Complex.fromReal d)
+
+translateComplex :: Ring.C a => Complex.T a -> T a -> T a
+translateComplex d f =
+   Cons
+      (Bell.translateComplex d $ bell f)
+      (Poly.translate d $ polynomial f)
+
+modulate :: Ring.C a => a -> T a -> T a
+modulate d f =
+   Cons
+      (Bell.modulate d $ bell f)
+      (polynomial f)
+
+turn :: Ring.C a => a -> T a -> T a
+turn d f =
+   Cons
+      (Bell.turn d $ bell f)
+      (polynomial f)
+
+reverse :: Additive.C a => T a -> T a
+reverse f =
+   Cons
+      (Bell.reverse $ bell f)
+      (Poly.reverse $ polynomial f)
+
+dilate :: Field.C a => a -> T a -> T a
+dilate k f =
+   Cons
+      (Bell.dilate k $ bell f)
+      (Poly.dilate (Complex.fromReal k) $ polynomial f)
+
+shrink :: Ring.C a => a -> T a -> T a
+shrink k f =
+   Cons
+      (Bell.shrink k $ bell f)
+      (Poly.shrink (Complex.fromReal k) $ polynomial f)
+
+{-
+We could also amplify the polynomial coefficients.
+-}
+amplify :: Ring.C a => a -> T a -> T a
+amplify k f =
+   Cons
+      (Bell.amplify k $ bell f)
+      (polynomial f)
+
+
+{- |
+Approximate a @T a@ using a linear combination of translated @Bell.T a@.
+The smaller the unit (e.g. 0.1, 0.01, 0.001)
+the better the approximation but the worse the numeric properties.
+
+We cannot put all information into @amp@ of @Bell@,
+since @amp@ must be real, but is complex here by construction.
+We really need at least signed amplitudes at this place,
+since we want to represent differences of Gaussians.
+-}
+approximateByBells ::
+   Field.C a =>
+   a -> T a -> [(Complex.T a, Bell.T a)]
+approximateByBells unit f =
+   let b = bell f
+       amps =
+          -- approximateByBellsByTranslation
+          approximateByBellsAtOnce
+             unit
+             (Complex.scale (recip (2 * Bell.c2 b)) (Bell.c1 b))
+             (recip (2*unit*Bell.c2 b))
+             (polynomial f)
+   in  zip (LPoly.coeffs amps) $
+       map
+          (\d -> Bell.translate d b)
+          (laurentAbscissas (unit/2) amps)
+
+approximateByBellsAtOnce ::
+   Field.C a =>
+   a -> Complex.T a -> a -> Poly.T (Complex.T a) -> LPoly.T (Complex.T a)
+approximateByBellsAtOnce unit d s p =
+   foldr
+      (\x amps0 ->
+         {-
+         Decompose (bell t * (t-d)) = bell t * t - bell t * d
+         -}
+         let y = fmap (Complex.scale s) amps0
+         in  -- \t -> bell t * t
+             --    ~   (translate unit bell - translate (-unit) bell) / unit
+             LPoly.shift 1 y -
+             LPoly.shift (-1) y +
+             -- bell t * d
+             zipWithAbscissas
+                (\t z -> (Complex.fromReal t - d) * z)
+                (unit/2) amps0 +
+             LPoly.const x)
+      (LPoly.fromCoeffs [])
+      (Poly.coeffs p)
+
+approximateByBellsByTranslation ::
+   Field.C a =>
+   a -> Complex.T a -> a -> Poly.T (Complex.T a) -> LPoly.T (Complex.T a)
+approximateByBellsByTranslation unit d s p =
+   foldr
+      (\x amps0 ->
+         {-
+         Decompose (bell t * (t-d)) = bell t * t - bell t * d
+         -}
+         let y = fmap (Complex.scale s) amps0
+         in  -- \t -> bell t * t
+             --    ~   (translate unit bell - translate (-unit) bell) / unit
+             LPoly.shift 1 y -
+             LPoly.shift (-1) y +
+             -- bell t * d
+             zipWithAbscissas Complex.scale (unit/2) amps0 +
+             LPoly.const x)
+      (LPoly.fromCoeffs [])
+      (Poly.coeffs $ Poly.translate d p)
+
+zipWithAbscissas ::
+   (Ring.C a) =>
+   (a -> b -> c) -> a -> LPoly.T b -> LPoly.T c
+zipWithAbscissas h unit y =
+   LPoly.fromShiftCoeffs (LPoly.expon y) $
+   zipWith h
+      (laurentAbscissas unit y)
+      (LPoly.coeffs y)
+
+laurentAbscissas :: Ring.C a => a -> LPoly.T c -> [a]
+laurentAbscissas unit =
+   map (\d -> fromIntegral d * unit) .
+   iterate (1+) . LPoly.expon
+
 
 {- No Ring instance for Gaussians
 instance (Ring.C a) => Differential.C (T a) where
diff --git a/src/MathObj/Gaussian/Variance.hs b/src/MathObj/Gaussian/Variance.hs
--- a/src/MathObj/Gaussian/Variance.hs
+++ b/src/MathObj/Gaussian/Variance.hs
@@ -1,12 +1,7 @@
 {-# LANGUAGE NoImplicitPrelude #-}
 {-
-Reciprocal of variance of a Gaussian bell curve.
-We describe the curve only in terms of its variance
-thus we represent a bell curve at the coordinate origin
-neglecting its amplitude.
-
-We could also define the amplitude as @root 4 c@,
-but then @dilate@ and @shrink@ also include an amplification.
+We represent a Gaussian bell curve in terms of the reciprocal of its variance
+and its value at the origin.
 
 We could do some projective geometry in the exponent
 in order to also have zero variance,
@@ -23,34 +18,38 @@
 import qualified Algebra.Ring           as Ring
 import qualified Algebra.Additive       as Additive
 
+{-
 import Algebra.Transcendental (pi, )
 import Algebra.Ring ((*), (^), )
 import Algebra.Additive ((+))
-
+-}
 import Test.QuickCheck (Arbitrary, arbitrary, coarbitrary, )
-
+import Control.Monad (liftM2, )
 
 -- import Prelude (($))
 import NumericPrelude
 import PreludeBase
 
 
-data T a = Cons {c :: a}
+data T a = Cons {amp, c :: a}
    deriving (Eq, Show)
 
 instance (Real.C a, Arbitrary a) => Arbitrary (T a) where
-   arbitrary = fmap (Cons . (1+) . abs) arbitrary
+   arbitrary =
+      liftM2 Cons
+         arbitrary
+         (fmap ((1+) . abs) arbitrary)
    coarbitrary = undefined
 
 
-constant :: Additive.C a => T a
-constant = Cons zero
+constant :: Ring.C a => T a
+constant = Cons one zero
 
 {-# INLINE evaluate #-}
 evaluate :: (Trans.C a) =>
    T a -> a -> a
 evaluate f x =
-   exp $ (-pi * c f * x^2)
+   sqrt (amp f) * exp (-pi * c f * x^2)
 
 exponentPolynomial :: (Additive.C a) =>
    T a -> Poly.T a
@@ -58,17 +57,17 @@
    Poly.fromCoeffs [zero, zero, c f]
 
 
-norm1 :: (Algebraic.C a) => T a -> a
+norm1 :: (Algebraic.C a, Real.C a) => T a -> a
 norm1 f =
-   recip $ sqrt $ c f
+   sqrt $ abs (amp f) / c f
 
-norm2 :: (Algebraic.C a) => T a -> a
+norm2 :: (Algebraic.C a, Real.C a) => T a -> a
 norm2 f =
-   recip $ sqrt $ sqrt $ 2 * c f
+   sqrt $ abs (amp f) / (sqrt $ 2 * c f)
 
-normP :: (Trans.C a) => a -> T a -> a
+normP :: (Trans.C a, Real.C a) => a -> T a -> a
 normP p f =
-   (p * c f) ^? (- recip (2*p))
+   sqrt (abs (amp f)) * (p * c f) ^? (- recip (2*p))
 
 
 variance :: (Trans.C a) =>
@@ -76,10 +75,10 @@
 variance f =
    recip $ c f * 2*pi
 
-multiply :: (Additive.C a) =>
+multiply :: (Ring.C a) =>
    T a -> T a -> T a
 multiply f g =
-   Cons $ c f + c g
+   Cons (amp f * amp g) (c f + c g)
 
 {- |
 > convolve x y t =
@@ -88,7 +87,9 @@
 convolve :: (Field.C a) =>
    T a -> T a -> T a
 convolve f g =
-   Cons $ recip $ recip (c f) + recip (c g)
+   Cons
+      (amp f * amp g / (c f + c g))
+      (recip $ recip (c f) + recip (c g))
 
 {- |
 > fourier x f =
@@ -97,18 +98,23 @@
 fourier :: (Field.C a) =>
    T a -> T a
 fourier f =
-   Cons $ recip $ c f
+   Cons (amp f / c f) (recip $ c f)
 {-
 fourier (t -> exp(-(a*t)^2))
 -}
 
 dilate :: (Field.C a) => a -> T a -> T a
 dilate k f =
-   Cons $ c f / k^2
+   Cons (amp f) $ c f / k^2
 
 shrink :: (Ring.C a) => a -> T a -> T a
 shrink k f =
-   Cons $ c f * k^2
+   Cons (amp f) $ c f * k^2
+
+amplify :: (Ring.C a) => a -> T a -> T a
+amplify k f =
+   Cons (k^2 * amp f) $ c f
+
 
 {- laws
 fourier (convolve f g) = multiply (fourier f) (fourier g)
diff --git a/src/MathObj/LaurentPolynomial.hs b/src/MathObj/LaurentPolynomial.hs
--- a/src/MathObj/LaurentPolynomial.hs
+++ b/src/MathObj/LaurentPolynomial.hs
@@ -69,8 +69,12 @@
 bounds :: T a -> (Int, Int)
 bounds (Cons xt x) = (xt, xt + length x - 1)
 
+shift :: Int -> T a -> T a
+shift t (Cons xt x) = Cons (xt+t) x
+
+{-# DEPRECATED translate "In order to avoid confusion with Polynomial.translate, use 'shift' instead" #-}
 translate :: Int -> T a -> T a
-translate t (Cons xt x) = Cons (xt+t) x
+translate = shift
 
 
 instance Functor T where
diff --git a/src/MathObj/Polynomial.hs b/src/MathObj/Polynomial.hs
--- a/src/MathObj/Polynomial.hs
+++ b/src/MathObj/Polynomial.hs
@@ -50,11 +50,12 @@
     compose, equal, add, sub, negate,
     horner, hornerCoeffVector, hornerArgVector,
     shift, unShift,
-    mul, scale, divMod,
+    mul, scale, divMod, divModRev,
     tensorProduct, tensorProductAlt,
     mulShear, mulShearTranspose,
     progression, differentiate, integrate, integrateInt,
-    fromRoots, alternate, reverse, )
+    fromRoots, alternate, reverse,
+    translate, dilate, shrink, )
 where
 
 import qualified Algebra.Differential         as Differential
@@ -76,8 +77,9 @@
 import Control.Monad (liftM, )
 import qualified Data.List as List
 import NumericPrelude.List (zipWithOverlap, )
+import Data.Tuple.HT (mapPair, mapFst, forcePair, )
 import Data.List.HT
-          (dropWhileRev, shear, shearTranspose, outerProduct, )
+          (dropWhileRev, switchL, shear, shearTranspose, outerProduct, )
 
 import Test.QuickCheck (Arbitrary(arbitrary,coarbitrary))
 
@@ -307,19 +309,27 @@
 
 divMod :: (ZeroTestable.C a, Field.C a) => [a] -> [a] -> ([a], [a])
 divMod x y =
-    let (y0:ys) = dropWhile isZero (List.reverse y)
-        aux l xs' =
-          if l < 0
-            then ([], xs')
-            else
+   mapPair (List.reverse, List.reverse) $
+   divModRev (List.reverse x) (List.reverse y)
+
+{-
+snd $ Poly.divMod (repeat (1::Double)) [1,1]
+-}
+divModRev :: (ZeroTestable.C a, Field.C a) => [a] -> [a] -> ([a], [a])
+divModRev x y =
+   let (y0:ys) = dropWhile isZero y
+       -- the second parameter represents lazily (length x - length y)
+       aux xs' =
+         forcePair .
+         switchL
+           ([], xs')
+           (P.const $
               let (x0:xs) = xs'
                   q0      = x0/y0
-                  (d',m') = aux (l-1) (sub xs (scale q0 ys))
-              in  (q0:d',m')
-        (d, m) = aux (length x - length y) (List.reverse x)
-    in  if isZero y
-          then error "MathObj.Polynomial: division by zero"
-          else (List.reverse d, List.reverse m)
+              in  mapFst (q0:) . aux (sub xs (scale q0 ys)))
+   in  if isZero y
+         then error "MathObj.Polynomial: division by zero"
+         else aux x (drop (length y - 1) x)
 
 instance (ZeroTestable.C a, Field.C a) => Integral.C (T a) where
   divMod (Cons x) (Cons y) =
@@ -351,7 +361,7 @@
 
 {-# INLINE differentiate #-}
 differentiate :: (Ring.C a) => [a] -> [a]
-differentiate = zipWith (*) progression . tail
+differentiate = zipWith (*) progression . drop 1
 
 {-# INLINE integrate #-}
 integrate :: (Field.C a) => a -> [a] -> [a]
@@ -388,6 +398,17 @@
 {-# INLINE reverse #-}
 reverse :: Additive.C a => T a -> T a
 reverse = lift1 alternate
+
+translate :: Ring.C a => a -> T a -> T a
+translate d =
+   lift1 $ foldr (\c p -> [c] + mulLinearFactor d p) []
+
+shrink :: Ring.C a => a -> T a -> T a
+shrink k =
+   lift1 $ zipWith (*) (iterate (k*) one)
+
+dilate :: Field.C a => a -> T a -> T a
+dilate = shrink . Field.recip
 
 {-
 see htam: Wavelet/DyadicResultant
diff --git a/src/Number/Complex.hs b/src/Number/Complex.hs
--- a/src/Number/Complex.hs
+++ b/src/Number/Complex.hs
@@ -34,6 +34,7 @@
         signum,
         toPolar,
         magnitude,
+        magnitudeSqr,
         phase,
         -- * Conjugate
         conjugate,
@@ -69,9 +70,16 @@
 import qualified Algebra.Indexable          as Indexable
 
 import Algebra.ZeroTestable(isZero)
-import Algebra.Module((*>))
-import Algebra.Algebraic((^/))
+import Algebra.Module((*>), (<*>.*>), )
+import Algebra.Algebraic((^/), )
 
+import qualified NumericPrelude.Elementwise as Elem
+import Algebra.Additive ((<*>.+), (<*>.-), (<*>.-$), )
+
+import Foreign.Storable (Storable (..), )
+import qualified Foreign.Storable.Record as Store
+import Control.Applicative (liftA2, )
+
 import Test.QuickCheck (Arbitrary, arbitrary, coarbitrary)
 import Control.Monad (liftM2)
 
@@ -120,7 +128,23 @@
    {-# INLINE coarbitrary #-}
    coarbitrary = undefined
 
+instance (Storable a) => Storable (T a) where
+   sizeOf    = Store.sizeOf store
+   alignment = Store.alignment store
+   peek      = Store.peek store
+   poke      = Store.poke store
 
+store ::
+   (Storable a) =>
+   Store.Dictionary (T a)
+store =
+   Store.run $
+   liftA2 (+:)
+      (Store.element real)
+      (Store.element imag)
+
+
+
 {- * Functions -}
 
 -- | Construct a complex number from real and imaginary part.
@@ -265,13 +289,13 @@
     {-# SPECULATE instance Additive.C (T Float) #-}
     {-# SPECULATE instance Additive.C (T Double) #-}
     {-# INLINE zero #-}
-    zero                        =  Cons zero zero
+    {-# INLINE negate #-}
     {-# INLINE (+) #-}
-    (Cons x y) + (Cons x' y')   =  Cons (x+x') (y+y')
     {-# INLINE (-) #-}
-    (Cons x y) - (Cons x' y')   =  Cons (x-x') (y-y')
-    {-# INLINE negate #-}
-    negate (Cons x y)           =  Cons (negate x) (negate y)
+    zero   = Cons zero zero
+    (+)    = Elem.run2 $ Elem.with Cons <*>.+  real <*>.+  imag
+    (-)    = Elem.run2 $ Elem.with Cons <*>.-  real <*>.-  imag
+    negate = Elem.run  $ Elem.with Cons <*>.-$ real <*>.-$ imag
 
 instance  (Ring.C a) => Ring.C (T a)  where
     {-# SPECULATE instance Ring.C (T Float) #-}
@@ -295,7 +319,8 @@
 --   because it requires the Algebra.Module constraint
 instance (Module.C a b) => Module.C a (T b) where
    {-# INLINE (*>) #-}
-   s *> (Cons x y)  = Cons (s *> x) (s *> y)
+   (*>) = Elem.run2 $ Elem.with Cons <*>.*> real <*>.*> imag
+   -- s *> (Cons x y)  = Cons (s *> x) (s *> y)
 
 instance (VectorSpace.C a b) => VectorSpace.C a (T b)
 
diff --git a/src/Number/Quaternion.hs b/src/Number/Quaternion.hs
--- a/src/Number/Quaternion.hs
+++ b/src/Number/Quaternion.hs
@@ -47,13 +47,15 @@
 import qualified Algebra.ZeroTestable as ZeroTestable
 
 import Algebra.ZeroTestable(isZero)
-import Algebra.Module((*>))
--- import Algebra.Algebraic((^/))
+import Algebra.Module((*>), (<*>.*>), )
 
 import qualified Number.Complex as Complex
 
 import Number.Complex ((+:))
 
+import qualified NumericPrelude.Elementwise as Elem
+import Algebra.Additive ((<*>.+), (<*>.-), (<*>.-$), )
+
 -- import qualified Data.Typeable as Ty
 import Data.Array (Array, (!))
 import qualified Data.Array as Array
@@ -256,10 +258,10 @@
 instance (Additive.C a) => Additive.C (T a)  where
    {-# SPECIALISE instance Additive.C (T Float) #-}
    {-# SPECIALISE instance Additive.C (T Double) #-}
-   zero			=  Cons zero zero
-   (Cons xr xi) + (Cons yr yi)	=  Cons (xr+yr) (xi+yi)
-   (Cons xr xi) - (Cons yr yi)	=  Cons (xr-yr) (xi-yi)
-   negate (Cons x y)		=  Cons (negate x) (negate y)
+   zero   = Cons zero zero
+   (+)    = Elem.run2 $ Elem.with Cons <*>.+  real <*>.+  imag
+   (-)    = Elem.run2 $ Elem.with Cons <*>.-  real <*>.-  imag
+   negate = Elem.run  $ Elem.with Cons <*>.-$ real <*>.-$ imag
 
 instance (Ring.C a) => Ring.C (T a)  where
    {-# SPECIALISE instance Ring.C (T Float) #-}
@@ -288,7 +290,7 @@
 -- | The '(*>)' method can't replace 'scale'
 --   because it requires the Algebra.Module constraint
 instance (Module.C a b) => Module.C a (T b) where
-   s *> (Cons r i)  = Cons (s *> r) (s *> i)
+   (*>) = Elem.run2 $ Elem.with Cons <*>.*> real <*>.*> imag
 
 instance (VectorSpace.C a b) => VectorSpace.C a (T b)
 
diff --git a/src/Number/Ratio.hs b/src/Number/Ratio.hs
--- a/src/Number/Ratio.hs
+++ b/src/Number/Ratio.hs
@@ -40,6 +40,10 @@
 
 import Control.Monad(liftM, liftM2, )
 
+import Foreign.Storable (Storable (..), )
+import qualified Foreign.Storable.Record as Store
+import Control.Applicative (liftA2, )
+
 import Test.QuickCheck (Arbitrary(arbitrary,coarbitrary))
 import System.Random (Random(..), RandomGen, )
 
@@ -157,6 +161,22 @@
       liftM2 (%) arbitrary
          (liftM (\x -> if isZero x then one else x) arbitrary)
    coarbitrary = undefined
+
+
+instance (Storable a, PID.C a) => Storable (T a) where
+   sizeOf    = Store.sizeOf store
+   alignment = Store.alignment store
+   peek      = Store.peek store
+   poke      = Store.poke store
+
+store ::
+   (Storable a, PID.C a) =>
+   Store.Dictionary (T a)
+store =
+   Store.run $
+   liftA2 (%)
+      (Store.element numerator)
+      (Store.element denominator)
 
 {-
 This instance may not be appropriate for mathematical objects other than numbers.
diff --git a/src/NumericPrelude/Elementwise.hs b/src/NumericPrelude/Elementwise.hs
new file mode 100644
--- /dev/null
+++ b/src/NumericPrelude/Elementwise.hs
@@ -0,0 +1,41 @@
+module NumericPrelude.Elementwise where
+
+import Control.Applicative (Applicative(pure, (<*>)), )
+
+{- |
+A reader monad for the special purpose
+of defining instances of certain operations on tuples and records.
+It does not add any new functionality to the common Reader monad,
+but it restricts the functions to the required ones
+and exports them from one module.
+That is you do not have to import
+both Control.Monad.Trans.Reader and Control.Applicative.
+The type also tells the user, for what the Reader monad is used.
+We can more easily replace or extend the implementation when needed.
+-}
+newtype T v a = Cons {run :: v -> a}
+
+{-# INLINE with #-}
+with :: a -> T v a
+with e = Cons $ \ _v -> e
+
+{-# INLINE element #-}
+element :: (v -> a) -> T v a
+element = Cons
+
+
+{-# INLINE run2 #-}
+run2 :: T (x,y) a -> x -> y -> a
+run2 = curry . run
+
+instance Functor (T v) where
+   {-# INLINE fmap #-}
+   fmap f (Cons e) =
+      Cons $ \v -> f $ e v
+
+instance Applicative (T v) where
+   {-# INLINE pure #-}
+   {-# INLINE (<*>) #-}
+   pure = with
+   Cons f <*> Cons e =
+      Cons $ \v -> f v $ e v
diff --git a/test/Gaussian.hs b/test/Gaussian.hs
new file mode 100644
--- /dev/null
+++ b/test/Gaussian.hs
@@ -0,0 +1,6 @@
+module Main where
+
+import qualified MathObj.Gaussian.Example as Example
+
+main :: IO ()
+main = Example.polyApprox
diff --git a/test/Test/MathObj/Gaussian/Bell.hs b/test/Test/MathObj/Gaussian/Bell.hs
--- a/test/Test/MathObj/Gaussian/Bell.hs
+++ b/test/Test/MathObj/Gaussian/Bell.hs
@@ -62,6 +62,8 @@
           simple $ \x -> nest 2 G.fourier x == G.reverse x) :
       ("reverse identity",
           simple $ \x -> nest 2 G.reverse x == x) :
+      ("fourier unit",
+          quickCheck $ G.fourier G.unit == (G.unit :: G.T Rational)) :
       ("translate additive",
           simple $ \x a b ->
              G.translate a (G.translate b x) == G.translate (a+b) x) :
@@ -90,5 +92,5 @@
              G.shrink a x == G.dilate (recip a) x) :
       ("fourier dilate",
           simple $ \x a -> a>0 ==>
-             G.fourier (G.dilate a x) == G.shrink a (G.fourier x)) :
+             G.fourier (G.dilate a x) == G.amplify a (G.shrink a (G.fourier x))) :
       []
diff --git a/test/Test/MathObj/Gaussian/Variance.hs b/test/Test/MathObj/Gaussian/Variance.hs
--- a/test/Test/MathObj/Gaussian/Variance.hs
+++ b/test/Test/MathObj/Gaussian/Variance.hs
@@ -58,5 +58,5 @@
              G.shrink a x == G.dilate (recip a) x) :
       ("fourier dilate",
           simple $ \x a -> a>0 ==>
-             G.fourier (G.dilate a x) == G.shrink a (G.fourier x)) :
+             G.fourier (G.dilate a x) == G.amplify a (G.shrink a (G.fourier x))) :
       []
diff --git a/test/Test/Run.hs b/test/Test/Run.hs
--- a/test/Test/Run.hs
+++ b/test/Test/Run.hs
@@ -1,5 +1,6 @@
 module Main where
 
+import qualified Test.MathObj.Gaussian.Polynomial as GaussPoly
 import qualified Test.MathObj.Gaussian.Variance as GaussVariance
 import qualified Test.MathObj.Gaussian.Bell as GaussBell
 import qualified Test.MathObj.PartialFraction as PartialFraction
@@ -13,6 +14,7 @@
    do HUnitText.runTestTT (HUnit.TestList $
          GaussVariance.tests :
          GaussBell.tests :
+         GaussPoly.tests :
          PartialFraction.tests :
          Polynomial.tests :
          PowerSeries.tests :