cl3 2.0.0.0 → 3.0.0.0

raw patch · 9 files changed

+1153/−102 lines, 9 filesdep +liquid-basedep +liquidhaskelldep +massivPVP ok

version bump matches the API change (PVP)

Dependencies added: liquid-base, liquidhaskell, massiv, time, vector, weigh

API changes (from Hackage documentation)

Files

• ChangeLog.md +24/−0
• README.md +5/−1
• benchmarks/NbodyGameCl3.hs +156/−7
• benchmarks/NbodyMassiv.hs +293/−0
• benchmarks/WeighStorableCl3.hs +40/−0
• cl3.cabal +65/−6
• src/Algebra/Geometric/Cl3.hs +548/−85
• src/Algebra/Geometric/Cl3/JonesCalculus.hs +4/−1
• stack.yaml +18/−2

@@ -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

@@ -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.

@@ -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 

@@ -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

@@ -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

@@ -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

@@ -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

@@ -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

@@ -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