diff --git a/ChangeLog.md b/ChangeLog.md new file mode 100644 --- /dev/null +++ b/ChangeLog.md @@ -0,0 +1,5 @@ +# Revision history for cl3 + +## 1.0.0.0 -- 2017-10-28 + +* First version. Released on an unsuspecting world. diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,30 @@ +Copyright (c) 2017-2018, Nathan Waivio + +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + * Redistributions in binary form must reproduce the above + copyright notice, this list of conditions and the following + disclaimer in the documentation and/or other materials provided + with the distribution. + + * Neither the name of Nathan Waivio nor the names of other + contributors may be used to endorse or promote products derived + from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/README.md b/README.md new file mode 100644 --- /dev/null +++ b/README.md @@ -0,0 +1,142 @@ + +```haskell +-- CCCCCCCCCC llllllllll 3333333333 +-- CCCCCCCCCCCCCCCCCC llllllllll 333333333333333333 +-- CCCCCCC CCCCCCCCCC llllllllll 3333333333 333333 +-- CCCCC CCC llllllllll 333 ## 33333 +-- CCCCCC llllllllll # 33333 +-- CCCCC llllllllll 33333 +-- CCCCC llllllllll / @ 333333 +-- CCCCC llllllllll /| + 33333 +-- CCCCC lllllllllll \| + 33333 +-- CCCCC lllllllllllll \ @ 333333 +-- CCCCC lllllll lllllll 33333 +-- CCCCCC lllllll lllllll # 33333 +-- CCCCC CCC lllllll lllllll 333 ## 33333 +-- CCCCCCC CCCCCCCCCC lllllll lllllll 3333333333 333333 +-- CCCCCCCCCCCCCCCCCC lllllll lllllll 333333333333333333 +-- CCCCCCCCCC lllllll lllllll 3333333333 +``` + + + +# Cl3 +Cl3 is a Haskell Library implementing standard functions for the [Algebra of Physical Space](https://en.wikipedia.org/wiki/Algebra_of_physical_space) Cl(3,0) + +The goal of the Cl3 library is to provide a specialized, safe, high performance, Algebra of Physical Space implementation. +This library is suitable for physics simulations. The library integrates into Haskell's standard prelude and has few dependencies. +The library uses a ADT data type to specialize to specific graded elements in the Algebra of Physical Space. + + +# ADT Interface +The constructors are specialized to single and double grade combinations and the general case of APS. +Using the specialized constructors helps the compiler to compile to code similar to that you would hand write. +The constructors follow the following conventions for basis. + +```haskell +scalar = R e0 +vector = V3 e1 e2 e3 +bivector = BV e23 e31 e12 +trivectorPseudoScalar = I e123 +paravector = PV e0 e1 e2 e3 +quarternion = H e0 e23 e31 e12 +complex = C e0 e123 +biparavector = BPV e1 e2 e3 e23 e31 e12 +oddparavector = ODD e1 e2 e3 e123 +triparavector = TPV e23 e31 e12 e123 +aps = APS e0 e1 e2 e3 e23 e31 e12 e123 +``` +# Usage +In MATLAB or Octave one can write: `sqrt(-25)` and get `5.0i` + +In standard Haskell `sqrt (-25)` will produce `NaN` + +But using the Cl3 library `sqrt (-25) :: Cl3` will produce `I 5.0`, and likewise `(I 5.0)^2` will produce `R (-25)` + +If the unit imaginary is defined as `i = I 1`, expressions very similar to MATLAB can be formed `1.2 + 2.3*i` will produce `C 1.2 2.3` + +Vector addition is also natural, two arbitrary vectors `v1 = V3 a1 a2 a3` and `v2 = V3 b1 b2 b3` can be added `v1 + v2` and scaled `2*(v1+v2)` + +The dot product (inner product) of two arbitrary vectors is `toR $ v1 * v2`, that is the scalar part of the geometric product of two vectors. + +The cross product is the Hodge Dual of the wedge product (outer product) `-i * toBV (v1*v2)` + +The multiplication of two unit vectors is related to the rotor rotating from `u_from` to `u_to` like so `rot = sqrt $ u_to * u_from` + +Any arbitrary vector can be rotated by a rotor with the equation of `v' = rot * v * dag rot` + +Rotors can also be formed with an axis unit vector `u` and real scalar angle `theta` in units of radians, it produces the versor (unit quaternion) `rot = exp $ (-i/2) * theta * u` + +For special relativity with the velocity vector `v` and speed of light scalar `c`: +* Beta is `beta = v / c` +* Rapidity is `rapidity = atanh beta` +* Gamma is `gamma = cosh rapidity` +* Composition of velocities is simply adding the two rapidities and converting back to velocity +* Proper Velocity is `w = c * sinh rapidity` or `w = gamma * v` +* Four Velocity is a paravector `u = exp rapidity` where the real scalar part is `gamma * c` and the vector part is `w / c` +* The Boost is `boost = exp $ rapidity / 2` + +# APS Basis +Where __e0__ is the scalar basis frequently refered to as "1", in other texts. + +__e1__, __e2__, and __e3__ are the vector basis of 3 orthonormal vectors. + +__e23__, __e31__, and __e12__ are the bivector basis, these are formed by the outer product of two vector basis. For instance in the case of __e23__, the outer product, or wedge product, is __e2__ /\ __e3__, but because this can be simplified to the geometric product of __e2__ * __e3__ because the scalar part is zero for orthoginal vector basis'. The geometric product of the two basis vectors is further shortened for brevity to __e23__. + +__e123__ is the trivector basis, and is formed by the wedge product of __e1__ /\ __e2__ /\ __e3__, and likewise shortened to __e123__ + + +# Multiplication of the basis elements +The basis vectors multiply with the following multiplication table: + +| Mult | e0 | e1 | e2 | e3 | e23 | e31 | e12 | e123 | +|:--------:|:----:|:----:|:----:|:----:|:----:|:----:|:----:|:----:| +| __e0__ | e0 | e1 | e2 | e3 | e23 | e31 | e12 | e123 | +| __e1__ | e1 | e0 | e12 | -e31 | e123 | -e3 | e2 | e23 | +| __e2__ | e2 | -e12 | e0 | e23 | e3 | e123 | -e1 | e31 | +| __e3__ | e3 | e31 | -e23 | e0 | -e2 | e1 | e123 | e12 | +| __e23__ | e23 | e123 | -e3 | e2 | -e0 | -e12 | e31 | -e1 | +| __e31__ | e31 | e3 | e123 | -e1 | e12 | -e0 | -e23 | -e2 | +| __e12__ | e12 | -e2 | e1 | e123 | -e31 | e23 | -e0 | -e3 | +| __e123__ | e123 | e23 | e31 | e12 | -e1 | -e2 | -e3 | -e0 | + + +# Multiplication of the ADT Constructors +The grade specialized type constructors multiply with the following multiplication table: + +| Mult | R | V3 | BV | I | PV | H | C | BPV | ODD | TPV | APS | +|:-------:|:---:|:---:|:---:|:---:|:---:|:---:|:---:|:---:|:---:|:---:|:---:| +| __R__ | R | V3 | BV | I | PV | H | C | BPV | ODD | TPV | APS | +| __V3__ | V3 | H | ODD | BV | APS | ODD | BPV | APS | ODD | APS | APS | +| __BV__ | BV | ODD | H | V3 | APS | H | BPV | APS | ODD | APS | APS | +| __I__ | I | BV | V3 | R | TPV | ODD | C | BPV | H | PV | APS | +| __PV__ | PV | APS | APS | TPV | APS | APS | APS | APS | APS | APS | APS | +| __H__ | H | ODD | H | ODD | APS | H | APS | APS | ODD | APS | APS | +| __C__ | C | BPV | BPV | C | APS | APS | C | BPV | APS | APS | APS | +| __BPV__ | BPV | APS | ODD | BPV | APS | APS | BPV | APS | APS | APS | APS | +| __ODD__ | ODD | ODD | TPV | H | APS | ODD | APS | APS | H | APS | APS | +| __TPV__ | TPV | APS | APS | PV | APS | APS | APS | APS | APS | APS | APS | +| __APS__ | APS | APS | APS | APS | APS | APS | APS | APS | APS | APS | APS | + + +# Performace Benchmarking +A benchmark has been developed based on the Haskell entry for the N-Body Benchmark in the [The Computer Language Benchmarks Game](https://benchmarksgame-team.pages.debian.net/benchmarksgame/) with some modifications to run with [Criterion](https://hackage.haskell.org/package/criterion). +On my machine the current fastest implementation completes 50M steps with a mean time of 13.35 seconds. The benchmark uses a hand rolled implementation of vector math. The Cl3 implementation completes 50M steps with a mean time of 15.73 seconds. This 2.38 second difference amounts to a 23.8 ns difference in the inner loop. + +# Design Philosophy +The design space for Clifford Algebra libraries was explored quite a bit before the development of this library. Initially the isomorphism of APS with 2x2 Complex Matrices was used, this had the draw back that multiplying the scalar 2 * 2 would incur all of the computational cost of multiplying two 2x2 complex matrices. +Then the design was changed to lists that contained the basis' values, but lists are computationally slow and do not produce well optimized code. +Then a single constructor data type for APS was developed, but this had all of the drawbacks of 2x2 complex matrices. +The specialized ADT Constructor version of the library was developed and it showed that it had some promise. +More of the design space was explored, a version of the Cl3 library was developed using Multi-parameter Type Classes and Functional Dependencies, this didn't appear to have much gained over the specialized ADT Syntax interface and it didn't use the standard Prelude classes like Num, Float, etc. It was also difficult for me to figure out how to code a `reduce` function. +So the specialized ADT Constructor design of the Cl3 library was finished and released. + +# How does this fit in with the existing Haskell ecosystem? +Cl3 is meant to be a [Linear](https://hackage.haskell.org/package/linear) killer based on Geometric Algebra. The linear package +consists of many different types that are not easily combinable using the Num Class, and require many specialized functions each to multiply a different combination of types. + +The [clifford](https://hackage.haskell.org/package/clifford) package uses the [Numeric Prelude](https://hackage.haskell.org/package/numeric-prelude), +for a Clifford Algebra of arbitrary signature that stores multivector blades in a list data structure. + +The [clif](https://hackage.haskell.org/package/clif) is for symbolic computing using symbolic and numeric computations +with finite and infinite-dimensional Clifford algebras arising from arbitrary bilinear forms. The libraries representation of a Cliffor also makes use of lists. diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/benchmarks/NbodyGameCl3.hs b/benchmarks/NbodyGameCl3.hs new file mode 100644 --- /dev/null +++ b/benchmarks/NbodyGameCl3.hs @@ -0,0 +1,740 @@ +{-# LANGUAGE ViewPatterns #-} +{-# LANGUAGE BangPatterns #-} + +------------------------------------------------------------------ +-- | +-- Description: Benchmark for Algebra.Geometric.Cl3 +-- Uses standard "The Computer language Benchmarks Game" implementation +-- with Criterion to compare the performance difference between +-- using the Cl3 library vs. the benchmarks implementation using doubles. +-- +-- +-- +-- This derivative work was derived from the excellent work of: +-- Branimir Maksimovic. +-- +------------------------------------------------------------------- + +module Main (main) where + +import Criterion.Main (defaultMain, bgroup, bench, nfIO) -- To add Criterion to the benchmark +import Algebra.Geometric.Cl3 (Cl3(..), toR, toV3) -- To add Cl3 + + +import Foreign.Ptr (Ptr, castPtr, plusPtr) +import Foreign.Storable (Storable(..)) +import Foreign.Marshal.Alloc (mallocBytes) +import Control.Monad (when) +import Text.Printf (printf) + +------------------------------------------------------------------- +-- | +-- An n-body solver for the Sun and gas giants using the semi-implicit Euler method. +-- The benchmark uses the Criterion tool to measure timing. +-- The benchmark advances 50 million steps. +-- +-- The @nbodyBaseline@ is derived from Branimir Maksimovic's work with some changes suggested by @hlint@. +-- @nbodyCl3@ aggressively uses the "to" functions like 'toR' and 'toV3' to take advantage of the case of +-- known constructor optimization. @nbodyAPS@ does not make aggressive use "to" functions. The +-- program would not finish compiling (memory would grow to 54GB while compiling) before INLINE pragmas +-- were added to the Cl3 library. +-- +-- * Baseline 20171007: +-- +-- Initial Commit of the Code, GHC 8.0.2 +-- +-- >>> ./benchmarks/NbodyGameCl3 +-- benchmarking nbodyBaseline/50000000 +-- time 37.99 s (37.92 s .. 38.08 s) +-- 1.000 R² (1.000 R² .. 1.000 R²) +-- mean 37.93 s (37.89 s .. 37.95 s) +-- std dev 36.25 ms (0.0 s .. 38.29 ms) +-- variance introduced by outliers: 19% (moderately inflated) +-- benchmarking nbodyCl3/50000000 +-- time 38.95 s (38.91 s .. 39.03 s) +-- 1.000 R² (1.000 R² .. 1.000 R²) +-- mean 38.95 s (38.94 s .. 38.96 s) +-- std dev 14.70 ms (0.0 s .. 16.66 ms) +-- variance introduced by outliers: 19% (moderately inflated) +-- benchmarking nbodyAPS/50000000 +-- time 179.4 s (177.4 s .. 181.4 s) +-- 1.000 R² (1.000 R² .. 1.000 R²) +-- mean 179.3 s (178.8 s .. 179.5 s) +-- std dev 395.0 ms (0.0 s .. 413.9 ms) +-- variance introduced by outliers: 19% (moderately inflated) +-- +-- * 20171016: +-- +-- Change of INLINE phase control to fix ghc simplifier ticks issue, GHC 8.0.2 +-- +-- >>> ./benchmarks/NbodyGameCl3 +-- benchmarking nbodyBaseline/50000000 +-- time 37.21 s (36.80 s .. 37.46 s) +-- 1.000 R² (1.000 R² .. 1.000 R²) +-- mean 37.02 s (36.89 s .. 37.10 s) +-- std dev 125.4 ms (0.0 s .. 144.2 ms) +-- variance introduced by outliers: 19% (moderately inflated) +-- benchmarking nbodyCl3/50000000 +-- time 38.57 s (38.07 s .. 39.20 s) +-- 1.000 R² (1.000 R² .. 1.000 R²) +-- mean 38.55 s (38.44 s .. 38.64 s) +-- std dev 128.8 ms (0.0 s .. 147.5 ms) +-- variance introduced by outliers: 19% (moderately inflated) +-- benchmarking nbodyAPS/50000000 +-- time 176.9 s (169.3 s .. 181.7 s) +-- 1.000 R² (0.999 R² .. 1.000 R²) +-- mean 177.1 s (176.1 s .. 177.9 s) +-- std dev 1.180 s (0.0 s .. 1.316 s) +-- variance introduced by outliers: 19% (moderately inflated) +-- +-- * 20171130: +-- +-- Though, thru, some unknown reason, an extreme increase of performance has occurred. +-- Not sure what changed but a huge improvement, dev-lang/ghc-8.0.2, +-- upgraded sys-devel/gcc-6.4.0, sys-devel/llvm-3.9.1-r1, linux-4.12.12-gentoo? +-- Spooky thing is that the parent compiled executable from before that ran ~35.6 Sec +-- is now ~10.7 Sec. +-- +-- >>> ./benchmarks/NbodyGameCl3 +-- benchmarking nbodyBaseline/50000000 +-- time 14.04 s (13.67 s .. 14.60 s) +-- 1.000 R² (1.000 R² .. 1.000 R²) +-- mean 13.74 s (13.59 s .. 13.85 s) +-- std dev 168.6 ms (0.0 s .. 190.4 ms) +-- variance introduced by outliers: 19% (moderately inflated) +-- benchmarking nbodyCl3/50000000 +-- time 15.68 s (14.80 s .. 16.12 s) +-- 1.000 R² (0.999 R² .. 1.000 R²) +-- mean 15.79 s (15.64 s .. 15.88 s) +-- std dev 141.7 ms (0.0 s .. 161.5 ms) +-- variance introduced by outliers: 19% (moderately inflated) +-- benchmarking nbodyAPS/50000000 +-- time 182.9 s (182.4 s .. 183.5 s) +-- 1.000 R² (1.000 R² .. 1.000 R²) +-- mean 182.7 s (182.5 s .. 182.8 s) +-- std dev 162.1 ms (0.0 s .. 183.2 ms) +-- variance introduced by outliers: 19% (moderately inflated) +-- +-- 20180309 : +-- +-- Well, while correcting spelling I found several issues re: Ord; +-- Ord is now more what you'd expect (now -2 is smaller than -1). +-- It appears to have sped up the performance a bit. +-- +-- benchmarking nbodyBaseline/50000000 +-- time 13.47 s (12.95 s .. 13.89 s) +-- 1.000 R² (1.000 R² .. NaN R²) +-- mean 13.54 s (13.42 s .. 13.59 s) +-- std dev 98.11 ms (0.0 s .. 103.4 ms) +-- variance introduced by outliers: 19% (moderately inflated) +-- benchmarking nbodyCl3/50000000 +-- time 15.50 s (15.45 s .. 15.59 s) +-- 1.000 R² (1.000 R² .. 1.000 R²) +-- mean 15.54 s (15.52 s .. 15.56 s) +-- std dev 27.35 ms (0.0 s .. 30.00 ms) +-- variance introduced by outliers: 19% (moderately inflated) +-- benchmarking nbodyAPS/50000000 +-- time 154.5 s (149.7 s .. 159.5 s) +-- 1.000 R² (0.999 R² .. 1.000 R²) +-- mean 153.9 s (152.8 s .. 154.6 s) +-- std dev 1.119 s (0.0 s .. 1.283 s) +-- variance introduced by outliers: 19% (moderately inflated) +-- +-- 20180518: +-- +-- Submitted the library to review on Haskell :: Reddit, implemented +-- several suggestions, the main one that may affect performance +-- would be removing the orphan instances and moving the instances +-- to Cl3 and using CPP gating of the optional instances. +-- Performance seems to stay the same. +-- +-- benchmarking nbodyBaseline/50000000 +-- time 13.38 s (12.37 s .. 14.02 s) +-- 0.999 R² (0.998 R² .. 1.000 R²) +-- mean 13.35 s (13.23 s .. 13.52 s) +-- std dev 167.1 ms (45.82 ms .. 225.3 ms) +-- variance introduced by outliers: 19% (moderately inflated) +-- +-- benchmarking nbodyCl3/50000000 +-- time 15.67 s (15.58 s .. 15.72 s) +-- 1.000 R² (1.000 R² .. 1.000 R²) +-- mean 15.73 s (15.70 s .. 15.79 s) +-- std dev 60.02 ms (583.2 μs .. 70.85 ms) +-- variance introduced by outliers: 19% (moderately inflated) +-- +-- benchmarking nbodyAPS/50000000 +-- time 156.2 s (154.5 s .. 159.3 s) +-- 1.000 R² (1.000 R² .. 1.000 R²) +-- mean 155.1 s (154.8 s .. 155.7 s) +-- std dev 560.8 ms (88.83 ms .. 714.2 ms) +-- variance introduced by outliers: 19% (moderately inflated) +-- +------------------------------------------------------------------- + +main :: IO () +main = do n <- return (50000000 :: Int) + defaultMain [bgroup "nbodyBaseline" [ bench (show n) $ nfIO (nbodyBaseline n)], + bgroup "nbodyCl3" [ bench (show n) $ nfIO (nbodyCl3 n)], + bgroup "nbodyAPS" [ bench (show n) $ nfIO (nbodyAPS n)] + ] + + + +nbodyBaseline :: (Num t, Eq t) => t -> IO () +nbodyBaseline n = do + pPlanets <- fromList planets -- load planets into a stack + nbodyInit pPlanets -- initialize the system momentum + energy pPlanets >>= printf "%.9f\n" -- calculate and print out the initial energy of the system + run n pPlanets -- solve the initial value problem, increment the system 'n' times + energy pPlanets >>= printf "%.9f\n" -- calculate and print out the final energy of the system + + +nbodyCl3 :: (Eq t, Num t) => t -> IO () +nbodyCl3 n = do + pPlanetsCl3 <- fromListCl3 planetsCl3 + nbodyInitCl3 pPlanetsCl3 + energyCl3 pPlanetsCl3 >>= print.toR -- printf "%0.9f\n" + runCl3 n pPlanetsCl3 + energyCl3 pPlanetsCl3 >>= print.toR -- printf "%0.9f\n" + + +nbodyAPS :: (Eq t, Num t) => t -> IO () +nbodyAPS n = do + pPlanetsAPS <- fromListCl3 planetsCl3 + nbodyInitAPS pPlanetsAPS + energyAPS pPlanetsAPS >>= print.toR -- printf "%0.9f\n" + runAPS n pPlanetsAPS + energyAPS pPlanetsAPS >>= print.toR -- printf "%0.9f\n" + + +-- | 'run' advances the system 'n' times by recursively calling 'run' +run :: (Eq t, Num t) => t -> Ptr Planet -> IO () +run 0 _ = return () +run i p = do + advance p + run (i-1) p + + +runCl3 :: (Eq t, Num t) => t -> Ptr PlanetCl3 -> IO () +runCl3 0 _ = return () +runCl3 i p = do + advanceCl3 p + runCl3 (i-1) p + + +runAPS :: (Eq t, Num t) => t -> Ptr PlanetCl3 -> IO () +runAPS 0 _ = return () +runAPS i p = do + advanceAPS p + runAPS (i-1) p + + +-- | 'offsetMomentum' is used to update the sun's momentum +-- so that the total momentum will be zero for the whole +-- system at initialization. +offsetMomentum :: Planet -> (Double, Double, Double) -> Planet +offsetMomentum p (px,py,pz) = p {vx = -px / solarMass, + vy = -py / solarMass, + vz = -pz / solarMass + } + + +offsetMomentumCl3 :: PlanetCl3 -> Cl3 -> PlanetCl3 +offsetMomentumCl3 (PlanetCl3 (pos') _ (mass')) (momentum) = PlanetCl3 pos' (toV3.negate $ toV3 momentum * (toR.recip.toR $ mass')) mass' + + +offsetMomentumAPS :: PlanetCl3 -> Cl3 -> PlanetCl3 +offsetMomentumAPS p momCl3 = p{velCl3 = negate momCl3 / massCl3 p} + + + +-- | 'nbodyInit' initializes the velocity of the sun so that +-- the total system momentum is zero. +nbodyInit :: Ptr Planet -> IO () +nbodyInit pPlanets = do + let initBl (px,py,pz) i = -- do "Is this do do?" -- init cycles through the stack and builds up the total momentum + if i < length planets + then do + p <- peekElemOff pPlanets i + initBl (px + vx p * mass p, py + vy p * mass p, pz + vz p * mass p) (i+1) -- requires velocity and mass + else return (px,py,pz) + s <- initBl (0,0,0) 0 -- total momentum of the system + p <- peek pPlanets -- grab the sun, don't get burned + poke pPlanets $ offsetMomentum p s -- update the sun's momentum +-- end of nbodyInit + + +nbodyInitCl3 :: Ptr PlanetCl3 -> IO () +nbodyInitCl3 pPlanetsCl3 = do + let initCl3 (toV3 -> pCl3) i = + if i < length planetsCl3 + then do + (PlanetCl3 _ (toV3 -> vel') (toR -> mass')) <- peekElemOff pPlanetsCl3 i + initCl3 (toV3 $! pCl3 + vel' * mass') (i+1) + else return pCl3 + + totalMomentum <- initCl3 (V3 0 0 0) 0 + sun0 <- peek pPlanetsCl3 + poke pPlanetsCl3 $ offsetMomentumCl3 sun0 totalMomentum +-- end of nbodyInitCl3 + + +nbodyInitAPS :: Ptr PlanetCl3 -> IO () +nbodyInitAPS pPlanetsCl3 = do + let initCl3 !pCl3 i = + if i < length planetsCl3 + then do + p <- peekElemOff pPlanetsCl3 i + initCl3 (pCl3 + velCl3 p * massCl3 p) (i+1) + else return pCl3 + + totalMomentum <- initCl3 (APS 0 0 0 0 0 0 0 0) 0 + sun0 <- peek pPlanetsCl3 + poke pPlanetsCl3 $ offsetMomentumAPS sun0 totalMomentum +-- end of nbodyInitAPS + + + +-- | 'squared' scalar product of a vector with its self +squared :: Num a => a -> a -> a -> a +squared x' y' z' = x' * x' + y' * y' + z' * z' + + +-- | 'energy' calculate the total energy of the system +-- this is the gravitational potential energy added to +-- the kinetic energy of all of the planets. +energy :: Ptr Planet -> IO Double +energy pPlanets = do + let energy' e i = if i < length planets -- '' + then do + p <- peekElemOff pPlanets i + e1 <- energy'' p (i+1) e + e2 <- energy' e (i+1) + return $ e + 0.5 * mass p * squared (vx p) (vy p) (vz p) + e1 + e2 -- requires mass and velocity + else return e + energy'' p j e = if j < length planets + then do + pj <- peekElemOff pPlanets j + let distance = sqrt $ squared dx dy dz + dx = x pj - x p + dy = y pj - y p + dz = z pj - z p + e1 <- energy'' p (j+1) e + return $ e - (mass p * mass pj) / distance + e1 -- requires mass and position + else return e + energy' 0.0 0 -- starts off the recursive calculation of energy + + +energyCl3 :: Ptr PlanetCl3 -> IO Cl3 +energyCl3 pPlanetsCl3 = do + let energy' e i = + if i < length planetsCl3 + then do + p@(PlanetCl3 _ (toV3 -> vel1) (toR -> mass1)) <- peekElemOff pPlanetsCl3 i + let !ke = toR $ 0.5 * mass1 * toR (vel1^(2 :: Int)) + e1 <- energy'' p (i+1) e + e2 <- energy' e (i+1) + return $! e + ke + e1 + e2 + else return e + + energy'' p@(PlanetCl3 (toV3 -> pos1) _ (toR -> mass1)) j e = + if j < length planetsCl3 + then do + (PlanetCl3 (toV3 -> pos2) _ (toR -> mass2)) <- peekElemOff pPlanetsCl3 j + let !dpos = toV3 $ pos2 - pos1 + !distance = abs dpos + !pe = toR $ (mass1 * mass2) / distance + e1 <- energy'' p (j+1) e + return $! e - pe + e1 -- requires mass and position + else return e + + energy' (R 0) 0 -- starts off the recursive calculation of energy +-- end of energyCl3 + + +energyAPS :: Ptr PlanetCl3 -> IO Cl3 +energyAPS pPlanetsCl3 = do + let energy' e i = + if i < length planetsCl3 + then do + p <- peekElemOff pPlanetsCl3 i + e1 <- energy'' p (i+1) e + e2 <- energy' e (i+1) + return $! e + 0.5 * massCl3 p * (velCl3 p)^(2 :: Int) + e1 + e2 + else return e + + energy'' p j e = + if j < length planetsCl3 + then do + pj <- peekElemOff pPlanetsCl3 j + let !distance = abs dpos + !dpos = posCl3 pj - posCl3 p + e1 <- energy'' p (j+1) e + return $! e - (massCl3 p * massCl3 pj) / distance + e1 + else return e + + energy' (APS 0 0 0 0 0 0 0 0) 0 -- starts off the recursive calculation of energy +-- end of energyAPS + + + +-- | 'advance' integrates the system of differential equations using a +-- semi-implicit Euler method, also called symplectic Euler method. +-- The first order method updates the velocities based on the acceleration +-- calculation, then it updates the positions. The velocity verlet +-- would be an improvement in accuracy because it is a second order +-- symplectic integrator. +advance :: Ptr Planet -> IO () +advance pPlanets = do + let advance' i = when (i < length planets) $ do -- loops through all the planets and updates the velocity + let loop j = when (j < length planets) $ do + ii <- peekElemOff pPlanets i + jj <- peekElemOff pPlanets j + let mag = dt / (dSquared * sqrt dSquared) + dSquared = squared dx dy dz + dx = x ii - x jj + dy = y ii - y jj + dz = z ii - z jj + pokeV pPlanets i ii{ -- requires position and mass to update the velocity + vx = vx ii - dx * mass jj * mag, + vy = vy ii - dy * mass jj * mag, + vz = vz ii - dz * mass jj * mag + } + pokeV pPlanets j jj{ -- requires position and mass to update the velocity + vx = vx jj + dx * mass ii * mag, + vy = vy jj + dy * mass ii * mag, + vz = vz jj + dz * mass ii * mag + } + loop (j+1) + loop (i+1) + advance' (i+1) + advance'' i = when (i < length planets) $ do -- loops through all of the planets and updates the position + p <- peekElemOff pPlanets i + pokeC pPlanets i p{ -- requires position and velocity to update the position + x = x p + dt * vx p, + y = y p + dt * vy p, + z = z p + dt * vz p + } + advance'' (i+1) + advance' 0 -- update all of the planets velocities + advance'' 0 -- update all of the planets positions +-- end of advance + + +advanceCl3 :: Ptr PlanetCl3 -> IO () +advanceCl3 pPlanetsCl3 = do + let advance' i = when (i < length planetsCl3) $ do -- loops through all the planets and updates the velocity + let loop j = when (j < length planetsCl3) $ do + (PlanetCl3 (toV3 -> posi) (toV3 -> veli) (toR -> massi)) <- peekElemOff pPlanetsCl3 i + (PlanetCl3 (toV3 -> posj) (toV3 -> velj) (toR -> massj)) <- peekElemOff pPlanetsCl3 j + let !dpos = posi - posj + !(R dMag) = abs dpos + !mag = R (dt / (dMag^(3 :: Int))) + !veli' = veli - dpos * massj * mag + !velj' = velj + dpos * massi * mag + pokeVCl3 pPlanetsCl3 i (PlanetCl3 posi veli' massi) + pokeVCl3 pPlanetsCl3 j (PlanetCl3 posj velj' massj) + loop (j+1) + loop (i+1) + advance' (i+1) + + advance'' i = when (i < length planetsCl3) $ do -- loops through all of the planets and updates the position + (PlanetCl3 (toV3 -> posi) (toV3 -> veli) massi) <- peekElemOff pPlanetsCl3 i + let !posi' = posi + R dt * veli + pokeCCl3 pPlanetsCl3 i (PlanetCl3 posi' veli massi) + advance'' (i+1) + + advance' 0 -- update all of the planets velocities + advance'' 0 -- update all of the planets positions +-- end of advanceCl3 + + +advanceAPS :: Ptr PlanetCl3 -> IO () +advanceAPS pPlanetsCl3 = do + let advance' i = when (i < length planetsCl3) $ do -- loops through all the planets and updates the velocity + let loop j = when (j < length planetsCl3) $ do + ii <- peekElemOff pPlanetsCl3 i + jj <- peekElemOff pPlanetsCl3 j + let !mag = R dt / (dSquared * abs dpos) + !dSquared = dpos * dpos + !dpos = posCl3 ii - posCl3 jj + pokeVCl3 pPlanetsCl3 i ii{velCl3 = velCl3 ii - dpos * massCl3 jj * mag} + pokeVCl3 pPlanetsCl3 j jj{velCl3 = velCl3 jj + dpos * massCl3 ii * mag} + loop (j+1) + loop (i+1) + advance' (i+1) + + advance'' i = when (i < length planetsCl3) $ do -- loops through all of the planets and updates the position + p <- peekElemOff pPlanetsCl3 i + pokeCCl3 pPlanetsCl3 i p{posCl3 = posCl3 p + R dt * velCl3 p} + advance'' (i+1) + + advance' 0 -- update all of the planets velocities + advance'' 0 -- update all of the planets positions +-- end of advanceAPS + + + + +data Planet = Planet {x :: !Double, + y :: !Double, + z :: !Double, + vx :: !Double, + vy :: !Double, + vz :: !Double, + mass :: !Double } deriving (Show) + +planets :: [Planet] +planets = [sun, jupiter, saturn, uranus, neptune] + + +-- | A planet can be described by three cliffors. +data PlanetCl3 = PlanetCl3 {posCl3 :: !Cl3, + velCl3 :: !Cl3, + massCl3 :: !Cl3} deriving (Show) + +planetsCl3 :: [PlanetCl3] +planetsCl3 = [PlanetCl3 sunPos sunVel sunMass, + PlanetCl3 jupiterPos jupiterVel jupiterMass, + PlanetCl3 saturnPos saturnVel saturnMass, + PlanetCl3 uranusPos uranusVel uranusMass, + PlanetCl3 neptunePos neptuneVel neptuneMass] + + + +sun :: Planet +sun = Planet {x = 0, y = 0, z = 0, + vx = 0, vy = 0, vz = 0, + mass = solarMass + } + + +sunPos :: Cl3 +sunPos = V3 0 0 0 +sunVel :: Cl3 +sunVel = V3 0 0 0 +sunMass :: Cl3 +sunMass = R solarMass + + +jupiter :: Planet +jupiter = Planet + {x = 4.84143144246472090e+00, y = -1.16032004402742839e+00, z= -1.03622044471123109e-01, + vx = 1.66007664274403694e-03*dp, vy = 7.69901118419740425e-03*dp, vz = -6.90460016972063023e-05*dp, + mass = 9.54791938424326609e-04 * solarMass + } + + +jupiterPos :: Cl3 +jupiterPos = V3 (4.84143144246472090e+00) (-1.16032004402742839e+00) (-1.03622044471123109e-01) +jupiterVel :: Cl3 +jupiterVel = V3 (1.66007664274403694e-03*dp) (7.69901118419740425e-03*dp) (-6.90460016972063023e-05*dp) +jupiterMass :: Cl3 +jupiterMass = R (9.54791938424326609e-04 * solarMass) + + +saturn :: Planet +saturn = Planet + {x = 8.34336671824457987e+00, y = 4.12479856412430479e+00, z = -4.03523417114321381e-01, + vx = -2.76742510726862411e-03*dp, vy = 4.99852801234917238e-03*dp, vz = 2.30417297573763929e-05*dp, + mass = 2.85885980666130812e-04 * solarMass + } + + +saturnPos :: Cl3 +saturnPos = V3 (8.34336671824457987e+00) (4.12479856412430479e+00) (-4.03523417114321381e-01) +saturnVel :: Cl3 +saturnVel = V3 (-2.76742510726862411e-03*dp) (4.99852801234917238e-03*dp) (2.30417297573763929e-05*dp) +saturnMass :: Cl3 +saturnMass = R (2.85885980666130812e-04 * solarMass) + + +uranus :: Planet +uranus = Planet + {x = 1.28943695621391310e+01,y = -1.51111514016986312e+01,z = -2.23307578892655734e-01, + vx = 2.96460137564761618e-03*dp,vy = 2.37847173959480950e-03*dp, vz = -2.96589568540237556e-05*dp, + mass = 4.36624404335156298e-05 * solarMass + } + + +uranusPos :: Cl3 +uranusPos = V3 (1.28943695621391310e+01) (-1.51111514016986312e+01) (-2.23307578892655734e-01) +uranusVel :: Cl3 +uranusVel = V3 (2.96460137564761618e-03*dp) (2.37847173959480950e-03*dp) (-2.96589568540237556e-05*dp) +uranusMass :: Cl3 +uranusMass = R (4.36624404335156298e-05 * solarMass) + + +neptune :: Planet +neptune = Planet + {x = 1.53796971148509165e+01,y = -2.59193146099879641e+01,z = 1.79258772950371181e-01, + vx = 2.68067772490389322e-03*dp,vy = 1.62824170038242295e-03*dp, vz = -9.51592254519715870e-05*dp, + mass = 5.15138902046611451e-05 * solarMass + } + + +neptunePos :: Cl3 +neptunePos = V3 (1.53796971148509165e+01) (-2.59193146099879641e+01) (1.79258772950371181e-01) +neptuneVel :: Cl3 +neptuneVel = V3 (2.68067772490389322e-03*dp) (1.62824170038242295e-03*dp) (-9.51592254519715870e-05*dp) +neptuneMass :: Cl3 +neptuneMass = R (5.15138902046611451e-05 * solarMass) + +daysPerYear :: Double +daysPerYear = 365.24 + +solarMass :: Double +solarMass = 4 * pi^(2 :: Int) +dp :: Double +dp = daysPerYear +dt :: Double +dt = 0.01 + + + +-- | 'Storable' instance of Planet +instance Storable Planet where + sizeOf _ = 8 * dblSz + alignment _ = dblSz + peekElemOff p i = peek (plusPtr p (i * sizeOf (undefined :: Planet))) + pokeElemOff p i e = poke (plusPtr p (i * sizeOf e)) e + peek p = do + x' <- peek (offset 0) + y' <- peek (offset 1) + z' <- peek (offset 2) + vx' <- peek (offset 3) + vy' <- peek (offset 4) + vz' <- peek (offset 5) + mass' <- peek (offset 6) + return Planet {x=x',y=y',z=z',vx=vx',vy=vy',vz=vz',mass=mass'} + where + offset i = plusPtr (castPtr p :: Ptr Double) (i*8) + poke p e = do + poke (offset 0) $ x e + poke (offset 1) $ y e + poke (offset 2) $ z e + poke (offset 3) $ vx e + poke (offset 4) $ vy e + poke (offset 5) $ vz e + poke (offset 6) $ mass e + where + offset i = plusPtr (castPtr p :: Ptr Double) (i*8) + +dblSz :: Int +dblSz = sizeOf (undefined :: Double) + + +-- | 'Storable' instance of PlanetCl3 +instance Storable PlanetCl3 where + sizeOf _ = 4 * cl3Sz + alignment _ = cl3Sz + peekElemOff p i = peek ( p `plusPtr` (i * sizeOf (undefined :: PlanetCl3))) + pokeElemOff p i e = poke (p `plusPtr` (i * sizeOf e)) e + peek p = do + pos' <- peek (offset 0) + vel' <- peek (offset 1) + mass' <- peek (offset 2) + return $ PlanetCl3 pos' vel' mass' + where + offset i = (castPtr p :: Ptr Cl3) `plusPtr` (i * cl3Sz) + poke p e = do + poke (offset 0) $! posCl3 e + poke (offset 1) $! velCl3 e + poke (offset 2) $! massCl3 e + where + offset i = (castPtr p :: Ptr Cl3) `plusPtr` (i * cl3Sz) + +cl3Sz :: Int +cl3Sz = sizeOf (undefined :: Cl3) + + +-- | 'pokeC' update the position +pokeC :: Ptr Planet -> Int -> Planet -> IO () +pokeC p i e = do + poke (offset 0) $ x e + poke (offset 1) $ y e + poke (offset 2) $ z e + where + offset o = (castPtr p :: Ptr Double) `plusPtr` (o*8+i*64) + + +pokeCCl3 :: Ptr PlanetCl3 -> Int -> PlanetCl3 -> IO () +pokeCCl3 p !i !e = + poke ((castPtr p :: Ptr Cl3) `plusPtr` (i * sizeOf (undefined :: PlanetCl3))) $! posCl3 e + + +-- | 'pokeV' update the velocity +pokeV :: Ptr Planet -> Int -> Planet -> IO () +pokeV p i e = do + poke (offset 3) $ vx e + poke (offset 4) $ vy e + poke (offset 5) $ vz e + where + offset o = (castPtr p::Ptr Double) `plusPtr` (o*8+i*64) + + +pokeVCl3 :: Ptr PlanetCl3 -> Int -> PlanetCl3 -> IO () +pokeVCl3 p !i !e = + poke ((castPtr p :: Ptr Cl3) `plusPtr` (cl3Sz + i * sizeOf (undefined :: PlanetCl3))) $! velCl3 e + + +-- | 'fromList' initialize the planets in the stack +fromList :: [Planet] -> IO (Ptr Planet) +fromList l = do + let len = length l + pa <- mallocBytes (len * sizeOf (undefined :: Planet)) + let loop [] _ = return () + loop (q:qs) i = do + poke (pa `plusPtr` (i * sizeOf (undefined :: Planet))) q + loop qs (i+1) + loop l 0 + return pa + + +fromListCl3 :: [PlanetCl3] -> IO (Ptr PlanetCl3) +fromListCl3 l = do + let len = length l + pa <- mallocBytes (len * sizeOf (undefined :: PlanetCl3)) + let loop [] _ = return () + loop (q:qs) i = do + poke (pa `plusPtr` (i * sizeOf (undefined :: PlanetCl3))) q + loop qs (i+1) + loop l 0 + return pa + + + + + + +------------------------------------------------------------------- +-- +-- The Computer Language +-- Benchmarks Game +-- Revised BSD license +-- +-- This is a specific instance of the Open Source Initiative (OSI) BSD license template. +-- +-- Copyright © 2004-2008 Brent Fulgham, 2005-2017 Isaac Gouy +-- +-- All rights reserved. +-- +-- Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: +-- +-- > Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. +-- +-- > Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. +-- +-- > Neither the name of "The Computer Language Benchmarks Game" nor the name of "The Computer Language Shootout Benchmarks" nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. +-- +-- THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +-- +-- +-- +-- The Computer Language Benchmarks Game +-- http://benchmarksgame.alioth.debian.org/ +-- +-- Contributed by Branimir Maksimovic +-- +------------------------------------------------------------------- + diff --git a/cl3.cabal b/cl3.cabal new file mode 100644 --- /dev/null +++ b/cl3.cabal @@ -0,0 +1,108 @@ +-- Initial cl3.cabal generated by cabal init. For further documentation, +-- see http://haskell.org/cabal/users-guide/ + +-- The name of the package. +name: cl3 + +-- The package version. See the Haskell package versioning policy (PVP) +-- for standards guiding when and how versions should be incremented. +-- https://wiki.haskell.org/Package_versioning_policy +-- PVP summary: +-+------- breaking API changes +-- | | +----- non-breaking API additions +-- | | | +--- code changes with no API change +version: 1.0.0.0 + +-- A short (one-line) description of the package. +synopsis: Clifford Algebra of three dimensional space. + +-- A longer description of the package. +description: Haskell Library implementing standard functions for the Algebra of Physical Space Cl(3,0) + +-- URL for the project homepage or repository. +homepage: https://github.com/waivio/cl3 + +bug-reports: https://github.com/waivio/cl3/issues + +-- The license under which the package is released. +license: BSD3 + +-- The file containing the license text. +license-file: LICENSE + +-- The package author(s). +author: Nathan Waivio + +-- An email address to which users can send suggestions, bug reports, and +-- patches. +maintainer: Nathan Waivio + +-- A copyright notice. +copyright: Copyright (C) 2017-2018 Nathan Waivio + +category: Math, Algebra + +build-type: Simple + +-- Extra files to be distributed with the package, such as examples or a +-- README. +extra-source-files: ChangeLog.md, README.md + +-- Constraint on the version of Cabal needed to build this package. +cabal-version: >=1.10 + +source-repository head + type: git + location: https://github.com/waivio/cl3.git + +library + -- Modules exported by the library. + exposed-modules: + Algebra.Geometric.Cl3, + Algebra.Geometric.Cl3.JonesCalculus + + -- Compiler options + ghc-options: -Wall -O2 + + -- LANGUAGE extensions used by modules in this package. + other-extensions: + Safe, + GADTSyntax, + ViewPatterns, + DeriveDataTypeable, + DeriveGeneric, + BangPatterns + + -- Other library packages from which modules are imported. + build-depends: + base >=4.8 && <5, + random >=1.0 && <2 + + -- Directories containing source files. + hs-source-dirs: src + + -- Base language which the package is written in. + default-language: Haskell2010 + +-- QuickCheck based test suite +test-suite test-cl3 + type: exitcode-stdio-1.0 + hs-source-dirs: tests + main-is: TestCl3.hs + ghc-options: -Wall -O2 + build-depends: + cl3, + base >=4.8 && <5, + QuickCheck >=2.9 && <3 + default-language: Haskell2010 + +-- Criterion based benchmark +benchmark bench-cl3-nbody + type: exitcode-stdio-1.0 + hs-source-dirs: benchmarks + main-is: NbodyGameCl3.hs + ghc-options: -Wall -O2 + build-depends: + cl3, + base >=4.8 && <5, + criterion >=1.1 && <2 + default-language: Haskell2010 diff --git a/src/Algebra/Geometric/Cl3.hs b/src/Algebra/Geometric/Cl3.hs new file mode 100644 --- /dev/null +++ b/src/Algebra/Geometric/Cl3.hs @@ -0,0 +1,1904 @@ +{-# LANGUAGE Safe #-} +{-# LANGUAGE GADTSyntax #-} +{-# LANGUAGE ViewPatterns #-} +{-# LANGUAGE DeriveDataTypeable #-} +{-# LANGUAGE DeriveGeneric #-} +{-# OPTIONS_GHC -fno-warn-type-defaults #-} +{-# OPTIONS_GHC -fno-warn-unused-top-binds #-} + + +-------------------------------------------------------------------------------------------- +-- | +-- Copyright : (C) 2018 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, + -- * 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, + -- * Helpful Functions + eigvals, hasNilpotent, + spectraldcmp, project +) where + + +import Data.Data (Typeable, Data) +import GHC.Generics (Generic) +import Foreign.Storable (Storable, sizeOf, alignment, peek, poke) +import Foreign.Ptr (Ptr, plusPtr, castPtr) +import System.Random (RandomGen, Random, randomR, random) + + +-- | 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 and live within the APS. +-- +-- * 'R' is the constructor for the Real Scalar Sub-algebra Grade-0 +-- +-- * 'V3' is the Vector constructor Grade-1 +-- +-- * 'BV' is the Bivector constructor Grade-2 +-- +-- * '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 +-- +-- * 'H' is the Quaternion constructor it is the Even Sub-algebra with Grade-0 and Grade-2 elements +-- +-- * 'C' is the Complex constructor it is the Scalar Sub-algebra with Grade-0 and Grade-3 elements +-- +-- * 'BPV' is the Biparavector constructor with Grade-1 and Grade-2 elements +-- +-- * 'ODD' is the Odd constructor with Grade-1 and Grade-3 elements +-- +-- * 'TPV' is the Triparavector constructor with Grade-2 and Grade-3 elements +-- +-- * '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 (G0) + V3 :: !Double -> !Double -> !Double -> Cl3 -- Vectors (G1) + BV :: !Double -> !Double -> !Double -> Cl3 -- Bivectors (G2) + I :: !Double -> Cl3 -- Trivector Imaginary Pseudo-Scalar (G3) + PV :: !Double -> !Double -> !Double -> !Double -> Cl3 -- Paravector (G0 + G1) + H :: !Double -> !Double -> !Double -> !Double -> Cl3 -- Quaternion Even Sub-algebra (G0 + G2) + C :: !Double -> !Double -> Cl3 -- Complex Sub-algebra (G0 + G3) + BPV :: !Double -> !Double -> !Double -> !Double -> !Double -> !Double -> Cl3 -- Biparavector (G1 + G2) + ODD :: !Double -> !Double -> !Double -> !Double -> Cl3 -- Odd (G1 + G3) + TPV :: !Double -> !Double -> !Double -> !Double -> Cl3 -- Triparavector (G2 + G3) + APS :: !Double -> !Double -> !Double -> !Double -> !Double -> !Double -> !Double -> !Double -> Cl3 -- Algebra of Physical Space (G0 + G1 + G2 + G3) + deriving (Show, Read, Typeable, Data, Generic) + + + +-- |'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. 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 + compare cliffor1 cliffor2 = + let (R a0) = abs cliffor1 + (R b0) = abs cliffor2 + in case compare a0 b0 of + EQ -> let (R a0') = lsv cliffor1 + (R b0') = lsv cliffor2 + in compare a0' b0' + LT -> LT + GT -> GT + + +-- |Cl3 has a "Num" instance. "Num" is addition, geometric product, negation, 'abs' the largest +-- singular value, and 'signum' like a unit vector of sorts. +-- +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 smallest singular value is zero then the element is not + -- invertable, we can see here that R, C, V3, BV, and H are all invertable. + 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 (sqrt (a0^2 + a1^2 + a2^2 + a3^2 + 2 * abs a0 * sqrt (a1^2 + a2^2 + a3^2))) + 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) = R (sqrt (a1^2 + a23^2 + a2^2 + a31^2 + a3^2 + a12^2 + + 2 * sqrt ((a1*a31 - a2*a23)^2 + (a1*a12 - a3*a23)^2 + (a2*a12 - a3*a31)^2))) + abs (ODD a1 a2 a3 a123) = R (sqrt (a1^2 + a2^2 + a3^2 + a123^2)) + abs (TPV a23 a31 a12 a123) = R (sqrt (a23^2 + a31^2 + a12^2 + a123^2 + 2 * abs a123 * sqrt (a23^2 + a31^2 + a12^2))) + abs (APS a0 a1 a2 a3 a23 a31 a12 a123) = R (sqrt (a0^2 + a1^2 + a2^2 + a3^2 + a23^2 + a31^2 + a12^2 + a123^2 + + 2 * 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))) + + + -- |'signum' satisfies the Law "abs x * signum x == x" + -- kind of cool: signum of a vector is the unit vector. + signum cliffor + | abs cliffor == 0 = 0 -- initially this was abs cliffor < tol, but this caused problems with 'spectraldcmp' + | otherwise = + let (R mag) = abs cliffor + in cliffor * R (recip mag) + + + -- |'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) + + +-- |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 v@(V3 a1 a2 a3) = + let (R mag) = abs v + sqmag = mag * mag :: Double + in V3 (a1 / sqmag) (a2 / sqmag) (a3 / sqmag) + recip bv@(BV a23 a31 a12) = + let (R mag) = abs bv + sqmag = mag * mag :: Double + in BV (negate $ a23 / sqmag) (negate $ a31 / sqmag) (negate $ a12 / sqmag) + recip i@(I a123) = + let (R mag) = abs i + sqmag = mag * mag :: Double + in I (negate $! a123 / sqmag) + recip pv@PV{} = + let mag = toR $! pv * bar pv + in recip mag * bar pv + recip h@(H a0 a23 a31 a12) = -- H is a division algebra + let (R mag) = abs h + sqmag = mag * mag :: Double + in H (a0 / sqmag) (negate $! a23 / sqmag) (negate $! a31 / sqmag) (negate $! a12 / sqmag) + recip z@(C a0 a123) = -- C is a division algebra + let (R mag) = abs z + sqmag = mag * mag :: Double + in C (a0 / sqmag) (negate $ a123 / sqmag) + recip bpv@BPV{} = reduce $! spectraldcmp recip recip' bpv + recip od@(ODD a1 a2 a3 a123) = + let (R mag) = abs od + sqmag = mag * mag :: Double + in ODD (a1 / sqmag) (a2 / sqmag) (a3 / sqmag) (negate $ a123 / sqmag) + recip tpv@TPV{} = + let mag = toR $! tpv * bar tpv + in recip mag * bar tpv + recip aps@APS{} = reduce $! spectraldcmp recip recip' aps + + -- |'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 = reduce $! spectraldcmp exp exp' cliffor + + -- + log (R a0) | a0 >= 0 = R (log a0) + | otherwise = C (log (negate a0)) pi + log (I a123) = C (log (abs a123)) (signum a123 * (pi/2)) + log (C a0 a123) = C (log (sqrt (a0^2 + a123^2))) (atan2 a123 a0) + log cliffor = reduce $! spectraldcmp log log' cliffor + + -- + sqrt (R a0) | a0 >= 0 = R (sqrt a0) + | otherwise = I (sqrt $ negate a0) + sqrt (I a123) = C u (if a123 < 0 then -v else v) + where v = if u < tol' then 0 else abs a123 / (2 * u) + u = sqrt (abs a123 / 2) + sqrt (C a0 a123) = C u (if a123 < 0 then -v else v) + where (u,v) = if a0 < 0 then (v',u') else (u',v') + v' = if u' < tol' then 0 else abs a123 / (u'*2) + u' = sqrt ((sqrt (a0^2 + a123^2) + abs a0) / 2) + sqrt cliffor = reduce $! spectraldcmp sqrt sqrt' cliffor + + -- + sin (R a0) = R (sin a0) + sin (I a123) = I (sinh a123) + sin (C a0 a123) = C (sin a0 * cosh a123) (cos a0 * sinh a123) + sin cliffor = reduce $! 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 = reduce $! spectraldcmp cos cos' cliffor + + -- + tan (R a0) = R (tan a0) + tan (I a123) = I (tanh a123) + tan (C a0 a123) = C (sinx*coshy) (cosx*sinhy) / C (cosx*coshy) (negate $ sinx*sinhy) + where sinx = sin a0 + cosx = cos a0 + sinhy = sinh a123 + coshy = cosh a123 + tan cliffor = reduce $! spectraldcmp tan tan' cliffor + + -- + asin (R a0) = if (-1) <= a0 && a0 <= 1 then R (asin a0) else asin $ C a0 0 + asin (I a123) = I (asinh a123) + asin (C a0 a123) = C a123' (-a0') + where (C a0' a123') = toC $ log (C (-a123) a0 + sqrt (1 - C a0 a123 * C a0 a123)) -- check this + asin cliffor = reduce $! spectraldcmp asin asin' cliffor + + -- + acos (R a0) = if (-1) <= a0 && a0 <= 1 then R (acos a0) else acos $ C a0 0 + acos (I a123) = C (pi/2) (negate $ asinh a123) + acos (C a0 a123) = C a123'' (-a0'') + where (C a0'' a123'') = log (C a0 a123 + C (-a123') a0') -- check this + (C a0' a123') = sqrt (1 - C a0 a123 * C a0 a123) -- check this + acos cliffor = reduce $! spectraldcmp acos acos' cliffor + + -- + atan (R a0) = R (atan a0) + atan (I a123) = C a123' (-a0') + where (C a0' a123') = toC.log $ ( R (1-a123) / sqrt (R (1 - a123^2))) -- check this + atan (C a0 a123) = C a123' (-a0') + where (C a0' a123') = toC $ log (C (1-a123) a0 / sqrt (1 + C a0 a123 * C a0 a123)) -- check this + atan cliffor = reduce $! 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 = reduce $! 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 = reduce $! spectraldcmp cosh cosh' cliffor + + -- + tanh (R a0) = R (tanh a0) + tanh (I a123) = I (tan a123) + tanh (C a0 a123) = C (cosy*sinhx) (siny*coshx) / C (cosy*coshx) (siny*sinhx) + where siny = sin a123 + cosy = cos a123 + sinhx = sinh a0 + coshx = cosh a0 + tanh cliffor = reduce $! spectraldcmp tanh tanh' cliffor + + -- + asinh (R a0) = R (asinh a0) + asinh (I a123) = log (I a123 + sqrt (R (1 - a123^2))) + asinh (C a0 a123) = log (C a0 a123 + sqrt (1 + C a0 a123 * C a0 a123)) + asinh cliffor = reduce $! spectraldcmp asinh asinh' cliffor + + -- + acosh (R a0) = log (R a0 + sqrt(R a0 - 1) * sqrt(R a0 + 1)) + acosh (I a123) = log (I a123 + sqrt(I a123 - 1) * sqrt(I a123 + 1)) + acosh (C a0 a123) = log (C a0 a123 + sqrt(C a0 a123 - 1) * sqrt(C a0 a123 + 1)) + acosh cliffor = reduce $! spectraldcmp acosh acosh' cliffor + + -- + atanh (R a0) = 0.5 * log (1 + R a0) - 0.5 * log (1 - R a0) + atanh (I a123) = 0.5 * log (1 + I a123) - 0.5 * log (1 - I a123) + atanh (C a0 a123) = 0.5 * log (1 + C a0 a123) - 0.5 * log (1 - C a0 a123) + atanh cliffor = reduce $! 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 (sqrt (a0^2 + a1^2 + a2^2 + a3^2 - + 2 * abs a0 * sqrt (a1^2 + a2^2 + a3^2))) +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) = R (sqrt (a1^2 + a23^2 + a2^2 + a31^2 + a3^2 + a12^2 - + 2 * sqrt ((a1*a31 - a2*a23)^2 + (a1*a12 - a3*a23)^2 + (a2*a12 - a3*a31)^2))) +lsv (ODD a1 a2 a3 a123) = R (sqrt (a1^2 + a2^2 + a3^2 + a123^2)) +lsv (TPV a23 a31 a12 a123) = R (sqrt (a23^2 + a31^2 + a12^2 + a123^2 - (abs a123 + abs a123) * sqrt (a23^2 + a31^2 + a12^2))) +lsv (APS a0 a1 a2 a3 a23 a31 a12 a123) = R (sqrt (a0^2 + a1^2 + a2^2 + a3^2 + a23^2 + a31^2 + a12^2 + a123^2 - + 2 * 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))) + + + +-- | 'spectraldcmp' the spectral decomposition of a function to calculate analytic functions of cliffors in Cl(3,0). +-- This function requires the desired function to be calculated and it'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. +-- 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 _ (reduce -> r@R{}) = fun r +spectraldcmp fun _ (reduce -> v@V3{}) = spectraldcmpSpecial toR fun v -- spectprojR fun v +spectraldcmp fun _ (reduce -> bv@BV{}) = spectraldcmpSpecial toI fun bv -- spectprojI fun bv +spectraldcmp fun _ (reduce -> i@I{}) = fun i +spectraldcmp fun _ (reduce -> pv@PV{}) = spectraldcmpSpecial toR fun pv -- spectprojR fun pv +spectraldcmp fun _ (reduce -> h@H{}) = spectraldcmpSpecial toC fun h -- spectprojC fun h +spectraldcmp fun _ (reduce -> c@C{}) = fun c +spectraldcmp fun fun' (reduce -> bpv@BPV{}) + | hasNilpotent bpv = jordan fun fun' bpv -- jordan normal form Cl3 style + | isColinear bpv = spectraldcmpSpecial toC fun bpv -- spectprojC fun bpv + | otherwise = -- transform it so it will be colinear + let (v,d,v_bar) = boost2colinear bpv + in v * spectraldcmpSpecial toC fun d * v_bar -- v * spectprojC fun d * v_bar +spectraldcmp fun _ (reduce -> od@ODD{}) = spectraldcmpSpecial toC fun od -- spectprojC fun od +spectraldcmp fun _ (reduce -> tpv@TPV{}) = spectraldcmpSpecial toI fun tpv -- spectprojI fun tpv +spectraldcmp fun fun' (reduce -> aps@APS{}) + | hasNilpotent aps = jordan fun fun' aps -- jordan normal form Cl3 style + | isColinear aps = spectraldcmpSpecial toC fun aps -- spectprojC fun aps + | otherwise = -- transform it so it will be colinear + let (v,d,v_bar) = boost2colinear aps + in v * spectraldcmpSpecial toC fun d * v_bar -- v * spectprojC fun d * v_bar +spectraldcmp _ _ _ = error "Major problems with 'spectraldcmp' or 'reduce'" +-- + +-- | '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 'spectproj'. 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 +jordan fun fun' cliffor = + let eigs = toC 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 -> r@R{}) = (r,r) +eigvals (reduce -> v@V3{}) = eigvalsSpecial toR v -- eigvalsR v +eigvals (reduce -> bv@BV{}) = eigvalsSpecial toI bv -- eigvalsI bv +eigvals (reduce -> i@I{}) = (i,i) +eigvals (reduce -> pv@PV{}) = eigvalsSpecial toR pv -- eigvalsR pv +eigvals (reduce -> h@H{}) = eigvalsSpecial toC h -- eigvalsC h +eigvals (reduce -> c@C{}) = (c,c) +eigvals (reduce -> bpv@BPV{}) + | hasNilpotent bpv = (0,0) -- this case is actually nilpotent + | isColinear bpv = eigvalsSpecial toC bpv -- eigvalsC bpv + | otherwise = -- transform it so it will be colinear + let (_,d,_) = boost2colinear bpv + in eigvalsSpecial toC d -- eigvalsC d +eigvals (reduce -> od@ODD{}) = eigvalsSpecial toC od -- eigvalsC od +eigvals (reduce -> tpv@TPV{}) = eigvalsSpecial toI tpv -- eigvalsI tpv +eigvals (reduce -> aps@APS{}) + | hasNilpotent aps = (toC aps,toC aps) -- a scalar plus nilpotent + | isColinear aps = eigvalsSpecial toC aps -- eigvalsC aps + | otherwise = -- transform it so it will be colinear + let (_,d,_) = boost2colinear aps + in eigvalsSpecial toC d -- eigvalsC d +eigvals _ = error "Major issues with 'eigvals' or 'reduce'" +-- + +-- | '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. +-- We have safty problem with unreduced values, so it calls reduce first, as a view pattern. +project :: Cl3 -> Cl3 +project (reduce -> R{}) = PV 0.5 0 0 0.5 -- default to e3 direction +project (reduce -> v@V3{}) = 0.5 * (1 + signum v) +project (reduce -> bv@BV{}) = 0.5 * (1 + signum (toV3 $ mI * toBV bv)) +project (reduce -> I{}) = PV 0.5 0 0 0.5 -- default to e3 direction +project (reduce -> pv@PV{}) = 0.5 * (1 + signum (toV3 pv)) +project (reduce -> h@H{}) = 0.5 * (1 + signum (toV3 $ mI * toBV h)) +project (reduce -> C{}) = PV 0.5 0 0 0.5 -- default to e3 direction +project (reduce -> bpv@BPV{}) + | abs (toV3 bpv + toV3 (mI * toBV bpv)) <= tol = 0.5 * (1 + signum (toV3 bpv)) -- gaurd for equal and opposite + | otherwise = 0.5 * (1 + signum (toV3 bpv + toV3 (mI * toBV bpv))) +project (reduce -> od@ODD{}) = 0.5 * (1 + signum (toV3 od)) +project (reduce -> tpv@TPV{}) = 0.5 * (1 + signum (toV3 $ mI * toBV tpv)) +project (reduce -> aps@APS{}) = project.toBPV $ aps +project (reduce -> _) = error "Error: Got some serious issues with 'project' and/or 'reduce'. Please Fix." + +-- | 'boost2colinear' calculates a boost that is perpendicular to both the vector and bivector +-- components, 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 = v * d * bar v and returns v, d, and v_bar as a tuple. First v must be calculated +-- and then d = bar v * cliffor * v. d will have colinear vector and bivector parts. +-- This is somewhat simular to finding the drift frame for an electromagnetic field. +boost2colinear :: Cl3 -> (Cl3, Cl3, Cl3) +boost2colinear cliffor = + let v = toV3 cliffor -- extract the vector + bv = toV3 $ mI * toBV cliffor -- extract the bivector and turn it into a vector + -- Find an orthonormal basis natural to the cliffor (eigen basis) + sum_direction = signum $ v + bv -- the natural basis is the sum of the vector and bivector + orthogonal_direction = signum.toV3 $ mI * toBV (v * bv) -- the natural basis is perpedicualr to both the vector and bivector + other_direction = signum.toV3 $ mI * toBV (sum_direction * orthogonal_direction) -- the natural basis is orthoganl to both the sum and ortho basis + -- Decompose the cliffor in our new basis via dot product + -- this decpomosition is in the plane of the biparavector + (C a1 a23) = toC $ other_direction * cliffor + (C a3 a12) = toC $ sum_direction * cliffor + -- Find the boost to make the vector and bivector parts colinear (a two page derivation) + sum_sq = a1^2 + a3^2 + a23^2 + a12^2 + numerator = 2 * (a1 * a12 - a3 * a23) + tanh4eta = numerator / sum_sq + _4eta = atanh tanh4eta + eta = _4eta / 4 + boost = exp (R eta * orthogonal_direction) + -- calculate the returned values + boost_bar = bar boost + d = boost_bar * cliffor * boost + in (boost, d, boost_bar) + + +-- | 'isColinear' takes a Cliffor and determines if the vector part and the bivector part are +-- not at all orthoganl and non-zero. +isColinear :: Cl3 -> Bool +isColinear cliffor = abs (toV3 cliffor) /= 0 && abs (mI * toBV cliffor) /= 0 && -- Non-Zero + abs (toBV $ signum (toV3 cliffor) * signum (mI * toBV cliffor)) <= tol -- Not Orthoganl + + +-- | '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 cliffor = abs (toV3 cliffor) /= 0 && abs (mI * toBV cliffor) /= 0 && -- Non-Zero + abs (toR $ signum (toV3 cliffor) * signum (mI * toBV cliffor)) <= tol && -- Orthoganl + abs (abs (toV3 cliffor) - abs (toBV cliffor)) <= tol -- Equal Magnitude + + +-- | '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 some are zero +reduce :: Cl3 -> Cl3 +reduce r@R{} = r +reduce v@V3{} + | abs v <= tol = R 0 + | otherwise = v +reduce bv@BV{} + | abs bv <= tol = R 0 + | otherwise = bv +reduce i@I{} + | abs i <= tol = R 0 + | otherwise = i +reduce pv@PV{} + | abs pv <= tol = R 0 + | abs (toR pv) <= tol = toV3 pv + | abs (toV3 pv) <= tol = toR pv + | otherwise = pv +reduce h@H{} + | abs h <= tol = R 0 + | abs (toR h) <= tol = toBV h + | abs (toBV h) <= tol = toR h + | otherwise = h +reduce c@C{} + | abs c <= tol = R 0 + | abs (toR c) <= tol = toI c + | abs (toI c) <= tol = toR c + | otherwise = c +reduce bpv@BPV{} + | abs bpv <= tol = R 0 + | abs (toV3 bpv) <= tol = toBV bpv + | abs (toBV bpv) <= tol = toV3 bpv + | otherwise = bpv +reduce od@ODD{} + | abs od <= tol = R 0 + | abs (toV3 od) <= tol = toI od + | abs (toI od) <= tol = toV3 od + | otherwise = od +reduce tpv@TPV{} + | abs tpv <= tol = R 0 + | abs (toBV tpv) <= tol = toI tpv + | abs (toI tpv) <= tol = toBV tpv + | otherwise = tpv +reduce aps@APS{} + | abs aps <= tol = R 0 + | abs (toC aps) <= tol = reduce (toBPV aps) + | abs (toBPV aps) <= tol = reduce (toC aps) + | abs (toH aps) <= tol = reduce (toODD aps) + | abs (toODD aps) <= tol = reduce (toH aps) + | abs (toPV aps) <= tol = reduce (toTPV aps) + | abs (toTPV aps) <= tol = reduce (toPV aps) + | otherwise = aps + +-- | 'mI' negative i +mI :: Cl3 +mI = I (-1) + +-- | 'tol' currently 128*eps +tol :: Cl3 +tol = R $ 128 * 1.1102230246251565e-16 + +tol' :: Double +tol' = 128 * 1.1102230246251565e-16 + + +-- | '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 Cliffor 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' x = negate.recip $ x * x -- pole at 0 + +exp' :: Cl3 -> Cl3 +exp' = exp + +log' :: Cl3 -> Cl3 +log' = recip -- pole at 0 + +sqrt' :: Cl3 -> Cl3 +sqrt' x = 0.5 * recip (sqrt x) -- pole at 0 + +sin' :: Cl3 -> Cl3 +sin' = cos + +cos' :: Cl3 -> Cl3 +cos' = negate.sin + +tan' :: Cl3 -> Cl3 +tan' x = recip (cos x) * recip (cos x) -- pole at pi/2*n for all integers + +asin' :: Cl3 -> Cl3 +asin' x = recip.sqrt $ 1 - (x * x) -- pole at +/-1 + +acos' :: Cl3 -> Cl3 +acos' x = negate.recip.sqrt $ 1 - (x * x) -- pole at +/-1 + +atan' :: Cl3 -> Cl3 +atan' x = recip $ 1 + (x * x) -- pole at +/-i + +sinh' :: Cl3 -> Cl3 +sinh' = cosh + +cosh' :: Cl3 -> Cl3 +cosh' = sinh + +tanh' :: Cl3 -> Cl3 +tanh' x = recip (cosh x) * recip (cosh x) + +asinh' :: Cl3 -> Cl3 +asinh' x = recip.sqrt $ (x * x) + 1 -- pole at +/-i + +acosh' :: Cl3 -> Cl3 +acosh' x = recip $ sqrt (x - 1) * sqrt (x + 1) -- pole at +/-1 + +atanh' :: Cl3 -> Cl3 +atanh' x = recip $ 1 - (x * x) -- 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" + + + + +------------------------------------------------------------------- +-- +-- 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' + + 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 + 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, g') = rangeR (lo, hi) g + (v3, g'') = rangeV3 (lo, hi) g' + in (r + v3, 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, g') = rangeR (lo, hi) g + (bv, g'') = rangeBV (lo, hi) g' + in (r + bv, 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, g') = rangeR (lo, hi) g + (i, g'') = rangeI (lo, hi) g' + in (r + i, 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 (v3, g') = rangeV3 (lo, hi) g + (bv, g'') = rangeBV (lo, hi) g' + in (v3 + bv, 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 (v3, g') = rangeV3 (lo, hi) g + (i, g'') = rangeI (lo, hi) g' + in (v3 + i, 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 (bv, g') = rangeBV (lo, hi) g + (i, g'') = rangeI (lo, hi) g' + in (bv + i, 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 (pv, g') = rangePV (lo, hi) g + (tpv, g'') = rangeTPV (lo, hi) g' + in (pv + tpv, g'') + + +------------------------------------------------------------------- +-- Additional Random generators +------------------------------------------------------------------- +-- | 'randUnitV3' a unit vector with a random direction +randUnitV3 :: RandomGen g => g -> (Cl3, g) +randUnitV3 g = + let (theta, g') = randomR (0,pi) g + (phi, g'') = randomR (0,2*pi) g' + in (V3 (sin theta * cos phi) (sin theta * sin phi) (cos theta), g'') + + +-- | 'randProjector' a projector with a random direction +randProjector :: RandomGen g => g -> (Cl3, g) +randProjector g = + let (v3, g') = randUnitV3 g + in (0.5 + 0.5 * v3, g') + + +-- | 'randNilpotent' a nilpotent element with a random orientation +randNilpotent :: RandomGen g => g -> (Cl3, g) +randNilpotent g = + let (p, g') = randProjector g + (v, g'') = randUnitV3 g' + vnormal = signum $ I (-1) * toBV ( toV3 p * v) -- unit vector normal to the projector + in (toBPV $ vnormal * p, 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 + (theta, g'') = randomR (0,pi) g' + (phi, g''') = randomR (0,2*pi) g'' + in (con (mag * sin theta * cos phi) (mag * sin theta * sin phi) (mag * cos theta), g''') + diff --git a/src/Algebra/Geometric/Cl3/JonesCalculus.hs b/src/Algebra/Geometric/Cl3/JonesCalculus.hs new file mode 100644 --- /dev/null +++ b/src/Algebra/Geometric/Cl3/JonesCalculus.hs @@ -0,0 +1,212 @@ +{-# LANGUAGE Safe #-} +{-# LANGUAGE ViewPatterns #-} +{-# OPTIONS_GHC -fno-warn-missing-signatures #-} + +-------------------------------------------------------------------------------------------- +-- | +-- Copyright : (C) 2018 Nathan Waivio +-- License : BSD3 +-- Maintainer : Nathan Waivio +-- Stability : Stable +-- Portability : unportable +-- +-- Library implementing standard functions for the Jones Calculus in the Cl3 Library. +-- This implementation of the Jones Calculus is based on the convensions of SPIE's Field Guide to Polarization (ϕ = ω t − k z). +-- +-- * E. Collett, Field Guide to Polarization, SPIE Field Guides vol. FG05, SPIE (2005). ISBN 0-8194-5868-6. +-- +-- +-- = Jones Vectors +-- +-- Within the system of the Bloch Sphere, the Jones Vectors in Cl3 are calculated +-- by generating the left ideal of the rotation of a unit vector to the 'e3' basis. +-- Standard form for for a versor is 'rot = exp $ (-i/2) * theta * u' for angle 'theta' +-- and the rotational axis unit vector 'u'. +-- +-- Bloch Sphere Coordinates: +-- +-- @ +-- e3 +-- | +-- |____e2 +-- / +-- / +-- e1 +-- @ +-- +-------------------------------------------- + +module Algebra.Geometric.Cl3.JonesCalculus +( + -- * Jones Vectors + hpv, vpv, + dpv, apv, + rpv, lpv, + jv, + -- * Jones Matrices + hpm, vpm, + dpm, apm, + rpm, lpm, + jm, + hpmRot, + -- * Wave Plates + qwp, hwp, + qwpRot, hwpRot, + wp, + wpRot, + -- * Reflection + refl, + -- * Random Jones Vectors + randJonesVec, + randOrthogonalJonesVec, + -- * Normalization Factorization + factorize +) where + + +import safe Algebra.Geometric.Cl3 (Cl3(..), dag, bar, toR, toV3, toC, project, randUnitV3) +import System.Random (RandomGen) + + +e0 = R 1 +e1 = V3 1 0 0 +e2 = V3 0 1 0 +e3 = V3 0 0 1 + +i = I 1 + +p1 = 0.5 * (e0 + e1) +p2 = 0.5 * (e0 + e2) +p3 = 0.5 * (e0 + e3) + + + + +-- | 'hpv' horizontally polarized Jones vector +hpv = signum $ e0 * p3 -- e0 == exp $ (-i/2) * 0 * e2, any vector orthoganl to e3 could have been selected as the rotational axis because the angle is zero + +-- | 'vpv' vertically polarized Jones vector +vpv = signum $ exp ((-i/2) * pi * e2) * p3 -- e2 is selected to obtain the standard form, e1 or any vector orthoganl to e3 could have been selected + +-- | 'dpv' diagonally polarized Jones vector +dpv = signum $ exp ((-i/2) * (pi/2) * e2) * p3 -- rotate -e1 to e3 around rotational axis e2, an angle of pi/2 + +-- | 'apv' anti-diagonally polarized Jones vector +apv = signum $ exp ((-i/2) * (pi/2) * (-e2)) * p3 -- rotate e1 to e3 around rotational axis -e2, an angle of pi/2 + +-- | 'rpv' right hand circularly polarized Jones vector +rpv = signum $ exp ((-i/2) * (pi/2) * (-e1)) * p3 -- rotate -e2 to e3 around rotational axis -e1, and angle of pi/2 + +-- | 'lpv' left hand circularly polarized Jones vector +lpv = signum $ exp ((-i/2) * (pi/2) * e1) * p3 -- rotate e2 to e3 around rotational axis e1, an angle of pi/2 + +-- | 'jv' function that returns Jones vector from input vector unit vector +-- currently converts the input to a unit vector +jv (signum.toV3 -> v) | v == e3 = hpv + | v == -e3 = vpv + | otherwise = signum $ sqrt (e3 * v) * p3 + + +-- | 'hpm' Horizontal Polarizer Jones Matrix +hpm = p3 + +-- | 'vpm' Vertical Polarizer Jones Matrix +vpm = bar p3 + +-- | 'dpm' Diagonal Polarizer Jones Matrix +dpm = p1 + +-- | 'apm' Anti-diagonal Polarizer Jones Matrix +apm = bar p1 + +-- | 'rpm' Right Hand Circular Polarizer Jones Matrix +rpm = p2 + +-- | 'lpm' Left Hand Circular Polarizer Jones Matrix +lpm = bar p2 + + +-- | 'jm' funciton that returns a Jones Matrix from an input Bloch Vector +-- currently converts the input to a unit vector +jm (signum.toV3 -> v) = project v + +-- | 'rot' will produce a versor that rotates a vector by an angle about +-- a unit vector axis. +rot (toR -> theta) (signum.toV3 -> axis) = exp $ (-i/2) * theta * axis + +-- | 'rotIsh' will produce a versor that rotates by double the input angle +-- for rotating polarizers and wave plates the axis is e2. +rotIsh (toR -> theta) = rot (2*theta) e2 + +-- | 'hpmRot' Jones matrix for a rotated ideal Linear Horizontal Polarizer. +-- Input value should be a scalar angle in Radians. +hpmRot (toR -> theta) = + let roted = rotIsh theta + in roted * hpm * dag roted + +-- | 'qwp' Quarter Wave Plate Jones Matrix +qwp = p3 - i * bar p3 + +-- | 'qwpRot' Rotated Quarter Wave Plate Jones Matrix. +-- Input value should be a scalar angle in Radians. +qwpRot (toR -> theta) = + let roted = rotIsh theta + in roted * qwp * dag roted + +-- | 'hwp' Half Wave Plate Jones Matrix +hwp = e3 + +-- | 'hwpRot' Rotated Half Wave Plate Jones Matrix. +-- Input value should be a scalar angle in Radians. +hwpRot (toR -> theta) = + let roted = rotIsh theta + in roted * hwp * dag roted + +-- | 'wp' a Wave Plate with phase shift of phi Jones Matrix. +-- Input value should be a scalar angle in Radians. +wp (toR -> phi) = exp $ (i/2) * phi * e3 + +-- | 'wpRot' a Rotated Wave Plate with phase shift of phi and rotation theta Jones Matrix. +-- The first input value is phi the phase shift as a scalar value in Radians. The second +-- input value is theta the rotation a scalar angle in Radians. +wpRot (toR -> phi) (toR -> theta) = + let roted = rotIsh theta + in roted * wp phi * dag roted + +-- | 'refl' a Refelection Jones Matrix +refl = e3 + + +-- | 'factorize' is a function that takes an Jones Vector after transformation by an +-- optical chain, and returns the amplitude (amp), phase (phi), and normalized Jones +-- Vector (vec), by the factorization of the input such that: @__amp * exp (i*phi/2) * vec__@ +factorize :: Cl3 -> (Cl3,Cl3,Cl3) +factorize jonesVec = + let c = toC jonesVec + jonesVec' = recip c * jonesVec + ampC = abs c + ampJonesVec' = abs jonesVec' + normJonesVec = recip ampJonesVec' * jonesVec' + amp = ampC * ampJonesVec' + normC = recip ampC * c + phi = 2 * (-i) * log normC + in (amp, phi, normJonesVec) + +------------------------------------------------------------------- +-- +-- Random Jones Vectors +-- +------------------------------------------------------------------- + +-- | 'randJonesVec' a Random Jones Vector. +randJonesVec :: RandomGen g => g -> (Cl3, g) +randJonesVec g = + let (v3, g') = randUnitV3 g + in (jv v3,g') + +-- | 'randOrthogonalJonesVec' a Random Orthogonal Complementary pair of Jones +-- Vectors. +randOrthogonalJonesVec :: RandomGen g => g -> ((Cl3, Cl3), g) +randOrthogonalJonesVec g = + let (v3, g') = randUnitV3 g + in ((jv v3, jv (bar v3)),g') diff --git a/tests/TestCl3.hs b/tests/TestCl3.hs new file mode 100644 --- /dev/null +++ b/tests/TestCl3.hs @@ -0,0 +1,432 @@ +{-# LANGUAGE ViewPatterns #-} +{-# OPTIONS_GHC -fno-warn-type-defaults #-} + + +------------------------------------------------------------------- +-- | +-- Copyright : (c) 2017 Nathan Waivio +-- License : BSD3 +-- Maintainer : Nathan Waivio +-- +-- A program to test Algebra.Geometric.Cl3 +-- The code runs tests on some standard test input and then +-- runs quckcheck for some trig identities. +-- +------------------------------------------------------------------- + +module Main (main) where + +import Test.QuickCheck (Arbitrary, arbitrary, oneof, suchThat, quickCheckWith, stdArgs, maxSuccess) +import Algebra.Geometric.Cl3 +import Control.Applicative ((<*>)) + + +------------------------------------------------------------------ +-- | +-- This program verifies the approximate equality of various trig +-- identities to the with the following limitations: +-- +-- * The magnitude of the cliffor is limited in some cases. +-- +-- * The imaginary part of the eigenvalues are unwrapped, due to the cyclical nature of some of the results, in a few cases. +-- +-- * The poles of the functions are excluded. +-- +-- * The poles of the derivatives of the functions are excluded when the cliffor is has a nilpotent component. +-- +-- * Approximate equivalence is tested due to limitations with respect to floating point math. +-- +-- * The implementation of Arbitrary for Cl3 limits the arbitrary cliffor such that the absolute value of cliff is less than 15 +-- +-- The following properties are verified in this module: +-- +-- * log.exp Identity +-- +-- * exp.log Identity +-- +-- * abs*signum law +-- +-- * The definition of recip +-- +-- * recip.recip Identity +-- +-- * sin.asin Identity +-- +-- * asin.sin Identity +-- +-- * cos.acos Identity +-- +-- * acos.cos Identity +-- +-- * sinh.asinh Identity +-- +-- * asinh.sinh Identity +-- +-- * cosh.acosh Identity +-- +-- * acosh.cosh Identity +-- +-- * Double Sin Identity +-- +-- * Double Cos Identity +-- +-- * Double Tan Identity +-- +-- * Double Sinh Identity +-- +-- * Double Cosh Identity +-- +-- * Double Tanh Identity +-- +-- * Positive Sin Shift Identity +-- +-- * Negative Sin Shift Identity +-- +-- * sin^2+cos^2 Identity +-- +-- * cosh^2-sinh^2 Identity +-- +-- * Symmetry of Cosh +-- +-- * Symmetry of Sinh +-- +-- * Double I Sin +-- +-- * Composition Algebra Tests +-- +------------------------------------------------------------------- + +main :: IO () +main = do moduleTests + print "Testing log.exp Identity:" + quickCheckWith stdArgs { maxSuccess = 30000 } prop_LogExp + print "Testing exp.log Identity:" + quickCheckWith stdArgs { maxSuccess = 30000 } prop_ExpLog + print "Testing abs*signum law:" + quickCheckWith stdArgs { maxSuccess = 30000 } prop_AbsSignum + print "Testing the definition of recip:" + quickCheckWith stdArgs { maxSuccess = 30000 } prop_RecipDef + print "Testing recip.recip Identity:" + quickCheckWith stdArgs { maxSuccess = 30000 } prop_RecipID + print "Testing sin.asin Identity:" + quickCheckWith stdArgs { maxSuccess = 30000 } prop_SinAsin + print "Testing asin.sin Identity:" + quickCheckWith stdArgs { maxSuccess = 30000 } prop_AsinSin + print "Testing cos.acos Identity:" + quickCheckWith stdArgs { maxSuccess = 30000 } prop_CosAcos + print "Testing acos.cos Identity:" + quickCheckWith stdArgs { maxSuccess = 30000 } prop_AcosCos + print "Testing sinh.asinh Identity:" + quickCheckWith stdArgs { maxSuccess = 30000 } prop_SinhAsinh + print "Testing asinh.sinh Identity:" + quickCheckWith stdArgs { maxSuccess = 30000 } prop_AsinhSinh + print "Testing cosh.acosh Identity:" + quickCheckWith stdArgs { maxSuccess = 30000 } prop_CoshAcosh + print "Testing acosh.cosh Identity:" + quickCheckWith stdArgs { maxSuccess = 30000 } prop_AcoshCosh + print "Testing acosh.cosh Identity2:" + quickCheckWith stdArgs { maxSuccess = 30000 } prop_AcoshCosh2 + print "Testing Double Sin Identity:" + quickCheckWith stdArgs { maxSuccess = 30000 } prop_DubSin + print "Testing Double Cos Identity:" + quickCheckWith stdArgs { maxSuccess = 30000 } prop_DubCos + print "Testing Double Tan Identity:" + quickCheckWith stdArgs { maxSuccess = 30000 } prop_DubTan + print "Testing Double Sinh Identity:" + quickCheckWith stdArgs { maxSuccess = 30000 } prop_DubSinh + print "Testing Double Cosh Identity:" + quickCheckWith stdArgs { maxSuccess = 30000 } prop_DubCosh + print "Testing Double Tanh Identity:" + quickCheckWith stdArgs { maxSuccess = 30000 } prop_DubTanh + print "Testing Positive Sin Shift Identity:" + quickCheckWith stdArgs { maxSuccess = 30000 } prop_PosSinShift + print "Testing Negative Sin Shift Identity:" + quickCheckWith stdArgs { maxSuccess = 30000 } prop_NegSinShift + print "Testing sin^2+cos^2 Identity:" + quickCheckWith stdArgs { maxSuccess = 30000 } prop_SinSqCosSq + print "Testing cosh^2-sinh^2 Identity:" + quickCheckWith stdArgs { maxSuccess = 30000 } prop_CoshSqmSinhSq + print "Testing Symmetry of Cosh:" + quickCheckWith stdArgs { maxSuccess = 30000 } prop_SymCosh + print "Testing Symmetry of Sinh:" + quickCheckWith stdArgs { maxSuccess = 30000 } prop_SymSinh + print "Testing Double I Sin:" + quickCheckWith stdArgs { maxSuccess = 30000 } prop_DoubleISin + print "Is has Composition Sub-Algebras:" + quickCheckWith stdArgs { maxSuccess = 30000 } prop_CompAlg + + + +---------------------------------------------------------- +-- |Start of Module Tests +moduleTests :: IO () +moduleTests = sequence_ $ tests <*> inputs + +inputs :: [Cl3] +inputs = [R 0 + ,APS 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 + ,PV 0.5 0 0 0.5 + ,PV 0.5 0 0 (-0.5) + ,BPV 0.5 0 0 0 (-0.5) 0 + ,BPV 0.5 0 0 0 0.5 0 + ,R 1 + ,V3 1 0 0 + ,APS 1 0.5 0 0 0 0.5 0 0 + ,APS 1 0.5 0 0 0 (-0.5) 0 0 + ,PV 1 1 0 0 + ,R 1 + ,R (-1) + ,V3 1 0 0 + ,V3 (-1) 0 0 + ,V3 0 1 0 + ,V3 0 (-1) 0 + ,V3 0 0 1 + ,V3 0 0 (-1) + ,BV 1 0 0 + ,BV (-1) 0 0 + ,BV 0 1 0 + ,BV 0 (-1) 0 + ,BV 0 0 1 + ,BV 0 0 (-1) + ,I 1 + ,I (-1) + ] + +-- | 'tests' is a list of tests +-- The out of bounds can be the poles of the function or if the cliffor has +-- nilpotent content then the poles of the derivative as well as the function +tests :: [Cl3 -> IO()] +tests = [runTest "Log.Exp Identity" (log.exp) id (const False) + ,runTest "Exp.Log Identity" (exp.log) id (\z -> lsv z < tol) -- singular inputs are out of bounds + ,runTest "Abs*Signum Identity" (\x->abs x * signum x) id (const False) + ,runTest "Reciprical Identity" (recip.recip) id (\z -> lsv z < tol) -- singular inputs are out of bounds + ,runTest "sin.asin" (sin.asin) id (\z -> hasNilpotent z && poles [R 1, R (-1)] z) + ,runTest "asin.sin" (asin.sin) (\z -> negate (I 1) * log (0.5 * (exp (I 1 * z) - exp (I (-1) * z)) + + sqrt (1+0.25*(exp (I (-1) * z) - exp (I 1 * z))^2))) (const False) + ,runTest "cos.acos" (cos.acos) id (\z -> hasNilpotent z && poles [R 1, R (-1)] z) + ,runTest "acos.cos" (acos.cos) (\z -> 0.5 * (pi - 2 * asin(cos z))) (\z -> hasNilpotent z && poles [R 0, pi, negate pi] z) + ,runTest "sinh.asinh" (sinh.asinh) id (const False) + ,runTest "asinh.sinh" (asinh.sinh) (\z -> log (0.5*(exp z - exp (negate z)) + sqrt (0.25 * (exp z - exp (negate z))^2 + 1))) (const False) + ,runTest "cosh.acosh" (cosh.acosh) id (\z -> hasNilpotent z && poles [R 1, R (-1)] z) + ,runTest "acosh.cosh" (acosh.cosh) (\z -> log (0.5*(exp z + exp (negate z)) + + sqrt (0.5*(exp z + exp (negate z)) - 1) * sqrt (0.5*(exp z + exp (negate z)) + 1))) (const False) + ,runTest "Double Angle sin" (\z -> sin (2 * z)) (\z -> 2 * sin z * cos z) (const False) + ,runTest "Double Angle cos" (\z -> cos (2 * z)) (\z -> cos z ^ 2 - sin z ^ 2) (const False) + ,runTest "Double Angle tan" (\z -> tan (2 * z)) (\z -> (2 * tan z) / (1 - tan z ^ 2)) (const False) + ,runTest "+Sin Shift" (\z -> sin (pi/2 + z)) cos (const False) + ,runTest "-Sin Shift" (\z -> sin (pi/2 - z)) cos (const False) + ,runTest "Double Angle sinh" (\z -> sinh (2 * z)) (\z -> 2 * sinh z * cosh z) (const False) + ,runTest "Double Angle cosh" (\z -> cosh (2 * z)) (\z -> 2 * cosh z ^ 2 - 1) (const False) + ,runTest "Double Angle tanh" (\z -> tanh (2 * z)) (\z -> (2 * tanh z) / (1 + tanh z ^ 2)) (const False) + ,runTest "sin^2+cos^2" (\z -> sin z ^ 2 + cos z ^ 2) (const $ R 1) (const False) + ,runTest "cosh^2-sinh^2" (\z -> cosh z ^ 2 - sinh z ^ 2) (const $ R 1) (const False) + ,runTest "Symetry of cosh" (cosh.negate) cosh (const False) + ,runTest "Symetry of sinh" (sinh.negate) (negate.sinh) (const False) + ,runTest "sin.acos" (sin.acos) (\z -> sqrt (1 - z^2)) (\z -> hasNilpotent z && poles [R 1, R (-1)] z) + ,runTest "sin.atan" (sin.atan) (\z -> z / sqrt (1 + z^2)) (poles [I 1, I (-1)]) + ,runTest "cos.atan" (cos.atan) (\z -> recip.sqrt $ 1 + z^2) (poles [I 1, I (-1)]) + ,runTest "cos.asin" (cos.asin) (\z -> sqrt (1 - z^2)) (\z -> hasNilpotent z && poles [R 1, R (-1)] z) + ,runTest "tan.asin" (tan.asin) (\z -> z / sqrt (1 - z^2)) (poles [R 1, R (-1)]) + ,runTest "tan.acos" (tan.acos) (\z -> sqrt (1 - z^2) / z) (\z -> if hasNilpotent z then poles [R 1, R 0, R (-1)] z else poles [R 0] z) + ] + +-- | The Properties +prop_LogExp :: ArbCl3 -> Bool +prop_LogExp (Arb cliffor) = (abs cliffor > 10) || ( + let cliffor' = unWrapIPartEigs cliffor -- imaginary part of log.exp repeats +-- round off errors get large for exp larger than 5 use spectproj (log.exp) for accuracy + in log (exp cliffor') ≈≈ cliffor') + +-- log 0 is -Inf, Infinite vectors don't play nice +-- spectproj (exp.log) doesn't have this issue +prop_ExpLog :: ArbCl3 -> Bool +prop_ExpLog (Arb cliffor) = (lsv cliffor < tol) || (exp (log cliffor) ≈≈ cliffor) + +prop_AbsSignum :: ArbCl3 -> Bool +prop_AbsSignum (Arb cliffor) = abs cliffor * signum cliffor ≈≈ cliffor + +prop_RecipDef :: ArbCl3 -> Bool +prop_RecipDef (Arb cliffor) = (lsv cliffor < tol) || (recip cliffor * cliffor ≈≈ 1) + +-- singular inputs don't recip also suffers from roundoff errors at large values +prop_RecipID :: ArbCl3 -> Bool +prop_RecipID (Arb cliffor) = (lsv cliffor < tol) || (recip (recip cliffor) ≈≈ cliffor) + +prop_SinAsin :: ArbCl3 -> Bool +prop_SinAsin (Arb cliffor) = if hasNilpotent cliffor + then poles [R 1, R (-1)] cliffor || (sin (asin cliffor) ≈≈ cliffor) + else sin (asin cliffor) ≈≈ cliffor + +prop_AsinSin :: ArbCl3 -> Bool +prop_AsinSin (Arb cliffor) = (abs cliffor > 10) || (asin (sin cliffor) ≈≈ (I (-1) * log (0.5 * (exp (I 1 * cliffor) - exp (I (-1) * cliffor)) + + sqrt (1+0.25*(exp (I (-1) * cliffor) - exp (I 1 * cliffor))^2)))) + +prop_CosAcos :: ArbCl3 -> Bool +prop_CosAcos (Arb cliffor) = if hasNilpotent cliffor + then poles [R 1, R (-1)] cliffor || (cos (acos cliffor) ≈≈ cliffor) + else cos (acos cliffor) ≈≈ cliffor + +prop_AcosCos :: ArbCl3 -> Bool +prop_AcosCos (Arb cliffor) = (abs cliffor > 10) || (if hasNilpotent cliffor + then poles [R 0, pi, negate pi] cliffor || (acos (cos cliffor) ≈≈ 0.5 * (pi - 2 * asin(cos cliffor))) + else acos (cos cliffor) ≈≈ 0.5 * (pi - 2 * asin(cos cliffor))) + +prop_SinhAsinh :: ArbCl3 -> Bool +prop_SinhAsinh (Arb cliffor) = sinh (asinh cliffor) ≈≈ cliffor + +prop_AsinhSinh :: ArbCl3 -> Bool +prop_AsinhSinh (Arb cliffor) = (abs cliffor > 10) || (asinh (sinh cliffor) ≈≈ log (0.5*(exp cliffor - exp (negate cliffor)) + + sqrt (0.25 * (exp cliffor - exp (negate cliffor))^2 + 1))) + +prop_CoshAcosh :: ArbCl3 -> Bool +prop_CoshAcosh (Arb cliffor) = if hasNilpotent cliffor + then poles [R 1, R (-1)] cliffor || (cosh (acosh cliffor) ≈≈ cliffor) + else cosh (acosh cliffor) ≈≈ cliffor + +prop_AcoshCosh :: ArbCl3 -> Bool +prop_AcoshCosh (Arb cliffor) = acosh (cosh cliffor) ≈≈ log (0.5*(exp cliffor + exp (negate cliffor)) + + sqrt (0.5*(exp cliffor + exp (negate cliffor)) - 1) * + sqrt (0.5*(exp cliffor + exp (negate cliffor)) + 1)) + +prop_AcoshCosh2 :: ArbCl3 -> Bool +prop_AcoshCosh2 (Arb cliffor) = acosh (cosh cliffor) ≈≈ log (cosh cliffor + sqrt (cosh cliffor - 1) * sqrt (cosh cliffor + 1)) + +prop_DubSin :: ArbCl3 -> Bool +prop_DubSin (Arb cliffor) = sin (2 * cliffor) ≈≈ 2 * sin cliffor * cos cliffor + +prop_DubCos :: ArbCl3 -> Bool +prop_DubCos (Arb cliffor) = cos (2 * cliffor) ≈≈ cos cliffor ^ 2 - sin cliffor ^ 2 + +prop_DubTan :: ArbCl3 -> Bool +prop_DubTan (Arb cliffor) = tan (2 * cliffor) ≈≈ (2 * tan cliffor) / (1 - tan cliffor ^ 2) + +prop_DubSinh :: ArbCl3 -> Bool +prop_DubSinh (Arb cliffor) = sinh (2 * cliffor) ≈≈ 2 * sinh cliffor * cosh cliffor + +prop_DubCosh :: ArbCl3 -> Bool +prop_DubCosh (Arb cliffor) = cosh (2 * cliffor) ≈≈ 2 * cosh cliffor ^ 2 - 1 + +prop_DubTanh :: ArbCl3 -> Bool +prop_DubTanh (Arb cliffor) = tanh (2 * cliffor) ≈≈ (2 * tanh cliffor) / (1 + tanh cliffor ^ 2) + +prop_PosSinShift :: ArbCl3 -> Bool +prop_PosSinShift (Arb cliffor) = sin (pi/2 + cliffor) ≈≈ cos cliffor + +prop_NegSinShift :: ArbCl3 -> Bool +prop_NegSinShift (Arb cliffor) = sin (pi/2 - cliffor) ≈≈ cos cliffor + +prop_SinSqCosSq :: ArbCl3 -> Bool +prop_SinSqCosSq (Arb cliffor) = (abs cliffor > 10) || (sin cliffor ^ 2 + cos cliffor ^ 2 ≈≈ 1) + +prop_CoshSqmSinhSq :: ArbCl3 -> Bool +prop_CoshSqmSinhSq (Arb cliffor) = (abs cliffor > 10) || (cosh cliffor ^ 2 - sinh cliffor ^ 2 ≈≈ 1) + +prop_SymCosh :: ArbCl3 -> Bool +prop_SymCosh (Arb cliffor) = cosh (negate cliffor) ≈≈ cosh cliffor + +prop_SymSinh :: ArbCl3 -> Bool +prop_SymSinh (Arb cliffor) = sinh (negate cliffor) ≈≈ negate (sinh cliffor) + +prop_DoubleISin :: ArbCl3 -> Bool +prop_DoubleISin (Arb cliffor) = 2 * I 1 * sin cliffor ≈≈ exp(I 1 * cliffor) - exp (I (-1) * cliffor) + +-- | Composition Sub-Algebras have a distributive norm over multiplication, like this: +-- +-- > norm $ clif * clif' = norm clif * norm clif' +-- +-- Strangly the constructor combinations with the "= True" don't play nice with 'abs' +-- they are the constructors with non-zero zero-divisors. +prop_CompAlg :: (ArbCl3, ArbCl3) -> Bool +prop_CompAlg (Arb PV{}, Arb PV{}) = True +prop_CompAlg (Arb PV{}, Arb BPV{}) = True +prop_CompAlg (Arb PV{}, Arb TPV{}) = True +prop_CompAlg (Arb PV{}, Arb APS{}) = True +prop_CompAlg (Arb BPV{}, Arb PV{}) = True +prop_CompAlg (Arb TPV{}, Arb PV{}) = True +prop_CompAlg (Arb APS{}, Arb PV{}) = True +prop_CompAlg (Arb BPV{}, Arb BPV{}) = True +prop_CompAlg (Arb BPV{}, Arb TPV{}) = True +prop_CompAlg (Arb BPV{}, Arb APS{}) = True +prop_CompAlg (Arb TPV{}, Arb BPV{}) = True +prop_CompAlg (Arb APS{}, Arb BPV{}) = True +prop_CompAlg (Arb TPV{}, Arb TPV{}) = True +prop_CompAlg (Arb TPV{}, Arb APS{}) = True +prop_CompAlg (Arb APS{}, Arb TPV{}) = True +prop_CompAlg (Arb APS{}, Arb APS{}) = True +prop_CompAlg (Arb cliffor, Arb cliffor') = abs ( cliffor * cliffor') ≈≈ abs cliffor * abs cliffor' + + +-- Run the test +-- compare the function under test (fUT) to a golden test funcion (gTF) +-- if the input is within bounds +runTest :: String -> (Cl3 -> Cl3) -> (Cl3 -> Cl3) -> (Cl3 -> Bool) -> Cl3 -> IO() +runTest testName fUT gTF outOB iVal = + let f = fUT iVal + g = gTF iVal + in if outOB iVal + then putStr (unlines [testName ++ ": Input Out of Bounds"]) + else if f ≈≈ g + then putStr (unlines [testName ++ ": Passed"]) + else putStr (unlines [testName ++ ": Failed" + ,"Expected: " ++ show g + ," got: " ++ show f + ,"on input: " ++ show iVal + ]) + +-- | '≈≈' aproximately equal +(≈≈) :: Cl3 -> Cl3 -> Bool +(reduce -> clifforA) ≈≈ (reduce -> clifforB) = + let ave = (abs clifforA + abs clifforB) / 2 + in abs (clifforA - clifforB) <= 1e-5*ave + tol +infix 4 ≈≈ + +-- | 'poles' a function that tests if a cliffor is one of the defined poles +poles :: [Cl3] -> Cl3 -> Bool +poles [] _ = False +poles [p] cliffor = eig1 ≈≈ p || eig2 ≈≈ p + where (eig1,eig2) = eigvals cliffor +poles (p:ps) cliffor = (eig1 ≈≈ p || eig2 ≈≈ p) || poles ps cliffor + where (eig1,eig2) = eigvals cliffor + +unWrapIPartEigs :: Cl3 -> Cl3 +unWrapIPartEigs cliffor = reduce $ spectraldcmp unWrapI id cliffor + where unWrapI (R a0) = R a0 + unWrapI (I a123) | a123 > pi = unWrapI $ I (a123 - 2*pi) + | a123 < (-pi) = unWrapI $ I (a123 + 2*pi) + | otherwise = I a123 + unWrapI (C a0 a123) | a123 > pi = unWrapI $ C a0 (a123 - 2*pi) + | a123 < (-pi) = unWrapI $ C a0 (a123 + 2*pi) + | otherwise = C a0 a123 + unWrapI _ = error "unWrapI should only be unWrapping R I and C" + +---------------------------------------------------------- + +------------------------------------------------------------------- +-- +-- Arbitrary Instance of Cl3 types, typically for use with the +-- "Test.QuickCheck" library. +-- +------------------------------------------------------------------- + +-- | 'ArbCl3' to provide a newtype wrapper to avoid the orphan instance +newtype ArbCl3 = Arb Cl3 deriving (Show) + +-- | 'Arbitrary' instance that has its largest singular value less than or equal to 15 +instance Arbitrary ArbCl3 where + arbitrary = + oneof [(Arb.)R <$> arbitrary, + ((Arb.).).V3 <$> arbitrary <*> arbitrary <*> arbitrary, + ((Arb.).).BV <$> arbitrary <*> arbitrary <*> arbitrary, + (Arb.)I <$> arbitrary, + (((Arb.).).).PV <$> arbitrary <*> arbitrary <*> arbitrary <*> arbitrary, + (((Arb.).).).H <$> arbitrary <*> arbitrary <*> arbitrary <*> arbitrary, + (Arb.).C <$> arbitrary <*> arbitrary, + (((((Arb.).).).).).BPV <$> arbitrary <*> arbitrary <*> arbitrary + <*> arbitrary <*> arbitrary <*> arbitrary, + (((Arb.).).).ODD <$> arbitrary <*> arbitrary <*> arbitrary <*> arbitrary, + (((Arb.).).).TPV <$> arbitrary <*> arbitrary <*> arbitrary <*> arbitrary, + (((((((Arb.).).).).).).).APS <$> arbitrary <*> arbitrary <*> arbitrary <*> arbitrary + <*> arbitrary <*> arbitrary <*> arbitrary <*> arbitrary + ] `suchThat` lessThan15 + where + lessThan15 (Arb cliffor) = abs cliffor <= 15