diff --git a/ChangeLog.md b/ChangeLog.md
--- a/ChangeLog.md
+++ b/ChangeLog.md
@@ -1,5 +1,29 @@
 # Revision history for cl3
 
+## 3.0.0.0  -- 2022-04-01
+
+* Added LiquidHaskell support
+* Added more memory efficient subtypes for R, V3, BV, I, PV, H, C, BPV, ODD, TPV
+* Added LiquidHaskell support: liquidhaskell >= 0.8.10
+* Sucessfully tested with: liquidhaskell-0.8.10.2 and Stackage lts-18.6
+* Added a LiquidHaskell flag to build: stack build --flag cl3:do-liquid
+* LiquidHaskell requires CPP Options: -DO_LIQUID -DO_NO_RANDOM -DO_NO_DERIVED -DO_NO_STORABLE
+* LiquidHaskell _|_ when: Random, Storable, or Derived instances are enabled
+* LiquidHaskell _|_ when: There is any number n.nn*e^-x see [Issue#1762](https://github.com/ucsd-progsys/liquidhaskell/issues/1762) for my 128*eps threshold, eventhough rationals are supposed to be supported
+* Added more memory efficient subtypes for R, V3, BV, I, PV, H, C, BPV, ODD, TPV, and APS
+* Namely: Cl3_R, Cl3_V3, Cl3_BV, Cl3_I, Cl3_PV, Cl3_H, Cl3_C, Cl3_BPV, Cl3_ODD, Cl3_TPV, and Cl3_APS
+* Added smart constructors to covert to/from the more memory efficient subtypes
+* Added a "weigh" based benchmark to verify that the more memory efficient subtypes are actually more memory efficent
+* Added a WeighStorableCl3 benchmark; command to run: stack bench cl3:bench-cl3-weigh
+* Added a NbodyMassiv benchmark; command to run: stack bench cl3:bench-cl3-massiv-nbody
+* Modified existing NbodyGameCl3 benchmark; command to run: stack bench cl3:bench-cl3-nbody
+* LiquidHaskell found some interesting unsaftey in the definition of 'tan' and 'tanh'
+* Hopefully it's solved with a judicious use of 'reduce' for instances above R and I, in the affected 'tan' and 'tanh' calls
+* Corrected spelling error sp: simi to semi, in the definition of a random unit vector
+* 'abssignum' why not calculate 'abs' once and use it twice
+* 'absolute' only constructs an R
+
+
 ## 2.0.0.0  -- 2020-06-20
 
 * Added work around for GHC 8.10 regression of Issue #15304 reproducing code changes from GHC MR 2608 in the source files
diff --git a/README.md b/README.md
--- a/README.md
+++ b/README.md
@@ -121,7 +121,11 @@
 
 # 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.
+On my machine with GHC-8.10.7 the current fastest implementation completes 50M steps with a mean time of 4.014 seconds.  The benchmark uses a hand rolled implementation of vector math.  The Cl3 implementation completes 50M steps with a mean time of 5.691 seconds.  This 1.67 second difference amounts to a 33.5 ns difference in the inner loop.  This performance has been degraded with GHC regressions in GHC-9.0.2 and GHC-9.2.2 by ~5x.
+In the 3.0 release a [massiv](https://hackage.haskell.org/package/massiv) benchmark was added in addition to a [weigh](https://hackage.haskell.org/package/weigh) based benchmark.
+
+# Saftey and Correctness
+In the 3.0 release [Liquid Haskell](https://hackage.haskell.org/package/liquidhaskell) support was added, Liquid Haskell did prove its worth by finding a couple of bugs in the implementation.  So far it is an initial release and not much has been done to fully integrate Liquid Haskell to the library.
 
 # 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.
diff --git a/benchmarks/NbodyGameCl3.hs b/benchmarks/NbodyGameCl3.hs
--- a/benchmarks/NbodyGameCl3.hs
+++ b/benchmarks/NbodyGameCl3.hs
@@ -1,6 +1,8 @@
 {-# LANGUAGE ViewPatterns #-}
 {-# LANGUAGE BangPatterns #-}
 {-# LANGUAGE CPP #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE InstanceSigs #-}
 
 #if __GLASGOW_HASKELL__ == 810
 -- Work around to fix GHC Issue #15304, issue popped up again in GHC 8.10, it should be fixed in GHC 8.12
@@ -25,7 +27,7 @@
 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 Algebra.Geometric.Cl3 (Cl3(..), toR, toV3, toCl3_V3, toCl3_R, fromCl3_V3, fromCl3_R, Cl3_V3(), Cl3_R()) -- To add Cl3
 
 
 import Foreign.Ptr (Ptr, castPtr, plusPtr)
@@ -205,10 +207,11 @@
 
 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)]
-                       ]
+          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)]
+                      ,bgroup "nbodyCl3_V3" [ bench (show n) $ nfIO (nbodyCl3_V3 n)]
+                      ]
 
 
 
@@ -216,9 +219,9 @@
 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
+  energy pPlanets >>= print -- 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
+  energy pPlanets >>= print -- printf "%.9f\n"   -- calculate and print out the final energy of the system
 
 
 nbodyCl3 :: (Eq t, Num t) => t -> IO ()
@@ -239,6 +242,14 @@
   energyAPS pPlanetsAPS >>= print.toR -- printf "%0.9f\n"
 
 
+nbodyCl3_V3 :: (Eq t, Num t) => t -> IO ()
+nbodyCl3_V3 n = do
+  pPlanetsCl3_V3 <- fromListCl3_V3 planetsCl3_V3
+  nbodyInitCl3_V3 pPlanetsCl3_V3
+  energyCl3_V3 pPlanetsCl3_V3 >>= print
+  runCl3_V3 n pPlanetsCl3_V3
+  energyCl3_V3 pPlanetsCl3_V3 >>= print
+
 -- | 'run' advances the system 'n' times by recursively calling 'run'
 run :: (Eq t, Num t) => t -> Ptr Planet -> IO ()
 run 0 _ = return ()
@@ -261,6 +272,14 @@
   runAPS (i-1) p
 
 
+-- runCl3_V3 is a tail recursive action to advance the system
+runCl3_V3 :: (Eq t, Num t) => t -> Ptr PlanetCl3_V3 -> IO ()
+runCl3_V3 0 _ = return ()
+runCl3_V3 i p = do
+  advanceCl3_V3 p
+  runCl3_V3 (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.
@@ -279,7 +298,11 @@
 offsetMomentumAPS p momCl3 = p{velCl3 = negate momCl3 / massCl3 p}
 
 
+offsetMomentumCl3_V3 :: PlanetCl3_V3 -> Cl3 -> PlanetCl3_V3
+offsetMomentumCl3_V3 p momCl3 = p{velCl3_V3 = toCl3_V3 (negate momCl3 / (fromCl3_R.massCl3_R $ p))}
 
+
+
 -- | 'nbodyInit' initializes the velocity of the sun so that
 -- the total system momentum is zero.
 nbodyInit :: Ptr Planet -> IO ()
@@ -326,7 +349,21 @@
 -- end of nbodyInitAPS
 
 
+nbodyInitCl3_V3 :: Ptr PlanetCl3_V3 -> IO ()
+nbodyInitCl3_V3 pPlanetsCl3_V3 = do
+  let initCl3_V3 !pCl3_V3 i = -- initCl3_V3 takes a momentum accumulator and an iteration count
+        if i < length planetsCl3_V3  -- if the iteration count is less than the length of the list of bodies; add more momentum
+        then do
+            p <- peekElemOff pPlanetsCl3_V3 i  -- pull off a body at the current iteration count
+            initCl3_V3 (pCl3_V3 + (fromCl3_V3.velCl3_V3 $ p) * (fromCl3_R.massCl3_R $ p)) (i+1)  -- accumulate the momentum in the accumulator, increase the iteration count and recursively call the function again
+        else return pCl3_V3  -- no more body's momentum left to accumulate, so return the total momentum of all of the bodies
+ 
+  totalMomentum <- initCl3_V3 (V3 0 0 0) 0  -- calculate the total momentum, start with a zero momentum vector in the accumulator, and an iteration count of zero, returns a Cl3
+  sun0 <- peek pPlanetsCl3_V3 -- grab the Sun, don't get burned
+  poke pPlanetsCl3_V3 $ offsetMomentumCl3_V3 sun0 totalMomentum  -- update the Sun's velocity with the total momentum, such that the system has zero total momentum
 
+
+
 -- | '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'
@@ -409,7 +446,31 @@
 -- end of energyAPS
 
 
+energyCl3_V3 :: Ptr PlanetCl3_V3 -> IO Cl3
+energyCl3_V3 pPlanetsCl3_V3 = do
+  let energy' e i =
+        if i < length planetsCl3_V3
+        then do
+          p <- peekElemOff pPlanetsCl3_V3 i
+          e1 <- energy'' p (i+1) e
+          e2 <- energy' e (i+1)
+          return $! e + 0.5 * (fromCl3_R.massCl3_R $ p) * (toR $ (fromCl3_V3.velCl3_V3 $ p)^(2 :: Int)) + e1 + e2 -- the Kenitic Energy and other accumulated energies
+        else return e  -- all of the bodies energies have been calculated, return the total energy
+     
+      energy'' p j e = -- a Potential Energy calculation
+        if j < length planetsCl3_V3
+        then do
+          pj <- peekElemOff pPlanetsCl3_V3 j
+          let !distance = abs dpos
+              !dpos = (fromCl3_V3.posCl3_V3 $ pj) - (fromCl3_V3.posCl3_V3 $ p)
+          e1 <- energy'' p (j+1) e
+          return $! e - (((fromCl3_R.massCl3_R $ p) * (fromCl3_R.massCl3_R $ pj)) / distance) + e1
+        else return e
+ 
+  energy' (R 0) 0  -- starts off the recursive calculation of energy, returns the total energy as a Cl3
 
+
+
 -- | '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
@@ -506,8 +567,34 @@
 -- end of advanceAPS
 
 
+advanceCl3_V3 :: Ptr PlanetCl3_V3 -> IO ()
+advanceCl3_V3 pPlanetsCl3_V3 = do
+  let advance' i = when (i < length planetsCl3_V3) $ do  -- loops through all the planets and updates the velocity
+                            let loop j = when (j < length planetsCl3_V3) $ do
+                                      ii <- peekElemOff pPlanetsCl3_V3 i
+                                      jj <- peekElemOff pPlanetsCl3_V3 j
+                                      let -- !mag = R dt / (toR dSquared * abs dpos)
+                                          -- !dSquared = dpos * dpos
+                                          !dpos = (fromCl3_V3.posCl3_V3 $ ii) - (fromCl3_V3.posCl3_V3 $ jj)
+                                          !(R !dMag) = abs dpos
+                                          !mag = R (dt / (dMag^(3 :: Int)))
+                                      pokeVCl3_V3 pPlanetsCl3_V3 i ii{velCl3_V3 = toCl3_V3 $ (fromCl3_V3.velCl3_V3 $ ii) - dpos * (fromCl3_R.massCl3_R $ jj) * mag}
+                                      pokeVCl3_V3 pPlanetsCl3_V3 j jj{velCl3_V3 = toCl3_V3 $ (fromCl3_V3.velCl3_V3 $ jj) + dpos * (fromCl3_R.massCl3_R $ ii) * mag}
+                                      loop (j+1)
+                            loop (i+1)
+                            advance' (i+1)
+     
+      advance'' i = when (i < length planetsCl3_V3) $ do  -- loops through all of the planets and updates the position
+                      p <- peekElemOff pPlanetsCl3_V3 i
+                      pokeCCl3_V3 pPlanetsCl3_V3 i p{posCl3_V3 = toCl3_V3 $ (fromCl3_V3.posCl3_V3 $ p) + R dt * (fromCl3_V3.velCl3_V3 $ p)}
+                      advance'' (i+1)
+ 
+  advance' 0   -- update all of the planets velocities
+  advance'' 0  -- update all of the planets positions
 
 
+
+
 data Planet = Planet {x :: !Double,
                       y :: !Double,
                       z :: !Double,
@@ -525,6 +612,7 @@
                             velCl3 :: !Cl3,
                             massCl3 :: !Cl3} deriving (Show)
 
+
 planetsCl3 :: [PlanetCl3]
 planetsCl3 = [PlanetCl3 sunPos sunVel sunMass,
               PlanetCl3 jupiterPos jupiterVel jupiterMass,
@@ -533,7 +621,21 @@
               PlanetCl3 neptunePos neptuneVel neptuneMass]
 
 
+data PlanetCl3_V3 = PlanetCl3_V3 {posCl3_V3 :: !Cl3_V3
+                                 ,velCl3_V3 :: !Cl3_V3
+                                 ,massCl3_R :: !Cl3_R}
 
+instance Show PlanetCl3_V3 where
+  show p = "PlanetCl3_V3 " ++ show (show.fromCl3_V3.posCl3_V3 $ p) ++ show (show.fromCl3_V3.velCl3_V3 $ p) ++ show (show.fromCl3_R.massCl3_R $ p)
+
+planetsCl3_V3 :: [PlanetCl3_V3]
+planetsCl3_V3 = [PlanetCl3_V3 (toCl3_V3 sunPos) (toCl3_V3 sunVel) (toCl3_R sunMass)
+                ,PlanetCl3_V3 (toCl3_V3 jupiterPos) (toCl3_V3 jupiterVel) (toCl3_R jupiterMass)
+                ,PlanetCl3_V3 (toCl3_V3 saturnPos) (toCl3_V3 saturnVel) (toCl3_R saturnMass)
+                ,PlanetCl3_V3 (toCl3_V3 uranusPos) (toCl3_V3 uranusVel) (toCl3_R uranusMass)
+                ,PlanetCl3_V3 (toCl3_V3 neptunePos) (toCl3_V3 neptuneVel) (toCl3_R neptuneMass)]
+
+
 sun :: Planet
 sun = Planet {x = 0, y = 0, z = 0,
               vx = 0, vy = 0, vz = 0,
@@ -630,6 +732,7 @@
   alignment _ = dblSz
   peekElemOff p i = peek (plusPtr p (i * sizeOf (undefined :: Planet)))
   pokeElemOff p i e = poke (plusPtr p (i * sizeOf e)) e
+  peek :: Ptr Planet -> IO Planet
   peek p = do
     x' <- peek (offset 0)
     y' <- peek (offset 1)
@@ -641,6 +744,7 @@
     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 :: Ptr Planet -> Planet -> IO ()
   poke p e = do
     poke (offset 0) $ x e
     poke (offset 1) $ y e
@@ -662,6 +766,7 @@
   alignment _ = cl3Sz
   peekElemOff p i = peek ( p `plusPtr` (i * sizeOf (undefined :: PlanetCl3)))
   pokeElemOff p i e = poke (p `plusPtr` (i * sizeOf e)) e
+  peek :: Ptr PlanetCl3 -> IO PlanetCl3
   peek p = do
     pos' <- peek (offset 0)
     vel' <- peek (offset 1)
@@ -669,6 +774,7 @@
     return $ PlanetCl3 pos' vel' mass'
       where
         offset i = (castPtr p :: Ptr Cl3) `plusPtr` (i * cl3Sz)
+  poke :: Ptr PlanetCl3 -> PlanetCl3 -> IO ()
   poke p e = do
     poke (offset 0) $! posCl3 e
     poke (offset 1) $! velCl3 e
@@ -680,6 +786,27 @@
 cl3Sz = sizeOf (undefined :: Cl3)
 
 
+-- | 'Storable' instance of PlanetCl3_V3
+instance Storable PlanetCl3_V3 where
+  sizeOf _ = sizeOf (undefined :: Cl3_V3) + sizeOf (undefined :: Cl3_V3) + sizeOf (undefined :: Cl3_R)
+  alignment _ = undefined  -- perhaps this should be some integer power of 2, sizeOf PlanetCl3_V3 is 54bytes, should I pick 64bytes? Does Haskell even support a mov64?, or should it be 8 because everything at the core is a Double?
+  peekElemOff ptr idx = peek (castPtr ptr `plusPtr` (idx * sizeOf (undefined :: PlanetCl3_V3)))
+  pokeElemOff ptr idx e = poke (castPtr ptr `plusPtr` (idx * sizeOf (undefined :: PlanetCl3_V3))) e
+  peek :: Ptr PlanetCl3_V3 -> IO PlanetCl3_V3 
+  peek ptr = do
+    pos :: Cl3_V3 <- peek (castPtr ptr `plusPtr` 0)
+    vel :: Cl3_V3 <- peek (castPtr ptr `plusPtr` (sizeOf pos))
+    m :: Cl3_R <- peek (castPtr ptr `plusPtr` (sizeOf pos) `plusPtr` (sizeOf vel))
+    return (PlanetCl3_V3 pos vel m)
+  poke :: Ptr PlanetCl3_V3 -> PlanetCl3_V3 -> IO ()
+  poke ptr (PlanetCl3_V3 pos vel m) = do
+    poke (castPtr ptr `plusPtr` 0) pos
+    poke (castPtr ptr `plusPtr` (sizeOf pos)) vel
+    poke (castPtr ptr `plusPtr` (sizeOf pos) `plusPtr` (sizeOf vel)) m
+
+
+
+
 -- | 'pokeC' update the position
 pokeC :: Ptr Planet -> Int -> Planet -> IO ()
 pokeC p i e = do
@@ -694,7 +821,12 @@
 pokeCCl3 p !i !e =
   poke ((castPtr p :: Ptr Cl3) `plusPtr` (i * sizeOf (undefined :: PlanetCl3))) $! posCl3 e
 
+-- Update the position of a planet at an index in the array of planets
+pokeCCl3_V3 :: Ptr PlanetCl3_V3 -> Int -> PlanetCl3_V3 -> IO ()
+pokeCCl3_V3 p idx e =
+  poke ((castPtr p) `plusPtr` (idx * sizeOf (undefined :: PlanetCl3_V3))) (posCl3_V3 e)
 
+
 -- | 'pokeV' update the velocity
 pokeV :: Ptr Planet -> Int -> Planet -> IO ()
 pokeV p i e = do
@@ -709,7 +841,12 @@
 pokeVCl3 p !i !e =
   poke ((castPtr p :: Ptr Cl3) `plusPtr` (cl3Sz + i * sizeOf (undefined :: PlanetCl3))) $! velCl3 e
 
+-- Update the velocity of a planet at an index in the array of planets
+pokeVCl3_V3 :: Ptr PlanetCl3_V3 -> Int -> PlanetCl3_V3 -> IO ()
+pokeVCl3_V3 p idx e =
+  poke ((castPtr p) `plusPtr` (idx * sizeOf (undefined :: PlanetCl3_V3)) `plusPtr` (sizeOf (undefined :: Cl3_V3))) (velCl3_V3 e)
 
+
 -- | 'fromList' initialize the planets in the stack
 fromList :: [Planet] -> IO (Ptr Planet)
 fromList l = do
@@ -731,6 +868,18 @@
       loop (q:qs) i = do
         poke (pa `plusPtr` (i * sizeOf (undefined :: PlanetCl3))) q
         loop qs (i+1)
+  loop l 0
+  return pa
+
+
+fromListCl3_V3 :: [PlanetCl3_V3] -> IO (Ptr PlanetCl3_V3)
+fromListCl3_V3 l = do
+  let len = length l
+  pa <- mallocBytes (len * sizeOf (undefined :: PlanetCl3_V3))
+  let loop [] _ = return ()
+      loop (q:qs) idx = do
+        poke (pa `plusPtr` (idx * sizeOf (undefined :: PlanetCl3_V3))) q
+        loop qs (idx+1)
   loop l 0
   return pa
 
diff --git a/benchmarks/NbodyMassiv.hs b/benchmarks/NbodyMassiv.hs
new file mode 100644
--- /dev/null
+++ b/benchmarks/NbodyMassiv.hs
@@ -0,0 +1,293 @@
+{-# LANGUAGE ViewPatterns #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE PartialTypeSignatures #-}
+{-# LANGUAGE InstanceSigs #-}
+
+
+
+import Algebra.Geometric.Cl3
+import Data.Massiv.Array as A
+import Control.Monad (replicateM_,when)
+import Data.Time.Clock.POSIX (getPOSIXTime)
+
+import Foreign.Ptr (Ptr, castPtr, plusPtr)
+import Foreign.Storable (Storable(..))
+
+import Data.IORef (IORef(),newIORef,modifyIORef',readIORef)
+
+
+main :: IO ()
+main = do
+  arrayProgrammingAlgorithm
+  bitTwiddlingAlgorithm
+  
+arrayProgrammingAlgorithm :: IO ()
+arrayProgrammingAlgorithm = 
+  do tStart <- getPOSIXTime
+     n <- return (50000000 :: Int)
+     planetsPos :: Array S Ix1 Cl3_V3 <- fromListsM Seq planetspos  -- generate immutible Massiv array of the planets position
+     planetsVel :: Array S Ix1 Cl3_V3 <- fromListsM Seq planetsvel  -- generate immutible Massiv array of the planets velocity
+     planetsMass :: Array S Ix1 Cl3_R <- fromListsM Seq planetsmass  -- generate immutible Massiv array of the planets mass
+     mutablePos <- thaw planetsPos  -- copy to a mutible Massiv array of the planets position
+     mutableVel <- thaw planetsVel  -- copy to a mutible Massiv array of the planets velocity
+     let (Sz1 c) = size planetsVel  -- Size of the velocity array
+         c' = (c^(2 :: Int) - c) `div` 2  -- size of the nChoose2With array
+     mutableScratch <- newMArray (Sz1 c') (toCl3_V3 $ V3 0 0 0) -- nChoose2With sized scratch pad
+     mutableScratch2D <- newMArray (Sz (c :. (c-1))) (toCl3_V3 $ V3 0 0 0) -- 2D scratch pad for holding acceleration components
+     offsetMomentum mutableVel planetsMass -- modify the Sun's velocity such that the total momentum is zero
+     printEnergy mutablePos mutableVel planetsMass
+     replicateM_ n (advance mutablePos mutableVel planetsMass mutableScratch mutableScratch2D)
+     printEnergy mutablePos mutableVel planetsMass
+     tEnd <- getPOSIXTime
+     print $ show (tEnd - tStart) ++ " (sec) for " ++ show n ++ " iterations." -- Relative time measurement
+
+-- | 'bitTwiddlingAlgorithm' with Ϟ symbolizing the bit twiddling algorithm versions
+bitTwiddlingAlgorithm :: IO ()
+bitTwiddlingAlgorithm = 
+  do tStart <- getPOSIXTime
+     n <- return (50000000 :: Int)
+     (solarSystem,sz) <- offsetMomentumϞ [sun,jupiter,saturn,uranus,neptune] -- returns a zero momentum solar system
+     printEnergyϞ solarSystem sz
+     replicateM_ n (advanceϞ solarSystem sz)
+     printEnergyϞ solarSystem sz
+     tEnd <- getPOSIXTime
+     print $ show (tEnd - tStart) ++ " (sec) for " ++ show n ++ " iterations." -- Relative time measurement
+
+offsetMomentumϞ :: [Body] -> IO (SolarSystem,Int) -- A SolarSystem and the Size : size: sz :: Int
+offsetMomentumϞ bodies = do
+  initBodies :: Array S Ix1 Body <- A.fromListsM Seq bodies
+  -- mutate the solar system so that the sun's velocity gives zero total momentum
+  let velocity :: Array D Ix1 Cl3 = A.map (fromCl3_V3.vel) initBodies
+      mass :: Array D Ix1 Cl3 = A.map (fromCl3_R.mas) initBodies
+  ss :: SolarSystem <- thaw initBodies
+  let sz = totalElem.sizeOfMArray $ ss
+  bs <- readM ss 0
+  let sunV = fromCl3_V3.vel $ bs
+      sunM = fromCl3_R.mas $ bs
+  writeM ss 0 (bs{vel = toCl3_V3 (sunV - ((velocity `dot_V3` mass) / sunM))})
+  return (ss,sz)
+
+
+printEnergyϞ :: SolarSystem -> Int -> IO ()
+printEnergyϞ ss sz = do
+  acc :: IORef Cl3 <- newIORef (R 0)
+  ssFroz <- freeze Seq ss
+  iforIO_ ssFroz $ \idx bi ->
+    do let bi_pos :: Cl3 = fromCl3_V3.pos $ bi
+           bi_vel :: Cl3 = fromCl3_V3.vel $ bi
+           bi_mas :: Cl3 = fromCl3_R.mas $ bi
+           ke :: Cl3 = 0.5 * bi_mas * (toR $ bi_vel^(2::Int))
+       modifyIORef' acc (+ ke)  -- accumulate the kenetic energy for each body
+       -- now calcualate the potential energy
+       let szMinus1 = sz-1
+       when (idx < szMinus1) $ do
+         brest <- extractM (idx+1) (Sz1 (szMinus1-idx)) ssFroz
+         forIO_ brest $ \bj ->
+           do let bj_pos :: Cl3 = fromCl3_V3.pos $ bj
+                  bj_mas :: Cl3 = fromCl3_R.mas $ bj
+                  distance :: Cl3 = sqrt.toR $ (bi_pos - bj_pos)^(2::Int)
+                  pe_ij :: Cl3 = (bi_mas * bj_mas) / distance
+              modifyIORef' acc (\ac -> ac - pe_ij)
+  energy <- readIORef acc
+  print energy
+
+advanceϞ :: SolarSystem -> Int -> IO ()
+advanceϞ ss sz = do
+  updateVelocityϞ ss sz
+  updatePositionϞ ss
+
+updateVelocityϞ :: SolarSystem -> Int -> IO ()
+updateVelocityϞ ss sz = do
+  ssFroz <- freeze Seq ss
+  iforIO_ ssFroz $ \idx bi ->
+    when (idx < sz-1) $ do
+      let bi_pos :: Cl3 = fromCl3_V3.pos $ bi
+          bi_mas :: Cl3 = fromCl3_R.mas $ bi
+      brest <- extractM (idx+1) (Sz1 (sz-idx-1)) ssFroz
+      iforIO_ brest $ \jdx bj ->
+        do let bj_pos :: Cl3 = fromCl3_V3.pos $ bj
+               bj_mas :: Cl3 = fromCl3_R.mas $ bj
+               (absMag,dposUnit) = abssignum $ bi_pos - bj_pos
+               mag = dt / ((absMag)^(2::Int))
+               aij = negate $ dposUnit * bj_mas * mag
+               aji = dposUnit * bi_mas * mag
+           modify_ ss (modVel aij) idx
+           modify_ ss (modVel aji) (idx+jdx+1)
+
+modVel :: PrimMonad m => Cl3 -> Body -> m Body
+modVel a bod = do
+  let vbod = fromCl3_V3 $ vel bod
+  return bod{vel = toCl3_V3 $ vbod + a}
+
+updatePositionϞ :: SolarSystem -> IO ()
+updatePositionϞ ss = do
+  iforPrimM_ ss $ \idx b ->
+    do let p :: Cl3 = fromCl3_V3.pos $ b
+           v :: Cl3 = fromCl3_V3.vel $ b
+       writeM ss idx b{pos = toCl3_V3 (p + dt * v)}
+        
+
+advance :: MArray RealWorld S Ix1 Cl3_V3 -> MArray RealWorld S Ix1 Cl3_V3 -> Array S Ix1 Cl3_R -> MArray RealWorld S Ix1 Cl3_V3 -> MArray RealWorld S Ix2 Cl3_V3 -> IO ()
+advance sysPos sysVel sysMass scratchPad scratchPad2D =
+  do updateVelocity sysPos sysVel sysMass scratchPad scratchPad2D
+     updatePosition sysPos sysVel
+
+updateVelocity :: MArray RealWorld S Ix1 Cl3_V3 -> MArray RealWorld S Ix1 Cl3_V3 -> Array S Ix1 Cl3_R -> MArray RealWorld S Ix1 Cl3_V3 -> MArray RealWorld S Ix2 Cl3_V3 -> IO ()
+updateVelocity sysPos sysVel sysMass scratchPad scratchPad2D =
+  do -- Calculate the inverse square distance and unit vector between the bodies
+     frozPos <- freeze Seq sysPos
+     let ps :: Array D Ix1 Cl3 = A.map fromCl3_V3 frozPos
+         scaledPij = nChoose2With (\ i j -> let (a,s) = abssignum $ i - j in (recip (a^(2 :: Int))) * s) ps -- probably should use a computeSource here, because scaledPij is used twice, or verify streams are still working.
+     computeInto scratchPad (A.map toCl3_V3 scaledPij)  -- Compute result to the scratch pad
+     frozScratch <- freeze Seq scratchPad  -- Freeze the result for future use
+     -- Scale the previous inverse square distance by the mass of the other body
+     let scratch :: Array D Ix1 Cl3 = A.map fromCl3_V3 frozScratch
+         ms :: Array D Ix1 Cl3 = A.map fromCl3_R sysMass
+         upper = nChoose2With (\ _ j -> negate j) ms !*! scratch
+         lower = nChoose2With (\ i _ -> i) ms !*! scratch
+     computeInto scratchPad (A.map toCl3_V3 upper)
+     frozScratchUpper <- freeze Seq scratchPad
+     computeInto scratchPad (A.map toCl3_V3 lower)
+     frozScratchLower <- freeze Seq scratchPad
+     -- Load each body acceleration contribution into a 2D array 
+     frozVel <- freeze Seq sysVel
+     let vs :: Array D Ix1 Cl3 = A.map fromCl3_V3 frozVel
+         rows = elemsCount vs
+     iforPrimM_ scratchPad2D (\ix _ -> write_ scratchPad2D ix (genAcc rows frozScratchUpper frozScratchLower ix))
+     frozScratch2D <- freeze Seq scratchPad2D
+     --  Sum each row to calcuate the total acceleration and then update the velocity
+     let accComponents :: Array D Ix2 Cl3 = A.map fromCl3_V3 frozScratch2D
+         as = makeArray Seq (size vs) (\ ix -> foldlS (+) (V3 0 0 0) (accComponents !> ix))
+         vs' = vs !+! (dt *. as)
+     computeInto sysVel (A.map toCl3_V3 vs')
+
+genAcc :: Int -> Array S Ix1 Cl3_V3 -> Array S Ix1 Cl3_V3 -> Ix2 -> Cl3_V3
+genAcc rows upper lower (row :. col) | row <= col = upper ! unSplitIdx rows (row,col+1)
+                                     | otherwise  = lower ! unSplitIdx rows (col,row)
+
+unSplitIdx :: Int -> (Int,Int) -> Int
+unSplitIdx c (row,col) = go 0 (c-1) 0 + (col-1 - row)
+  where
+    go n c1 acc | n == row = acc
+                | otherwise = go (n+1) (c1-1) (acc+c1)
+
+
+
+updatePosition :: MArray RealWorld S Ix1 Cl3_V3 -> MArray RealWorld S Ix1 Cl3_V3 -> IO ()
+updatePosition sysPos sysVel =
+  do frozPos <- freeze Seq sysPos
+     frozVel <- freeze Seq sysVel
+     let ps :: Array D Ix1 Cl3 = A.map fromCl3_V3 frozPos
+         vs :: Array D Ix1 Cl3 = A.map fromCl3_V3 frozVel
+         ps' = ps !+! (dt *. vs)
+     computeInto sysPos (A.map toCl3_V3 ps')
+
+
+offsetMomentum :: MArray RealWorld S Ix1 Cl3_V3 -> Array S Ix1 Cl3_R -> IO ()
+offsetMomentum sysVel sysMass =
+  do frozVel <- freeze Seq sysVel
+     planetsvs :: Array D Ix1 Cl3 <- extractM 1 4 (A.map fromCl3_V3 frozVel)  -- select just the planets and not the sun
+     planetms :: Array D Ix1 Cl3 <- extractM 1 4 (A.map fromCl3_R sysMass)  -- select just the planets and not the sun
+     (fromCl3_R -> sunMassCl3) <- sysMass !? 0
+     writeM sysVel 0 (toCl3_V3 $ negate (planetsvs `dot_V3` planetms) / sunMassCl3)
+     
+
+-- dot_V3 is to work around the typical sum in the dot product that has an aculumulator of a Real 0 value
+dot_V3 :: Array D Ix1 Cl3 -> Array D Ix1 Cl3 -> Cl3
+dot_V3 v1 v2 = foldlS (+) (V3 0 0 0) (A.zipWith (*) v1 v2)
+
+
+printEnergy :: MArray RealWorld S Ix1 Cl3_V3 -> MArray RealWorld S Ix1 Cl3_V3 -> Array S Ix1 Cl3_R -> IO ()
+printEnergy sysPos sysVel sysMass =
+  do frozPos <- freeze Seq sysPos
+     frozVel <- freeze Seq sysVel
+     let ps :: Array D Ix1 Cl3 = A.map fromCl3_V3 frozPos
+         vs :: Array D Ix1 Cl3 = A.map fromCl3_V3 frozVel
+         ms :: Array D Ix1 Cl3 = A.map fromCl3_R sysMass
+         ke = 0.5 * (ms !.! A.map toR (vs.^2)) -- Kenetic Energy O(n)
+         pe = negate.A.sum $ nChoose2With (*) ms !/! nChoose2With (\ x y -> abs $ x - y) ps -- Gravitational Potential Energy O((n^2 - n)/2)
+     print $ ke + pe -- print the total energy
+
+-- This function is a bit of a combinatorial explosion so it is marked to be evaluated in Seqallel, probably should convert to a sgenerate :: Sz1 -> (Ix1 -> e) -> Vector DS e 
+nChoose2With :: (Cl3 -> Cl3 -> Cl3) -> Array D Ix1 Cl3 -> Array D Ix1 Cl3
+nChoose2With f bs =
+  let (Sz1 c) = size bs
+      c' = (c^(2 :: Int) - c) `div` 2
+      bs' :: Array BN Ix1 Cl3 = computeSource bs
+  in makeArray Seq (Sz1 c') (genElements f bs' c)
+
+genElements :: (Cl3 -> Cl3 -> Cl3) -> Array BN Ix1 Cl3 -> Int -> Ix1 -> Cl3
+genElements f bs c idx =
+  let (ix1,ix2) = splitIdx c idx
+  in f (bs ! ix1) (bs ! ix2)
+
+splitIdx :: Int -> Int -> (Int,Int)
+splitIdx c = go (0,1) (c-2)
+  where
+    go (row,col) x ix | ix - x <= 0 = (row, col + ix)
+                      | otherwise = go (row+1,col+1) (x-1) (ix-x-1)
+
+planetspos :: [Cl3_V3]
+planetspos = pos <$> [sun,jupiter,saturn,uranus,neptune]
+
+planetsvel :: [Cl3_V3]
+planetsvel = vel <$> [sun,jupiter,saturn,uranus,neptune]
+
+planetsmass :: [Cl3_R]
+planetsmass = mas <$> [sun,jupiter,saturn,uranus,neptune]
+
+type SolarSystem = MArray RealWorld S Ix1 Body
+
+data Body = Body
+  {pos :: Cl3_V3
+  ,vel :: Cl3_V3
+  ,mas :: Cl3_R}
+
+instance Storable Body where
+  sizeOf _ = sizeOf (undefined :: Cl3_V3) + sizeOf (undefined :: Cl3_V3) + sizeOf (undefined :: Cl3_R)
+  alignment _ = sizeOf (undefined :: Double) -- alignment must be a 2^x
+  peekElemOff ptr idx = peek (castPtr ptr `plusPtr` (idx * sizeOf (undefined :: Body)))
+  pokeElemOff ptr idx e = poke (castPtr ptr `plusPtr` (idx * sizeOf (undefined :: Body))) e
+  peek :: Ptr Body -> IO Body
+  peek ptr = do
+    p :: Cl3_V3 <- peek (castPtr ptr `plusPtr` 0)
+    v :: Cl3_V3 <- peek (castPtr ptr `plusPtr` (sizeOf p))
+    m :: Cl3_R <- peek (castPtr ptr `plusPtr` (sizeOf p) `plusPtr` (sizeOf v))
+    return (Body p v m)
+  poke :: Ptr Body -> Body -> IO ()
+  poke ptr (Body p v m) = do
+    poke (castPtr ptr `plusPtr` 0) p
+    poke (castPtr ptr `plusPtr` (sizeOf p)) v
+    poke (castPtr ptr `plusPtr` (sizeOf p) `plusPtr` (sizeOf v)) m
+
+sun :: Body
+sun = Body (toCl3_V3 $ V3 0 0 0) -- Sun Position
+           (toCl3_V3 $ V3 0 0 0) -- Sun Velocity 
+           (toCl3_R $ R solarMass) -- Sun Mass
+jupiter :: Body
+jupiter = Body (toCl3_V3 $ V3 (4.84143144246472090e+00) (-1.16032004402742839e+00) (-1.03622044471123109e-01))
+               (toCl3_V3 $ V3 (1.66007664274403694e-03*dp) (7.69901118419740425e-03*dp) (-6.90460016972063023e-05*dp))
+               (toCl3_R $ R (9.54791938424326609e-04 * solarMass))
+saturn :: Body
+saturn = Body (toCl3_V3 $ V3 (8.34336671824457987e+00) (4.12479856412430479e+00) (-4.03523417114321381e-01))
+              (toCl3_V3 $ V3 (-2.76742510726862411e-03*dp) (4.99852801234917238e-03*dp) (2.30417297573763929e-05*dp))
+              (toCl3_R $ R (2.85885980666130812e-04 * solarMass))
+uranus :: Body
+uranus = Body (toCl3_V3 $ V3 (1.28943695621391310e+01) (-1.51111514016986312e+01) (-2.23307578892655734e-01))
+              (toCl3_V3 $ V3 (2.96460137564761618e-03*dp) (2.37847173959480950e-03*dp) (-2.96589568540237556e-05*dp))
+              (toCl3_R $ R (4.36624404335156298e-05 * solarMass))
+neptune :: Body
+neptune = Body (toCl3_V3 $ V3 (1.53796971148509165e+01) (-2.59193146099879641e+01) (1.79258772950371181e-01))
+               (toCl3_V3 $ V3 (2.68067772490389322e-03*dp) (1.62824170038242295e-03*dp) (-9.51592254519715870e-05*dp))
+               (toCl3_R $ R (5.15138902046611451e-05 * solarMass))
+
+
+daysPerYear :: Double
+daysPerYear = 365.24
+
+solarMass :: Double
+solarMass = 4 * pi^(2 :: Int)
+dp :: Double
+dp = daysPerYear
+dt :: Cl3
+dt = R 0.01
diff --git a/benchmarks/WeighStorableCl3.hs b/benchmarks/WeighStorableCl3.hs
new file mode 100644
--- /dev/null
+++ b/benchmarks/WeighStorableCl3.hs
@@ -0,0 +1,40 @@
+
+import Weigh
+import Data.Vector.Storable as V
+import Data.Massiv.Array as A
+import Algebra.Geometric.Cl3
+import Foreign.Ptr (Ptr(), plusPtr)
+import Foreign.Storable (Storable(..))
+import Foreign.Marshal.Alloc (mallocBytes)
+
+main :: IO ()
+main = mainWith $ do func' "Cl3.R in Vector" apsR unitR
+                     func' "Cl3.Cl3_R in Vector" compactR unitR
+                     func' "Cl3.R in Massiv Vector" apsRMassiv unitR
+                     func' "Cl3.Cl3_R in Massiv Vector" compactRMassiv unitR
+                     action "malloc Storable Cl3_R by (1024*1024)" $ mallocR unitR
+
+apsR :: Cl3 -> V.Vector Cl3
+apsR cl3 = V.replicate (1024*1204) cl3
+
+compactR :: Cl3 -> V.Vector Cl3_R
+compactR cl3 = V.replicate (1024*1204) (toCl3_R cl3)
+
+apsRMassiv :: Cl3 -> A.Array A.S A.Ix1 Cl3
+apsRMassiv cl3 = A.makeVectorR A.S A.Seq (A.Sz1 (1024*1204)) (const cl3)
+
+compactRMassiv :: Cl3 -> A.Array A.S A.Ix1 Cl3_R
+compactRMassiv cl3 = A.makeVectorR A.S A.Seq (A.Sz1 (1024*1204)) (const $ toCl3_R cl3)
+
+mallocR :: Cl3 -> IO (Ptr Cl3_R)
+mallocR cl3 = do
+  pa <- mallocBytes ((1024*1204) * sizeOf (undefined :: Cl3_R))
+  let fill idx | idx < (1024*1204) = do
+        poke (pa `plusPtr` (idx * sizeOf (undefined :: Cl3_R))) cl3
+        fill (idx+1)
+      fill _ | otherwise = return ()
+  fill 0
+  peekElemOff pa (1024*1204)
+
+unitR :: Cl3
+unitR = R 1
diff --git a/cl3.cabal b/cl3.cabal
--- a/cl3.cabal
+++ b/cl3.cabal
@@ -10,7 +10,7 @@
 -- PVP summary:      +-+------- breaking API changes
 --                   | | +----- non-breaking API additions
 --                   | | | +--- code changes with no API change
-version:             2.0.0.0
+version:             3.0.0.0
 
 -- A short (one-line) description of the package.
 synopsis:            Clifford Algebra of three dimensional space.
@@ -37,7 +37,7 @@
 maintainer:          Nathan Waivio <nathan.waivio@gmail.com>
 
 -- A copyright notice.
-copyright:           Copyright (C) 2017-2020 Nathan Waivio
+copyright:           Copyright (C) 2017-2022 Nathan Waivio
 
 category:            Math, Algebra
 
@@ -50,7 +50,12 @@
                      GHC == 8.4.2,
                      GHC == 8.4.4,
                      GHC == 8.6.5,
-                     GHC == 8.8.3
+                     GHC == 8.8.3,
+                     GHC == 8.10.2,
+                     GHC == 8.10.4,
+                     GHC == 8.10.7,
+                     GHC == 9.0.2,
+                     GHC == 9.2.2
 
 -- Extra files to be distributed with the package, such as examples or a 
 -- README.
@@ -75,6 +80,16 @@
   manual:      True
   default:     False
 
+flag do-no-storable
+  description: Build without Storable Class support
+  manual:      True
+  default:     False
+
+flag do-liquid
+  description: Build with Liquid Haskell checking
+  manual:      True
+  default:     False
+
 library
   -- Modules exported by the library.
   exposed-modules:
@@ -84,12 +99,21 @@
   -- Compiler options
   ghc-options: -Wall -O2
   
+  if flag(do-liquid)
+    ghc-options: -fplugin=LiquidHaskell -fplugin-opt=LiquidHaskell:--fast -fplugin-opt=LiquidHaskell:--max-case-expand=0 -fplugin-opt=LiquidHaskell:--no-termination -fplugin-opt=LiquidHaskell:--short-names
+  
   if flag(do-no-derived-instances)
     cpp-options: -DO_NO_DERIVED
   
   if flag(do-no-random)
     cpp-options: -DO_NO_RANDOM
   
+  if flag(do-no-storable)
+    cpp-options: -DO_NO_STORABLE
+  
+  if flag(do-liquid)
+    cpp-options: -DO_LIQUID -DO_NO_RANDOM -DO_NO_DERIVED -DO_NO_STORABLE
+  
   -- LANGUAGE extensions used by modules in this package.
   other-extensions:
     Safe,
@@ -102,14 +126,22 @@
     MultiWayIf
   
   -- Other library packages from which modules are imported.
-  build-depends:       
-    base >=4.7 && <5,
+  build-depends:
     deepseq >=1.1 && <2
   
   if !flag(do-no-random)
     build-depends:
       random >=1.0 && <2
   
+  if !flag(do-liquid)
+    build-depends:
+      base >=4.7 && <5
+  
+  if flag(do-liquid)
+    build-depends:
+      liquid-base,
+      liquidhaskell >= 0.8.10
+  
   -- Directories containing source files.
   hs-source-dirs:      src
   
@@ -129,7 +161,7 @@
     random >=1.0 && <2
   default-language: Haskell2010
 
--- Criterion based benchmark
+-- Criterion based benchmark -threaded -rtsopts "-with-rtsopts=-N"
 benchmark bench-cl3-nbody
   type: exitcode-stdio-1.0
   hs-source-dirs: benchmarks
@@ -139,4 +171,31 @@
     cl3,
     base >=4.7 && <5,
     criterion >=1.1 && <2
+  default-language: Haskell2010
+
+-- Weigh based benchmark for Vector and Massiv
+benchmark bench-cl3-weigh
+  type: exitcode-stdio-1.0
+  hs-source-dirs: benchmarks
+  main-is: WeighStorableCl3.hs
+  ghc-options: -Wall -O2
+  build-depends:
+    cl3,
+    base >=4.7 && <5,
+    vector,
+    massiv >= 1.0,
+    weigh
+  default-language: Haskell2010
+
+-- Cl3/Massiv N-body benchmark -threaded -rtsopts "-with-rtsopts=-N"
+benchmark bench-cl3-massiv-nbody
+  type: exitcode-stdio-1.0
+  hs-source-dirs: benchmarks
+  main-is: NbodyMassiv.hs
+  ghc-options: -Wall -O2
+  build-depends:
+    cl3,
+    base >=4.7 && <5,
+    time,
+    massiv >= 1.0
   default-language: Haskell2010
diff --git a/src/Algebra/Geometric/Cl3.hs b/src/Algebra/Geometric/Cl3.hs
--- a/src/Algebra/Geometric/Cl3.hs
+++ b/src/Algebra/Geometric/Cl3.hs
@@ -1,4 +1,3 @@
-{-# LANGUAGE Safe #-}
 {-# LANGUAGE GADTSyntax #-}
 {-# LANGUAGE ViewPatterns #-}
 {-# LANGUAGE DeriveDataTypeable #-}
@@ -8,15 +7,23 @@
 {-# LANGUAGE CPP #-}
 {-# LANGUAGE MultiWayIf #-}
 
+
+#ifndef O_LIQUID
+-- Turn off Safe Haskell language extension due to liquid-base re-exports
+{-# LANGUAGE Safe #-}
+#endif
+
 #if __GLASGOW_HASKELL__ == 810
 -- Work around to fix GHC Issue #15304, issue popped up again in GHC 8.10, it should be fixed in GHC 8.12
 -- This code is meant to reproduce MR 2608 for GHC 8.10
 {-# OPTIONS_GHC -funfolding-keeness-factor=1 -funfolding-use-threshold=80 #-}
 #endif
 
+
+
 --------------------------------------------------------------------------------------------
 -- |
--- Copyright   :  (C) 2017-2020 Nathan Waivio
+-- Copyright   :  (C) 2017-2022 Nathan Waivio
 -- License     :  BSD3
 -- Maintainer  :  Nathan Waivio <nathan.waivio@gmail.com>
 -- Stability   :  Stable
@@ -43,6 +50,20 @@
  showOctave,
  -- * Eliminate grades that are less than 'tol' to use a simpler Constructor
  reduce, tol,
+#ifndef O_NO_STORABLE
+ -- * Compact Storable types for the Cl3 Constructors with smart constructors
+ Cl3_R, toCl3_R, fromCl3_R,
+ Cl3_V3, toCl3_V3, fromCl3_V3,
+ Cl3_BV, toCl3_BV, fromCl3_BV,
+ Cl3_I, toCl3_I, fromCl3_I,
+ Cl3_PV, toCl3_PV, fromCl3_PV,
+ Cl3_H, toCl3_H, fromCl3_H,
+ Cl3_C, toCl3_C, fromCl3_C,
+ Cl3_BPV, toCl3_BPV, fromCl3_BPV,
+ Cl3_ODD, toCl3_ODD, fromCl3_ODD,
+ Cl3_TPV, toCl3_TPV, fromCl3_TPV,
+ Cl3_APS, toCl3_APS, fromCl3_APS,
+#endif
 #ifndef O_NO_RANDOM
  -- * Random Instances
  randR, rangeR,
@@ -64,19 +85,22 @@
  -- * Helpful Functions
  eigvals, hasNilpotent,
  spectraldcmp, project,
- mIx, timesI
+ mIx, timesI,
+ abssignum
 ) where
 
 #ifndef O_NO_DERIVED
 import Data.Data (Typeable, Data)
 import GHC.Generics (Generic)
+import Text.Read (Read,readPrec)
 #endif
 
-
 import Control.DeepSeq (NFData,rnf)
+
+#ifndef O_NO_STORABLE
 import Foreign.Storable (Storable, sizeOf, alignment, peek, poke)
 import Foreign.Ptr (Ptr, plusPtr, castPtr)
-
+#endif
 
 #ifndef O_NO_RANDOM
 import System.Random (RandomGen, Random, randomR, random)
@@ -136,6 +160,7 @@
 #endif
 
 
+-- | Cl3 can be reduced to a normal form.
 instance NFData Cl3 where
   rnf !_ = ()
 
@@ -955,90 +980,20 @@
                                                                                 (a0*b123 + a123*b0 + a1*b23 + a23*b1 + a2*b31 + a31*b2 + a3*b12 + a12*b3)
 
 
+
+
   -- |'abs' is the spectral norm aka the spectral radius
   -- it is the largest singular value. This function may need to be fiddled with
   -- to make the math a bit safer wrt overflows.  This makes use of the largest
   -- singular value, if the littlest singular value is zero then the element is not
   -- invertable, we can see here that R, C, V3, BV, and H are all invertable, and
   -- by implication R, C, and H are division algebras.
-  abs (R a0) = R (abs a0) -- absolute value of a real number
-  abs (V3 a1 a2 a3) = R (sqrt (a1^2 + a2^2 + a3^2)) -- magnitude of a vector
-  abs (BV a23 a31 a12) = R (sqrt (a23^2 + a31^2 + a12^2)) -- magnitude of a bivector
-  abs (I a123) = R (abs a123) -- magnitude of a Imaginary number
-  abs (PV a0 a1 a2 a3) = R (reimMag a0 a1 a2 a3)
-  abs (TPV a23 a31 a12 a123) = R (reimMag a123 a23 a31 a12)
-  abs (H a0 a23 a31 a12) = R (sqrt (a0^2 + a23^2 + a31^2 + a12^2)) -- largest singular value
-  abs (C a0 a123) = R (sqrt (a0^2 + a123^2)) -- magnitude of a complex number
-  abs (BPV a1 a2 a3 a23 a31 a12) =
-    let x = sqrt ((a1*a31 - a2*a23)^2 + (a1*a12 - a3*a23)^2 + (a2*a12 - a3*a31)^2) -- core was duplicating this computation added let to hopefully reduce the duplication
-    in R (sqrt (a1^2 + a23^2 + a2^2 + a31^2 + a3^2 + a12^2 + x + x))
-  abs (ODD a1 a2 a3 a123) = R (sqrt (a1^2 + a2^2 + a3^2 + a123^2))
-  abs (APS a0 a1 a2 a3 a23 a31 a12 a123) =
-    let x = sqrt ((a0*a1 + a123*a23)^2 + (a0*a2 + a123*a31)^2 + (a0*a3 + a123*a12)^2 +
-                  (a2*a12 - a3*a31)^2 + (a3*a23 - a1*a12)^2 + (a1*a31 - a2*a23)^2) -- core was duplicating this computation added let to hopefully reduce the duplication
-    in R (sqrt (a0^2 + a1^2 + a2^2 + a3^2 + a23^2 + a31^2 + a12^2 + a123^2 + x + x))
+  abs cl3 = fst $ abssignum cl3
 
 
   -- |'signum' satisfies the Law "abs x * signum x == x"
   -- kind of cool: signum of a vector is it's unit vector.
-  signum (R a0) = R (signum a0)
-  signum (V3 a1 a2 a3) =
-    let mag = sqrt (a1^2 + a2^2 + a3^2)
-        invMag = recip mag
-    in if mag == 0
-       then R 0
-       else V3 (invMag * a1) (invMag * a2) (invMag * a3)
-  signum (BV a23 a31 a12) =
-    let mag = sqrt (a23^2 + a31^2 + a12^2)
-        invMag = recip mag
-    in if mag == 0
-       then R 0
-       else BV (invMag * a23) (invMag * a31) (invMag * a12)
-  signum (I a123) = I (signum a123)
-  signum (PV a0 a1 a2 a3) =
-    let mag = reimMag a0 a1 a2 a3
-        invMag = recip mag
-    in if mag == 0
-       then R 0
-       else PV (invMag * a0) (invMag * a1) (invMag * a2) (invMag * a3)
-  signum (H a0 a23 a31 a12) =
-    let mag = sqrt (a0^2 + a23^2 + a31^2 + a12^2)
-        invMag = recip mag
-    in if mag == 0
-       then R 0
-       else H (invMag * a0) (invMag * a23) (invMag * a31) (invMag * a12)
-  signum (C a0 a123) =
-    let mag = sqrt (a0^2 + a123^2)
-        invMag = recip mag
-    in if mag == 0
-       then R 0
-       else C (invMag * a0) (invMag * a123)
-  signum (BPV a1 a2 a3 a23 a31 a12) =
-    let x = sqrt ((a1*a31 - a2*a23)^2 + (a1*a12 - a3*a23)^2 + (a2*a12 - a3*a31)^2)
-        mag = sqrt (a1^2 + a23^2 + a2^2 + a31^2 + a3^2 + a12^2 + x + x)
-        invMag = recip mag
-    in if mag == 0
-       then R 0
-       else BPV (invMag * a1) (invMag * a2) (invMag * a3) (invMag * a23) (invMag * a31) (invMag * a12)
-  signum (ODD a1 a2 a3 a123) =
-    let mag = sqrt (a1^2 + a2^2 + a3^2 + a123^2)
-        invMag = recip mag
-    in if mag == 0
-       then R 0
-       else ODD (invMag * a1) (invMag * a2) (invMag * a3) (invMag * a123)
-  signum (TPV a23 a31 a12 a123) =
-    let mag = reimMag a123 a23 a31 a12
-        invMag = recip mag
-    in if mag == 0
-       then R 0
-       else TPV (invMag * a23) (invMag * a31) (invMag * a12) (invMag * a123)
-  signum (APS a0 a1 a2 a3 a23 a31 a12 a123) =
-    let x = sqrt ((a0*a1 + a123*a23)^2 + (a0*a2 + a123*a31)^2 + (a0*a3 + a123*a12)^2 + (a2*a12 - a3*a31)^2 + (a3*a23 - a1*a12)^2 + (a1*a31 - a2*a23)^2)
-        mag = sqrt (a0^2 + a1^2 + a2^2 + a3^2 + a23^2 + a31^2 + a12^2 + a123^2 + x + x)
-        invMag = recip mag
-    in if mag == 0
-       then R 0
-       else APS (invMag * a0) (invMag * a1) (invMag * a2) (invMag * a3) (invMag * a23) (invMag * a31) (invMag * a12) (invMag * a123)
+  signum cl3 = snd $ abssignum cl3
 
 
   -- |'fromInteger'
@@ -1072,7 +1027,7 @@
 reimMag v0 v1 v2 v3 =
   let sumsqs = v1^2 + v2^2 + v3^2
       x = abs v0 * sqrt sumsqs
-  in sqrt (v0^2 + sumsqs + x + x)
+  in sqrt (v0^2 + sumsqs + 2*x)
 
 -- |Cl(3,0) has a Fractional instance
 instance Fractional Cl3 where
@@ -1111,6 +1066,8 @@
     in C (expa0 * cos a123) (expa0 * sin a123)
   exp cliffor = spectraldcmp exp exp' cliffor
 
+
+
   --
   log (R a0)
     | a0 >= 0 = R (log a0)
@@ -1123,6 +1080,7 @@
   log (C a0 a123) = C (log (a0^2 + a123^2) / 2) (atan2 a123 a0)
   log cliffor = spectraldcmp log log' cliffor
 
+
   --
   sqrt (R a0)
     | a0 >= 0 = R (sqrt a0)
@@ -1158,7 +1116,7 @@
   tan (I a123)
     | a123 == 0 = R 0
     | otherwise = I (tanh a123)
-  tan (C a0 a123) =
+  tan (reduce -> C a0 a123) =
     let
       m = x2^2 + y2^2
       x1 = sinx*coshy
@@ -1172,6 +1130,8 @@
     in C ((x1*x2 + y1*y2)/m) ((x2*y1 - x1*y2)/m)
   tan cliffor = spectraldcmp tan tan' cliffor
 
+
+
   --
   asin (R a0)
       -- asin (R a0) = I (-1) * log (I 1 * R a0 + sqrt (1 - (R a0)^2))
@@ -1340,7 +1300,7 @@
   --
   tanh (R a0) = R (tanh a0)
   tanh (I a123) = I (tan a123)
-  tanh (C a0 a123) =
+  tanh (reduce -> C a0 a123) =
     let
       m = x2^2 + y2^2
       x1 = cosy*sinhx
@@ -1963,6 +1923,49 @@
 timesI (TPV a23 a31 a12 a123) = PV (negate a123) (negate a23) (negate a31) (negate a12)
 timesI (APS a0 a1 a2 a3 a23 a31 a12 a123) = APS (negate a123) (negate a23) (negate a31) (negate a12) a1 a2 a3 a0
 
+-- | 'abssignum' is a more effecient '\cl3 -> (abs cl3, signum cl3)'
+-- So 'abs' is always R and 'signum' is the same type of constructor as the input
+-- 'signum' is the element divided by its largest singular value 'abs'
+abssignum :: Cl3 -> (Cl3,Cl3)
+abssignum cl3 =
+  let (R m0) = absolute cl3
+  in if m0 == 0
+     then (R 0, R 0) -- (abs 0 == 0, signum 0 == 0)
+     else (R m0, cl3/(R m0))
+
+absolute :: Cl3 -> Cl3
+absolute (R a0) = R (abs a0)
+absolute (V3 a1 a2 a3) = let m = rss3 a1 a2 a3 in R m
+absolute (BV a23 a31 a12) = let m = rss3 a23 a31 a12 in R m
+absolute (I a123) = R (abs a123)
+absolute (PV a0 a1 a2 a3) = let m = reimMag a0 a1 a2 a3 in R m
+absolute (H a0 a23 a31 a12) = let m = rss4 a0 a23 a31 a12 in R m
+absolute (C a0 a123) = let m = rss2 a0 a123 in R m
+absolute (BPV a1 a2 a3 a23 a31 a12) = let mag0 = rss3 (a1*a31 - a2*a23) (a1*a12 - a3*a23) (a2*a12 - a3*a31)
+                                          m = sqrt $ a1^2 + a23^2 + a2^2 + a31^2 + a3^2 + a12^2 + 2*mag0 in R m
+absolute (ODD a1 a2 a3 a123) = let m = rss4 a1 a2 a3 a123 in R m
+absolute (TPV a23 a31 a12 a123) = let m = reimMag a123 a23 a31 a12 in R m
+absolute (APS a0 a1 a2 a3 a23 a31 a12 a123) = let mag0 = sqrt $ (a0*a1 + a123*a23)^2 + (a0*a2 + a123*a31)^2 + (a0*a3 + a123*a12)^2 + (a2*a12 - a3*a31)^2 + (a3*a23 - a1*a12)^2 + (a1*a31 - a2*a23)^2
+                                                  m = sqrt $ a0^2 + a1^2 + a2^2 + a3^2 + a23^2 + a31^2 + a12^2 + a123^2 + 2*mag0 in R m
+
+rss2 :: Double -> Double -> Double
+rss2 a0 a123 = sqrt $ a0^2 + a123^2
+
+rss3 :: Double -> Double -> Double -> Double
+rss3 x y z = sqrt $ x^2 + y^2 + z^2
+
+rss4 :: Double -> Double -> Double -> Double -> Double
+rss4 t x y z = sqrt $ t^2 + x^2 + y^2 + z^2
+
+
+
+#ifdef O_LIQUID
+tol :: Cl3
+tol = R 0
+
+tol' :: Double
+tol' = 0
+#else
 -- | 'tol' currently 128*eps
 tol :: Cl3
 {-# INLINE tol #-}
@@ -1971,7 +1974,7 @@
 tol' :: Double
 {-# INLINE tol' #-}
 tol' = 1.4210854715202004e-14
-
+#endif
 
 -- | 'bar' is a Clifford Conjugate, the vector grades are negated
 bar :: Cl3 -> Cl3
@@ -2208,6 +2211,7 @@
 atanh' = recip.(1-).(^2)  -- pole at +/-1
 
 
+#ifndef O_NO_STORABLE
 -------------------------------------------------------------------
 -- 
 -- Instance of Cl3 types with the "Foreign.Storable" library.
@@ -2219,6 +2223,9 @@
 
 -- | Cl3 instance of Storable uses the APS constructor as its standard interface.
 -- "peek" returns a cliffor constructed with APS. "poke" converts a cliffor to APS.
+-- For a more compact storing of constructors other than APS use the storable
+-- subtypes Cl3_R, Cl3_V3, Cl3_BV, Cl3_I, Cl3_PV, Cl3_H, Cl3_C, Cl3_BPV,
+-- Cl3_ODD, Cl3_TPV.
 instance Storable Cl3 where
   sizeOf _ = 8 * sizeOf (undefined :: Double)
   alignment _ = sizeOf (undefined :: Double)
@@ -2249,8 +2256,464 @@
   poke _ _ = error "Serious Issues with poke in Cl3.Storable"
 
 
+-- | 'Cl3_R' a compact storable data type for R.
+data Cl3_R where
+  Cl3_R :: !Double -> Cl3_R
 
+-- | 'toCl3_R' converts a Cl3 value constructed with R to its compact form.
+toCl3_R :: Cl3 -> Cl3_R
+toCl3_R (R a0) = Cl3_R a0
+toCl3_R err = error $ "Please don't try and cast something that's not R to Cl3_R, Got: " ++ show err
 
+-- | 'fromCl3_R' converts the compact Cl3_R type back to a Cl3 type.
+fromCl3_R :: Cl3_R -> Cl3
+fromCl3_R (Cl3_R a0) = R a0
+
+instance Show Cl3_R where
+  show = show.fromCl3_R
+
+#ifndef O_NO_DERIVED
+instance Read Cl3_R where
+  readPrec = toCl3_R <$> readPrec
+#endif
+
+instance Storable Cl3_R where
+  sizeOf _ = sizeOf (undefined :: Double)
+  alignment _ = sizeOf (undefined :: Double)
+  peek ptr = do
+    a0 <- peek offset
+    return $ Cl3_R a0
+      where
+        offset = (castPtr ptr :: Ptr Double)
+
+  poke ptr (Cl3_R a0) = do
+    poke offset a0
+      where
+        offset = (castPtr ptr :: Ptr Double)
+
+
+-- | 'Cl3_V3' a compact storable data type for V3.
+data Cl3_V3 where
+  Cl3_V3 :: !Double -> !Double -> !Double -> Cl3_V3
+
+-- | 'toCl3_V3' converts a Cl3 value constructed with V3 to its compact form.
+toCl3_V3 :: Cl3 -> Cl3_V3
+toCl3_V3 (V3 a1 a2 a3) = Cl3_V3 a1 a2 a3
+toCl3_V3 err = error $ "Please don't try and cast something that's not V3 to Cl3_V3, Got: " ++ show err
+
+-- | 'fromCl3_V3' converts the compact Cl3_V3 type back to a Cl3 type.
+fromCl3_V3 :: Cl3_V3 -> Cl3
+fromCl3_V3 (Cl3_V3 a1 a2 a3) = V3 a1 a2 a3
+
+instance Show Cl3_V3 where
+  show = show.fromCl3_V3
+
+#ifndef O_NO_DERIVED
+instance Read Cl3_V3 where
+  readPrec = toCl3_V3 <$> readPrec
+#endif
+
+instance Storable Cl3_V3 where
+  sizeOf _ = 3 * sizeOf (undefined :: Double)
+  alignment _ = sizeOf (undefined :: Double)
+  peek ptr = do
+    a1 <- peek (offset 0)
+    a2 <- peek (offset 1)
+    a3 <- peek (offset 2)
+    return $ Cl3_V3 a1 a2 a3
+      where
+        offset i = (castPtr ptr :: Ptr Double) `plusPtr` (i*8)
+
+  poke ptr (Cl3_V3 a1 a2 a3) = do
+    poke (offset 0) a1
+    poke (offset 1) a2
+    poke (offset 2) a3
+      where
+        offset i = (castPtr ptr :: Ptr Double) `plusPtr` (i*8)
+
+
+-- | 'Cl3_BV' a compact storable data type for BV.
+data Cl3_BV where
+  Cl3_BV :: !Double -> !Double -> !Double -> Cl3_BV
+
+-- | 'toCl3_BV' converts a Cl3 value constructed with BV to its compact form.
+toCl3_BV :: Cl3 -> Cl3_BV
+toCl3_BV (BV a23 a31 a12) = Cl3_BV a23 a31 a12
+toCl3_BV err = error $ "Please don't try and cast something that's not BV to Cl3_BV, Got: " ++ show err
+
+-- | 'fromCl3_BV' converts the compact Cl3_BV type back to a Cl3 type.
+fromCl3_BV :: Cl3_BV -> Cl3
+fromCl3_BV (Cl3_BV a23 a31 a12) = BV a23 a31 a12
+
+instance Show Cl3_BV where
+  show = show.fromCl3_BV
+
+#ifndef O_NO_DERIVED
+instance Read Cl3_BV where
+  readPrec = toCl3_BV <$> readPrec
+#endif
+
+instance Storable Cl3_BV where
+  sizeOf _ = 3 * sizeOf (undefined :: Double)
+  alignment _ = sizeOf (undefined :: Double)
+  peek ptr = do
+    a23 <- peek (offset 0)
+    a31 <- peek (offset 1)
+    a12 <- peek (offset 2)
+    return $ Cl3_BV a23 a31 a12
+      where
+        offset i = (castPtr ptr :: Ptr Double) `plusPtr` (i*8)
+
+  poke ptr (Cl3_BV a23 a31 a12) = do
+    poke (offset 0) a23
+    poke (offset 1) a31
+    poke (offset 2) a12
+      where
+        offset i = (castPtr ptr :: Ptr Double) `plusPtr` (i*8)
+
+
+-- | 'Cl3_I' a compact storable data type for I.
+data Cl3_I where
+  Cl3_I :: !Double -> Cl3_I
+
+-- | 'toCl3_I' converts a Cl3 value constructed with I to its compact form.
+toCl3_I :: Cl3 -> Cl3_I
+toCl3_I (I a123) = Cl3_I a123
+toCl3_I err = error $ "Please don't try and cast something that's not R to Cl3_R, Got: " ++ show err
+
+-- | 'fromCl3_I' converts the compact Cl3_I type back to a Cl3 type.
+fromCl3_I :: Cl3_I -> Cl3
+fromCl3_I (Cl3_I a123) = I a123
+
+instance Show Cl3_I where
+  show = show.fromCl3_I
+
+#ifndef O_NO_DERIVED
+instance Read Cl3_I where
+  readPrec = toCl3_I <$> readPrec
+#endif
+
+instance Storable Cl3_I where
+  sizeOf _ = sizeOf (undefined :: Double)
+  alignment _ = sizeOf (undefined :: Double)
+  peek ptr = do
+    a123 <- peek offset
+    return $ Cl3_I a123
+      where
+        offset = (castPtr ptr :: Ptr Double)
+
+  poke ptr (Cl3_I a123) = do
+    poke offset a123
+      where
+        offset = (castPtr ptr :: Ptr Double)
+
+
+-- | 'Cl3_PV' a compact storable data type for PV.
+data Cl3_PV where
+  Cl3_PV :: !Double -> !Double -> !Double -> !Double -> Cl3_PV
+
+-- | 'toCl3_PV' converts a Cl3 value constructed with PV to its compact form.
+toCl3_PV :: Cl3 -> Cl3_PV
+toCl3_PV (PV a0 a1 a2 a3) = Cl3_PV a0 a1 a2 a3
+toCl3_PV err = error $ "Please don't try and cast something that's not PV to Cl3_PV, Got: " ++ show err
+
+-- | 'fromCl3_PV' converts the compact Cl3_PV type back to a Cl3 type.
+fromCl3_PV :: Cl3_PV -> Cl3
+fromCl3_PV (Cl3_PV a0 a1 a2 a3) = PV a0 a1 a2 a3
+
+instance Show Cl3_PV where
+  show = show.fromCl3_PV
+
+#ifndef O_NO_DERIVED
+instance Read Cl3_PV where
+  readPrec = toCl3_PV <$> readPrec
+#endif
+
+instance Storable Cl3_PV where
+  sizeOf _ = 4 * sizeOf (undefined :: Double)
+  alignment _ = sizeOf (undefined :: Double)
+  peek ptr = do
+    a0 <- peek (offset 0)
+    a1 <- peek (offset 1)
+    a2 <- peek (offset 2)
+    a3 <- peek (offset 4)
+    return $ Cl3_PV a0 a1 a2 a3
+      where
+        offset i = (castPtr ptr :: Ptr Double) `plusPtr` (i*8)
+
+  poke ptr (Cl3_PV a0 a1 a2 a3) = do
+    poke (offset 0) a0
+    poke (offset 1) a1
+    poke (offset 2) a2
+    poke (offset 3) a3
+      where
+        offset i = (castPtr ptr :: Ptr Double) `plusPtr` (i*8)
+
+
+-- | 'Cl3_H' a compact storable data type for H.
+data Cl3_H where
+  Cl3_H :: !Double -> !Double -> !Double -> !Double -> Cl3_H
+
+-- | 'toCl3_H' converts a Cl3 value constructed with H to its compact form.
+toCl3_H :: Cl3 -> Cl3_H
+toCl3_H (H a0 a23 a31 a12) = Cl3_H a0 a23 a31 a12
+toCl3_H err = error $ "Please don't try and cast something that's not H to Cl3_H, Got: " ++ show err
+
+-- | 'fromCl3_H' converts the compact Cl3_H type back to a Cl3 type.
+fromCl3_H :: Cl3_H -> Cl3
+fromCl3_H (Cl3_H a0 a23 a31 a12) = H a0 a23 a31 a12
+
+instance Show Cl3_H where
+  show = show.fromCl3_H
+
+#ifndef O_NO_DERIVED
+instance Read Cl3_H where
+  readPrec = toCl3_H <$> readPrec
+#endif
+
+instance Storable Cl3_H where
+  sizeOf _ = 4 * sizeOf (undefined :: Double)
+  alignment _ = sizeOf (undefined :: Double)
+  peek ptr = do
+    a0 <- peek (offset 0)
+    a23 <- peek (offset 1)
+    a31 <- peek (offset 2)
+    a12 <- peek (offset 3)
+    return $ Cl3_H a0 a23 a31 a12
+      where
+        offset i = (castPtr ptr :: Ptr Double) `plusPtr` (i*8)
+
+  poke ptr (Cl3_H a0 a23 a31 a12) = do
+    poke (offset 0) a0
+    poke (offset 1) a23
+    poke (offset 2) a31
+    poke (offset 3) a12
+      where
+        offset i = (castPtr ptr :: Ptr Double) `plusPtr` (i*8)
+
+
+-- | 'Cl3_C' a compact storable data type for C.
+data Cl3_C where
+  Cl3_C :: !Double -> !Double -> Cl3_C
+
+-- | 'toCl3_C' converts a Cl3 value constructed with C to its compact form.
+toCl3_C :: Cl3 -> Cl3_C
+toCl3_C (C a0 a123) = Cl3_C a0 a123
+toCl3_C err = error $ "Please don't try and cast something that's not C to Cl3_C, Got: " ++ show err
+
+-- | 'fromCl3_C' converts the compact Cl3_C type back to a Cl3 type.
+fromCl3_C :: Cl3_C -> Cl3
+fromCl3_C (Cl3_C a0 a123) = C a0 a123
+
+instance Show Cl3_C where
+  show = show.fromCl3_C
+
+#ifndef O_NO_DERIVED
+instance Read Cl3_C where
+  readPrec = toCl3_C <$> readPrec
+#endif
+
+instance Storable Cl3_C where
+  sizeOf _ = 2 * sizeOf (undefined :: Double)
+  alignment _ = sizeOf (undefined :: Double)
+  peek ptr = do
+    a0 <- peek (offset 0)
+    a123 <- peek (offset 1)
+    return $ Cl3_C a0 a123
+      where
+        offset i = (castPtr ptr :: Ptr Double) `plusPtr` (i*8)
+
+  poke ptr (Cl3_C a0 a123) = do
+    poke (offset 0) a0
+    poke (offset 1) a123
+      where
+        offset i = (castPtr ptr :: Ptr Double) `plusPtr` (i*8)
+
+
+-- | 'Cl3_BPV' a compact storable data type for BPV.
+data Cl3_BPV where
+  Cl3_BPV :: !Double -> !Double -> !Double -> !Double -> !Double -> !Double -> Cl3_BPV
+
+-- | 'toCl3_BPV' converts a Cl3 value constructed with BPV to its compact form.
+toCl3_BPV :: Cl3 -> Cl3_BPV
+toCl3_BPV (BPV a1 a2 a3 a23 a31 a12) = Cl3_BPV a1 a2 a3 a23 a31 a12
+toCl3_BPV err = error $ "Please don't try and cast something that's not BPV to Cl3_BPV, Got: " ++ show err
+
+-- | 'fromCl3_BPV' converts the compact Cl3_BPV type back to a Cl3 type.
+fromCl3_BPV :: Cl3_BPV -> Cl3
+fromCl3_BPV (Cl3_BPV a1 a2 a3 a23 a31 a12) = BPV a1 a2 a3 a23 a31 a12
+
+instance Show Cl3_BPV where
+  show = show.fromCl3_BPV
+
+#ifndef O_NO_DERIVED
+instance Read Cl3_BPV where
+  readPrec = toCl3_BPV <$> readPrec
+#endif
+
+instance Storable Cl3_BPV where
+  sizeOf _ = 6 * sizeOf (undefined :: Double)
+  alignment _ = sizeOf (undefined :: Double)
+  peek ptr = do
+    a1 <- peek (offset 0)
+    a2 <- peek (offset 1)
+    a3 <- peek (offset 2)
+    a23 <- peek (offset 3)
+    a31 <- peek (offset 4)
+    a12 <- peek (offset 5)
+    return $ Cl3_BPV a1 a2 a3 a23 a31 a12
+      where
+        offset i = (castPtr ptr :: Ptr Double) `plusPtr` (i*8)
+
+  poke ptr (Cl3_BPV a1 a2 a3 a23 a31 a12) = do
+    poke (offset 0) a1
+    poke (offset 1) a2
+    poke (offset 2) a3
+    poke (offset 3) a23
+    poke (offset 4) a31
+    poke (offset 5) a12
+      where
+        offset i = (castPtr ptr :: Ptr Double) `plusPtr` (i*8)
+
+
+-- | 'Cl3_ODD' a compact storable data type for ODD.
+data Cl3_ODD where
+  Cl3_ODD :: !Double -> !Double -> !Double -> !Double -> Cl3_ODD
+
+-- | 'toCl3_ODD' converts a Cl3 value constructed with ODD to its compact form.
+toCl3_ODD :: Cl3 -> Cl3_ODD
+toCl3_ODD (ODD a1 a2 a3 a123) = Cl3_ODD a1 a2 a3 a123
+toCl3_ODD err = error $ "Please don't try and cast something that's not ODD to Cl3_ODD, Got: " ++ show err
+
+-- | 'fromCl3_ODD' converts the compact Cl3_ODD type back to a Cl3 type.
+fromCl3_ODD :: Cl3_ODD -> Cl3
+fromCl3_ODD (Cl3_ODD a1 a2 a3 a123) = ODD a1 a2 a3 a123
+
+instance Show Cl3_ODD where
+  show = show.fromCl3_ODD
+
+#ifndef O_NO_DERIVED
+instance Read Cl3_ODD where
+  readPrec = toCl3_ODD <$> readPrec
+#endif
+
+instance Storable Cl3_ODD where
+  sizeOf _ = 4 * sizeOf (undefined :: Double)
+  alignment _ = sizeOf (undefined :: Double)
+  peek ptr = do
+    a1 <- peek (offset 0)
+    a2 <- peek (offset 1)
+    a3 <- peek (offset 2)
+    a123 <- peek (offset 3)
+    return $ Cl3_ODD a1 a2 a3 a123
+      where
+        offset i = (castPtr ptr :: Ptr Double) `plusPtr` (i*8)
+
+  poke ptr (Cl3_ODD a1 a2 a3 a123) = do
+    poke (offset 0) a1
+    poke (offset 1) a2
+    poke (offset 2) a3
+    poke (offset 3) a123
+      where
+        offset i = (castPtr ptr :: Ptr Double) `plusPtr` (i*8)
+
+
+-- | 'Cl3_TPV' a compact storable data type for TPV.
+data Cl3_TPV where
+  Cl3_TPV :: !Double -> !Double -> !Double -> !Double -> Cl3_TPV
+
+-- | 'toCl3_TPV' converts a Cl3 value constructed with TPV to its compact form.
+toCl3_TPV :: Cl3 -> Cl3_TPV
+toCl3_TPV (TPV a23 a31 a12 a123) = Cl3_TPV a23 a31 a12 a123
+toCl3_TPV err = error $ "Please don't try and cast something that's not TPV to Cl3_TPV, Got: " ++ show err
+
+-- | 'fromCl3_TPV' converts the compact Cl3_TPV type back to a Cl3 type.
+fromCl3_TPV :: Cl3_TPV -> Cl3
+fromCl3_TPV (Cl3_TPV a23 a31 a12 a123) = TPV a23 a31 a12 a123
+
+instance Show Cl3_TPV where
+  show = show.fromCl3_TPV
+
+#ifndef O_NO_DERIVED
+instance Read Cl3_TPV where
+  readPrec = toCl3_TPV <$> readPrec
+#endif
+
+instance Storable Cl3_TPV where
+  sizeOf _ = 4 * sizeOf (undefined :: Double)
+  alignment _ = sizeOf (undefined :: Double)
+  peek ptr = do
+    a23 <- peek (offset 0)
+    a31 <- peek (offset 1)
+    a12 <- peek (offset 2)
+    a123 <- peek (offset 3)
+    return $ Cl3_TPV a23 a31 a12 a123
+      where
+        offset i = (castPtr ptr :: Ptr Double) `plusPtr` (i*8)
+
+  poke ptr (Cl3_TPV a23 a31 a12 a123) = do
+    poke (offset 0) a23
+    poke (offset 1) a31
+    poke (offset 2) a12
+    poke (offset 3) a123
+      where
+        offset i = (castPtr ptr :: Ptr Double) `plusPtr` (i*8)
+
+
+-- | 'Cl3_APS' a compact storable data type for APS.
+data Cl3_APS where
+  Cl3_APS :: !Double -> !Double -> !Double -> !Double -> !Double -> !Double -> !Double -> !Double -> Cl3_APS
+
+-- | 'toCl3_APS' converts a Cl3 value constructed with APS to its compact form.
+toCl3_APS :: Cl3 -> Cl3_APS
+toCl3_APS (APS a0 a1 a2 a3 a23 a31 a12 a123) = Cl3_APS a0 a1 a2 a3 a23 a31 a12 a123
+toCl3_APS err = error $ "Please don't try and cast something that's not APS to Cl3_APS, Got: " ++ show err
+
+-- | 'fromCl3_APS' converts the compact Cl3_APS type back to a Cl3 type.
+fromCl3_APS :: Cl3_APS -> Cl3
+fromCl3_APS (Cl3_APS a0 a1 a2 a3 a23 a31 a12 a123) = APS a0 a1 a2 a3 a23 a31 a12 a123
+
+instance Show Cl3_APS where
+  show = show.fromCl3_APS
+
+#ifndef O_NO_DERIVED
+instance Read Cl3_APS where
+  readPrec = toCl3_APS <$> readPrec
+#endif
+
+instance Storable Cl3_APS where
+  sizeOf _ = 8 * sizeOf (undefined :: Double)
+  alignment _ = sizeOf (undefined :: Double)
+  peek ptr = do
+    a0 <- peek (offset 0)
+    a1 <- peek (offset 1)
+    a2 <- peek (offset 2)
+    a3 <- peek (offset 3)
+    a23 <- peek (offset 4)
+    a31 <- peek (offset 5)
+    a12 <- peek (offset 6)
+    a123 <- peek (offset 7)
+    return $ Cl3_APS a0 a1 a2 a3 a23 a31 a12 a123
+      where
+        offset i = (castPtr ptr :: Ptr Double) `plusPtr` (i*8)
+
+  poke ptr (Cl3_APS a0 a1 a2 a3 a23 a31 a12 a123) = do
+    poke (offset 0) a0
+    poke (offset 1) a1
+    poke (offset 2) a2
+    poke (offset 3) a3
+    poke (offset 4) a23
+    poke (offset 5) a31
+    poke (offset 6) a12
+    poke (offset 7) a123
+      where
+        offset i = (castPtr ptr :: Ptr Double) `plusPtr` (i*8)
+
+
+#endif
+
+
+
 #ifndef O_NO_RANDOM
 -------------------------------------------------------------------
 -- 
@@ -2476,8 +2939,8 @@
 randUnitV3 g =
   let (theta, g') = randomR (0,2*pi) g
       (u, g'') = randomR (-1,1) g'
-      simicircle = sqrt (1-u^2)
-  in (V3 (simicircle * cos theta) (simicircle * sin theta) u, g'')
+      semicircle = sqrt (1-u^2)
+  in (V3 (semicircle * cos theta) (semicircle * sin theta) u, g'')
 
 
 -- | 'randProjector' a projector with a random direction
diff --git a/src/Algebra/Geometric/Cl3/JonesCalculus.hs b/src/Algebra/Geometric/Cl3/JonesCalculus.hs
--- a/src/Algebra/Geometric/Cl3/JonesCalculus.hs
+++ b/src/Algebra/Geometric/Cl3/JonesCalculus.hs
@@ -1,7 +1,10 @@
-{-# LANGUAGE Safe #-}
 {-# LANGUAGE ViewPatterns #-}
 {-# OPTIONS_GHC -fno-warn-missing-signatures #-}
 {-# LANGUAGE CPP #-}
+
+#ifndef O_LIQUID
+{-# LANGUAGE Safe #-}
+#endif
 
 #if __GLASGOW_HASKELL__ == 810
 -- Work around to fix GHC Issue #15304, issue popped up again in GHC 8.10, it should be fixed in GHC 8.12
diff --git a/stack.yaml b/stack.yaml
--- a/stack.yaml
+++ b/stack.yaml
@@ -1,3 +1,19 @@
+# resolver: nightly-2022-03-23 # ghc-9.2.2
+# resolver: lts-19.0 # ghc-9.0.2
+# resolver: lts-18.28 # ghc-8.10.7
+resolver: lts-18.6 # ghc-8.10.4 # Currently the only version that seems to work with LiquidHaskell
 packages:
-- .
-resolver: nightly-2020-07-19
+  - .
+allow-newer: true
+extra-deps:
+  # For Massiv benchmarks:
+  - massiv-1.0.1.1
+  - scheduler-2.0.0
+  # For LiquidHaskell:
+  - text-format-0.3.2
+  - Diff-0.3.4
+  - optparse-applicative-0.16.1.0
+  - liquid-fixpoint-0.8.10.2
+  - liquidhaskell-0.8.10.2
+  - liquid-base-4.14.1.0
+  - liquid-ghc-prim-0.6.1