content: {-# LANGUAGE GADTSyntax #-}
  {-# LANGUAGE ViewPatterns #-}
  {-# LANGUAGE DeriveDataTypeable #-}
  {-# LANGUAGE DeriveGeneric #-}
  {-# OPTIONS_GHC -fno-warn-type-defaults #-}
  {-# LANGUAGE BangPatterns #-}
  {-# LANGUAGE CPP #-}
  {-# LANGUAGE MultiWayIf #-}
  
  
  #ifndef O_LIQUID
  -- Turn off Safe Haskell language extension due to liquid-base re-exports
  {-# LANGUAGE Safe #-}
  #endif
  
  #if __GLASGOW_HASKELL__ == 810
  -- Work around to fix GHC Issue #15304, issue popped up again in GHC 8.10, it should be fixed in GHC 8.12
  -- This code is meant to reproduce MR 2608 for GHC 8.10
  {-# OPTIONS_GHC -funfolding-keeness-factor=1 -funfolding-use-threshold=80 #-}
  #endif
  
  
  
  --------------------------------------------------------------------------------------------
  -- |
  -- Copyright   :  (C) 2017-2022 Nathan Waivio
  -- License     :  BSD3
  -- Maintainer  :  Nathan Waivio <nathan.waivio@gmail.com>
  -- Stability   :  Stable
  -- Portability :  unportable
  --
  -- Library implementing standard functions for the <https://en.wikipedia.org/wiki/Algebra_of_physical_space Algebra of Physical Space> Cl(3,0)
  -- 
  ---------------------------------------------------------------------------------------------
  
  
  module Algebra.Geometric.Cl3
  (-- * The type for the Algebra of Physical Space
   Cl3(..),
   -- * Clifford Conjugate and Complex Conjugate
   bar, dag,
   -- * The littlest singular value
   lsv,
   -- * Constructor Selectors - For optimizing and simplifying calculations
   toR, toV3, toBV, toI,
   toPV, toH, toC,
   toBPV, toODD, toTPV,
   toAPS,
   -- * Pretty Printing for use with Octave
   showOctave,
   -- * Eliminate grades that are less than 'tol' to use a simpler Constructor
   reduce, tol,
  #ifndef O_NO_STORABLE
   -- * Compact Storable types for the Cl3 Constructors with smart constructors
   Cl3_R, toCl3_R, fromCl3_R,
   Cl3_V3, toCl3_V3, fromCl3_V3,
   Cl3_BV, toCl3_BV, fromCl3_BV,
   Cl3_I, toCl3_I, fromCl3_I,
   Cl3_PV, toCl3_PV, fromCl3_PV,
   Cl3_H, toCl3_H, fromCl3_H,
   Cl3_C, toCl3_C, fromCl3_C,
   Cl3_BPV, toCl3_BPV, fromCl3_BPV,
   Cl3_ODD, toCl3_ODD, fromCl3_ODD,
   Cl3_TPV, toCl3_TPV, fromCl3_TPV,
   Cl3_APS, toCl3_APS, fromCl3_APS,
  #endif
  #ifndef O_NO_RANDOM
   -- * Random Instances
   randR, rangeR,
   randV3, rangeV3,
   randBV, rangeBV,
   randI, rangeI,
   randPV, rangePV,
   randH, rangeH,
   randC, rangeC,
   randBPV, rangeBPV,
   randODD, rangeODD,
   randTPV, rangeTPV,
   randAPS, rangeAPS,
   randUnitV3,
   randProjector,
   randNilpotent,
   randUnitary,
  #endif
   -- * Helpful Functions
   eigvals, hasNilpotent,
   spectraldcmp, project,
   mIx, timesI,
   abssignum
  ) where
  
  #ifndef O_NO_DERIVED
  import Data.Data (Typeable, Data)
  import GHC.Generics (Generic)
  import Text.Read (Read,readPrec)
  #endif
  
  import Control.DeepSeq (NFData,rnf)
  
  #ifndef O_NO_STORABLE
  import Foreign.Storable (Storable, sizeOf, alignment, peek, poke)
  import Foreign.Ptr (Ptr, plusPtr, castPtr)
  #endif
  
  #ifndef O_NO_RANDOM
  import System.Random (RandomGen, Random, randomR, random)
  #endif
  
  
  -- | Cl3 provides specialized constructors for sub-algebras and other geometric objects
  -- contained in the algebra.  Cl(3,0), abbreviated to Cl3, is a Geometric Algebra
  -- of 3 dimensional space known as the Algebra of Physical Space (APS).  Geometric Algebras are Real
  -- Clifford Algebras, double precision floats are used to approximate real numbers in this
  -- library.  Single and Double grade combinations are specialized using algebraic datatypes
  -- and live within the APS.
  --
  --   * 'R' is the constructor for the Real Scalar Sub-algebra Grade-0
  --
  --   * 'V3' is the Three Dimensional Real Vector constructor Grade-1
  --
  --   * 'BV' is the Bivector constructor Grade-2 an Imaginary Three Dimensional Vector
  --
  --   * 'I' is the Imaginary constructor Grade-3 and is the Pseudo-Scalar for APS
  --
  --   * 'PV' is the Paravector constructor with Grade-0 and Grade-1 elements, a Real Scalar plus Vector, (R + V3)
  --
  --   * 'H' is the Quaternion constructor it is the Even Sub-algebra with Grade-0 and Grade-2 elements, a Real Scalar plus Bivector, (R + BV)
  --
  --   * 'C' is the Complex constructor it is the Scalar Sub-algebra with Grade-0 and Grade-3 elements, a Real Scalar plus Imaginar Scalar, (R + I)
  --
  --   * 'BPV' is the Biparavector constructor with Grade-1 and Grade-2 elements, a Real Vector plus Bivector, (V3 + BV)
  --
  --   * 'ODD' is the Odd constructor with Grade-1 and Grade-3 elements, a Vector plus Imaginary Scalar, (V3 + I)
  --
  --   * 'TPV' is the Triparavector constructor with Grade-2 and Grade-3 elements, a Bivector plus Imaginary, (BV + I)
  --
  --   * 'APS' is the constructor for an element in the Algebra of Physical Space with Grade-0 through Grade-3 elements
  --
  data Cl3 where
    R   :: !Double -> Cl3 -- Real Scalar Sub-algebra
    V3  :: !Double -> !Double -> !Double -> Cl3 -- Three Dimensional Vectors
    BV  :: !Double -> !Double -> !Double -> Cl3 -- Bivectors, Imaginary Three Dimenstional Vectors
    I   :: !Double -> Cl3 -- Trivector Imaginary Pseudo-Scalar, Imaginary Scalar
    PV  :: !Double -> !Double -> !Double -> !Double -> Cl3 -- Paravector, Real Scalar plus Three Dimensional Real Vector, (R + V3)
    H   :: !Double -> !Double -> !Double -> !Double -> Cl3 -- Quaternion Even Sub-algebra, Real Scalar plus Bivector, (R + BV)
    C   :: !Double -> !Double -> Cl3 -- Complex Sub-algebra, Real Scalar plus Imaginary Scalar, (R + I)
    BPV :: !Double -> !Double -> !Double -> !Double -> !Double -> !Double -> Cl3 -- Biparavector, Vector plus Bivector, (V3 + BV)
    ODD :: !Double -> !Double -> !Double -> !Double -> Cl3 -- Odd, Vector plus Imaginary, (V3 + I)
    TPV :: !Double -> !Double -> !Double -> !Double -> Cl3 -- Triparavector, Bivector plus Imaginary Scalar, (BV + I)
    APS :: !Double -> !Double -> !Double -> !Double -> !Double -> !Double -> !Double -> !Double -> Cl3 -- Algebra of Physical Space
  #ifndef O_NO_DERIVED
      deriving (Show, Read, Typeable, Data, Generic)
  
  #else
  
  -- | In case we don't derive Show, provide 'showOctave' as the Show instance
  instance Show Cl3 where
    show = showOctave
  
  #endif
  
  
  -- | Cl3 can be reduced to a normal form.
  instance NFData Cl3 where
    rnf !_ = ()
  
  
  -- |'showOctave' for useful for debug purposes.
  -- The additional octave definition is needed:  
  -- 
  -- > e0 = [1,0;0,1]; e1=[0,1;1,0]; e2=[0,-i;i,0]; e3=[1,0;0,-1];
  --
  -- This allows one to take advantage of the isomorphism between Cl3 and M(2,C)
  showOctave :: Cl3 -> String
  showOctave (R a0) = show a0 ++ "*e0"
  showOctave (V3 a1 a2 a3) = show a1 ++ "*e1 + " ++ show a2 ++ "*e2 + " ++ show a3 ++ "*e3"
  showOctave (BV a23 a31 a12) = show a23 ++ "i*e1 + " ++ show a31 ++ "i*e2 + " ++ show a12 ++ "i*e3"
  showOctave (I a123) = show a123 ++ "i*e0"
  showOctave (PV a0 a1 a2 a3) = show a0 ++ "*e0 + " ++ show a1 ++ "*e1 + " ++ show a2 ++ "*e2 + " ++ show a3 ++ "*e3"
  showOctave (H a0 a23 a31 a12) = show a0 ++ "*e0 + " ++ show a23 ++ "i*e1 + " ++ show a31 ++ "i*e2 + " ++ show a12 ++ "i*e3"
  showOctave (C a0 a123) = show a0 ++ "*e0 + " ++ show a123 ++ "i*e0"
  showOctave (BPV a1 a2 a3 a23 a31 a12) = show a1 ++ "*e1 + " ++ show a2 ++ "*e2 + " ++ show a3 ++ "*e3 + " ++
                                          show a23 ++ "i*e1 + " ++ show a31 ++ "i*e2 + " ++ show a12 ++ "i*e3"
  showOctave (ODD a1 a2 a3 a123) = show a1 ++ "*e1 + " ++ show a2 ++ "*e2 + " ++ show a3 ++ "*e3 + " ++ show a123 ++ "i*e0"
  showOctave (TPV a23 a31 a12 a123) = show a23 ++ "i*e1 + " ++ show a31 ++ "i*e2 + " ++ show a12 ++ "i*e3 + " ++ show a123 ++ "i*e0"
  showOctave (APS a0 a1 a2 a3 a23 a31 a12 a123) = show a0 ++ "*e0 + " ++ show a1 ++ "*e1 + " ++ show a2 ++ "*e2 + " ++ show a3 ++ "*e3 + " ++
                                                  show a23 ++ "i*e1 + " ++ show a31 ++ "i*e2 + " ++ show a12 ++ "i*e3 + " ++ show a123 ++ "i*e0"
  
  
  -- |Cl(3,0) has the property of equivalence.  "Eq" is "True" when all of the grade elements are equivalent.
  instance Eq Cl3 where
    (R a0) == (R b0) = a0 == b0
  
    (R a0) == (V3 b1 b2 b3) = a0 == 0 && b1 == 0 && b2 == 0 && b3 == 0
    (R a0) == (BV b23 b31 b12) = a0 == 0 && b23 == 0 && b31 == 0 && b12 == 0
    (R a0) == (I b123) = a0 == 0 && b123 == 0
    (R a0) == (PV b0 b1 b2 b3) = a0 == b0 && b1 == 0 && b2 == 0 && b3 == 0
    (R a0) == (H b0 b23 b31 b12) = a0 == b0 && b23 == 0 && b31 == 0 && b12 == 0
    (R a0) == (C b0 b123) = a0 == b0 && b123 == 0
    (R a0) == (BPV b1 b2 b3 b23 b31 b12) = a0 == 0 && b1 == 0 && b2 == 0 && b3 == 0 && b23 == 0 && b31 == 0 && b12 == 0
    (R a0) == (ODD b1 b2 b3 b123) = a0 == 0 && b1 == 0 && b2 == 0 && b3 == 0 && b123 == 0
    (R a0) == (TPV b23 b31 b12 b123) = a0 == 0 && b23 == 0 && b31 == 0 && b12 == 0 && b123 == 0
    (R a0) == (APS b0 b1 b2 b3 b23 b31 b12 b123) = a0 == b0 && b1 == 0 && b2 == 0 && b3 == 0 && b23 == 0 && b31 == 0 && b12 == 0 && b123 == 0
  
    (V3 a1 a2 a3) == (R b0) = a1 == 0 && a2 == 0 && a3 == 0 && b0 == 0
    (BV a23 a31 a12) == (R b0) = a23 == 0 && a31 == 0 && a12 == 0 && b0 == 0
    (I a123) == (R b0) = a123 == 0 && b0 == 0
    (PV a0 a1 a2 a3) == (R b0) = a0 == b0 && a1 == 0 && a2 == 0 && a3 == 0
    (H a0 a23 a31 a12) == (R b0) = a0 == b0 && a23 == 0 && a31 == 0 && a12 == 0
    (C a0 a123) == (R b0) = a0 == b0 && a123 == 0
    (BPV a1 a2 a3 a23 a31 a12) == (R b0) = a1 == 0 && a2 == 0 && a3 == 0 && a23 == 0 && a31 == 0 && a12 == 0 && b0 == 0
    (ODD a1 a2 a3 a123) == (R b0) = a1 == 0 && a2 == 0 && a3 == 0 && a123 == 0 && b0 == 0
    (TPV a23 a31 a12 a123) == (R b0) = a23 == 0 && a31 == 0 && a12 == 0 && a123 == 0 && b0 == 0
    (APS a0 a1 a2 a3 a23 a31 a12 a123) == (R b0) = a0 == b0 && a1 == 0 && a2 == 0 && a3 == 0 && a23 == 0 && a31 == 0 && a12 == 0 && a123 == 0
  
    (V3 a1 a2 a3) == (V3 b1 b2 b3) = a1 == b1 && a2 == b2 && a3 == b3
  
    (V3 a1 a2 a3) == (BV b23 b31 b12) = a1 == 0 && a2 == 0 && a3 == 0 && b23 == 0 && b31 == 0 && b12 == 0
    (V3 a1 a2 a3) == (I b123) = a1 == 0 && a2 == 0 && a3 == 0 && b123 == 0
    (V3 a1 a2 a3) == (PV b0 b1 b2 b3) = a1 == b1 && a2 == b2 && a3 == b3 && b0 == 0
    (V3 a1 a2 a3) == (H b0 b23 b31 b12) = a1 == 0 && a2 == 0 && a3 == 0 && b0 == 0 && b23 == 0 && b31 == 0 && b12 == 0
    (V3 a1 a2 a3) == (C b0 b123) = a1 == 0 && a2 == 0 && a3 == 0 && b0 == 0 && b123 == 0
    (V3 a1 a2 a3) == (BPV b1 b2 b3 b23 b31 b12) = a1 == b1 && a2 == b2 && a3 == b3 && b23 == 0 && b31 == 0 && b12 == 0
    (V3 a1 a2 a3) == (ODD b1 b2 b3 b123) = a1 == b1 && a2 == b2 && a3 == b3 && b123 == 0
    (V3 a1 a2 a3) == (TPV b23 b31 b12 b123) = a1 == 0 && a2 == 0 && a3 == 0 && b23 == 0 && b31 == 0 && b12 == 0 && b123 == 0
    (V3 a1 a2 a3) == (APS b0 b1 b2 b3 b23 b31 b12 b123) = a1 == b1 && a2 == b2 && a3 == b3 && b0 == 0 && b23 == 0 && b31 == 0 && b12 == 0 && b123 == 0
  
    (BV a23 a31 a12) == (V3 b1 b2 b3) = a23 == 0 && a31 == 0 && a12 == 0 && b1 == 0 && b2 == 0 && b3 == 0
    (I a123) == (V3 b1 b2 b3) = a123 == 0 && b1 == 0 && b2 == 0 && b3 == 0
    (PV a0 a1 a2 a3) == (V3 b1 b2 b3) = a0 == 0 && a1 == b1 && a2 == b2 && a3 == b3
    (H a0 a23 a31 a12) == (V3 b1 b2 b3) = a0 == 0 && a23 == 0 && a31 == 0 && a12 == 0 && b1 == 0 && b2 == 0 && b3 == 0
    (C a0 a123) == (V3 b1 b2 b3) = a0 == 0 && a123 == 0 && b1 == 0 && b2 == 0 && b3 == 0
    (BPV a1 a2 a3 a23 a31 a12) == (V3 b1 b2 b3) = a1 == b1 && a2 == b2 && a3 == b3 && a23 == 0 && a31 == 0 && a12 == 0
    (ODD a1 a2 a3 a123) == (V3 b1 b2 b3) = a1 == b1 && a2 == b2 && a3 == b3 && a123 == 0
    (TPV a23 a31 a12 a123) == (V3 b1 b2 b3) = b1 == 0 && b2 == 0 && b3 == 0 && a23 == 0 && a31 == 0 && a12 == 0 && a123 == 0
    (APS a0 a1 a2 a3 a23 a31 a12 a123) == (V3 b1 b2 b3) = a0 == 0 && a1 == b1 && a2 == b2 && a3 == b3 && a23 == 0 && a31 == 0 && a12 == 0 && a123 == 0
  
    (BV a23 a31 a12) == (BV b23 b31 b12) = a23 == b23 && a31 == b31 && a12 == b12
  
    (BV a23 a31 a12) == (I b123) = a23 == 0 && a31 == 0 && a12 == 0 && b123 == 0
    (BV a23 a31 a12) == (PV b0 b1 b2 b3) = a23 == 0 && a31 == 0 && a12 == 0 && b0 == 0 && b1 == 0 && b2 == 0 && b3 == 0
    (BV a23 a31 a12) == (H b0 b23 b31 b12) = a23 == b23 && a31 == b31 && a12 == b12 && b0 == 0
    (BV a23 a31 a12) == (C b0 b123) = a23 == 0 && a31 == 0 && a12 == 0 && b0 == 0 && b123 == 0
    (BV a23 a31 a12) == (BPV b1 b2 b3 b23 b31 b12) = a23 == b23 && a31 == b31 && a12 == b12 && b1 == 0 && b2 == 0 && b3 == 0
    (BV a23 a31 a12) == (ODD b1 b2 b3 b123) = a23 == 0 && a31 == 0 && a12 == 0 && b1 == 0 && b2 == 0 && b3 == 0 && b123 == 0
    (BV a23 a31 a12) == (TPV b23 b31 b12 b123) = a23 == b23 && a31 == b31 && a12 == b12 && b123 == 0
    (BV a23 a31 a12) == (APS b0 b1 b2 b3 b23 b31 b12 b123) = a23 == b23 && a31 == b31 && a12 == b12 && b0 == 0 && b1 == 0 && b2 == 0 && b3 == 0 && b123 == 0
  
    (I a123) == (BV b23 b31 b12) = a123 == 0 && b23 == 0 && b31 == 0 && b12 == 0
    (PV a0 a1 a2 a3) == (BV b23 b31 b12) = a0 == 0 && a1 == 0 && a2 == 0 && a3 == 0 && b23 == 0 && b31 == 0 && b12 == 0
    (H a0 a23 a31 a12) == (BV b23 b31 b12) = a0 == 0 && a23 == b23 && a31 == b31 && a12 == b12
    (C a0 a123) == (BV b23 b31 b12) = a0 == 0 && a123 == 0 && b23 == 0 && b31 == 0 && b12 == 0
    (BPV a1 a2 a3 a23 a31 a12) == (BV b23 b31 b12) = a1 == 0 && a2 == 0 && a3 == 0 && a23 == b23 && a31 == b31 && a12 == b12
    (ODD a1 a2 a3 a123) == (BV b23 b31 b12) = a1 == 0 && a2 == 0 && a3 == 0 && a123 == 0 && b23 == 0 && b31 == 0 && b12 == 0
    (TPV a23 a31 a12 a123) == (BV b23 b31 b12) = a23 == b23 && a31 == b31 && a12 == b12 && a123 == 0
    (APS a0 a1 a2 a3 a23 a31 a12 a123) == (BV b23 b31 b12) = a0 == 0 && a1 == 0 && a2 == 0 && a3 == 0 && a23 == b23 && a31 == b31 && a12 == b12 && a123 == 0
  
    (I a123) == (I b123) = a123 == b123
  
    (I a123) == (PV b0 b1 b2 b3) = a123 == 0 && b0 == 0 && b1 == 0 && b2 == 0 && b3 == 0
    (I a123) == (H b0 b23 b31 b12) = a123 == 0 && b0 == 0 && b23 == 0 && b31 == 0 && b12 == 0
    (I a123) == (C b0 b123) = a123 == b123 && b0 == 0
    (I a123) == (BPV b1 b2 b3 b23 b31 b12) = a123 == 0 && b1 == 0 && b2 == 0 && b3 == 0 && b23 == 0 && b31 == 0 && b12 == 0
    (I a123) == (ODD b1 b2 b3 b123) = a123 == b123 && b1 == 0 && b2 == 0 && b3 == 0
    (I a123) == (TPV b23 b31 b12 b123) = a123 == b123 && b23 == 0 && b31 == 0 && b12 == 0
    (I a123) == (APS b0 b1 b2 b3 b23 b31 b12 b123) = a123 == b123 && b0 == 0 && b1 == 0 && b2 == 0 && b3 == 0 && b23 == 0 && b31 == 0 && b12 == 0
  
    (PV a0 a1 a2 a3) == (I b123) = b123 == 0 && a0 == 0 && a1 == 0 && a2 == 0 && a3 == 0
    (H a0 a23 a31 a12) == (I b123) = b123 == 0 && a0 == 0 && a23 == 0 && a31 == 0 && a12 == 0
    (C a0 a123) == (I b123) = a123 == b123 && a0 == 0
    (BPV a1 a2 a3 a23 a31 a12) == (I b123) = b123 == 0 && a1 == 0 && a2 == 0 && a3 == 0 && a23 == 0 && a31 == 0 && a12 == 0
    (ODD a1 a2 a3 a123) == (I b123) = a123 == b123 && a1 == 0 && a2 == 0 && a3 == 0
    (TPV a23 a31 a12 a123) == (I b123) = a123 == b123 && a23 == 0 && a31 == 0 && a12 == 0
    (APS a0 a1 a2 a3 a23 a31 a12 a123) == (I b123) = a123 == b123 && a0 == 0 && a1 == 0 && a2 == 0 && a3 == 0 && a23 == 0 && a31 == 0 && a12 == 0
  
    (PV a0 a1 a2 a3) == (PV b0 b1 b2 b3) = a0 == b0 && a1 == b1 && a2 == b2 && a3 == b3
  
    (PV a0 a1 a2 a3) == (H b0 b23 b31 b12) = a0 == b0 && a1 == 0 && a2 == 0 && a3 == 0 && b23 == 0 && b31 == 0 && b12 == 0
    (PV a0 a1 a2 a3) == (C b0 b123) = a0 == b0 && a1 == 0 && a2 == 0 && a3 == 0 && b123 == 0
    (PV a0 a1 a2 a3) == (BPV b1 b2 b3 b23 b31 b12) = a0 == 0 && a1 == b1 && a2 == b2 && a3 == b3 && b23 == 0 && b31 == 0 && b12 == 0
    (PV a0 a1 a2 a3) == (ODD b1 b2 b3 b123) = a0 == 0 && a1 == b1 && a2 == b2 && a3 == b3 && b123 == 0
    (PV a0 a1 a2 a3) == (TPV b23 b31 b12 b123) = a0 == 0 && a1 == 0 && a2 == 0 && a3 == 0 && b23 == 0 && b31 == 0 && b12 == 0 && b123 == 0
    (PV a0 a1 a2 a3) == (APS b0 b1 b2 b3 b23 b31 b12 b123) = a0 == b0 && a1 == b1 && a2 == b2 && a3 == b3 && b23 == 0 && b31 == 0 && b12 == 0 && b123 == 0
  
    (H a0 a23 a31 a12) == (PV b0 b1 b2 b3) = a0 == b0 && a23 == 0 && a31 == 0 && a12 == 0 && b1 == 0 && b2 == 0 && b3 == 0
    (C a0 a123) == (PV b0 b1 b2 b3) = a0 == b0 && a123 == 0 && b1 == 0 && b2 == 0 && b3 == 0
    (BPV a1 a2 a3 a23 a31 a12) == (PV b0 b1 b2 b3) = a1 == b1 && a2 == b2 && a3 == b3 && a23 == 0 && a31 == 0 && a12 == 0 && b0 == 0
    (ODD a1 a2 a3 a123) == (PV b0 b1 b2 b3) = a1 == b1 && a2 == b2 && a3 == b3 && a123 == 0 && b0 == 0
    (TPV a23 a31 a12 a123) == (PV b0 b1 b2 b3) = a23 == 0 && a31 == 0 && a12 == 0 && b0 == 0 && a123 == 0 && b1 == 0 && b2 == 0 && b3 == 0
    (APS a0 a1 a2 a3 a23 a31 a12 a123) == (PV b0 b1 b2 b3) = a0 == b0 && a1 == b1 && a2 == b2 && a3 == b3 && a23 == 0 && a31 == 0 && a12 == 0 && a123 == 0
  
    (H a0 a23 a31 a12) == (H b0 b23 b31 b12) = a0 == b0 && a23 == b23 && a31 == b31 && a12 == b12
  
    (H a0 a23 a31 a12) == (C b0 b123) = a0 == b0 && a23 == 0 && a31 == 0 && a12 == 0 && b123 == 0
    (H a0 a23 a31 a12) == (BPV b1 b2 b3 b23 b31 b12) = a0 == 0 && a23 == b23 && a31 == b31 && a12 == b12 && b1 == 0 && b2 == 0 && b3 == 0
    (H a0 a23 a31 a12) == (ODD b1 b2 b3 b123) = a0 == 0 && b1 == 0 && b2 == 0 && b3 == 0 && a23 == 0 && a31 == 0 && a12 == 0 && b123 == 0
    (H a0 a23 a31 a12) == (TPV b23 b31 b12 b123) = a0 == 0 && a23 == b23 && a31 == b31 && a12 == b12 && b123 == 0
    (H a0 a23 a31 a12) == (APS b0 b1 b2 b3 b23 b31 b12 b123) = a0 == b0 && a23 == b23 && a31 == b31 && a12 == b12 && b1 == 0 && b2 == 0 && b3 == 0 && b123 == 0
  
    (C a0 a123) == (H b0 b23 b31 b12) = a0 == b0 && a123 == 0 && b23 == 0 && b31 == 0 && b12 == 0
    (BPV a1 a2 a3 a23 a31 a12) == (H b0 b23 b31 b12) = a1 == 0 && a2 == 0 && a3 == 0 && a23 == b23 && a31 == b31 && a12 == b12 && b0 == 0
    (ODD a1 a2 a3 a123) == (H b0 b23 b31 b12) = a1 == 0 && a2 == 0 && a3 == 0 && a123 == 0 && b23 == 0 && b31 == 0 && b12 == 0 && b0 == 0
    (TPV a23 a31 a12 a123) == (H b0 b23 b31 b12) = a23 == b23 && a31 == b31 && a12 == b12 && b0 == 0 && a123 == 0
    (APS a0 a1 a2 a3 a23 a31 a12 a123) == (H b0 b23 b31 b12) = a0 == b0 && a1 == 0 && a2 == 0 && a3 == 0 && a23 == b23 && a31 == b31 && a12 == b12 && a123 == 0
  
    (C a0 a123) == (C b0 b123) = a0 == b0 && a123 == b123
  
    (C a0 a123) == (BPV b1 b2 b3 b23 b31 b12) = a0 == 0 && a123 == 0 && b1 == 0 && b2 == 0 && b3 == 0 && b23 == 0 && b31 == 0 && b12 == 0
    (C a0 a123) == (ODD b1 b2 b3 b123) = a0 == 0 && a123 == b123 && b1 == 0 && b2 == 0 && b3 == 0
    (C a0 a123) == (TPV b23 b31 b12 b123) = a0 == 0 && a123 == b123 && b23 == 0 && b31 == 0 && b12 == 0
    (C a0 a123) == (APS b0 b1 b2 b3 b23 b31 b12 b123) = a0 == b0 && a123 == b123 && b1 == 0 && b2 == 0 && b3 == 0 && b23 == 0 && b31 == 0 && b12 == 0
  
    (BPV a1 a2 a3 a23 a31 a12) == (C b0 b123) = a1 == 0 && a2 == 0 && a3 == 0 && a23 == 0 && a31 == 0 && a12 == 0 && b0 == 0 && b123 == 0
    (ODD a1 a2 a3 a123) == (C b0 b123) = b0 == 0 && a123 == b123 && a1 == 0 && a2 == 0 && a3 == 0
    (TPV a23 a31 a12 a123) == (C b0 b123) = b0 == 0 && a123 == b123 && a23 == 0 && a31 == 0 && a12 == 0
    (APS a0 a1 a2 a3 a23 a31 a12 a123) == (C b0 b123) = a0 == b0 && a123 == b123 && a1 == 0 && a2 == 0 && a3 == 0 && a23 == 0 && a31 == 0 && a12 == 0
  
    (BPV a1 a2 a3 a23 a31 a12) == (BPV b1 b2 b3 b23 b31 b12) = a1 == b1 && a2 == b2 && a3 == b3 && a23 == b23 && a31 == b31 && a12 == b12
  
    (BPV a1 a2 a3 a23 a31 a12) == (ODD b1 b2 b3 b123) = a1 == b1 && a2 == b2 && a3 == b3 && b123 == 0 && a23 == 0 && a31 == 0 && a12 == 0
    (BPV a1 a2 a3 a23 a31 a12) == (TPV b23 b31 b12 b123) = a23 == b23 && a31 == b31 && a12 == b12 && b123 == 0 && a1 == 0 && a2 == 0 && a3 == 0
    (BPV a1 a2 a3 a23 a31 a12) == (APS b0 b1 b2 b3 b23 b31 b12 b123) = a1 == b1 && a2 == b2 && a3 == b3 && a23 == b23 && a31 == b31 && a12 == b12
                                                                                && b0 == 0 && b123 == 0
  
    (ODD a1 a2 a3 a123) == (BPV b1 b2 b3 b23 b31 b12) = a1 == b1 && a2 == b2 && a3 == b3 && a123 == 0 && b23 == 0 && b31 == 0 && b12 == 0
    (TPV a23 a31 a12 a123) == (BPV b1 b2 b3 b23 b31 b12) = a23 == b23 && a31 == b31 && a12 == b12 && a123 == 0 && b1 == 0 && b2 == 0 && b3 == 0
    (APS a0 a1 a2 a3 a23 a31 a12 a123) == (BPV b1 b2 b3 b23 b31 b12) = a0 == 0 && a1 == b1 && a2 == b2 && a3 == b3 && a23 == b23 && a31 == b31
                                                                               && a12 == b12 && a123 == 0
  
    (ODD a1 a2 a3 a123) == (ODD b1 b2 b3 b123) = a1 == b1 && a2 == b2 && a3 == b3 && a123 == b123
  
    (ODD a1 a2 a3 a123) == (TPV b23 b31 b12 b123) = a123 == b123 && a1 == 0 && a2 == 0 && a3 == 0 && b23 == 0 && b31 == 0 && b12 == 0
    (ODD a1 a2 a3 a123) == (APS b0 b1 b2 b3 b23 b31 b12 b123) = a1 == b1 && a2 == b2 && a3 == b3 && a123 == b123 && b0 == 0 && b23 == 0 && b31 == 0 && b12 == 0
  
    (TPV a23 a31 a12 a123) == (ODD b1 b2 b3 b123) = a123 == b123 && b1 == 0 && b2 == 0 && b3 == 0 && a23 == 0 && a31 == 0 && a12 == 0
    (APS a0 a1 a2 a3 a23 a31 a12 a123) == (ODD b1 b2 b3 b123) = a1 == b1 && a2 == b2 && a3 == b3 && a123 == b123 && a0 == 0 && a23 == 0 && a31 == 0 && a12 == 0
  
    (TPV a23 a31 a12 a123) == (TPV b23 b31 b12 b123) = a23 == b23 && a31 == b31 && a12 == b12 && a123 == b123
  
    (TPV a23 a31 a12 a123) == (APS b0 b1 b2 b3 b23 b31 b12 b123) = a23 == b23 && a31 == b31 && a12 == b12 && a123 == b123
                                                                              && b0 == 0 && b1 == 0 && b2 == 0 && b3 == 0
  
    (APS a0 a1 a2 a3 a23 a31 a12 a123) == (TPV b23 b31 b12 b123) = a23 == b23 && a31 == b31 && a12 == b12 && a123 == b123
                                                                              && a0 == 0 && a1 == 0 && a2 == 0 && a3 == 0
  
    (APS a0 a1 a2 a3 a23 a31 a12 a123) == (APS b0 b1 b2 b3 b23 b31 b12 b123) = a0 == b0 && a1 == b1 && a2 == b2 && a3 == b3 && a23 == b23
                                                                                        && a31 == b31 && a12 == b12 && a123 == b123
  
  
  -- |Cl3 has a total preorder ordering in which all pairs are comparable by two real valued functions.
  -- Comparison of two reals is just the typical real compare function.  Comparison of to imaginary numbers
  -- is just the typical comparison function.  When reals are compared to anything else it will compare the
  -- absolute value of the reals to the magnitude of the other cliffor.  Compare of two complex values
  -- compares the polar magnitude of the complex numbers.  Compare of two vectors compares the vector
  -- magnitudes.  The Ord instance for the general case is based on the singular values of each cliffor and
  -- this Ordering compares the largest singular value 'abs' and then the littlest singular value 'lsv'.
  -- Some arbitrary cliffors may return EQ for Ord but not be exactly '==' equivalent, but they are related
  -- by a right and left multiplication of two unitary elements.  For instance for the Cliffors A and B,
  -- A == B could be False, but compare A B is EQ, because A * V = U * B, where V and U are unitary.  
  instance Ord Cl3 where
    compare (R a0) (R b0) = compare a0 b0 -- Real Numbers have a total order within the limitations of Double Precision comparison
    compare (I a123) (I b123) = compare a123 b123 -- Imaginary Numbers have a total order within the limitations of Double Precision comparison
    compare cliffor1 cliffor2 =
       let (R a0) = abs cliffor1
           (R b0) = abs cliffor2
           (R a0') = lsv cliffor1
           (R b0') = lsv cliffor2
       in case compare a0 b0 of
            LT -> LT
            GT -> GT
            EQ -> compare a0' b0'
  
  
  
  -- |Cl3 has a "Num" instance.  "Num" is addition, geometric product, negation, 'abs' the largest
  -- singular value, and 'signum'.
  -- 
  instance Num Cl3 where
    -- | Cl3 can be added
    (R a0) + (R b0) = R (a0 + b0)
  
    (R a0) + (V3 b1 b2 b3) = PV a0 b1 b2 b3
    (R a0) + (BV b23 b31 b12) = H a0 b23 b31 b12
    (R a0) + (I b123) = C a0 b123
    (R a0) + (PV b0 b1 b2 b3) = PV (a0 + b0) b1 b2 b3
    (R a0) + (H b0 b23 b31 b12) = H (a0 + b0) b23 b31 b12
    (R a0) + (C b0 b123) = C (a0 + b0) b123
    (R a0) + (BPV b1 b2 b3 b23 b31 b12) = APS a0 b1 b2 b3 b23 b31 b12 0
    (R a0) + (ODD b1 b2 b3 b123) = APS a0 b1 b2 b3 0 0 0 b123
    (R a0) + (TPV b23 b31 b12 b123) = APS a0 0 0 0 b23 b31 b12 b123
    (R a0) + (APS b0 b1 b2 b3 b23 b31 b12 b123) = APS (a0 + b0) b1 b2 b3 b23 b31 b12 b123
  
    (V3 a1 a2 a3) + (R b0) = PV b0 a1 a2 a3
    (BV a23 a31 a12) + (R b0) = H b0 a23 a31 a12
    (I a123) + (R b0) = C b0 a123
    (PV a0 a1 a2 a3) + (R b0) = PV (a0 + b0) a1 a2 a3
    (H a0 a23 a31 a12) + (R b0) = H (a0 + b0) a23 a31 a12
    (C a0 a123) + (R b0) = C (a0 + b0) a123
    (BPV a1 a2 a3 a23 a31 a12) + (R b0) = APS b0 a1 a2 a3 a23 a31 a12 0
    (ODD a1 a2 a3 a123) + (R b0) = APS b0 a1 a2 a3 0 0 0 a123
    (TPV a23 a31 a12 a123) + (R b0) = APS b0 0 0 0 a23 a31 a12 a123
    (APS a0 a1 a2 a3 a23 a31 a12 a123) + (R b0) = APS (a0 + b0) a1 a2 a3 a23 a31 a12 a123
  
    (V3 a1 a2 a3) + (V3 b1 b2 b3) = V3 (a1 + b1) (a2 + b2) (a3 + b3)
  
    (V3 a1 a2 a3) + (BV b23 b31 b12) = BPV a1 a2 a3 b23 b31 b12
    (V3 a1 a2 a3) + (I b123) = ODD a1 a2 a3 b123
    (V3 a1 a2 a3) + (PV b0 b1 b2 b3) = PV b0 (a1 + b1) (a2 + b2) (a3 + b3)
    (V3 a1 a2 a3) + (H b0 b23 b31 b12) = APS b0 a1 a2 a3 b23 b31 b12 0
    (V3 a1 a2 a3) + (C b0 b123) = APS b0 a1 a2 a3 0 0 0 b123
    (V3 a1 a2 a3) + (BPV b1 b2 b3 b23 b31 b12) = BPV (a1 + b1) (a2 + b2) (a3 + b3) b23 b31 b12
    (V3 a1 a2 a3) + (ODD b1 b2 b3 b123) = ODD (a1 + b1) (a2 + b2) (a3 + b3) b123
    (V3 a1 a2 a3) + (TPV b23 b31 b12 b123) = APS 0 a1 a2 a3 b23 b31 b12 b123
    (V3 a1 a2 a3) + (APS b0 b1 b2 b3 b23 b31 b12 b123) = APS b0 (a1 + b1) (a2 + b2) (a3 + b3) b23 b31 b12 b123
  
    (BV a23 a31 a12) + (V3 b1 b2 b3) = BPV b1 b2 b3 a23 a31 a12
    (I a123) + (V3 b1 b2 b3) = ODD b1 b2 b3 a123
    (PV a0 a1 a2 a3) + (V3 b1 b2 b3) = PV a0 (a1 + b1) (a2 + b2) (a3 + b3)
    (H a0 a23 a31 a12) + (V3 b1 b2 b3) = APS a0 b1 b2 b3 a23 a31 a12 0
    (C a0 a123) + (V3 b1 b2 b3) = APS a0 b1 b2 b3 0 0 0 a123
    (BPV a1 a2 a3 a23 a31 a12) + (V3 b1 b2 b3) = BPV (a1 + b1) (a2 + b2) (a3 + b3) a23 a31 a12
    (ODD a1 a2 a3 a123) + (V3 b1 b2 b3) = ODD (a1 + b1) (a2 + b2) (a3 + b3) a123
    (TPV a23 a31 a12 a123) + (V3 b1 b2 b3) = APS 0 b1 b2 b3 a23 a31 a12 a123
    (APS a0 a1 a2 a3 a23 a31 a12 a123) + (V3 b1 b2 b3) = APS a0 (a1 + b1) (a2 + b2) (a3 + b3) a23 a31 a12 a123
  
    (BV a23 a31 a12) + (BV b23 b31 b12) = BV (a23 + b23) (a31 + b31) (a12 + b12)
  
    (BV a23 a31 a12) + (I b123) = TPV a23 a31 a12 b123
    (BV a23 a31 a12) + (PV b0 b1 b2 b3) = APS b0 b1 b2 b3 a23 a31 a12 0
    (BV a23 a31 a12) + (H b0 b23 b31 b12) = H b0 (a23 + b23) (a31 + b31) (a12 + b12)
    (BV a23 a31 a12) + (C b0 b123) = APS b0 0 0 0 a23 a31 a12 b123
    (BV a23 a31 a12) + (BPV b1 b2 b3 b23 b31 b12) = BPV b1 b2 b3 (a23 + b23) (a31 + b31) (a12 + b12)
    (BV a23 a31 a12) + (ODD b1 b2 b3 b123) = APS 0 b1 b2 b3 a23 a31 a12 b123
    (BV a23 a31 a12) + (TPV b23 b31 b12 b123) = TPV (a23 + b23) (a31 + b31) (a12 + b12) b123
    (BV a23 a31 a12) + (APS b0 b1 b2 b3 b23 b31 b12 b123) = APS b0 b1 b2 b3 (a23 + b23) (a31 + b31) (a12 + b12) b123
  
    (I a123) + (BV b23 b31 b12) = TPV b23 b31 b12 a123
    (PV a0 a1 a2 a3) + (BV b23 b31 b12) = APS a0 a1 a2 a3 b23 b31 b12 0
    (H a0 a23 a31 a12) + (BV b23 b31 b12) = H a0 (a23 + b23) (a31 + b31) (a12 + b12)
    (C a0 a123) + (BV b23 b31 b12) = APS a0 0 0 0 b23 b31 b12 a123
    (BPV a1 a2 a3 a23 a31 a12) + (BV b23 b31 b12) = BPV a1 a2 a3 (a23 + b23) (a31 + b31) (a12 + b12)
    (ODD a1 a2 a3 a123) + (BV b23 b31 b12) = APS 0 a1 a2 a3 b23 b31 b12 a123
    (TPV a23 a31 a12 a123) + (BV b23 b31 b12) = TPV (a23 + b23) (a31 + b31) (a12 + b12) a123
    (APS a0 a1 a2 a3 a23 a31 a12 a123) + (BV b23 b31 b12) = APS a0 a1 a2 a3 (a23 + b23) (a31 + b31) (a12 + b12) a123
  
    (I a123) + (I b123) = I (a123 + b123)
  
    (I a123) + (PV b0 b1 b2 b3) = APS b0 b1 b2 b3 0 0 0 a123
    (I a123) + (H b0 b23 b31 b12) = APS b0 0 0 0 b23 b31 b12 a123
    (I a123) + (C b0 b123) = C b0 (a123 + b123)
    (I a123) + (BPV b1 b2 b3 b23 b31 b12) = APS 0 b1 b2 b3 b23 b31 b12 a123
    (I a123) + (ODD b1 b2 b3 b123) = ODD b1 b2 b3 (a123 + b123)
    (I a123) + (TPV b23 b31 b12 b123) = TPV b23 b31 b12 (a123 + b123)
    (I a123) + (APS b0 b1 b2 b3 b23 b31 b12 b123) = APS b0 b1 b2 b3 b23 b31 b12 (a123 + b123)
  
    (PV a0 a1 a2 a3) + (I b123) = APS a0 a1 a2 a3 0 0 0 b123
    (H a0 a23 a31 a12) + (I b123) = APS a0 0 0 0 a23 a31 a12 b123
    (C a0 a123) + (I b123) = C a0 (a123 + b123)
    (BPV a1 a2 a3 a23 a31 a12) + (I b123) = APS 0 a1 a2 a3 a23 a31 a12 b123
    (ODD a1 a2 a3 a123) + (I b123) = ODD a1 a2 a3 (a123 + b123)
    (TPV a23 a31 a12 a123) + (I b123) = TPV a23 a31 a12 (a123 + b123)
    (APS a0 a1 a2 a3 a23 a31 a12 a123) + (I b123) = APS a0 a1 a2 a3 a23 a31 a12 (a123 + b123)
  
    (PV a0 a1 a2 a3) + (PV b0 b1 b2 b3) = PV (a0 + b0) (a1 + b1) (a2 + b2) (a3 + b3)
  
    (PV a0 a1 a2 a3) + (H b0 b23 b31 b12) = APS (a0 + b0) a1 a2 a3 b23 b31 b12 0
    (PV a0 a1 a2 a3) + (C b0 b123) = APS (a0 + b0) a1 a2 a3 0 0 0 b123
    (PV a0 a1 a2 a3) + (BPV b1 b2 b3 b23 b31 b12) = APS a0 (a1 + b1) (a2 + b2) (a3 + b3) b23 b31 b12 0
    (PV a0 a1 a2 a3) + (ODD b1 b2 b3 b123) = APS a0 (a1 + b1) (a2 + b2) (a3 + b3) 0 0 0 b123
    (PV a0 a1 a2 a3) + (TPV b23 b31 b12 b123) = APS a0 a1 a2 a3 b23 b31 b12 b123
    (PV a0 a1 a2 a3) + (APS b0 b1 b2 b3 b23 b31 b12 b123) = APS (a0 + b0) (a1 + b1) (a2 + b2) (a3 + b3) b23 b31 b12 b123
  
    (H a0 a23 a31 a12) + (PV b0 b1 b2 b3) = APS (a0 + b0) b1 b2 b3 a23 a31 a12 0
    (C a0 a123) + (PV b0 b1 b2 b3) = APS (a0 + b0) b1 b2 b3 0 0 0 a123
    (BPV a1 a2 a3 a23 a31 a12) + (PV b0 b1 b2 b3) = APS b0 (a1 + b1) (a2 + b2) (a3 + b3) a23 a31 a12 0
    (ODD a1 a2 a3 a123) + (PV b0 b1 b2 b3) = APS b0 (a1 + b1) (a2 + b2) (a3 + b3) 0 0 0 a123
    (TPV a23 a31 a12 a123) + (PV b0 b1 b2 b3) = APS b0 b1 b2 b3 a23 a31 a12 a123
    (APS a0 a1 a2 a3 a23 a31 a12 a123) + (PV b0 b1 b2 b3) = APS (a0 + b0) (a1 + b1) (a2 + b2) (a3 + b3) a23 a31 a12 a123
  
    (H a0 a23 a31 a12) + (H b0 b23 b31 b12) = H (a0 + b0) (a23 + b23) (a31 + b31) (a12 + b12)
  
    (H a0 a23 a31 a12) + (C b0 b123) = APS (a0 + b0) 0 0 0 a23 a31 a12 b123
    (H a0 a23 a31 a12) + (BPV b1 b2 b3 b23 b31 b12) = APS a0 b1 b2 b3 (a23 + b23) (a31 + b31) (a12 + b12) 0
    (H a0 a23 a31 a12) + (ODD b1 b2 b3 b123) = APS a0 b1 b2 b3 a23 a31 a12 b123
    (H a0 a23 a31 a12) + (TPV b23 b31 b12 b123) = APS a0 0 0 0 (a23 + b23) (a31 + b31) (a12 + b12) b123
    (H a0 a23 a31 a12) + (APS b0 b1 b2 b3 b23 b31 b12 b123) = APS (a0 + b0) b1 b2 b3 (a23 + b23) (a31 + b31) (a12 + b12) b123
  
    (C a0 a123) + (H b0 b23 b31 b12) = APS (a0 + b0) 0 0 0 b23 b31 b12 a123
    (BPV a1 a2 a3 a23 a31 a12) + (H b0 b23 b31 b12) = APS b0 a1 a2 a3 (a23 + b23) (a31 + b31) (a12 + b12) 0
    (ODD a1 a2 a3 a123) + (H b0 b23 b31 b12) = APS b0 a1 a2 a3 b23 b31 b12 a123
    (TPV a23 a31 a12 a123) + (H b0 b23 b31 b12) = APS b0 0 0 0 (a23 + b23) (a31 + b31) (a12 + b12) a123
    (APS a0 a1 a2 a3 a23 a31 a12 a123) + (H b0 b23 b31 b12) = APS (a0 + b0) a1 a2 a3 (a23 + b23) (a31 + b31) (a12 + b12) a123
  
    (C a0 a123) + (C b0 b123) = C (a0 + b0) (a123 + b123)
  
    (C a0 a123) + (BPV b1 b2 b3 b23 b31 b12) = APS a0 b1 b2 b3 b23 b31 b12 a123
    (C a0 a123) + (ODD b1 b2 b3 b123) = APS a0 b1 b2 b3 0 0 0 (a123 + b123)
    (C a0 a123) + (TPV b23 b31 b12 b123) = APS a0 0 0 0 b23 b31 b12 (a123 + b123)
    (C a0 a123) + (APS b0 b1 b2 b3 b23 b31 b12 b123) = APS (a0 + b0) b1 b2 b3 b23 b31 b12 (a123 + b123)
  
    (BPV a1 a2 a3 a23 a31 a12) + (C b0 b123) = APS b0 a1 a2 a3 a23 a31 a12 b123
    (ODD a1 a2 a3 a123) + (C b0 b123) = APS b0 a1 a2 a3 0 0 0 (a123 + b123)
    (TPV a23 a31 a12 a123) + (C b0 b123) = APS b0 0 0 0 a23 a31 a12 (a123 + b123)
    (APS a0 a1 a2 a3 a23 a31 a12 a123) + (C b0 b123) = APS (a0 + b0) a1 a2 a3 a23 a31 a12 (a123 + b123)
  
    (BPV a1 a2 a3 a23 a31 a12) + (BPV b1 b2 b3 b23 b31 b12) = BPV (a1 + b1) (a2 + b2) (a3 + b3) (a23 + b23) (a31 + b31) (a12 + b12)
  
    (BPV a1 a2 a3 a23 a31 a12) + (ODD b1 b2 b3 b123) = APS 0 (a1 + b1) (a2 + b2) (a3 + b3) a23 a31 a12 b123
    (BPV a1 a2 a3 a23 a31 a12) + (TPV b23 b31 b12 b123) = APS 0 a1 a2 a3 (a23 + b23) (a31 + b31) (a12 + b12) b123
    (BPV a1 a2 a3 a23 a31 a12) + (APS b0 b1 b2 b3 b23 b31 b12 b123) = APS b0 (a1 + b1) (a2 + b2) (a3 + b3) (a23 + b23) (a31 + b31) (a12 + b12) b123
  
    (ODD a1 a2 a3 a123) + (BPV b1 b2 b3 b23 b31 b12) = APS 0 (a1 + b1) (a2 + b2) (a3 + b3) b23 b31 b12 a123
    (TPV a23 a31 a12 a123) + (BPV b1 b2 b3 b23 b31 b12) = APS 0 b1 b2 b3 (a23 + b23) (a31 + b31) (a12 + b12) a123
    (APS a0 a1 a2 a3 a23 a31 a12 a123) + (BPV b1 b2 b3 b23 b31 b12) = APS a0 (a1 + b1) (a2 + b2) (a3 + b3) (a23 + b23) (a31 + b31) (a12 + b12) a123
  
    (ODD a1 a2 a3 a123) + (ODD b1 b2 b3 b123) = ODD (a1 + b1) (a2 + b2) (a3 + b3) (a123 + b123)
  
    (ODD a1 a2 a3 a123) + (TPV b23 b31 b12 b123) = APS 0 a1 a2 a3 b23 b31 b12 (a123 + b123)
    (ODD a1 a2 a3 a123) + (APS b0 b1 b2 b3 b23 b31 b12 b123) = APS b0 (a1 + b1) (a2 + b2) (a3 + b3) b23 b31 b12 (a123 + b123)
  
    (TPV a23 a31 a12 a123) + (ODD b1 b2 b3 b123) = APS 0 b1 b2 b3 a23 a31 a12 (a123 + b123)
    (APS a0 a1 a2 a3 a23 a31 a12 a123) + (ODD b1 b2 b3 b123) = APS a0 (a1 + b1) (a2 + b2) (a3 + b3) a23 a31 a12 (a123 + b123)
  
    (TPV a23 a31 a12 a123) + (TPV b23 b31 b12 b123) = TPV (a23 + b23) (a31 + b31) (a12 + b12) (a123 + b123)
  
    (TPV a23 a31 a12 a123) + (APS b0 b1 b2 b3 b23 b31 b12 b123) = APS b0 b1 b2 b3 (a23 + b23) (a31 + b31) (a12 + b12) (a123 + b123)
  
    (APS a0 a1 a2 a3 a23 a31 a12 a123) + (TPV b23 b31 b12 b123) = APS a0 a1 a2 a3 (a23 + b23) (a31 + b31) (a12 + b12) (a123 + b123)
  
    (APS a0 a1 a2 a3 a23 a31 a12 a123) + (APS b0 b1 b2 b3 b23 b31 b12 b123) = APS (a0 + b0)
                                                                                  (a1 + b1) (a2 + b2) (a3 + b3)
                                                                                  (a23 + b23) (a31 + b31) (a12 + b12)
                                                                                  (a123 + b123)
  
    -- | Multiplication Instance implementing a Geometric Product
    (R a0) * (R b0) = R (a0*b0)
  
    (R a0) * (V3 b1 b2 b3) = V3 (a0*b1) (a0*b2) (a0*b3)
    (R a0) * (BV b23 b31 b12) = BV (a0*b23) (a0*b31) (a0*b12)
    (R a0) * (I b123) = I (a0*b123)
    (R a0) * (PV b0 b1 b2 b3) = PV (a0*b0)
                                   (a0*b1) (a0*b2) (a0*b3)
    (R a0) * (H b0 b23 b31 b12) = H (a0*b0)
                                    (a0*b23) (a0*b31) (a0*b12)
    (R a0) * (C b0 b123) = C (a0*b0)
                             (a0*b123)
    (R a0) * (BPV b1 b2 b3 b23 b31 b12) = BPV (a0*b1) (a0*b2) (a0*b3)
                                              (a0*b23) (a0*b31) (a0*b12)
    (R a0) * (ODD b1 b2 b3 b123) = ODD (a0*b1) (a0*b2) (a0*b3)
                                       (a0*b123)
    (R a0) * (TPV b23 b31 b12 b123) = TPV (a0*b23) (a0*b31) (a0*b12)
                                          (a0*b123)
    (R a0) * (APS b0 b1 b2 b3 b23 b31 b12 b123) = APS (a0*b0)
                                                      (a0*b1) (a0*b2) (a0*b3)
                                                      (a0*b23) (a0*b31) (a0*b12)
                                                      (a0*b123)
  
    (V3 a1 a2 a3) * (R b0) = V3 (a1*b0) (a2*b0) (a3*b0)
    (BV a23 a31 a12) * (R b0) = BV (a23*b0) (a31*b0) (a12*b0)
    (I a123) * (R b0) = I (a123*b0)
    (PV a0 a1 a2 a3) * (R b0) = PV (a0*b0)
                                   (a1*b0) (a2*b0) (a3*b0)
    (H a0 a23 a31 a12) * (R b0) = H (a0*b0)
                                    (a23*b0) (a31*b0) (a12*b0)
    (C a0 a123) * (R b0) = C (a0*b0)
                             (a123*b0)
    (BPV a1 a2 a3 a23 a31 a12) * (R b0) = BPV (a1*b0) (a2*b0) (a3*b0)
                                              (a23*b0) (a31*b0) (a12*b0)
    (ODD a1 a2 a3 a123) * (R b0) = ODD (a1*b0) (a2*b0) (a3*b0)
                                       (a123*b0)
    (TPV a23 a31 a12 a123) * (R b0) = TPV (a23*b0) (a31*b0) (a12*b0)
                                          (a123*b0)
    (APS a0 a1 a2 a3 a23 a31 a12 a123) * (R b0) = APS (a0*b0)
                                                      (a1*b0) (a2*b0) (a3*b0)
                                                      (a23*b0) (a31*b0) (a12*b0)
                                                      (a123*b0)
  
    (V3 a1 a2 a3) * (V3 b1 b2 b3) = H (a1*b1 + a2*b2 + a3*b3)
                                      (a2*b3 - a3*b2) (a3*b1 - a1*b3) (a1*b2 - a2*b1)
  
    (V3 a1 a2 a3) * (BV b23 b31 b12) = ODD (a3*b31 - a2*b12) (a1*b12 - a3*b23) (a2*b23 - a1*b31)
                                           (a1*b23 + a2*b31 + a3*b12)
    (V3 a1 a2 a3) * (I b123) = BV (a1*b123) (a2*b123) (a3*b123)
    (V3 a1 a2 a3) * (PV b0 b1 b2 b3) = APS (a1*b1 + a2*b2 + a3*b3)
                                           (a1*b0) (a2*b0) (a3*b0)
                                           (a2*b3 - a3*b2) (a3*b1 - a1*b3) (a1*b2 - a2*b1)
                                           0
    (V3 a1 a2 a3) * (H b0 b23 b31 b12) = ODD (a1*b0 - a2*b12 + a3*b31) (a2*b0 + a1*b12 - a3*b23) (a3*b0 - a1*b31 + a2*b23)
                                             (a1*b23 + a2*b31 + a3*b12)
    (V3 a1 a2 a3) * (C b0 b123) = BPV (a1*b0) (a2*b0) (a3*b0)
                                      (a1*b123) (a2*b123) (a3*b123)
    (V3 a1 a2 a3) * (BPV b1 b2 b3 b23 b31 b12) = APS (a1*b1 + a2*b2 + a3*b3)
                                                     (a3*b31 - a2*b12) (a1*b12 - a3*b23) (a2*b23 - a1*b31)
                                                     (a2*b3 - a3*b2) (a3*b1 - a1*b3) (a1*b2 - a2*b1)
                                                     (a1*b23 + a2*b31 + a3*b12)
    (V3 a1 a2 a3) * (ODD b1 b2 b3 b123) = H (a1*b1 + a2*b2 + a3*b3)
                                            (a1*b123 + a2*b3 - a3*b2) (a2*b123 - a1*b3 + a3*b1) (a3*b123 + a1*b2 - a2*b1)
    (V3 a1 a2 a3) * (TPV b23 b31 b12 b123) = APS 0
                                                 (a3*b31 - a2*b12) (a1*b12 - a3*b23) (a2*b23 - a1*b31)
                                                 (a1*b123) (a2*b123) (a3*b123)
                                                 (a1*b23 + a2*b31 + a3*b12)
    (V3 a1 a2 a3) * (APS b0 b1 b2 b3 b23 b31 b12 b123) = APS (a1*b1 + a2*b2 + a3*b3)
                                                             (a1*b0 - a2*b12 + a3*b31) (a2*b0 + a1*b12 - a3*b23) (a3*b0 - a1*b31 + a2*b23)
                                                             (a1*b123 + a2*b3 - a3*b2) (a3*b1 - a1*b3 + a2*b123) (a1*b2 - a2*b1 + a3*b123)
                                                             (a1*b23 + a2*b31 + a3*b12)
  
    (BV a23 a31 a12) * (V3 b1 b2 b3) = ODD (a12*b2  - a31*b3) (a23*b3 - a12*b1) (a31*b1  - a23*b2)
                                           (a23*b1  + a31*b2  + a12*b3)
    (I a123) * (V3 b1 b2 b3) = BV (a123*b1) (a123*b2) (a123*b3)
    (PV a0 a1 a2 a3) * (V3 b1 b2 b3) = APS (a1*b1 + a2*b2 + a3*b3)
                                           (a0*b1) (a0*b2) (a0*b3)
                                           (a2*b3 - a3*b2) (a3*b1 - a1*b3) (a1*b2 - a2*b1)
                                           0
    (H a0 a23 a31 a12) * (V3 b1 b2 b3) = ODD (a0*b1 + a12*b2 - a31*b3) (a0*b2 - a12*b1 + a23*b3) (a0*b3 + a31*b1 - a23*b2)
                                             (a23*b1 + a31*b2 + a12*b3)
    (C a0 a123) * (V3 b1 b2 b3) = BPV (a0*b1) (a0*b2) (a0*b3)
                                      (a123*b1) (a123*b2) (a123*b3)
    (BPV a1 a2 a3 a23 a31 a12) * (V3 b1 b2 b3) = APS (a1*b1 + a2*b2 + a3*b3)
                                                     (a12*b2 - a31*b3) (a23*b3 - a12*b1) (a31*b1 - a23*b2)
                                                     (a2*b3 - a3*b2) (a3*b1 - a1*b3) (a1*b2 - a2*b1)
                                                     (a23*b1 + a31*b2 + a12*b3)
    (ODD a1 a2 a3 a123) * (V3 b1 b2 b3) = H (a1*b1 + a2*b2 + a3*b3)
                                            (a123*b1 + a2*b3 - a3*b2) (a123*b2 - a1*b3 + a3*b1) (a123*b3 + a1*b2 - a2*b1)
    (TPV a23 a31 a12 a123) * (V3 b1 b2 b3) = APS 0
                                                 (a12*b2 - a31*b3) (a23*b3 - a12*b1) (a31*b1 - a23*b2)
                                                 (a123*b1) (a123*b2) (a123*b3)
                                                 (a23*b1 + a31*b2 + a12*b3)
    (APS a0 a1 a2 a3 a23 a31 a12 a123) * (V3 b1 b2 b3) = APS (a1*b1 + a2*b2 + a3*b3)
                                                             (a0*b1 + a12*b2 - a31*b3) (a0*b2 - a12*b1 + a23*b3) (a0*b3 + a31*b1 - a23*b2)
                                                             (a123*b1 + a2*b3 - a3*b2) (a3*b1 - a1*b3 + a123*b2) (a1*b2 - a2*b1 + a123*b3)
                                                             (a23*b1 + a31*b2 + a12*b3)
  
    (BV a23 a31 a12) * (BV b23 b31 b12) = H (negate $ a23*b23 + a31*b31 + a12*b12)
                                            (a12*b31 - a31*b12) (a23*b12 - a12*b23) (a31*b23 - a23*b31)
  
    (BV a23 a31 a12) * (I b123) = V3 (negate $ a23*b123) (negate $ a31*b123) (negate $ a12*b123)
    (BV a23 a31 a12) * (PV b0 b1 b2 b3) = APS 0
                                              (a12*b2 - a31*b3) (a23*b3 - a12*b1) (a31*b1 - a23*b2)
                                              (a23*b0) (a31*b0) (a12*b0)
                                              (a23*b1 + a31*b2 + a12*b3)
    (BV a23 a31 a12) * (H b0 b23 b31 b12) = H (negate $ a23*b23 + a31*b31 + a12*b12)
                                              (a23*b0 - a31*b12 + a12*b31) (a31*b0 + a23*b12 - a12*b23) (a12*b0 - a23*b31 + a31*b23)
    (BV a23 a31 a12) * (C b0 b123) = BPV (negate $ a23*b123) (negate $ a31*b123) (negate $ a12*b123)
                                         (a23*b0) (a31*b0) (a12*b0)
    (BV a23 a31 a12) * (BPV b1 b2 b3 b23 b31 b12) = APS (negate $ a23*b23 + a31*b31 + a12*b12)
                                                        (a12*b2 - a31*b3) (a23*b3 - a12*b1) (a31*b1 - a23*b2)  
                                                        (a12*b31 - a31*b12) (a23*b12 - a12*b23) (a31*b23 - a23*b31)
                                                        (a23*b1 + a31*b2 + a12*b3)
    (BV a23 a31 a12) * (ODD b1 b2 b3 b123) = ODD (a12*b2 - a31*b3 - a23*b123) (a23*b3 - a12*b1 - a31*b123) (a31*b1 - a23*b2 - a12*b123)
                                                 (a23*b1 + a31*b2 + a12*b3)
    (BV a23 a31 a12) * (TPV b23 b31 b12 b123) = APS (negate $ a23*b23 + a31*b31 + a12*b12)
                                                    (negate $ a23*b123) (negate $ a31*b123) (negate $ a12*b123)
                                                    (a12*b31 - a31*b12) (a23*b12 - a12*b23) (a31*b23 - a23*b31)
                                                    0
    (BV a23 a31 a12) * (APS b0 b1 b2 b3 b23 b31 b12 b123) = APS (negate $ a23*b23 + a31*b31 + a12*b12)
                                                                (a12*b2 - a31*b3 - a23*b123) (a23*b3 - a31*b123 - a12*b1) (a31*b1 - a23*b2 - a12*b123)
                                                                (a23*b0 - a31*b12 + a12*b31) (a31*b0 + a23*b12 - a12*b23) (a12*b0 - a23*b31 + a31*b23)
                                                                (a23*b1 + a31*b2 + a12*b3)
  
    (I a123) * (BV b23 b31 b12) = V3 (negate $ a123*b23) (negate $ a123*b31) (negate $ a123*b12)
    (PV a0 a1 a2 a3) * (BV b23 b31 b12) = APS 0
                                              (a3*b31 - a2*b12) (a1*b12 - a3*b23) (a2*b23 - a1*b31)
                                              (a0*b23) (a0*b31) (a0*b12)
                                              (a1*b23 + a2*b31 + a3*b12)
    (H a0 a23 a31 a12) * (BV b23 b31 b12) = H (negate $ a23*b23 + a31*b31 + a12*b12)
                                              (a0*b23 - a31*b12 + a12*b31) (a0*b31 + a23*b12 - a12*b23) (a0*b12 - a23*b31 + a31*b23)
    (C a0 a123) * (BV b23 b31 b12) = BPV (negate $ a123*b23) (negate $ a123*b31) (negate $ a123*b12)
                                         (a0*b23) (a0*b31) (a0*b12)
    (BPV a1 a2 a3 a23 a31 a12) * (BV b23 b31 b12) = APS (negate $ a23*b23 + a31*b31 + a12*b12)
                                                        (a3*b31 - a2*b12) (a1*b12 - a3*b23) (a2*b23 - a1*b31)    
                                                        (a12*b31 - a31*b12) (a23*b12 - a12*b23) (a31*b23 - a23*b31)
                                                        (a1*b23 + a2*b31 + a3*b12)
    (ODD a1 a2 a3 a123) * (BV b23 b31 b12) = ODD (negate $ a123*b23 + a2*b12 - a3*b31)
                                                 (negate $ a123*b31 - a1*b12 + a3*b23)
                                                 (negate $ a123*b12 + a1*b31 - a2*b23)
                                                 (a1*b23 + a2*b31 + a3*b12)
    (TPV a23 a31 a12 a123) * (BV b23 b31 b12) = APS (negate $ a23*b23 + a31*b31 + a12*b12)
                                                    (negate $ a123*b23) (negate $ a123*b31) (negate $ a123*b12)
                                                    (negate $ a31*b12 - a12*b31) (negate $ a12*b23 - a23*b12) (negate $ a23*b31 - a31*b23)
                                                    0
    (APS a0 a1 a2 a3 a23 a31 a12 a123) * (BV b23 b31 b12) = APS (negate $ a23*b23 + a31*b31 + a12*b12)  
                                                                (a3*b31 - a123*b23 - a2*b12) (a1*b12 - a3*b23 - a123*b31) (a2*b23 - a123*b12 - a1*b31)
                                                                (a0*b23 - a31*b12 + a12*b31) (a0*b31 + a23*b12 - a12*b23) (a0*b12 - a23*b31 + a31*b23)
                                                                (a1*b23 + a2*b31 + a3*b12)
  
    (I a123) * (I b123) = R (negate $ a123*b123)
  
    (I a123) * (PV b0 b1 b2 b3) = TPV (a123*b1) (a123*b2) (a123*b3)
                                      (a123*b0)
    (I a123) * (H b0 b23 b31 b12) = ODD (negate $ a123*b23) (negate $ a123*b31) (negate $ a123*b12)
                                        (a123*b0)
    (I a123) * (C b0 b123) = C (negate $ a123*b123)
                               (a123*b0)
    (I a123) * (BPV b1 b2 b3 b23 b31 b12) = BPV (negate $ a123*b23) (negate $ a123*b31) (negate $ a123*b12)
                                                (a123*b1) (a123*b2) (a123*b3)
    (I a123) * (ODD b1 b2 b3 b123) = H (negate $ a123*b123)
                                       (a123*b1) (a123*b2) (a123*b3)
    (I a123) * (TPV b23 b31 b12 b123) = PV (negate $ a123*b123)
                                           (negate $ a123*b23) (negate $ a123*b31) (negate $ a123*b12)
    (I a123) * (APS b0 b1 b2 b3 b23 b31 b12 b123) = APS (negate $ a123*b123)
                                                        (negate $ a123*b23) (negate $ a123*b31) (negate $ a123*b12)
                                                        (a123*b1) (a123*b2) (a123*b3)
                                                        (a123*b0)
  
    (PV a0 a1 a2 a3) * (I b123) = TPV (a1*b123) (a2*b123) (a3*b123)
                                      (a0*b123)
    (H a0 a23 a31 a12) * (I b123) = ODD (negate $ a23*b123) (negate $ a31*b123) (negate $ a12*b123)
                                        (a0*b123)
    (C a0 a123) * (I b123) = C (negate $ a123*b123)
                               (a0*b123)
    (BPV a1 a2 a3 a23 a31 a12) * (I b123) = BPV (negate $ a23*b123) (negate $ a31*b123) (negate $ a12*b123)
                                                (a1*b123) (a2*b123) (a3*b123)
    (ODD a1 a2 a3 a123) * (I b123) = H (negate $ a123*b123)
                                       (a1*b123) (a2*b123) (a3*b123)
    (TPV a23 a31 a12 a123) * (I b123) = PV (negate $ a123*b123)
                                           (negate $ a23*b123) (negate $ a31*b123) (negate $ a12*b123)
    (APS a0 a1 a2 a3 a23 a31 a12 a123) * (I b123) = APS (negate $ a123*b123)
                                                        (negate $ a23*b123) (negate $ a31*b123) (negate $ a12*b123)
                                                        (a1*b123) (a2*b123) (a3*b123)
                                                        (a0*b123)
  
  
    (PV a0 a1 a2 a3) * (PV b0 b1 b2 b3) = APS (a0*b0 + a1*b1 + a2*b2 + a3*b3)
                                              (a0*b1 + a1*b0) (a0*b2 + a2*b0) (a0*b3 + a3*b0)
                                              (a2*b3 - a3*b2) (a3*b1 - a1*b3) (a1*b2 - a2*b1)
                                              0
  
    (PV a0 a1 a2 a3) * (H b0 b23 b31 b12) = APS (a0*b0)
                                                (a1*b0 - a2*b12 + a3*b31) (a2*b0 + a1*b12 - a3*b23) (a3*b0 - a1*b31 + a2*b23)
                                                (a0*b23) (a0*b31) (a0*b12)
                                                (a1*b23 + a2*b31 + a3*b12)
    (PV a0 a1 a2 a3) * (C b0 b123) = APS (a0*b0)
                                         (a1*b0) (a2*b0) (a3*b0)
                                         (a1*b123) (a2*b123) (a3*b123)
                                         (a0*b123)
    (PV a0 a1 a2 a3) * (BPV b1 b2 b3 b23 b31 b12) = APS (a1*b1 + a2*b2 + a3*b3)
                                                        (a0*b1 - a2*b12 + a3*b31) (a0*b2 + a1*b12 - a3*b23) (a0*b3 - a1*b31 + a2*b23)
                                                        (a0*b23 + a2*b3 - a3*b2) (a0*b31 - a1*b3 + a3*b1) (a0*b12 + a1*b2 - a2*b1)
                                                        (a1*b23 + a2*b31 + a3*b12)
    (PV a0 a1 a2 a3) * (ODD b1 b2 b3 b123) = APS (a1*b1 + a2*b2 + a3*b3)
                                                 (a0*b1) (a0*b2) (a0*b3)
                                                 (a1*b123 + a2*b3 - a3*b2) (a2*b123 - a1*b3 + a3*b1) (a3*b123 + a1*b2 - a2*b1)
                                                 (a0*b123)
    (PV a0 a1 a2 a3) * (TPV b23 b31 b12 b123) = APS 0
                                                    (a3*b31 - a2*b12) (a1*b12 - a3*b23) (a2*b23 - a1*b31)
                                                    (a0*b23 + a1*b123) (a0*b31 + a2*b123) (a0*b12 + a3*b123)
                                                    (a0*b123 + a1*b23 + a2*b31 + a3*b12)
    (PV a0 a1 a2 a3) * (APS b0 b1 b2 b3 b23 b31 b12 b123) = APS (a0*b0 + a1*b1 + a2*b2 + a3*b3)
                                                                (a0*b1 + a1*b0 - a2*b12 + a3*b31)
                                                                (a0*b2 + a2*b0 + a1*b12 - a3*b23)
                                                                (a0*b3 + a3*b0 - a1*b31 + a2*b23)
                                                                (a0*b23 + a1*b123 + a2*b3 - a3*b2)
                                                                (a0*b31 - a1*b3 + a3*b1 + a2*b123)
                                                                (a0*b12 + a1*b2 - a2*b1 + a3*b123)
                                                                (a0*b123 + a1*b23 + a2*b31 + a3*b12)
  
    (H a0 a23 a31 a12) * (PV b0 b1 b2 b3) = APS (a0*b0)
                                                (a0*b1 + a12*b2 - a31*b3) (a0*b2 - a12*b1 + a23*b3) (a0*b3 + a31*b1 - a23*b2)
                                                (a23*b0) (a31*b0) (a12*b0)
                                                (a23*b1 + a31*b2 + a12*b3)
    (C a0 a123) * (PV b0 b1 b2 b3) = APS (a0*b0)
                                         (a0*b1) (a0*b2) (a0*b3)
                                         (a123*b1) (a123*b2) (a123*b3)
                                         (a123*b0)
    (BPV a1 a2 a3 a23 a31 a12) * (PV b0 b1 b2 b3) = APS (a1*b1 + a2*b2 + a3*b3)
                                                        (a1*b0 + a12*b2 - a31*b3) (a2*b0 - a12*b1 + a23*b3) (a3*b0 + a31*b1 - a23*b2)
                                                        (a23*b0 + a2*b3 - a3*b2) (a31*b0 - a1*b3 + a3*b1) (a12*b0 + a1*b2 - a2*b1)
                                                        (a23*b1 + a31*b2 + a12*b3)
    (ODD a1 a2 a3 a123) * (PV b0 b1 b2 b3) = APS (a1*b1 + a2*b2 + a3*b3)
                                                 (a1*b0) (a2*b0) (a3*b0)
                                                 (a123*b1 + a2*b3 - a3*b2)
                                                 (a123*b2 - a1*b3 + a3*b1)
                                                 (a123*b3 + a1*b2 - a2*b1)
                                                 (a123*b0)
    (TPV a23 a31 a12 a123) * (PV b0 b1 b2 b3) = APS 0
                                                    (a12*b2 - a31*b3) (a23*b3 - a12*b1) (a31*b1 - a23*b2)
                                                    (a23*b0 + a123*b1) (a31*b0 + a123*b2) (a12*b0 + a123*b3)
                                                    (a123*b0 + a23*b1 + a31*b2 + a12*b3)
    (APS a0 a1 a2 a3 a23 a31 a12 a123) * (PV b0 b1 b2 b3) = APS (a0*b0 + a1*b1 + a2*b2 + a3*b3)
                                                                (a0*b1 + a1*b0 + a12*b2 - a31*b3)
                                                                (a0*b2 + a2*b0 - a12*b1 + a23*b3)
                                                                (a0*b3 + a3*b0 + a31*b1 - a23*b2)
                                                                (a23*b0 + a123*b1 + a2*b3 - a3*b2)
                                                                (a31*b0 - a1*b3 + a3*b1 + a123*b2)
                                                                (a12*b0 + a1*b2 - a2*b1 + a123*b3)
                                                                (a123*b0 + a23*b1 + a31*b2 + a12*b3)
  
    (H a0 a23 a31 a12) * (H b0 b23 b31 b12) = H (a0*b0 - a23*b23 - a31*b31 - a12*b12)
                                                (a0*b23 + a23*b0 - a31*b12 + a12*b31)
                                                (a0*b31 + a31*b0 + a23*b12 - a12*b23)
                                                (a0*b12 + a12*b0 - a23*b31 + a31*b23)
  
    (H a0 a23 a31 a12) * (C b0 b123) = APS (a0*b0)
                                           (negate $ a23*b123) (negate $ a31*b123) (negate $ a12*b123)
                                           (a23*b0) (a31*b0) (a12*b0)
                                           (a0*b123)
    (H a0 a23 a31 a12) * (BPV b1 b2 b3 b23 b31 b12) = APS (negate $ a23*b23 + a31*b31 + a12*b12)
                                                          (a0*b1 + a12*b2 - a31*b3) (a0*b2 - a12*b1 + a23*b3) (a0*b3 + a31*b1 - a23*b2)
                                                          (a0*b23 - a31*b12 + a12*b31) (a0*b31 + a23*b12 - a12*b23) (a0*b12 - a23*b31 + a31*b23)
                                                          (a23*b1 + a31*b2  + a12*b3)
    (H a0 a23 a31 a12) * (ODD b1 b2 b3 b123) = ODD (a0*b1 + a12*b2 - a31*b3 - a23*b123)
                                                   (a0*b2 - a12*b1 + a23*b3 - a31*b123)
                                                   (a0*b3 + a31*b1 - a23*b2 - a12*b123)
                                                   (a0*b123 + a23*b1 + a31*b2 + a12*b3)
    (H a0 a23 a31 a12) * (TPV b23 b31 b12 b123) = APS (negate $ a23*b23 + a31*b31 + a12*b12)
                                                      (negate $ a23*b123) (negate $ a31*b123) (negate $ a12*b123)
                                                      (a0*b23 - a31*b12 + a12*b31) (a0*b31 + a23*b12 - a12*b23) (a0*b12 - a23*b31 + a31*b23)
                                                      (a0*b123)
    (H a0 a23 a31 a12) * (APS b0 b1 b2 b3 b23 b31 b12 b123) = APS (a0*b0 - a23*b23 - a31*b31 - a12*b12)
                                                                  (a0*b1 + a12*b2 - a31*b3 - a23*b123)
                                                                  (a0*b2 - a12*b1 + a23*b3 - a31*b123)
                                                                  (a0*b3 + a31*b1 - a23*b2 - a12*b123)
                                                                  (a0*b23 + a23*b0 - a31*b12 + a12*b31)
                                                                  (a0*b31 + a31*b0 + a23*b12 - a12*b23)
                                                                  (a0*b12 + a12*b0 - a23*b31 + a31*b23)
                                                                  (a0*b123 + a23*b1 + a31*b2 + a12*b3)
  
    (C a0 a123) * (H b0 b23 b31 b12) = APS (a0*b0)
                                           (negate $ a123*b23) (negate $ a123*b31) (negate $ a123*b12)
                                           (a0*b23) (a0*b31) (a0*b12)
                                           (a123*b0)
    (BPV a1 a2 a3 a23 a31 a12) * (H b0 b23 b31 b12) = APS (negate $ a23*b23 + a31*b31 + a12*b12)
                                                          (a1*b0 - a2*b12 + a3*b31) (a2*b0 + a1*b12 - a3*b23) (a3*b0 - a1*b31 + a2*b23)
                                                          (a23*b0 - a31*b12 + a12*b31) (a31*b0 + a23*b12 - a12*b23) (a12*b0 - a23*b31 + a31*b23)
                                                          (a1*b23 + a2*b31 + a3*b12)
    (ODD a1 a2 a3 a123) * (H b0 b23 b31 b12) = ODD (a1*b0 - a2*b12 + a3*b31 - a123*b23)
                                                   (a2*b0 + a1*b12 - a3*b23 - a123*b31)
                                                   (a3*b0 - a1*b31 + a2*b23 - a123*b12)
                                                   (a123*b0 + a1*b23 + a2*b31 + a3*b12)
    (TPV a23 a31 a12 a123) * (H b0 b23 b31 b12) = APS (negate $ a23*b23 + a31*b31 + a12*b12)
                                                      (negate $ a123*b23) (negate $ a123*b31) (negate $ a123*b12)
                                                      (a23*b0 - a31*b12 + a12*b31) (a31*b0 + a23*b12 - a12*b23) (a12*b0 - a23*b31 + a31*b23)
                                                      (a123*b0)
    (APS a0 a1 a2 a3 a23 a31 a12 a123) * (H b0 b23 b31 b12) = APS (a0*b0 - a23*b23 - a31*b31 - a12*b12)
                                                                  (a1*b0 - a2*b12 + a3*b31 - a123*b23)
                                                                  (a2*b0 + a1*b12 - a3*b23 - a123*b31)
                                                                  (a3*b0 - a1*b31 + a2*b23 - a123*b12)
                                                                  (a0*b23 + a23*b0 - a31*b12 + a12*b31)
                                                                  (a0*b31 + a31*b0 + a23*b12 - a12*b23)
                                                                  (a0*b12 + a12*b0 - a23*b31 + a31*b23)
                                                                  (a123*b0 + a1*b23 + a2*b31 + a3*b12)
  
    (C a0 a123) * (C b0 b123) = C (a0*b0 - a123*b123)
                                  (a0*b123 + a123*b0)
  
    (C a0 a123) * (BPV b1 b2 b3 b23 b31 b12) = BPV (a0*b1 - a123*b23) (a0*b2 - a123*b31) (a0*b3 - a123*b12)
                                                   (a0*b23 + a123*b1) (a0*b31 + a123*b2) (a0*b12 + a123*b3)
    (C a0 a123) * (ODD b1 b2 b3 b123) = APS (negate $ a123*b123)
                                            (a0*b1) (a0*b2) (a0*b3)
                                            (a123*b1) (a123*b2) (a123*b3)
                                            (a0*b123)
    (C a0 a123) * (TPV b23 b31 b12 b123) = APS (negate $ a123*b123)
                                               (negate $ a123*b23) (negate $ a123*b31) (negate $ a123*b12)
                                               (a0*b23) (a0*b31) (a0*b12)
                                               (a0*b123)
    (C a0 a123) * (APS b0 b1 b2 b3 b23 b31 b12 b123) = APS (a0*b0 - a123*b123)
                                                           (a0*b1 - a123*b23) (a0*b2 - a123*b31) (a0*b3 - a123*b12)
                                                           (a0*b23 + a123*b1) (a0*b31 + a123*b2) (a0*b12 + a123*b3)
                                                           (a0*b123 + a123*b0)
  
    (BPV a1 a2 a3 a23 a31 a12) * (C b0 b123) = BPV (a1*b0 - a23*b123) (a2*b0 - a31*b123) (a3*b0 - a12*b123)
                                                   (a23*b0 + a1*b123) (a31*b0 + a2*b123) (a12*b0 + a3*b123)
    (ODD a1 a2 a3 a123) * (C b0 b123) = APS (negate $ a123*b123)
                                            (a1*b0) (a2*b0) (a3*b0)
                                            (a1*b123) (a2*b123) (a3*b123)
                                            (a123*b0)
    (TPV a23 a31 a12 a123) * (C b0 b123) = APS (negate $ a123*b123)
                                               (negate $ a23*b123) (negate $ a31*b123) (negate $ a12*b123)
                                               (a23*b0) (a31*b0) (a12*b0)
                                               (a123*b0)
    (APS a0 a1 a2 a3 a23 a31 a12 a123) * (C b0 b123) = APS (a0*b0 - a123*b123)
                                                           (a1*b0 - a23*b123) (a2*b0 - a31*b123) (a3*b0 - a12*b123)
                                                           (a23*b0 + a1*b123) (a31*b0 + a2*b123) (a12*b0 + a3*b123)
                                                           (a0*b123 + a123*b0)
  
    (BPV a1 a2 a3 a23 a31 a12) * (BPV b1 b2 b3 b23 b31 b12) = APS (a1*b1 + a2*b2 + a3*b3 - a23*b23 - a31*b31 - a12*b12)
                                                                  (a12*b2 - a2*b12 + a3*b31 - a31*b3)
                                                                  (a1*b12 - a12*b1 - a3*b23 + a23*b3)
                                                                  (a31*b1 - a1*b31 + a2*b23 - a23*b2)
                                                                  (a2*b3 - a3*b2 - a31*b12 + a12*b31)
                                                                  (a3*b1 - a1*b3 + a23*b12 - a12*b23)
                                                                  (a1*b2 - a2*b1 - a23*b31 + a31*b23)
                                                                  (a1*b23 + a23*b1 + a2*b31 + a31*b2 + a3*b12 + a12*b3)
  
    (BPV a1 a2 a3 a23 a31 a12) * (ODD b1 b2 b3 b123) = APS (a1*b1 + a2*b2 + a3*b3)
                                                           (a12*b2 - a31*b3 - a23*b123) (a23*b3 - a12*b1 - a31*b123) (a31*b1 - a23*b2 - a12*b123)
                                                           (a1*b123 + a2*b3 - a3*b2) (a2*b123 - a1*b3 + a3*b1) (a3*b123 + a1*b2 - a2*b1)
                                                           (a23*b1 + a31*b2 + a12*b3)
    (BPV a1 a2 a3 a23 a31 a12) * (TPV b23 b31 b12 b123) = APS (negate $ a23*b23 + a31*b31 + a12*b12)
                                                              (a3*b31 - a2*b12 - a23*b123) (a1*b12 - a3*b23 - a31*b123) (a2*b23 - a1*b31 - a12*b123)
                                                              (a1*b123 - a31*b12 + a12*b31) (a2*b123 + a23*b12 - a12*b23) (a3*b123 - a23*b31 + a31*b23)
                                                              (a1*b23 + a2*b31 + a3*b12)
    (BPV a1 a2 a3 a23 a31 a12) * (APS b0 b1 b2 b3 b23 b31 b12 b123) = APS (a1*b1 + a2*b2 + a3*b3 - a23*b23 - a31*b31 - a12*b12)
                                                                          (a1*b0 - a2*b12 + a12*b2 + a3*b31 - a31*b3 - a23*b123)
                                                                          (a2*b0 + a1*b12 - a12*b1 - a3*b23 + a23*b3 - a31*b123)
                                                                          (a3*b0 - a1*b31 + a31*b1 + a2*b23 - a23*b2 - a12*b123)
                                                                          (a23*b0 + a1*b123 + a2*b3 - a3*b2 - a31*b12 + a12*b31)
                                                                          (a31*b0 - a1*b3 + a3*b1 + a2*b123 + a23*b12 - a12*b23)
                                                                          (a12*b0 + a1*b2 - a2*b1 + a3*b123 - a23*b31 + a31*b23)
                                                                          (a1*b23 + a23*b1 + a2*b31 + a31*b2 + a3*b12 + a12*b3)
  
    (ODD a1 a2 a3 a123) * (BPV b1 b2 b3 b23 b31 b12) = APS (a1*b1 + a2*b2 + a3*b3)
                                                           (a3*b31 - a2*b12 - a123*b23) (a1*b12 - a3*b23 - a123*b31) (a2*b23 - a1*b31 - a123*b12)
                                                           (a123*b1 + a2*b3 - a3*b2) (a123*b2 - a1*b3 + a3*b1) (a123*b3 + a1*b2 - a2*b1)
                                                           (a1*b23 + a2*b31 + a3*b12)
    (TPV a23 a31 a12 a123) * (BPV b1 b2 b3 b23 b31 b12) = APS (negate $ a23*b23 + a31*b31 + a12*b12)
                                                              (a12*b2 - a31*b3 - a123*b23) (a23*b3 - a12*b1 - a123*b31) (a31*b1 - a23*b2 - a123*b12)
                                                              (a123*b1 - a31*b12 + a12*b31) (a123*b2 + a23*b12 - a12*b23) (a123*b3 - a23*b31 + a31*b23)
                                                              (a23*b1 + a31*b2 + a12*b3)
    (APS a0 a1 a2 a3 a23 a31 a12 a123) * (BPV b1 b2 b3 b23 b31 b12) = APS (a1*b1 + a2*b2 + a3*b3 - a23*b23 - a31*b31 - a12*b12)
                                                                          (a0*b1 - a2*b12 + a12*b2 + a3*b31 - a31*b3 - a123*b23)
                                                                          (a0*b2 + a1*b12 - a12*b1 - a3*b23 + a23*b3 - a123*b31)
                                                                          (a0*b3 - a1*b31 + a31*b1 + a2*b23 - a23*b2 - a123*b12)
                                                                          (a0*b23 + a123*b1 + a2*b3 - a3*b2 - a31*b12 + a12*b31)
                                                                          (a0*b31 - a1*b3 + a3*b1 + a123*b2 + a23*b12 - a12*b23)
                                                                          (a0*b12 + a1*b2 - a2*b1 + a123*b3 - a23*b31 + a31*b23)
                                                                          (a1*b23 + a23*b1 + a2*b31 + a31*b2 + a3*b12 + a12*b3)
  
    (ODD a1 a2 a3 a123) * (ODD b1 b2 b3 b123) = H (a1*b1 + a2*b2 + a3*b3 - a123*b123)
                                                  (a1*b123 + a123*b1 + a2*b3 - a3*b2)
                                                  (a2*b123 + a123*b2 - a1*b3 + a3*b1)
                                                  (a3*b123 + a123*b3 + a1*b2 - a2*b1)
  
    (ODD a1 a2 a3 a123) * (TPV b23 b31 b12 b123) = APS (negate $ a123*b123)
                                                       (a3*b31 - a2*b12 - a123*b23) (a1*b12 - a3*b23 - a123*b31) (a2*b23 - a1*b31 - a123*b12)
                                                       (a1*b123) (a2*b123) (a3*b123)
                                                       (a1*b23 + a2*b31 + a3*b12)
    (ODD a1 a2 a3 a123) * (APS b0 b1 b2 b3 b23 b31 b12 b123) = APS (a1*b1 + a2*b2 + a3*b3 - a123*b123)
                                                                   (a1*b0 - a2*b12 + a3*b31 - a123*b23)
                                                                   (a2*b0 + a1*b12 - a3*b23 - a123*b31)
                                                                   (a3*b0 - a1*b31 + a2*b23 - a123*b12)
                                                                   (a1*b123 + a123*b1 + a2*b3 - a3*b2)
                                                                   (a2*b123 + a123*b2 - a1*b3 + a3*b1)
                                                                   (a3*b123 + a123*b3 + a1*b2 - a2*b1)
                                                                   (a123*b0 + a1*b23 + a2*b31 + a3*b12)
  
    (TPV a23 a31 a12 a123) * (ODD b1 b2 b3 b123) = APS (negate $ a123*b123)
                                                       (a12*b2 - a31*b3 - a23*b123) (a23*b3 - a12*b1 - a31*b123) (a31*b1 - a23*b2 - a12*b123)
                                                       (a123*b1) (a123*b2) (a123*b3)
                                                       (a23*b1 + a31*b2 + a12*b3)
    (APS a0 a1 a2 a3 a23 a31 a12 a123) * (ODD b1 b2 b3 b123) = APS (a1*b1 + a2*b2 + a3*b3 - a123*b123)
                                                                   (a0*b1 + a12*b2 - a31*b3 - a23*b123)
                                                                   (a0*b2 - a12*b1 + a23*b3 - a31*b123)
                                                                   (a0*b3 + a31*b1 - a23*b2 - a12*b123)
                                                                   (a1*b123 + a123*b1 + a2*b3 - a3*b2)
                                                                   (a2*b123 + a123*b2 - a1*b3 + a3*b1)
                                                                   (a3*b123 + a123*b3 + a1*b2 - a2*b1)
                                                                   (a0*b123 + a23*b1 + a31*b2 + a12*b3)
  
    (TPV a23 a31 a12 a123) * (TPV b23 b31 b12 b123) = APS (negate $ a23*b23 + a31*b31 + a12*b12 + a123*b123)
                                                          (negate $ a23*b123 + a123*b23) (negate $ a31*b123 + a123*b31) (negate $ a12*b123 + a123*b12)
                                                          (a12*b31 - a31*b12) (a23*b12 - a12*b23) (a31*b23 - a23*b31)
                                                          0
  
    (TPV a23 a31 a12 a123) * (APS b0 b1 b2 b3 b23 b31 b12 b123) = APS (negate $ a23*b23 + a31*b31 + a12*b12 + a123*b123)
                                                                      (a12*b2 - a31*b3 - a23*b123 - a123*b23)
                                                                      (a23*b3 - a12*b1 - a31*b123 - a123*b31)
                                                                      (a31*b1 - a23*b2 - a12*b123 - a123*b12)
                                                                      (a23*b0 + a123*b1 - a31*b12 + a12*b31)
                                                                      (a31*b0 + a123*b2 + a23*b12 - a12*b23)
                                                                      (a12*b0 + a123*b3 - a23*b31 + a31*b23)
                                                                      (a123*b0 + a23*b1 + a31*b2 + a12*b3)
  
    (APS a0 a1 a2 a3 a23 a31 a12 a123) * (TPV b23 b31 b12 b123) = APS (negate $ a23*b23 + a31*b31 + a12*b12 + a123*b123)
                                                                      (a3*b31 - a2*b12 - a23*b123 - a123*b23)
                                                                      (a1*b12 - a3*b23 - a31*b123 - a123*b31)
                                                                      (a2*b23 - a1*b31 - a12*b123 - a123*b12)
                                                                      (a0*b23 + a1*b123 - a31*b12 + a12*b31)
                                                                      (a0*b31 + a2*b123 + a23*b12 - a12*b23)
                                                                      (a0*b12 + a3*b123 - a23*b31 + a31*b23)
                                                                      (a0*b123 + a1*b23 + a2*b31 + a3*b12)
  
    (APS a0 a1 a2 a3 a23 a31 a12 a123) * (APS b0 b1 b2 b3 b23 b31 b12 b123) = APS (a0*b0 + a1*b1 + a2*b2 + a3*b3 - a23*b23 - a31*b31 - a12*b12 - a123*b123)
                                                                                  (a0*b1 + a1*b0 - a2*b12 + a12*b2 + a3*b31 - a31*b3 - a23*b123 - a123*b23)
                                                                                  (a0*b2 + a2*b0 + a1*b12 - a12*b1 - a3*b23 + a23*b3 - a31*b123 - a123*b31)
                                                                                  (a0*b3 + a3*b0 - a1*b31 + a31*b1 + a2*b23 - a23*b2 - a12*b123 - a123*b12)
                                                                                  (a0*b23 + a23*b0 + a1*b123 + a123*b1 + a2*b3 - a3*b2 - a31*b12 + a12*b31)
                                                                                  (a0*b31 + a31*b0 - a1*b3 + a3*b1 + a2*b123 + a123*b2 + a23*b12 - a12*b23)
                                                                                  (a0*b12 + a12*b0 + a1*b2 - a2*b1 + a3*b123 + a123*b3 - a23*b31 + a31*b23)
                                                                                  (a0*b123 + a123*b0 + a1*b23 + a23*b1 + a2*b31 + a31*b2 + a3*b12 + a12*b3)
  
  
  
  
    -- |'abs' is the spectral norm aka the spectral radius
    -- it is the largest singular value. This function may need to be fiddled with
    -- to make the math a bit safer wrt overflows.  This makes use of the largest
    -- singular value, if the littlest singular value is zero then the element is not
    -- invertable, we can see here that R, C, V3, BV, and H are all invertable, and
    -- by implication R, C, and H are division algebras.
    abs cl3 = fst $ abssignum cl3
  
  
    -- |'signum' satisfies the Law "abs x * signum x == x"
    -- kind of cool: signum of a vector is it's unit vector.
    signum cl3 = snd $ abssignum cl3
  
  
    -- |'fromInteger'
    fromInteger int = R (fromInteger int)
  
  
    -- |'negate' simply distributes into the grade components
    negate (R a0) = R (negate a0)
    negate (V3 a1 a2 a3) = V3 (negate a1) (negate a2) (negate a3)
    negate (BV a23 a31 a12) = BV (negate a23) (negate a31) (negate a12)
    negate (I a123) = I (negate a123)
    negate (PV a0 a1 a2 a3) =  PV (negate a0)
                                  (negate a1) (negate a2) (negate a3)
    negate (H a0 a23 a31 a12) = H (negate a0)
                                  (negate a23) (negate a31) (negate a12)
    negate (C a0 a123) = C (negate a0)
                           (negate a123)
    negate (BPV a1 a2 a3 a23 a31 a12) = BPV (negate a1) (negate a2) (negate a3)
                                            (negate a23) (negate a31) (negate a12)
    negate (ODD a1 a2 a3 a123) = ODD (negate a1) (negate a2) (negate a3)
                                     (negate a123)
    negate (TPV a23 a31 a12 a123) = TPV (negate a23) (negate a31) (negate a12)
                                        (negate a123)
    negate (APS a0 a1 a2 a3 a23 a31 a12 a123) = APS (negate a0)
                                                    (negate a1) (negate a2) (negate a3)
                                                    (negate a23) (negate a31) (negate a12)
                                                    (negate a123)
  
  -- | 'reimMag' small helper function to calculate magnitude for PV and TPV
  reimMag :: Double -> Double -> Double -> Double -> Double
  reimMag v0 v1 v2 v3 =
    let sumsqs = v1^2 + v2^2 + v3^2
        x = abs v0 * sqrt sumsqs
    in sqrt (v0^2 + sumsqs + 2*x)
  
  -- |Cl(3,0) has a Fractional instance
  instance Fractional Cl3 where
    -- |Some of the sub algebras are division algebras but APS is not a division algebra
    recip (R a0) = R (recip a0)   -- R is a division algebra
    recip cliff = 
      let (R mag) = abs cliff
          recipsqmag = recip mag^2
          negrecipsqmag = negate recipsqmag
          recipmag2 = recip.toR $ cliff * bar cliff
          go_recip (V3 a1 a2 a3) = V3 (recipsqmag * a1) (recipsqmag * a2) (recipsqmag * a3)
          go_recip (BV a23 a31 a12) = BV (negrecipsqmag * a23) (negrecipsqmag * a31) (negrecipsqmag * a12)
          go_recip (I a123) = I (negrecipsqmag * a123)
          go_recip (H a0 a23 a31 a12) = H (recipsqmag * a0) (negrecipsqmag * a23) (negrecipsqmag * a31) (negrecipsqmag * a12)  -- H is a division algebra
          go_recip (C a0 a123) = C (recipsqmag * a0) (negrecipsqmag * a123)   -- C is a division algebra
          go_recip (ODD a1 a2 a3 a123) = ODD (recipsqmag * a1) (recipsqmag * a2) (recipsqmag * a3) (negrecipsqmag * a123)
          go_recip pv@PV{} = recipmag2 * bar pv
          go_recip tpv@TPV{} = recipmag2 * bar tpv
          go_recip cliffor = reduce $ spectraldcmp recip recip' cliffor
      in go_recip cliff
  
  
    -- |'fromRational'
    fromRational rat = R (fromRational rat)
  
  
  -- |Cl(3,0) has a "Floating" instance.
  instance Floating Cl3 where
    pi = R pi
  
    --
    exp (R a0) = R (exp a0)
    exp (I a123) = C (cos a123) (sin a123)
    exp (C a0 a123) =
      let expa0 = exp a0
      in C (expa0 * cos a123) (expa0 * sin a123)
    exp cliffor = spectraldcmp exp exp' cliffor
  
  
  
    --
    log (R a0)
      | a0 >= 0 = R (log a0)
      | a0 == (-1) = I pi
      | otherwise = C (log.negate $ a0) pi
    log (I a123)
      | a123 == 1 = I (pi/2)
      | a123 == (-1) = I (-pi/2)
      | otherwise = C (log.abs $ a123) (signum a123 * (pi/2))
    log (C a0 a123) = C (log (a0^2 + a123^2) / 2) (atan2 a123 a0)
    log cliffor = spectraldcmp log log' cliffor
  
  
    --
    sqrt (R a0)
      | a0 >= 0 = R (sqrt a0)
      | otherwise = I (sqrt.negate $ a0)
    sqrt (I a123)
      | a123 == 0 = R 0
      | otherwise =
          let sqrtr = sqrt.abs $ a123
              phiby2 = signum a123 * (pi/4) -- evaluated: atan2 a123 0 / 2
          in C (sqrtr * cos phiby2) (sqrtr * sin phiby2)
    sqrt (C a0 a123) =
      let sqrtr = sqrt.sqrt $ a0^2 + a123^2
          phiby2 = atan2 a123 a0 / 2
      in C (sqrtr * cos phiby2) (sqrtr * sin phiby2)
    sqrt cliffor = spectraldcmp sqrt sqrt' cliffor
  
    --
    sin (R a0) = R (sin a0)
    sin (I a123)
      | a123 == 0 = R 0
      | otherwise = I (sinh a123)
    sin (C a0 a123) = C (sin a0 * cosh a123) (cos a0 * sinh a123)
    sin cliffor = spectraldcmp sin sin' cliffor
  
    --
    cos (R a0) = R (cos a0)
    cos (I a123) = R (cosh a123)
    cos (C a0 a123) = C (cos a0 * cosh a123) (negate $ sin a0 * sinh a123)
    cos cliffor = spectraldcmp cos cos' cliffor
  
    --
    tan (R a0) = R (tan a0)
    tan (I a123)
      | a123 == 0 = R 0
      | otherwise = I (tanh a123)
    tan (reduce -> C a0 a123) =
      let
        m = x2^2 + y2^2
        x1 = sinx*coshy
        y1 = cosx*sinhy
        x2 = cosx*coshy
        y2 = negate $ sinx*sinhy
        sinx  = sin a0
        cosx  = cos a0
        sinhy = sinh a123
        coshy = cosh a123
      in C ((x1*x2 + y1*y2)/m) ((x2*y1 - x1*y2)/m)
    tan cliffor = spectraldcmp tan tan' cliffor
  
  
  
    --
    asin (R a0)
        -- asin (R a0) = I (-1) * log (I 1 * R a0 + sqrt (1 - (R a0)^2))
        -- I (-1) * log (I a0 + sqrt (R 1 - (R a0)^2))
        -- I (-1) * log (I a0 + sqrt (R (1 - a0^2)))
        -- I (-1) * log (I a0 + (I (sqrt.negate $ 1 - a0^2)))
        -- I (-1) * log (I a0 + (sqrt.negate $ 1 - a0^2))
        -- Def ==> log (I a123) = C (log.abs $ a123) (signum a123 * (pi/2))
        -- I (-1) * C (log.abs $ (a0 + (sqrt.negate $ 1 - a0^2))) (signum (a0 + (sqrt.negate $ 1 - a0^2)) * (pi/2))
        -- C (signum (a0 + (sqrt.negate $ 1 - a0^2)) * (pi/2)) (negate.log.abs $ (a0 + (sqrt.negate $ 1 - a0^2)))
      | a0 > 1 = C (pi/2) (negate.log $ (a0 + sqrt (a0^2 - 1)))
        -- I (-1) * log (I a0 + R (sqrt $ 1 - a0^2))
        -- I (-1) * log (C (sqrt $ 1 - a0^2) a0)
        -- Def ==> log (C a0 a123) = C (log.sqrt $ a0^2 + a123^2) (atan2 a123 a0)
        -- I (-1) * C (log.sqrt $ (sqrt $ 1 - a0^2)^2 + a0^2) (atan2 a0 (sqrt $ 1 - a0^2))
        -- C (atan2 a0 (sqrt $ 1 - a0^2)) (negate.log.sqrt $ (sqrt $ 1 - a0^2)^2 + a0^2)
        -- C (atan(a0/(sqrt $ 1 - a0^2))) (negate.log.sqrt $ 1)
        -- Apply sqrt 1 == 1, Apply log 1 == 0, reduce
        -- R (atan(a0/(sqrt $ 1 - a0^2)))
        -- Identity: tan(asin x) = x / (sqrt (1 - x^2))
        -- R (asin a0)
      | a0 >= (-1) = R (asin a0)
        -- I (-1) * log (I a0 + sqrt (R (1 - a0^2)))
        -- I (-1) * log (I (a0 + (sqrt.negate $ 1 - a0^2)))
        -- Def ==> log (I a123) = C (log.abs $ a123) (signum a123 * (pi/2))
        -- I (-1) * C (log.abs $ (a0 + (sqrt.negate $ 1 - a0^2))) (signum (a0 + (sqrt.negate $ 1 - a0^2)) * (pi/2))
        -- C (signum (a0 + (sqrt.negate $ 1 - a0^2)) * (pi/2)) (negate.log.abs $ (a0 + (sqrt.negate $ 1 - a0^2)))
        -- For the negative branch signum is -1
        -- C (-pi/2) (negate.log.abs $ (a0 + (sqrt $ a0^2 - 1)))
      | otherwise = C (-pi/2) (negate.log.abs $ (a0 + sqrt (a0^2 - 1)))
        --
        -- For I:
        -- I (-1) * log (I (1) * I a123 + sqrt (R 1 - (I a123)^2))
        -- I (-1) * log (R (-a123) + sqrt (R 1 - (I a123)^2))
        -- I (-1) * log (R (-a123) + sqrt (R 1 - R (-a123^2)))
        -- I (-1) * log (R (-a123) + sqrt (R (1 + a123^2)))
        -- I (-1) * log (R (-a123) + R (sqrt $ 1 + a123^2))
        -- I (-1) * log (R ((sqrt $ 1 + a123^2) - a123))
        -- ((sqrt $ 1 + a123^2) - a123)) is always positive
        -- Def ==> log (R a0) | a0 >= 0 = R (log a0)
        -- I (-1) * (R (log $ (sqrt $ 1 + a123^2) - a123))
        -- I (negate.log $ (sqrt $ 1 + a123^2) - a123)
        -- I (negate.log $ (sqrt $ 1 + a123^2) - a123)
        -- because ((sqrt $ 1 + a123^2) - a123)) is always positive: negate.log == log.recip
        -- I (log.recip $ (sqrt $ 1 + a123^2) - a123)
        -- recip $ (sqrt $ 1 + a123^2) - a123) == (sqrt $ 1 + a123^2) + a123)
        -- I (log $ (sqrt $ 1 + a123^2) + a123)
        -- I (asinh a123)
    asin (I a123)
      | a123 == 0 = R 0
      | otherwise = I (asinh a123)
      --
    asin (C a0 a123) =
        -- For C:
        -- I (-1) * log (I 1 * C a0 a123 + sqrt (R 1 - (C a0 a123)^2))
        -- I (-1) * log (C (-a123) a0 + sqrt (R 1 - (C a0 a123)^2))
        -- I (-1) * log (C (-a123) a0 + sqrt (C (1 - a0^2 + a123^2) (-2*a0*a123)))
        -- Def ==> sqrt (C a0 a123) = C ((sqrt.sqrt $ a0^2 + a123^2) * cos (atan2 a123 a0 / 2)) ((sqrt.sqrt $ a0^2 + a123^2) * sin (atan2 a123 a0 / 2))
        -- I (-1) * log (C (-a123) a0 + C ((sqrt.sqrt $ (1 - a0^2 + a123^2)^2 + (-2*a0*a123)^2) * cos (atan2 (-2*a0*a123) (1 - a0^2 + a123^2) / 2)) ((sqrt.sqrt $ (1 - a0^2 + a123^2)^2 + (-2*a0*a123)^2) * sin (atan2 (-2*a0*a123) (1 - a0^2 + a123^2) / 2)))
        -- I (-1) * log (C (((sqrt.sqrt $ (1 - a0^2 + a123^2)^2 + (-2*a0*a123)^2) * cos (atan2 (-2*a0*a123) (1 - a0^2 + a123^2) / 2)) - a123) (((sqrt.sqrt $ (1 - a0^2 + a123^2)^2 + (-2*a0*a123)^2) * sin (atan2 (-2*a0*a123) (1 - a0^2 + a123^2) / 2)) + a0))
        -- Def ==> log (C a0 a123) = C (log.sqrt $ a0^2 + a123^2) (atan2 a123 a0)
        -- C (atan2 (((sqrt.sqrt $ (1 - a0^2 + a123^2)^2 + (-2*a0*a123)^2) * sin (atan2 (-2*a0*a123) (1 - a0^2 + a123^2) / 2)) + a0)
        --          (((sqrt.sqrt $ (1 - a0^2 + a123^2)^2 + (-2*a0*a123)^2) * cos (atan2 (-2*a0*a123) (1 - a0^2 + a123^2) / 2)) - a123))
        --   (negate.log.sqrt $ (((sqrt.sqrt $ (1 - a0^2 + a123^2)^2 + (-2*a0*a123)^2) * cos (atan2 (-2*a0*a123) (1 - a0^2 + a123^2) / 2)) - a123)^2 +
        --                      (((sqrt.sqrt $ (1 - a0^2 + a123^2)^2 + (-2*a0*a123)^2) * sin (atan2 (-2*a0*a123) (1 - a0^2 + a123^2) / 2)) + a0)^2)
        -- Collect like terms:
      let theta = atan2 (-2*a0*a123) (1 - a0^2 + a123^2)
          rho = sqrt.sqrt $ (1 - a0^2 + a123^2)^2 + (-2*a0*a123)^2
          b0 = rho * cos (theta/2) - a123
          b123 = rho * sin (theta/2) + a0
      in C (atan2 b123 b0) (log (b0^2 + b123^2) / (-2))
      --
    asin cliffor = spectraldcmp asin asin' cliffor
  
    --
    acos (R a0)
        -- acos x == (pi/2) - asin x so just subistute
        -- For R a0 > 1:
        -- R (pi/2) - C (pi/2) (negate.log $ (a0 + (sqrt $ a0^2 - 1)))
        -- C 0 (negate.negate.log $ (a0 + (sqrt $ a0^2 - 1)))
        -- I (log $ (a0 + (sqrt $ a0^2 - 1)))
      | a0 > 1 = I (log (a0 + sqrt (a0^2 - 1)))
        -- For R a0 > (-1)
        -- R (pi/2) - R (asin a0) == R (acos a0)
      | a0 >= (-1) = R (acos a0)
        -- For R otherwise:
        -- R (pi/2) - C (-pi/2) (negate.log.abs $ (a0 + (sqrt $ a0^2 - 1)))
        -- C pi (negate.negate.log.abs $ (a0 + (sqrt $ a0^2 - 1)))
        -- C pi (log.abs $ (a0 + (sqrt $ a0^2 - 1)))
      | otherwise = C pi (log.abs $ (a0 + sqrt (a0^2 - 1)))
        --
        -- For I:
        -- asin (I a123)  = I (asinh a123) -- so,
        -- acos x == R (pi/2) - I (asinh a123)
        -- C (pi/2) (negate $ asinh a123)
    acos (I a123)
      | a123 == 0 = R (pi/2)
      | otherwise = C (pi/2) (negate $ asinh a123)
    --
    acos (C a0 a123) =
        -- For C:
        -- asin (C a0 a123) = C (atan2 (((sqrt.sqrt $ (1 - a0^2 + a123^2)^2 + (-2*a0*a123)^2) * sin (atan2 (-2*a0*a123) (1 - a0^2 + a123^2) / 2)) + a0) (((sqrt.sqrt $ (1 - a0^2 + a123^2)^2 + (-2*a0*a123)^2) * cos (atan2 (-2*a0*a123) (1 - a0^2 + a123^2) / 2)) - a123)) (negate.log.sqrt $ (((sqrt.sqrt $ (1 - a0^2 + a123^2)^2 + (-2*a0*a123)^2) * cos (atan2 (-2*a0*a123) (1 - a0^2 + a123^2) / 2)) - a123)^2 + (((sqrt.sqrt $ (1 - a0^2 + a123^2)^2 + (-2*a0*a123)^2) * sin (atan2 (-2*a0*a123) (1 - a0^2 + a123^2) / 2)) + a0)^2)
        -- acos x == (pi/2) - asin x so just subistute
        -- R (pi/2) - C (atan2 (((sqrt.sqrt $ (1 - a0^2 + a123^2)^2 + (-2*a0*a123)^2) * sin (atan2 (-2*a0*a123) (1 - a0^2 + a123^2) / 2)) + a0) (((sqrt.sqrt $ (1 - a0^2 + a123^2)^2 + (-2*a0*a123)^2) * cos (atan2 (-2*a0*a123) (1 - a0^2 + a123^2) / 2)) - a123)) (negate.log.sqrt $ (((sqrt.sqrt $ (1 - a0^2 + a123^2)^2 + (-2*a0*a123)^2) * cos (atan2 (-2*a0*a123) (1 - a0^2 + a123^2) / 2)) - a123)^2 + (((sqrt.sqrt $ (1 - a0^2 + a123^2)^2 + (-2*a0*a123)^2) * sin (atan2 (-2*a0*a123) (1 - a0^2 + a123^2) / 2)) + a0)^2)
      let theta = atan2 (-2*a0*a123) (1 - a0^2 + a123^2)
          rho = sqrt.sqrt $ (1 - a0^2 + a123^2)^2 + (-2*a0*a123)^2
          b0 = rho * cos (theta/2) - a123
          b123 = rho * sin (theta/2) + a0
      in C ((pi/2) - atan2 b123 b0) (log (b0^2 + b123^2) / 2)
      --
    acos cliffor = spectraldcmp acos acos' cliffor
  
    --
    atan (R a0) = R (atan a0)
    --
    atan (I a123)
        -- I (0.5) * (log (R 1 - (I 1 * I a123)) - log (R 1 + (I 1 * I a123)))
        -- I (0.5) * (log (R 1 - (R (-a123))) - log (R 1 + (R (-a123))))
        -- I (0.5) * ((log (R (1 + a123))) - log (R (1 - a123)))
        -- Def ==> C (log.negate $ a0) pi for negative a123
        -- I (0.5) * ((log (R (1 + a123))) - (C (log.negate $ (1 - a123)) pi))
        -- Def ==> R (log a0) for positive a123
        -- I (0.5) * ((R (log (1 + a123))) - (C (log.negate $ (1 - a123)) pi))
        -- I (0.5) * (C (log (1 + a123) - (log.negate $ (1 - a123))) (-pi))
        -- C (pi/2) ((log (1 + a123) - (log.negate $ (1 - a123)))/2)
      | a123 > 1 = C (pi/2) (0.5*(log (1 + a123) - log (a123 - 1)))
        -- I (0.5) * (log (R 1 - (I 1 * I a123)) - log (R 1 + (I 1 * I a123)))
        -- I (0.5) * (log (R 1 - (R (-a123))) - log (R 1 + (R (-a123))))
        -- I (0.5) * ((log (R (1 + a123))) - log (R (1 - a123)))
        -- I (0.5) * ((R (log (1 + a123))) - R (log (1 - a123)))
        -- I (0.5) * (R ((log (1 + a123)) - (log (1 - a123))))
        -- I (((log (1 + a123)) - (log (1 - a123)))/2)
        -- I (atanh a123)
      | a123 == 0 = R 0
      | a123 >= (-1) = I (atanh a123)
        -- I (0.5) * (log (R 1 - (I 1 * I a123)) - log (R 1 + (I 1 * I a123)))
        -- I (0.5) * (log (R 1 - (R (-a123))) - log (R 1 + (R (-a123))))
        -- I (0.5) * ((log (R (1 + a123))) - R (log (1 - a123)))
        -- C (-pi/2) (((log.negate $ (1 + a123)) - (log (1 - a123)))/2)
      | otherwise = C (-pi/2) (((log.negate $ (1 + a123)) - log (1 - a123))/2)
        --
        -- I (0.5) * (log (R 1 - (I 1 * C a0 a123)) - log (R 1 + (I 1 * C a0 a123)))
        -- I (0.5) * (log (C (1 + a123) (-a0)) - log (C (1 - a123) a0))
        -- Def ==> log (C a0 a123) = C (log.sqrt $ a0^2 + a123^2) (atan2 a123 a0)
        -- I (0.5) * ((C (log.sqrt $ (1 + a123)^2 + (-a0)^2) (atan2 (-a0) (1 + a123))) - (C (log.sqrt $ (1 - a123)^2 + a0^2) (atan2 a0 (1 - a123))))
        -- I (0.5) * C ((log.sqrt $ (1 + a123)^2 + (-a0)^2) - (log.sqrt $ (1 - a123)^2 + a0^2)) ((atan2 (-a0) (1 + a123)) - (atan2 a0 (1 - a123)))
        -- I (0.5) * C (0.5*((log $ (1 + a123)^2 + a0^2) - (log $ (1 - a123)^2 + a0^2))) ((atan2 (-a0) (1 + a123)) - (atan2 a0 (1 - a123)))
        -- C (((atan2 a0 (1 - a123)) + (atan2 a0 (1 + a123)))/2) (((log $ (1 + a123)^2 + a0^2) - (log $ (1 - a123)^2 + a0^2))/4)
    atan (C a0 a123) = C ((atan2 a0 (1 - a123) + atan2 a0 (1 + a123))/2)
                         ((log ((1 + a123)^2 + a0^2) - log ((1 - a123)^2 + a0^2))/4)
      --
    atan cliffor = spectraldcmp atan atan' cliffor
  
    --
    sinh (R a0) = R (sinh a0)
    sinh (I a123) = I (sin a123)
    sinh (C a0 a123) = C (cos a123 * sinh a0) (sin a123 * cosh a0)
    sinh cliffor = spectraldcmp sinh sinh' cliffor
  
    --
    cosh (R a0) = R (cosh a0)
    cosh (I a123) = R (cos a123)
    cosh (C a0 a123) = C (cos a123 * cosh a0) (sin a123 * sinh a0)
    cosh cliffor = spectraldcmp cosh cosh' cliffor
  
    --
    tanh (R a0) = R (tanh a0)
    tanh (I a123) = I (tan a123)
    tanh (reduce -> C a0 a123) =
      let
        m = x2^2 + y2^2
        x1 = cosy*sinhx
        y1 = siny*coshx
        x2 = cosy*coshx
        y2 = siny*sinhx
        siny  = sin a123
        cosy  = cos a123
        sinhx = sinh a0
        coshx = cosh a0
      in C ((x1*x2 + y1*y2)/m) ((x2*y1 - x1*y2)/m)
    tanh cliffor = spectraldcmp tanh tanh' cliffor
  
    --
    asinh (R a0) = R (asinh a0)
    --
    asinh (I a123)
        -- log (I a123 + sqrt (R (1 - a123^2)))
        -- 3 branches where between -1 and 1 it is just asin
        -- For a123 > 1:
        -- log (I a123 + I (sqrt.negate $ (1 - a123^2)))
        -- log (I (a123 + (sqrt (a123^2 - 1))))
        -- Def ==> log (I a123) = C (log.abs $ a123) (signum a123 * (pi/2))
        -- C (log.abs $ (a123 + (sqrt (a123^2 - 1)))) (signum (a123 + (sqrt (a123^2 - 1))) * (pi/2))
        -- a123 is positive so signum evaluates to 1
        -- C (log.abs $ (a123 + (sqrt (a123^2 - 1)))) (pi/2)
      | a123 > 1 = C (log.abs $ (a123 + sqrt (a123^2 - 1))) (pi/2)
        -- log (I a123 + sqrt (R (1 - a123^2)))
        -- log (I a123 + R (sqrt (1 - a123^2)))
        -- log (C (sqrt (1 - a123^2)) a123)
        -- Def ==> log (C a0 a123) = C (log.sqrt $ a0^2 + a123^2) (atan2 a123 a0)
        -- C (log.sqrt $ (sqrt (1 - a123^2))^2 + a123^2) (atan2 a123 (sqrt (1 - a123^2)))
        -- (sqrt (1 - a123^2))^2 + a123^2 == 1
        -- sqrt 1 == 1
        -- log 1 == 0
        -- I (atan2 a123 (sqrt (1 - a123^2)))
        -- I (atan (a123 / (sqrt (1 - a123^2))))
        -- Identity: tan(asin x) = x / (sqrt (1 - x^2))
        -- asin a123 = atan (a123 / (sqrt (1 - a123^2)))
      | a123 == 0 = R 0
      | a123 >= (-1) = I (asin a123)
        -- log (I a123 + sqrt (R (1 - a123^2)))
        -- For a123 < (-1):
        -- log (I a123 + I (sqrt.negate $ (1 - a123^2)))
        -- log (I (a123 + (sqrt (a123^2 - 1))))
        -- Def ==> log (I a123) = C (log.abs $ a123) (signum a123 * (pi/2))
        -- C (log.abs $ (a123 + (sqrt (a123^2 - 1)))) (signum (a123 + (sqrt (a123^2 - 1))) * (pi/2))
        -- for a123 lt (-1) signum evaluates to -1
      | otherwise = C (log.abs $ (a123 + sqrt (a123^2 - 1))) (-pi/2)
      --
    asinh (C a0 a123) =
        -- For C:
        -- log (C a0 a123 + sqrt (C (a0^2 - a123^2 +1) (2*a0*a123)))
        -- Def ==> sqrt (C a0 a123) = C ((sqrt.sqrt $ a0^2 + a123^2) * cos (atan2 a123 a0 / 2)) ((sqrt.sqrt $ a0^2 + a123^2) * sin (atan2 a123 a0 / 2))
        -- log (C a0 a123 + C ((sqrt.sqrt $ (a0^2 - a123^2 +1)^2 + (2*a0*a123)^2) * cos (atan2 (2*a0*a123) (a0^2 - a123^2 +1) / 2)) ((sqrt.sqrt $ (a0^2 - a123^2 +1)^2 + (2*a0*a123)^2) * sin (atan2 (2*a0*a123) (a0^2 - a123^2 +1) / 2)))
        -- log (C (a0 + ((sqrt.sqrt $ (a0^2 - a123^2 +1)^2 + (2*a0*a123)^2) * cos (atan2 (2*a0*a123) (a0^2 - a123^2 +1) / 2))) (a123 + ((sqrt.sqrt $ (a0^2 - a123^2 +1)^2 + (2*a0*a123)^2) * sin (atan2 (2*a0*a123) (a0^2 - a123^2 +1) / 2))))
        -- Def ==> log (C a0 a123) = C (log.sqrt $ a0^2 + a123^2) (atan2 a123 a0)
        -- C (log.sqrt $ (a0 + ((sqrt.sqrt $ (a0^2 - a123^2 +1)^2 + (2*a0*a123)^2) * cos (atan2 (2*a0*a123) (a0^2 - a123^2 +1) / 2)))^2 +
        --               (a123 + ((sqrt.sqrt $ (a0^2 - a123^2 +1)^2 + (2*a0*a123)^2) * sin (atan2 (2*a0*a123) (a0^2 - a123^2 +1) / 2)))^2)
        --   (atan2 (a123 + ((sqrt.sqrt $ (a0^2 - a123^2 +1)^2 + (2*a0*a123)^2) * sin (atan2 (2*a0*a123) (a0^2 - a123^2 +1) / 2)))
        --          (a0 + ((sqrt.sqrt $ (a0^2 - a123^2 +1)^2 + (2*a0*a123)^2) * cos (atan2 (2*a0*a123) (a0^2 - a123^2 +1) / 2))))
        -- Collect like terms:
      let theta = atan2 (2*a0*a123) (a0^2 - a123^2 +1)
          rho = sqrt.sqrt $ (a0^2 - a123^2 +1)^2 + (2*a0*a123)^2
          b0 = a0 + rho * cos (theta/2)
          b123 = a123 + rho * sin (theta/2)
      in C (log (b0^2 + b123^2) / 2) (atan2 b123 b0)
      --
    asinh cliffor = spectraldcmp asinh asinh' cliffor
  
    --
    acosh (R a0)
      -- log (R a0 + sqrt(R (a0+1)) * sqrt(R (a0-1)))
      | a0 >= 1 = R (acosh a0)
        -- log (R a0 + sqrt(R (a0+1)) * sqrt(R (a0-1)))
        -- log (R a0 + R (sqrt $ a0+1) * R (sqrt $ a0-1))
        -- log (R a0 + R ((sqrt $ a0+1) * (sqrt $ a0-1)))
        -- log (R (a0 + (sqrt $ a0+1) * (sqrt $ a0-1)))
        -- R (log $ a0 + (sqrt $ a0+1) * (sqrt $ a0-1))
        -- R (acosh a0)
        -- Strangely ghc substitutes 'acosh a0' with something like:
        -- R (log $ a0 + (a0 + 1 ) * (sqrt $ (a0 - 1)/(a0 + 1)))
      | a0 >= (-1) = I (atan2 (sqrt $ 1-a0^2) a0) -- This is I because of cancelation of the real component
        -- log (R a0 + sqrt(R (a0+1)) * sqrt(R (a0-1)))
        -- log (R a0 + R (sqrt $ a0+1) * I (sqrt.negate $ a0-1))
        -- log (R a0 + I ((sqrt $ a0+1) * (sqrt.negate $ a0-1)))
        -- log (R a0 + I ((sqrt $ a0+1) * (sqrt $ 1-a0)))
        -- log $ C (a0) ((sqrt $ a0+1) * (sqrt $ 1-a0))
        -- Def log ==> log (C b0 b123) = C (log.sqrt $ b0^2 + b123^2) (atan2 b123 b0)
        -- let b0 = a0
        --     b123 = (sqrt $ a0+1) * (sqrt $ 1-a0) = sqrt $ 1-a0^2
        -- in C (log.sqrt $ b0^2 + b123^2) (atan2 b123 b0)
        -- b123^2 = 1-a0^2
        -- C (log.sqrt $ a0^2 + 1-a0^2) (atan2 (sqrt $ 1-a0^2) a0)
        -- C (log.sqrt $ 1) (atan2 (sqrt $ 1-a0^2) a0)
        -- C 0 (atan2 (sqrt $ 1-a0^2) a0)
        -- I (atan2 (sqrt $ 1-a0^2) a0)
      | otherwise = C (acosh.negate $ a0) pi
        -- log (R a0 + sqrt(R (a0+1)) * sqrt(R (a0-1)))
        -- log (R a0 + I (sqrt.negate $ a0+1) * I (sqrt.negate $ a0-1))
        -- Def ==> (I a123) * (I b123) = R (negate $ a123*b123)
        -- log (R a0 + R (negate $ (sqrt.negate $ a0+1) * (sqrt.negate $ a0-1))
        -- log (R (a0 + (negate $ (sqrt.negate $ a0+1) * (sqrt.negate $ a0-1))))
        -- C (log.negate $ (a0 + (negate $ (sqrt.negate $ a0+1) * (sqrt.negate $ a0-1)))) pi
        -- C (log $ (negate a0 + ((sqrt $ (negate a0)+1) * (sqrt $ (negate a0)-1)))) pi
        -- C (acosh (negate a0)) pi
        --
    acosh (I a123)
        -- log (I a123 + sqrt(C 1 a123) * sqrt(C (-1) a123))
        -- Def ==> sqrt (C a0 a123) =
        --   C ((sqrt.sqrt $ a0^2 + a123^2) * cos (atan2 a123 a0 / 2))
        --      ((sqrt.sqrt $ a0^2 + a123^2) * sin (atan2 a123 a0 / 2))
        -- log (I a123 +
        --      C ((sqrt.sqrt $ 1 + a123^2) * cos (atan2 a123 1 / 2))
        --        ((sqrt.sqrt $ 1 + a123^2) * sin (atan2 a123 1 / 2)) *
        --      C ((sqrt.sqrt $ 1 + a123^2) * cos (atan2 a123 (-1) / 2))
        --        ((sqrt.sqrt $ 1 + a123^2) * sin (atan2 a123 (-1) / 2)) )
        -- Factor out "(sqrt.sqrt $ 1 + a123^2)*"
        -- log (I a123 + R (sqrt.sqrt $ 1 + a123^2) *
        --               C (cos (atan2 a123 1 / 2)) (sin (atan2 a123 1 / 2)) *
        --               R (sqrt.sqrt $ 1 + a123^2) *
        --               C (cos (atan2 a123 (-1) / 2)) (sin (atan2 a123 (-1) / 2)))
        -- Collect both R's and simplify
        -- log (I a123 + (R (sqrt $ 1 + a123^2)) *
        --                C (cos (atan2 a123 1 / 2)) (sin (atan2 a123 1 / 2)) *
        --                C (cos (atan2 a123 (-1) / 2)) (sin (atan2 a123 (-1) / 2)))
        -- Def ==> (C a0 a123) * (C b0 b123) = C (a0*b0 - a123*b123) (a0*b123 + a123*b0)
        -- log (I a123 + R (sqrt $ 1 + a123^2) *
        --               C ((cos (atan2 a123 1 / 2))*(cos (atan2 a123 (-1) / 2)) - (sin (atan2 a123 1 / 2))*(sin (atan2 a123 (-1) / 2)))
        --                 ((cos (atan2 a123 1 / 2))*(sin (atan2 a123 (-1) / 2)) + (sin (atan2 a123 1 / 2))*(cos (atan2 a123 (-1) / 2))) )
        --
        -- Solution now branches for positive and negative a123
        --
        -- For a123 > 0 Substitute (cos (atan2 a123 (-1) / 2)) == (sin (atan2 a123 1 / 2)) AND
        --                         (sin (atan2 a123 (-1) / 2)) == (cos (atan2 a123 1 / 2)) AND
        --                         atan2 a123 1 == atan a123
        -- log (I a123 + R (sqrt $ 1 + a123^2) *
        --               C ((cos (atan a123 / 2))*(sin (atan a123 / 2)) - (sin (atan a123 / 2))*(cos (atan a123 / 2)))
        --                 ((cos (atan a123 / 2))*(cos (atan a123 / 2)) + (sin (atan a123 / 2))*(sin (atan a123 / 2))) )
        -- sin^2 + cos^2 == 1 AND cos*sin - sin*cos == 0 AND Reduce C 0 1 to I 1 AND apply (*) AND apply (+)
        -- log (I (a123 + sqrt (1 + a123^2)))
        -- Def ==> log (I a123) = C (log.abs $ a123) (signum a123 * (pi/2))
        -- C (log.abs $ (a123 + sqrt (1 + a123^2))) (signum (a123 + sqrt (1 + a123^2)) * (pi/2))
        -- With a123 positive Apply signum:
        -- C (log.abs $ (a123 + sqrt (1 + a123^2))) (pi/2)
      | a123 > 0 = C (log.abs $ (a123 + sqrt (1 + a123^2))) (pi/2)
        -- With a123 == 0:
        -- reduce C 0 (pi/2)
        -- I (pi/2)
      | a123 == 0 = I (pi/2)
        -- For a123 < 0 Substitute (cos (atan2 a123 (-1) / 2)) == (negate.sin $ (atan2 a123 1 / 2)) AND
        --                         (sin (atan2 a123 (-1) / 2)) == (negate.cos $ (atan2 a123 1 / 2)) AND
        --                         atan2 a123 1 == atan a123
        -- log (I a123 + R (sqrt $ 1 + a123^2) *
        --               C ((cos (atan2 a123 1 / 2))*(negate.sin $ (atan2 a123 1 / 2)) - (sin (atan2 a123 1 / 2))*(negate.cos $ (atan2 a123 1 / 2)))
        --                 ((cos (atan2 a123 1 / 2))*(negate.cos $ (atan2 a123 1 / 2)) + (sin (atan2 a123 1 / 2))*(negate.sin $ (atan2 a123 1 / 2))) )
        -- Factor negate out AND sin^2 + cos^2 == 1 AND cos*sin - sin*cos == 0 AND Reduce C 0 (-1) to I (-1) AND apply (*) AND apply (+)
        -- log (I (a123 - sqrt (1 + a123^2)))
        -- Def ==> log (I a123) = C (log.abs $ a123) (signum a123 * (pi/2))
        -- C (log.abs $ (a123 - sqrt (1 + a123^2))) (signum (a123 - sqrt (1 + a123^2)) * (pi/2))
        -- With a123 negateive Apply signum:
        -- C (log.abs $ (a123 - sqrt (1 + a123^2))) (-pi/2)
      | otherwise = C (log.abs $ (a123 - sqrt (1 + a123^2))) (-pi/2)
      --
    acosh (C a0 a123) =
        -- log (C a0 a123 + sqrt(C (a0+1) a123) * sqrt(C (a0-1) a123))
        -- Def ==> sqrt (C a0 a123) =
        --   C ((sqrt.sqrt $ a0^2 + a123^2) * cos (atan2 a123 a0 / 2))
        --      ((sqrt.sqrt $ a0^2 + a123^2) * sin (atan2 a123 a0 / 2))
        -- log (C a0 a123 +
        --      C ((sqrt.sqrt $ (a0+1)^2 + a123^2) * cos (atan2 a123 (a0+1) / 2))
        --        ((sqrt.sqrt $ (a0+1)^2 + a123^2) * sin (atan2 a123 (a0+1) / 2)) *
        --      C ((sqrt.sqrt $ (a0-1)^2 + a123^2) * cos (atan2 a123 (a0-1) / 2))
        --        ((sqrt.sqrt $ (a0-1)^2 + a123^2) * sin (atan2 a123 (a0-1) / 2)) )
        -- Factor out the scalar in both Complex numbers
        -- log (C a0 a123 +
        --      R (sqrt.sqrt $ (a0+1)^2 + a123^2) *
        --      C (cos (atan2 a123 (a0+1) / 2)) (sin (atan2 a123 (a0+1) / 2)) *
        --      R (sqrt.sqrt $ (a0-1)^2 + a123^2) *
        --      C (cos (atan2 a123 (a0-1) / 2)) (sin (atan2 a123 (a0-1) / 2)) )
        -- Combine the R terms
        -- log (C a0 a123 +
        --      R (sqrt.sqrt $ ((a0+1)^2 + a123^2) * ((a0-1)^2 + a123^2)) *
        --      C (cos (atan2 a123 (a0+1) / 2)) (sin (atan2 a123 (a0+1) / 2)) *
        --      C (cos (atan2 a123 (a0-1) / 2)) (sin (atan2 a123 (a0-1) / 2)) )
        -- Def ==> (C a0 a123) * (C b0 b123) = C (a0*b0 - a123*b123)
        --                                       (a0*b123 + a123*b0)
        -- log (C a0 a123 +
        --      R (sqrt.sqrt $ ((a0+1)^2 + a123^2) * ((a0-1)^2 + a123^2)) *
        --      C (((cos (atan2 a123 (a0+1) / 2))*(cos (atan2 a123 (a0-1) / 2))) - ((sin (atan2 a123 (a0+1) / 2))*(sin (atan2 a123 (a0-1) / 2))))
        --        (((cos (atan2 a123 (a0+1) / 2))*(sin (atan2 a123 (a0-1) / 2))) + ((sin (atan2 a123 (a0+1) / 2))*(cos (atan2 a123 (a0-1) / 2)))) )
        -- =
        -- log (C a0 a123 +
        --      R (sqrt.sqrt $ ((a0+1)^2 + a123^2) * ((a0-1)^2 + a123^2)) *
        --      C (cos(0.5*(atan2 a123 (a0+1) + atan2 a123 (a0-1))))
        --        (sin(0.5*(atan2 a123 (a0-1) + atan2 a123 (a0+1)))) )
        -- Apply (*)
        -- log (C a0 a123 +
        --      C ((sqrt.sqrt $ ((a0+1)^2 + a123^2) * ((a0-1)^2 + a123^2)) *(cos(0.5*(atan2 a123 (a0+1) + atan2 a123 (a0-1)))))
        --        ((sqrt.sqrt $ ((a0+1)^2 + a123^2) * ((a0-1)^2 + a123^2)) *(sin(0.5*(atan2 a123 (a0-1) + atan2 a123 (a0+1))))) )
        -- Apply (+)
        -- log (C (a0 + (sqrt.sqrt $ ((a0+1)^2 + a123^2) * ((a0-1)^2 + a123^2)) * ((cos(0.5*(atan2 a123 (a0+1) + atan2 a123 (a0-1))))))
        --        (a123 + (sqrt.sqrt $ ((a0+1)^2 + a123^2) * ((a0-1)^2 + a123^2)) * ((sin(0.5*(atan2 a123 (a0-1) + atan2 a123 (a0+1)))))) )
        -- Def ==>  log (C a0 a123) = C (log.sqrt $ a0^2 + a123^2) (atan2 a123 a0)
        -- = C (log.sqrt $ (a0 + (sqrt.sqrt $ ((a0+1)^2 + a123^2) * ((a0-1)^2 + a123^2)) * ((cos(0.5*(atan2 a123 (a0+1) + atan2 a123 (a0-1))))))^2 + (a123 + (sqrt.sqrt $ ((a0+1)^2 + a123^2) * ((a0-1)^2 + a123^2)) * ((sin(0.5*(atan2 a123 (a0-1) + atan2 a123 (a0+1))))))^2) 
        --     (atan2 (a123 + (sqrt.sqrt $ ((a0+1)^2 + a123^2) * ((a0-1)^2 + a123^2)) * ((sin(0.5*(atan2 a123 (a0-1) + atan2 a123 (a0+1)))))) (a0 + (sqrt.sqrt $ ((a0+1)^2 + a123^2) * ((a0-1)^2 + a123^2)) * ((cos(0.5*(atan2 a123 (a0+1) + atan2 a123 (a0-1)))))))
        -- Collect like terms:
      let theta = atan2 a123 (a0+1) + atan2 a123 (a0-1)
          rho = sqrt.sqrt $ ((a0+1)^2 + a123^2) * ((a0-1)^2 + a123^2)
          b0 = a0 + rho * cos(theta/2)
          b123 = a123 + rho * sin(theta/2)
      in C (log (b0^2 + b123^2) / 2) (atan2 b123 b0)
      --
    acosh cliffor = spectraldcmp acosh acosh' cliffor
  
    --
    atanh (R a0)
        -- = 0.5*log (R (1+a0)) - 0.5*log (R (1-a0))
        -- = (R ((0.5*).log $ 1+a0)) - (C ((0.5*).log.negate $ 1-a0) (pi/2))
        -- = C (((0.5*).log $ 1+a0) - ((0.5*).log.negate $ 1-a0)) (-pi/2)
        -- = C (0.5*((log $ 1+a0) - (log $ a0-1))) (-pi/2)
      | a0 > 1 = C ((log (1+a0) - log (a0-1))/2) (-pi/2)
        -- = 0.5 * (log (R (1+a0)) - log (R (1-a0)))
        -- = 0.5*(R (log $ 1+a0) - R (log $ 1-a0))
        -- = R (0.5*(log $ 1+a0) - 0.5*(log $ 1-a0))
        -- = R (atanh a0)
      | a0 >= (-1) = R (atanh a0)
        -- = 0.5 * (log (R (1+a0)) - log (R (1-a0)))
        -- = (C ((0.5*).log.negate $ 1+a0) (pi/2)) - (R ((0.5*).log $ 1-a0))
        -- = C (((0.5*).log.negate $ 1+a0) - ((0.5*).log $ 1-a0)) (pi/2)
        -- = C (0.5*((log.negate $ 1+a0) - (log $ 1-a0))) (pi/2)
      | otherwise = C (((log.negate $ 1+a0) - log (1-a0))/2) (pi/2)
      --
      -- For I:
      -- = 0.5 * (log (C 1 a123) - log (C 1 (-a123)))
      -- = I (atan a123)
    atanh (I a123)
      | a123 == 0 = R 0
      | otherwise = I (atan a123)
      -- = 0.5 * (log (C (1+a0) a123) - log (C (1-a0) (-a123)))
      -- Def log ==> log (C a0 a123) = C (log.sqrt $ a0^2 + a123^2) (atan2 a123 a0)
      -- log (C (1+a0) a123) = C (log.sqrt $ (1+a0)^2 + a123^2) (atan2 a123 (1+a0))
      -- log (C (1-a0) (-a123)) = C (log.sqrt $ (1-a0)^2 + (-a123)^2) (atan2 (-a123) (1-a0))
      -- = C (((0.5*).log.sqrt $ (1+a0)^2 + a123^2) - ((0.5*).log.sqrt $ (1-a0)^2 + a123^2)) (0.5*((atan2 a123 (1+a0)) - (atan2 (-a123) (1-a0))))
      -- C (((log $ (1+a0)^2 + a123^2) - (log $ (1-a0)^2 + a123^2))/4) (((atan2 a123 (1-a0)) + (atan2 a123 (1+a0)))/2)
    atanh (C a0 a123) = C ((log ((1+a0)^2 + a123^2) - log ((1-a0)^2 + a123^2))/4) ((atan2 a123 (1-a0) + atan2 a123 (1+a0))/2)
    --
    atanh cliffor = spectraldcmp atanh atanh' cliffor
  
  
  
  -- |'lsv' the littlest singular value. Useful for testing for invertability.
  lsv :: Cl3 -> Cl3
  lsv (R a0) = R (abs a0) -- absolute value of a real number
  lsv (V3 a1 a2 a3) = R (sqrt (a1^2 + a2^2 + a3^2)) -- magnitude of a vector
  lsv (BV a23 a31 a12) = R (sqrt (a23^2 + a31^2 + a12^2)) -- magnitude of a bivector
  lsv (I a123) = R (abs a123)
  lsv (PV a0 a1 a2 a3) = R (loDisc a0 a1 a2 a3)
  lsv (TPV a23 a31 a12 a123) = R (loDisc a123 a23 a31 a12)
  lsv (H a0 a23 a31 a12) = R (sqrt (a0^2 + a23^2 + a31^2 + a12^2))
  lsv (C a0 a123) = R (sqrt (a0^2 + a123^2)) -- magnitude of a complex number
  lsv (BPV a1 a2 a3 a23 a31 a12) =
    let x = negate.sqrt $ (a1*a31 - a2*a23)^2 + (a1*a12 - a3*a23)^2 + (a2*a12 - a3*a31)^2 -- core was duplicating this computation added let to hopefully reduce the duplication
        y = a1^2 + a23^2 + a2^2 + x + a31^2 + a3^2 + a12^2 + x -- attempted to balance out the sum of several positives with a negitive before the next sum of positives and negitive
    in if y <= tol' -- gaurd for numerical errors, y could be negative with large enough biparavectors
       then R 0
       else R (sqrt y)
  lsv (ODD a1 a2 a3 a123) = R (sqrt (a1^2 + a2^2 + a3^2 + a123^2))
  lsv (APS a0 a1 a2 a3 a23 a31 a12 a123) =
    let x = negate.sqrt $ (a0*a1 + a123*a23)^2 + (a0*a2 + a123*a31)^2 + (a0*a3 + a123*a12)^2 +
                          (a2*a12 - a3*a31)^2 + (a3*a23 - a1*a12)^2 + (a1*a31 - a2*a23)^2 -- core was duplicating this computation added let to hopefully reduce the duplication
        y = a0^2 + a1^2 + a2^2 + a3^2 + x + a23^2 + a31^2 + a12^2 + a123^2 + x -- attempted to balance out the sum of several positives with a negitive before the next sum of positives and negitive
    in if y <= tol' -- gaurd for numerical errors, y could be negative with large enough cliffors
       then R 0
       else R (sqrt y)
  
  
  -- | 'loDisc' The Lower Discriminant for Paravectors and Triparavectors, real and imagninary portions of APS
  loDisc :: Double -> Double -> Double -> Double -> Double
  loDisc v0 v1 v2 v3 =
    let sumsqs = v1^2 + v2^2 + v3^2
        x = negate $ abs v0 * sqrt sumsqs
        y = v0^2 + x + sumsqs + x
    in if y <= tol' -- gaurd for numerical errors, y could be negative with large enough paravectors
       then 0
       else sqrt y
  
  
  -- | 'spectraldcmp' the spectral decomposition of a function to calculate analytic functions of cliffors in Cl(3,0).
  -- This function requires the desired function's R, I, and C instances to be calculated and the function's derivative.
  -- If multiple functions are being composed, its best to pass the composition of the funcitons
  -- to this function and the derivative to this function.  Any function with a Taylor Series
  -- approximation should be able to be used.  A real, imaginary, and complex version of the function to be decomposed
  -- must be provided and spectraldcmp will handle the case for an arbitrary Cliffor.
  -- 
  -- It may be possible to add, in the future, a RULES pragma like:
  --
  -- > "spectral decomposition function composition"
  -- > forall f f' g g' cliff.
  -- > spectraldcmp f f' (spectraldcmp g g' cliff) = spectraldcmp (f.g) (f'.g') cliff
  -- 
  -- 
  spectraldcmp :: (Cl3 -> Cl3) -> (Cl3 -> Cl3) -> Cl3 -> Cl3
  spectraldcmp fun fun' (reduce -> cliffor) = dcmp cliffor
    where
      dcmp r@R{} = fun r
      dcmp i@I{} = fun i
      dcmp c@C{} = fun c
      dcmp v@V3{} = spectraldcmpSpecial toR fun v -- spectprojR fun v
      dcmp pv@PV{} = spectraldcmpSpecial toR fun pv -- spectprojR fun pv
      dcmp bv@BV{} = spectraldcmpSpecial toI fun bv -- spectprojI fun bv
      dcmp tpv@TPV{} = spectraldcmpSpecial toI fun tpv -- spectprojI fun tpv
      dcmp h@H{} = spectraldcmpSpecial toC fun h -- spectprojC fun h
      dcmp od@ODD{} = spectraldcmpSpecial toC fun od -- spectprojC fun od
      dcmp cliff
        | hasNilpotent cliff = jordan toC fun fun' cliff  -- jordan normal form Cl3 style
        | isColinear cliff = spectraldcmpSpecial toC fun cliff -- spectprojC fun bpv
        | otherwise =                               -- transform it so it will be colinear
            let (BPV a1 a2 a3 a23 a31 a12) = toBPV cliff
                boost = boost2colinear a1 a2 a3 a23 a31 a12
            in boost * spectraldcmpSpecial toC fun (bar boost * cliff * boost) * bar boost -- v * spectprojC fun d * v_bar
  --
  
  
  -- | 'jordan' does a Cl(3,0) version of the decomposition into Jordan Normal Form and Matrix Function Calculation
  -- The intended use is for calculating functions for cliffors with vector parts simular to Nilpotent.
  -- It is a helper function for 'spectraldcmp'.  It is fortunate because eigen decomposition doesn't
  -- work with elements with nilpotent content, so it fills the gap.
  jordan :: (Cl3 -> Cl3) -> (Cl3 -> Cl3) -> (Cl3 -> Cl3) -> Cl3 -> Cl3
  jordan toSpecial fun fun' cliffor =
    let eigs = toSpecial cliffor
    in fun eigs + fun' eigs * toBPV cliffor
  
  -- | 'spectraldcmpSpecial' helper function for with specialization for real, imaginary, or complex eigenvalues.
  -- To specialize for Reals pass 'toR', to specialize for Imaginary pass 'toI', to specialize for Complex pass 'toC'
  spectraldcmpSpecial :: (Cl3 -> Cl3) -> (Cl3 -> Cl3) -> Cl3 -> Cl3
  spectraldcmpSpecial toSpecial function cliffor =
    let (p,p_bar,eig1,eig2) = projEigs toSpecial cliffor
    in function eig1 * p + function eig2 * p_bar
  
  
  
  -- | 'eigvals' calculates the eignenvalues of the cliffor.
  -- This is useful for determining if a cliffor is the pole
  -- of a function.
  eigvals :: Cl3 -> (Cl3,Cl3)
  eigvals (reduce -> cliffor) = eigv cliffor
    where
      eigv r@R{} = dup r
      eigv i@I{} = dup i
      eigv c@C{} = dup c
      eigv v@V3{} = eigvalsSpecial toR v -- eigvalsR v
      eigv pv@PV{} = eigvalsSpecial toR pv -- eigvalsR pv
      eigv bv@BV{} = eigvalsSpecial toI bv -- eigvalsI bv
      eigv tpv@TPV{} = eigvalsSpecial toI tpv -- eigvalsI tpv
      eigv h@H{} = eigvalsSpecial toC h -- eigvalsC h
      eigv od@ODD{} = eigvalsSpecial toC od -- eigvalsC od
      eigv cliff
        | hasNilpotent cliff = dup.reduce.toC $ cliff  -- this case is actually nilpotent
        | isColinear cliff = eigvalsSpecial toC cliff  -- eigvalsC bpv
        | otherwise =                           -- transform it so it will be colinear
            let (BPV a1 a2 a3 a23 a31 a12) = toBPV cliff
                boost = boost2colinear a1 a2 a3 a23 a31 a12
            in eigvalsSpecial toC (bar boost * cliff * boost) -- eigvalsC d
  --
  
  
  dup :: Cl3 -> (Cl3,Cl3)
  dup cliff = (cliff, cliff)
  
  -- | 'eigvalsSpecial' helper function to calculate Eigenvalues
  eigvalsSpecial :: (Cl3 -> Cl3) -> Cl3 -> (Cl3,Cl3)
  eigvalsSpecial toSpecial cliffor =
    let (_,_,eig1,eig2) = projEigs toSpecial cliffor
    in (eig1,eig2)
  
  
  -- | 'project' makes a projector based off of the vector content of the Cliffor.
  project :: Cl3 -> Cl3  -- PV<:Cl3
  project R{} = PV 0.5 0 0 0.5   -- default to e3 direction
  project I{} = PV 0.5 0 0 0.5   -- default to e3 direction
  project C{} = PV 0.5 0 0 0.5   -- default to e3 direction
  project (V3 a1 a2 a3) = triDProj a1 a2 a3   -- proj v@V3{} = 0.5 + 0.5*signum v
  project (PV _ a1 a2 a3) = triDProj a1 a2 a3   -- proj pv@PV{} = 0.5 + 0.5*(signum.toV3 $ pv)
  project (ODD a1 a2 a3 _) = triDProj a1 a2 a3   -- od@ODD{} = 0.5 + 0.5*(signum.toV3 $ od)
  project (BV a23 a31 a12) = triDProj a23 a31 a12   -- bv@BV{} = 0.5 + 0.5*(mIx.signum $ bv)
  project (H _ a23 a31 a12) = triDProj a23 a31 a12   -- h@H{} = 0.5 + 0.5*(mIx.signum.toBV $ h)
  project (TPV a23 a31 a12 _) = triDProj a23 a31 a12   -- tpv@TPV{} = 0.5 + 0.5*(mIx.signum.toBV $ tpv)
  project (BPV a1 a2 a3 a23 a31 a12) = biTriDProj a1 a2 a3 a23 a31 a12
  project (APS _ a1 a2 a3 a23 a31 a12 _) = biTriDProj a1 a2 a3 a23 a31 a12
  
  
  
  -- If Dot product is negative or zero we have a problem, if it is zero
  -- it either the vector or bivector par is zero or they are orthognal
  -- if the dot product is negative the vectors could be antiparallel
  biTriDProj :: Double -> Double -> Double -> Double -> Double -> Double -> Cl3  -- PV<:Cl3
  biTriDProj a1 a2 a3 a23 a31 a12 =
    let v3Mag = sqrt $ a1^2 + a2^2 + a3^2
        v3MagltTol = v3Mag < tol'
        halfInvV3Mag = recip v3Mag / 2
        bvMag = sqrt $ a23^2 + a31^2 + a12^2
        bvMagltTol = bvMag < tol'
        halfInvBVMag = recip bvMag / 2
        dotPos = (a1*a23) + (a2*a31) + (a3*a12) >= 0
        b1 = a1 + a23
        b2 = a2 + a31
        b3 = a3 + a12
        bHalfInvMag = (/2).recip.sqrt $ b1^2 + b2^2 + b3^2
        c1 = a1 - a23
        c2 = a2 - a31
        c3 = a3 - a12
        cHalfInvMag = (/2).recip.sqrt $ c1^2 + c2^2 + c3^2
    in if | v3MagltTol && bvMagltTol -> PV 0.5 0 0 0.5
          | bvMagltTol -> PV 0.5 (halfInvV3Mag * a1) (halfInvV3Mag * a2) (halfInvV3Mag * a3)
          | v3MagltTol -> PV 0.5 (halfInvBVMag * a23) (halfInvBVMag * a31) (halfInvBVMag * a12)
          | dotPos -> PV 0.5 (bHalfInvMag * b1) (bHalfInvMag * b2) (bHalfInvMag * b3)
          | otherwise -> PV 0.5 (cHalfInvMag * c1) (cHalfInvMag * c2) (cHalfInvMag * c3)
  
  
  -- | 'triDProj' a single 3 dimensional vector grade to a projector
  triDProj :: Double -> Double -> Double -> Cl3  -- PV<:Cl3
  triDProj v1 v2 v3 =
    let mag = sqrt $ v1^2 + v2^2 + v3^2
        halfInvMag = recip mag / 2
    in if mag == 0
       then PV 0.5 0 0 0.5
       else PV 0.5 (halfInvMag * v1) (halfInvMag * v2) (halfInvMag * v3)
  
  
  -- | 'boost2colinear' calculates a boost that is perpendicular to both the vector and bivector
  -- components of the cliffor, that will mix the vector and bivector parts such that the vector and bivector
  -- parts become colinear. This function is a simularity transform such that:
  --
  -- > cliffor = boost * colinear * bar boost
  --
  -- and returns the boost given the inputs.  First the boost must be calculated
  -- and then
  --
  -- > colinear = bar boost * cliffor * boost
  --
  -- and colinear will have colinear vector and bivector parts of the cliffor.
  -- This is somewhat simular to finding the drift frame for a static electromagnetic field.
  --
  -- > v = toV3 cliffor  -- extract the vector
  -- > bv = mIx.toBV $ cliffor  -- extract the bivector and turn it into a vector
  -- > invariant = ((2*).mIx.toBV $ v * bv) / (toR (v^2) + toR (bv^2))
  -- > boost = spectraldcmpSpecial toR (exp.(/4).atanh) invariant
  --
  boost2colinear :: Double -> Double -> Double -> Double -> Double -> Double -> Cl3  -- PV<:Cl3
  boost2colinear a1 a2 a3 a23 a31 a12 =
    let scale = recip $ a1^2 + a2^2 + a3^2 + a23^2 + a31^2 + a12^2
        b1 = scale * (a2*a12 - a3*a31)
        b2 = scale * (a3*a23 - a1*a12)
        b3 = scale * (a1*a31 - a2*a23)
        eig1 = (2*).sqrt $ b1^2 + b2^2 + b3^2
        eig2 = negate eig1
        transEig1 = exp.(/4).atanh $ eig1
        transEig2 = exp.(/4).atanh $ eig2
        sumTransEigs = (transEig1 - transEig2) * recip eig1
    in PV (0.5 * (transEig1 + transEig2)) (sumTransEigs * b1) (sumTransEigs * b2) (sumTransEigs * b3)
  
  
  -- | 'isColinear' takes a Cliffor and determines if either the vector part or the bivector part are
  -- zero or both aligned in the same direction.
  isColinear :: Cl3 -> Bool
  isColinear R{} = True
  isColinear V3{} = True
  isColinear BV{} = True
  isColinear I{} = True
  isColinear PV{} = True
  isColinear H{} = True
  isColinear C{} = True
  isColinear ODD{} = True
  isColinear TPV{} = True
  isColinear (BPV a1 a2 a3 a23 a31 a12) = colinearHelper a1 a2 a3 a23 a31 a12
  isColinear (APS _ a1 a2 a3 a23 a31 a12 _) = colinearHelper a1 a2 a3 a23 a31 a12
  
  colinearHelper :: Double -> Double -> Double -> Double -> Double -> Double -> Bool
  colinearHelper a1 a2 a3 a23 a31 a12 =
    let magV3 = sqrt $ a1^2 + a2^2 + a3^2
        invMagV3 = recip magV3
        magBV = sqrt $ a23^2 + a31^2 + a12^2
        invMagBV = recip magBV
        crss = sqrt (((invMagV3 * a2)*(invMagBV * a12) - (invMagV3 * a3)*(invMagBV * a31))^2 +
                     ((invMagV3 * a3)*(invMagBV * a23) - (invMagV3 * a1)*(invMagBV * a12))^2 +
                     ((invMagV3 * a1)*(invMagBV * a31) - (invMagV3 * a2)*(invMagBV * a23))^2)
    in magV3 == 0 ||     -- Zero Vector
       magBV == 0 ||     -- Zero Bivector
       crss <= tol'      -- Orthoganl part is zero-ish
  
  
  -- | 'hasNilpotent' takes a Cliffor and determines if the vector part and the bivector part are
  -- orthoganl and equal in magnitude, i.e. that it is simular to a nilpotent BPV.
  hasNilpotent :: Cl3 -> Bool
  hasNilpotent R{} = False
  hasNilpotent V3{} = False
  hasNilpotent BV{} = False
  hasNilpotent I{} = False
  hasNilpotent PV{} = False
  hasNilpotent H{} = False
  hasNilpotent C{} = False
  hasNilpotent ODD{} = False
  hasNilpotent TPV{} = False
  hasNilpotent (BPV a1 a2 a3 a23 a31 a12) = nilpotentHelper a1 a2 a3 a23 a31 a12
  hasNilpotent (APS _ a1 a2 a3 a23 a31 a12 _) = nilpotentHelper a1 a2 a3 a23 a31 a12
  
  nilpotentHelper :: Double -> Double -> Double -> Double -> Double -> Double -> Bool
  nilpotentHelper a1 a2 a3 a23 a31 a12 =
    let magV3 = sqrt $ a1^2 + a2^2 + a3^2
        invMagV3 = recip magV3
        magBV = sqrt $ a23^2 + a31^2 + a12^2
        invMagBV = recip magV3
        magDiff = abs (magV3 - magBV)
        b1 = invMagV3 * a1
        b2 = invMagV3 * a2
        b3 = invMagV3 * a3
        b23 = invMagBV * a23
        b31 = invMagBV * a31
        b12 = invMagBV * a12
        c0 = b1*b1 + b2*b2 + b3*b3 - b23*b23 - b31*b31 - b12*b12
        c1 = b12*b2 - b2*b12 + b3*b31 - b31*b3
        c2 = b1*b12 - b12*b1 - b3*b23 + b23*b3
        c3 = b31*b1 - b1*b31 + b2*b23 - b23*b2
        c23 = b2*b3 - b3*b2 - b31*b12 + b12*b31
        c31 = b3*b1 - b1*b3 + b23*b12 - b12*b23
        c12 = b1*b2 - b2*b1 - b23*b31 + b31*b23
        c123 = b1*b23 + b23*b1 + b2*b31 + b31*b2 + b3*b12 + b12*b3
        x = sqrt ((c0*c1 + c123*c23)^2 + (c0*c2 + c123*c31)^2 + (c0*c3 + c123*c12)^2 +
                  (c2*c12 - c3*c31)^2 + (c3*c23 - c1*c12)^2 + (c1*c31 - c2*c23)^2)
        sqMag = sqrt (c0^2 + c1^2 + c2^2 + c3^2 + c23^2 + c31^2 + c12^2 + c123^2 + x + x)
    in magV3 /= 0 &&          -- Non-Zero Vector Part
       magBV /= 0 &&          -- Non-Zero Bivector Part
       magDiff <= tol' &&     -- Vector and Bivector are Equal Magnitude
       sqMag <= tol'          -- It's non-zero but squares to zero
  
  
  -- | 'projEigs' function returns complementary projectors and eigenvalues for a Cliffor with specialization.
  -- The Cliffor at this point is allready colinear and the Eigenvalue is known to be real, imaginary, or complex.
  projEigs :: (Cl3 -> Cl3) -> Cl3 -> (Cl3,Cl3,Cl3,Cl3)
  projEigs toSpecial cliffor =
    let p = project cliffor
        p_bar = bar p
        eig1 = 2 * toSpecial (p * cliffor * p)
        eig2 = 2 * toSpecial (p_bar * cliffor * p_bar)
    in (p,p_bar,eig1,eig2)
  
  -- | 'reduce' function reduces the number of grades in a specialized Cliffor if they
  -- are zero-ish
  reduce :: Cl3 -> Cl3
  reduce cliff
    | abs cliff <= tol = R 0
    | otherwise = go_reduce cliff
      where
        go_reduce r@R{} = r
        go_reduce v@V3{} = v
        go_reduce bv@BV{} = bv
        go_reduce i@I{} = i
        go_reduce pv@PV{}
          | abs (toV3 pv) <= tol = toR pv
          | abs (toR pv) <= tol = toV3 pv
          | otherwise = pv
        go_reduce h@H{}
          | abs (toBV h) <= tol = toR h
          | abs (toR h) <= tol = toBV h
          | otherwise = h
        go_reduce c@C{}
          | abs (toI c) <= tol = toR c
          | abs (toR c) <= tol = toI c
          | otherwise = c
        go_reduce bpv@BPV{}
          | abs (toBV bpv) <= tol = toV3 bpv
          | abs (toV3 bpv) <= tol = toBV bpv
          | otherwise = bpv
        go_reduce od@ODD{}
          | abs (toI od) <= tol = toV3 od
          | abs (toV3 od) <= tol = toI od
          | otherwise = od
        go_reduce tpv@TPV{}
          | abs (toBV tpv) <= tol = toI tpv
          | abs (toI tpv) <= tol = toBV tpv
          | otherwise = tpv
        go_reduce aps@APS{}
          | abs (toBPV aps) <= tol = go_reduce (toC aps)
          | abs (toODD aps) <= tol = go_reduce (toH aps)
          | abs (toTPV aps) <= tol = go_reduce (toPV aps)
          | abs (toC aps) <= tol = go_reduce (toBPV aps)
          | abs (toH aps) <= tol = go_reduce (toODD aps)
          | abs (toPV aps) <= tol = go_reduce (toTPV aps)
          | otherwise = aps
  
  
  -- | 'mIx' a more effecient '\x -> I (-1) * x' typically useful for converting a
  -- Bivector to a Vector in the same direction. Related to Hodge Dual and/or
  -- Inverse Hodge Star.
  mIx :: Cl3 -> Cl3
  mIx (R a0) = I (negate a0)
  mIx (V3 a1 a2 a3) = BV (negate a1) (negate a2) (negate a3)
  mIx (BV a23 a31 a12) = V3 a23 a31 a12
  mIx (I a123) = R a123
  mIx (PV a0 a1 a2 a3) = TPV (negate a1) (negate a2) (negate a3) (negate a0)
  mIx (H a0 a23 a31 a12) = ODD a23 a31 a12 (negate a0)
  mIx (C a0 a123) = C a123 (negate a0)
  mIx (BPV a1 a2 a3 a23 a31 a12) = BPV a23 a31 a12 (negate a1) (negate a2) (negate a3)
  mIx (ODD a1 a2 a3 a123) = H a123 (negate a1) (negate a2) (negate a3)
  mIx (TPV a23 a31 a12 a123) = PV a123 a23 a31 a12
  mIx (APS a0 a1 a2 a3 a23 a31 a12 a123) = APS a123 a23 a31 a12 (negate a1) (negate a2) (negate a3) (negate a0)
  
  -- | 'timesI' is a more effecient '\x -> I 1 * x'
  timesI :: Cl3 -> Cl3
  timesI (R a0) = I a0
  timesI (V3 a1 a2 a3) = BV a1 a2 a3
  timesI (BV a23 a31 a12) = V3 (negate a23) (negate a31) (negate a12)
  timesI (I a123) = R (negate a123)
  timesI (PV a0 a1 a2 a3) = TPV a1 a2 a3 a0
  timesI (H a0 a23 a31 a12) = ODD (negate a23) (negate a31) (negate a12) a0
  timesI (C a0 a123) = C (negate a123) a0
  timesI (BPV a1 a2 a3 a23 a31 a12) = BPV (negate a23) (negate a31) (negate a12) a1 a2 a3
  timesI (ODD a1 a2 a3 a123) = H (negate a123) a1 a2 a3
  timesI (TPV a23 a31 a12 a123) = PV (negate a123) (negate a23) (negate a31) (negate a12)
  timesI (APS a0 a1 a2 a3 a23 a31 a12 a123) = APS (negate a123) (negate a23) (negate a31) (negate a12) a1 a2 a3 a0
  
  -- | 'abssignum' is a more effecient '\cl3 -> (abs cl3, signum cl3)'
  -- So 'abs' is always R and 'signum' is the same type of constructor as the input
  -- 'signum' is the element divided by its largest singular value 'abs'
  abssignum :: Cl3 -> (Cl3,Cl3)
  abssignum cl3 =
    let (R m0) = absolute cl3
    in if m0 == 0
       then (R 0, R 0) -- (abs 0 == 0, signum 0 == 0)
       else (R m0, cl3/(R m0))
  
  absolute :: Cl3 -> Cl3
  absolute (R a0) = R (abs a0)
  absolute (V3 a1 a2 a3) = let m = rss3 a1 a2 a3 in R m
  absolute (BV a23 a31 a12) = let m = rss3 a23 a31 a12 in R m
  absolute (I a123) = R (abs a123)
  absolute (PV a0 a1 a2 a3) = let m = reimMag a0 a1 a2 a3 in R m
  absolute (H a0 a23 a31 a12) = let m = rss4 a0 a23 a31 a12 in R m
  absolute (C a0 a123) = let m = rss2 a0 a123 in R m
  absolute (BPV a1 a2 a3 a23 a31 a12) = let mag0 = rss3 (a1*a31 - a2*a23) (a1*a12 - a3*a23) (a2*a12 - a3*a31)
                                            m = sqrt $ a1^2 + a23^2 + a2^2 + a31^2 + a3^2 + a12^2 + 2*mag0 in R m
  absolute (ODD a1 a2 a3 a123) = let m = rss4 a1 a2 a3 a123 in R m
  absolute (TPV a23 a31 a12 a123) = let m = reimMag a123 a23 a31 a12 in R m
  absolute (APS a0 a1 a2 a3 a23 a31 a12 a123) = let mag0 = sqrt $ (a0*a1 + a123*a23)^2 + (a0*a2 + a123*a31)^2 + (a0*a3 + a123*a12)^2 + (a2*a12 - a3*a31)^2 + (a3*a23 - a1*a12)^2 + (a1*a31 - a2*a23)^2
                                                    m = sqrt $ a0^2 + a1^2 + a2^2 + a3^2 + a23^2 + a31^2 + a12^2 + a123^2 + 2*mag0 in R m
  
  rss2 :: Double -> Double -> Double
  rss2 a0 a123 = sqrt $ a0^2 + a123^2
  
  rss3 :: Double -> Double -> Double -> Double
  rss3 x y z = sqrt $ x^2 + y^2 + z^2
  
  rss4 :: Double -> Double -> Double -> Double -> Double
  rss4 t x y z = sqrt $ t^2 + x^2 + y^2 + z^2
  
  
  
  #ifdef O_LIQUID
  tol :: Cl3
  tol = R 0
  
  tol' :: Double
  tol' = 0
  #else
  -- | 'tol' currently 128*eps
  tol :: Cl3
  {-# INLINE tol #-}
  tol = R 1.4210854715202004e-14
  
  tol' :: Double
  {-# INLINE tol' #-}
  tol' = 1.4210854715202004e-14
  #endif
  
  -- | 'bar' is a Clifford Conjugate, the vector grades are negated
  bar :: Cl3 -> Cl3
  bar (R a0) = R a0
  bar (V3 a1 a2 a3) = V3 (negate a1) (negate a2) (negate a3)
  bar (BV a23 a31 a12) = BV (negate a23) (negate a31) (negate a12)
  bar (I a123) = I a123
  bar (PV a0 a1 a2 a3) = PV a0 (negate a1) (negate a2) (negate a3)
  bar (H a0 a23 a31 a12) = H a0 (negate a23) (negate a31) (negate a12)
  bar (C a0 a123) = C a0 a123
  bar (BPV a1 a2 a3 a23 a31 a12) = BPV (negate a1) (negate a2) (negate a3) (negate a23) (negate a31) (negate a12)
  bar (ODD a1 a2 a3 a123) = ODD (negate a1) (negate a2) (negate a3) a123
  bar (TPV a23 a31 a12 a123) = TPV (negate a23) (negate a31) (negate a12) a123
  bar (APS a0 a1 a2 a3 a23 a31 a12 a123) = APS a0 (negate a1) (negate a2) (negate a3) (negate a23) (negate a31) (negate a12) a123
  
  -- | 'dag' is the Complex Conjugate, the imaginary grades are negated
  dag :: Cl3 -> Cl3
  dag (R a0) = R a0
  dag (V3 a1 a2 a3) = V3 a1 a2 a3
  dag (BV a23 a31 a12) = BV (negate a23) (negate a31) (negate a12)
  dag (I a123) = I (negate a123)
  dag (PV a0 a1 a2 a3) =  PV a0 a1 a2 a3
  dag (H a0 a23 a31 a12) = H a0 (negate a23) (negate a31) (negate a12)
  dag (C a0 a123) = C a0 (negate a123)
  dag (BPV a1 a2 a3 a23 a31 a12) = BPV a1 a2 a3 (negate a23) (negate a31) (negate a12)
  dag (ODD a1 a2 a3 a123) = ODD a1 a2 a3 (negate a123)
  dag (TPV a23 a31 a12 a123) = TPV (negate a23) (negate a31) (negate a12) (negate a123)
  dag (APS a0 a1 a2 a3 a23 a31 a12 a123) = APS a0 a1 a2 a3 (negate a23) (negate a31) (negate a12) (negate a123)
  
  ----------------------------------------------------------------------------------------------------------------
  -- the to... functions provide a lossy cast from one Cl3 constructor to another
  ---------------------------------------------------------------------------------------------------------------
  -- | 'toR' takes any Cliffor and returns the R portion
  toR :: Cl3 -> Cl3
  toR (R a0) = R a0
  toR V3{} = R 0
  toR BV{} = R 0
  toR I{} = R 0
  toR (PV a0 _ _ _) = R a0
  toR (H a0 _ _ _) = R a0
  toR (C a0 _) = R a0
  toR BPV{} = R 0
  toR ODD{} = R 0
  toR TPV{} = R 0
  toR (APS a0 _ _ _ _ _ _ _) = R a0
  
  -- | 'toV3' takes any Cliffor and returns the V3 portion
  toV3 :: Cl3 -> Cl3
  toV3 R{} = V3 0 0 0
  toV3 (V3 a1 a2 a3) = V3 a1 a2 a3
  toV3 BV{} = V3 0 0 0
  toV3 I{} = V3 0 0 0
  toV3 (PV _ a1 a2 a3) = V3 a1 a2 a3
  toV3 H{} = V3 0 0 0
  toV3 C{} = V3 0 0 0
  toV3 (BPV a1 a2 a3 _ _ _) = V3 a1 a2 a3
  toV3 (ODD a1 a2 a3 _) = V3 a1 a2 a3
  toV3 TPV{} = V3 0 0 0
  toV3 (APS _ a1 a2 a3 _ _ _ _) = V3 a1 a2 a3
  
  -- | 'toBV' takes any Cliffor and returns the BV portion
  toBV :: Cl3 -> Cl3
  toBV R{} = BV 0 0 0
  toBV V3{} = BV 0 0 0
  toBV (BV a23 a31 a12) = BV a23 a31 a12
  toBV I{} = BV 0 0 0
  toBV PV{} = BV 0 0 0
  toBV (H _ a23 a31 a12) = BV a23 a31 a12
  toBV C{} = BV 0 0 0
  toBV (BPV _ _ _ a23 a31 a12) = BV a23 a31 a12
  toBV ODD{} = BV 0 0 0
  toBV (TPV a23 a31 a12 _) = BV a23 a31 a12
  toBV (APS _ _ _ _ a23 a31 a12 _) = BV a23 a31 a12
  
  -- | 'toI' takes any Cliffor and returns the I portion
  toI :: Cl3 -> Cl3
  toI R{} = I 0
  toI V3{} = I 0
  toI BV{} = I 0
  toI (I a123) = I a123
  toI PV{} = I 0
  toI H{} = I 0
  toI (C _ a123) = I a123
  toI BPV{} = I 0
  toI (ODD _ _ _ a123) = I a123
  toI (TPV _ _ _ a123) = I a123
  toI (APS _ _ _ _ _ _ _ a123) = I a123
  
  -- | 'toPV' takes any Cliffor and returns the PV poriton
  toPV :: Cl3 -> Cl3
  toPV (R a0) = PV a0 0 0 0
  toPV (V3 a1 a2 a3) = PV 0 a1 a2 a3
  toPV BV{} = PV 0 0 0 0
  toPV I{} = PV 0 0 0 0
  toPV (PV a0 a1 a2 a3) = PV a0 a1 a2 a3
  toPV (H a0 _ _ _) = PV a0 0 0 0
  toPV (C a0 _) = PV a0 0 0 0
  toPV (BPV a1 a2 a3 _ _ _) = PV 0 a1 a2 a3
  toPV (ODD a1 a2 a3 _) = PV a1 a2 a3 0
  toPV TPV{} = PV 0 0 0 0
  toPV (APS a0 a1 a2 a3 _ _ _ _) = PV a0 a1 a2 a3
  
  -- | 'toH' takes any Cliffor and returns the H portion
  toH :: Cl3 -> Cl3
  toH (R a0) = H a0 0 0 0
  toH V3{} = H 0 0 0 0
  toH (BV a23 a31 a12) = H 0 a23 a31 a12
  toH (I _) = H 0 0 0 0
  toH (PV a0 _ _ _) = H a0 0 0 0
  toH (H a0 a23 a31 a12) = H a0 a23 a31 a12
  toH (C a0 _) = H a0 0 0 0
  toH (BPV _ _ _ a23 a31 a12) = H 0 a23 a31 a12
  toH ODD{} = H 0 0 0 0
  toH (TPV a23 a31 a12 _) = H 0 a23 a31 a12
  toH (APS a0 _ _ _ a23 a31 a12 _) = H a0 a23 a31 a12
  
  -- | 'toC' takes any Cliffor and returns the C portion
  toC :: Cl3 -> Cl3
  toC (R a0) = C a0 0
  toC V3{} = C 0 0
  toC BV{} = C 0 0
  toC (I a123) = C 0 a123
  toC (PV a0 _ _ _) = C a0 0
  toC (H a0 _ _ _) = C a0 0
  toC (C a0 a123) = C a0 a123
  toC BPV{} = C 0 0
  toC (ODD _ _ _ a123) = C 0 a123
  toC (TPV _ _ _ a123) = C 0 a123
  toC (APS a0 _ _ _ _ _ _ a123) = C a0 a123
  
  -- | 'toBPV' takes any Cliffor and returns the BPV portion
  toBPV :: Cl3 -> Cl3
  toBPV R{} = BPV 0 0 0 0 0 0
  toBPV (V3 a1 a2 a3) = BPV a1 a2 a3 0 0 0
  toBPV (BV a23 a31 a12) = BPV 0 0 0 a23 a31 a12
  toBPV I{} = BPV 0 0 0 0 0 0
  toBPV (PV _ a1 a2 a3) = BPV a1 a2 a3 0 0 0
  toBPV (H _ a23 a31 a12) = BPV 0 0 0 a23 a31 a12
  toBPV C{} = BPV 0 0 0 0 0 0
  toBPV (BPV a1 a2 a3 a23 a31 a12) = BPV a1 a2 a3 a23 a31 a12
  toBPV (ODD a1 a2 a3 _) = BPV a1 a2 a3 0 0 0
  toBPV (TPV a23 a31 a12 _) = BPV 0 0 0 a23 a31 a12
  toBPV (APS _ a1 a2 a3 a23 a31 a12 _) = BPV a1 a2 a3 a23 a31 a12
  
  -- | 'toODD' takes any Cliffor and returns the ODD portion
  toODD :: Cl3 -> Cl3
  toODD R{} = ODD 0 0 0 0
  toODD (V3 a1 a2 a3) = ODD a1 a2 a3 0
  toODD BV{} = ODD 0 0 0 0
  toODD (I a123) = ODD 0 0 0 a123
  toODD (PV _ a1 a2 a3) = ODD a1 a2 a3 0
  toODD H{} = ODD 0 0 0 0
  toODD (C _ a123) = ODD 0 0 0 a123
  toODD (BPV a1 a2 a3 _ _ _) = ODD a1 a2 a3 0
  toODD (ODD a1 a2 a3 a123) = ODD a1 a2 a3 a123
  toODD (TPV _ _ _ a123) = ODD 0 0 0 a123
  toODD (APS _ a1 a2 a3 _ _ _ a123) = ODD a1 a2 a3 a123
  
  -- | 'toTPV' takes any Cliffor and returns the TPV portion
  toTPV :: Cl3 -> Cl3
  toTPV R{} = TPV 0 0 0 0
  toTPV V3{} = TPV 0 0 0 0
  toTPV (BV a23 a31 a12) = TPV a23 a31 a12 0
  toTPV (I a123) = TPV 0 0 0 a123
  toTPV PV{} = TPV 0 0 0 0
  toTPV (H _ a23 a31 a12) = TPV a23 a31 a12 0
  toTPV (C _ a123) = TPV 0 0 0 a123
  toTPV (BPV _ _ _ a23 a31 a12) = TPV a23 a31 a12 0
  toTPV (ODD _ _ _ a123) = TPV 0 0 0 a123
  toTPV (TPV a23 a31 a12 a123) = TPV a23 a31 a12 a123
  toTPV (APS _ _ _ _ a23 a31 a12 a123) = TPV a23 a31 a12 a123
  
  -- | 'toAPS' takes any Cliffor and returns the APS portion
  toAPS :: Cl3 -> Cl3
  toAPS (R a0) = APS a0 0 0 0 0 0 0 0
  toAPS (V3 a1 a2 a3) = APS 0 a1 a2 a3 0 0 0 0
  toAPS (BV a23 a31 a12) = APS 0 0 0 0 a23 a31 a12 0
  toAPS (I a123) = APS 0 0 0 0 0 0 0 a123
  toAPS (PV a0 a1 a2 a3) = APS a0 a1 a2 a3 0 0 0 0
  toAPS (H a0 a23 a31 a12) = APS a0 0 0 0 a23 a31 a12 0
  toAPS (C a0 a123) = APS a0 0 0 0 0 0 0 a123
  toAPS (BPV a1 a2 a3 a23 a31 a12) = APS 0 a1 a2 a3 a23 a31 a12 0
  toAPS (ODD a1 a2 a3 a123) = APS 0 a1 a2 a3 0 0 0 a123
  toAPS (TPV a23 a31 a12 a123) = APS 0 0 0 0 a23 a31 a12 a123
  toAPS (APS a0 a1 a2 a3 a23 a31 a12 a123) = APS a0 a1 a2 a3 a23 a31 a12 a123
  
  -- derivatives of the functions in the Fractional Class for use in Jordan NF functon implemetnation
  recip' :: Cl3 -> Cl3
  recip' = negate.recip.(^2)   -- pole at 0
  
  exp' :: Cl3 -> Cl3
  exp' = exp
  
  log' :: Cl3 -> Cl3
  log' = recip  -- pole at 0
  
  sqrt' :: Cl3 -> Cl3
  sqrt' = (/2).recip.sqrt   -- pole at 0
  
  sin' :: Cl3 -> Cl3
  sin' = cos
  
  cos' :: Cl3 -> Cl3
  cos' = negate.sin
  
  tan' :: Cl3 -> Cl3
  tan' = recip.(^2).cos  -- pole at pi/2*n for all integers
  
  asin' :: Cl3 -> Cl3
  asin' = recip.sqrt.(1-).(^2)  -- pole at +/-1
  
  acos' :: Cl3 -> Cl3
  acos' = negate.recip.sqrt.(1-).(^2)  -- pole at +/-1
  
  atan' :: Cl3 -> Cl3
  atan' = recip.(1+).(^2)  -- pole at +/-i
  
  sinh' :: Cl3 -> Cl3
  sinh' = cosh
  
  cosh' :: Cl3 -> Cl3
  cosh' = sinh
  
  tanh' :: Cl3 -> Cl3
  tanh' = recip.(^2).cosh
  
  asinh' :: Cl3 -> Cl3
  asinh' = recip.sqrt.(1+).(^2)  -- pole at +/-i
  
  acosh' :: Cl3 -> Cl3
  acosh' x = recip $ sqrt (x - 1) * sqrt (x + 1)  -- pole at +/-1
  
  atanh' :: Cl3 -> Cl3
  atanh' = recip.(1-).(^2)  -- pole at +/-1
  
  
  #ifndef O_NO_STORABLE
  -------------------------------------------------------------------
  -- 
  -- Instance of Cl3 types with the "Foreign.Storable" library.
  --  
  -- For use with high performance data structures like Data.Vector.Storable
  -- or Data.Array.Storable
  -- 
  -------------------------------------------------------------------
  
  -- | Cl3 instance of Storable uses the APS constructor as its standard interface.
  -- "peek" returns a cliffor constructed with APS. "poke" converts a cliffor to APS.
  -- For a more compact storing of constructors other than APS use the storable
  -- subtypes Cl3_R, Cl3_V3, Cl3_BV, Cl3_I, Cl3_PV, Cl3_H, Cl3_C, Cl3_BPV,
  -- Cl3_ODD, Cl3_TPV.
  instance Storable Cl3 where
    sizeOf _ = 8 * sizeOf (undefined :: Double)
    alignment _ = sizeOf (undefined :: Double)
    peek ptr = do
      a0 <- peek (offset 0)
      a1 <- peek (offset 1)
      a2 <- peek (offset 2)
      a3 <- peek (offset 3)
      a23 <- peek (offset 4)
      a31 <- peek (offset 5)
      a12 <- peek (offset 6)
      a123 <- peek (offset 7)
      return $ APS a0 a1 a2 a3 a23 a31 a12 a123
        where
          offset i = (castPtr ptr :: Ptr Double) `plusPtr` (i*8)
    
    poke ptr (toAPS -> APS a0 a1 a2 a3 a23 a31 a12 a123) = do
      poke (offset 0) a0
      poke (offset 1) a1
      poke (offset 2) a2
      poke (offset 3) a3
      poke (offset 4) a23
      poke (offset 5) a31
      poke (offset 6) a12
      poke (offset 7) a123
        where
          offset i = (castPtr ptr :: Ptr Double) `plusPtr` (i*8)
    poke _ _ = error "Serious Issues with poke in Cl3.Storable"
  
  
  -- | 'Cl3_R' a compact storable data type for R.
  data Cl3_R where
    Cl3_R :: !Double -> Cl3_R
  
  -- | 'toCl3_R' converts a Cl3 value constructed with R to its compact form.
  toCl3_R :: Cl3 -> Cl3_R
  toCl3_R (R a0) = Cl3_R a0
  toCl3_R err = error $ "Please don't try and cast something that's not R to Cl3_R, Got: " ++ show err
  
  -- | 'fromCl3_R' converts the compact Cl3_R type back to a Cl3 type.
  fromCl3_R :: Cl3_R -> Cl3
  fromCl3_R (Cl3_R a0) = R a0
  
  instance Show Cl3_R where
    show = show.fromCl3_R
  
  #ifndef O_NO_DERIVED
  instance Read Cl3_R where
    readPrec = toCl3_R <$> readPrec
  #endif
  
  instance Storable Cl3_R where
    sizeOf _ = sizeOf (undefined :: Double)
    alignment _ = sizeOf (undefined :: Double)
    peek ptr = do
      a0 <- peek offset
      return $ Cl3_R a0
        where
          offset = (castPtr ptr :: Ptr Double)
  
    poke ptr (Cl3_R a0) = do
      poke offset a0
        where
          offset = (castPtr ptr :: Ptr Double)
  
  
  -- | 'Cl3_V3' a compact storable data type for V3.
  data Cl3_V3 where
    Cl3_V3 :: !Double -> !Double -> !Double -> Cl3_V3
  
  -- | 'toCl3_V3' converts a Cl3 value constructed with V3 to its compact form.
  toCl3_V3 :: Cl3 -> Cl3_V3
  toCl3_V3 (V3 a1 a2 a3) = Cl3_V3 a1 a2 a3
  toCl3_V3 err = error $ "Please don't try and cast something that's not V3 to Cl3_V3, Got: " ++ show err
  
  -- | 'fromCl3_V3' converts the compact Cl3_V3 type back to a Cl3 type.
  fromCl3_V3 :: Cl3_V3 -> Cl3
  fromCl3_V3 (Cl3_V3 a1 a2 a3) = V3 a1 a2 a3
  
  instance Show Cl3_V3 where
    show = show.fromCl3_V3
  
  #ifndef O_NO_DERIVED
  instance Read Cl3_V3 where
    readPrec = toCl3_V3 <$> readPrec
  #endif
  
  instance Storable Cl3_V3 where
    sizeOf _ = 3 * sizeOf (undefined :: Double)
    alignment _ = sizeOf (undefined :: Double)
    peek ptr = do
      a1 <- peek (offset 0)
      a2 <- peek (offset 1)
      a3 <- peek (offset 2)
      return $ Cl3_V3 a1 a2 a3
        where
          offset i = (castPtr ptr :: Ptr Double) `plusPtr` (i*8)
  
    poke ptr (Cl3_V3 a1 a2 a3) = do
      poke (offset 0) a1
      poke (offset 1) a2
      poke (offset 2) a3
        where
          offset i = (castPtr ptr :: Ptr Double) `plusPtr` (i*8)
  
  
  -- | 'Cl3_BV' a compact storable data type for BV.
  data Cl3_BV where
    Cl3_BV :: !Double -> !Double -> !Double -> Cl3_BV
  
  -- | 'toCl3_BV' converts a Cl3 value constructed with BV to its compact form.
  toCl3_BV :: Cl3 -> Cl3_BV
  toCl3_BV (BV a23 a31 a12) = Cl3_BV a23 a31 a12
  toCl3_BV err = error $ "Please don't try and cast something that's not BV to Cl3_BV, Got: " ++ show err
  
  -- | 'fromCl3_BV' converts the compact Cl3_BV type back to a Cl3 type.
  fromCl3_BV :: Cl3_BV -> Cl3
  fromCl3_BV (Cl3_BV a23 a31 a12) = BV a23 a31 a12
  
  instance Show Cl3_BV where
    show = show.fromCl3_BV
  
  #ifndef O_NO_DERIVED
  instance Read Cl3_BV where
    readPrec = toCl3_BV <$> readPrec
  #endif
  
  instance Storable Cl3_BV where
    sizeOf _ = 3 * sizeOf (undefined :: Double)
    alignment _ = sizeOf (undefined :: Double)
    peek ptr = do
      a23 <- peek (offset 0)
      a31 <- peek (offset 1)
      a12 <- peek (offset 2)
      return $ Cl3_BV a23 a31 a12
        where
          offset i = (castPtr ptr :: Ptr Double) `plusPtr` (i*8)
  
    poke ptr (Cl3_BV a23 a31 a12) = do
      poke (offset 0) a23
      poke (offset 1) a31
      poke (offset 2) a12
        where
          offset i = (castPtr ptr :: Ptr Double) `plusPtr` (i*8)
  
  
  -- | 'Cl3_I' a compact storable data type for I.
  data Cl3_I where
    Cl3_I :: !Double -> Cl3_I
  
  -- | 'toCl3_I' converts a Cl3 value constructed with I to its compact form.
  toCl3_I :: Cl3 -> Cl3_I
  toCl3_I (I a123) = Cl3_I a123
  toCl3_I err = error $ "Please don't try and cast something that's not R to Cl3_R, Got: " ++ show err
  
  -- | 'fromCl3_I' converts the compact Cl3_I type back to a Cl3 type.
  fromCl3_I :: Cl3_I -> Cl3
  fromCl3_I (Cl3_I a123) = I a123
  
  instance Show Cl3_I where
    show = show.fromCl3_I
  
  #ifndef O_NO_DERIVED
  instance Read Cl3_I where
    readPrec = toCl3_I <$> readPrec
  #endif
  
  instance Storable Cl3_I where
    sizeOf _ = sizeOf (undefined :: Double)
    alignment _ = sizeOf (undefined :: Double)
    peek ptr = do
      a123 <- peek offset
      return $ Cl3_I a123
        where
          offset = (castPtr ptr :: Ptr Double)
  
    poke ptr (Cl3_I a123) = do
      poke offset a123
        where
          offset = (castPtr ptr :: Ptr Double)
  
  
  -- | 'Cl3_PV' a compact storable data type for PV.
  data Cl3_PV where
    Cl3_PV :: !Double -> !Double -> !Double -> !Double -> Cl3_PV
  
  -- | 'toCl3_PV' converts a Cl3 value constructed with PV to its compact form.
  toCl3_PV :: Cl3 -> Cl3_PV
  toCl3_PV (PV a0 a1 a2 a3) = Cl3_PV a0 a1 a2 a3
  toCl3_PV err = error $ "Please don't try and cast something that's not PV to Cl3_PV, Got: " ++ show err
  
  -- | 'fromCl3_PV' converts the compact Cl3_PV type back to a Cl3 type.
  fromCl3_PV :: Cl3_PV -> Cl3
  fromCl3_PV (Cl3_PV a0 a1 a2 a3) = PV a0 a1 a2 a3
  
  instance Show Cl3_PV where
    show = show.fromCl3_PV
  
  #ifndef O_NO_DERIVED
  instance Read Cl3_PV where
    readPrec = toCl3_PV <$> readPrec
  #endif
  
  instance Storable Cl3_PV where
    sizeOf _ = 4 * sizeOf (undefined :: Double)
    alignment _ = sizeOf (undefined :: Double)
    peek ptr = do
      a0 <- peek (offset 0)
      a1 <- peek (offset 1)
      a2 <- peek (offset 2)
      a3 <- peek (offset 4)
      return $ Cl3_PV a0 a1 a2 a3
        where
          offset i = (castPtr ptr :: Ptr Double) `plusPtr` (i*8)
  
    poke ptr (Cl3_PV a0 a1 a2 a3) = do
      poke (offset 0) a0
      poke (offset 1) a1
      poke (offset 2) a2
      poke (offset 3) a3
        where
          offset i = (castPtr ptr :: Ptr Double) `plusPtr` (i*8)
  
  
  -- | 'Cl3_H' a compact storable data type for H.
  data Cl3_H where
    Cl3_H :: !Double -> !Double -> !Double -> !Double -> Cl3_H
  
  -- | 'toCl3_H' converts a Cl3 value constructed with H to its compact form.
  toCl3_H :: Cl3 -> Cl3_H
  toCl3_H (H a0 a23 a31 a12) = Cl3_H a0 a23 a31 a12
  toCl3_H err = error $ "Please don't try and cast something that's not H to Cl3_H, Got: " ++ show err
  
  -- | 'fromCl3_H' converts the compact Cl3_H type back to a Cl3 type.
  fromCl3_H :: Cl3_H -> Cl3
  fromCl3_H (Cl3_H a0 a23 a31 a12) = H a0 a23 a31 a12
  
  instance Show Cl3_H where
    show = show.fromCl3_H
  
  #ifndef O_NO_DERIVED
  instance Read Cl3_H where
    readPrec = toCl3_H <$> readPrec
  #endif
  
  instance Storable Cl3_H where
    sizeOf _ = 4 * sizeOf (undefined :: Double)
    alignment _ = sizeOf (undefined :: Double)
    peek ptr = do
      a0 <- peek (offset 0)
      a23 <- peek (offset 1)
      a31 <- peek (offset 2)
      a12 <- peek (offset 3)
      return $ Cl3_H a0 a23 a31 a12
        where
          offset i = (castPtr ptr :: Ptr Double) `plusPtr` (i*8)
  
    poke ptr (Cl3_H a0 a23 a31 a12) = do
      poke (offset 0) a0
      poke (offset 1) a23
      poke (offset 2) a31
      poke (offset 3) a12
        where
          offset i = (castPtr ptr :: Ptr Double) `plusPtr` (i*8)
  
  
  -- | 'Cl3_C' a compact storable data type for C.
  data Cl3_C where
    Cl3_C :: !Double -> !Double -> Cl3_C
  
  -- | 'toCl3_C' converts a Cl3 value constructed with C to its compact form.
  toCl3_C :: Cl3 -> Cl3_C
  toCl3_C (C a0 a123) = Cl3_C a0 a123
  toCl3_C err = error $ "Please don't try and cast something that's not C to Cl3_C, Got: " ++ show err
  
  -- | 'fromCl3_C' converts the compact Cl3_C type back to a Cl3 type.
  fromCl3_C :: Cl3_C -> Cl3
  fromCl3_C (Cl3_C a0 a123) = C a0 a123
  
  instance Show Cl3_C where
    show = show.fromCl3_C
  
  #ifndef O_NO_DERIVED
  instance Read Cl3_C where
    readPrec = toCl3_C <$> readPrec
  #endif
  
  instance Storable Cl3_C where
    sizeOf _ = 2 * sizeOf (undefined :: Double)
    alignment _ = sizeOf (undefined :: Double)
    peek ptr = do
      a0 <- peek (offset 0)
      a123 <- peek (offset 1)
      return $ Cl3_C a0 a123
        where
          offset i = (castPtr ptr :: Ptr Double) `plusPtr` (i*8)
  
    poke ptr (Cl3_C a0 a123) = do
      poke (offset 0) a0
      poke (offset 1) a123
        where
          offset i = (castPtr ptr :: Ptr Double) `plusPtr` (i*8)
  
  
  -- | 'Cl3_BPV' a compact storable data type for BPV.
  data Cl3_BPV where
    Cl3_BPV :: !Double -> !Double -> !Double -> !Double -> !Double -> !Double -> Cl3_BPV
  
  -- | 'toCl3_BPV' converts a Cl3 value constructed with BPV to its compact form.
  toCl3_BPV :: Cl3 -> Cl3_BPV
  toCl3_BPV (BPV a1 a2 a3 a23 a31 a12) = Cl3_BPV a1 a2 a3 a23 a31 a12
  toCl3_BPV err = error $ "Please don't try and cast something that's not BPV to Cl3_BPV, Got: " ++ show err
  
  -- | 'fromCl3_BPV' converts the compact Cl3_BPV type back to a Cl3 type.
  fromCl3_BPV :: Cl3_BPV -> Cl3
  fromCl3_BPV (Cl3_BPV a1 a2 a3 a23 a31 a12) = BPV a1 a2 a3 a23 a31 a12
  
  instance Show Cl3_BPV where
    show = show.fromCl3_BPV
  
  #ifndef O_NO_DERIVED
  instance Read Cl3_BPV where
    readPrec = toCl3_BPV <$> readPrec
  #endif
  
  instance Storable Cl3_BPV where
    sizeOf _ = 6 * sizeOf (undefined :: Double)
    alignment _ = sizeOf (undefined :: Double)
    peek ptr = do
      a1 <- peek (offset 0)
      a2 <- peek (offset 1)
      a3 <- peek (offset 2)
      a23 <- peek (offset 3)
      a31 <- peek (offset 4)
      a12 <- peek (offset 5)
      return $ Cl3_BPV a1 a2 a3 a23 a31 a12
        where
          offset i = (castPtr ptr :: Ptr Double) `plusPtr` (i*8)
  
    poke ptr (Cl3_BPV a1 a2 a3 a23 a31 a12) = do
      poke (offset 0) a1
      poke (offset 1) a2
      poke (offset 2) a3
      poke (offset 3) a23
      poke (offset 4) a31
      poke (offset 5) a12
        where
          offset i = (castPtr ptr :: Ptr Double) `plusPtr` (i*8)
  
  
  -- | 'Cl3_ODD' a compact storable data type for ODD.
  data Cl3_ODD where
    Cl3_ODD :: !Double -> !Double -> !Double -> !Double -> Cl3_ODD
  
  -- | 'toCl3_ODD' converts a Cl3 value constructed with ODD to its compact form.
  toCl3_ODD :: Cl3 -> Cl3_ODD
  toCl3_ODD (ODD a1 a2 a3 a123) = Cl3_ODD a1 a2 a3 a123
  toCl3_ODD err = error $ "Please don't try and cast something that's not ODD to Cl3_ODD, Got: " ++ show err
  
  -- | 'fromCl3_ODD' converts the compact Cl3_ODD type back to a Cl3 type.
  fromCl3_ODD :: Cl3_ODD -> Cl3
  fromCl3_ODD (Cl3_ODD a1 a2 a3 a123) = ODD a1 a2 a3 a123
  
  instance Show Cl3_ODD where
    show = show.fromCl3_ODD
  
  #ifndef O_NO_DERIVED
  instance Read Cl3_ODD where
    readPrec = toCl3_ODD <$> readPrec
  #endif
  
  instance Storable Cl3_ODD where
    sizeOf _ = 4 * sizeOf (undefined :: Double)
    alignment _ = sizeOf (undefined :: Double)
    peek ptr = do
      a1 <- peek (offset 0)
      a2 <- peek (offset 1)
      a3 <- peek (offset 2)
      a123 <- peek (offset 3)
      return $ Cl3_ODD a1 a2 a3 a123
        where
          offset i = (castPtr ptr :: Ptr Double) `plusPtr` (i*8)
  
    poke ptr (Cl3_ODD a1 a2 a3 a123) = do
      poke (offset 0) a1
      poke (offset 1) a2
      poke (offset 2) a3
      poke (offset 3) a123
        where
          offset i = (castPtr ptr :: Ptr Double) `plusPtr` (i*8)
  
  
  -- | 'Cl3_TPV' a compact storable data type for TPV.
  data Cl3_TPV where
    Cl3_TPV :: !Double -> !Double -> !Double -> !Double -> Cl3_TPV
  
  -- | 'toCl3_TPV' converts a Cl3 value constructed with TPV to its compact form.
  toCl3_TPV :: Cl3 -> Cl3_TPV
  toCl3_TPV (TPV a23 a31 a12 a123) = Cl3_TPV a23 a31 a12 a123
  toCl3_TPV err = error $ "Please don't try and cast something that's not TPV to Cl3_TPV, Got: " ++ show err
  
  -- | 'fromCl3_TPV' converts the compact Cl3_TPV type back to a Cl3 type.
  fromCl3_TPV :: Cl3_TPV -> Cl3
  fromCl3_TPV (Cl3_TPV a23 a31 a12 a123) = TPV a23 a31 a12 a123
  
  instance Show Cl3_TPV where
    show = show.fromCl3_TPV
  
  #ifndef O_NO_DERIVED
  instance Read Cl3_TPV where
    readPrec = toCl3_TPV <$> readPrec
  #endif
  
  instance Storable Cl3_TPV where
    sizeOf _ = 4 * sizeOf (undefined :: Double)
    alignment _ = sizeOf (undefined :: Double)
    peek ptr = do
      a23 <- peek (offset 0)
      a31 <- peek (offset 1)
      a12 <- peek (offset 2)
      a123 <- peek (offset 3)
      return $ Cl3_TPV a23 a31 a12 a123
        where
          offset i = (castPtr ptr :: Ptr Double) `plusPtr` (i*8)
  
    poke ptr (Cl3_TPV a23 a31 a12 a123) = do
      poke (offset 0) a23
      poke (offset 1) a31
      poke (offset 2) a12
      poke (offset 3) a123
        where
          offset i = (castPtr ptr :: Ptr Double) `plusPtr` (i*8)
  
  
  -- | 'Cl3_APS' a compact storable data type for APS.
  data Cl3_APS where
    Cl3_APS :: !Double -> !Double -> !Double -> !Double -> !Double -> !Double -> !Double -> !Double -> Cl3_APS
  
  -- | 'toCl3_APS' converts a Cl3 value constructed with APS to its compact form.
  toCl3_APS :: Cl3 -> Cl3_APS
  toCl3_APS (APS a0 a1 a2 a3 a23 a31 a12 a123) = Cl3_APS a0 a1 a2 a3 a23 a31 a12 a123
  toCl3_APS err = error $ "Please don't try and cast something that's not APS to Cl3_APS, Got: " ++ show err
  
  -- | 'fromCl3_APS' converts the compact Cl3_APS type back to a Cl3 type.
  fromCl3_APS :: Cl3_APS -> Cl3
  fromCl3_APS (Cl3_APS a0 a1 a2 a3 a23 a31 a12 a123) = APS a0 a1 a2 a3 a23 a31 a12 a123
  
  instance Show Cl3_APS where
    show = show.fromCl3_APS
  
  #ifndef O_NO_DERIVED
  instance Read Cl3_APS where
    readPrec = toCl3_APS <$> readPrec
  #endif
  
  instance Storable Cl3_APS where
    sizeOf _ = 8 * sizeOf (undefined :: Double)
    alignment _ = sizeOf (undefined :: Double)
    peek ptr = do
      a0 <- peek (offset 0)
      a1 <- peek (offset 1)
      a2 <- peek (offset 2)
      a3 <- peek (offset 3)
      a23 <- peek (offset 4)
      a31 <- peek (offset 5)
      a12 <- peek (offset 6)
      a123 <- peek (offset 7)
      return $ Cl3_APS a0 a1 a2 a3 a23 a31 a12 a123
        where
          offset i = (castPtr ptr :: Ptr Double) `plusPtr` (i*8)
  
    poke ptr (Cl3_APS a0 a1 a2 a3 a23 a31 a12 a123) = do
      poke (offset 0) a0
      poke (offset 1) a1
      poke (offset 2) a2
      poke (offset 3) a3
      poke (offset 4) a23
      poke (offset 5) a31
      poke (offset 6) a12
      poke (offset 7) a123
        where
          offset i = (castPtr ptr :: Ptr Double) `plusPtr` (i*8)
  
  
  #endif
  
  
  
  #ifndef O_NO_RANDOM
  -------------------------------------------------------------------
  -- 
  -- Random Instance of Cl3 types with the "System.Random" library.
  -- 
  --
  -- Random helper functions will be based on the "abs x * signum x" decomposition
  -- for the single grade elements. The "abs x" will be the random magnitude that
  -- is by the default [0,1), and the "signum x" part will be a random direction
  -- of a vector or the sign of a scalar. The multi-grade elements will be constructed from
  -- a combination of the single grade generators.  Each grade will be evenly
  -- distributed across the range.
  -- 
  -------------------------------------------------------------------
  
  -- | 'Random' instance for the 'System.Random' library
  instance Random Cl3 where
    randomR (minAbs,maxAbs) g =
      case randomR (fromEnum (minBound :: ConCl3), fromEnum (maxBound :: ConCl3)) g of
        (r, g') -> case toEnum r of
          ConR -> rangeR (minAbs,maxAbs) g'
          ConV3 -> rangeV3 (minAbs,maxAbs) g'
          ConBV -> rangeBV (minAbs,maxAbs) g'
          ConI -> rangeI (minAbs,maxAbs) g'
          ConPV -> rangePV (minAbs,maxAbs) g'
          ConH -> rangeH (minAbs,maxAbs) g'
          ConC -> rangeC (minAbs,maxAbs) g'
          ConBPV -> rangeBPV (minAbs,maxAbs) g'
          ConODD -> rangeODD (minAbs,maxAbs) g'
          ConTPV -> rangeTPV (minAbs,maxAbs) g'
          ConAPS -> rangeAPS (minAbs,maxAbs) g'
          ConProj -> rangeProjector (minAbs,maxAbs) g'
          ConNilpotent -> rangeNilpotent (minAbs,maxAbs) g'
          ConUnitary -> rangeUnitary (minAbs,maxAbs) g'
  
    random = randomR (0,1)
  
  
  
  -- | 'ConCl3' Bounded Enum Algebraic Data Type of constructors of Cl3
  data ConCl3 = ConR
              | ConV3
              | ConBV
              | ConI
              | ConPV
              | ConH
              | ConC
              | ConBPV
              | ConODD
              | ConTPV
              | ConAPS
              | ConProj
              | ConNilpotent
              | ConUnitary
    deriving (Bounded, Enum)
  
  
  
  
  -- | 'randR' random Real Scalar (Grade 0) with random magnitude and random sign
  randR :: RandomGen g => g -> (Cl3, g)
  randR = rangeR (0,1)
  
  
  -- | 'rangeR' random Real Scalar (Grade 0) with random magnitude within a range and a random sign
  rangeR :: RandomGen g => (Cl3, Cl3) -> g -> (Cl3, g)
  rangeR = scalarHelper R
  
  
  -- | 'randV3' random Vector (Grade 1) with random magnitude and random direction
  -- the direction is using spherical coordinates
  randV3 :: RandomGen g => g -> (Cl3, g)
  randV3 = rangeV3 (0,1)
  
  
  -- | 'rangeV3' random Vector (Grade 1) with random magnitude within a range and a random direction
  rangeV3 :: RandomGen g => (Cl3, Cl3) -> g -> (Cl3, g)
  rangeV3 = vectorHelper V3
  
  
  -- | 'randBV' random Bivector (Grade 2) with random magnitude and random direction
  -- the direction is using spherical coordinates
  randBV :: RandomGen g => g -> (Cl3, g)
  randBV = rangeBV (0,1)
  
  
  -- | 'rangeBV' random Bivector (Grade 2) with random magnitude in a range and a random direction
  rangeBV :: RandomGen g => (Cl3, Cl3) -> g -> (Cl3, g)
  rangeBV = vectorHelper BV
  
  
  -- | 'randI' random Imaginary Scalar (Grade 3) with random magnitude and random sign
  randI :: RandomGen g => g -> (Cl3, g)
  randI = rangeI (0,1)
  
  
  -- | 'rangeI' random Imaginary Scalar (Grade 3) with random magnitude within a range and random sign
  rangeI :: RandomGen g => (Cl3, Cl3) -> g -> (Cl3, g)
  rangeI = scalarHelper I
  
  
  -- | 'randPV' random Paravector made from random Grade 0 and Grade 1 elements
  randPV :: RandomGen g => g -> (Cl3, g)
  randPV = rangePV (0,1)
  
  
  -- | 'rangePV' random Paravector made from random Grade 0 and Grade 1 elements within a range
  rangePV :: RandomGen g => (Cl3, Cl3) -> g -> (Cl3, g)
  rangePV (lo, hi) g =
    let (R scale, g') = rangeR (lo, hi) g
        (R a0, g'') = randR g'
        (V3 a1 a2 a3, g''') = randV3 g''
        sumsqs = a1^2 + a2^2 + a3^2
        x = abs a0 * sqrt sumsqs
        invMag = recip.sqrt $ a0^2 + sumsqs + x + x
        mag = scale * invMag
    in (PV (mag * a0) (mag * a1) (mag * a2) (mag * a3), g''')
  
  
  -- | 'randH' random Quaternion made from random Grade 0 and Grade 2 elements
  randH :: RandomGen g => g -> (Cl3, g)
  randH = rangeH (0,1)
  
  
  -- | 'rangeH' random Quaternion made from random Grade 0 and Grade 2 elements within a range
  rangeH :: RandomGen g => (Cl3, Cl3) -> g -> (Cl3, g)
  rangeH (lo, hi) g =
    let (R scale, g') = rangeR (lo, hi) g
        (R a0, g'') = randR g'
        (BV a23 a31 a12, g''') = randBV g''
        invMag = recip.sqrt $ a0^2 + a23^2 + a31^2 + a12^2
        mag = scale * invMag
    in (H (mag * a0) (mag * a23) (mag * a31) (mag * a12), g''')
  
  
  -- | 'randC' random combination of Grade 0 and Grade 3
  randC :: RandomGen g => g -> (Cl3, g)
  randC = rangeC (0,1)
  
  
  -- | 'rangeC' random combination of Grade 0 and Grade 3 within a range
  rangeC :: RandomGen g => (Cl3, Cl3) -> g -> (Cl3, g)
  rangeC (lo, hi) g =
    let (R scale, g') = rangeR (lo, hi) g
        (phi, g'') = randomR (0, 2*pi) g'
    in (C (scale * cos phi) (scale * sin phi), g'')
  
  
  -- | 'randBPV' random combination of Grade 1 and Grade 2
  randBPV :: RandomGen g => g -> (Cl3, g)
  randBPV = rangeBPV (0,1)
  
  
  -- | 'rangeBPV' random combination of Grade 1 and Grade 2 within a range
  rangeBPV :: RandomGen g => (Cl3, Cl3) -> g -> (Cl3, g)
  rangeBPV (lo, hi) g =
    let (R scale, g') = rangeR (lo, hi) g
        (V3 a1 a2 a3, g'') = randV3 g'
        (BV a23 a31 a12, g''') = randBV g''
        x = sqrt $ (a1*a31 - a2*a23)^2 + (a1*a12 - a3*a23)^2 + (a2*a12 - a3*a31)^2
        invMag = recip.sqrt $ a1^2 + a23^2 + a2^2 + a31^2 + a3^2 + a12^2 + x + x
        mag = scale * invMag
    in (BPV (mag * a1) (mag * a2) (mag * a3) (mag * a23) (mag * a31) (mag * a12), g''')
  
  
  -- | 'randODD' random combination of Grade 1 and Grade 3
  randODD :: RandomGen g => g -> (Cl3, g)
  randODD = rangeODD (0,1)
  
  
  -- | 'rangeODD' random combination of Grade 1 and Grade 3 within a range
  rangeODD :: RandomGen g => (Cl3, Cl3) -> g -> (Cl3, g)
  rangeODD (lo, hi) g =
    let (R scale, g') = rangeR (lo, hi) g
        (V3 a1 a2 a3, g'') = randV3 g'
        (I a123, g''') = randI g''
        invMag = recip.sqrt $ a1^2 + a2^2 + a3^2 + a123^2
        mag = scale * invMag
    in (ODD (mag * a1) (mag * a2) (mag * a3) (mag * a123), g''')
  
  
  -- | 'randTPV' random combination of Grade 2 and Grade 3
  randTPV :: RandomGen g => g -> (Cl3, g)
  randTPV = rangeTPV (0,1)
  
  
  -- | 'rangeTPV' random combination of Grade 2 and Grade 3 within a range
  rangeTPV :: RandomGen g => (Cl3, Cl3) -> g -> (Cl3, g)
  rangeTPV (lo, hi) g =
    let (R scale, g') = rangeR (lo, hi) g
        (BV a23 a31 a12, g'') = randBV g'
        (I a123, g''') = randI g''
        sumsqs = a23^2 + a31^2 + a12^2
        x = abs a123 * sqrt sumsqs
        invMag = recip.sqrt $ sumsqs + a123^2 + x + x
        mag = scale * invMag
    in (TPV (mag * a23) (mag * a31) (mag * a12) (mag * a123), g''')
  
  
  -- | 'randAPS' random combination of all 4 grades
  randAPS :: RandomGen g => g -> (Cl3, g)
  randAPS = rangeAPS (0,1)
  
  
  -- | 'rangeAPS' random combination of all 4 grades within a range
  rangeAPS :: RandomGen g => (Cl3, Cl3) -> g -> (Cl3, g)
  rangeAPS (lo, hi) g =
    let (R scale, g') = rangeR (lo, hi) g
        (C a0 a123, g'') = randC g'
        (V3 a1 a2 a3, g''') = randV3 g''
        (BV a23 a31 a12, g'v) = randBV g'''
        x = sqrt $ (a0*a1 + a123*a23)^2 + (a0*a2 + a123*a31)^2 + (a0*a3 + a123*a12)^2 + (a2*a12 - a3*a31)^2 + (a3*a23 - a1*a12)^2 + (a1*a31 - a2*a23)^2
        invMag = recip.sqrt $ a0^2 + a1^2 + a2^2 + a3^2 + a23^2 + a31^2 + a12^2 + a123^2 + x + x
        mag = scale * invMag
    in (APS (mag * a0) (mag * a1) (mag * a2) (mag * a3) (mag * a23) (mag * a31) (mag * a12) (mag * a123), g'v)
  
  
  -------------------------------------------------------------------
  -- Additional Random generators
  -------------------------------------------------------------------
  -- | 'randUnitV3' a unit vector with a random direction
  randUnitV3 :: RandomGen g => g -> (Cl3, g)
  randUnitV3 g =
    let (theta, g') = randomR (0,2*pi) g
        (u, g'') = randomR (-1,1) g'
        semicircle = sqrt (1-u^2)
    in (V3 (semicircle * cos theta) (semicircle * sin theta) u, g'')
  
  
  -- | 'randProjector' a projector with a random direction
  randProjector :: RandomGen g => g -> (Cl3, g)
  randProjector g =
    let (V3 a1 a2 a3, g') = randUnitV3 g
    in (PV 0.5 (0.5 * a1) (0.5 * a2) (0.5 * a3), g')
  
  
  -- | 'rangeProjector' a projector with a range of random magnitudes and directions
  rangeProjector :: RandomGen g => (Cl3, Cl3) -> g -> (Cl3, g)
  rangeProjector (lo, hi) g =
    let (R mag, g') = rangeR (lo, hi) g
        (PV a0 a1 a2 a3, g'') = randProjector g'
    in (PV (mag * a0) (mag * a1) (mag * a2) (mag * a3), g'')
  
  
  -- | 'randNilpotent' a nilpotent element with a random orientation
  randNilpotent :: RandomGen g => g -> (Cl3, g)
  randNilpotent g =
    let (PV a0 a1 a2 a3, g') = randProjector g
        (V3 b1 b2 b3, g'') = randUnitV3 g'
        c1 = a2*b3 - a3*b2
        c2 = a3*b1 - a1*b3
        c3 = a1*b2 - a2*b1 -- (V3 c1 c2 c3) vector normal to the projector: mIx.toBV $ toV3 p * v
        invMag = recip.sqrt $ c1^2 + c2^2 + c3^2
        d1 = invMag * c1
        d2 = invMag * c2
        d3 = invMag * c3  -- (V3 d1 d2 d3) unit vector normal to the projector
    in (BPV (d1*a0) (d2*a0) (d3*a0) (d2*a3 - d3*a2) (d3*a1 - d1*a3) (d1*a2 - d2*a1), g'')
  
  
  -- | 'rangeNilpotent' a nilpotent with a range of random magnitudes and orientations
  rangeNilpotent :: RandomGen g => (Cl3, Cl3) -> g -> (Cl3, g)
  rangeNilpotent (lo, hi)  g =
    let (R mag, g') = rangeR (lo, hi) g
        (BPV a1 a2 a3 a23 a31 a12, g'') = randNilpotent g'
    in (BPV (mag * a1) (mag * a2) (mag * a3) (mag * a23) (mag * a31) (mag * a12), g'')
  
  
  -- | 'randUnitary' a unitary element with a random orientation
  randUnitary :: RandomGen g => g -> (Cl3, g)
  randUnitary g =
    let (tpv,g') = randTPV g
    in (exp tpv,g')
  
  
  -- | 'rangeUnitary' a unitary element with a range of random magnitudes and orientations, the exponential of a triparavector
  rangeUnitary :: RandomGen g => (Cl3, Cl3) -> g -> (Cl3, g)
  rangeUnitary (lo, hi) g =
    let (tpv, g') = rangeTPV (lo, hi) g
    in (exp tpv, g')
  
  
  -------------------------------------------------------------------
  -- helper functions
  -------------------------------------------------------------------
  magHelper :: RandomGen g => (Cl3, Cl3) -> g -> (Double, g)
  magHelper (lo, hi) g =
    let R lo' = abs lo
        R hi' = abs hi
    in randomR (lo', hi') g
  
  
  scalarHelper :: RandomGen g => (Double -> Cl3) -> (Cl3, Cl3) -> g -> (Cl3, g)
  scalarHelper con rng g =
    let (mag, g') = magHelper rng g
        (sign, g'') = random g'
    in if sign
       then (con mag, g'')
       else (con (negate mag), g'')
  
  
  vectorHelper :: RandomGen g => (Double -> Double -> Double -> Cl3) -> (Cl3, Cl3) -> g -> (Cl3, g)
  vectorHelper con rng g =
    let (mag, g') = magHelper rng g
        (V3 x y z, g'') = randUnitV3 g'
    in (con (mag * x) (mag * y) (mag * z), g'')
  
  
  #endif
  
  -- End of File
  