geomancy 0.3.1.0 → 0.4.0.0

raw patch · 95 files changed

+23777/−2151 lines, 95 filesdep +tastydep +tasty-benchdep +tasty-hedgehogdep −criterionPVP ok

version bump matches the API change (PVP)

Dependencies added: tasty, tasty-bench, tasty-hedgehog, vector

Dependencies removed: criterion

API changes (from Hackage documentation)

Files

• ChangeLog.md +57/−0
• LICENSE +1/−1
• README.md +36/−38
• bench/Bench/Common.hs +15/−0
• bench/Bench/DMat4.hs +69/−0
• bench/Bench/DQuaternion.hs +70/−0
• bench/Bench/Mat4.hs +69/−0
• bench/Bench/Quaternion.hs +49/−0
• bench/Main.hs +15/−80
• cbits/dmat4.c +153/−0
• cbits/mat4.c +10/−10
• geomancy.cabal +83/−6
• src/Geomancy.hs +77/−1
• src/Geomancy/DMat4.hs +351/−0
• src/Geomancy/DMat4/Internal.hs +392/−0
• src/Geomancy/DQuaternion.hs +358/−0
• src/Geomancy/DTransform.hs +226/−0
• src/Geomancy/DTransform/Tree.hs +132/−0
• src/Geomancy/DVec2.hs +340/−0
• src/Geomancy/DVec3.hs +510/−0
• src/Geomancy/DVec4.hs +516/−0
• src/Geomancy/I64Vec2.hs +255/−0
• src/Geomancy/I64Vec3.hs +339/−0
• src/Geomancy/I64Vec4.hs +288/−0
• src/Geomancy/IVec2.hs +11/−0
• src/Geomancy/IVec3.hs +31/−3
• src/Geomancy/IVec4.hs +21/−1
• src/Geomancy/Mat4/Internal.hs +63/−28
• src/Geomancy/Point.hs +37/−4
• src/Geomancy/Quaternion.hs +14/−4
• src/Geomancy/Swizzle.hs +83/−1513
• src/Geomancy/Swizzle/DVec/Down.hs +369/−0
• src/Geomancy/Swizzle/DVec/Extend.hs +45/−0
• src/Geomancy/Swizzle/DVec/Shuffle.hs +1161/−0
• src/Geomancy/Swizzle/DVec/Up.hs +433/−0
• src/Geomancy/Swizzle/I64Vec/Down.hs +369/−0
• src/Geomancy/Swizzle/I64Vec/Extend.hs +45/−0
• src/Geomancy/Swizzle/I64Vec/Shuffle.hs +1161/−0
• src/Geomancy/Swizzle/I64Vec/Up.hs +433/−0
• src/Geomancy/Swizzle/IVec/Down.hs +369/−0
• src/Geomancy/Swizzle/IVec/Extend.hs +45/−0
• src/Geomancy/Swizzle/IVec/Shuffle.hs +1161/−0
• src/Geomancy/Swizzle/IVec/Up.hs +433/−0
• src/Geomancy/Swizzle/U64Vec/Down.hs +369/−0
• src/Geomancy/Swizzle/U64Vec/Extend.hs +45/−0
• src/Geomancy/Swizzle/U64Vec/Shuffle.hs +1161/−0
• src/Geomancy/Swizzle/U64Vec/Up.hs +433/−0
• src/Geomancy/Swizzle/UVec/Down.hs +369/−0
• src/Geomancy/Swizzle/UVec/Extend.hs +45/−0
• src/Geomancy/Swizzle/UVec/Shuffle.hs +1161/−0
• src/Geomancy/Swizzle/UVec/Up.hs +433/−0
• src/Geomancy/Swizzle/Vec/Down.hs +369/−0
• src/Geomancy/Swizzle/Vec/Extend.hs +45/−0
• src/Geomancy/Swizzle/Vec/Shuffle.hs +1161/−0
• src/Geomancy/Swizzle/Vec/Up.hs +433/−0
• src/Geomancy/Transform.hs +1/−1
• src/Geomancy/U64Vec2.hs +255/−0
• src/Geomancy/U64Vec3.hs +349/−0
• src/Geomancy/U64Vec4.hs +288/−0
• src/Geomancy/UVec2.hs +11/−0
• src/Geomancy/UVec3.hs +25/−0
• src/Geomancy/UVec4.hs +19/−0
• src/Geomancy/Vec2.hs +10/−0
• src/Geomancy/Vec3.hs +25/−0
• src/Geomancy/Vec4.hs +86/−31
• src/Geomancy/Vulkan/DProjection.hs +51/−0
• src/Geomancy/Vulkan/DView.hs +38/−0
• test/Main.hs +38/−0
• test/Spec.hs +0/−430
• test/Spec/DMat4.hs +159/−0
• test/Spec/DQuaternion.hs +180/−0
• test/Spec/DTransform.hs +138/−0
• test/Spec/DVec.hs +86/−0
• test/Spec/Fusion.hs +258/−0
• test/Spec/Mat4.hs +88/−0
• test/Spec/Point.hs +30/−0
• test/Spec/Quaternion.hs +143/−0
• test/Spec/Storable.hs +87/−0
• test/Spec/Swizzle.hs +3114/−0
• test/Spec/Transform.hs +105/−0
• test/Spec/Vec4.hs +57/−0
• test/Spec/Vec64.hs +94/−0
• test/Spec/Vulkan.hs +130/−0
• test/Support/Ball.hs +265/−0
• test/Support/Fixtures.hs +116/−0
• test/Support/Gen.hs +168/−0
• test/Support/Linear.hs +69/−0
• test/Support/Property.hs +30/−0
• test/Support/Reference/Linear.hs +20/−0
• test/Support/Reference/Mat4.hs +140/−0
• test/Support/Reference/Quaternion.hs +62/−0
• test/Support/Reference/Transform.hs +107/−0
• test/Support/Reference/Vec3.hs +31/−0
• test/Support/Reference/Vulkan.hs +58/−0
• test/Support/Swizzle.hs +127/−0

@@ -1,5 +1,62 @@ # Changelog for geomancy +## 0.4.0.0++A generally safe upgrade, but quite a bit of changes across the board.++- Added double-precision `DVec2`, `DVec3` (with `Packed`), `DVec4`, `DMat4` and `DQuaternion`.+  Includes SSE2/NEON kernels for `DMat4` products and `DMat4 !* DVec4`.+- Added `DTransform` (with `DTransform.Tree`), `Vulkan.DProjection` and `Vulkan.DView`.+  Compose in double and convert with `fromTransform`/`toTransform` to keep precision far from the origin.+- Added `DPoint2`, `DPoint3`, `DPoint3P` and `DPoint4`.+  `Point`'s `AffineSpace`, `qd`, `distance` and `lerp` now work over any scalar, not just `Float`.+- Added 64-bit integer vectors `I64Vec2`, `I64Vec3` (with `Packed`), `I64Vec4`,+  `U64Vec2`, `U64Vec3` (with `Packed`), `U64Vec4`.+- Added `fromVecN`/`toVecN`, `fromMat4`/`toMat4` and `fromQuaternion`/`toQuaternion`+  for single-precision conversions.+- Added `HasField` instances for the lanes of every vector type, for `v.x`, `v.y`, `v.z`, `v.w` with `OverloadedRecordDot`.+  * `Packed` vectors have them too, and a `Point` forwards each lane to the vector it wraps.+  * Optimized code is the same as with `withVecN`: one unpack, and only the lanes used get read.+- Added rewrite rules that fuse `Vec4` and `DVec4` producers into their consumers.+  * Intermediate vectors in expressions like `dot (a + b) c`, `lerp t a b` or `(a + b) * c - d` are no longer allocated.+  * Folds over such expressions got 20-30% faster.+  * When the rules fire, only the lanes the consumer uses are evaluated, so a fused expression may avoid forcing an undefined lane.+- Added `HasField` swizzles for every vector family, in `Geomancy.Swizzle.<Family>.Shuffle`, `.Down` and `.Up`+  (families: `Vec`, `DVec`, `IVec`, `UVec`, `I64Vec`, `U64Vec`).+  * `Shuffle` keeps the size (`v.zyx`, and the identity `v.xyz`), `Down` narrows (`v.xy` of a `Vec4`), `Up` widens (`v.xyxy` of a `Vec2`).+  * The instances are orphans. `import Geomancy.Swizzle ()` brings all of them,+    or import just the modules you need, like `import Geomancy.Swizzle.Vec.Down ()`.+  * Chains like `v.wzyx.zyx.yx.y` compile to a single lane read.+- Added `Extend` swizzles that widen a vector with constant lanes, in `Geomancy.Swizzle.<Family>.Extend`.+  * The source lanes stay in order and get `0` or `1` appended: `flatPos.xy01`, `pos.xyz1`, `uv.xy0`.+  * Masks compose from extending and narrowing: `white.xyz0.wyz` is `(0, y, z)`.+- Changed `Geomancy.Swizzle` labels to use the `HasField` swizzles, so `v .@ #xzy` works for every family.++- Modified Quaternion rotation to use a unit shortcut.+  * The contract was already "use normalized" quaternions for rotations.+  * Float got a modest speedup, Double is twice-ish as fast. Both are faster than Linear.++- Added `Eq` for `Mat4`, `DMat4`, `Transform` and `DTransform`.+  Elements are compared pointwise with IEEE semantics, like the derived `Eq` of `Vec3`:+  `0 == -0` and a matrix containing `NaN` is not equal to itself.+- Changed `Eq` and `Ord` for `Vec4` from bytewise to pointwise, matching `Vec3`.+  * `Ord` is now lexicographic over x, y, z, w. Previously it compared the raw bytes, which is not a numeric order.+  * `Map`/`Set` keyed by `Vec4` will order differently, and keys differing only in the sign of zero now collapse.+    + Keys containing `NaN` are no longer found by lookups, since `Ord` is no longer total for them.+    + Basically, don't use floats for keys.++- Fixed `abs` for `IVec3` and `IVec4`, which returned the vector unchanged.+- Fixed `Block` layout for `IVec3.Packed`, whose `alignment140` and `sizeOf140` were swapped in 0.3.1.0 (now 16 and 12).+- Fixed `peek` for `Vec4` and `Mat4` returning the wrong element under optimization.+  * `Data.Vector.Storable` indexing runs `peek` with a shared `realWorld#` token, which let GHC merge the+    buffer allocations of neighbouring reads, so every read in a strict context got the last element.+  * `peek` now loads on the caller's token and builds the value in a private buffer.+  * `Vec4` `peek` and `poke` are about twice as fast, since the 16-byte copies they replace were done bytewise.++- Changed `Vec4` and `Mat4` storage from aligned pinned to regular unpinned byte arrays, like the new `DVec4` and `DMat4`.+  * SIMD kernels now use unaligned loads and stores.+  * Code reaching into the exported constructors must not take `byteArrayContents#` or pass the array to `safe` foreign calls, since the GC may now move it.+ ## 0.3.1.0  - Fixed Block.sizeOf40 for Packed 3-vectors (was 16, now 12).

@@ -1,4 +1,4 @@-Copyright Alexander Bondarenko (c) 2020+Copyright IC Rainbow (c) 2020-2026  All rights reserved. 

@@ -8,56 +8,54 @@  ## The Numbers -Storing a list of 1000 transformations (e.g. rendering instance data):+Storing a list of 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)+    poke+      10+        geomancy: OK+          42.3 ns ± 1.4 ns, 0.47x+        linear:   OK+          90.5 ns ± 3.1 ns+        linear/T: OK+          81.6 ns ± 2.7 ns+...+      1000+        geomancy: OK+          2.85 μs ±  37 ns, 0.36x+        linear:   OK+          7.93 μs ± 132 ns+        linear/T: OK+          7.08 μs ± 348 ns ``` -If you're willing to adjust your shaders, it's only 2.4 times slower.+Folding down a `gloss`-style scene tree is where it is all started:  ```-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)+    multiply+      2+        geomancy: OK+          12.3 ns ± 164 ps, 0.34x+        linear:   OK+          36.1 ns ± 1.3 ns+...+      1000+        geomancy: OK+          7.71 μs ±  93 ns, 0.23x+        linear:   OK+          34.0 μs ± 1.3 μs ``` -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)-```+Add that time to the poking that'll follow. -Folding down a `gloss`-style scene graph is where it is all started:+Sure, it is in the lower microseconds range, but this budget can be used elsewhere. -```-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 and comparing -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)+```sh+stack bench --ba '--csv baseline.csv --stdev 2 --timeout 120'+stack bench --ba '--baseline baseline.csv --fail-if-slower 5' ```--Add that time to the poking that'll follow.--Sure, it is in the lower microseconds range, but this budget can be used elsewhere.  ## Conventions 

@@ -0,0 +1,15 @@+module Bench.Common+  ( clones+  , chain+  ) where++import qualified Data.List as List++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++chain :: (a -> a) -> a -> Int -> a+chain f x0 n = List.foldl' (\x _ -> f x) x0 [1 .. n]

@@ -0,0 +1,69 @@+{-# LANGUAGE BlockArguments #-}++module Bench.DMat4 (benchmarks) where++import Test.Tasty.Bench++import Data.List (foldl1')+import qualified Foreign++import qualified Geomancy+import qualified Geomancy.DMat4++import qualified Linear+import qualified Linear.Matrix++import Bench.Common (clones)++benchmarks :: Benchmark+benchmarks = bgroup "dmat4"+  [ bgroup "transpose" dmat4transpose++  , bgroup "multiply" $ flip map [2, 3, 10, 100, 1000] \size ->+      env (clones size mempty Linear.identity) $+        dmat4multiply size++  , bgroup "poke" $ flip map [10, 100, 1000, 10000] \size ->+      env (clones size mempty Linear.identity) $+        dmat4poke size+  ]++dmat4transpose :: [Benchmark]+dmat4transpose =+  [ bench "geomancy" $ whnf Geomancy.DMat4.transpose mempty+  , bench "linear"   $ whnf Linear.Matrix.transpose (Linear.identity :: Linear.M44 Double)+  ]++dmat4multiply+  :: Int+  -> ([Geomancy.DMat4], [Linear.M44 Double])+  -> Benchmark+dmat4multiply size ~(dmat4s, m44s) = bgroup (show size)+  [ bench "geomancy" $ whnf multGeomancy dmat4s+  , bench "linear"   $ whnf multLinear m44s+  ]+  where+    multGeomancy = foldl1' (<>)+    multLinear = foldl1' (Linear.!*!)++dmat4poke+  :: Int+  -> ([Geomancy.DMat4], [Linear.M44 Double])+  -> Benchmark+dmat4poke size ~(dmat4s, m44s) = bgroup (show size)+  [ bench "geomancy" $ nfIO pokeGeomancy+  , bench "linear"   $ nfIO pokeLinear+  , bench "linear/T" $ nfIO pokeLinearTranspose+  ]+  where+    pokeGeomancy =+      Foreign.withArray dmat4s \_ptr ->+        pure ()++    pokeLinear =+      Foreign.withArray m44s \_ptr ->+        pure ()++    pokeLinearTranspose =+      Foreign.withArray (map Linear.transpose m44s) \_ptr ->+        pure ()

@@ -0,0 +1,70 @@+{-# LANGUAGE BlockArguments #-}++module Bench.DQuaternion (benchmarks) where++import Test.Tasty.Bench++import Data.List (foldl1')++import qualified Geomancy+import qualified Geomancy.DQuaternion++import qualified Linear+import qualified Linear.Quaternion++import Bench.Common (chain, clones)++benchmarks :: Benchmark+benchmarks = bgroup "dquaternion"+  [ bgroup "euler" dquatFromEuler+  , bgroup "axisAngle" dquatAxisAngle+  , bgroup "multiply" $ flip map [10, 100, 1000] \size ->+      env (clones size benchDQuat benchLQuat) $+        dquatMultiply size+  , bgroup "rotate" $ flip map [10, 100, 1000] \size ->+      env (clones size (Geomancy.dvec3 1 2 3) (Linear.V3 1 2 3)) $+        dquatRotate size+  , bgroup "rotate chain" dquatRotateChain+  ]++dquatFromEuler :: [Benchmark]+dquatFromEuler =+  [ bench "fused" $ whnf (\a -> Geomancy.DQuaternion.intrinsic a a a) (pi/3)+  , bench "composed" $ whnf (\a -> Geomancy.DQuaternion.extrinsic a a a) (pi/3)+  ]++dquatAxisAngle :: [Benchmark]+dquatAxisAngle =+  [ bench "geomancy" $ whnf (Geomancy.DQuaternion.axisAngle (Geomancy.dvec3 1 2 3)) (pi/7)+  , bench "linear"   $ whnf (Linear.Quaternion.axisAngle (Linear.V3 1 2 3 :: Linear.V3 Double)) (pi/7)+  ]++benchDQuat :: Geomancy.DQuaternion+benchDQuat = Geomancy.DQuaternion.axisAngle (Geomancy.dvec3 1 2 3) (pi/7)++benchLQuat :: Linear.Quaternion.Quaternion Double+benchLQuat = Linear.Quaternion.axisAngle (Linear.V3 1 2 3) (pi/7)++dquatMultiply+  :: Int+  -> ([Geomancy.DQuaternion], [Linear.Quaternion.Quaternion Double])+  -> Benchmark+dquatMultiply size ~(dquats, lquats) = bgroup (show size)+  [ bench "geomancy" $ whnf (foldl1' (*)) dquats+  , bench "linear"   $ whnf (foldl1' (*)) lquats+  ]++dquatRotate+  :: Int+  -> ([Geomancy.DVec3], [Linear.V3 Double])+  -> Benchmark+dquatRotate size ~(dvecs, lvecs) = bgroup (show size)+  [ bench "geomancy" $ nf (map (Geomancy.DQuaternion.rotate benchDQuat)) dvecs+  , bench "linear"   $ nf (map (Linear.Quaternion.rotate benchLQuat)) lvecs+  ]++dquatRotateChain :: [Benchmark]+dquatRotateChain =+  [ bench "geomancy" $ whnf (chain (Geomancy.DQuaternion.rotate benchDQuat) (Geomancy.dvec3 1 2 3)) 1000+  , bench "linear"   $ whnf (chain (Linear.Quaternion.rotate benchLQuat) (Linear.V3 1 2 3)) 1000+  ]

@@ -0,0 +1,69 @@+{-# LANGUAGE BlockArguments #-}++module Bench.Mat4 (benchmarks) where++import Test.Tasty.Bench++import Data.List (foldl1')+import qualified Foreign++import qualified Geomancy+import qualified Geomancy.Mat4++import qualified Linear+import qualified Linear.Matrix++import Bench.Common (clones)++benchmarks :: Benchmark+benchmarks = bgroup "mat4"+  [ bgroup "transpose" mat4transpose++  , bgroup "multiply" $ flip map [2, 3, 10, 100, 1000] \size ->+      env (clones size mempty Linear.identity) $+        mat4multiply size++  , bgroup "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.!*!)++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 ()

@@ -0,0 +1,49 @@+{-# LANGUAGE BlockArguments #-}++module Bench.Quaternion (benchmarks) where++import Test.Tasty.Bench++import qualified Geomancy+import qualified Geomancy.Quaternion++import qualified Linear+import qualified Linear.Quaternion++import Bench.Common (chain, clones)++benchmarks :: Benchmark+benchmarks = bgroup "quaternion"+  [ bgroup "euler" quatFromEuler+  , bgroup "rotate" $ flip map [10, 100, 1000] \size ->+      env (clones size (Geomancy.vec3 1 2 3) (Linear.V3 1 2 3)) $+        quatRotate size+  , bgroup "rotate chain" quatRotateChain+  ]++quatFromEuler :: [Benchmark]+quatFromEuler =+  [ bench "fused" $ whnf (\a -> Geomancy.Quaternion.intrinsic a a a) (pi/3)+  , bench "composed" $ whnf (\a -> Geomancy.Quaternion.extrinsic a a a) (pi/3)+  ]++benchQuat :: Geomancy.Quaternion+benchQuat = Geomancy.Quaternion.axisAngle (Geomancy.vec3 1 2 3) (pi/7)++benchLQuatF :: Linear.Quaternion.Quaternion Float+benchLQuatF = Linear.Quaternion.axisAngle (Linear.V3 1 2 3) (pi/7)++quatRotate+  :: Int+  -> ([Geomancy.Vec3], [Linear.V3 Float])+  -> Benchmark+quatRotate size ~(vecs, lvecs) = bgroup (show size)+  [ bench "geomancy" $ nf (map (Geomancy.Quaternion.rotate benchQuat)) vecs+  , bench "linear"   $ nf (map (Linear.Quaternion.rotate benchLQuatF)) lvecs+  ]++quatRotateChain :: [Benchmark]+quatRotateChain =+  [ bench "geomancy" $ whnf (chain (Geomancy.Quaternion.rotate benchQuat) (Geomancy.vec3 1 2 3)) 1000+  , bench "linear"   $ whnf (chain (Linear.Quaternion.rotate benchLQuatF) (Linear.V3 1 2 3)) 1000+  ]

@@ -1,86 +1,21 @@-{-# LANGUAGE BlockArguments #-}-{-# LANGUAGE TypeApplications #-}--module Main where--import Criterion.Main--import Data.List (foldl1')-import qualified Foreign+module Main (main) where -import qualified Geomancy-import qualified Geomancy.Mat4-import qualified Geomancy.Quaternion+import Test.Tasty.Bench (Benchmark, bcompare, defaultMain, locateBenchmark, mapLeafBenchmarks)+import Test.Tasty.Patterns.Printer (printAwkExpr) -import qualified Linear-import qualified Linear.Matrix+import qualified Bench.DMat4+import qualified Bench.DQuaternion+import qualified Bench.Mat4+import qualified Bench.Quaternion  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--  , bgroup "quaternion"-      [ bgroup "euler" quatFromEuler-      ]-  ]--mat4transpose :: [Benchmark]-mat4transpose =-  [ bench "geomancy" $ whnf Geomancy.Mat4.transpose mempty-  , bench "linear"   $ whnf Linear.Matrix.transpose (Linear.identity :: Linear.M44 Float)+main = defaultMain $ map (mapLeafBenchmarks vsLinear)+  [ Bench.Mat4.benchmarks+  , Bench.DMat4.benchmarks+  , Bench.Quaternion.benchmarks+  , Bench.DQuaternion.benchmarks   ] -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 ()--quatFromEuler :: [Benchmark]-quatFromEuler =-  [ bench "fused" $ whnf (\a -> Geomancy.Quaternion.intrinsic a a a) (pi/3)-  , bench "composed" $ whnf (\a -> Geomancy.Quaternion.extrinsic a a a) (pi/3)-  ]+vsLinear :: [String] -> Benchmark -> Benchmark+vsLinear ("geomancy" : path) = bcompare (printAwkExpr $ locateBenchmark ("linear" : path))+vsLinear _ = id

@@ -0,0 +1,153 @@+#if defined(__SSE2__)++#include <emmintrin.h>++void DMat4xDMat4_SIMD(double *A, double *B, double *O) {+  __m128d row0lo = _mm_loadu_pd(&A[0]);+  __m128d row0hi = _mm_loadu_pd(&A[2]);+  __m128d row1lo = _mm_loadu_pd(&A[4]);+  __m128d row1hi = _mm_loadu_pd(&A[6]);+  __m128d row2lo = _mm_loadu_pd(&A[8]);+  __m128d row2hi = _mm_loadu_pd(&A[10]);+  __m128d row3lo = _mm_loadu_pd(&A[12]);+  __m128d row3hi = _mm_loadu_pd(&A[14]);+  for(int i=0; i<4; i++) {+    __m128d brod0 = _mm_set1_pd(B[4*i + 0]);+    __m128d brod1 = _mm_set1_pd(B[4*i + 1]);+    __m128d brod2 = _mm_set1_pd(B[4*i + 2]);+    __m128d brod3 = _mm_set1_pd(B[4*i + 3]);+    __m128d rowlo = _mm_add_pd(+      _mm_add_pd(+        _mm_mul_pd(brod0, row0lo),+        _mm_mul_pd(brod1, row1lo)),+      _mm_add_pd(+        _mm_mul_pd(brod2, row2lo),+        _mm_mul_pd(brod3, row3lo)));+    __m128d rowhi = _mm_add_pd(+      _mm_add_pd(+        _mm_mul_pd(brod0, row0hi),+        _mm_mul_pd(brod1, row1hi)),+      _mm_add_pd(+        _mm_mul_pd(brod2, row2hi),+        _mm_mul_pd(brod3, row3hi)));+    _mm_storeu_pd(&O[4*i], rowlo);+    _mm_storeu_pd(&O[4*i + 2], rowhi);+  }+}++void DMat4xDVec4_SIMD(double *M, double *V, double *O) {+  __m128d row0lo = _mm_loadu_pd(&M[0]);+  __m128d row0hi = _mm_loadu_pd(&M[2]);+  __m128d row1lo = _mm_loadu_pd(&M[4]);+  __m128d row1hi = _mm_loadu_pd(&M[6]);+  __m128d row2lo = _mm_loadu_pd(&M[8]);+  __m128d row2hi = _mm_loadu_pd(&M[10]);+  __m128d row3lo = _mm_loadu_pd(&M[12]);+  __m128d row3hi = _mm_loadu_pd(&M[14]);++  __m128d vec00 = _mm_set1_pd(V[0]);+  __m128d vec01 = _mm_set1_pd(V[1]);+  __m128d vec02 = _mm_set1_pd(V[2]);+  __m128d vec03 = _mm_set1_pd(V[3]);++  __m128d reslo = _mm_mul_pd(row0lo, vec00);+  reslo = _mm_add_pd(reslo, _mm_mul_pd(row1lo, vec01));+  reslo = _mm_add_pd(reslo, _mm_mul_pd(row2lo, vec02));+  reslo = _mm_add_pd(reslo, _mm_mul_pd(row3lo, vec03));++  __m128d reshi = _mm_mul_pd(row0hi, vec00);+  reshi = _mm_add_pd(reshi, _mm_mul_pd(row1hi, vec01));+  reshi = _mm_add_pd(reshi, _mm_mul_pd(row2hi, vec02));+  reshi = _mm_add_pd(reshi, _mm_mul_pd(row3hi, vec03));++  _mm_storeu_pd(&O[0], reslo);+  _mm_storeu_pd(&O[2], reshi);+}++#elif defined(__aarch64__)++#include <arm_neon.h>++void DMat4xDMat4_SIMD(double *A, double *B, double *O) {+  float64x2_t A0lo = vld1q_f64(&A[0]);+  float64x2_t A0hi = vld1q_f64(&A[2]);+  float64x2_t A1lo = vld1q_f64(&A[4]);+  float64x2_t A1hi = vld1q_f64(&A[6]);+  float64x2_t A2lo = vld1q_f64(&A[8]);+  float64x2_t A2hi = vld1q_f64(&A[10]);+  float64x2_t A3lo = vld1q_f64(&A[12]);+  float64x2_t A3hi = vld1q_f64(&A[14]);++  for(int i=0; i<4; i++) {+    float64x2_t Blo = vld1q_f64(&B[4*i]);+    float64x2_t Bhi = vld1q_f64(&B[4*i + 2]);++    float64x2_t Clo = vmovq_n_f64(0);+    Clo = vfmaq_laneq_f64(Clo, A0lo, Blo, 0);+    Clo = vfmaq_laneq_f64(Clo, A1lo, Blo, 1);+    Clo = vfmaq_laneq_f64(Clo, A2lo, Bhi, 0);+    Clo = vfmaq_laneq_f64(Clo, A3lo, Bhi, 1);++    float64x2_t Chi = vmovq_n_f64(0);+    Chi = vfmaq_laneq_f64(Chi, A0hi, Blo, 0);+    Chi = vfmaq_laneq_f64(Chi, A1hi, Blo, 1);+    Chi = vfmaq_laneq_f64(Chi, A2hi, Bhi, 0);+    Chi = vfmaq_laneq_f64(Chi, A3hi, Bhi, 1);++    vst1q_f64(&O[4*i], Clo);+    vst1q_f64(&O[4*i + 2], Chi);+  }+}++void DMat4xDVec4_SIMD(double *M, double *V, double *O) {+  float64x2_t M0lo = vld1q_f64(&M[0]);+  float64x2_t M0hi = vld1q_f64(&M[2]);+  float64x2_t M1lo = vld1q_f64(&M[4]);+  float64x2_t M1hi = vld1q_f64(&M[6]);+  float64x2_t M2lo = vld1q_f64(&M[8]);+  float64x2_t M2hi = vld1q_f64(&M[10]);+  float64x2_t M3lo = vld1q_f64(&M[12]);+  float64x2_t M3hi = vld1q_f64(&M[14]);++  float64x2_t Vlo = vld1q_f64(&V[0]);+  float64x2_t Vhi = vld1q_f64(&V[2]);++  float64x2_t Llo = vmovq_n_f64(0);+  Llo = vfmaq_laneq_f64(Llo, M0lo, Vlo, 0);+  Llo = vfmaq_laneq_f64(Llo, M1lo, Vlo, 1);+  Llo = vfmaq_laneq_f64(Llo, M2lo, Vhi, 0);+  Llo = vfmaq_laneq_f64(Llo, M3lo, Vhi, 1);++  float64x2_t Lhi = vmovq_n_f64(0);+  Lhi = vfmaq_laneq_f64(Lhi, M0hi, Vlo, 0);+  Lhi = vfmaq_laneq_f64(Lhi, M1hi, Vlo, 1);+  Lhi = vfmaq_laneq_f64(Lhi, M2hi, Vhi, 0);+  Lhi = vfmaq_laneq_f64(Lhi, M3hi, Vhi, 1);++  vst1q_f64(&O[0], Llo);+  vst1q_f64(&O[2], Lhi);+}++#else++void DMat4xDMat4_SIMD(double *A, double *B, double *O) {+  for(int i=0; i<4; i++) {+    for(int j=0; j<4; j++) {+      O[4*i + j] =+        (B[4*i + 0] * A[0 + j] + B[4*i + 1] * A[4 + j]) ++        (B[4*i + 2] * A[8 + j] + B[4*i + 3] * A[12 + j]);+    }+  }+}++void DMat4xDVec4_SIMD(double *M, double *V, double *O) {+  for(int j=0; j<4; j++) {+    O[j] =+      M[0 + j] * V[0] ++      M[4 + j] * V[1] ++      M[8 + j] * V[2] ++      M[12 + j] * V[3];+  }+}++#endif

@@ -3,10 +3,10 @@ #include <xmmintrin.h>  void Mat4xMat4_SIMD(float *A, float *B, float *O) {-  __m128 row0 = _mm_load_ps(&A[0]);-  __m128 row1 = _mm_load_ps(&A[4]);-  __m128 row2 = _mm_load_ps(&A[8]);-  __m128 row3 = _mm_load_ps(&A[12]);+  __m128 row0 = _mm_loadu_ps(&A[0]);+  __m128 row1 = _mm_loadu_ps(&A[4]);+  __m128 row2 = _mm_loadu_ps(&A[8]);+  __m128 row3 = _mm_loadu_ps(&A[12]);   for(int i=0; i<4; i++) {     __m128 brod0 = _mm_set1_ps(B[4*i + 0]);     __m128 brod1 = _mm_set1_ps(B[4*i + 1]);@@ -19,17 +19,17 @@       _mm_add_ps(         _mm_mul_ps(brod2, row2),         _mm_mul_ps(brod3, row3)));-    _mm_store_ps(&O[4*i], row);+    _mm_storeu_ps(&O[4*i], row);   } }  void Mat4xVec4_SIMD(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 row0 = _mm_loadu_ps(&M[0]);+  __m128 row1 = _mm_loadu_ps(&M[4]);+  __m128 row2 = _mm_loadu_ps(&M[8]);+  __m128 row3 = _mm_loadu_ps(&M[12]); -  __m128 vec0 = _mm_load_ps(V);+  __m128 vec0 = _mm_loadu_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));

@@ -1,18 +1,18 @@ cabal-version: 1.12 --- This file has been generated from package.yaml by hpack version 0.39.1.+-- This file has been generated from package.yaml by hpack version 0.39.6. -- -- see: https://github.com/sol/hpack  name:           geomancy-version:        0.3.1.0+version:        0.4.0.0 synopsis:       Vectors 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:         IC Rainbow maintainer:     aenor.realm@gmail.com-copyright:      2021-2025 IC Rainbow+copyright:      2020-2026 IC Rainbow license:        BSD3 license-file:   LICENSE build-type:     Simple@@ -35,8 +35,19 @@ library   exposed-modules:       Geomancy+      Geomancy.DMat4+      Geomancy.DMat4.Internal+      Geomancy.DQuaternion+      Geomancy.DTransform+      Geomancy.DTransform.Tree+      Geomancy.DVec2+      Geomancy.DVec3+      Geomancy.DVec4       Geomancy.Elementwise       Geomancy.Gl.Funs+      Geomancy.I64Vec2+      Geomancy.I64Vec3+      Geomancy.I64Vec4       Geomancy.Interpolate       Geomancy.IVec2       Geomancy.IVec3@@ -46,9 +57,36 @@       Geomancy.Point       Geomancy.Quaternion       Geomancy.Swizzle+      Geomancy.Swizzle.DVec.Down+      Geomancy.Swizzle.DVec.Extend+      Geomancy.Swizzle.DVec.Shuffle+      Geomancy.Swizzle.DVec.Up+      Geomancy.Swizzle.I64Vec.Down+      Geomancy.Swizzle.I64Vec.Extend+      Geomancy.Swizzle.I64Vec.Shuffle+      Geomancy.Swizzle.I64Vec.Up+      Geomancy.Swizzle.IVec.Down+      Geomancy.Swizzle.IVec.Extend+      Geomancy.Swizzle.IVec.Shuffle+      Geomancy.Swizzle.IVec.Up+      Geomancy.Swizzle.U64Vec.Down+      Geomancy.Swizzle.U64Vec.Extend+      Geomancy.Swizzle.U64Vec.Shuffle+      Geomancy.Swizzle.U64Vec.Up+      Geomancy.Swizzle.UVec.Down+      Geomancy.Swizzle.UVec.Extend+      Geomancy.Swizzle.UVec.Shuffle+      Geomancy.Swizzle.UVec.Up+      Geomancy.Swizzle.Vec.Down+      Geomancy.Swizzle.Vec.Extend+      Geomancy.Swizzle.Vec.Shuffle+      Geomancy.Swizzle.Vec.Up       Geomancy.Transform       Geomancy.Transform.Tree       Geomancy.Tree+      Geomancy.U64Vec2+      Geomancy.U64Vec3+      Geomancy.U64Vec4       Geomancy.UVec2       Geomancy.UVec3       Geomancy.UVec4@@ -56,6 +94,8 @@       Geomancy.Vec3       Geomancy.Vec4       Geomancy.Vector+      Geomancy.Vulkan.DProjection+      Geomancy.Vulkan.DView       Geomancy.Vulkan.Projection       Geomancy.Vulkan.View   other-modules:@@ -65,6 +105,7 @@   ghc-options: -Wall   c-sources:       cbits/mat4.c+      cbits/dmat4.c   build-depends:       base >=4.7 && <5     , containers@@ -86,8 +127,34 @@  test-suite geomancy-test   type: exitcode-stdio-1.0-  main-is: Spec.hs+  main-is: Main.hs   other-modules:+      Spec.DMat4+      Spec.DQuaternion+      Spec.DTransform+      Spec.DVec+      Spec.Fusion+      Spec.Mat4+      Spec.Point+      Spec.Quaternion+      Spec.Storable+      Spec.Swizzle+      Spec.Transform+      Spec.Vec4+      Spec.Vec64+      Spec.Vulkan+      Support.Ball+      Support.Fixtures+      Support.Gen+      Support.Linear+      Support.Property+      Support.Reference.Linear+      Support.Reference.Mat4+      Support.Reference.Quaternion+      Support.Reference.Transform+      Support.Reference.Vec3+      Support.Reference.Vulkan+      Support.Swizzle       Paths_geomancy   hs-source-dirs:       test@@ -96,24 +163,34 @@       base >=4.7 && <5     , deepseq     , geomancy+    , gl-block     , hedgehog     , linear     , simple-affine-space+    , tasty+    , tasty-hedgehog+    , vector   default-language: Haskell2010  benchmark geomancy-bench   type: exitcode-stdio-1.0   main-is: Main.hs   other-modules:+      Bench.Common+      Bench.DMat4+      Bench.DQuaternion+      Bench.Mat4+      Bench.Quaternion       Paths_geomancy   hs-source-dirs:       bench-  ghc-options: -Wall -O2 -threaded -rtsopts "-with-rtsopts=-N1 -A128m -s"+  ghc-options: -Wall -rtsopts "-with-rtsopts=-A32m -T"   build-depends:       base >=4.7 && <5-    , criterion     , deepseq     , geomancy     , linear     , simple-affine-space+    , tasty+    , tasty-bench   default-language: Haskell2010

@@ -25,6 +25,27 @@   , pattern WithVec4   , Point4 +  -- ** Double-precision / float64s++  , DVec2+  , dvec2+  , withDVec2+  , pattern WithDVec2+  , DPoint2++  , DVec3+  , dvec3+  , withDVec3+  , pattern WithDVec3+  , DPoint3+  , DPoint3P++  , DVec4+  , dvec4+  , withDVec4+  , pattern WithDVec4+  , DPoint4+   -- ** Signed / int32s    , IVec2@@ -59,6 +80,40 @@   , withUVec4   , pattern WithUVec4 +  -- ** Signed / int64s++  , I64Vec2+  , i64vec2+  , withI64Vec2+  , pattern WithI64Vec2++  , I64Vec3+  , i64vec3+  , withI64Vec3+  , pattern WithI64Vec3++  , I64Vec4+  , i64vec4+  , withI64Vec4+  , pattern WithI64Vec4++  -- ** Unsigned / word64s++  , U64Vec2+  , u64vec2+  , withU64Vec2+  , pattern WithU64Vec2++  , U64Vec3+  , u64vec3+  , withU64Vec3+  , pattern WithU64Vec3++  , U64Vec4+  , u64vec4+  , withU64Vec4+  , pattern WithU64Vec4+   -- ** Generic maps    , Elementwise(..)@@ -66,7 +121,9 @@   -- * Matrices    , Mat4+  , DMat4   , Transform(..)+  , DTransform(..)    -- * Other beasts @@ -74,6 +131,10 @@   , quaternion   , withQuaternion +  , DQuaternion+  , dquaternion+  , withDQuaternion+   -- * Spaces    , Point(..)@@ -88,6 +149,10 @@ import Geomancy.Vec3 (Vec3, vec3, withVec3, pattern WithVec3) import Geomancy.Vec4 (Vec4, vec4, withVec4, pattern WithVec4) +import Geomancy.DVec2 (DVec2, dvec2, withDVec2, pattern WithDVec2)+import Geomancy.DVec3 (DVec3, dvec3, withDVec3, pattern WithDVec3)+import Geomancy.DVec4 (DVec4, dvec4, withDVec4, pattern WithDVec4)+ import Geomancy.IVec2 (IVec2, ivec2, withIVec2, pattern WithIVec2) import Geomancy.IVec3 (IVec3, ivec3, withIVec3, pattern WithIVec3) import Geomancy.IVec4 (IVec4, ivec4, withIVec4, pattern WithIVec4)@@ -96,12 +161,23 @@ import Geomancy.UVec3 (UVec3, uvec3, withUVec3, pattern WithUVec3) import Geomancy.UVec4 (UVec4, uvec4, withUVec4, pattern WithUVec4) +import Geomancy.I64Vec2 (I64Vec2, i64vec2, withI64Vec2, pattern WithI64Vec2)+import Geomancy.I64Vec3 (I64Vec3, i64vec3, withI64Vec3, pattern WithI64Vec3)+import Geomancy.I64Vec4 (I64Vec4, i64vec4, withI64Vec4, pattern WithI64Vec4)++import Geomancy.U64Vec2 (U64Vec2, u64vec2, withU64Vec2, pattern WithU64Vec2)+import Geomancy.U64Vec3 (U64Vec3, u64vec3, withU64Vec3, pattern WithU64Vec3)+import Geomancy.U64Vec4 (U64Vec4, u64vec4, withU64Vec4, pattern WithU64Vec4)+ import Geomancy.Elementwise (Elementwise(..))  import Geomancy.Mat4 (Mat4)+import Geomancy.DMat4 (DMat4) import Geomancy.Transform (Transform(..))+import Geomancy.DTransform (DTransform(..))  import Geomancy.Quaternion (Quaternion, quaternion, withQuaternion)+import Geomancy.DQuaternion (DQuaternion, dquaternion, withDQuaternion) -import Geomancy.Point (AffineSpace(..), Point(..), Point2, Point3, Point3P, Point4)+import Geomancy.Point (AffineSpace(..), Point(..), Point2, Point3, Point3P, Point4, DPoint2, DPoint3, DPoint3P, DPoint4) import Geomancy.Vector (VectorSpace(..), lerp, quadrance, (^*))

@@ -0,0 +1,351 @@+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE BlockArguments #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE MagicHash #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TypeApplications #-}+{-# LANGUAGE UnboxedTuples #-}+{-# LANGUAGE UnliftedFFITypes #-}+{-# LANGUAGE ViewPatterns #-}++{-# OPTIONS_GHC -Wno-orphans #-}++module Geomancy.DMat4+  ( DMat4+    -- * Operations+  , identity+  , transpose+  , inverse+  , det+  , pointwise+  , zipWith+  , matrixProduct+  , scalarMultiply+  , (!*)++    -- * Construction++    -- ** Col-major+  , colMajor+  , withColMajor+  , toListColMajor+  , toListColMajor2d+  , showColumns++    -- ** Row-major+  , rowMajor+  , withRowMajor+  , toListRowMajor+  , toListRowMajor2d+  , fromRowMajor2d+  , showRows++    -- * Precision conversion+  , fromMat4+  , toMat4+  ) where++import Prelude hiding (zipWith)+import GHC.Exts hiding (VecCount(..), toList)+import GHC.Float (double2Float, float2Double)+import Geomancy.DMat4.Internal++import System.IO.Unsafe (unsafePerformIO)+import Text.Printf (printf)++import qualified Data.Foldable as Foldable++import Geomancy.DVec4 (DVec4(..), newDVec4)+import Geomancy.Mat4.Internal (Mat4)+import qualified Geomancy.Mat4.Internal as Mat4++{-# INLINE fromMat4 #-}+fromMat4 :: Coercible a Mat4 => a -> DMat4+fromMat4 m = Mat4.withMemory (coerce m)+  \ e0 e1 e2 e3+    e4 e5 e6 e7+    e8 e9 eA eB+    eC eD eE eF ->+  fromMemory+    (float2Double e0) (float2Double e1) (float2Double e2) (float2Double e3)+    (float2Double e4) (float2Double e5) (float2Double e6) (float2Double e7)+    (float2Double e8) (float2Double e9) (float2Double eA) (float2Double eB)+    (float2Double eC) (float2Double eD) (float2Double eE) (float2Double eF)++{-# INLINE toMat4 #-}+toMat4 :: DMat4 -> Mat4+toMat4 m = withMemory m+  \ e0 e1 e2 e3+    e4 e5 e6 e7+    e8 e9 eA eB+    eC eD eE eF ->+  Mat4.fromMemory+    (double2Float e0) (double2Float e1) (double2Float e2) (double2Float e3)+    (double2Float e4) (double2Float e5) (double2Float e6) (double2Float e7)+    (double2Float e8) (double2Float e9) (double2Float eA) (double2Float eB)+    (double2Float eC) (double2Float eD) (double2Float eE) (double2Float eF)++toListColMajor :: Coercible a DMat4 => a -> [Double]+toListColMajor = toListTrans . coerce++toListColMajor2d :: Coercible a DMat4 => a -> [[Double]]+toListColMajor2d = toList2dTrans . coerce++rowMajor+  :: Coercible DMat4 a+  => Double -> Double -> Double -> Double+  -> Double -> Double -> Double -> Double+  -> Double -> Double -> Double -> Double+  -> Double -> Double -> Double -> Double+  -> a+rowMajor+  e0 e1 e2 e3+  e4 e5 e6 e7+  e8 e9 eA eB+  eC eD eE eF =+    coerce $ fromMemory+      e0 e1 e2 e3+      e4 e5 e6 e7+      e8 e9 eA eB+      eC eD eE eF++fromRowMajor2d+  :: forall t a+  .  ( Foldable t+     , Coercible DMat4 a+     )+  => t (t Double)+  -> Maybe a+fromRowMajor2d rows =+  case Foldable.toList rows of+    [r0, r1, r2, r3] ->+      withRow r0 \e0 e1 e2 e3 ->+      withRow r1 \e4 e5 e6 e7 ->+      withRow r2 \e8 e9 eA eB ->+      withRow r3 \eC eD eE eF ->+        Just . coerce @a $ rowMajor+          e0 e1 e2 e3+          e4 e5 e6 e7+          e8 e9 eA eB+          eC eD eE eF+    _ ->+      Nothing+  where+    withRow row f =+      case Foldable.toList row of+        [c0, c1, c2, c3] ->+          f c0 c1 c2 c3+        _ ->+          Nothing++{-# INLINE withRowMajor #-}+withRowMajor+  :: Coercible a DMat4+  => a+  ->+    ( Double -> Double -> Double -> Double ->+      Double -> Double -> Double -> Double ->+      Double -> Double -> Double -> Double ->+      Double -> Double -> Double -> Double ->+      r+    )+  -> r+withRowMajor m f = withMemory (coerce m)+  \ e0 e1 e2 e3+    e4 e5 e6 e7+    e8 e9 eA eB+    eC eD eE eF ->+  f+    e0 e1 e2 e3+    e4 e5 e6 e7+    e8 e9 eA eB+    eC eD eE eF++{-# INLINE colMajor #-}+colMajor+  :: forall a . Coercible DMat4 a+  => Double -> Double -> Double -> Double+  -> Double -> Double -> Double -> Double+  -> Double -> Double -> Double -> Double+  -> Double -> Double -> Double -> Double+  -> a+colMajor+  e0 e4 e8 eC+  e1 e5 e9 eD+  e2 e6 eA eE+  e3 e7 eB eF =+    coerce $ fromMemory+      e0 e1 e2 e3+      e4 e5 e6 e7+      e8 e9 eA eB+      eC eD eE eF++{-# INLINE withColMajor #-}+withColMajor+  :: Coercible a DMat4+  => a+  ->+    ( Double -> Double -> Double -> Double ->+      Double -> Double -> Double -> Double ->+      Double -> Double -> Double -> Double ->+      Double -> Double -> Double -> Double ->+      r+    )+  -> r+withColMajor m f = withMemory (coerce m)+  \ e0 e4 e8 eC+    e1 e5 e9 eD+    e2 e6 eA eE+    e3 e7 eB eF ->+  f+    e0 e1 e2 e3+    e4 e5 e6 e7+    e8 e9 eA eB+    eC eD eE eF++toListRowMajor :: Coercible a DMat4 => a -> [Double]+toListRowMajor = toListMemory . coerce++toListRowMajor2d :: Coercible a DMat4 => a -> [[Double]]+toListRowMajor2d = toList2dMemory . coerce++det :: Coercible DMat4 a => a -> Double+det m =+  withRowMajor m+    \ m00 m01 m02 m03+      m10 m11 m12 m13+      m20 m21 m22 m23+      m30 m31 m32 m33 ->+        let+          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+        in+            s0 * c5+          - s1 * c4+          + s2 * c3+          + s3 * c2+          - s4 * c1+          + s5 * c0++inverse :: Coercible DMat4 a => a -> a+inverse m =+  withRowMajor m+    \ m00 m01 m02 m03+      m10 m11 m12 m13+      m20 m21 m22 m23+      m30 m31 m32 m33 ->+        let+          invDet = recip $+              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+          rowMajor+            i00 i01 i02 i03+            i10 i11 i12 i13+            i20 i21 i22 i23+            i30 i31 i32 i33++foreign import ccall unsafe "DMat4xDMat4_SIMD" m4m4simd :: ByteArray# -> ByteArray# -> ByteArray# -> IO ()++{-# INLINE matrixProduct #-}+matrixProduct :: DMat4 -> DMat4 -> DMat4+matrixProduct (DMat4 l) (DMat4 r) = unsafePerformIO do+  result@(DMat4 m) <- newDMat4+  m4m4simd l r m+  pure result++foreign import ccall unsafe "DMat4xDVec4_SIMD" m4v4simd :: ByteArray# -> ByteArray# -> ByteArray# -> IO ()++(!*) :: Coercible a DMat4 => a -> DVec4 -> DVec4+(!*) (coerce -> DMat4 m) (DVec4 v) = unsafePerformIO do+  result@(DVec4 o) <- newDVec4+  !() <- m4v4simd m v o+  pure result++infixr 5 !*++instance Semigroup DMat4 where+  {-# INLINE (<>) #-}+  (<>) = matrixProduct++instance Monoid DMat4 where+  {-# INLINE mempty #-}+  mempty = identity++instance Show DMat4 where+  show = showColumns++showColumns :: DMat4 -> String+showColumns cm = withColMajor cm+  \ e0 e4 e8 eC+    e1 e5 e9 eD+    e2 e6 eA eE+    e3 e7 eB eF ->+  unlines+    [ printf  "/ %.4f %.4f %.4f %.4f \\" e0 e4 e8 eC+    , printf  "| %.4f %.4f %.4f %.4f |"  e1 e5 e9 eD+    , printf  "| %.4f %.4f %.4f %.4f |"  e2 e6 eA eE+    , printf "\\ %.4f %.4f %.4f %.4f /"  e3 e7 eB eF+    ]++showRows :: DMat4 -> String+showRows rm = withRowMajor rm+  \ e0 e1 e2 e3+    e4 e5 e6 e7+    e8 e9 eA eB+    eC eD eE eF ->+  unlines+    [ printf "[ %.4f %.4f %.4f %.4f |" e0 e1 e2 e3+    , printf "| %.4f %.4f %.4f %.4f |" e4 e5 e6 e7+    , printf "| %.4f %.4f %.4f %.4f |" e8 e9 eA eB+    , printf "| %.4f %.4f %.4f %.4f ]" eC eD eE eF+    ]

@@ -0,0 +1,392 @@+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE BlockArguments #-}+{-# LANGUAGE CPP #-}+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE MagicHash #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE UnboxedTuples #-}++#ifdef TH_LIFT+{-# LANGUAGE TemplateHaskell #-}+#endif++module Geomancy.DMat4.Internal+  ( DMat4(..)+  , fromMemory+  , withMemory+  , newDMat4+    -- * Order-independent operations+  , identity+  , transpose+  , pointwise+  , scalarMultiply+    -- * List operations in memory order+  , toListMemory+  , toList2dMemory+  , toListTrans+  , toList2dTrans+  , zipWith+  ) where++import Prelude hiding (zipWith)+import GHC.Exts hiding (VecCount(..), toList)++#ifdef TH_LIFT+import Language.Haskell.TH.Syntax (Lift(..))+#endif++#ifdef SERIALISE+import Codec.Serialise (Serialise(..))+import qualified Codec.Serialise.Encoding as Serialise+import qualified Codec.Serialise.Decoding as Serialise+#endif++import Control.DeepSeq (NFData(rnf))+import Foreign (Storable(..))+import Foreign.Ptr.Diff (peekDiffOff, pokeDiffOff)+import GHC.IO (IO(..))++import qualified Data.List as List++import Graphics.Gl.Block (Block(..))++data DMat4 = DMat4 ByteArray#++instance Eq DMat4 where+  DMat4 a == DMat4 b = go 0#+    where+      go i+        | isTrue# (i ==# 16#) = True+        | isTrue# (indexDoubleArray# a i ==## indexDoubleArray# b i) = go (i +# 1#)+        | otherwise = False++instance NFData DMat4 where+  rnf DMat4{} = ()++{-# INLINE fromMemory #-}+fromMemory+  :: Double -> Double -> Double -> Double+  -> Double -> Double -> Double -> Double+  -> Double -> Double -> Double -> Double+  -> Double -> Double -> Double -> Double+  -> DMat4+fromMemory+  e0 e1 e2 e3+  e4 e5 e6 e7+  e8 e9 ea eb+  ec ed ee ef =+  runRW# \world ->+    let+      !(# world_, arr #) = newByteArray# 128# world++      write i (D# e) = writeDoubleArray# arr i e++      world0 = write 0x0# e0 world_+      world1 = write 0x1# e1 world0+      world2 = write 0x2# e2 world1+      world3 = write 0x3# e3 world2++      world4 = write 0x4# e4 world3+      world5 = write 0x5# e5 world4+      world6 = write 0x6# e6 world5+      world7 = write 0x7# e7 world6++      world8 = write 0x8# e8 world7+      world9 = write 0x9# e9 world8+      worldA = write 0xA# ea world9+      worldB = write 0xB# eb worldA++      worldC = write 0xC# ec worldB+      worldD = write 0xD# ed worldC+      worldE = write 0xE# ee worldD+      worldF = write 0xF# ef worldE++      !(# _world', arr' #) = unsafeFreezeByteArray# arr worldF+    in+      DMat4 arr'++{-# INLINE withMemory #-}+withMemory+  :: DMat4+  ->+    ( Double -> Double -> Double -> Double ->+      Double -> Double -> Double -> Double ->+      Double -> Double -> Double -> Double ->+      Double -> Double -> Double -> Double ->+      r+    )+  -> r+withMemory (DMat4 arr) f =+  f+    (D# (indexDoubleArray# arr 0x0#))+    (D# (indexDoubleArray# arr 0x1#))+    (D# (indexDoubleArray# arr 0x2#))+    (D# (indexDoubleArray# arr 0x3#))++    (D# (indexDoubleArray# arr 0x4#))+    (D# (indexDoubleArray# arr 0x5#))+    (D# (indexDoubleArray# arr 0x6#))+    (D# (indexDoubleArray# arr 0x7#))++    (D# (indexDoubleArray# arr 0x8#))+    (D# (indexDoubleArray# arr 0x9#))+    (D# (indexDoubleArray# arr 0xA#))+    (D# (indexDoubleArray# arr 0xB#))++    (D# (indexDoubleArray# arr 0xC#))+    (D# (indexDoubleArray# arr 0xD#))+    (D# (indexDoubleArray# arr 0xE#))+    (D# (indexDoubleArray# arr 0xF#))++#ifdef TH_LIFT+instance Lift DMat4 where+  lift m = withMemory m+    \ e0 e1 e2 e3+      e4 e5 e6 e7+      e8 e9 eA eB+      eC eD eE eF ->+      [|+        fromMemory+          $(lift e0) $(lift e1) $(lift e2) $(lift e3)+          $(lift e4) $(lift e5) $(lift e6) $(lift e7)+          $(lift e8) $(lift e9) $(lift eA) $(lift eB)+          $(lift eC) $(lift eD) $(lift eE) $(lift eF)+      |]+  liftTyped m = withMemory m+    \ e0 e1 e2 e3+      e4 e5 e6 e7+      e8 e9 eA eB+      eC eD eE eF ->+      [||+        fromMemory+          $$(liftTyped e0) $$(liftTyped e1) $$(liftTyped e2) $$(liftTyped e3)+          $$(liftTyped e4) $$(liftTyped e5) $$(liftTyped e6) $$(liftTyped e7)+          $$(liftTyped e8) $$(liftTyped e9) $$(liftTyped eA) $$(liftTyped eB)+          $$(liftTyped eC) $$(liftTyped eD) $$(liftTyped eE) $$(liftTyped eF)+      ||]+#endif+++{-# INLINE identity #-}+identity :: DMat4+identity = fromMemory+  1 0 0 0+  0 1 0 0+  0 0 1 0+  0 0 0 1++{-# INLINE transpose #-}+transpose :: DMat4 -> DMat4+transpose =+  flip withMemory+    \ e0 e1 e2 e3+      e4 e5 e6 e7+      e8 e9 eA eB+      eC eD eE eF ->+    fromMemory+      e0 e4 e8 eC+      e1 e5 e9 eD+      e2 e6 eA eE+      e3 e7 eB eF++pointwise :: DMat4 -> DMat4 -> (Double -> Double -> Double) -> DMat4+pointwise a b f =+  withMemory a+    \ a0 a1 a2 a3+      a4 a5 a6 a7+      a8 a9 aA aB+      aC aD aE aF ->+  withMemory b+    \ b0 b1 b2 b3+      b4 b5 b6 b7+      b8 b9 bA bB+      bC bD bE bF ->+  fromMemory+    (f a0 b0) (f a1 b1) (f a2 b2) (f a3 b3)+    (f a4 b4) (f a5 b5) (f a6 b6) (f a7 b7)+    (f a8 b8) (f a9 b9) (f aA bA) (f aB bB)+    (f aC bC) (f aD bD) (f aE bE) (f aF bF)++{-# INLINE newDMat4 #-}+newDMat4 :: IO DMat4+newDMat4 =+  IO \world ->+    let+      !(# world_, arr_ #) = newByteArray# 128# world+      !(# world', arr #) = unsafeFreezeByteArray# arr_ world_+    in+      (# world', DMat4 arr #)++{-# INLINE scalarMultiply #-}+scalarMultiply :: Double -> DMat4 -> DMat4+scalarMultiply x m =+  withMemory m+    \ e0 e1 e2 e3+      e4 e5 e6 e7+      e8 e9 eA eB+      eC eD eE eF ->+      fromMemory+      (e0 * x) (e1 * x) (e2 * x) (e3 * x)+      (e4 * x) (e5 * x) (e6 * x) (e7 * x)+      (e8 * x) (e9 * x) (eA * x) (eB * x)+      (eC * x) (eD * x) (eE * x) (eF * x)++instance Storable DMat4 where+  sizeOf _mat4 = 128++  alignment _mat4 = 32++  {-# INLINE poke #-}+  poke (Ptr addr) (DMat4 arr) = IO \world ->+    let+      world' = copyByteArrayToAddr# arr 0# addr 128# world+    in+      (# world', () #)++  {-# INLINE peek #-}+  peek (Ptr addr) = IO \world ->+    let+      load i = readDoubleOffAddr# addr i++      !(# world0, e0 #) = load 0x0# world+      !(# world1, e1 #) = load 0x1# world0+      !(# world2, e2 #) = load 0x2# world1+      !(# world3, e3 #) = load 0x3# world2++      !(# world4, e4 #) = load 0x4# world3+      !(# world5, e5 #) = load 0x5# world4+      !(# world6, e6 #) = load 0x6# world5+      !(# world7, e7 #) = load 0x7# world6++      !(# world8, e8 #) = load 0x8# world7+      !(# world9, e9 #) = load 0x9# world8+      !(# worldA, eA #) = load 0xA# world9+      !(# worldB, eB #) = load 0xB# worldA++      !(# worldC, eC #) = load 0xC# worldB+      !(# worldD, eD #) = load 0xD# worldC+      !(# worldE, eE #) = load 0xE# worldD+      !(# world', eF #) = load 0xF# worldE++      !mat = fromMemory+        (D# e0) (D# e1) (D# e2) (D# e3)+        (D# e4) (D# e5) (D# e6) (D# e7)+        (D# e8) (D# e9) (D# eA) (D# eB)+        (D# eC) (D# eD) (D# eE) (D# eF)+    in+      (# world', mat #)++instance Block DMat4 where+  type PackedSize DMat4 = 128+  alignment140 _  = 32+  sizeOf140       = sizeOfPacked+  alignment430    = alignment140+  sizeOf430       = sizeOf140+  isStruct _      = False+  read140     = peekDiffOff+  write140    = pokeDiffOff+  read430     = read140+  write430    = write140+  readPacked  = read140+  writePacked = write140+  {-# INLINE alignment140 #-}+  {-# INLINE sizeOf140 #-}+  {-# INLINE alignment430 #-}+  {-# INLINE sizeOf430 #-}+  {-# INLINE isStruct #-}+  {-# INLINE read140 #-}+  {-# INLINE write140 #-}+  {-# INLINE read430 #-}+  {-# INLINE write430 #-}+  {-# INLINE readPacked #-}+  {-# INLINE writePacked #-}++#ifdef SERIALISE+instance Serialise DMat4 where+  encode m =+    withMemory m+      \ e0 e1 e2 e3+        e4 e5 e6 e7+        e8 e9 eA eB+        eC eD eE eF ->+          Serialise.encodeListLen 16 <>+          encode e0 <>+          encode e1 <> encode e2 <> encode e3 <>+          encode e4 <> encode e5 <> encode e6  <> encode e7 <>+          encode e8 <> encode e9 <> encode eA  <> encode eB <>+          encode eC <> encode eD <> encode eE  <> encode eF+  decode = do+    Serialise.decodeListLenOf 16+    e0 <- decode+    e1 <- decode+    e2 <- decode+    e3 <- decode+    e4 <- decode+    e5 <- decode+    e6 <- decode+    e7 <- decode+    e8 <- decode+    e9 <- decode+    eA <- decode+    eB <- decode+    eC <- decode+    eD <- decode+    eE <- decode+    eF <- decode+    pure $! fromMemory+      e0 e1 e2 e3+      e4 e5 e6 e7+      e8 e9 eA eB+      eC eD eE eF+#endif++toListMemory :: DMat4 -> [Double]+toListMemory = flip withMemory+    \ e0 e1 e2 e3+      e4 e5 e6 e7+      e8 e9 eA eB+      eC eD eE eF ->+    [ e0, e1, e2, e3+    , e4, e5, e6, e7+    , e8, e9, eA, eB+    , eC, eD, eE, eF+    ]++toList2dMemory :: DMat4 -> [[Double]]+toList2dMemory = flip withMemory+    \ e0 e1 e2 e3+      e4 e5 e6 e7+      e8 e9 eA eB+      eC eD eE eF ->+    [ [e0, e1, e2, e3]+    , [e4, e5, e6, e7]+    , [e8, e9, eA, eB]+    , [eC, eD, eE, eF]+    ]++toListTrans :: DMat4 -> [Double]+toListTrans = flip withMemory+    \ e0 e1 e2 e3+      e4 e5 e6 e7+      e8 e9 eA eB+      eC eD eE eF ->+    [ e0, e4, e8, eC+    , e1, e5, e9, eD+    , e2, e6, eA, eE+    , e3, e7, eB, eF+    ]++toList2dTrans :: DMat4 -> [[Double]]+toList2dTrans = flip withMemory+    \ e0 e1 e2 e3+      e4 e5 e6 e7+      e8 e9 eA eB+      eC eD eE eF ->+    [ [e0, e4, e8, eC]+    , [e1, e5, e9, eD]+    , [e2, e6, eA, eE]+    , [e3, e7, eB, eF]+    ]++zipWith :: (Double -> Double -> c) -> DMat4 -> DMat4 -> [c]+zipWith f a b = List.zipWith f (toListMemory a) (toListMemory b)

@@ -0,0 +1,358 @@+{-# LANGUAGE BlockArguments #-}+{-# LANGUAGE CPP #-}+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE TypeFamilies #-}++#ifdef TH_LIFT+{-# LANGUAGE DeriveLift #-}+{-# LANGUAGE DerivingStrategies #-}+#endif++module Geomancy.DQuaternion+  ( DQuaternion+  , dquaternion+  , withDQuaternion++  , axisAngle+  , intrinsic+  , extrinsic+  , rotate+  , rotatePoint+  , rotationBetween+  , lookAtUp++  , (^*)+  , (^/)+  , slerp++  , conjugate+  , norm+  , quadrance+  , dot+  , normalize+  , qNaN++  , fromQuaternion+  , toQuaternion+  ) where++import Control.DeepSeq (NFData(rnf))+import Foreign (Storable(..), castPtr)+import Foreign.Ptr.Diff (peekDiffOff, pokeDiffOff)+import GHC.Float (double2Float, float2Double)++#ifdef TH_LIFT+import Language.Haskell.TH.Syntax (Lift)+#endif++#ifdef SERIALISE+import Codec.Serialise (Serialise(..))+import qualified Codec.Serialise.Encoding as Serialise+import qualified Codec.Serialise.Decoding as Serialise+#endif++import Graphics.Gl.Block (Block(..))+import Geomancy.DVec3 (DVec3, dvec3, withDVec3)+import Geomancy.Quaternion (Quaternion, quaternion, withQuaternion)++import qualified Geomancy.DVec3 as DVec3++data DQuaternion = DQuaternion+  {-# UNPACK #-} !Double+  {-# UNPACK #-} !Double+  {-# UNPACK #-} !Double+  {-# UNPACK #-} !Double+  deriving (Eq, Ord, Show)+#ifdef TH_LIFT+  deriving Lift+#endif++{-# INLINE dquaternion #-}+dquaternion :: Double -> Double -> Double -> Double -> DQuaternion+dquaternion = DQuaternion++{-# INLINE withDQuaternion #-}+withDQuaternion+  :: DQuaternion+  -> (Double -> Double -> Double -> Double -> r)+  -> r+withDQuaternion (DQuaternion a b c d) f = f a b c d++{-# INLINE fromQuaternion #-}+fromQuaternion :: Quaternion -> DQuaternion+fromQuaternion q =+  withQuaternion q \a b c d ->+    dquaternion (float2Double a) (float2Double b) (float2Double c) (float2Double d)++{-# INLINE toQuaternion #-}+toQuaternion :: DQuaternion -> Quaternion+toQuaternion q =+  withDQuaternion q \a b c d ->+    quaternion (double2Float a) (double2Float b) (double2Float c) (double2Float d)++{-# INLINE (^*) #-}+(^*) :: DQuaternion -> Double -> DQuaternion+DQuaternion a b c d ^* x =+  DQuaternion+    (a * x)+    (b * x)+    (c * x)+    (d * x)++{-# INLINE (^/) #-}+(^/) :: DQuaternion -> Double -> DQuaternion+DQuaternion a b c d ^/ x =+  DQuaternion+    (a / x)+    (b / x)+    (c / x)+    (d / x)++slerp :: DQuaternion -> DQuaternion -> Double -> DQuaternion+slerp q p t+  | 1.0 - cosphi < 1e-12 =+      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 :: DQuaternion -> DQuaternion+conjugate (DQuaternion e x y z) = DQuaternion e (-x) (-y) (-z)++{-# INLINE norm #-}+norm :: DQuaternion -> Double+norm = sqrt . quadrance++{-# INLINE quadrance #-}+quadrance :: DQuaternion -> Double+quadrance q = dot q q++{-# INLINE dot #-}+dot :: DQuaternion -> DQuaternion -> Double+dot (DQuaternion a b c d) (DQuaternion e f g h) =+  a * e ++  b * f ++  c * g ++  d * h++{-# INLINE normalize #-}+normalize :: DQuaternion -> DQuaternion+normalize v =+  if nearZero q || nearZero (1-q) then+    v+  else+    let+      DQuaternion e i j k = v+    in+      DQuaternion (e / l) (i / l) (j / l) (k / l)++  where+    q = dot v v+    l = sqrt q++    nearZero a = abs a <= 1e-12++instance NFData DQuaternion where+  rnf DQuaternion{} = ()++instance Num DQuaternion where+  {-# INLINE (+) #-}+  DQuaternion a b c d + DQuaternion e f g h =+    DQuaternion+      (a + e)+      (b + f)+      (c + g)+      (d + h)++  {-# INLINE (-) #-}+  DQuaternion a b c d - DQuaternion e f g h =+    DQuaternion+      (a - e)+      (b - f)+      (c - g)+      (d - h)++  {-# INLINE (*) #-}+  DQuaternion as ax ay az * DQuaternion bs bx by bz =+    withDVec3 v \x y z ->+      DQuaternion s x y z+    where+      s = as * bs - DVec3.dot v1 v2+      v = DVec3.cross v1 v2 + v2 DVec3.^* as + v1 DVec3.^* bs+      v1 = dvec3 ax ay az+      v2 = dvec3 bx by bz++  {-# INLINE fromInteger #-}+  fromInteger x = DQuaternion (fromInteger x) 0 0 0++  {-# INLINE abs #-}+  abs z = DQuaternion (norm z) 0 0 0++  {-# INLINE signum #-}+  signum q@(DQuaternion e i j k)+    | m == 0 =+        q+    | not (isInfinite m || isNaN m) =+        DQuaternion (e * misqrt) (i * misqrt) (j * misqrt) (k * misqrt)+    | any isNaN [e, i, j, k] = qNaN+    | not (ii || ij || ik) = DQuaternion 1 0 0 0+    | not (ie || ij || ik) = DQuaternion 0 1 0 0+    | not (ie || ii || ik) = DQuaternion 0 0 1 0+    | not (ie || ii || ij) = DQuaternion 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 :: DQuaternion+qNaN = DQuaternion dNaN dNaN dNaN dNaN+  where+    dNaN = 0/0++instance Storable DQuaternion where+  {-# INLINE sizeOf #-}+  sizeOf _ = 32++  {-# INLINE alignment #-}+  alignment _ = 32++  {-# INLINE poke #-}+  poke ptr (DQuaternion 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 = DQuaternion+    <$> peek ptr'+    <*> peekElemOff ptr' 1+    <*> peekElemOff ptr' 2+    <*> peekElemOff ptr' 3+    where+      ptr' = castPtr ptr++{-# INLINE axisAngle #-}+axisAngle :: DVec3 -> Double -> DQuaternion+axisAngle axis rads =+  withDVec3 (DVec3.normalize axis DVec3.^* sin half) $+    dquaternion (cos half)+  where+    half = rads / 2++intrinsic :: Double -> Double -> Double -> DQuaternion+intrinsic roll pitch yaw =+  dquaternion+    (cx * cy * cz + sx * sy * sz)+    (sx * cy * cz + cx * sy * sz)+    (cx * sy * cz - sx * cy * sz)+    (cx * cy * sz - sx * sy * cz)+  where+    cx = cos (pitch * 0.5)+    sx = sin (pitch * 0.5)+    cy = cos (yaw * 0.5)+    sy = sin (yaw * 0.5)+    cz = cos (roll * 0.5)+    sz = sin (roll * 0.5)++extrinsic :: Double -> Double -> Double -> DQuaternion+extrinsic heading elevation tilt =+  axisAngle (dvec3 0 0 1) tilt *+  axisAngle (dvec3 1 0 0) elevation *+  axisAngle (dvec3 0 1 0) heading++{-# INLINE rotate #-}+rotate :: DQuaternion -> DVec3 -> DVec3+rotate (DQuaternion w x y z) v =+  withDVec3 v \vx vy vz ->+    let+      x2 = x + x+      y2 = y + y+      z2 = z + z+      tx = y2 * vz - z2 * vy+      ty = z2 * vx - x2 * vz+      tz = x2 * vy - y2 * vx+    in+      dvec3+        (vx + w * tx + y * tz - z * ty)+        (vy + w * ty + z * tx - x * tz)+        (vz + w * tz + x * ty - y * tx)++{-# INLINE rotatePoint #-}+rotatePoint :: DQuaternion -> DVec3 -> DVec3 -> DVec3+rotatePoint q origin point =+  origin + rotate q (point - origin)++rotationBetween :: DVec3 -> DVec3 -> DQuaternion+rotationBetween v1 v2 = axisAngle axis angle+  where+    axis = DVec3.cross v1 v2+    angle = acos cosAngle+    cosAngle =+      max (-1) . min 1 $+        DVec3.dot (DVec3.normalize v1) (DVec3.normalize v2)++lookAtUp :: DVec3 -> DVec3 -> DVec3 -> DQuaternion+lookAtUp src dst up = rot2 * rot1+  where+    dir3 = dst - src++    rot1 = rotationBetween (dvec3 0 0 1) dir3++    rot2 = rotationBetween newUp fixedUp++    newUp = rotate rot1 up+    fixedUp = DVec3.cross (DVec3.cross dir3 up) dir3++instance Block DQuaternion where+  type PackedSize DQuaternion = 32+  alignment140 _  = 32+  sizeOf140       = sizeOfPacked+  alignment430    = alignment140+  sizeOf430       = sizeOf140+  isStruct _      = True+  read140     = peekDiffOff+  write140    = pokeDiffOff+  read430     = read140+  write430    = write140+  readPacked  = read140+  writePacked = write140+  {-# INLINE alignment140 #-}+  {-# INLINE sizeOf140 #-}+  {-# INLINE alignment430 #-}+  {-# INLINE sizeOf430 #-}+  {-# INLINE isStruct #-}+  {-# INLINE read140 #-}+  {-# INLINE write140 #-}+  {-# INLINE read430 #-}+  {-# INLINE write430 #-}+  {-# INLINE readPacked #-}+  {-# INLINE writePacked #-}++#ifdef SERIALISE+instance Serialise DQuaternion where+  encode v = withDQuaternion v \x y z w ->+    Serialise.encodeListLen 4 <> encode x <> encode y <> encode z <> encode w+  decode =+    Serialise.decodeListLenOf 4 *> (dquaternion <$> decode <*> decode <*> decode <*> decode)+#endif

@@ -0,0 +1,226 @@+{-# LANGUAGE BlockArguments #-}+{-# LANGUAGE CPP #-}+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE DerivingStrategies #-}+{-# LANGUAGE GeneralisedNewtypeDeriving #-}+{-# LANGUAGE ImportQualifiedPost #-}+{-# LANGUAGE TypeFamilies #-}++#ifdef TH_LIFT+{-# LANGUAGE DeriveLift #-}+#endif++module Geomancy.DTransform+  ( DTransform(..)+  , DMat4.inverse++  , apply+  , (!.)++  , translate+  , translateV++  , rotateX+  , rotateY+  , rotateZ+  , rotateQ++  , scale+  , scaleX+  , scaleY+  , scaleZ+  , scaleXY+  , scale3+  , scaleV++  , dirPos+  , node++  , fromTransform+  , toTransform+  ) where++import Foreign (Storable(..))+import Foreign.Ptr.Diff (peekDiffOff, pokeDiffOff)++#ifdef TH_LIFT+import Language.Haskell.TH.Syntax (Lift)+#endif++#ifdef SERIALISE+import Codec.Serialise (Serialise)+#endif++import Geomancy.DMat4 (DMat4, colMajor)+import Geomancy.DQuaternion (DQuaternion, withDQuaternion)+import Geomancy.DVec3 (DVec3, dvec3, withDVec3)+import Geomancy.DVec4 (fromDVec3, withDVec4)+import Geomancy.DMat4 qualified as DMat4+import Geomancy.Transform (Transform(..))++import Graphics.Gl.Block (Block(..))++newtype DTransform = DTransform { unDTransform :: DMat4 }+  deriving newtype (Eq, Show, Semigroup, Monoid, Storable)+#ifdef TH_LIFT+  deriving Lift+#endif+#ifdef SERIALISE+  deriving Serialise+#endif++instance Block DTransform where+  type PackedSize DTransform = 128+  alignment140 _  = 32+  sizeOf140       = sizeOfPacked+  alignment430    = alignment140+  sizeOf430       = sizeOf140+  isStruct _      = False+  read140     = peekDiffOff+  write140    = pokeDiffOff+  read430     = read140+  write430    = write140+  readPacked  = read140+  writePacked = write140+  {-# INLINE alignment140 #-}+  {-# INLINE sizeOf140 #-}+  {-# INLINE alignment430 #-}+  {-# INLINE sizeOf430 #-}+  {-# INLINE isStruct #-}+  {-# INLINE read140 #-}+  {-# INLINE write140 #-}+  {-# INLINE read430 #-}+  {-# INLINE write430 #-}+  {-# INLINE readPacked #-}+  {-# INLINE writePacked #-}++{-# INLINE fromTransform #-}+fromTransform :: Transform -> DTransform+fromTransform = DTransform . DMat4.fromMat4++{-# INLINE toTransform #-}+toTransform :: DTransform -> Transform+toTransform = Transform . DMat4.toMat4 . unDTransform++apply :: DVec3 -> DTransform -> DVec3+apply = flip (!.)++(!.) :: DTransform -> DVec3 -> DVec3+(!.) mat vec =+  withDVec4 res \x y z w ->+    dvec3 (x / w) (y / w) (z / w)+  where+    res = mat DMat4.!* fromDVec3 vec 1.0++infixr 5 !.++{-# INLINE translate #-}+translate :: Double -> Double -> Double -> DTransform+translate x y z = colMajor+  1 0 0 x+  0 1 0 y+  0 0 1 z+  0 0 0 1++{-# INLINE translateV #-}+translateV :: DVec3 -> DTransform+translateV vec = withDVec3 vec translate++{-# INLINE scale3 #-}+scale3 :: Double -> Double -> Double -> DTransform+scale3 x y z = colMajor+  x 0 0 0+  0 y 0 0+  0 0 z 0+  0 0 0 1++{-# INLINE scale #-}+scale :: Double -> DTransform+scale s = scale3 s s s++{-# INLINE scaleX #-}+scaleX :: Double -> DTransform+scaleX x = scale3 x 1 1++{-# INLINE scaleY #-}+scaleY :: Double -> DTransform+scaleY y = scale3 1 y 1++{-# INLINE scaleZ #-}+scaleZ :: Double -> DTransform+scaleZ z = scale3 1 1 z++{-# INLINE scaleXY #-}+scaleXY :: Double -> Double -> DTransform+scaleXY x y = scale3 x y 1++{-# INLINE scaleV #-}+scaleV :: DVec3 -> DTransform+scaleV s = withDVec3 s scale3++{-# INLINE rotateX #-}+rotateX :: Double -> DTransform+rotateX rads =+  colMajor+    1 0   0  0+    0 c (-s) 0+    0 s   c  0+    0 0   0  1+  where+    c = cos rads+    s = sin rads++{-# INLINE rotateY #-}+rotateY :: Double -> DTransform+rotateY rads =+  colMajor+    c  0 s 0+    0  1 0 0+  (-s) 0 c 0+    0  0 0 1+  where+    c = cos rads+    s = sin rads++{-# INLINE rotateZ #-}+rotateZ :: Double -> DTransform+rotateZ rads =+  colMajor+    c (-s) 0 0+    s   c  0 0+    0   0  1 0+    0   0  0 1+  where+   c = cos rads+   s = sin rads++{-# INLINE rotateQ #-}+rotateQ :: DQuaternion -> DTransform+rotateQ dir = dirPos dir 0++{-# INLINE dirPos #-}+dirPos :: DQuaternion -> DVec3 -> DTransform+dirPos rs t =+  withDQuaternion rs \w x y z ->+  withDVec3 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++node :: DVec3 -> DQuaternion -> DVec3 -> DTransform+node t r s+ | s == 1.0 = dirPos r t+ | otherwise = dirPos r t <> scaleV s

@@ -0,0 +1,132 @@+module Geomancy.DTransform.Tree+  ( Transformed+  , apply+  , node_+  , leaf_+  , collect_++  , translate_+  , translateV_++  , rotateX_+  , rotateY_+  , rotateZ_+  , rotateQ_++  , scale_+  , scaleX_+  , scaleY_+  , scaleZ_+  , scaleXY_+  , scale3_+  ) where++import Geomancy.DQuaternion (DQuaternion)+import Geomancy.DTransform (DTransform)+import Geomancy.Tree (Tree, apply, node_, leaf_, collect_)+import Geomancy.DVec3 (DVec3)++import qualified Geomancy.DTransform as DTransform++type Transformed a = Tree (DTransform, Maybe a)++{-# INLINE translate_ #-}+translate_+  :: Double+  -> Double+  -> Double+  -> [Transformed a]+  -> Transformed a+translate_ x y z =+  node_ (DTransform.translate x y z)++{-# INLINE translateV_ #-}+translateV_+  :: DVec3+  -> [Transformed a]+  -> Transformed a+translateV_ xyz =+  node_ (DTransform.translateV xyz)++{-# INLINE rotateX_ #-}+rotateX_+  :: Double+  -> [Transformed a]+  -> Transformed a+rotateX_ rads =+  node_ (DTransform.rotateX rads)++{-# INLINE rotateY_ #-}+rotateY_+  :: Double+  -> [Transformed a]+  -> Transformed a+rotateY_ rads =+  node_ (DTransform.rotateY rads)++{-# INLINE rotateZ_ #-}+rotateZ_+  :: Double+  -> [Transformed a]+  -> Transformed a+rotateZ_ rads =+  node_ (DTransform.rotateZ rads)++{-# INLINE rotateQ_ #-}+rotateQ_+  :: DQuaternion+  -> [Transformed a]+  -> Transformed a+rotateQ_ rot =+  node_ (DTransform.rotateQ rot)++{-# INLINE scale_ #-}+scale_+  :: Double+  -> [Transformed a]+  -> Transformed a+scale_ x =+  node_ (DTransform.scale x)++{-# INLINE scaleX_ #-}+scaleX_+  :: Double+  -> [Transformed a]+  -> Transformed a+scaleX_ rads =+  node_ (DTransform.scaleX rads)++{-# INLINE scaleY_ #-}+scaleY_+  :: Double+  -> [Transformed a]+  -> Transformed a+scaleY_ rads =+  node_ (DTransform.scaleY rads)++{-# INLINE scaleZ_ #-}+scaleZ_+  :: Double+  -> [Transformed a]+  -> Transformed a+scaleZ_ rads =+  node_ (DTransform.scaleZ rads)++{-# INLINE scaleXY_ #-}+scaleXY_+  :: Double+  -> Double+  -> [Transformed a]+  -> Transformed a+scaleXY_ x y =+  node_ (DTransform.scaleXY x y)++{-# INLINE scale3_ #-}+scale3_+  :: Double+  -> Double+  -> Double+  -> [Transformed a]+  -> Transformed a+scale3_ x y z =+  node_ (DTransform.scale3 x y z)

@@ -0,0 +1,340 @@+{-# LANGUAGE BlockArguments #-}+{-# LANGUAGE CPP #-}+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE PatternSynonyms #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE ViewPatterns #-}++#ifdef TH_LIFT+{-# LANGUAGE DeriveLift #-}+{-# LANGUAGE DerivingStrategies #-}+#endif++module Geomancy.DVec2+  ( DVec2+  , dvec2+  , withDVec2+  , pattern WithDVec2+  , fromTuple+  , convert++  , (^*)+  , (^/)+  , lerp++  , dot+  , normalize++  , fromVec2+  , toVec2+  ) where++import Control.DeepSeq (NFData(rnf))+import Data.Coerce (Coercible, coerce)+import Data.MonoTraversable (Element, MonoFunctor(..), MonoPointed(..))+import Data.VectorSpace (VectorSpace)+import Foreign (Storable(..))+import Foreign.Ptr.Diff (peekDiffOff, pokeDiffOff)+import GHC.Float (double2Float, float2Double)+import GHC.Records (HasField(..))+import qualified Data.VectorSpace as VectorSpace++#ifdef TH_LIFT+import Language.Haskell.TH.Syntax (Lift)+#endif++#ifdef SERIALISE+import Codec.Serialise (Serialise(..))+import qualified Codec.Serialise.Encoding as Serialise+import qualified Codec.Serialise.Decoding as Serialise+#endif++import Geomancy.Elementwise (Elementwise(..))+import Graphics.Gl.Block (Block(..))+import Geomancy.Vec2 (Vec2)+import qualified Geomancy.Vec2 as Vec2++data DVec2 = DVec2+  {-# UNPACK #-} !Double+  {-# UNPACK #-} !Double+  deriving (Eq, Ord, Show)+#ifdef TH_LIFT+  deriving Lift+#endif++{-# INLINE dvec2 #-}+dvec2 :: Double -> Double -> DVec2+dvec2 = DVec2++{-# INLINE withDVec2 #-}+withDVec2+  :: DVec2+  -> (Double -> Double -> r)+  -> r+withDVec2 (DVec2 a b) f = f a b++{-# INLINE convert #-}+convert :: Coercible v DVec2 => (Double -> a) -> (a -> a -> r) -> v -> r+convert f t v =+  withDVec2 (coerce v) \a b ->+    t (f a) (f b)++pattern WithDVec2 :: Double -> Double -> DVec2+pattern WithDVec2 a b <- ((`withDVec2` (,)) -> (a, b))+{-# COMPLETE WithDVec2 #-}++instance HasField "x" DVec2 Double where+  {-# INLINE getField #-}+  getField v = withDVec2 v \x _ -> x++instance HasField "y" DVec2 Double where+  {-# INLINE getField #-}+  getField v = withDVec2 v \_ y -> y++{-# INLINE fromTuple #-}+fromTuple :: (Double, Double) -> DVec2+fromTuple (x, y) = dvec2 x y++{-# INLINE fromVec2 #-}+fromVec2 :: Vec2 -> DVec2+fromVec2 = Vec2.convert float2Double dvec2++{-# INLINE toVec2 #-}+toVec2 :: DVec2 -> Vec2+toVec2 = convert double2Float Vec2.vec2++instance NFData DVec2 where+  rnf DVec2{} = ()++type instance Element DVec2 = Double++instance MonoFunctor DVec2 where+  {-# INLINE omap #-}+  omap f v =+    withDVec2 v \x y ->+      dvec2 (f x) (f y)++instance MonoPointed DVec2 where+  opoint x = dvec2 x x++instance Elementwise DVec2 where+  {-# INLINE emap2 #-}+  emap2 f p0 p1 =+    withDVec2 p0 \x0 y0 ->+    withDVec2 p1 \x1 y1 ->+      dvec2+        (f x0 x1)+        (f y0 y1)++  {-# INLINE emap3 #-}+  emap3 f p0 p1 p2 =+    withDVec2 p0 \x0 y0 ->+    withDVec2 p1 \x1 y1 ->+    withDVec2 p2 \x2 y2 ->+      dvec2+        (f x0 x1 x2)+        (f y0 y1 y2)++  {-# INLINE emap4 #-}+  emap4 f p0 p1 p2 p3 =+    withDVec2 p0 \x0 y0 ->+    withDVec2 p1 \x1 y1 ->+    withDVec2 p2 \x2 y2 ->+    withDVec2 p3 \x3 y3 ->+      dvec2+        (f x0 x1 x2 x3)+        (f y0 y1 y2 y3)++  {-# INLINE emap5 #-}+  emap5 f p0 p1 p2 p3 p4 =+    withDVec2 p0 \x0 y0 ->+    withDVec2 p1 \x1 y1 ->+    withDVec2 p2 \x2 y2 ->+    withDVec2 p3 \x3 y3 ->+    withDVec2 p4 \x4 y4 ->+      dvec2+        (f x0 x1 x2 x3 x4)+        (f y0 y1 y2 y3 y4)++instance Num DVec2 where+  {-# INLINE (+) #-}+  DVec2 l1 l2 + DVec2 r1 r2 =+    DVec2+      (l1 + r1)+      (l2 + r2)++  {-# INLINE (-) #-}+  DVec2 l1 l2 - DVec2 r1 r2 =+    DVec2+      (l1 - r1)+      (l2 - r2)++  {-# INLINE (*) #-}+  DVec2 l1 l2 * DVec2 r1 r2 =+    DVec2+      (l1 * r1)+      (l2 * r2)++  {-# INLINE abs #-}+  abs (DVec2 a b) =+    DVec2 (abs a) (abs b)++  {-# INLINE signum #-}+  signum (DVec2 a b) =+    DVec2 (signum a) (signum b)++  {-# INLINE fromInteger #-}+  fromInteger x = DVec2 x' x'+    where+      x' = fromInteger x++instance Fractional DVec2 where+  {-# INLINE (/) #-}+  DVec2 l1 l2 / DVec2 r1 r2 =+    DVec2 (l1 / r1) (l2 / r2)++  {-# INLINE recip #-}+  recip (DVec2 a b) =+    DVec2 (recip a) (recip b)++  {-# INLINE fromRational #-}+  fromRational x = DVec2 x' x'+    where+      x' = fromRational x++instance Floating DVec2 where+  pi = opoint pi++  exp = omap exp+  log = omap log+  sqrt = omap sqrt+  sin = omap sin+  cos = omap cos+  asin = omap asin+  acos = omap acos+  atan = omap atan+  sinh = omap sinh+  cosh = omap cosh+  asinh = omap asinh+  acosh = omap acosh+  atanh = omap atanh++  (**) = emap2 (**)++{-# INLINE (^*) #-}+(^*) :: DVec2 -> Double -> DVec2+DVec2 a b ^* x =+  DVec2+    (a * x)+    (b * x)++{-# INLINE (^/) #-}+(^/) :: DVec2 -> Double -> DVec2+DVec2 a b ^/ x =+  DVec2+    (a / x)+    (b / x)++{-# INLINE lerp #-}+lerp :: Double -> DVec2 -> DVec2 -> DVec2+lerp alpha u v = u ^* alpha + v ^* (1 - alpha)++{-# INLINE dot #-}+dot :: DVec2 -> DVec2 -> Double+dot (DVec2 l1 l2) (DVec2 r1 r2) =+  l1 * r1 + l2 * r2++{-# INLINE normalize #-}+normalize :: DVec2 -> DVec2+normalize v =+  if nearZero q || nearZero (1 - q) then+    v+  else+    let+      DVec2 x y = v+    in+      DVec2 (x / l) (y / l)++  where+    q = dot v v+    l = sqrt q++    nearZero a = abs a <= 1e-12++instance Storable DVec2 where+  {-# INLINE sizeOf #-}+  sizeOf _ = 16++  {-# INLINE alignment #-}+  alignment _ = 16++  {-# INLINE poke #-}+  poke ptr v2 =+    withDVec2 v2 \a b -> do+      pokeByteOff ptr 0 a+      pokeByteOff ptr 8 b++  {-# INLINE peek #-}+  peek ptr = dvec2+    <$> peekByteOff ptr 0+    <*> peekByteOff ptr 8++instance Block DVec2 where+  type PackedSize DVec2 = 16+  alignment140 _  = 16+  sizeOf140 _     = 16+  alignment430    = alignment140+  sizeOf430       = sizeOf140+  isStruct _      = False+  read140     = peekDiffOff+  write140    = pokeDiffOff+  read430     = read140+  write430    = write140+  readPacked  = read140+  writePacked = write140+  {-# INLINE alignment140 #-}+  {-# INLINE sizeOf140 #-}+  {-# INLINE alignment430 #-}+  {-# INLINE sizeOf430 #-}+  {-# INLINE isStruct #-}+  {-# INLINE read140 #-}+  {-# INLINE write140 #-}+  {-# INLINE read430 #-}+  {-# INLINE write430 #-}+  {-# INLINE readPacked #-}+  {-# INLINE writePacked #-}++#ifdef SERIALISE+instance Serialise DVec2 where+  encode v = withDVec2 v \x y ->+    Serialise.encodeListLen 2 <> encode x <> encode y+  decode =+    Serialise.decodeListLenOf 2 *> (dvec2 <$> decode <*> decode)+#endif++instance VectorSpace DVec2 Double where+  zeroVector = epoint 0++  {-# INLINE (*^) #-}+  a *^ v = v Geomancy.DVec2.^* a++  {-# INLINE (^/) #-}+  v ^/ a = v Geomancy.DVec2.^/ a++  {-# INLINE (^+^) #-}+  (^+^) = emap2 (+)++  {-# INLINE (^-^) #-}+  (^-^) = emap2 (-)++  {-# INLINE negateVector #-}+  negateVector = emap negate++  {-# INLINE dot #-}+  dot = Geomancy.DVec2.dot++  {-# INLINE normalize #-}+  normalize = Geomancy.DVec2.normalize

@@ -0,0 +1,510 @@+{-# LANGUAGE BlockArguments #-}+{-# LANGUAGE CPP #-}+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE DerivingStrategies #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE GeneralizedNewtypeDeriving #-}+{-# LANGUAGE ImportQualifiedPost #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE PatternSynonyms #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TypeApplications #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE ViewPatterns #-}++#ifdef TH_LIFT+{-# LANGUAGE DeriveLift #-}+#endif++module Geomancy.DVec3+  ( DVec3+  , dvec3+  , withDVec3+  , pattern WithDVec3+  , convert+  , fromDVec2+  , fromTuple++  , (^*)+  , (^/)+  , lerp++  , cross+  , dot+  , normalize++  , Packed(..)+  , packed++  , emap2+  , emap3+  , emap4++  , fromVec3+  , toVec3+  ) where++import Control.DeepSeq (NFData(rnf))+import Data.Coerce (Coercible, coerce)+import Data.MonoTraversable (Element, MonoFunctor(..), MonoPointed(..))+import Data.VectorSpace (VectorSpace)+import Data.VectorSpace qualified as VectorSpace+import Foreign (Storable(..), castPtr)+import Foreign.Ptr.Diff (peekDiffOff, pokeDiffOff)+import GHC.Float (double2Float, float2Double)+import GHC.OverloadedLabels (IsLabel(..))+import GHC.Records (HasField(..))+import WebColor.Labels (IsWebColor(..))++#ifdef TH_LIFT+import Language.Haskell.TH.Syntax (Lift(..))+#endif++#ifdef SERIALISE+import Codec.Serialise (Serialise(..))+import qualified Codec.Serialise.Encoding as Serialise+import qualified Codec.Serialise.Decoding as Serialise+#endif++import Geomancy.DVec2 (DVec2, withDVec2)+import Geomancy.Elementwise (Elementwise(..))+import Graphics.Gl.Block (Block(..))+import Geomancy.Vec3 (Vec3)+import Geomancy.Vec3 qualified as Vec3++data DVec3 = DVec3+  {-# UNPACK #-} !Double+  {-# UNPACK #-} !Double+  {-# UNPACK #-} !Double+  deriving (Eq, Ord, Show)+#ifdef TH_LIFT+  deriving Lift+#endif++{-# INLINE dvec3 #-}+dvec3 :: Double -> Double -> Double -> DVec3+dvec3 = DVec3++{-# INLINE withDVec3 #-}+withDVec3+  :: DVec3+  -> (Double -> Double -> Double -> r)+  -> r+withDVec3 (DVec3 a b c) f = f a b c++{-# INLINE convert #-}+convert :: Coercible v DVec3 => (Double -> a) -> (a -> a -> a -> r) -> v -> r+convert f t v =+  withDVec3 (coerce v) \a b c ->+    t (f a) (f b) (f c)++pattern WithDVec3 :: Double -> Double -> Double -> DVec3+pattern WithDVec3 a b c <- ((`withDVec3` (,,)) -> (a, b, c))+{-# COMPLETE WithDVec3 #-}++instance HasField "x" DVec3 Double where+  {-# INLINE getField #-}+  getField v = withDVec3 v \x _ _ -> x++instance HasField "y" DVec3 Double where+  {-# INLINE getField #-}+  getField v = withDVec3 v \_ y _ -> y++instance HasField "z" DVec3 Double where+  {-# INLINE getField #-}+  getField v = withDVec3 v \_ _ z -> z++{-# INLINE fromDVec2 #-}+fromDVec2 :: Coercible DVec3 a => DVec2 -> Double -> a+fromDVec2 xy z =+  withDVec2 xy \x y ->+    coerce (dvec3 x y z)++{-# INLINE fromTuple #-}+fromTuple :: Coercible DVec3 a => (Double, Double, Double) -> a+fromTuple (x, y, z) = coerce (dvec3 x y z)++{-# INLINE fromVec3 #-}+fromVec3 :: Vec3 -> DVec3+fromVec3 = Vec3.convert float2Double dvec3++{-# INLINE toVec3 #-}+toVec3 :: DVec3 -> Vec3+toVec3 = convert double2Float Vec3.vec3++instance IsWebColor s => IsLabel s DVec3 where+  {-# INLINE fromLabel #-}+  fromLabel =+    webColor @s \r g b ->+      dvec3+        (fromIntegral r / 255)+        (fromIntegral g / 255)+        (fromIntegral b / 255)++instance NFData DVec3 where+  rnf DVec3{} = ()++type instance Element DVec3 = Double++instance MonoFunctor DVec3 where+  {-# INLINE omap #-}+  omap f v =+    withDVec3 v \x y z ->+      dvec3 (f x) (f y) (f z)++instance MonoPointed DVec3 where+  opoint x = dvec3 x x x++instance Elementwise DVec3 where+  {-# INLINE emap2 #-}+  emap2 f p0 p1 =+    withDVec3 p0 \x0 y0 z0 ->+    withDVec3 p1 \x1 y1 z1 ->+      dvec3+        (f x0 x1)+        (f y0 y1)+        (f z0 z1)++  {-# INLINE emap3 #-}+  emap3 f p0 p1 p2 =+    withDVec3 p0 \x0 y0 z0 ->+    withDVec3 p1 \x1 y1 z1 ->+    withDVec3 p2 \x2 y2 z2 ->+      dvec3+        (f x0 x1 x2)+        (f y0 y1 y2)+        (f z0 z1 z2)++  {-# INLINE emap4 #-}+  emap4 f p0 p1 p2 p3 =+    withDVec3 p0 \x0 y0 z0 ->+    withDVec3 p1 \x1 y1 z1 ->+    withDVec3 p2 \x2 y2 z2 ->+    withDVec3 p3 \x3 y3 z3 ->+      dvec3+        (f x0 x1 x2 x3)+        (f y0 y1 y2 y3)+        (f z0 z1 z2 z3)++  {-# INLINE emap5 #-}+  emap5 f p0 p1 p2 p3 p4 =+    withDVec3 p0 \x0 y0 z0 ->+    withDVec3 p1 \x1 y1 z1 ->+    withDVec3 p2 \x2 y2 z2 ->+    withDVec3 p3 \x3 y3 z3 ->+    withDVec3 p4 \x4 y4 z4 ->+      dvec3+        (f x0 x1 x2 x3 x4)+        (f y0 y1 y2 y3 y4)+        (f z0 z1 z2 z3 z4)++instance Num DVec3 where+  {-# INLINE (+) #-}+  DVec3 a b c + DVec3 d e f =+    DVec3+      (a + d)+      (b + e)+      (c + f)++  {-# INLINE (-) #-}+  DVec3 a b c - DVec3 d e f =+    DVec3+      (a - d)+      (b - e)+      (c - f)++  {-# INLINE (*) #-}+  DVec3 a b c * DVec3 d e f =+    DVec3+      (a * d)+      (b * e)+      (c * f)++  {-# INLINE abs #-}+  abs (DVec3 a b c) =+    DVec3 (abs a) (abs b) (abs c)++  {-# INLINE signum #-}+  signum (DVec3 a b c) =+    DVec3 (signum a) (signum b) (signum c)++  {-# INLINE fromInteger #-}+  fromInteger x = DVec3 x' x' x'+    where+      x' = fromInteger x++instance Fractional DVec3 where+  {-# INLINE (/) #-}+  DVec3 l1 l2 l3 / DVec3 r1 r2 r3 =+    DVec3 (l1 / r1) (l2 / r2) (l3 / r3)++  {-# INLINE recip #-}+  recip (DVec3 a b c) =+    DVec3 (recip a) (recip b) (recip c)++  {-# INLINE fromRational #-}+  fromRational x = DVec3 x' x' x'+    where+      x' = fromRational x++instance Floating DVec3 where+  pi = opoint pi++  exp = omap exp+  log = omap log+  sqrt = omap sqrt+  sin = omap sin+  cos = omap cos+  asin = omap asin+  acos = omap acos+  atan = omap atan+  sinh = omap sinh+  cosh = omap cosh+  asinh = omap asinh+  acosh = omap acosh+  atanh = omap atanh++  a ** b =+    withDVec3 a \ax ay az ->+    withDVec3 b \bx by bz ->+      dvec3+        (ax ** bx)+        (ay ** by)+        (az ** bz)++instance Storable DVec3 where+  {-# INLINE sizeOf #-}+  sizeOf _ = 32++  {-# INLINE alignment #-}+  alignment _ = 8++  {-# INLINE poke #-}+  poke ptr v3 =+    withDVec3 v3 \a b c -> do+      poke ptr' a+      pokeElemOff ptr' 1 b+      pokeElemOff ptr' 2 c+      pokeElemOff ptr' 3 (1.0 :: Double)+    where+      ptr' = castPtr ptr++  {-# INLINE peek #-}+  peek ptr =+    dvec3 <$> peek ptr' <*> peekElemOff ptr' 1 <*> peekElemOff ptr' 2+    where+      ptr' = castPtr ptr++instance Block DVec3 where+  type PackedSize DVec3 = 24+  alignment140 _  = 32+  sizeOf140 _     = 32+  alignment430    = alignment140+  sizeOf430       = sizeOf140+  isStruct _      = False+  read140     = peekDiffOff+  write140    = pokeDiffOff+  read430     = read140+  write430    = write140+  readPacked  = read140+  writePacked ptr diff = pokeDiffOff ptr (coerce diff) . Packed+  {-# INLINE alignment140 #-}+  {-# INLINE sizeOf140 #-}+  {-# INLINE alignment430 #-}+  {-# INLINE sizeOf430 #-}+  {-# INLINE isStruct #-}+  {-# INLINE read140 #-}+  {-# INLINE write140 #-}+  {-# INLINE read430 #-}+  {-# INLINE write430 #-}+  {-# INLINE readPacked #-}+  {-# INLINE writePacked #-}++#ifdef SERIALISE+instance Serialise DVec3 where+  encode v = withDVec3 v \x y z ->+    Serialise.encodeListLen 3 <> encode x <> encode y <> encode z+  decode =+    Serialise.decodeListLenOf 3 *> (dvec3 <$> decode <*> decode <*> decode)+#endif++{-# INLINE (^*) #-}+(^*) :: DVec3 -> Double -> DVec3+DVec3 a b c ^* x =+  DVec3+    (a * x)+    (b * x)+    (c * x)++{-# INLINE (^/) #-}+(^/) :: DVec3 -> Double -> DVec3+DVec3 a b c ^/ x =+  DVec3+    (a / x)+    (b / x)+    (c / x)++{-# INLINE lerp #-}+lerp :: Double -> DVec3 -> DVec3 -> DVec3+lerp alpha u v = u ^* alpha + v ^* (1 - alpha)++{-# INLINE cross #-}+cross :: DVec3 -> DVec3 -> DVec3+cross (DVec3 a b c) (DVec3 d e f) =+  DVec3+    (b * f - c * e)+    (c * d - a * f)+    (a * e - b * d)++{-# INLINE dot #-}+dot :: DVec3 -> DVec3 -> Double+dot (DVec3 a b c) (DVec3 d e f) =+  a * d ++  b * e ++  c * f++{-# INLINE normalize #-}+normalize :: DVec3 -> DVec3+normalize v =+  if nearZero q || nearZero (1-q) then+    v+  else+    let+      DVec3 x y z = v+    in+      DVec3 (x / l) (y / l) (z / l)++  where+    q = dot v v+    l = sqrt q++    nearZero a = abs a <= 1e-12++instance VectorSpace DVec3 Double where+  zeroVector = epoint 0++  {-# INLINE (*^) #-}+  (*^) = flip (Geomancy.DVec3.^*)++  {-# INLINE (^/) #-}+  (^/) = (Geomancy.DVec3.^/)++  {-# INLINE (^+^) #-}+  (^+^) = emap2 (+)++  {-# INLINE (^-^) #-}+  (^-^) = emap2 (-)++  {-# INLINE negateVector #-}+  negateVector = emap negate++  {-# INLINE dot #-}+  dot = Geomancy.DVec3.dot++  {-# INLINE normalize #-}+  normalize = Geomancy.DVec3.normalize++type instance Element Packed = Double++newtype Packed = Packed { unPacked :: DVec3 }+  deriving stock+    ( Eq, Ord, Show+    )+  deriving newtype+    ( NFData, Num, Fractional, Floating+    , MonoFunctor, MonoPointed+    , Elementwise+    )+#ifdef TH_LIFT+  deriving Lift+#endif++{-# INLINE packed #-}+packed :: Double -> Double -> Double -> Packed+packed x y z = Packed (dvec3 x y z)++instance HasField "x" Packed Double where+  {-# INLINE getField #-}+  getField (Packed v) = withDVec3 v \x _ _ -> x++instance HasField "y" Packed Double where+  {-# INLINE getField #-}+  getField (Packed v) = withDVec3 v \_ y _ -> y++instance HasField "z" Packed Double where+  {-# INLINE getField #-}+  getField (Packed v) = withDVec3 v \_ _ z -> z++instance Storable Packed where+  {-# INLINE sizeOf #-}+  sizeOf _ = 24++  {-# INLINE alignment #-}+  alignment _ = 8++  {-# INLINE poke #-}+  poke ptr (Packed v3) =+    withDVec3 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++instance Block Packed where+  type PackedSize Packed = 24+  alignment140 _  = 32+  sizeOf140 _     = 24+  alignment430    = alignment140+  sizeOf430       = sizeOf140+  isStruct _      = False+  read140     = peekDiffOff+  write140    = pokeDiffOff+  read430     = read140+  write430    = write140+  readPacked  = read140+  writePacked = write140+  {-# INLINE alignment140 #-}+  {-# INLINE sizeOf140 #-}+  {-# INLINE alignment430 #-}+  {-# INLINE sizeOf430 #-}+  {-# INLINE isStruct #-}+  {-# INLINE read140 #-}+  {-# INLINE write140 #-}+  {-# INLINE read430 #-}+  {-# INLINE write430 #-}+  {-# INLINE readPacked #-}+  {-# INLINE writePacked #-}++instance VectorSpace Packed Double where+  zeroVector = epoint 0++  {-# INLINE (*^) #-}+  (*^) = flip $ coerce (Geomancy.DVec3.^*)++  {-# INLINE (^/) #-}+  (^/) = coerce (Geomancy.DVec3.^/)++  {-# INLINE (^+^) #-}+  v1 ^+^ v2 = v1 + v2++  {-# INLINE (^-^) #-}+  v1 ^-^ v2 = v1 - v2++  {-# INLINE negateVector #-}+  negateVector = emap negate++  {-# INLINE dot #-}+  dot = coerce Geomancy.DVec3.dot++  {-# INLINE normalize #-}+  normalize = coerce Geomancy.DVec3.normalize

@@ -0,0 +1,516 @@+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE BlockArguments #-}+{-# LANGUAGE CPP #-}+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE ImportQualifiedPost #-}+{-# LANGUAGE MagicHash #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE PatternSynonyms #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TypeApplications #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE UnboxedTuples #-}+{-# LANGUAGE ViewPatterns #-}++#ifdef TH_LIFT+{-# LANGUAGE TemplateHaskell #-}+#endif++module Geomancy.DVec4+  ( DVec4(..)+  , dvec4+  , withDVec4+  , pattern WithDVec4+  , convert+  , fromDVec2+  , fromDVec22+  , fromDVec3+  , fromTuple++  , (^*)+  , (^/)+  , lerp++  , dot+  , normalize++  , newDVec4++  , fromVec4+  , toVec4+  ) where++import GHC.Exts hiding (VecCount(..), toList)++import Control.DeepSeq (NFData(rnf))+import Data.MonoTraversable (Element, MonoFunctor(..), MonoPointed(..))+import Data.VectorSpace (VectorSpace)+import Data.VectorSpace qualified as VectorSpace+import Foreign (Storable(..))+import Foreign.Ptr.Diff (peekDiffOff, pokeDiffOff)+import GHC.Float (double2Float, float2Double)+import GHC.IO (IO(..))+import GHC.OverloadedLabels (IsLabel(..))+import GHC.Records (HasField(..))+import Text.Printf (printf)+import WebColor.Labels (IsWebColorAlpha(..))++#ifdef TH_LIFT+import Language.Haskell.TH.Syntax (Lift(..))+#endif++#ifdef SERIALISE+import Codec.Serialise (Serialise(..))+import qualified Codec.Serialise.Encoding as Serialise+import qualified Codec.Serialise.Decoding as Serialise+#endif++import Geomancy.DVec2 (DVec2, withDVec2)+import Geomancy.DVec3 (DVec3, withDVec3)+import Geomancy.Elementwise (Elementwise(..))+import Graphics.Gl.Block (Block(..))+import Geomancy.Vec4 (Vec4)+import Geomancy.Vec4 qualified as Vec4++data DVec4 = DVec4 ByteArray#++{-# INLINE [1] dvec4 #-}+dvec4 :: Double -> Double -> Double -> Double -> DVec4+dvec4 v0 v1 v2 v3 =+  runRW# \world ->+    let+      !(# world_, arr #) = newByteArray# 32# world++      write i (D# v) = writeDoubleArray# arr i v++      world0 = write 0x0# v0 world_+      world1 = write 0x1# v1 world0+      world2 = write 0x2# v2 world1+      world3 = write 0x3# v3 world2+      !(# _world', arr' #) = unsafeFreezeByteArray# arr world3+    in+      DVec4 arr'++{-# INLINE [1] withDVec4 #-}+withDVec4+  :: DVec4+  -> (Double -> Double -> Double -> Double -> r)+  -> r+withDVec4 (DVec4 arr) f =+  f+    (D# (indexDoubleArray# arr 0x0#))+    (D# (indexDoubleArray# arr 0x1#))+    (D# (indexDoubleArray# arr 0x2#))+    (D# (indexDoubleArray# arr 0x3#))++{-# RULES+"withDVec4/dvec4" forall a b c d f.+  withDVec4 (dvec4 a b c d) f = f a b c d+"withDVec4/withDVec4" forall v k f.+  withDVec4 (withDVec4 v k) f = withDVec4 v \a b c d -> withDVec4 (k a b c d) f+  #-}++{-# INLINE convert #-}+convert :: Coercible v DVec4 => (Double -> a) -> (a -> a -> a -> a -> r) -> v -> r+convert f t v =+  withDVec4 (coerce v) \a b c d->+    t (f a) (f b) (f c) (f d)++instance Eq DVec4 where+  DVec4 a == DVec4 b = go 0#+    where+      go i+        | isTrue# (i ==# 4#) = True+        | isTrue# (indexDoubleArray# a i ==## indexDoubleArray# b i) = go (i +# 1#)+        | otherwise = False++instance Ord DVec4 where+  compare = lexicographic LT GT compare+  a < b = lexicographic True False (<) a b+  a <= b = not (b < a)+  a > b = b < a+  a >= b = not (a < b)++{-# INLINE lexicographic #-}+lexicographic :: r -> r -> (Double -> Double -> r) -> DVec4 -> DVec4 -> r+lexicographic lt gt lastElement (DVec4 a) (DVec4 b) = go 0#+  where+    go 3# = lastElement (D# (indexDoubleArray# a 3#)) (D# (indexDoubleArray# b 3#))+    go i =+      case compare (D# (indexDoubleArray# a i)) (D# (indexDoubleArray# b i)) of+        LT -> lt+        EQ -> go (i +# 1#)+        GT -> gt++instance Show DVec4 where+  show v =+    withDVec4 v $+      printf "DVec4 %.4f %.4f %.4f %.4f"++pattern WithDVec4 :: Double -> Double -> Double -> Double -> DVec4+pattern WithDVec4 a b c d <- ((`withDVec4` (,,,)) -> (a, b, c, d))+{-# COMPLETE WithDVec4 #-}++instance HasField "x" DVec4 Double where+  {-# INLINE getField #-}+  getField v = withDVec4 v \x _ _ _ -> x++instance HasField "y" DVec4 Double where+  {-# INLINE getField #-}+  getField v = withDVec4 v \_ y _ _ -> y++instance HasField "z" DVec4 Double where+  {-# INLINE getField #-}+  getField v = withDVec4 v \_ _ z _ -> z++instance HasField "w" DVec4 Double where+  {-# INLINE getField #-}+  getField v = withDVec4 v \_ _ _ w -> w++{-# INLINE fromDVec2 #-}+fromDVec2 :: DVec2 -> Double -> Double -> DVec4+fromDVec2 xy z w =+  withDVec2 xy \x y ->+    dvec4 x y z w++{-# INLINE fromDVec22 #-}+fromDVec22 :: DVec2 -> DVec2 -> DVec4+fromDVec22 xy zw =+  withDVec2 xy \x y ->+  withDVec2 zw \z w ->+    dvec4 x y z w++{-# INLINE fromDVec3 #-}+fromDVec3 :: Coercible a DVec3 => a -> Double -> DVec4+fromDVec3 xyz w =+  withDVec3 (coerce xyz) \x y z ->+    dvec4 x y z w++{-# INLINE fromTuple #-}+fromTuple :: (Double, Double, Double, Double) -> DVec4+fromTuple (x, y, z, w) = dvec4 x y z w++{-# INLINE fromVec4 #-}+fromVec4 :: Vec4 -> DVec4+fromVec4 = Vec4.convert float2Double dvec4++{-# INLINE toVec4 #-}+toVec4 :: DVec4 -> Vec4+toVec4 = convert double2Float Vec4.vec4++instance IsWebColorAlpha s => IsLabel s DVec4 where+  {-# INLINE fromLabel #-}+  fromLabel =+    webColorAlpha @s \r g b a ->+      dvec4+        (fromIntegral r / 255)+        (fromIntegral g / 255)+        (fromIntegral b / 255)+        (fromIntegral a / 255)++instance NFData DVec4 where+  rnf DVec4{} = ()++#ifdef TH_LIFT+instance Lift DVec4 where+  lift (WithDVec4 a b c d) = [| dvec4 $(lift a) $(lift b) $(lift c) $(lift d) |]+  liftTyped (WithDVec4 a b c d) = [|| dvec4 $$(liftTyped a) $$(liftTyped b) $$(liftTyped c) $$(liftTyped d) ||]+#endif++type instance Element DVec4 = Double++instance MonoFunctor DVec4 where+  {-# INLINE omap #-}+  omap f v =+    withDVec4 v \x y z w ->+      dvec4 (f x) (f y) (f z) (f w)++instance MonoPointed DVec4 where+  {-# INLINE opoint #-}+  opoint x = dvec4 x x x x++instance Elementwise DVec4 where+  {-# INLINE emap2 #-}+  emap2 f p0 p1 =+    withDVec4 p0 \x0 y0 z0 w0 ->+    withDVec4 p1 \x1 y1 z1 w1 ->+      dvec4+        (f x0 x1)+        (f y0 y1)+        (f z0 z1)+        (f w0 w1)++  {-# INLINE emap3 #-}+  emap3 f p0 p1 p2 =+    withDVec4 p0 \x0 y0 z0 w0 ->+    withDVec4 p1 \x1 y1 z1 w1 ->+    withDVec4 p2 \x2 y2 z2 w2 ->+      dvec4+        (f x0 x1 x2)+        (f y0 y1 y2)+        (f z0 z1 z2)+        (f w0 w1 w2)++  {-# INLINE emap4 #-}+  emap4 f p0 p1 p2 p3 =+    withDVec4 p0 \x0 y0 z0 w0 ->+    withDVec4 p1 \x1 y1 z1 w1 ->+    withDVec4 p2 \x2 y2 z2 w2 ->+    withDVec4 p3 \x3 y3 z3 w3 ->+      dvec4+        (f x0 x1 x2 x3)+        (f y0 y1 y2 y3)+        (f z0 z1 z2 z3)+        (f w0 w1 w2 w3)++  {-# INLINE emap5 #-}+  emap5 f p0 p1 p2 p3 p4 =+    withDVec4 p0 \x0 y0 z0 w0 ->+    withDVec4 p1 \x1 y1 z1 w1 ->+    withDVec4 p2 \x2 y2 z2 w2 ->+    withDVec4 p3 \x3 y3 z3 w3 ->+    withDVec4 p4 \x4 y4 z4 w4 ->+      dvec4+        (f x0 x1 x2 x3 x4)+        (f y0 y1 y2 y3 y4)+        (f z0 z1 z2 z3 z4)+        (f w0 w1 w2 w3 w4)++instance Num DVec4 where+  {-# INLINE (+) #-}+  (+) l r =+    withDVec4 l \l1 l2 l3 l4 ->+      withDVec4 r \r1 r2 r3 r4 ->+        dvec4+          (l1 + r1)+          (l2 + r2)+          (l3 + r3)+          (l4 + r4)++  {-# INLINE (-) #-}+  (-) l r =+    withDVec4 l \l1 l2 l3 l4 ->+      withDVec4 r \r1 r2 r3 r4 ->+        dvec4+          (l1 - r1)+          (l2 - r2)+          (l3 - r3)+          (l4 - r4)++  {-# INLINE (*) #-}+  (*) l r =+    withDVec4 l \l1 l2 l3 l4 ->+      withDVec4 r \r1 r2 r3 r4 ->+        dvec4+          (l1 * r1)+          (l2 * r2)+          (l3 * r3)+          (l4 * r4)++  {-# INLINE abs #-}+  abs v =+    withDVec4 v \a b c d ->+      dvec4 (abs a) (abs b) (abs c) (abs d)++  {-# INLINE signum #-}+  signum v =+    withDVec4 v \a b c d ->+      dvec4 (signum a) (signum b) (signum c) (signum d)++  {-# INLINE fromInteger #-}+  fromInteger x = dvec4 x' x' x' x'+    where+      x' = fromInteger x++instance Fractional DVec4 where+  {-# INLINE (/) #-}+  (/) l r =+    withDVec4 l \l1 l2 l3 l4 ->+      withDVec4 r \r1 r2 r3 r4 ->+        dvec4 (l1 / r1) (l2 / r2) (l3 / r3) (l4 / r4)++  {-# INLINE recip #-}+  recip v =+    withDVec4 v \a b c d ->+      dvec4 (recip a) (recip b) (recip c) (recip d)++  {-# INLINE fromRational #-}+  fromRational x = dvec4 x' x' x' x'+    where+      x' = fromRational x++instance Floating DVec4 where+  pi = opoint pi++  {-# INLINE exp #-}+  exp = omap exp+  {-# INLINE log #-}+  log = omap log+  {-# INLINE sqrt #-}+  sqrt = omap sqrt+  {-# INLINE sin #-}+  sin = omap sin+  {-# INLINE cos #-}+  cos = omap cos+  {-# INLINE asin #-}+  asin = omap asin+  {-# INLINE acos #-}+  acos = omap acos+  {-# INLINE atan #-}+  atan = omap atan+  {-# INLINE sinh #-}+  sinh = omap sinh+  {-# INLINE cosh #-}+  cosh = omap cosh+  {-# INLINE asinh #-}+  asinh = omap asinh+  {-# INLINE acosh #-}+  acosh = omap acosh+  {-# INLINE atanh #-}+  atanh = omap atanh++  {-# INLINE (**) #-}+  a ** b =+    withDVec4 a \ax ay az aw ->+    withDVec4 b \bx by bz bw ->+      dvec4+        (ax ** bx)+        (ay ** by)+        (az ** bz)+        (aw ** bw)++instance Storable DVec4 where+  {-# INLINE sizeOf #-}+  sizeOf _ = 32++  {-# INLINE alignment #-}+  alignment _ = 32++  {-# INLINE poke #-}+  poke (Ptr addr) (DVec4 arr) = IO \world ->+    let+      world' = copyByteArrayToAddr# arr 0# addr 32# world+    in+      (# world', () #)++  {-# INLINE peek #-}+  peek (Ptr addr) = IO \world ->+    let+      !(# world0, v0 #) = readDoubleOffAddr# addr 0x0# world+      !(# world1, v1 #) = readDoubleOffAddr# addr 0x1# world0+      !(# world2, v2 #) = readDoubleOffAddr# addr 0x2# world1+      !(# world3, v3 #) = readDoubleOffAddr# addr 0x3# world2+      !v = dvec4 (D# v0) (D# v1) (D# v2) (D# v3)+    in+      (# world3, v #)++instance Block DVec4 where+  type PackedSize DVec4 = 32+  alignment140 _  = 32+  sizeOf140       = sizeOfPacked+  alignment430    = alignment140+  sizeOf430       = sizeOf140+  isStruct _      = False+  read140     = peekDiffOff+  write140    = pokeDiffOff+  read430     = read140+  write430    = write140+  readPacked  = read140+  writePacked = write140+  {-# INLINE alignment140 #-}+  {-# INLINE sizeOf140 #-}+  {-# INLINE alignment430 #-}+  {-# INLINE sizeOf430 #-}+  {-# INLINE isStruct #-}+  {-# INLINE read140 #-}+  {-# INLINE write140 #-}+  {-# INLINE read430 #-}+  {-# INLINE write430 #-}+  {-# INLINE readPacked #-}+  {-# INLINE writePacked #-}++#ifdef SERIALISE+instance Serialise DVec4 where+  encode v = withDVec4 v \x y z w ->+    Serialise.encodeListLen 4 <> encode x <> encode y <> encode z <> encode w+  decode =+    Serialise.decodeListLenOf 4 *> (dvec4 <$> decode <*> decode <*> decode <*> decode)+#endif++instance VectorSpace DVec4 Double where+  zeroVector = epoint 0++  {-# INLINE (*^) #-}+  (*^) = flip (Geomancy.DVec4.^*)++  {-# INLINE (^/) #-}+  (^/) = (Geomancy.DVec4.^/)++  {-# INLINE (^+^) #-}+  (^+^) = emap2 (+)++  {-# INLINE (^-^) #-}+  (^-^) = emap2 (-)++  {-# INLINE negateVector #-}+  negateVector = emap negate++  {-# INLINE dot #-}+  dot = Geomancy.DVec4.dot++  {-# INLINE normalize #-}+  normalize = Geomancy.DVec4.normalize++{-# INLINE (^*) #-}+(^*) :: DVec4 -> Double -> DVec4+(^*) v x =+  withDVec4 v \a b c d ->+    dvec4+      (a * x)+      (b * x)+      (c * x)+      (d * x)++{-# INLINE (^/) #-}+(^/) :: DVec4 -> Double -> DVec4+v ^/ x = v ^* recip x++{-# INLINE lerp #-}+lerp :: Double -> DVec4 -> DVec4 -> DVec4+lerp alpha u v = u ^* alpha + v ^* (1 - alpha)++{-# INLINE dot #-}+dot :: DVec4 -> DVec4 -> Double+dot a b =+  withDVec4 a \a1 a2 a3 a4 ->+    withDVec4 b \b1 b2 b3 b4 ->+      a1 * b1 ++      a2 * b2 ++      a3 * b3 ++      a4 * b4++{-# INLINE normalize #-}+normalize :: DVec4 -> DVec4+normalize v =+  if nearZero q || nearZero (1-q) then+    v+  else+    v ^/ l++  where+    q = dot v v+    l = sqrt q++    nearZero a = abs a <= 1e-12++{-# INLINE newDVec4 #-}+newDVec4 :: IO DVec4+newDVec4 =+  IO \world ->+    let+      !(# world_, arr_ #) = newByteArray# 32# world+      !(# world', arr #) = unsafeFreezeByteArray# arr_ world_+    in+      (# world', DVec4 arr #)