packages feed

geomancy (empty) → 0.2.2.3

raw patch · 23 files changed

+2896/−0 lines, 23 filesdep +basedep +criteriondep +deepseqsetup-changed

Dependencies added: base, criterion, deepseq, geomancy, hedgehog, linear

Files

+ ChangeLog.md view
@@ -0,0 +1,76 @@+# Changelog for geomancy++## 0.2.2.3+++ Add IVec and UVec 32-bit integer vectors.++## 0.2.2.2+++ Fixed using Transform.apply with projection inverses.++## 0.2.2.1+++ Add Vec2.++ Add pattern synonyms to `Geomancy` re-exports.++ Add Fractional instances.++ Add conversions from tuples and lower-dimension vecs.++## 0.2.2.0+++ Add pattern synonym alternatives to `withVecN`.++ Add Mat4 converstion from `Linear.M44`.+- Hide `toList` and `toListTrans`.++ Add `toList2d` and publish its element-order wrappers.++## 0.2.1.0++Transform rewrite++* Change perspective FoV to radians.++ Add `infinitePerspective`.++ Add `Transform.inverse`.+- Hide `mat4`, `withMat4`.++## 0.2.0.0++Mat4 rewrite+++ Add `rowMajor`, `withRowMajor`, `toListRowMajor`.++ Add `colMajor`, `withColMajor`, `toListColMajor`.++ Add `Mat4.inverse`.+* Expose `matrixProduct`.++ Add `scalarMultiply`.++ Add `Mat4.(!*)` to use with `Vec4`.+* Extract transformations to `Geomancy.Transform` and use column notation.++ Add `Transform.(!.)` and `apply` to use with `Vec3`.+* Rename `scale` to `scale3`.++ Add uniform `scale`.++ Add `scaleXY` for flat meshes.+* Rename `mkTransformation` to `dirPos`.++ Add `rotateQ` via `dirPos` with empty translation.++ Add `Vec3.Packed` newtype without extra padding.++## 0.1.3.0++* Update tests++ Add `zipWith`+* Rename `elementwise` to `pointwise`+* Rename `colMajor` to `toList`+* Rename `rowMajor` to `toListTrans`++## 0.1.2.1+++ Add Mat4 multiplication test via `linear`++ Add `elementwise`, `colMajor`, `rowMajor`++## 0.1.2.0++* Move projections and views to Vulkan namespace.++## 0.1.1.2++* Fix bug in quaternion rotationBetween.++## 0.1.1.1++* Add lookAtUp and rotationBetween.
+ LICENSE view
@@ -0,0 +1,30 @@+Copyright Alexander Bondarenko (c) 2020++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions are met:++    * Redistributions of source code must retain the above copyright+      notice, this list of conditions and the following disclaimer.++    * Redistributions in binary form must reproduce the above+      copyright notice, this list of conditions and the following+      disclaimer in the documentation and/or other materials provided+      with the distribution.++    * Neither the name of Alexander Bondarenko nor the names of other+      contributors may be used to endorse or promote products derived+      from this software without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ README.md view
@@ -0,0 +1,62 @@+# Geomancy++Linear is nice, but slow. Those are naughty, but a bit faster.++* All data types are monomorphic, unpacked and specialized to `Float`.+* `Mat4` is a `ByteArray#`.+* `Mat4` multiplication with SSE.+* Matrix construction states their element order.+* Transforms don't require transposition for GLSL.++### The Numbers++Storing a list of 1000 transformations (e.g. rendering instance data):++```+benchmarking 4x4 poke/1000/geomancy+time                 11.76 μs   (11.66 μs .. 11.92 μs)+                     0.999 R²   (0.998 R² .. 1.000 R²)+mean                 11.75 μs   (11.69 μs .. 11.86 μs)+std dev              283.4 ns   (199.0 ns .. 399.0 ns)+variance introduced by outliers: 26% (moderately inflated)+```++If you're willing to adjust your shaders, it's only 2.4 times slower.++```+benchmarking 4x4 poke/1000/linear+time                 28.29 μs   (28.21 μs .. 28.38 μs)+                     1.000 R²   (1.000 R² .. 1.000 R²)+mean                 28.40 μs   (28.34 μs .. 28.50 μs)+std dev              267.4 ns   (145.5 ns .. 419.9 ns)+```++Keeping your shaders straight make the affair 6.1x slower.++```+benchmarking 4x4 poke/1000/linear/T+time                 73.70 μs   (73.06 μs .. 74.49 μs)+                     1.000 R²   (0.999 R² .. 1.000 R²)+mean                 72.77 μs   (72.50 μs .. 73.22 μs)+std dev              1.129 μs   (793.5 ns .. 1.580 μs)+```++Folding down a `gloss`-style scene graph is where it is all started:++```+benchmarking 4x4 multiply/1000/geomancy+time                 20.79 μs   (20.77 μs .. 20.83 μs)+                     1.000 R²   (1.000 R² .. 1.000 R²)+mean                 20.80 μs   (20.78 μs .. 20.83 μs)+std dev              76.71 ns   (60.01 ns .. 99.06 ns)++benchmarking 4x4 multiply/1000/linear+time                 173.9 μs   (173.6 μs .. 174.4 μs)+                     1.000 R²   (1.000 R² .. 1.000 R²)+mean                 173.5 μs   (173.2 μs .. 174.4 μs)+std dev              1.733 μs   (727.8 ns .. 3.422 μs)+```++Add that time to the poking that'll follow.++Sure, it is in the lower microseconds range, but this budget can be used elsewhere.
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ bench/Main.hs view
@@ -0,0 +1,75 @@+{-# LANGUAGE BlockArguments #-}+{-# LANGUAGE TypeApplications #-}++module Main where++import Criterion.Main++import Data.List (foldl1')+import qualified Foreign++import qualified Geomancy+import qualified Geomancy.Mat4++import qualified Linear+import qualified Linear.Matrix++main :: IO ()+main = defaultMain+  [ bgroup "4x4 transpose"+      mat4transpose++  , bgroup "4x4 multiply" $ flip map [2, 3, 10, 100, 1000] \size ->+      env (clones size mempty Linear.identity) $+        mat4multiply size++  , bgroup "4x4 poke" $ flip map [10, 100, 1000, 10000] \size ->+      env (clones size mempty Linear.identity) $+        mat4poke size+  ]++mat4transpose :: [Benchmark]+mat4transpose =+  [ bench "geomancy" $ whnf Geomancy.Mat4.transpose mempty+  , bench "linear"   $ whnf Linear.Matrix.transpose (Linear.identity :: Linear.M44 Float)+  ]++mat4multiply+  :: Int+  -> ([Geomancy.Mat4], [Linear.M44 Float])+  -> Benchmark+mat4multiply size ~(mat4s, m44s) = bgroup (show size)+  [ bench "geomancy" $ whnf multGeomancy mat4s+  , bench "linear"   $ whnf multLinear m44s+  ]+  where+    multGeomancy = foldl1' (<>)+    multLinear = foldl1' (Linear.!*!)++clones :: Applicative f => Int -> a -> b -> f ([a], [b])+clones n geo lin = (,) <$> genMat4s <*> genM44s+  where+    genMat4s = pure $! replicate n geo+    genM44s = pure $! replicate n lin++mat4poke+  :: Int+  -> ([Geomancy.Mat4], [Linear.M44 Float])+  -> Benchmark+mat4poke size ~(mat4s, m44s) = bgroup (show size)+  [ bench "geomancy" $ nfIO pokeGeomancy --  mat4s+  , bench "linear"   $ nfIO pokeLinear --  m44s+  , bench "linear/T" $ nfIO pokeLinearTranspose+  ]+  where+    pokeGeomancy =+      Foreign.withArray mat4s \_ptr ->+        pure ()++    pokeLinear =+      Foreign.withArray m44s \_ptr ->+        pure ()++    pokeLinearTranspose =+      Foreign.withArray (map Linear.transpose m44s) \_ptr ->+        pure ()
+ cbits/mat4.c view
@@ -0,0 +1,22 @@+#include <xmmintrin.h>++void M4x4_SSE(float *A, float *B, float *C) {+    __m128 row1 = _mm_load_ps(&B[0]);+    __m128 row2 = _mm_load_ps(&B[4]);+    __m128 row3 = _mm_load_ps(&B[8]);+    __m128 row4 = _mm_load_ps(&B[12]);+    for(int i=0; i<4; i++) {+        __m128 brod1 = _mm_set1_ps(A[4*i + 0]);+        __m128 brod2 = _mm_set1_ps(A[4*i + 1]);+        __m128 brod3 = _mm_set1_ps(A[4*i + 2]);+        __m128 brod4 = _mm_set1_ps(A[4*i + 3]);+        __m128 row = _mm_add_ps(+                    _mm_add_ps(+                        _mm_mul_ps(brod1, row1),+                        _mm_mul_ps(brod2, row2)),+                    _mm_add_ps(+                        _mm_mul_ps(brod3, row3),+                        _mm_mul_ps(brod4, row4)));+        _mm_store_ps(&C[4*i], row);+    }+}
+ geomancy.cabal view
@@ -0,0 +1,88 @@+cabal-version: 1.12++-- This file has been generated from package.yaml by hpack version 0.33.0.+--+-- see: https://github.com/sol/hpack+--+-- hash: 0bbbcdde42a68176d0c54cb31c5798701541c1ce73877ad066e6aec522450641++name:           geomancy+version:        0.2.2.3+synopsis:       Geometry and matrix manipulation+description:    Sometimes it is unavoidable you have to do stuff on CPU.+                Let's at least do it faster.+category:       Graphics+author:         Alexander Bondarenko+maintainer:     aenor.realm@gmail.com+copyright:      2021 Alexander Bondarenko+license:        BSD3+license-file:   LICENSE+build-type:     Simple+extra-source-files:+    README.md+    ChangeLog.md++source-repository head+  type: git+  location: https://gitlab.com/dpwiz/geomancy++library+  exposed-modules:+      Geomancy+      Geomancy.IVec2+      Geomancy.IVec3+      Geomancy.IVec4+      Geomancy.Mat4+      Geomancy.Quaternion+      Geomancy.Transform+      Geomancy.UVec2+      Geomancy.UVec3+      Geomancy.UVec4+      Geomancy.Vec2+      Geomancy.Vec3+      Geomancy.Vec4+      Geomancy.Vulkan.Projection+      Geomancy.Vulkan.View+  other-modules:+      Paths_geomancy+  hs-source-dirs:+      src+  ghc-options: -Wall+  c-sources:+      cbits/mat4.c+  build-depends:+      base >=4.7 && <5+    , deepseq+  default-language: Haskell2010++test-suite geomancy-test+  type: exitcode-stdio-1.0+  main-is: Spec.hs+  other-modules:+      Paths_geomancy+  hs-source-dirs:+      test+  ghc-options: -Wall -threaded -rtsopts -with-rtsopts=-N -with-rtsopts=-qn1 -with-rtsopts=-A128m+  build-depends:+      base >=4.7 && <5+    , deepseq+    , geomancy+    , hedgehog+    , linear+  default-language: Haskell2010++benchmark geomancy-bench+  type: exitcode-stdio-1.0+  main-is: Main.hs+  other-modules:+      Paths_geomancy+  hs-source-dirs:+      bench+  ghc-options: -Wall -threaded -rtsopts -with-rtsopts=-N -with-rtsopts=-qn1 -with-rtsopts=-A16m+  build-depends:+      base >=4.7 && <5+    , criterion+    , deepseq+    , geomancy+    , linear+  default-language: Haskell2010
+ src/Geomancy.hs view
@@ -0,0 +1,85 @@+{-# LANGUAGE PatternSynonyms #-}++module Geomancy+  (+  -- * Vectors++  -- ** Single-precision / float32s++    Vec2+  , vec2+  , withVec2+  , pattern WithVec2++  , Vec3+  , vec3+  , withVec3+  , pattern WithVec3++  , Vec4+  , vec4+  , withVec4+  , pattern WithVec4++  -- ** Signed / int32s++  , IVec2+  , ivec2+  , withIVec2+  , pattern WithIVec2++  , IVec3+  , ivec3+  , withIVec3+  , pattern WithIVec3++  , IVec4+  , ivec4+  , withIVec4+  , pattern WithIVec4++  -- ** Unsigned / word32s++  , UVec2+  , uvec2+  , withUVec2+  , pattern WithUVec2++  , UVec3+  , uvec3+  , withUVec3+  , pattern WithUVec3++  , UVec4+  , uvec4+  , withUVec4+  , pattern WithUVec4++  -- * Matrices++  , Mat4+  , Transform(..)++  -- * Other beasts++  , Quaternion+  , quaternion+  , withQuaternion+  ) where++import Geomancy.Vec2 (Vec2, vec2, withVec2, pattern WithVec2)+import Geomancy.Vec3 (Vec3, vec3, withVec3, pattern WithVec3)+import Geomancy.Vec4 (Vec4, vec4, withVec4, pattern WithVec4)++import Geomancy.IVec2 (IVec2, ivec2, withIVec2, pattern WithIVec2)+import Geomancy.IVec3 (IVec3, ivec3, withIVec3, pattern WithIVec3)+import Geomancy.IVec4 (IVec4, ivec4, withIVec4, pattern WithIVec4)++import Geomancy.UVec2 (UVec2, uvec2, withUVec2, pattern WithUVec2)+import Geomancy.UVec3 (UVec3, uvec3, withUVec3, pattern WithUVec3)+import Geomancy.UVec4 (UVec4, uvec4, withUVec4, pattern WithUVec4)++import Geomancy.Mat4 (Mat4)+import Geomancy.Transform (Transform(..))++import Geomancy.Quaternion (Quaternion, quaternion, withQuaternion)
+ src/Geomancy/IVec2.hs view
@@ -0,0 +1,95 @@+{-# LANGUAGE BlockArguments #-}+{-# LANGUAGE PatternSynonyms #-}+{-# LANGUAGE ViewPatterns #-}++-- | Specialized and inlined @V2 Int32@.++module Geomancy.IVec2+  ( IVec2+  , ivec2+  , withIVec2+  , pattern WithIVec2+  , fromTuple+  ) where++import Control.DeepSeq (NFData(rnf))+import Data.Int (Int32)+import Foreign (Storable(..))++data IVec2 = IVec2+  {-# UNPACK #-} !Int32+  {-# UNPACK #-} !Int32+  deriving (Eq, Ord, Show)++{-# INLINE ivec2 #-}+ivec2 :: Int32 -> Int32 -> IVec2+ivec2 = IVec2++{-# INLINE withIVec2 #-}+withIVec2+  :: IVec2+  -> (Int32 -> Int32 -> r)+  -> r+withIVec2 (IVec2 a b) f = f a b++pattern WithIVec2 :: Int32 -> Int32 -> IVec2+pattern WithIVec2 a b <- ((`withIVec2` (,)) -> (a, b))+{-# COMPLETE WithIVec2 #-}++{-# INLINE fromTuple #-}+fromTuple :: (Int32, Int32) -> IVec2+fromTuple (x, y) = ivec2 x y++instance NFData IVec2 where+  rnf IVec2{} = ()++-- XXX: That's one nasty instance...+instance Num IVec2 where+  {-# INLINE (+) #-}+  IVec2 l1 l2 + IVec2 r1 r2 =+    IVec2+      (l1 + r1)+      (l2 + r2)++  {-# INLINE (-) #-}+  IVec2 l1 l2 - IVec2 r1 r2 =+    IVec2+      (l1 - r1)+      (l2 - r2)++  {-# INLINE (*) #-}+  IVec2 l1 l2 * IVec2 r1 r2 =+    IVec2+      (l1 * r1)+      (l2 * r2)++  {-# INLINE abs #-}+  abs (IVec2 a b) =+    IVec2 (abs a) (abs b)++  {-# INLINE signum #-}+  signum v2 = withIVec2 v2 \a b ->+    ivec2 (signum a) (signum b)++  {-# INLINE fromInteger #-}+  fromInteger x = IVec2 x' x'+    where+      x' = fromInteger x++instance Storable IVec2 where+  {-# INLINE sizeOf #-}+  sizeOf _ = 8++  {-# INLINE alignment #-}+  alignment _ = 8++  {-# INLINE poke #-}+  poke ptr v4 =+    withIVec2 v4 \a b -> do+      pokeByteOff ptr 0 a+      pokeByteOff ptr 4 b++  {-# INLINE peek #-}+  peek ptr = ivec2+    <$> peekByteOff ptr 0+    <*> peekByteOff ptr 4
+ src/Geomancy/IVec3.hs view
@@ -0,0 +1,131 @@+{-# LANGUAGE BlockArguments #-}+{-# LANGUAGE GeneralizedNewtypeDeriving #-}+{-# LANGUAGE PatternSynonyms #-}+{-# LANGUAGE ViewPatterns #-}++-- | Specialized and inlined @V2 Int32@.++module Geomancy.IVec3+  ( IVec3+  , ivec3+  , withIVec3+  , pattern WithIVec3+  , fromTuple++  , Packed(..)+  , packed+  ) where++import Control.DeepSeq (NFData(rnf))+import Data.Int (Int32)+import Foreign (Storable(..))++data IVec3 = IVec3+  {-# UNPACK #-} !Int32+  {-# UNPACK #-} !Int32+  {-# UNPACK #-} !Int32+  deriving (Eq, Ord, Show)++{-# INLINE ivec3 #-}+ivec3 :: Int32 -> Int32 -> Int32 -> IVec3+ivec3 = IVec3++{-# INLINE withIVec3 #-}+withIVec3+  :: IVec3+  -> (Int32 -> Int32 -> Int32 -> r)+  -> r+withIVec3 (IVec3 a b c) f = f a b c++pattern WithIVec3 :: Int32 -> Int32 -> Int32 -> IVec3+pattern WithIVec3 a b c <- ((`withIVec3` (,,)) -> (a, b, c))+{-# COMPLETE WithIVec3 #-}++{-# INLINE fromTuple #-}+fromTuple :: (Int32, Int32, Int32) -> IVec3+fromTuple (a, b, c) = ivec3 a b c++instance NFData IVec3 where+  rnf IVec3{} = ()++-- XXX: That's another nasty instance...+instance Num IVec3 where+  {-# INLINE (+) #-}+  IVec3 l1 l2 l3 + IVec3 r1 r2 r3 =+    IVec3+      (l1 + r1)+      (l2 + r2)+      (l3 + r3)++  {-# INLINE (-) #-}+  IVec3 l1 l2 l3 - IVec3 r1 r2 r3 =+    IVec3+      (l1 - r1)+      (l2 - r2)+      (l3 - r3)++  {-# INLINE (*) #-}+  IVec3 l1 l2 l3 * IVec3 r1 r2 r3 =+    IVec3+      (l1 * r1)+      (l2 * r2)+      (l3 * r3)++  {-# INLINE abs #-}+  abs x = x++  {-# INLINE signum #-}+  signum v3 = withIVec3 v3 \a b c ->+    ivec3 (signum a) (signum b) (signum c)++  {-# INLINE fromInteger #-}+  fromInteger x = IVec3 x' x' x'+    where+      x' = fromInteger x++instance Storable IVec3 where+  {-# INLINE sizeOf #-}+  sizeOf _ = 16++  {-# INLINE alignment #-}+  alignment _ = 8++  {-# INLINE poke #-}+  poke ptr v3 =+    withIVec3 v3 \a b c -> do+      pokeByteOff ptr  0 a+      pokeByteOff ptr  4 b+      pokeByteOff ptr  8 c++  {-# INLINE peek #-}+  peek ptr = ivec3+    <$> peekByteOff ptr  0+    <*> peekByteOff ptr  4+    <*> peekByteOff ptr  8++newtype Packed = Packed { unPacked :: IVec3 }+  deriving (Eq, Ord, Show, NFData, Num)++{-# INLINE packed #-}+packed :: Int32 -> Int32 -> Int32 -> Packed+packed a b c = Packed (ivec3 a b c)++instance Storable Packed where+  {-# INLINE sizeOf #-}+  sizeOf _ = 12++  {-# INLINE alignment #-}+  alignment _ = 8++  {-# INLINE poke #-}+  poke ptr (Packed v3) =+    withIVec3 v3 \a b c -> do+      pokeByteOff ptr 0 a+      pokeByteOff ptr 4 b+      pokeByteOff ptr 8 c++  {-# INLINE peek #-}+  peek ptr = packed+    <$> peekByteOff ptr 0+    <*> peekByteOff ptr 4+    <*> peekByteOff ptr 8
+ src/Geomancy/IVec4.hs view
@@ -0,0 +1,106 @@+{-# LANGUAGE BlockArguments #-}+{-# LANGUAGE PatternSynonyms #-}+{-# LANGUAGE ViewPatterns #-}++-- | Specialized and inlined @V2 Int32@.++module Geomancy.IVec4+  ( IVec4+  , ivec4+  , withIVec4+  , pattern WithIVec4+  , fromTuple+  ) where++import Control.DeepSeq (NFData(rnf))+import Data.Int (Int32)+import Foreign (Storable(..))++data IVec4 = IVec4+  {-# UNPACK #-} !Int32+  {-# UNPACK #-} !Int32+  {-# UNPACK #-} !Int32+  {-# UNPACK #-} !Int32+  deriving (Eq, Ord, Show)++{-# INLINE ivec4 #-}+ivec4 :: Int32 -> Int32 -> Int32 -> Int32 -> IVec4+ivec4 = IVec4++{-# INLINE withIVec4 #-}+withIVec4+  :: IVec4+  -> (Int32 -> Int32 -> Int32 -> Int32 -> r)+  -> r+withIVec4 (IVec4 a b c d) f = f a b c d++pattern WithIVec4 :: Int32 -> Int32 -> Int32 -> Int32 -> IVec4+pattern WithIVec4 a b c d <- ((`withIVec4` (,,,)) -> (a, b, c, d))+{-# COMPLETE WithIVec4 #-}++{-# INLINE fromTuple #-}+fromTuple :: (Int32, Int32, Int32, Int32) -> IVec4+fromTuple (x, y, z, w) = ivec4 x y z w++instance NFData IVec4 where+  rnf IVec4{} = ()++-- XXX: That's another nasty instance...+instance Num IVec4 where+  {-# INLINE (+) #-}+  IVec4 l1 l2 l3 l4 + IVec4 r1 r2 r3 r4 =+    IVec4+      (l1 + r1)+      (l2 + r2)+      (l3 + r3)+      (l4 + r4)++  {-# INLINE (-) #-}+  IVec4 l1 l2 l3 l4 - IVec4 r1 r2 r3 r4 =+    IVec4+      (l1 - r1)+      (l2 - r2)+      (l3 - r3)+      (l4 - r4)++  {-# INLINE (*) #-}+  IVec4 l1 l2 l3 l4 * IVec4 r1 r2 r3 r4 =+    IVec4+      (l1 * r1)+      (l2 * r2)+      (l3 * r3)+      (l4 * r4)++  {-# INLINE abs #-}+  abs x = x++  {-# INLINE signum #-}+  signum v4 = withIVec4 v4 \a b c d ->+    ivec4 (signum a) (signum b) (signum c) (signum d)++  {-# INLINE fromInteger #-}+  fromInteger x = IVec4 x' x' x' x'+    where+      x' = fromInteger x++instance Storable IVec4 where+  {-# INLINE sizeOf #-}+  sizeOf _ = 16++  {-# INLINE alignment #-}+  alignment _ = 8++  {-# INLINE poke #-}+  poke ptr v4 =+    withIVec4 v4 \a b c d -> do+      pokeByteOff ptr  0 a+      pokeByteOff ptr  4 b+      pokeByteOff ptr  8 c+      pokeByteOff ptr 12 d++  {-# INLINE peek #-}+  peek ptr = ivec4+    <$> peekByteOff ptr  0+    <*> peekByteOff ptr  4+    <*> peekByteOff ptr  8+    <*> peekByteOff ptr 12
+ src/Geomancy/Mat4.hs view
@@ -0,0 +1,498 @@+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE BlockArguments #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE MagicHash #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE UnboxedTuples #-}+{-# LANGUAGE UnliftedFFITypes #-}++-- | General matrix storage and operations.++module Geomancy.Mat4+  ( Mat4++  , rowMajor+  , withRowMajor+  , toListRowMajor+  , toListRowMajor2d+  , fromRowMajor2d++  , colMajor+  , withColMajor+  , toListColMajor+  , toListColMajor2d++  , identity+  , transpose+  , inverse+  , pointwise+  , zipWith+  , matrixProduct+  , scalarMultiply+  , (!*)+  ) where++import Prelude hiding (zipWith)+import GHC.Exts hiding (toList)++import Control.DeepSeq (NFData(rnf))+import Foreign (Storable(..))+import GHC.IO (IO(..))+import System.IO.Unsafe (unsafePerformIO)+import Text.Printf (printf)++import qualified Data.Foldable as Foldable+import qualified Data.List as List++import Geomancy.Vec4 (Vec4, vec4, withVec4)++data Mat4 = Mat4 ByteArray#++{- | Construct 'Mat4' from @row@ notation.+-}+rowMajor+  :: Coercible Mat4 a+  => Float -> Float -> Float -> Float+  -> Float -> Float -> Float -> Float+  -> Float -> Float -> Float -> Float+  -> Float -> Float -> Float -> Float+  -> a+rowMajor = coerce mat4++{- | Reduce 'Mat4' with a function with @row@ notation of arguments.+-}+withRowMajor+  :: Coercible a Mat4+  => a+  ->+    ( Float -> Float -> Float -> Float ->+      Float -> Float -> Float -> Float ->+      Float -> Float -> Float -> Float ->+      Float -> Float -> Float -> Float ->+      r+    )+  -> r+withRowMajor m = withMat4 (coerce m)++toListRowMajor :: Coercible a Mat4 => a -> [Float]+toListRowMajor = toList . coerce++toListRowMajor2d :: Coercible a Mat4 => a -> [[Float]]+toListRowMajor2d = toList2d . coerce++{- |+  Build a Mat4 from a list-of-lists kind of container+  with row-major ordering of elements.++@+  fromRowMajor2d (Linear.mkTransformation dir pos) :: Transform+@+-}+fromRowMajor2d+  :: forall t a+  .  ( Foldable t+     , Coercible Mat4 a+     )+  => t (t Float)+  -> Maybe a+fromRowMajor2d rows =+  case Foldable.toList rows of+    [r0, r1, r2, r3] ->+      withRow r0 \m00 m01 m02 m03 ->+      withRow r1 \m10 m11 m12 m13 ->+      withRow r2 \m20 m21 m22 m23 ->+      withRow r3 \m30 m31 m32 m33 ->+        Just . coerce $ mat4+          m00 m01 m02 m03+          m10 m11 m12 m13+          m20 m21 m22 m23+          m30 m31 m32 m33+    _ ->+      Nothing+  where+    withRow row f =+      case Foldable.toList row of+        [c0, c1, c2, c3] ->+          f c0 c1 c2 c3+        _ ->+          Nothing++{- | Construct a 'Mat4' from @column@ notation.+-}+{-# INLINE colMajor #-}+colMajor+  :: Coercible Mat4 a+  => Float -> Float -> Float -> Float+  -> Float -> Float -> Float -> Float+  -> Float -> Float -> Float -> Float+  -> Float -> Float -> Float -> Float+  -> a+colMajor+  m00 m01 m02 m03+  m10 m11 m12 m13+  m20 m21 m22 m23+  m30 m31 m32 m33 =+    coerce $ mat4+      m00 m10 m20 m30+      m01 m11 m21 m31+      m02 m12 m22 m32+      m03 m13 m23 m33++{- | Reduce 'Mat4' with a function with @column@ notation for arguments.+-}+{-# INLINE withColMajor #-}+withColMajor+  :: Coercible a Mat4+  => a+  ->+    ( Float -> Float -> Float -> Float ->+      Float -> Float -> Float -> Float ->+      Float -> Float -> Float -> Float ->+      Float -> Float -> Float -> Float ->+      r+    )+  -> r+withColMajor m f = withMat4 (coerce m)+  \ m00 m01 m02 m03+    m10 m11 m12 m13+    m20 m21 m22 m23+    m30 m31 m32 m33 ->+  f+    m00 m10 m20 m30+    m01 m11 m21 m31+    m02 m12 m22 m32+    m03 m13 m23 m33++toListColMajor :: Coercible a Mat4 => a -> [Float]+toListColMajor = toListTrans . coerce++toListColMajor2d :: Coercible a Mat4 => a -> [[Float]]+toListColMajor2d = toList2dTrans . coerce++{- | Construct 'Mat4' from elements in memory order.+-}+{-# INLINE mat4 #-}+mat4+  :: Float -> Float -> Float -> Float+  -> Float -> Float -> Float -> Float+  -> Float -> Float -> Float -> Float+  -> Float -> Float -> Float -> Float+  -> Mat4+mat4+  (F# m00) (F# m01) (F# m02) (F# m03)+  (F# m10) (F# m11) (F# m12) (F# m13)+  (F# m20) (F# m21) (F# m22) (F# m23)+  (F# m30) (F# m31) (F# m32) (F# m33) =+  runRW# \world ->+    let+      !(# world_, arr #) = newAlignedPinnedByteArray# 64# 16# world++      world00 = writeFloatArray# arr 0x0# m00 world_+      world01 = writeFloatArray# arr 0x1# m01 world00+      world02 = writeFloatArray# arr 0x2# m02 world01+      world03 = writeFloatArray# arr 0x3# m03 world02++      world10 = writeFloatArray# arr 0x4# m10 world03+      world11 = writeFloatArray# arr 0x5# m11 world10+      world12 = writeFloatArray# arr 0x6# m12 world11+      world13 = writeFloatArray# arr 0x7# m13 world12++      world20 = writeFloatArray# arr 0x8# m20 world13+      world21 = writeFloatArray# arr 0x9# m21 world20+      world22 = writeFloatArray# arr 0xA# m22 world21+      world23 = writeFloatArray# arr 0xB# m23 world22++      world30 = writeFloatArray# arr 0xC# m30 world23+      world31 = writeFloatArray# arr 0xD# m31 world30+      world32 = writeFloatArray# arr 0xE# m32 world31+      world33 = writeFloatArray# arr 0xF# m33 world32++      !(# _world', arr' #) = unsafeFreezeByteArray# arr world33+    in+      Mat4 arr'++{- | Reduce 'Mat4' with a function with @memory@ notation for arguments.+-}+{-# INLINE withMat4 #-}+withMat4+  :: Mat4+  ->+    ( Float -> Float -> Float -> Float ->+      Float -> Float -> Float -> Float ->+      Float -> Float -> Float -> Float ->+      Float -> Float -> Float -> Float ->+      r+    )+  -> r+withMat4 (Mat4 arr) f =+  f+    (F# (indexFloatArray# arr 0x0#))+    (F# (indexFloatArray# arr 0x1#))+    (F# (indexFloatArray# arr 0x2#))+    (F# (indexFloatArray# arr 0x3#))++    (F# (indexFloatArray# arr 0x4#))+    (F# (indexFloatArray# arr 0x5#))+    (F# (indexFloatArray# arr 0x6#))+    (F# (indexFloatArray# arr 0x7#))++    (F# (indexFloatArray# arr 0x8#))+    (F# (indexFloatArray# arr 0x9#))+    (F# (indexFloatArray# arr 0xA#))+    (F# (indexFloatArray# arr 0xB#))++    (F# (indexFloatArray# arr 0xC#))+    (F# (indexFloatArray# arr 0xD#))+    (F# (indexFloatArray# arr 0xE#))+    (F# (indexFloatArray# arr 0xF#))++{- | @I@, the identity matrix.++Neutral element of its monoid, so you can use 'mempty'.+-}+{-# INLINE identity #-}+identity :: Mat4+identity = mat4+  1 0 0 0+  0 1 0 0+  0 0 1 0+  0 0 0 1++{-# INLINE transpose #-}+transpose :: Mat4 -> Mat4+transpose =+  flip withMat4+    \ m00 m01 m02 m03+      m10 m11 m12 m13+      m20 m21 m22 m23+      m30 m31 m32 m33 ->+    mat4+      m00 m10 m20 m30+      m01 m11 m21 m31+      m02 m12 m22 m32+      m03 m13 m23 m33++{- | Compute an inverse matrix, slowly.+-}+inverse :: (Coercible Mat4 a, Coercible Mat4 a) => a -> a+inverse m =+  coerce $ withMat4 (coerce m)+    \ m00 m01 m02 m03+      m10 m11 m12 m13+      m20 m21 m22 m23+      m30 m31 m32 m33 ->+        let+          invDet = recip det++          det+            = s0 * c5+            - s1 * c4+            + s2 * c3+            + s3 * c2+            - s4 * c1+            + s5 * c0++          s0 = m00 * m11 - m10 * m01+          s1 = m00 * m12 - m10 * m02+          s2 = m00 * m13 - m10 * m03+          s3 = m01 * m12 - m11 * m02+          s4 = m01 * m13 - m11 * m03+          s5 = m02 * m13 - m12 * m03++          c5 = m22 * m33 - m32 * m23+          c4 = m21 * m33 - m31 * m23+          c3 = m21 * m32 - m31 * m22+          c2 = m20 * m33 - m30 * m23+          c1 = m20 * m32 - m30 * m22+          c0 = m20 * m31 - m30 * m21++          i00 = ( m11 * c5 - m12 * c4 + m13 * c3) * invDet+          i01 = (-m01 * c5 + m02 * c4 - m03 * c3) * invDet+          i02 = ( m31 * s5 - m32 * s4 + m33 * s3) * invDet+          i03 = (-m21 * s5 + m22 * s4 - m23 * s3) * invDet++          i10 = (-m10 * c5 + m12 * c2 - m13 * c1) * invDet+          i11 = ( m00 * c5 - m02 * c2 + m03 * c1) * invDet+          i12 = (-m30 * s5 + m32 * s2 - m33 * s1) * invDet+          i13 = ( m20 * s5 - m22 * s2 + m23 * s1) * invDet++          i20 = ( m10 * c4 - m11 * c2 + m13 * c0) * invDet+          i21 = (-m00 * c4 + m01 * c2 - m03 * c0) * invDet+          i22 = ( m30 * s4 - m31 * s2 + m33 * s0) * invDet+          i23 = (-m20 * s4 + m21 * s2 - m23 * s0) * invDet++          i30 = (-m10 * c3 + m11 * c1 - m12 * c0) * invDet+          i31 = ( m00 * c3 - m01 * c1 + m02 * c0) * invDet+          i32 = (-m30 * s3 + m31 * s1 - m32 * s0) * invDet+          i33 = ( m20 * s3 - m21 * s1 + m22 * s0) * invDet+        in+          mat4+            i00 i01 i02 i03+            i10 i11 i12 i13+            i20 i21 i22 i23+            i30 i31 i32 i33++pointwise :: Mat4 -> Mat4 -> (Float -> Float -> Float) -> Mat4+pointwise a b f =+  withMat4 a+    \ a00 a01 a02 a03+      a10 a11 a12 a13+      a20 a21 a22 a23+      a30 a31 a32 a33 ->+  withMat4 b+    \ b00 b01 b02 b03+      b10 b11 b12 b13+      b20 b21 b22 b23+      b30 b31 b32 b33 ->+  mat4+    (f a00 b00) (f a01 b01) (f a02 b02) (f a03 b03)+    (f a10 b10) (f a11 b11) (f a12 b12) (f a13 b13)+    (f a20 b20) (f a21 b21) (f a22 b22) (f a23 b23)+    (f a30 b30) (f a31 b31) (f a32 b32) (f a33 b33)++toList :: Mat4 -> [Float]+toList = flip withMat4+    \ a00 a01 a02 a03+      a10 a11 a12 a13+      a20 a21 a22 a23+      a30 a31 a32 a33 ->+    [ a00, a01, a02, a03+    , a10, a11, a12, a13+    , a20, a21, a22, a23+    , a30, a31, a32, a33+    ]++toList2d :: Mat4 -> [[Float]]+toList2d = flip withMat4+    \ a00 a01 a02 a03+      a10 a11 a12 a13+      a20 a21 a22 a23+      a30 a31 a32 a33 ->+    [ [a00, a01, a02, a03]+    , [a10, a11, a12, a13]+    , [a20, a21, a22, a23]+    , [a30, a31, a32, a33]+    ]++toListTrans :: Mat4 -> [Float]+toListTrans = flip withMat4+    \ a00 a01 a02 a03+      a10 a11 a12 a13+      a20 a21 a22 a23+      a30 a31 a32 a33 ->+    [ a00, a10, a20, a30+    , a01, a11, a21, a31+    , a02, a12, a22, a32+    , a03, a13, a23, a33+    ]++toList2dTrans :: Mat4 -> [[Float]]+toList2dTrans = flip withMat4+    \ a00 a01 a02 a03+      a10 a11 a12 a13+      a20 a21 a22 a23+      a30 a31 a32 a33 ->+    [ [a00, a10, a20, a30]+    , [a01, a11, a21, a31]+    , [a02, a12, a22, a32]+    , [a03, a13, a23, a33]+    ]++zipWith :: (Float -> Float -> c) -> Mat4 -> Mat4 -> [c]+zipWith f a b = List.zipWith f (toList a) (toList b)++{-# INLINE scalarMultiply #-}+scalarMultiply :: Float -> Mat4 -> Mat4+scalarMultiply x m =+  withMat4 m+    \ m00 m01 m02 m03+      m10 m11 m12 m13+      m20 m21 m22 m23+      m30 m31 m32 m33 ->+      mat4+        (m00 * x) (m10 * x) (m20 * x) (m30 * x)+        (m01 * x) (m11 * x) (m21 * x) (m31 * x)+        (m02 * x) (m12 * x) (m22 * x) (m32 * x)+        (m03 * x) (m13 * x) (m23 * x) (m33 * x)++-- | Matrix - column vector multiplication+(!*) :: Coercible a Mat4 => a -> Vec4 -> Vec4+(!*) mat vec =+  withVec4 vec \v1 v2 v3 v4 ->+    withColMajor mat+      \ m11 m12 m13 m14+        m21 m22 m23 m24+        m31 m32 m33 m34+        m41 m42 m43 m44 ->+          vec4+            (m11 * v1 + m12 * v2 + m13 * v3 + m14 * v4)+            (m21 * v1 + m22 * v2 + m23 * v3 + m24 * v4)+            (m31 * v1 + m32 * v2 + m33 * v3 + m34 * v4)+            (m41 * v1 + m42 * v2 + m43 * v3 + m44 * v4)++foreign import ccall unsafe "M4x4_SSE" mm4sse :: Addr# -> Addr# -> Addr# -> IO ()++{-# INLINE matrixProduct #-}+matrixProduct :: Mat4 -> Mat4 -> Mat4+matrixProduct (Mat4 l) (Mat4 r) = unsafePerformIO do+  result@(Mat4 m) <- unsafeNewMat4+  mm4sse+    (byteArrayContents# l)+    (byteArrayContents# r)+    (byteArrayContents# m)+  pure result++{-# INLINE unsafeNewMat4 #-}+unsafeNewMat4 :: IO Mat4+unsafeNewMat4 =+  IO \world ->+    let+      !(# world_, arr_ #) = newAlignedPinnedByteArray# 64# 16# world+      !(# _world', arr #) = unsafeFreezeByteArray# arr_ world_+    in+      (# world, Mat4 arr #)++instance NFData Mat4 where+  rnf Mat4{} = ()++instance Semigroup Mat4 where+  {-# INLINE (<>) #-}+  (<>) = matrixProduct++instance Monoid Mat4 where+  {-# INLINE mempty #-}+  mempty = identity++instance Show Mat4 where+  show = flip withMat4+    \ m00 m01 m02 m03+      m10 m11 m12 m13+      m20 m21 m22 m23+      m30 m31 m32 m33 ->+    unlines+      [ printf "| %.4f %.4f %.4f %.4f |" m00 m01 m02 m03+      , printf "| %.4f %.4f %.4f %.4f |" m10 m11 m12 m13+      , printf "| %.4f %.4f %.4f %.4f |" m20 m21 m22 m23+      , printf "| %.4f %.4f %.4f %.4f |" m30 m31 m32 m33+      ]++instance Storable Mat4 where+  sizeOf _mat4 = 64++  alignment _mat4 = 16++  {-# INLINE poke #-}+  poke (Ptr addr) (Mat4 arr) = IO \world ->+    let+      world' = copyByteArrayToAddr# arr 0# addr 64# world+    in+      (# world', () #)++  {-# INLINE peek #-}+  peek (Ptr addr) = IO \world ->+    let+      !(# world0, arr #)  = newAlignedPinnedByteArray# 64# 16# world+      world1              = copyAddrToByteArray# addr arr 0# 64# world0+      !(# world', arr' #) = unsafeFreezeByteArray# arr world1+    in+      (# world', Mat4 arr' #)
+ src/Geomancy/Quaternion.hs view
@@ -0,0 +1,265 @@+{-# LANGUAGE BlockArguments #-}++-- | Specialized and inlined @Quaternion Float@.++module Geomancy.Quaternion+  ( Quaternion+  , quaternion+  , withQuaternion++  , axisAngle+  , rotate+  , rotatePoint+  , rotationBetween+  , lookAtUp++  , (^*)+  , (^/)+  , slerp++  , conjugate+  , norm+  , quadrance+  , dot+  , normalize+  , qNaN+  ) where++import Control.DeepSeq (NFData(rnf))+import Foreign (Storable(..), castPtr)++import Geomancy.Vec3 (Vec3, vec3, withVec3)++import qualified Geomancy.Vec3 as Vec3++data Quaternion = Quaternion+  {-# UNPACK #-} !Float+  {-# UNPACK #-} !Float+  {-# UNPACK #-} !Float+  {-# UNPACK #-} !Float+  deriving (Eq, Ord, Show)++{-# INLINE quaternion #-}+quaternion :: Float -> Float -> Float -> Float -> Quaternion+quaternion = Quaternion++{-# INLINE withQuaternion #-}+withQuaternion+  :: Quaternion+  -> (Float -> Float -> Float -> Float -> r)+  -> r+withQuaternion (Quaternion a b c d) f = f a b c d++{-# INLINE (^*) #-}+(^*) :: Quaternion -> Float -> Quaternion+Quaternion a b c d ^* x =+  Quaternion+    (a * x)+    (b * x)+    (c * x)+    (d * x)++{-# INLINE (^/) #-}+(^/) :: Quaternion -> Float -> Quaternion+Quaternion a b c d ^/ x =+  Quaternion+    (a / x)+    (b / x)+    (c / x)+    (d / x)++slerp :: Quaternion -> Quaternion -> Float -> Quaternion+slerp q p t+  | 1.0 - cosphi < 1e-8 =+      q+  | otherwise =+      ( (q   ^* sin ((1 - t) * phi)) ++         f p ^* sin (t * phi)+      ) ^/ sin phi+  where+    phi = acos cosphi++    (cosphi, f) =+      if dqp < 0 then+        (-dqp, negate)+      else+        (dqp, id)++    dqp = dot q p++{-# INLINE conjugate #-}+conjugate :: Quaternion -> Quaternion+conjugate (Quaternion e x y z) = Quaternion e (-x) (-y) (-z)++{-# INLINE norm #-}+norm :: Quaternion -> Float+norm = sqrt . quadrance++{-# INLINE quadrance #-}+quadrance :: Quaternion -> Float+quadrance q = dot q q++{-# INLINE dot #-}+dot :: Quaternion -> Quaternion -> Float+dot (Quaternion a b c d) (Quaternion e f g h) =+  a * e ++  b * f ++  c * g ++  d * h -- XXX: SIMD time!++{-# INLINE normalize #-}+normalize :: Quaternion -> Quaternion+normalize v =+  if nearZero q || nearZero (1-q) then+    v+  else+    let+      Quaternion e i j k = v+    in+      Quaternion (e / l) (i / l) (j / l) (k / l)++  where+    q = dot v v+    l = sqrt q++    nearZero a = abs a <= 1e-6++instance NFData Quaternion where+  rnf Quaternion{} = ()++instance Num Quaternion where+  {-# INLINE (+) #-}+  Quaternion a b c d + Quaternion e f g h =+    Quaternion+      (a + e)+      (b + f)+      (c + g)+      (d + h)++  {-# INLINE (-) #-}+  Quaternion a b c d - Quaternion e f g h =+    Quaternion+      (a - e)+      (b - f)+      (c - g)+      (d - h)++  {-# INLINE (*) #-}+  Quaternion a b c d * Quaternion e f g h =+    withVec3 v \y z w ->+      Quaternion x y z w+    where+      x = a * e - Vec3.dot v1 v2+      v = Vec3.cross v1 v2 + v2 Vec3.^* a + v1 Vec3.^* e+      v1 = vec3 b c d+      v2 = vec3 f g h++  {-# INLINE fromInteger #-}+  fromInteger x = Quaternion (fromInteger x) 0 0 0++  {-# INLINE abs #-}+  abs z = Quaternion (norm z) 0 0 0++  {-# INLINE signum #-}+  signum q@(Quaternion e i j k)+    | m == 0 =+        q+    | not (isInfinite m || isNaN m) =+        Quaternion (e * misqrt) (i * misqrt) (j * misqrt) (k * misqrt)+    | any isNaN [e, i, j, k] = qNaN+    | not (ii || ij || ik) = Quaternion 1 0 0 0+    | not (ie || ij || ik) = Quaternion 0 1 0 0+    | not (ie || ii || ik) = Quaternion 0 0 1 0+    | not (ie || ii || ij) = Quaternion 0 0 0 1+    | otherwise = qNaN+    where+      m = quadrance q+      misqrt = recip (sqrt m)++      ie = isInfinite e+      ii = isInfinite i+      ij = isInfinite j+      ik = isInfinite k++{-# INLINE qNaN #-}+qNaN :: Quaternion+qNaN = Quaternion fNaN fNaN fNaN fNaN+  where+    fNaN = 0/0++-- XXX: GPU layouts call for some padding.+instance Storable Quaternion where+  {-# INLINE sizeOf #-}+  sizeOf _ = 16++  {-# INLINE alignment #-}+  alignment _ = 16++  {-# INLINE poke #-}+  poke ptr (Quaternion a b c d) = do+    poke ptr' a+    pokeElemOff ptr' 1 b+    pokeElemOff ptr' 2 c+    pokeElemOff ptr' 3 d+    where+      ptr' = castPtr ptr++  {-# INLINE peek #-}+  peek ptr = Quaternion+    <$> peek ptr'+    <*> peekElemOff ptr' 1+    <*> peekElemOff ptr' 2+    <*> peekElemOff ptr' 3+    where+      ptr' = castPtr ptr++-- | Quaternion construction from axis and angle.+{-# INLINE axisAngle #-}+axisAngle :: Vec3 -> Float -> Quaternion+axisAngle axis rads =+  withVec3 (Vec3.normalize axis Vec3.^* sin half) $+    quaternion (cos half)+  where+    half = rads / 2++{-# INLINE rotate #-}+rotate :: Quaternion -> Vec3 -> Vec3+rotate q v = withQuaternion q' \_a b c d -> vec3 b c d+  where+    q' = withVec3 v \x y z ->+      q * quaternion 0 x y z * conjugate q++{-# INLINE rotatePoint #-}+rotatePoint :: Quaternion -> Vec3 -> Vec3 -> Vec3+rotatePoint q origin point =+  origin + rotate q (point - origin)++{- | Rotation between vectors.++(in other words: the quaternion needed to rotate @v1@ so that it matches @v2@)+-}+rotationBetween :: Vec3 -> Vec3 -> Quaternion+rotationBetween v1 v2 = axisAngle axis angle+  where+    axis = Vec3.cross v1 v2+    angle = acos cosAngle+    cosAngle =+      max (-1) . min 1 $+        Vec3.dot (Vec3.normalize v1) (Vec3.normalize v2)++{- | Orient towards a point.++Use "rotationBetween" if you don't need to keep the object upright.+-}+lookAtUp :: Vec3 -> Vec3 -> Vec3 -> Quaternion+lookAtUp src dst up = rot2 * rot1+  where+    dir3 = dst - src++    -- XXX: turn "eye"+    rot1 = rotationBetween (vec3 0 0 1) dir3++    rot2 = rotationBetween newUp fixedUp++    newUp = rotate rot1 up+    fixedUp = Vec3.cross (Vec3.cross dir3 up) dir3
+ src/Geomancy/Transform.hs view
@@ -0,0 +1,178 @@+{-# LANGUAGE BlockArguments #-}+{-# LANGUAGE DerivingStrategies #-}+{-# LANGUAGE GeneralisedNewtypeDeriving #-}++module Geomancy.Transform+  ( Transform(..)+  , inverse++  , apply+  , (!.)++  , translate+  , translateV++  , rotateX+  , rotateY+  , rotateZ+  , rotateQ++  , scale+  , scaleX+  , scaleY+  , scaleZ+  , scaleXY+  , scale3++  , dirPos+  ) where++import Foreign (Storable(..))++import Geomancy.Mat4 (Mat4, colMajor, withColMajor, inverse)+import Geomancy.Quaternion (Quaternion, withQuaternion)+import Geomancy.Vec3 (Vec3, vec3, withVec3)++newtype Transform = Transform { unTransform :: Mat4 }+  deriving newtype (Show, Semigroup, Monoid, Storable)++-- | Apply transformation to a vector, then normalize with perspective division+apply :: Vec3 -> Transform -> Vec3+apply = flip (!.)++-- | Matrix - column vector multiplication with perspective division+(!.) :: Transform -> Vec3 -> Vec3+(!.) mat vec =+  withVec3 vec \v1 v2 v3 ->+    withColMajor mat+      \ m11 m12 m13 m14+        m21 m22 m23 m24+        m31 m32 m33 m34+        m41 m42 m43 m44 ->+          let+            px = m11 * v1 + m12 * v2 + m13 * v3 + m14+            py = m21 * v1 + m22 * v2 + m23 * v3 + m24+            pz = m31 * v1 + m32 * v2 + m33 * v3 + m34+            p  = m41 * v1 + m42 * v2 + m43 * v3 + m44+          in+            vec3 (px / p) (py / p) (pz / p)++-- ** Translation++{-# INLINE translate #-}+translate :: Float -> Float -> Float -> Transform+translate x y z = colMajor+  1 0 0 x+  0 1 0 y+  0 0 1 z+  0 0 0 1++{-# INLINE translateV #-}+translateV :: Vec3 -> Transform+translateV vec = withVec3 vec translate++-- ** Scaling++{-# INLINE scale3 #-}+scale3 :: Float -> Float -> Float -> Transform+scale3 x y z = colMajor+  x 0 0 0+  0 y 0 0+  0 0 z 0+  0 0 0 1++{-# INLINE scale #-}+scale :: Float -> Transform+scale s = scale3 s s s++{-# INLINE scaleX #-}+scaleX :: Float -> Transform+scaleX x = scale3 x 1 1++{-# INLINE scaleY #-}+scaleY :: Float -> Transform+scaleY y = scale3 1 y 1++{-# INLINE scaleZ #-}+scaleZ :: Float -> Transform+scaleZ z = scale3 1 1 z++{-# INLINE scaleXY #-}+scaleXY :: Float -> Float -> Transform+scaleXY x y = scale3 x y 1++-- ** Euler angle rotations++{-# INLINE rotateX #-}+rotateX :: Float -> Transform+rotateX rads = colMajor+  1 0   0   0+  0 t11 t21 0+  0 t12 t22 0+  0 0   0   1+  where+    t11 = cost+    t12 = -sint+    t21 = sint+    t22 = cost++    cost = cos rads+    sint = sin rads++{-# INLINE rotateY #-}+rotateY :: Float -> Transform+rotateY rads = colMajor+  t00 0 t20 0+  0   1 0   0+  t02 0 t22 0+  0   0 0   1+  where+    cost = cos rads+    sint = sin rads++    t00 = cost+    t02 = sint+    t20 = -sint+    t22 = cost++{-# INLINE rotateZ #-}+rotateZ :: Float -> Transform+rotateZ rads = colMajor+  t00 t10 0 0+  t01 t11 0 0+  0   0   1 0+  0   0   0 1+  where+   t00 = cost+   t01 = -sint+   t10 = sint+   t11 = cost++   cost = cos rads+   sint = sin rads++{-# INLINE rotateQ #-}+rotateQ :: Quaternion -> Transform+rotateQ dir = dirPos dir 0++{-# INLINE dirPos #-}+dirPos :: Quaternion -> Vec3 -> Transform+dirPos rs t =+  withQuaternion rs \w x y z ->+  withVec3 t \tx ty tz ->+    let+      x2 = x * x+      y2 = y * y+      z2 = z * z+      xy = x * y+      xz = x * z+      xw = x * w+      yz = y * z+      yw = y * w+      zw = z * w+    in+      colMajor+        (1 - 2 * (y2 + z2)) (    2 * (xy - zw)) (    2 * (xz + yw)) tx+        (    2 * (xy + zw)) (1 - 2 * (x2 + z2)) (    2 * (yz - xw)) ty+        (    2 * (xz - yw)) (    2 * (yz + xw)) (1 - 2 * (x2 + y2)) tz+         0                   0                   0                  1
+ src/Geomancy/UVec2.hs view
@@ -0,0 +1,95 @@+{-# LANGUAGE BlockArguments #-}+{-# LANGUAGE PatternSynonyms #-}+{-# LANGUAGE ViewPatterns #-}++-- | Specialized and inlined @V2 Word32@.++module Geomancy.UVec2+  ( UVec2+  , uvec2+  , withUVec2+  , pattern WithUVec2+  , fromTuple+  ) where++import Control.DeepSeq (NFData(rnf))+import Data.Word (Word32)+import Foreign (Storable(..))++data UVec2 = UVec2+  {-# UNPACK #-} !Word32+  {-# UNPACK #-} !Word32+  deriving (Eq, Ord, Show)++{-# INLINE uvec2 #-}+uvec2 :: Word32 -> Word32 -> UVec2+uvec2 = UVec2++{-# INLINE withUVec2 #-}+withUVec2+  :: UVec2+  -> (Word32 -> Word32 -> r)+  -> r+withUVec2 (UVec2 a b) f = f a b++pattern WithUVec2 :: Word32 -> Word32 -> UVec2+pattern WithUVec2 a b <- ((`withUVec2` (,)) -> (a, b))+{-# COMPLETE WithUVec2 #-}++{-# INLINE fromTuple #-}+fromTuple :: (Word32, Word32) -> UVec2+fromTuple (x, y) = uvec2 x y++instance NFData UVec2 where+  rnf UVec2{} = ()++-- XXX: That's one nasty instance...+instance Num UVec2 where+  {-# INLINE (+) #-}+  UVec2 l1 l2 + UVec2 r1 r2 =+    UVec2+      (l1 + r1)+      (l2 + r2)++  {-# INLINE (-) #-}+  UVec2 l1 l2 - UVec2 r1 r2 =+    UVec2+      (l1 - r1)+      (l2 - r2)++  {-# INLINE (*) #-}+  UVec2 l1 l2 * UVec2 r1 r2 =+    UVec2+      (l1 * r1)+      (l2 * r2)++  {-# INLINE abs #-}+  abs (UVec2 a b) =+    UVec2 (abs a) (abs b)++  {-# INLINE signum #-}+  signum v2 = withUVec2 v2 \a b ->+    uvec2 (signum a) (signum b)++  {-# INLINE fromInteger #-}+  fromInteger x = UVec2 x' x'+    where+      x' = fromInteger x++instance Storable UVec2 where+  {-# INLINE sizeOf #-}+  sizeOf _ = 8++  {-# INLINE alignment #-}+  alignment _ = 8++  {-# INLINE poke #-}+  poke ptr v4 =+    withUVec2 v4 \a b -> do+      pokeByteOff ptr 0 a+      pokeByteOff ptr 4 b++  {-# INLINE peek #-}+  peek ptr = uvec2+    <$> peekByteOff ptr 0+    <*> peekByteOff ptr 4
+ src/Geomancy/UVec3.hs view
@@ -0,0 +1,131 @@+{-# LANGUAGE BlockArguments #-}+{-# LANGUAGE GeneralizedNewtypeDeriving #-}+{-# LANGUAGE PatternSynonyms #-}+{-# LANGUAGE ViewPatterns #-}++-- | Specialized and inlined @V2 Word32@.++module Geomancy.UVec3+  ( UVec3+  , uvec3+  , withUVec3+  , pattern WithUVec3+  , fromTuple++  , Packed(..)+  , packed+  ) where++import Control.DeepSeq (NFData(rnf))+import Data.Word (Word32)+import Foreign (Storable(..))++data UVec3 = UVec3+  {-# UNPACK #-} !Word32+  {-# UNPACK #-} !Word32+  {-# UNPACK #-} !Word32+  deriving (Eq, Ord, Show)++{-# INLINE uvec3 #-}+uvec3 :: Word32 -> Word32 -> Word32 -> UVec3+uvec3 = UVec3++{-# INLINE withUVec3 #-}+withUVec3+  :: UVec3+  -> (Word32 -> Word32 -> Word32 -> r)+  -> r+withUVec3 (UVec3 a b c) f = f a b c++pattern WithUVec3 :: Word32 -> Word32 -> Word32 -> UVec3+pattern WithUVec3 a b c <- ((`withUVec3` (,,)) -> (a, b, c))+{-# COMPLETE WithUVec3 #-}++{-# INLINE fromTuple #-}+fromTuple :: (Word32, Word32, Word32) -> UVec3+fromTuple (a, b, c) = uvec3 a b c++instance NFData UVec3 where+  rnf UVec3{} = ()++-- XXX: That's another nasty instance...+instance Num UVec3 where+  {-# INLINE (+) #-}+  UVec3 l1 l2 l3 + UVec3 r1 r2 r3 =+    UVec3+      (l1 + r1)+      (l2 + r2)+      (l3 + r3)++  {-# INLINE (-) #-}+  UVec3 l1 l2 l3 - UVec3 r1 r2 r3 =+    UVec3+      (l1 - r1)+      (l2 - r2)+      (l3 - r3)++  {-# INLINE (*) #-}+  UVec3 l1 l2 l3 * UVec3 r1 r2 r3 =+    UVec3+      (l1 * r1)+      (l2 * r2)+      (l3 * r3)++  {-# INLINE abs #-}+  abs x = x++  {-# INLINE signum #-}+  signum v3 = withUVec3 v3 \a b c ->+    uvec3 (signum a) (signum b) (signum c)++  {-# INLINE fromInteger #-}+  fromInteger x = UVec3 x' x' x'+    where+      x' = fromInteger x++instance Storable UVec3 where+  {-# INLINE sizeOf #-}+  sizeOf _ = 16++  {-# INLINE alignment #-}+  alignment _ = 8++  {-# INLINE poke #-}+  poke ptr v3 =+    withUVec3 v3 \a b c -> do+      pokeByteOff ptr  0 a+      pokeByteOff ptr  4 b+      pokeByteOff ptr  8 c++  {-# INLINE peek #-}+  peek ptr = uvec3+    <$> peekByteOff ptr  0+    <*> peekByteOff ptr  4+    <*> peekByteOff ptr  8++newtype Packed = Packed { unPacked :: UVec3 }+  deriving (Eq, Ord, Show, NFData, Num)++{-# INLINE packed #-}+packed :: Word32 -> Word32 -> Word32 -> Packed+packed a b c = Packed (uvec3 a b c)++instance Storable Packed where+  {-# INLINE sizeOf #-}+  sizeOf _ = 12++  {-# INLINE alignment #-}+  alignment _ = 8++  {-# INLINE poke #-}+  poke ptr (Packed v3) =+    withUVec3 v3 \a b c -> do+      pokeByteOff ptr 0 a+      pokeByteOff ptr 4 b+      pokeByteOff ptr 8 c++  {-# INLINE peek #-}+  peek ptr = packed+    <$> peekByteOff ptr 0+    <*> peekByteOff ptr 4+    <*> peekByteOff ptr 8
+ src/Geomancy/UVec4.hs view
@@ -0,0 +1,106 @@+{-# LANGUAGE BlockArguments #-}+{-# LANGUAGE PatternSynonyms #-}+{-# LANGUAGE ViewPatterns #-}++-- | Specialized and inlined @V2 Word32@.++module Geomancy.UVec4+  ( UVec4+  , uvec4+  , withUVec4+  , pattern WithUVec4+  , fromTuple+  ) where++import Control.DeepSeq (NFData(rnf))+import Data.Word (Word32)+import Foreign (Storable(..))++data UVec4 = UVec4+  {-# UNPACK #-} !Word32+  {-# UNPACK #-} !Word32+  {-# UNPACK #-} !Word32+  {-# UNPACK #-} !Word32+  deriving (Eq, Ord, Show)++{-# INLINE uvec4 #-}+uvec4 :: Word32 -> Word32 -> Word32 -> Word32 -> UVec4+uvec4 = UVec4++{-# INLINE withUVec4 #-}+withUVec4+  :: UVec4+  -> (Word32 -> Word32 -> Word32 -> Word32 -> r)+  -> r+withUVec4 (UVec4 a b c d) f = f a b c d++pattern WithUVec4 :: Word32 -> Word32 -> Word32 -> Word32 -> UVec4+pattern WithUVec4 a b c d <- ((`withUVec4` (,,,)) -> (a, b, c, d))+{-# COMPLETE WithUVec4 #-}++{-# INLINE fromTuple #-}+fromTuple :: (Word32, Word32, Word32, Word32) -> UVec4+fromTuple (x, y, z, w) = uvec4 x y z w++instance NFData UVec4 where+  rnf UVec4{} = ()++-- XXX: That's another nasty instance...+instance Num UVec4 where+  {-# INLINE (+) #-}+  UVec4 l1 l2 l3 l4 + UVec4 r1 r2 r3 r4 =+    UVec4+      (l1 + r1)+      (l2 + r2)+      (l3 + r3)+      (l4 + r4)++  {-# INLINE (-) #-}+  UVec4 l1 l2 l3 l4 - UVec4 r1 r2 r3 r4 =+    UVec4+      (l1 - r1)+      (l2 - r2)+      (l3 - r3)+      (l4 - r4)++  {-# INLINE (*) #-}+  UVec4 l1 l2 l3 l4 * UVec4 r1 r2 r3 r4 =+    UVec4+      (l1 * r1)+      (l2 * r2)+      (l3 * r3)+      (l4 * r4)++  {-# INLINE abs #-}+  abs x = x++  {-# INLINE signum #-}+  signum v4 = withUVec4 v4 \a b c d ->+    uvec4 (signum a) (signum b) (signum c) (signum d)++  {-# INLINE fromInteger #-}+  fromInteger x = UVec4 x' x' x' x'+    where+      x' = fromInteger x++instance Storable UVec4 where+  {-# INLINE sizeOf #-}+  sizeOf _ = 16++  {-# INLINE alignment #-}+  alignment _ = 8++  {-# INLINE poke #-}+  poke ptr v4 =+    withUVec4 v4 \a b c d -> do+      pokeByteOff ptr  0 a+      pokeByteOff ptr  4 b+      pokeByteOff ptr  8 c+      pokeByteOff ptr 12 d++  {-# INLINE peek #-}+  peek ptr = uvec4+    <$> peekByteOff ptr  0+    <*> peekByteOff ptr  4+    <*> peekByteOff ptr  8+    <*> peekByteOff ptr 12
+ src/Geomancy/Vec2.hs view
@@ -0,0 +1,154 @@+{-# LANGUAGE BlockArguments #-}+{-# LANGUAGE PatternSynonyms #-}+{-# LANGUAGE ViewPatterns #-}++-- | Specialized and inlined @V2 Float@.++module Geomancy.Vec2+  ( Vec2+  , vec2+  , withVec2+  , pattern WithVec2+  , fromTuple++  , (^*)+  , (^/)+  , lerp++  , dot+  , normalize+  ) where++import Control.DeepSeq (NFData(rnf))+import Foreign (Storable(..))++data Vec2 = Vec2+  {-# UNPACK #-} !Float+  {-# UNPACK #-} !Float+  deriving (Eq, Ord, Show)++{-# INLINE vec2 #-}+vec2 :: Float -> Float -> Vec2+vec2 = Vec2++{-# INLINE withVec2 #-}+withVec2+  :: Vec2+  -> (Float -> Float -> r)+  -> r+withVec2 (Vec2 a b) f = f a b++pattern WithVec2 :: Float -> Float -> Vec2+pattern WithVec2 a b <- ((`withVec2` (,)) -> (a, b))+{-# COMPLETE WithVec2 #-}++{-# INLINE fromTuple #-}+fromTuple :: (Float, Float) -> Vec2+fromTuple (x, y) = vec2 x y++instance NFData Vec2 where+  rnf Vec2{} = ()++instance Num Vec2 where+  {-# INLINE (+) #-}+  Vec2 l1 l2 + Vec2 r1 r2 =+    Vec2+      (l1 + r1)+      (l2 + r2)++  {-# INLINE (-) #-}+  Vec2 l1 l2 - Vec2 r1 r2 =+    Vec2+      (l1 - r1)+      (l2 - r2)++  {-# INLINE (*) #-}+  Vec2 l1 l2 * Vec2 r1 r2 =+    Vec2+      (l1 * r1)+      (l2 * r2)++  {-# INLINE abs #-}+  abs (Vec2 a b) =+    Vec2 (abs a) (abs b)++  {-# INLINE signum #-}+  signum (Vec2 a b) =+    Vec2 (signum a) (signum b)++  {-# INLINE fromInteger #-}+  fromInteger x = Vec2 x' x'+    where+      x' = fromInteger x++instance Fractional Vec2 where+  {-# INLINE (/) #-}+  Vec2 l1 l2 / Vec2 r1 r2 =+    Vec2 (l1 / r1) (l2 / r2)++  {-# INLINE recip #-}+  recip (Vec2 a b) =+    Vec2 (recip a) (recip b)++  {-# INLINE fromRational #-}+  fromRational x = Vec2 x' x'+    where+      x' = fromRational x++{-# INLINE (^*) #-}+(^*) :: Vec2 -> Float -> Vec2+Vec2 a b ^* x =+  Vec2+    (a * x)+    (b * x)++{-# INLINE (^/) #-}+(^/) :: Vec2 -> Float -> Vec2+Vec2 a b ^/ x =+  Vec2+    (a / x)+    (b / x)++{-# INLINE lerp #-}+lerp :: Float -> Vec2 -> Vec2 -> Vec2+lerp alpha u v = u ^* alpha + v ^* (1 - alpha)++{-# INLINE dot #-}+dot :: Vec2 -> Vec2 -> Float+dot (Vec2 l1 l2) (Vec2 r1 r2) =+  l1 * r1 + l2 * r2++{-# INLINE normalize #-}+normalize :: Vec2 -> Vec2+normalize v =+  if nearZero q || nearZero (1 - q) then+    v+  else+    let+      Vec2 x y = v+    in+      Vec2 (x / l) (y / l)++  where+    q = dot v v+    l = sqrt q++    nearZero a = abs a <= 1e-6++instance Storable Vec2 where+  {-# INLINE sizeOf #-}+  sizeOf _ = 8++  {-# INLINE alignment #-}+  alignment _ = 8++  {-# INLINE poke #-}+  poke ptr v4 =+    withVec2 v4 \a b -> do+      pokeByteOff ptr 0 a+      pokeByteOff ptr 4 b++  {-# INLINE peek #-}+  peek ptr = vec2+    <$> peekByteOff ptr 0+    <*> peekByteOff ptr 4
+ src/Geomancy/Vec3.hs view
@@ -0,0 +1,229 @@+{-# LANGUAGE BlockArguments #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE GeneralizedNewtypeDeriving #-}+{-# LANGUAGE PatternSynonyms #-}+{-# LANGUAGE ViewPatterns #-}++-- | Specialized and inlined @V3 Float@.++module Geomancy.Vec3+  ( Vec3+  , vec3+  , withVec3+  , pattern WithVec3+  , fromVec2+  , fromTuple++  , (^*)+  , (^/)+  , lerp++  , cross+  , dot+  , normalize++  , Packed(..)+  , packed+  ) where++import Control.DeepSeq (NFData(rnf))+import Data.Coerce (Coercible, coerce)+import Foreign (Storable(..), castPtr)++import Geomancy.Vec2 (Vec2, withVec2)++data Vec3 = Vec3+  {-# UNPACK #-} !Float+  {-# UNPACK #-} !Float+  {-# UNPACK #-} !Float+  deriving (Eq, Ord, Show)++{-# INLINE vec3 #-}+vec3 :: Float -> Float -> Float -> Vec3+vec3 = Vec3++{-# INLINE withVec3 #-}+withVec3+  :: Vec3+  -> (Float -> Float -> Float -> r)+  -> r+withVec3 (Vec3 a b c) f = f a b c++pattern WithVec3 :: Float -> Float -> Float -> Vec3+pattern WithVec3 a b c <- ((`withVec3` (,,)) -> (a, b, c))+{-# COMPLETE WithVec3 #-}++{-# INLINE fromVec2 #-}+fromVec2 :: Coercible Vec3 a => Vec2 -> Float -> a+fromVec2 xy z =+  withVec2 xy \x y ->+    coerce (vec3 x y z)++{-# INLINE fromTuple #-}+fromTuple :: Coercible Vec3 a => (Float, Float, Float) -> a+fromTuple (x, y, z) = coerce (vec3 x y z)++instance NFData Vec3 where+  rnf Vec3{} = ()++instance Num Vec3 where+  {-# INLINE (+) #-}+  Vec3 a b c + Vec3 d e f =+    Vec3+      (a + d)+      (b + e)+      (c + f)++  {-# INLINE (-) #-}+  Vec3 a b c - Vec3 d e f =+    Vec3+      (a - d)+      (b - e)+      (c - f)++  {-# INLINE (*) #-}+  Vec3 a b c * Vec3 d e f =+    Vec3+      (a * d)+      (b * e)+      (c * f)++  {-# INLINE abs #-}+  abs (Vec3 a b c) =+    Vec3 (abs a) (abs b) (abs c)++  {-# INLINE signum #-}+  signum (Vec3 a b c) =+    Vec3 (signum a) (signum b) (signum c)++  {-# INLINE fromInteger #-}+  fromInteger x = Vec3 x' x' x'+    where+      x' = fromInteger x++instance Fractional Vec3 where+  {-# INLINE (/) #-}+  Vec3 l1 l2 l3 / Vec3 r1 r2 r3 =+    Vec3 (l1 / r1) (l2 / r2) (l3 / r3)++  {-# INLINE recip #-}+  recip (Vec3 a b c) =+    Vec3 (recip a) (recip b) (recip c)++  {-# INLINE fromRational #-}+  fromRational x = Vec3 x' x' x'+    where+      x' = fromRational x++{-+  XXX: GPU layouts call for some padding.++  Maybe it would be worth it to flip the sizeOf-s.+-}+instance Storable Vec3 where+  {-# INLINE sizeOf #-}+  sizeOf _ = 16++  {-# INLINE alignment #-}+  alignment _ = 4++  {-# INLINE poke #-}+  poke ptr v3 =+    withVec3 v3 \a b c -> do+      poke ptr' a+      pokeElemOff ptr' 1 b+      pokeElemOff ptr' 2 c+      pokeElemOff ptr' 3 (1.0 :: Float)+    where+      ptr' = castPtr ptr++  {-# INLINE peek #-}+  peek ptr =+    vec3 <$> peek ptr' <*> peekElemOff ptr' 1 <*> peekElemOff ptr' 2+    where+      ptr' = castPtr ptr++{-# INLINE (^*) #-}+(^*) :: Vec3 -> Float -> Vec3+Vec3 a b c ^* x =+  Vec3+    (a * x)+    (b * x)+    (c * x)++{-# INLINE (^/) #-}+(^/) :: Vec3 -> Float -> Vec3+Vec3 a b c ^/ x =+  Vec3+    (a / x)+    (b / x)+    (c / x)++{-# INLINE lerp #-}+lerp :: Float -> Vec3 -> Vec3 -> Vec3+lerp alpha u v = u ^* alpha + v ^* (1 - alpha)++{-# INLINE cross #-}+cross :: Vec3 -> Vec3 -> Vec3+cross (Vec3 a b c) (Vec3 d e f) =+  Vec3+    (b * f - c * e)+    (c * d - a * f)+    (a * e - b * d)++{-# INLINE dot #-}+dot :: Vec3 -> Vec3 -> Float+dot (Vec3 a b c) (Vec3 d e f) =+  a * d ++  b * e ++  c * f++{-# INLINE normalize #-}+normalize :: Vec3 -> Vec3+normalize v =+  if nearZero q || nearZero (1-q) then+    v+  else+    let+      Vec3 x y z = v+    in+      Vec3 (x / l) (y / l) (z / l)++  where+    q = dot v v+    l = sqrt q++    nearZero a = abs a <= 1e-6++-- * Unpadded++newtype Packed = Packed { unPacked :: Vec3 }+  deriving (Eq, Ord, Show, NFData, Num, Fractional)++{-# INLINE packed #-}+packed :: Float -> Float -> Float -> Packed+packed x y z = Packed (vec3 x y z)++instance Storable Packed where+  {-# INLINE sizeOf #-}+  sizeOf _ = 12++  {-# INLINE alignment #-}+  alignment _ = 4++  {-# INLINE poke #-}+  poke ptr (Packed v3) =+    withVec3 v3 \a b c -> do+      poke ptr' a+      pokeElemOff ptr' 1 b+      pokeElemOff ptr' 2 c+    where+      ptr' = castPtr ptr++  {-# INLINE peek #-}+  peek ptr = packed+    <$> peek ptr'+    <*> peekElemOff ptr' 1+    <*> peekElemOff ptr' 2+    where+      ptr' = castPtr ptr
+ src/Geomancy/Vec4.hs view
@@ -0,0 +1,204 @@+{-# LANGUAGE BlockArguments #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE PatternSynonyms #-}+{-# LANGUAGE ViewPatterns #-}++-- | Specialized and inlined @V4 Float@.++module Geomancy.Vec4+  ( Vec4+  , vec4+  , withVec4+  , pattern WithVec4+  , fromVec2+  , fromVec22+  , fromVec3+  , fromTuple++  , (^*)+  , (^/)+  , lerp++  , dot+  , normalize+  ) where++import Control.DeepSeq (NFData(rnf))+import Data.Coerce (Coercible, coerce)+import Foreign (Storable(..), castPtr)++import Geomancy.Vec2 (Vec2, withVec2)+import Geomancy.Vec3 (Vec3, withVec3)++data Vec4 = Vec4+  {-# UNPACK #-} !Float+  {-# UNPACK #-} !Float+  {-# UNPACK #-} !Float+  {-# UNPACK #-} !Float+  deriving (Eq, Ord, Show)++{-# INLINE vec4 #-}+vec4 :: Float -> Float -> Float -> Float -> Vec4+vec4 = Vec4++{-# INLINE withVec4 #-}+withVec4+  :: Vec4+  -> (Float -> Float -> Float -> Float -> r)+  -> r+withVec4 (Vec4 a b c d) f = f a b c d++pattern WithVec4 :: Float -> Float -> Float -> Float -> Vec4+pattern WithVec4 a b c d <- ((`withVec4` (,,,)) -> (a, b, c, d))+{-# COMPLETE WithVec4 #-}++{-# INLINE fromVec2 #-}+fromVec2 :: Vec2 -> Float -> Float -> Vec4+fromVec2 xy z w =+  withVec2 xy \x y ->+    vec4 x y z w++{-# INLINE fromVec22 #-}+fromVec22 :: Vec2 -> Vec2 -> Vec4+fromVec22 xy zw =+  withVec2 xy \x y ->+  withVec2 zw \z w ->+    vec4 x y z w++{-# INLINE fromVec3 #-}+fromVec3 :: Coercible a Vec3 => a -> Float -> Vec4+fromVec3 xyz w =+  withVec3 (coerce xyz) \x y z ->+    vec4 x y z w++{-# INLINE fromTuple #-}+fromTuple :: (Float, Float, Float, Float) -> Vec4+fromTuple (x, y, z, w) = vec4 x y z w++instance NFData Vec4 where+  rnf Vec4{} = ()++instance Num Vec4 where+  {-# INLINE (+) #-}+  Vec4 l1 l2 l3 l4 + Vec4 r1 r2 r3 r4 =+    Vec4+      (l1 + r1)+      (l2 + r2)+      (l3 + r3)+      (l4 + r4)++  {-# INLINE (-) #-}+  Vec4 l1 l2 l3 l4 - Vec4 r1 r2 r3 r4 =+    Vec4+      (l1 - r1)+      (l2 - r2)+      (l3 - r3)+      (l4 - r4)++  {-# INLINE (*) #-}+  Vec4 l1 l2 l3 l4 * Vec4 r1 r2 r3 r4 =+    Vec4+      (l1 * r1)+      (l2 * r2)+      (l3 * r3)+      (l4 * r4)++  {-# INLINE abs #-}+  abs (Vec4 a b c d) =+    Vec4 (abs a) (abs b) (abs c) (abs d)++  {-# INLINE signum #-}+  signum (Vec4 a b c d) =+    Vec4 (signum a) (signum b) (signum c) (signum d)++  {-# INLINE fromInteger #-}+  fromInteger x = Vec4 x' x' x' x'+    where+      x' = fromInteger x++instance Fractional Vec4 where+  {-# INLINE (/) #-}+  Vec4 l1 l2 l3 l4 / Vec4 r1 r2 r3 r4 =+    Vec4 (l1 / r1) (l2 / r2) (l3 / r3) (l4 / r4)++  {-# INLINE recip #-}+  recip (Vec4 a b c d) =+    Vec4 (recip a) (recip b) (recip c) (recip d)++  {-# INLINE fromRational #-}+  fromRational x = Vec4 x' x' x' x'+    where+      x' = fromRational x++instance Storable Vec4 where+  {-# INLINE sizeOf #-}+  sizeOf _ = 16++  {-# INLINE alignment #-}+  alignment _ = 16++  {-# INLINE poke #-}+  poke ptr v4 =+    withVec4 v4 \a b c d -> do+      poke ptr' a+      pokeElemOff ptr' 1 b+      pokeElemOff ptr' 2 c+      pokeElemOff ptr' 3 d+    where+      ptr' = castPtr ptr++  {-# INLINE peek #-}+  peek ptr = vec4+    <$> peek ptr'+    <*> peekElemOff ptr' 1+    <*> peekElemOff ptr' 2+    <*> peekElemOff ptr' 3+    where+      ptr' = castPtr ptr++{-# INLINE (^*) #-}+(^*) :: Vec4 -> Float -> Vec4+Vec4 a b c d ^* x =+  Vec4+    (a * x)+    (b * x)+    (c * x)+    (d * x)++{-# INLINE (^/) #-}+(^/) :: Vec4 -> Float -> Vec4+Vec4 a b c d ^/ x =+  Vec4+    (a / x)+    (b / x)+    (c / x)+    (d / x)++{-# INLINE lerp #-}+lerp :: Float -> Vec4 -> Vec4 -> Vec4+lerp alpha u v = u ^* alpha + v ^* (1 - alpha)++{-# INLINE dot #-}+dot :: Vec4 -> Vec4 -> Float+dot (Vec4 l1 l2 l3 l4) (Vec4 r1 r2 r3 r4) =+  l1 * r1 ++  l2 * r2 ++  l3 * r3 ++  l4 * r4++{-# INLINE normalize #-}+normalize :: Vec4 -> Vec4+normalize v =+  if nearZero q || nearZero (1-q) then+    v+  else+    let+      Vec4 x y z w = v+    in+      Vec4 (x / l) (y / l) (z / l) (w / l)++  where+    q = dot v v+    l = sqrt q++    nearZero a = abs a <= 1e-6
+ src/Geomancy/Vulkan/Projection.hs view
@@ -0,0 +1,72 @@+module Geomancy.Vulkan.Projection+  ( perspective+  , infinitePerspective+  , orthoOffCenter+  ) where++import Geomancy.Mat4 (colMajor)+import Geomancy.Transform (Transform(..))++perspective+  :: Integral side+  => Float+  -> Float -> Float+  -> side -> side+  -> Transform+perspective fovRads near far width height = colMajor+  x 0   0   0+  0 y   0   0+  0 0   z w23+  0 0 w32   1++  where+    x = cotFoV * aspectX+    y = -cotFoV+    z = far / (near - far)+    w23 = near * far / (near - far)+    w32 = -1++    cotFoV = recip . tan $ 0.5 * fovRads++    aspectX = fromIntegral height / fromIntegral width++infinitePerspective+  :: Integral side+  => Float+  -> side+  -> side+  -> Transform+infinitePerspective fovRads width height = colMajor+  x   0   0  0+  0 (-y)  0  0+  0   0 (-1) w+  0   0 (-1) 0+  where+    (x, y) =+      if width > height then+        ( cotFoV / camAspect+        , cotFoV+        )+      else+        ( cotFoV+        , cotFoV * camAspect+        )+    camAspect = fromIntegral width / fromIntegral height+    cotFoV = recip . tan $ 0.5 * fovRads++    w = -2 * near+    near = 1/128 -- 2048++orthoOffCenter :: Integral side => Float -> Float -> side -> side -> Transform+orthoOffCenter near far width height = colMajor+  x 0 0 0+  0 y 0 0+  0 0 z w+  0 0 0 1++  where+    x = 2 / fromIntegral width+    y = 2 / fromIntegral height+    z = 1 / (far - near)++    w = near * (near - far)
+ src/Geomancy/Vulkan/View.hs view
@@ -0,0 +1,41 @@+{-# LANGUAGE BlockArguments #-}++module Geomancy.Vulkan.View+  ( orthoFitScreen+  , lookAt+  ) where++import Geomancy.Mat4 (Mat4, rowMajor)+import Geomancy.Vec3 (Vec3, withVec3)+import Geomancy.Transform (Transform(..))++import qualified Geomancy.Vec3 as Vec3++lookAt :: Vec3 -> Vec3 -> Vec3 -> Transform+lookAt eye center up =+  withVec3 xa \xaX xaY xaZ ->+  withVec3 ya \yaX yaY yaZ ->+  withVec3 za \zaX zaY zaZ ->+  rowMajor+    xaX yaX (-zaX) 0+    xaY yaY (-zaY) 0+    xaZ yaZ (-zaZ) 0+    xd  yd    zd   1+  where+    xa = Vec3.normalize $ Vec3.cross za up+    ya = Vec3.cross xa za+    za = Vec3.normalize $ center - eye++    xd = - Vec3.dot xa eye+    yd = - Vec3.dot ya eye+    zd =   Vec3.dot za eye++orthoFitScreen :: Float -> Float -> Float -> Float -> Mat4+orthoFitScreen screenWidth screenHeight targetWidth targetHeight =+  rowMajor+    s 0 0 0+    0 s 0 0+    0 0 1 0+    0 0 0 1+  where+    s = min (screenWidth / targetWidth) (screenHeight / targetHeight)
+ test/Spec.hs view
@@ -0,0 +1,151 @@+{-# LANGUAGE BlockArguments #-}+{-# LANGUAGE NumericUnderscores #-}+{-# LANGUAGE PatternSynonyms #-}+{-# LANGUAGE TemplateHaskell #-}++import Hedgehog+import qualified Hedgehog.Gen as Gen+import qualified Hedgehog.Range as Range++import Control.Monad (unless)+import Control.Monad.IO.Class (MonadIO(..))+import Data.Foldable (toList)+import Data.Maybe (catMaybes)+import GHC.Stack (withFrozenCallStack)+import System.Exit (exitFailure, exitSuccess)+import Text.Printf (printf)++import Geomancy.Transform (Transform(..))++import qualified Foreign+import qualified Linear+import qualified Geomancy+import qualified Geomancy.Mat4+import qualified Geomancy.Transform as Transform++import Linear ((!*!))++main :: IO ()+main = do+  passed <- checkParallel discovered+  if passed then+    exitSuccess+  else+    exitFailure++-- | Enough tests to cover the principles+pattern PROP_TESTS :: TestLimit+pattern PROP_TESTS = 10_000++-- | Try harder to catch FP precision errors+pattern PROP_TESTS_BRUTAL :: TestLimit+pattern PROP_TESTS_BRUTAL = 10_000_000++prop_assoc_multiply :: Property+prop_assoc_multiply = withTests PROP_TESTS $ property do+  (a, a_) <- forAllTransform+  (b, b_) <- forAllTransform+  (c, c_) <- forAllTransform++  let+    ab'c = (a <> b) <> c+    a'bc = a <> (b <> c)+    delta' = nearlyEqualMat4 ab'c a'bc+  annotateShow delta'++  let+    ab_c = (a_ !*! b_) !*! c_+    a_bc = a_ !*! (b_ !*! c_)+    delta_ = nearlyEqualM44 ab_c a_bc+  annotateShow delta_++  -- Intra-library transitivity+  unless (null $ catMaybes delta') do+    -- XXX: check only if there is some outstanding error+    delta' === delta_++  -- Inter-library calculated values nearlyEqual+  ab'c_ <- toM44 ab'c+  catMaybes (nearlyEqualM44 ab_c ab'c_) === []++  a'bc_ <- toM44 a'bc+  catMaybes (nearlyEqualM44 a_bc a'bc_) === []++forAllTransform :: PropertyT IO (Geomancy.Mat4.Mat4, Linear.M44 Float)+forAllTransform = withFrozenCallStack do+  (_name, Transform g) <- forAllWith fst genTransform+  l <- toM44 g+  pure (g, l)++genTransform :: Gen ([Char], Transform)+genTransform = Gen.choice+  [ genIdentity+  , genTranslate+  , genRotate+  , genScale+  ]+  where+    genIdentity = pure ("identity", mempty)++    genTranslate = do+      x <- Gen.float (Range.linearFracFrom 0.0 (-1e6) 1e6)+      y <- Gen.float (Range.linearFracFrom 0.0 (-1e6) 1e6)+      z <- Gen.float (Range.linearFracFrom 0.0 (-1e6) 1e6)+      pure+        ( printf "translate %0.4f %0.4f %0.4f" x y z+        , Transform.translate x y z+        )++    genRotate = do+      (name, axis) <- Gen.element+        [ ("rotate/x", Transform.rotateX)+        , ("rotate/y", Transform.rotateY)+        , ("rotate/z", Transform.rotateZ)+        ]+      angle <- Gen.float (Range.linearFracFrom 0.0 (-4 * pi) (4 * pi))+      pure+        ( printf "%s %0.4f" name angle+        , axis angle+        )++    genScale = do+      x <- Gen.float (Range.linearFracFrom 1.0 1e-6 1e6)+      y <- Gen.float (Range.linearFracFrom 1.0 1e-6 1e6)+      z <- Gen.float (Range.linearFracFrom 1.0 1e-6 1e6)+      pure+        ( printf "scale %0.4f %0.4f %0.4f" x y z+        , Transform.scale3 x y z+        )++-- toVulkan :: Linear.M44 Float -> Linear.M44 Float+-- toVulkan = {- Linear.transpose . -} (correction Linear.!*!)+--   where+--     correction = Linear.V4+--       (Linear.V4 1   0  0   0)+--       (Linear.V4 0 (-1) 0   0)+--       (Linear.V4 0   0  0.5 0.5)+--       (Linear.V4 0   0  0   1)++toM44 :: MonadIO io => Geomancy.Mat4 -> io (Linear.M44 Float)+toM44 mat4 =+  liftIO $+    Foreign.with mat4 $+      Foreign.peek . Foreign.castPtr++type NearlyEqual = Maybe (Float, Float, Float)++nearlyEqualMat4 :: Geomancy.Mat4.Mat4 -> Geomancy.Mat4.Mat4 -> [NearlyEqual]+nearlyEqualMat4 a b = Geomancy.Mat4.zipWith nearlyEqual a b++nearlyEqualM44 :: Linear.M44 Float -> Linear.M44 Float -> [NearlyEqual]+nearlyEqualM44 a b = zipWith nearlyEqual (concatMap toList a) (concatMap toList b)++nearlyEqual :: Float -> Float -> NearlyEqual+nearlyEqual x y =+  if x == y || abs (1 - x / y) < 0.001 then+    Nothing+  else+    Just (x, y, abs (1 - x / y))++discovered :: Group+discovered = $$(discover)