Geometric algebra split in even and odd subalgebras

geoalg.hs

{-# LANGUAGE
    GADTs
  , DataKinds
  , InstanceSigs
  , ViewPatterns
  , TypeFamilies
  , TypeOperators
  , TypeApplications
  , FlexibleInstances
  , StandaloneDeriving
  , AllowAmbiguousTypes
  , ScopedTypeVariables
  , UndecidableInstances
  #-}

import Prelude hiding (or, reverse)

import Linear.Epsilon
import Linear.Metric
import Linear.Quaternion (Quaternion(..))
import Linear.Vector hiding (el)
import Linear.V hiding (int)
import Linear.V3 (V3(..))
import qualified Data.Vector as V

import Control.Lens.Iso
import Data.Complex (Complex(..))

import Data.List (sortOn)
import Data.Monoid ((<>), Sum(..))
import GHC.TypeLits (type (*))


data Nat = Z | S Nat
type D1 = 'S 'Z
type D2 = 'S D1
type D3 = 'S D2
type D4 = 'S D3

-- ! A multivector of dimension `d` over field `a`
data MV d a = MV !(MVE d a) !(MVO d a)

-- | The even grades of a multivector
data MVE (d :: Nat) a where
  Sc :: !a -> MVE d a                             -- scalar
  ZE :: MVE d a                                   -- zero
  BE :: !(MVE d a) -> !(MVO d a) -> MVE ('S d) a  -- BE a b == a + b * e_(d+1)

-- | The odd grades of a multivector
data MVO (d :: Nat) a where
  ZO :: MVO d a                                   -- zero
  BO :: !(MVO d a) -> !(MVE d a) -> MVO ('S d) a  -- BO a b == a + b * e_(d+1)


mkSc :: (Num a, Eq a) => a -> MVE d a
mkSc 0 = ZE
mkSc a = Sc a
mkBE :: MVE d a -> MVO d a -> MVE ('S d) a
mkBE ZE ZO = ZE
mkBE (Sc a) ZO = Sc a
mkBE e o = BE e o
mkBO :: MVO d a -> MVE d a -> MVO ('S d) a
mkBO ZO ZE = ZO
mkBO o e = BO o e

isBlade :: (Num a, Eq a, IsDim d) => MV d a -> Bool
isBlade = (<= (1 :: Int)) . getSum . foldMap (Sum . const 1) . mkSparse

liftDim :: (Num a, Eq a) => MV d a -> MV ('S d) a
liftDim (MV e o) = MV (mkBE e ZO) (mkBO o ZE)


grade :: (Num a, Eq a) => Int -> MV d a -> MV d a
grade d (MV mve mvo) | even d    = MV (gradeE d mve) ZO
                     | otherwise = MV ZE (gradeO d mvo)
  where
    gradeE :: (Num a, Eq a) => Int -> MVE d a -> MVE d a
    gradeE 0 (Sc a) = mkSc a
    gradeE g (BE e o) = mkBE (gradeE g e) (gradeO (g - 1) o)
    gradeE _ _ = ZE
    gradeO :: (Num a, Eq a) => Int -> MVO d a -> MVO d a
    gradeO g (BO o e) = mkBO (gradeO g o) (gradeE (g - 1) e)
    gradeO _ _ = ZO

reverse :: (Eq a, Floating a) => MV d a -> MV d a
reverse (MV mve mvo) = MV (reverseE mve) (reverseO mvo)
  where
    reverseE :: (Eq a, Floating a) => MVE d a -> MVE d a
    reverseE (BE e o) = BE (reverseE e) (-reverseO o)
    reverseE e = e
    reverseO :: (Eq a, Floating a) => MVO d a -> MVO d a
    reverseO (BO o e) = BO (reverseO o) (reverseE e)
    reverseO o = o

involute :: (Eq a, Floating a) => MV d a -> MV d a
involute (MV e o) = MV e (-o)

cliffordConjugate :: (Eq a, Floating a) => MV d a -> MV d a
cliffordConjugate = involute . reverse


mult :: (Eq a, Floating a) => MV d a -> MV d a -> MV d a
mult (MV mel mol) (MV mer mor) = MV (multEE mel mer + multOO mol mor) (multEO mel mor + multOE mol mer)

multEE :: (Eq a, Floating a) => MVE d a -> MVE d a -> MVE d a
multEE ZE _ = ZE
multEE _ ZE = ZE
multEE (Sc a) v = a *^ v
multEE v (Sc a) = v ^* a
multEE (BE el ol) (BE er or) = mkBE (multEE el er - multOO ol or) (multEO el or + multOE ol er)
multEO :: (Eq a, Floating a) => MVE d a -> MVO d a -> MVO d a
multEO ZE _ = ZO
multEO _ ZO = ZO
multEO (Sc a) v = a *^ v
multEO (BE el ol) (BO or er) = mkBO (multEO el or + multOE ol er) (multEE el er - multOO ol or)
multOE :: (Eq a, Floating a) => MVO d a -> MVE d a -> MVO d a
multOE ZO _ = ZO
multOE _ ZE = ZO
multOE v (Sc a) = v ^* a
multOE (BO ol el) (BE er or) = mkBO (multOE ol er - multEO el or) (multOO ol or + multEE el er)
multOO :: (Eq a, Floating a) => MVO d a -> MVO d a -> MVE d a
multOO ZO _ = ZE
multOO _ ZO = ZE
multOO (BO ol el) (BO or er) = mkBE (multOO ol or + multEE el er) (multOE ol er - multEO el or)

filteredMult :: (Eq a, Floating a) => (Int -> Int -> Int) -> MV d a -> MV d a -> MV d a
filteredMult op (MV mel mol) (MV mer mor) = MV (fmultEE 0 0 0 mel mer + fmultOO 0 0 0 mol mor) (fmultEO 0 0 0 mel mor + fmultOE 0 0 0 mol mer)
  where
    fmultEE :: (Eq a, Floating a) => Int -> Int -> Int -> MVE d a -> MVE d a -> MVE d a
    fmultEE _ _ _ ZE _ = ZE
    fmultEE _ _ _ _ ZE = ZE
    fmultEE r s o (Sc a) (Sc b) = if r `op` s == o then mkSc (a * b) else ZE
    fmultEE r s o (Sc a) (BE er or) = mkBE (fmultEE r s o (Sc a) er) (fmultEO r (s + 1) (o + 1) (Sc a) or)
    fmultEE r s o (BE el ol) (Sc a) = mkBE (fmultEE r s o el (Sc a)) (fmultOE (r + 1) s (o + 1) ol (Sc a))
    fmultEE r s o (BE el ol) (BE er or) =
      mkBE (fmultEE r s o el er - fmultOO (r + 1) (s + 1) o ol or) (fmultEO r (s + 1) (o + 1) el or + fmultOE (r + 1) s (o + 1) ol er)
    fmultEO :: (Eq a, Floating a) => Int -> Int -> Int -> MVE d a -> MVO d a -> MVO d a
    fmultEO _ _ _ ZE _ = ZO
    fmultEO _ _ _ _ ZO = ZO
    fmultEO r s o (Sc a) (BO or er) = mkBO (fmultEO r s o (Sc a) or) (fmultEE r (s + 1) (o + 1) (Sc a) er)
    fmultEO r s o (BE el ol) (BO or er) =
      mkBO (fmultEO r s o el or + fmultOE (r + 1) (s + 1) o ol er) (fmultEE r (s + 1) (o + 1) el er - fmultOO (r + 1) s (o + 1) ol or)
    fmultOE :: (Eq a, Floating a) => Int -> Int -> Int -> MVO d a -> MVE d a -> MVO d a
    fmultOE _ _ _ ZO _ = ZO
    fmultOE _ _ _ _ ZE = ZO
    fmultOE r s o (BO ol el) (Sc a) = mkBO (fmultOE r s o ol (Sc a)) (fmultEE (r + 1) s (o + 1) el (Sc a))
    fmultOE r s o (BO ol el) (BE er or) =
      mkBO (fmultOE r s o ol er - fmultEO (r + 1) (s + 1) o el or) (fmultOO r (s + 1) (o + 1) ol or + fmultEE (r + 1) s (o + 1) el er)
    fmultOO :: (Eq a, Floating a) => Int -> Int -> Int -> MVO d a -> MVO d a -> MVE d a
    fmultOO _ _ _ ZO _ = ZE
    fmultOO _ _ _ _ ZO = ZE
    fmultOO r s o (BO ol el) (BO or er) =
      mkBE (fmultOO r s o ol or + fmultEE (r + 1) (s + 1) o el er) (fmultOE r (s + 1) (o + 1) ol er - fmultEO (r + 1) s (o + 1) el or)

(/\) :: (Eq a, Floating a) => MV d a -> MV d a -> MV d a
(/\) = filteredMult (\r s -> r + s)

contractL :: (Eq a, Floating a) => MV d a -> MV d a -> MV d a
contractL = filteredMult (\r s -> s - r)

contractR :: (Eq a, Floating a) => MV d a -> MV d a -> MV d a
contractR = filteredMult (\r s -> r - s)

dotProduct :: (Eq a, Floating a) => MV d a -> MV d a -> MV d a
dotProduct = filteredMult (\r s -> abs (r - s))

scalarProduct :: (Eq a, Floating a) => MV d a -> MV d a -> MV d a
scalarProduct = filteredMult (\_ _ -> 0)


instance Functor (MV d) where
  fmap f (MV e o) = MV (fmap f e) (fmap f o)
  {-# INLINE fmap #-}
  a <$ MV e o = MV (a <$ e) (a <$ o)
  {-# INLINE (<$) #-}
instance Functor (MVE d) where
  fmap f (Sc a) = Sc (f a)
  fmap _ ZE = ZE
  fmap f (BE e o) = BE (fmap f e) (fmap f o)
  {-# INLINE fmap #-}
  a <$ (Sc _) = Sc a
  _ <$ ZE = ZE
  a <$ BE e o = BE (a <$ e) (a <$ o)
  {-# INLINE (<$) #-}
instance Functor (MVO d) where
  fmap _ ZO = ZO
  fmap f (BO o e) = BO (fmap f o) (fmap f e)
  {-# INLINE fmap #-}
  _ <$ ZO = ZO
  a <$ BO o e = BO (a <$ o) (a <$ e)
  {-# INLINE (<$) #-}
instance Foldable (MV d) where
  foldMap f (MV e o) = foldMap f e <> foldMap f o
  {-# INLINE foldMap #-}
instance Foldable (MVE d) where
  foldMap f (Sc a) = f a
  foldMap _ ZE = mempty
  foldMap f (BE e o) = foldMap f e <> foldMap f o
  {-# INLINE foldMap #-}
instance Foldable (MVO d) where
  foldMap _ ZO = mempty
  foldMap f (BO o e) = foldMap f o <> foldMap f e
  {-# INLINE foldMap #-}
instance Traversable (MV d) where
  traverse f (MV e o) = MV <$> traverse f e <*> traverse f o
  {-# INLINE traverse #-}
instance Traversable (MVE d) where
  traverse f (Sc a) = Sc <$> f a
  traverse _ ZE = pure ZE
  traverse f (BE e o) = BE <$> traverse f e <*> traverse f o
  {-# INLINE traverse #-}
instance Traversable (MVO d) where
  traverse _ ZO = pure ZO
  traverse f (BO o e) = BO <$> traverse f o <*> traverse f e
  {-# INLINE traverse #-}


class IsDim (n :: Nat) where
  int :: Int
  showAsListE :: (Show a, Eq a, Num a) => MVE n a -> [(a, String)]
  showAsListO :: (Show a, Eq a, Num a) => MVO n a -> [(a, String)]
  toPseudo :: Num a => a -> MV ('S n) a
  eO :: Num a => Int -> MVO ('S n) a
  mkDense :: (Num a, Eq a) => MV n a -> MV n a
  mkSparse :: (Num a, Eq a) => MV n a -> MV n a
  mkDenseE :: (Num a, Eq a) => MVE n a -> MVE n a
  mkSparseE :: (Num a, Eq a) => MVE n a -> MVE n a
  mkDenseO :: (Num a, Eq a) => MVO n a -> MVO n a
  mkSparseO :: (Num a, Eq a) => MVO n a -> MVO n a
instance IsDim 'Z where
  int = 0
  showAsListE ZE = []
  showAsListE (Sc 0) = []
  showAsListE (Sc a) = [(a, "")]
  showAsListO ZO = []
  toPseudo o = MV ZE (BO ZO (Sc o))
  eO 1 = BO ZO (Sc 1)
  eO _ = ZO
  mkDense (MV e o) = MV (mkDenseE e) o
  mkSparse (MV e o) = MV (mkSparseE e) o
  mkDenseE ZE = Sc 0
  mkDenseE v = v
  mkSparseE (Sc a) = mkSc a
  mkSparseE ZE = ZE
  mkDenseO ZO = ZO
  mkSparseO ZO = ZO
instance IsDim n => IsDim ('S n) where
  int = 1 + int @n
  showAsListE ZE = []
  showAsListE (Sc 0) = []
  showAsListE (Sc a) = [(a, "")]
  showAsListE (BE e o) = showAsListE e ++ fmap (fmap (++ 'e':subscript (1 + int @n))) (showAsListO o)
  showAsListO ZO = []
  showAsListO (BO o e) = showAsListO o ++ fmap (fmap (++ 'e':subscript (1 + int @n))) (showAsListE e)
  eO i | int @n + 2 == i = BO ZO (Sc 1)
       | otherwise       = BO (eO i) ZE
  toPseudo a = case toPseudo a of
    MV mve ZO -> MV ZE (BO ZO mve)
    MV ZE mvo -> MV (BE ZE mvo) ZO
    _ -> error "unexpected toPseudo return value"
  mkDense (MV e o) = MV (mkDenseE e) (mkDenseO o)
  mkSparse (MV e o) = MV (mkSparseE e) (mkSparseO o)
  mkDenseE (Sc a) = BE (mkDenseE (Sc a)) (mkDenseO ZO)
  mkDenseE ZE = BE (mkDenseE ZE) (mkDenseO ZO)
  mkDenseE (BE e o) = BE (mkDenseE e) (mkDenseO o)
  mkSparseE (Sc a) = mkSc a
  mkSparseE ZE = ZE
  mkSparseE (BE e o) = mkBE (mkSparseE e) (mkSparseO o)
  mkDenseO ZO = BO (mkDenseO ZO) (mkDenseE ZE)
  mkDenseO (BO o e) = BO (mkDenseO o) (mkDenseE e)
  mkSparseO ZO = ZO
  mkSparseO (BO o e) = mkBO (mkSparseO o) (mkSparseE e)



e_ :: (Num a, IsDim d) => Int -> MV ('S d) a
e_ i = MV ZE (eO i)

e1, e2, e3 :: MV D3 Double
e1 = e_ 1
e2 = e_ 2
e3 = e_ 3

qi, qj, qk :: MV D3 Double
qi = e3 * e2
qj = e1 * e3
qk = e2 * e1

complexIso :: (Eq a, Num a) => Iso' (Complex a) (MVE D2 a)
complexIso = iso to from
  where
    to :: (Eq a, Num a) => Complex a -> MVE D2 a
    to (a :+ b) = mkBE (mkSc a) (mkBO ZO (mkSc b))
    from :: (Eq a, Num a) => MVE D2 a -> Complex a
    from (mkDenseE -> BE (BE (Sc a) ZO) (BO ZO (Sc b))) = a :+ b
    from _ = error "not possible after mkDenseE"

quaternionIso :: (Eq a, Num a) => Iso' (Quaternion a) (MVE D3 a)
quaternionIso = iso to from
  where
    to :: (Eq a, Num a) => Quaternion a -> MVE D3 a
    to (Quaternion a (V3 b c d)) = mkBE (mkBE (mkSc a) (BO ZO (mkSc b))) (mkBO (mkBO ZO (mkSc c)) (mkSc d))
    from :: (Eq a, Num a) => MVE D3 a -> Quaternion a
    from (mkDenseE -> BE (BE (BE (Sc a) ZO) (BO ZO (Sc b))) (BO (BO ZO (Sc c)) (BE (Sc d) ZO))) = Quaternion a (V3 b c d)
    from _ = error "not possible after mkDenseE"


instance Finite (MV 'Z) where
  type Size (MV 'Z) = 1
  fromV (V v) = MV (Sc (v V.! 0)) ZO
instance Finite (MV d) => Finite (MV ('S d)) where
  type Size (MV ('S d)) = 2 * Size (MV d)
  fromV (V v) = MV (BE el ol) (BO or er)
    where
      MV el or = fromV (V v1)
      MV er ol = fromV (V v2)
      (v1, v2) = V.splitAt (V.length v `div` 2) v

instance Additive (MV d) where
  zero = MV ZE ZO
  {-# INLINE zero #-}
  liftU2 f (MV ae ao) (MV be bo) = MV (liftU2 f ae be) (liftU2 f ao bo)
  {-# INLINE liftU2 #-}
  liftI2 f (MV ae ao) (MV be bo) = MV (liftI2 f ae be) (liftI2 f ao bo)
  {-# INLINE liftI2 #-}
instance Additive (MVE d) where
  zero = ZE
  {-# INLINE zero #-}
  liftU2 _ ZE v = v
  liftU2 _ v ZE = v
  liftU2 f (Sc a) (Sc b) = Sc (f a b)
  liftU2 f (Sc a) (BE be bo) = BE (liftU2 f (Sc a) be) bo
  liftU2 f (BE ae ao) (Sc b) = BE (liftU2 f ae (Sc b)) ao
  liftU2 f (BE ae ao) (BE be bo) = BE (liftU2 f ae be) (liftU2 f ao bo)
  {-# INLINE liftU2 #-}
  liftI2 _ ZE _ = ZE
  liftI2 _ _ ZE = ZE
  liftI2 f (Sc a) (Sc b) = Sc (f a b)
  liftI2 f (Sc a) (BE be _) = BE (liftI2 f (Sc a) be) ZO
  liftI2 f (BE ae _) (Sc b) = BE (liftI2 f ae (Sc b)) ZO
  liftI2 f (BE ae ao) (BE be bo) = BE (liftI2 f ae be) (liftI2 f ao bo)
  {-# INLINE liftI2 #-}
instance Additive (MVO d) where
  zero = ZO
  {-# INLINE zero #-}
  liftU2 _ ZO v = v
  liftU2 _ v ZO = v
  liftU2 f (BO ao ae) (BO bo be) = BO (liftU2 f ao bo) (liftU2 f ae be)
  {-# INLINE liftU2 #-}
  liftI2 _ ZO _ = ZO
  liftI2 _ _ ZO = ZO
  liftI2 f (BO ao ae) (BO bo be) = BO (liftI2 f ao bo) (liftI2 f ae be)
  {-# INLINE liftI2 #-}

instance Metric (MV d)
instance Metric (MVE d)
instance Metric (MVO d)

instance (Eq a, Floating a, Epsilon a) => Epsilon (MV d a) where
  nearZero = nearZero . quadrance
  {-# INLINE nearZero #-}
instance (Eq a, Floating a, Epsilon a) => Epsilon (MVE d a) where
  nearZero = nearZero . quadrance
  {-# INLINE nearZero #-}
instance (Eq a, Floating a, Epsilon a) => Epsilon (MVO d a) where
  nearZero = nearZero . quadrance
  {-# INLINE nearZero #-}


instance (Eq a, Floating a) => Num (MV d a) where
  (+) = (^+^)
  (-) = (^-^)
  (*) = mult
  negate = fmap negate
  fromInteger a = MV (fromInteger a) ZO
  abs mv = MV (Sc (norm mv)) ZO
  signum mv = mv ^/ norm mv

instance (Eq a, Floating a) => Num (MVE d a) where
  (+) = (^+^)
  (-) = (^-^)
  (*) = multEE
  negate = fmap negate
  fromInteger = Sc . fromInteger
  abs v = Sc (norm v)
  signum v = v ^/ norm v

instance (Eq a, Floating a) => Num (MVO d a) where
  (+) = (^+^)
  (-) = (^-^)
  negate = fmap negate
  (*) = undefined
  fromInteger = undefined
  abs = undefined
  signum = undefined

-- instance (Eq a, Floating a) => Fractional (MV d a)
-- instance (Eq a, Floating a) => Fractional (MVE d a)

-- The following is wrong, since exp (a + b) != exp a * exp b, and similar for cos and sin, but it does work for blades
-- instance (Eq a, Floating a) => Floating (MV d a) where
--   exp (MV e o) = MV (exp e) ZO * expO o
--   sin (MV e o) = MV (multEE (sin e) (cosO o)) (multEO (cos e) (sinO o))
--   cos (MV e o) = MV (multEE (cos e) (cosO o)) (negate (multEO (sin e) (sinO o)))
-- instance (Eq a, Floating a) => Floating (MVE d a) where
--   exp (Sc a) = mkSc (exp a)
--   exp ZE = Sc 1
--   exp (BE e o) = let expe = exp e in mkBE (multEE expe (cosO o)) (multEO expe (sinO o))
--   sin (Sc a) = mkSc (sin a)
--   sin ZE = ZE
--   sin (BE e o) = let MV cosho sinho = expO o in mkBE (multEE (sin e) cosho) (multEO (cos e) sinho)
--   cos (Sc a) = mkSc (cos a)
--   cos ZE = Sc 1
--   cos (BE e o) = let MV cosho sinho = expO o in BE (multEE (cos e) cosho) (negate (multEO (sin e) sinho))
--   cosh x = (exp x + exp (negate x)) ^/ 2
--   sinh x = (exp x - exp (negate x)) ^/ 2
--
-- expO :: (Eq a, Floating a) => MVO d a -> MV d a
-- expO ZO = 1
-- expO (BO o e) = liftDim (expO o) * MV (mkBE (cosh e) ZO) (mkBO ZO (sinh e))
-- cosO :: (Eq a, Floating a) => MVO d a -> MVE d a
-- cosO ZO = Sc 1
-- cosO (BO o e) = mkBE (multEE (cosO o) (cos e)) (multOE (sinO o) (sin e))
-- sinO :: (Eq a, Floating a) => MVO d a -> MVO d a
-- sinO ZO = ZO
-- sinO (BO o e) = mkBO (multOE (sinO o) (cos e)) (multEE (cosO o) (sin e))


expSer :: (Eq a, Floating a) => MV d a -> MV d a
expSer v = 1 + v * (1 + (v ^/ 2) * (1 + (v ^/ 3) * (1 + (v ^/ 4) * (1 + (v ^/ 5) * (1 + v ^/ 6)))))

cosSer :: (Eq a, Floating a) => MV d a -> MV d a
cosSer v = 1 + v * (0 - (v ^/ 2) * (1 + (v ^/ 3) * (0 - (v ^/ 4) * (1 + (v ^/ 5) * (0 - v ^/ 6)))))

sinSer :: (Eq a, Floating a) => MV d a -> MV d a
sinSer v = 0 + v * (1 - (v ^/ 2) * (0 + (v ^/ 3) * (1 - (v ^/ 4) * (0 + (v ^/ 5)))))

-- instance Floating a => Fractional (MVE 'Z a) where
--   recip ZE = Sc (1/0)
--   recip (Sc a) = Sc (1/a)
--
-- instance (Floating a, Fractional (MVE d a)) => Fractional (MVE ('S d) a) where
--   recip ZE = Sc (1/0)
--   recip (Sc a) = Sc (1/a)
--   recip (BE e o) = BE (e*a) (negate (multOE o b))
--     where
--       a :: Floating a => MVE d a
--       a = recip (multEE e e + multOO o o)
--       b :: Floating a => MVE d a
--       b = multEE (multOO (multOE (recip o) (recip e)) (multOE o e)) a


showFactors :: (Num a, Eq a, Show a) => [(a, String)] -> String
showFactors [] = "0"
showFactors l = fixFirst . concatMap (\(a, es) -> (if abs a /= a then " - " ++ show (-a) else " + " ++ show a) ++ es) . sortOn (length . snd) $ l
  where
    fixFirst (' ':'+':' ':s) = s
    fixFirst (' ':'-':' ':s) = '-':s
    fixFirst s = s

subscript :: Int -> String
subscript 0 = "\x2080"
subscript a = s' a where
  s' 0 = ""
  s' ((`divMod` 10) -> (d, n)) = s' d ++ [toEnum (0x2080 + n)]

instance (Show a, Eq a, Num a, IsDim d) => Show (MV d a) where
  show = pretty

pretty :: (IsDim n, Show a, Eq a, Num a) => MV n a -> String
pretty (MV e o) = showFactors (showAsListE e ++ showAsListO o)

-- deriving instance Show a => Show (MV d a)
deriving instance Show a => Show (MVE d a)
deriving instance Show a => Show (MVO d a)


testE :: MV D3 Double
testE = 1 + 2*e1*e2 + 3*e1*e3 + 4*e2*e3

testO :: MV D3 Double
testO = e1 + 2*e2 + 3*e3 + 4*e1*e2*e3

test :: MV D3 Double
test = (1 + e1) * (1 + 3 * e2) * (1 + 5 * e3)