{-# LANGUAGE Safe #-} {-# LANGUAGE GADTSyntax #-} {-# LANGUAGE ViewPatterns #-} {-# LANGUAGE DeriveDataTypeable #-} {-# LANGUAGE DeriveGeneric #-} {-# OPTIONS_GHC -fno-warn-type-defaults #-} {-# LANGUAGE BangPatterns #-} {-# LANGUAGE CPP #-} {-# LANGUAGE MultiWayIf #-} #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-2020 Nathan Waivio -- License : BSD3 -- Maintainer : Nathan Waivio -- Stability : Stable -- Portability : unportable -- -- Library implementing standard functions for the 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_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 ) where #ifndef O_NO_DERIVED import Data.Data (Typeable, Data) import GHC.Generics (Generic) #endif import Control.DeepSeq (NFData,rnf) import Foreign.Storable (Storable, sizeOf, alignment, peek, poke) import Foreign.Ptr (Ptr, plusPtr, castPtr) #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 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 (R a0) = R (abs a0) -- absolute value of a real number abs (V3 a1 a2 a3) = R (sqrt (a1^2 + a2^2 + a3^2)) -- magnitude of a vector abs (BV a23 a31 a12) = R (sqrt (a23^2 + a31^2 + a12^2)) -- magnitude of a bivector abs (I a123) = R (abs a123) -- magnitude of a Imaginary number abs (PV a0 a1 a2 a3) = R (reimMag a0 a1 a2 a3) abs (TPV a23 a31 a12 a123) = R (reimMag a123 a23 a31 a12) abs (H a0 a23 a31 a12) = R (sqrt (a0^2 + a23^2 + a31^2 + a12^2)) -- largest singular value abs (C a0 a123) = R (sqrt (a0^2 + a123^2)) -- magnitude of a complex number abs (BPV a1 a2 a3 a23 a31 a12) = let x = 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 in R (sqrt (a1^2 + a23^2 + a2^2 + a31^2 + a3^2 + a12^2 + x + x)) abs (ODD a1 a2 a3 a123) = R (sqrt (a1^2 + a2^2 + a3^2 + a123^2)) abs (APS a0 a1 a2 a3 a23 a31 a12 a123) = let 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) -- core was duplicating this computation added let to hopefully reduce the duplication in R (sqrt (a0^2 + a1^2 + a2^2 + a3^2 + a23^2 + a31^2 + a12^2 + a123^2 + x + x)) -- |'signum' satisfies the Law "abs x * signum x == x" -- kind of cool: signum of a vector is it's unit vector. signum (R a0) = R (signum a0) signum (V3 a1 a2 a3) = let mag = sqrt (a1^2 + a2^2 + a3^2) invMag = recip mag in if mag == 0 then R 0 else V3 (invMag * a1) (invMag * a2) (invMag * a3) signum (BV a23 a31 a12) = let mag = sqrt (a23^2 + a31^2 + a12^2) invMag = recip mag in if mag == 0 then R 0 else BV (invMag * a23) (invMag * a31) (invMag * a12) signum (I a123) = I (signum a123) signum (PV a0 a1 a2 a3) = let mag = reimMag a0 a1 a2 a3 invMag = recip mag in if mag == 0 then R 0 else PV (invMag * a0) (invMag * a1) (invMag * a2) (invMag * a3) signum (H a0 a23 a31 a12) = let mag = sqrt (a0^2 + a23^2 + a31^2 + a12^2) invMag = recip mag in if mag == 0 then R 0 else H (invMag * a0) (invMag * a23) (invMag * a31) (invMag * a12) signum (C a0 a123) = let mag = sqrt (a0^2 + a123^2) invMag = recip mag in if mag == 0 then R 0 else C (invMag * a0) (invMag * a123) signum (BPV a1 a2 a3 a23 a31 a12) = let x = sqrt ((a1*a31 - a2*a23)^2 + (a1*a12 - a3*a23)^2 + (a2*a12 - a3*a31)^2) mag = sqrt (a1^2 + a23^2 + a2^2 + a31^2 + a3^2 + a12^2 + x + x) invMag = recip mag in if mag == 0 then R 0 else BPV (invMag * a1) (invMag * a2) (invMag * a3) (invMag * a23) (invMag * a31) (invMag * a12) signum (ODD a1 a2 a3 a123) = let mag = sqrt (a1^2 + a2^2 + a3^2 + a123^2) invMag = recip mag in if mag == 0 then R 0 else ODD (invMag * a1) (invMag * a2) (invMag * a3) (invMag * a123) signum (TPV a23 a31 a12 a123) = let mag = reimMag a123 a23 a31 a12 invMag = recip mag in if mag == 0 then R 0 else TPV (invMag * a23) (invMag * a31) (invMag * a12) (invMag * a123) signum (APS a0 a1 a2 a3 a23 a31 a12 a123) = let 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) mag = sqrt (a0^2 + a1^2 + a2^2 + a3^2 + a23^2 + a31^2 + a12^2 + a123^2 + x + x) invMag = recip mag in if mag == 0 then R 0 else APS (invMag * a0) (invMag * a1) (invMag * a2) (invMag * a3) (invMag * a23) (invMag * a31) (invMag * a12) (invMag * a123) -- |'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 + x + 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 (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 (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 -- | 'tol' currently 128*eps tol :: Cl3 {-# INLINE tol #-} tol = R 1.4210854715202004e-14 tol' :: Double {-# INLINE tol' #-} tol' = 1.4210854715202004e-14 -- | '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 ------------------------------------------------------------------- -- -- 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. 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" #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' simicircle = sqrt (1-u^2) in (V3 (simicircle * cos theta) (simicircle * 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