geomancy 0.2.2.4 → 0.2.3.0
raw patch · 6 files changed
+233/−143 lines, 6 filesPVP ok
version bump matches the API change (PVP)
API changes (from Hackage documentation)
+ Geomancy.Vec4: Vec4 :: ByteArray# -> Vec4
+ Geomancy.Vec4: unsafeNewVec4 :: IO Vec4
Files
- ChangeLog.md +5/−0
- cbits/mat4.c +40/−19
- geomancy.cabal +2/−2
- src/Geomancy/Mat4.hs +43/−33
- src/Geomancy/Transform.hs +8/−14
- src/Geomancy/Vec4.hs +135/−75
ChangeLog.md view
@@ -1,5 +1,10 @@ # Changelog for geomancy +## 0.2.3.0++- Vec4 moved to ByteArray#.+- Mat4 !* Vec4 is now SIMD too.+ ## 0.2.2.4 + Add Transform trees.
cbits/mat4.c view
@@ -1,22 +1,43 @@ #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);- }+void Mat4xMat4_SSE(float *A, float *B, float *O) {+ __m128 row0 = _mm_load_ps(&B[0]);+ __m128 row1 = _mm_load_ps(&B[4]);+ __m128 row2 = _mm_load_ps(&B[8]);+ __m128 row3 = _mm_load_ps(&B[12]);+ for(int i=0; i<4; i++) {+ __m128 brod0 = _mm_set1_ps(A[4*i + 0]);+ __m128 brod1 = _mm_set1_ps(A[4*i + 1]);+ __m128 brod2 = _mm_set1_ps(A[4*i + 2]);+ __m128 brod3 = _mm_set1_ps(A[4*i + 3]);+ __m128 row = _mm_add_ps(+ _mm_add_ps(+ _mm_mul_ps(brod0, row0),+ _mm_mul_ps(brod1, row1)),+ _mm_add_ps(+ _mm_mul_ps(brod2, row2),+ _mm_mul_ps(brod3, row3)));+ _mm_store_ps(&O[4*i], row);+ }+}++void Mat4xVec4_SSE(float *M, float *V, float *O) {+ __m128 row0 = _mm_load_ps(&M[0]);+ __m128 row1 = _mm_load_ps(&M[4]);+ __m128 row2 = _mm_load_ps(&M[8]);+ __m128 row3 = _mm_load_ps(&M[12]);++ __m128 vec0 = _mm_load_ps(V);++ __m128 vec00 = _mm_shuffle_ps(vec0, vec0, _MM_SHUFFLE(0, 0, 0, 0));+ __m128 vec01 = _mm_shuffle_ps(vec0, vec0, _MM_SHUFFLE(1, 1, 1, 1));+ __m128 vec02 = _mm_shuffle_ps(vec0, vec0, _MM_SHUFFLE(2, 2, 2, 2));+ __m128 vec03 = _mm_shuffle_ps(vec0, vec0, _MM_SHUFFLE(3, 3, 3, 3));++ __m128 res = _mm_mul_ps(row0, vec00);+ res = _mm_add_ps(res, _mm_mul_ps(row1, vec01));+ res = _mm_add_ps(res, _mm_mul_ps(row2, vec02));+ res = _mm_add_ps(res, _mm_mul_ps(row3, vec03));++ _mm_storeu_ps(O, res); }
geomancy.cabal view
@@ -4,10 +4,10 @@ -- -- see: https://github.com/sol/hpack ----- hash: 1480b5bc5e86691e54a70d648214d4951a3685c03bc3eeb8741a8429eaa0b13c+-- hash: aa93c14151d6ba606892efc64f6795f6f0ec8ebfc5b9f24e104d5b4676cf8bb8 name: geomancy-version: 0.2.2.4+version: 0.2.3.0 synopsis: Geometry and matrix manipulation description: Sometimes it is unavoidable you have to do stuff on CPU. Let's at least do it faster.
src/Geomancy/Mat4.hs view
@@ -5,6 +5,7 @@ {-# LANGUAGE ScopedTypeVariables #-} {-# LANGUAGE UnboxedTuples #-} {-# LANGUAGE UnliftedFFITypes #-}+{-# LANGUAGE ViewPatterns #-} -- | General matrix storage and operations. @@ -33,7 +34,7 @@ ) where import Prelude hiding (zipWith)-import GHC.Exts hiding (toList)+import GHC.Exts hiding (VecCount(..), toList) import Control.DeepSeq (NFData(rnf)) import Foreign (Storable(..))@@ -44,7 +45,7 @@ import qualified Data.Foldable as Foldable import qualified Data.List as List -import Geomancy.Vec4 (Vec4, vec4, withVec4)+import Geomancy.Vec4 (Vec4(..), unsafeNewVec4) data Mat4 = Mat4 ByteArray# @@ -401,42 +402,13 @@ 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 ()+foreign import ccall unsafe "Mat4xMat4_SSE" m4m4sse :: Addr# -> Addr# -> Addr# -> IO () {-# INLINE matrixProduct #-} matrixProduct :: Mat4 -> Mat4 -> Mat4 matrixProduct (Mat4 l) (Mat4 r) = unsafePerformIO do result@(Mat4 m) <- unsafeNewMat4- mm4sse+ m4m4sse (byteArrayContents# l) (byteArrayContents# r) (byteArrayContents# m)@@ -451,6 +423,44 @@ !(# _world', arr #) = unsafeFreezeByteArray# arr_ world_ in (# world, Mat4 arr #)++{-# 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)++foreign import ccall unsafe "Mat4xVec4_SSE" m4v4sse :: Addr# -> Addr# -> Addr# -> IO ()++-- | Matrix - column vector multiplication+(!*) :: Coercible a Mat4 => a -> Vec4 -> Vec4+(!*) (coerce -> Mat4 m) (Vec4 v) = unsafePerformIO do+ result@(Vec4 o) <- unsafeNewVec4+ m4v4sse+ (byteArrayContents# m)+ (byteArrayContents# v)+ (byteArrayContents# o)+ pure result++ -- 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) instance NFData Mat4 where rnf Mat4{} = ()
src/Geomancy/Transform.hs view
@@ -29,10 +29,13 @@ import Foreign (Storable(..)) -import Geomancy.Mat4 (Mat4, colMajor, withColMajor, inverse)+import Geomancy.Mat4 (Mat4, colMajor, inverse) import Geomancy.Quaternion (Quaternion, withQuaternion) import Geomancy.Vec3 (Vec3, vec3, withVec3)+import Geomancy.Vec4 (fromVec3, withVec4) +import qualified Geomancy.Mat4 as Mat4+ newtype Transform = Transform { unTransform :: Mat4 } deriving newtype (Show, Semigroup, Monoid, Storable) @@ -43,19 +46,10 @@ -- | 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)+ withVec4 res \x y z w ->+ vec3 (x / w) (y / w) (z / w)+ where+ res = mat Mat4.!* fromVec3 vec 1.0 -- ** Translation
src/Geomancy/Vec4.hs view
@@ -1,12 +1,16 @@+{-# LANGUAGE BangPatterns #-} {-# LANGUAGE BlockArguments #-} {-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE MagicHash #-} {-# LANGUAGE PatternSynonyms #-}+{-# LANGUAGE UnboxedTuples #-}+{-# LANGUAGE UnliftedFFITypes #-} {-# LANGUAGE ViewPatterns #-} -- | Specialized and inlined @V4 Float@. module Geomancy.Vec4- ( Vec4+ ( Vec4(..) , vec4 , withVec4 , pattern WithVec4@@ -21,33 +25,74 @@ , dot , normalize++ , unsafeNewVec4 ) where +import GHC.Exts hiding (VecCount(..), toList)+ import Control.DeepSeq (NFData(rnf))-import Data.Coerce (Coercible, coerce)-import Foreign (Storable(..), castPtr)+import Foreign (Storable(..))+import GHC.IO (IO(..))+-- import System.IO.Unsafe (unsafePerformIO)+import Text.Printf (printf) 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)+data Vec4 = Vec4 ByteArray# {-# INLINE vec4 #-} vec4 :: Float -> Float -> Float -> Float -> Vec4-vec4 = Vec4+vec4 (F# v0) (F# v1) (F# v2) (F# v3) =+ runRW# \world ->+ let+ !(# world_, arr #) = newAlignedPinnedByteArray# 16# 16# world + world0 = writeFloatArray# arr 0x0# v0 world_+ world1 = writeFloatArray# arr 0x1# v1 world0+ world2 = writeFloatArray# arr 0x2# v2 world1+ world3 = writeFloatArray# arr 0x3# v3 world2+ !(# _world', arr' #) = unsafeFreezeByteArray# arr world3+ in+ Vec4 arr'+ {-# INLINE withVec4 #-} withVec4 :: Vec4 -> (Float -> Float -> Float -> Float -> r) -> r-withVec4 (Vec4 a b c d) f = f a b c d+withVec4 (Vec4 arr) f =+ f+ (F# (indexFloatArray# arr 0x0#))+ (F# (indexFloatArray# arr 0x1#))+ (F# (indexFloatArray# arr 0x2#))+ (F# (indexFloatArray# arr 0x3#)) +{-# INLINE compareVec4 #-}+compareVec4 :: Vec4 -> Vec4 -> Ordering+compareVec4 (Vec4 src1) (Vec4 src2) =+ compare (I# (compareByteArrays# src1 0# src2 0# 16#)) 0++instance Eq Vec4 where+ (==) a b =+ case compareVec4 a b of+ EQ -> True+ _ -> False++ (/=) a b =+ case compareVec4 a b of+ EQ -> False+ _ -> True++instance Ord Vec4 where+ compare = compareVec4++instance Show Vec4 where+ show v =+ withVec4 v $+ printf "Vec4 %.4f %.4f %.4f %.4f"+ pattern WithVec4 :: Float -> Float -> Float -> Float -> Vec4 pattern WithVec4 a b c d <- ((`withVec4` (,,,)) -> (a, b, c, d)) {-# COMPLETE WithVec4 #-}@@ -80,53 +125,64 @@ instance Num Vec4 where {-# INLINE (+) #-}- Vec4 l1 l2 l3 l4 + Vec4 r1 r2 r3 r4 =- Vec4- (l1 + r1)- (l2 + r2)- (l3 + r3)- (l4 + r4)+ (+) l r =+ withVec4 l \l1 l2 l3 l4 ->+ withVec4 r \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)+ (-) l r =+ withVec4 l \l1 l2 l3 l4 ->+ withVec4 r \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)+ (*) l r =+ withVec4 l \l1 l2 l3 l4 ->+ withVec4 r \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)+ abs v =+ withVec4 v \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)+ signum v =+ withVec4 v \a b c d ->+ vec4 (signum a) (signum b) (signum c) (signum d) {-# INLINE fromInteger #-}- fromInteger x = Vec4 x' x' x' x'+ 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)+ (/) l r =+ withVec4 l \l1 l2 l3 l4 ->+ withVec4 r \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)+ recip v =+ withVec4 v \a b c d ->+ vec4 (recip a) (recip b) (recip c) (recip d) {-# INLINE fromRational #-}- fromRational x = Vec4 x' x' x' x'+ fromRational x = vec4 x' x' x' x' where x' = fromRational x @@ -138,53 +194,50 @@ 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+ poke (Ptr addr) (Vec4 arr) = IO \world ->+ let+ world' = copyByteArrayToAddr# arr 0# addr 16# world+ in+ (# world', () #) {-# INLINE peek #-}- peek ptr = vec4- <$> peek ptr'- <*> peekElemOff ptr' 1- <*> peekElemOff ptr' 2- <*> peekElemOff ptr' 3- where- ptr' = castPtr ptr+ peek (Ptr addr) = IO \world ->+ let+ !(# world0, arr #) = newAlignedPinnedByteArray# 16# 16# world+ world1 = copyAddrToByteArray# addr arr 0# 16# world0+ !(# world', arr' #) = unsafeFreezeByteArray# arr world1+ in+ (# world', Vec4 arr' #) +-- TODO: SIMD {-# INLINE (^*) #-} (^*) :: Vec4 -> Float -> Vec4-Vec4 a b c d ^* x =- Vec4- (a * x)- (b * x)- (c * x)- (d * x)+(^*) v x =+ withVec4 v \a b c d ->+ 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)+v ^/ x = v ^* recip x {-# INLINE lerp #-} lerp :: Float -> Vec4 -> Vec4 -> Vec4 lerp alpha u v = u ^* alpha + v ^* (1 - alpha) +-- TODO: SIMD {-# 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+dot a b =+ withVec4 a \a1 a2 a3 a4 ->+ withVec4 b \b1 b2 b3 b4 ->+ a1 * b1 ++ a2 * b2 ++ a3 * b3 ++ a4 * b4 {-# INLINE normalize #-} normalize :: Vec4 -> Vec4@@ -192,13 +245,20 @@ 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)+ v ^/ l where q = dot v v l = sqrt q nearZero a = abs a <= 1e-6++{-# INLINE unsafeNewVec4 #-}+unsafeNewVec4 :: IO Vec4+unsafeNewVec4 =+ IO \world ->+ let+ !(# world_, arr_ #) = newAlignedPinnedByteArray# 16# 16# world+ !(# _world', arr #) = unsafeFreezeByteArray# arr_ world_+ in+ (# world, Vec4 arr #)