linear 0.9.2 → 1.23.3
raw patch · 41 files changed
Files
- .ghci +0/−1
- .gitignore +19/−0
- .hlint.yaml +7/−0
- .travis.yml +0/−28
- CHANGELOG.markdown +417/−0
- LICENSE +1/−1
- README.markdown +2/−0
- Setup.lhs +4/−41
- linear.cabal +105/−31
- src/Linear.hs +14/−9
- src/Linear/Affine.hs +303/−0
- src/Linear/Algebra.hs +136/−0
- src/Linear/Binary.hs +27/−0
- src/Linear/Conjugate.hs +39/−2
- src/Linear/Core.hs +0/−25
- src/Linear/Covector.hs +73/−0
- src/Linear/Epsilon.hs +14/−2
- src/Linear/Instances.hs +4/−67
- src/Linear/Matrix.hs +573/−78
- src/Linear/Metric.hs +52/−5
- src/Linear/Plucker.hs +456/−42
- src/Linear/Plucker/Coincides.hs +38/−0
- src/Linear/Projection.hs +260/−0
- src/Linear/Quaternion.hs +306/−57
- src/Linear/Trace.hs +116/−0
- src/Linear/V.hs +594/−0
- src/Linear/V0.hs +244/−15
- src/Linear/V1.hs +400/−0
- src/Linear/V2.hs +306/−40
- src/Linear/V3.hs +332/−23
- src/Linear/V4.hs +463/−26
- src/Linear/Vector.hs +251/−46
- tests/Prop/Quaternion.hs +28/−0
- tests/Prop/V3.hs +8/−0
- tests/Test.hs +24/−0
- tests/Unit/Binary.hs +20/−0
- tests/Unit/Plucker.hs +36/−0
- tests/Unit/V.hs +14/−0
- tests/doctests.hs +16/−27
- travis/cabal-apt-install +0/−16
- travis/config +0/−16
− .ghci
@@ -1,1 +0,0 @@-:set -isrc -idist/build/autogen -optP-include -optPdist/build/autogen/cabal_macros.h
.gitignore view
@@ -1,4 +1,5 @@ dist+dist-newstyle docs wiki TAGS@@ -11,3 +12,21 @@ *.hi *~ *#+.stack-work/+cabal-dev+*.chi+*.chs.h+*.dyn_o+*.dyn_hi+.hpc+.hsenv+.cabal-sandbox/+cabal.sandbox.config+*.prof+*.aux+*.hp+*.eventlog+cabal.project.local+cabal.project.local~+.HTF/+.ghc.environment.*
+ .hlint.yaml view
@@ -0,0 +1,7 @@+- arguments: [-XCPP]++- ignore: {name: Use fmap}+- ignore: {name: Avoid lambda}+- ignore: {name: Redundant lambda}+- ignore: {name: Unused LANGUAGE pragma}+- ignore: {name: Eta reduce, within: [Linear.Plucker, Linear.Quaternion, Linear.V, Linear.V0, Linear.V1, Linear.V2, Linear.V3, Linear.V4]}
− .travis.yml
@@ -1,28 +0,0 @@-language: haskell-before_install:- # Uncomment whenever hackage is down.- # - mkdir -p ~/.cabal && cp travis/config ~/.cabal/config && cabal update-- # Try installing some of the build-deps with apt-get for speed.- - travis/cabal-apt-install --only-dependencies --force-reinstall $mode-- - sudo apt-get -q -y install hlint || cabal install hlint--install:- - cabal configure $mode- - cabal build--script:- - $script- - hlint src --cpp-define HLINT--notifications:- irc:- channels:- - "irc.freenode.org#haskell-lens"- skip_join: true- template:- - "\x0313linear\x03/\x0306%{branch}\x03 \x0314%{commit}\x03 %{build_url} %{message}"--env:- - mode="--enable-tests" script="cabal test"
CHANGELOG.markdown view
@@ -1,3 +1,420 @@+1.23.3 [2026.01.10]+-------------------+* Remove unused `ghc-prim`, `tagged, `transformers-compat`, and `void`+ dependencies.++1.23.2 [2025.06.17]+-------------------+* Replace `test-framework` with `tasty` in the test suite.++1.23.1 [2025.03.03]+-------------------+* Add `Uniform` and `UniformRange` instances for `Plucker`, `Quaternion`, `V`,+ and `V{0,1,2,3,4}`.++1.23 [2024.04.15]+-----------------+* The direction of interpolation of `lerp` has been reversed;+ now `lerp 0 a b == a` and `lerp 1 a b == b`.+ This brings `lerp` in line not only with its implementation+ in other languages and frameworks, but also with `slerp` in this package.++1.22 [2022.11.30]+-----------------+* The types of `_Point` and `lensP` have been generalized:++ ```diff+ -_Point :: Iso' (Point f a) (f a)+ +_Point :: Iso (Point f a) (Point g b) (f a) (g b)++ -lensP :: Lens' (Point g a) (g a)+ +lensP :: Lens (Point f a) (Point g b) (f a) (g b)+ ```++ There is a chance that existing uses of `_Point` or `lensP` will fail to+ typecheck due to their more general types. You can use `_Point.simple` or+ `lensP.simple` to restore their old, more restricted types (where `simple`+ comes from `Control.Lens` in the `lens` library).++1.21.10 [2022.06.21]+--------------------+* Allow building with `vector-0.13.*`.++1.21.9 [2022.05.18]+-------------------+* Allow building with `transformers-0.6.*`.++1.21.8 [2021.11.15]+-------------------+* Allow building with `hashable-1.4.*`.+* Drop support for pre-8.0 versions of GHC.++1.21.7 [2021.09.20]+-------------------+* Fix a build error when using `random-1.2.1` or later.++1.21.6 [2021.07.05]+-------------------+* Fix a build error when configured with `-template-haskell`.++1.21.5 [2021.02.18]+-------------------+* Allow building with `lens-5.*`.++1.21.4 [2021.01.29]+-------------------+* Allow building with `vector-0.12.2` or later.+* The build-type has been changed from `Custom` to `Simple`.+ To achieve this, the `doctests` test suite has been removed in favor of using+ [`cabal-docspec`](https://github.com/phadej/cabal-extras/tree/master/cabal-docspec)+ to run the doctests.++1.21.3 [2020.10.03]+-------------------+* Allow building with GHC 9.0.++1.21.2 [2020.09.30]+-------------------+* Use `base-orphans-0.8.3` or later. This means that the `Linear.Instances`+ module no longer defines any orphan instances of its own, and the module is+ now a simple shim on top of `Data.Orphans` from `base-orphans`.++1.21.1 [2020.06.25]+-------------------+* Allow building with `random-1.2.*`.++1.21 [2020.02.03]+-----------------+* Add instances for direct sums (`Product`) and tensor products (`Compose`) of+ other vector spaces. This makes is much more convenient to do things like treat+ a matrix temporarily as a vector through Compose, or to consider things like+ Gauss-Jordan elimination, which wants augmented structures.+* Add `frobenius` for computing the Frobenius norm of a matrix.+* Added `Random` instances for `System.Random`. We had an indirect dependency+ through `vector` anyways.+* Add "obvious" zipping `Semigroup` and `Monoid` instances to all the+ representable vector spaces.+* Add `R1`..`R4` instances to `Quaternion`. `_w` is the scalar component so that+ `_x`,`_y`,`_z` can be directional.+* Add more solvers to `Linear.Matrix`, available with `base-4.8` or later.+* Add `unangle` function to `Linear.V2`.++1.20.9 [2019.05.02]+-------------------+* Derive `Lift` instances for `Plucker`, `Quaternion`, and `V{0,1,2,3,4}`.++1.20.8 [2018.07.03]+-------------------+* Add instances of the `Field` classes from `lens`.+* Add `Epsilon` instance for `Complex`.+* Use specialized implementations of the `null` and `length` methods in+ `Foldable` instances.+* Add `Hashable1` instances for data types in `linear`. Also add a+ `Hashable` instance for `V`.+* Fix a bug in which `Quaternion`s were incorrectly exponentiated.++1.20.7+------+* Support `semigroupoids-5.2.1` and `doctest-0.12`++1.20.6+------+* Revamp `Setup.hs` to use `cabal-doctest`. This makes it build+ with `Cabal-2.0`, and makes the `doctest`s work with `cabal new-build` and+ sandboxes.+* Make `(1 / x)` and `recip x` agree in the `Fractional` instance for `Quaternion`+* Use newtype instances for `Point` vectors in `Linear.Affine`+* Enable `PolyKinds` in `Linear.Trace`. Also enable `PolyKinds` when GHC 7.6 or+ later is used (previously, it was GHC 7.8 or later).+* Fix a segfault arising from the `MVector` instance for `V`+* Add `Finite` class for conversion between `V` and fixed-size vector types++1.20.5+------+* GHC 8 compatibility+* Fixed the `perspective` calculation.++1.20.4+------+* Compatibility with `base-orphans` 0.5++1.20.3+------+* Support `vector` 0.11.0.0.+* Support `cereal` 0.5+* You can now unboxed vectors of `V n` vectors.++1.20.2+------+* Modified the `doctest` machinery to work with `stack` and builds to non-standard locations.+* Removed the local `.ghci` file.+* Various numerical stability improvements were made to the quaternion and projection functions.++1.20.1+------+* Fixed doctests broken by the previous change.+* Unboxed vector instances for various linear data types now use unpacked integers even on older GHCs.++1.20+----+* `inv22`, `inv33` and `inv44` no longer attempt an epsilon check. They no longer return a `Maybe` result as a consequence.+ You should filter for the 0 determinant case yourself.++1.19.1.3+--------+* `vector` 0.11.0.0 support++1.19.1.2+--------+* Fix GHC 7.4.++1.19.1.1+--------+* Proper `reflection` 2 support++1.19.1+------+* `reflection` 2 support++1.19+----+* Change the Ixed instance for `Linear.V` to use `Int` as the index type. This makes `V n` a _lot_ easier to use.++1.18.3+------+* Compile warning-free on GHC 7.10.+++1.18.2+------+* Added `NFData` instance for `Point`++1.18.1+------+* Added an `-f-template-haskell` option to allow disabling `template-haskell` support. This is an unsupported configuration but may be useful for expert users in sandbox configurations.+* Added lenses for extracting corner various sub-matrices e.g. `_m22`, `_m33`++1.18.0.2+--------+* Fixed builds on even older GHCs.++1.18.0.1+--------+* Fixed the test suite.+* Fixed builds on older GHCs.++1.18+----+* Consolidated `eye2` .. `eye4` into a single `identity` combinator.+* Fixed the `Data` instance `V n a` for GHC 7.10-RC3.++1.17.1.1+--------+* `filepath` 1.4 support++1.17.1+------+* Added support for `Data.Functor.Classes` from `transformers` 0.5 via `transformers-compat`.+* Added missing support for `binary`, `bytes` and `cereal` for `Point`++1.17+----+* Better support for `binary`. Added support for `bytes` and `cereal`++1.16.4+------+* `ortho` and `inverseOrtho` now only require a `Fractional` constraint.+* Added missing `Floating` instances.++1.16.3+----+* Improve the performance of `fromQuaternion`, `mkTransformation`,+ `mkTransformationMat`, `basisFor`, `scaled` by using implementations+ that inline well for functions that were previously reference+ implementations.++1.16.2+----+* Added `NFData` instances for the various vector types.+* Added `!!/` operator for matrix division by scalar.++1.16.1+----+* Added `Trace` instance for `V1`.++1.16+----+* Renamed `kronecker` to `scaled`.++1.15.5+------+* Added `Metric` instances for `[]`, `ZipList`, `Maybe`+* Added `det44` and `inv44` to `Linear.Matrix`+* Added `Data` instance for `Point`++1.15.4+------+* Added Typeable and Data instances for V++1.15.3+------+* Added missing `FunctorWithIndex`, `FoldableWithIndex` and `TraversableWithIndex Int (V n)` instances for `V`++1.15.2+------+* Added `frustum`, analogous to the old `glFrustum` call.+* Added `inverseInfinitePerspective`, `inverseOrtho`, `inverseFrustum`.++1.15.1+------+* Added `inversePerspective`. It is much more accurate to compute it directly than to compute an inverse.++1.15.0.1+--------+* Fixed build failures caused by `Linear` re-exporting the old name.++1.15+----+* Renamed `Linear.Perspective` to `Linear.Projection`.+* Fixed a build issue with GHC HEAD.++1.14.0.1+--------+* Fixed test failures caused by 1.14++1.14+----+* Moved `Coincides` to `Linear.Plucker.Coincides`. The constructors `Line` and `Ray` oft collided with user code.++1.13+----+* Switched 'ortho' to follow the OpenGL handedness.++1.12.1+------+* Added "swizzle" lenses **e.g.** `_yzx`, which are useful for working with libraries like `gl`.++1.12+------+* Added 'transpose'+* Added missing 'Mxy' matrices up to 4 dimensions -- they were commonly reimplemented by users.++1.11.3+------+* Fixed an issue with `UndecidableInstances` on GHC 7.6.3++1.11.2+------+* Added `Linear.Perspective`.++1.11.1+------+* Added `_Point`, `relative` and a few instances for `Point`.++1.11+----+* Changed the 'representation' of `V n` from `E (V n)`, which was hard to use, to `Int`, which is a bit too permissive, but is easy to use.++1.10.1+------+* Added `Linear.V2.angle`.++1.10+----+* Added `Hashable` instances.++1.9.1+-----+* Added a role annotation to `V n a` to prevent users from using GHC 7.8's `Coercible` machinery to violate invariants.++1.9.0.1+-----+* Fixed a broken build++1.9+---+* Added `MonadZip` instances.+* Added `MonadFix` instances.+* Added `Control.Lens.Each.Each` instances++1.8.1+-----+* Bugfixed `slerp`++1.8+---+* Added missing `Unbox` instances for working with unboxed vectors of `linear` data types.++1.7+---+* Fixed `axisAngle`+* `unit` now has a rank 1 type.++1.5+---+* `lens` 4 compatibility++1.4+---+* Renamed `incore` to `column` and added an example.++1.3.1.1+-------+* Build bugfix++1.3.1+---+* Better implementations of `basis` and `basisFor`.+* Derived Generic instances.++1.2+---+* Improved matrix multiplication to properly support the sparse/sparse case.++1.1.4+-----+* Marked modules `Trustworthy` as necessary.++1.1.2+-----+* Dependency bump for `reflection` compatibility++1.1.1+-----+* Fixed an infinite loop in the default definition of `liftI2`.++1.1+---+* Added `Additive` instances for `[]`, `Maybe` and `Vector`.++1.0+---+* Strict vectors+* Exported `mkTransformationMat`+* Bumped dependency bounds++0.9.1 [bug fix]+-----+* Exported `Linear.V0`!++0.9+---+* Added sparse vector support.++0.8+---+* Added `Linear.V0`++0.7+---+* Added `Linear.Instances`+* More documentation++0.6+---+* Removed the direct dependency on `lens`.+* Added `Linear.Core` to cover vector spaces as corepresentable functors.+ 0.5 ------- * Added `Ix` instances for `V2`, `V3`, and `V4`
LICENSE view
@@ -1,4 +1,4 @@-Copyright 2011-13 Edward Kmett+Copyright 2011-2015 Edward Kmett All rights reserved.
README.markdown view
@@ -1,6 +1,8 @@ linear ====== +[](https://hackage.haskell.org/package/linear) [](https://github.com/ekmett/linear/actions?query=workflow%3AHaskell-CI)+ Highly polymorphic vector space operations on sparse and free vector spaces. Contact Information
Setup.lhs view
@@ -1,44 +1,7 @@ #!/usr/bin/runhaskell-\begin{code}-{-# OPTIONS_GHC -Wall #-}-module Main (main) where--import Data.List ( nub )-import Data.Version ( showVersion )-import Distribution.Package ( PackageName(PackageName), PackageId, InstalledPackageId, packageVersion, packageName )-import Distribution.PackageDescription ( PackageDescription(), TestSuite(..) )-import Distribution.Simple ( defaultMainWithHooks, UserHooks(..), simpleUserHooks )-import Distribution.Simple.Utils ( rewriteFile, createDirectoryIfMissingVerbose )-import Distribution.Simple.BuildPaths ( autogenModulesDir )-import Distribution.Simple.Setup ( BuildFlags(buildVerbosity), fromFlag )-import Distribution.Simple.LocalBuildInfo ( withLibLBI, withTestLBI, LocalBuildInfo(), ComponentLocalBuildInfo(componentPackageDeps) )-import Distribution.Verbosity ( Verbosity )-import System.FilePath ( (</>) )--main :: IO ()-main = defaultMainWithHooks simpleUserHooks- { buildHook = \pkg lbi hooks flags -> do- generateBuildModule (fromFlag (buildVerbosity flags)) pkg lbi- buildHook simpleUserHooks pkg lbi hooks flags- }--generateBuildModule :: Verbosity -> PackageDescription -> LocalBuildInfo -> IO ()-generateBuildModule verbosity pkg lbi = do- let dir = autogenModulesDir lbi- createDirectoryIfMissingVerbose verbosity True dir- withLibLBI pkg lbi $ \_ libcfg -> do- withTestLBI pkg lbi $ \suite suitecfg -> do- rewriteFile (dir </> "Build_" ++ testName suite ++ ".hs") $ unlines- [ "module Build_" ++ testName suite ++ " where"- , "deps :: [String]"- , "deps = " ++ (show $ formatdeps (testDeps libcfg suitecfg))- ]- where- formatdeps = map (formatone . snd)- formatone p = case packageName p of- PackageName n -> n ++ "-" ++ showVersion (packageVersion p)+> module Main (main) where -testDeps :: ComponentLocalBuildInfo -> ComponentLocalBuildInfo -> [(InstalledPackageId, PackageId)]-testDeps xs ys = nub $ componentPackageDeps xs ++ componentPackageDeps ys+> import Distribution.Simple -\end{code}+> main :: IO ()+> main = defaultMain
linear.cabal view
@@ -1,75 +1,149 @@ name: linear category: Math, Algebra-version: 0.9.2+version: 1.23.3 license: BSD3-cabal-version: >= 1.8+cabal-version: >= 1.10 license-file: LICENSE author: Edward A. Kmett maintainer: Edward A. Kmett <ekmett@gmail.com> stability: provisional homepage: http://github.com/ekmett/linear/ bug-reports: http://github.com/ekmett/linear/issues-copyright: Copyright (C) 2012-2013 Edward A. Kmett+copyright: Copyright (C) 2012-2015 Edward A. Kmett synopsis: Linear Algebra description: Types and combinators for linear algebra on free vector spaces-build-type: Custom-tested-with: GHC == 7.4.1, GHC == 7.4.2, GHC == 7.6.1+build-type: Simple+tested-with: GHC == 8.0.2+ , GHC == 8.2.2+ , GHC == 8.4.4+ , GHC == 8.6.5+ , GHC == 8.8.4+ , GHC == 8.10.7+ , GHC == 9.0.2+ , GHC == 9.2.8+ , GHC == 9.4.8+ , GHC == 9.6.7+ , GHC == 9.8.4+ , GHC == 9.10.3+ , GHC == 9.12.2+ , GHC == 9.14.1 extra-source-files:- .ghci .gitignore- .travis.yml+ .hlint.yaml .vim.custom- travis/cabal-apt-install- travis/config CHANGELOG.markdown README.markdown +flag template-haskell+ description:+ You can disable the use of the `template-haskell` package using `-f-template-haskell`.+ .+ Disabling this is an unsupported configuration, but it may be useful for accelerating builds in sandboxes for expert users.+ default: True+ manual: True++flag herbie+ description: Enable `herbie`.+ default: False+ manual: True+ source-repository head type: git- location: git://github.com/ekmett/linear.git+ location: https://github.com/ekmett/linear library build-depends:- base >= 4.5 && < 5,- containers >= 0.4 && < 0.6,- distributive >= 0.2.2,- hashable >= 1.1 && < 1.3,- semigroups >= 0.9,- semigroupoids >= 3,- transformers >= 0.2 && < 0.4,- unordered-containers >= 0.2.3+ adjunctions >= 4 && < 5,+ base >= 4.9 && < 5,+ base-orphans >= 0.8.3 && < 1,+ binary >= 0.5 && < 0.9,+ bytes >= 0.15 && < 1,+ cereal >= 0.4.1.1 && < 0.6,+ containers >= 0.4 && < 0.9,+ deepseq >= 1.1 && < 1.6,+ distributive >= 0.5.1 && < 1,+ hashable >= 1.2.7.0 && < 1.6,+ indexed-traversable >= 0.1.1 && < 0.2,+ lens >= 4.15.2 && < 6,+ random >= 1.2 && < 1.4,+ reflection >= 2 && < 3,+ semigroupoids >= 5.2.1 && < 7,+ transformers >= 0.5 && < 0.7,+ unordered-containers >= 0.2.3 && < 0.3,+ vector >= 0.12.1.2 && < 0.14 + if flag(template-haskell) && impl(ghc)+ build-depends: template-haskell >= 2.11.1.0 && < 3.0++ if flag(herbie)+ build-depends: HerbiePlugin >= 0.1 && < 0.2+ ghc-options: -fplugin=Herbie+ cpp-options: -DHERBIE+ exposed-modules: Linear+ Linear.Affine+ Linear.Algebra+ Linear.Binary Linear.Conjugate- Linear.Core+ Linear.Covector Linear.Epsilon Linear.Instances Linear.Matrix Linear.Metric Linear.Plucker+ Linear.Plucker.Coincides+ Linear.Projection Linear.Quaternion+ Linear.Trace+ Linear.V Linear.V0+ Linear.V1 Linear.V2 Linear.V3 Linear.V4 Linear.Vector - ghc-options: -Wall -fwarn-tabs -O2 -fdicts-cheap -funbox-strict-fields+ ghc-options: -Wall -Wtabs -O2 -fdicts-cheap -funbox-strict-fields -Wno-trustworthy-safe hs-source-dirs: src --- Verify the results of the examples+ default-language: Haskell2010++ x-docspec-extra-packages: simple-reflect++-- We need this dummy test-suite to add simple-reflect to the install plan+--+-- When cabal-install's extra-packages support becomes widely available+-- (i.e. after 3.4 release), we can remove this test-suite. test-suite doctests- type: exitcode-stdio-1.0- main-is: doctests.hs- build-depends:- base,- directory >= 1.0 && < 1.3,- doctest >= 0.8 && < 0.10,- filepath >= 1.3 && < 1.4,- lens >= 3.7,- simple-reflect >= 0.3.1+ type: exitcode-stdio-1.0+ main-is: doctests.hs+ hs-source-dirs: tests+ default-language: Haskell2010 - ghc-options: -Wall -Werror -threaded+ build-depends: base < 5, simple-reflect >= 0.3.1++test-suite test+ type: exitcode-stdio-1.0+ main-is: Test.hs+ other-modules: Prop.Quaternion+ Prop.V3+ Unit.Binary+ Unit.Plucker+ Unit.V+ ghc-options: -Wall -threaded hs-source-dirs: tests+ build-depends:+ base,+ binary,+ bytestring,+ deepseq,+ tasty >= 1.4 && < 1.6,+ tasty-hunit >= 0.10 && < 0.11,+ tasty-quickcheck >= 0.10 && < 0.12,+ linear,+ QuickCheck >= 2.5,+ reflection,+ vector+ default-language: Haskell2010
src/Linear.hs view
@@ -1,7 +1,6 @@ ----------------------------------------------------------------------------- -- |--- Module : Linear--- Copyright : (C) 2012-2013 Edward Kmett,+-- Copyright : (C) 2012-2015 Edward Kmett -- License : BSD-style (see the file LICENSE) -- -- Maintainer : Edward Kmett <ekmett@gmail.com>@@ -12,32 +11,38 @@ -- that make up the linear package. ---------------------------------------------------------------------------- module Linear- ( module Linear.Conjugate- , module Linear.Core+ ( module Linear.Algebra+ , module Linear.Binary+ , module Linear.Conjugate+ , module Linear.Covector , module Linear.Epsilon , module Linear.Matrix , module Linear.Metric- , module Linear.Plucker+ , module Linear.Projection , module Linear.Quaternion+ , module Linear.Trace , module Linear.V0+ , module Linear.V1 , module Linear.V2 , module Linear.V3 , module Linear.V4 , module Linear.Vector ) where +import Linear.Algebra+import Linear.Binary import Linear.Conjugate-import Linear.Core+import Linear.Covector import Linear.Epsilon import Linear.Instances () import Linear.Matrix import Linear.Metric-import Linear.Plucker+import Linear.Projection import Linear.Quaternion+import Linear.Trace import Linear.V0+import Linear.V1 import Linear.V2 import Linear.V3 import Linear.V4 import Linear.Vector--{-# ANN module "Hlint: ignore Use import/export shortcut" #-}
+ src/Linear/Affine.hs view
@@ -0,0 +1,303 @@+{-# LANGUAGE CPP #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE GeneralizedNewtypeDeriving #-}+{-# LANGUAGE DeriveTraversable #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE UndecidableInstances #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE RankNTypes #-}+{-# LANGUAGE Trustworthy #-}+{-# LANGUAGE DeriveGeneric #-}+{-# LANGUAGE DeriveDataTypeable #-}+{-# LANGUAGE StandaloneDeriving #-}+{-# LANGUAGE ScopedTypeVariables #-}++-----------------------------------------------------------------------------+-- |+-- License : BSD-style (see the file LICENSE)+-- Maintainer : Edward Kmett <ekmett@gmail.com>+-- Stability : provisional+-- Portability : portable+--+-- Operations on affine spaces.+-----------------------------------------------------------------------------+module Linear.Affine where++import Control.Applicative+import Control.DeepSeq+import Control.Monad (liftM)+import Control.Lens+import Data.Binary as Binary+import Data.Bytes.Serial+import Data.Coerce+import Data.Complex (Complex)+import Data.Data+import Data.Distributive+import Data.Foldable as Foldable+import Data.Functor.Bind+import Data.Functor.Classes+import Data.Functor.Product+import Data.Functor.Rep as Rep+import Data.HashMap.Lazy (HashMap)+import Data.Hashable+import Data.Hashable.Lifted+import Data.IntMap (IntMap)+import Data.Ix+import Data.Kind+import Data.Map (Map)+#if !(MIN_VERSION_base(4,11,0))+import Data.Semigroup (Semigroup)+#endif+import Data.Serialize as Cereal+import Data.Vector (Vector)+import qualified Data.Vector.Generic.Mutable as M+import qualified Data.Vector.Generic as G+import qualified Data.Vector.Unboxed.Base as U+import Foreign.Storable+import GHC.Generics (Generic, Generic1)+import Linear.Epsilon+import Linear.Metric+import Linear.Plucker+import Linear.Quaternion+import Linear.V+import Linear.V0+import Linear.V1+import Linear.V2+import Linear.V3+import Linear.V4+import Linear.Vector+import System.Random (Random(..))++-- | An affine space is roughly a vector space in which we have+-- forgotten or at least pretend to have forgotten the origin.+--+-- > a .+^ (b .-. a) = b@+-- > (a .+^ u) .+^ v = a .+^ (u ^+^ v)@+-- > (a .-. b) ^+^ v = (a .+^ v) .-. q@+class Additive (Diff p) => Affine p where+ type Diff p :: Type -> Type++ infixl 6 .-.+ -- | Get the difference between two points as a vector offset.+ (.-.) :: Num a => p a -> p a -> Diff p a++ infixl 6 .+^+ -- | Add a vector offset to a point.+ (.+^) :: Num a => p a -> Diff p a -> p a++ infixl 6 .-^+ -- | Subtract a vector offset from a point.+ (.-^) :: Num a => p a -> Diff p a -> p a+ p .-^ v = p .+^ negated v+ {-# INLINE (.-^) #-}++instance (Affine f, Affine g) => Affine (Product f g) where+ type Diff (Product f g) = Product (Diff f) (Diff g)+ Pair a b .-. Pair c d = Pair (a .-. c) (b .-. d)+ Pair a b .+^ Pair c d = Pair (a .+^ c) (b .+^ d)+ Pair a b .-^ Pair c d = Pair (a .+^ c) (b .+^ d)++-- | Compute the quadrance of the difference (the square of the distance)+qdA :: (Affine p, Foldable (Diff p), Num a) => p a -> p a -> a+qdA a b = Foldable.sum (fmap (join (*)) (a .-. b))+{-# INLINE qdA #-}++-- | Distance between two points in an affine space+distanceA :: (Floating a, Foldable (Diff p), Affine p) => p a -> p a -> a+distanceA a b = sqrt (qdA a b)+{-# INLINE distanceA #-}++#define ADDITIVEC(CTX,T) instance CTX => Affine T where type Diff T = T ; \+ (.-.) = (^-^) ; {-# INLINE (.-.) #-} ; (.+^) = (^+^) ; {-# INLINE (.+^) #-} ; \+ (.-^) = (^-^) ; {-# INLINE (.-^) #-}+#define ADDITIVE(T) ADDITIVEC((), T)++ADDITIVE([])+ADDITIVE(Complex)+ADDITIVE(ZipList)+ADDITIVE(Maybe)+ADDITIVE(IntMap)+ADDITIVE(Identity)+ADDITIVE(Vector)+ADDITIVE(V0)+ADDITIVE(V1)+ADDITIVE(V2)+ADDITIVE(V3)+ADDITIVE(V4)+ADDITIVE(Plucker)+ADDITIVE(Quaternion)+ADDITIVE(((->) b))+ADDITIVEC(Ord k, (Map k))+ADDITIVEC((Eq k, Hashable k), (HashMap k))+ADDITIVEC(Dim n, (V n))++-- | A handy wrapper to help distinguish points from vectors at the+-- type level+newtype Point f a = P (f a)+ deriving ( Eq, Ord, Show, Read, Monad, Functor, Applicative, Foldable+ , Eq1, Ord1, Show1, Read1+ , Traversable, Apply, Additive, Metric+ , Fractional , Num, Ix, Storable, Epsilon+ , Semigroup, Monoid+ , Random, Hashable+ , Generic, Generic1, Data+ )++instance Finite f => Finite (Point f) where+ type Size (Point f) = Size f+ toV (P v) = toV v+ fromV v = P (fromV v)++instance NFData (f a) => NFData (Point f a) where+ rnf (P x) = rnf x++instance Serial1 f => Serial1 (Point f) where+ serializeWith f (P p) = serializeWith f p+ deserializeWith m = P `liftM` deserializeWith m++instance Serial (f a) => Serial (Point f a) where+ serialize (P p) = serialize p+ deserialize = P `liftM` deserialize++instance Binary (f a) => Binary (Point f a) where+ put (P p) = Binary.put p+ get = P `liftM` Binary.get++instance Serialize (f a) => Serialize (Point f a) where+ put (P p) = Cereal.put p+ get = P `liftM` Cereal.get++instance Hashable1 f => Hashable1 (Point f) where+ liftHashWithSalt h s (P f) = liftHashWithSalt h s f+ {-# INLINE liftHashWithSalt #-}++lensP :: Lens (Point f a) (Point g b) (f a) (g b)+lensP afb (P a) = P <$> afb a+{-# INLINE lensP #-}++_Point :: Iso (Point f a) (Point g b) (f a) (g b)+_Point = iso (\(P a) -> a) P+{-# INLINE _Point #-}++instance (t ~ Point g b) => Rewrapped (Point f a) t+instance Wrapped (Point f a) where+ type Unwrapped (Point f a) = f a+ _Wrapped' = _Point+ {-# INLINE _Wrapped' #-}++-- These are stolen from Data.Profunctor.Unsafe+(.#) :: Coercible b a => (b -> c) -> (a -> b) -> a -> c+f .# _ = coerce f+{-# INLINE (.#) #-}++(#.) :: Coercible c b => (b -> c) -> (a -> b) -> a -> c+(#.) _ = coerce (\x -> x :: b) :: forall a b. Coercible b a => a -> b+{-# INLINE (#.) #-}++unP :: Point f a -> f a+unP (P x) = x+{-# INLINE unP #-}++-- We can't use GND to derive 'Bind' because 'join' causes+-- role troubles. However, GHC 7.8 and above let us use+-- explicit coercions for (>>-).+instance Bind f => Bind (Point f) where+ (>>-) = ((P .) . (. (unP .))) #. (>>-) .# unP+ join (P m) = P $ m >>- \(P m') -> m'++instance Distributive f => Distributive (Point f) where+ distribute = P . collect (\(P p) -> p)+ collect = (P .) #. collect .# (unP .)++instance Representable f => Representable (Point f) where+ type Rep (Point f) = Rep f+ tabulate = P #. tabulate+ {-# INLINE tabulate #-}+ index = Rep.index .# unP+ {-# INLINE index #-}++type instance Index (Point f a) = Index (f a)+type instance IxValue (Point f a) = IxValue (f a)++instance Ixed (f a) => Ixed (Point f a) where+ ix l = lensP . ix l+ {-# INLINE ix #-}++instance Traversable f => Each (Point f a) (Point f b) a b where+ each = traverse+ {-# INLINE each #-}++instance R1 f => R1 (Point f) where+ _x = lensP . _x+ {-# INLINE _x #-}++instance R2 f => R2 (Point f) where+ _y = lensP . _y+ {-# INLINE _y #-}+ _xy = lensP . _xy+ {-# INLINE _xy #-}++instance R3 f => R3 (Point f) where+ _z = lensP . _z+ {-# INLINE _z #-}+ _xyz = lensP . _xyz+ {-# INLINE _xyz #-}++instance R4 f => R4 (Point f) where+ _w = lensP . _w+ {-# INLINE _w #-}+ _xyzw = lensP . _xyzw+ {-# INLINE _xyzw #-}++instance Additive f => Affine (Point f) where+ type Diff (Point f) = f+ (.-.) = (. unP) #. (^-^) .# unP+ {-# INLINE (.-.) #-}+ (.+^) = (P .) #. (^+^) .# unP+ {-# INLINE (.+^) #-}+ (.-^) = (P .) #. (^-^) .# unP+ {-# INLINE (.-^) #-}++-- | Vector spaces have origins.+origin :: (Additive f, Num a) => Point f a+origin = P zero++-- | An isomorphism between points and vectors, given a reference+-- point.+relative :: (Additive f, Num a) => Point f a -> Iso' (Point f a) (f a)+relative p0 = iso (.-. p0) (p0 .+^)+{-# INLINE relative #-}++newtype instance U.Vector (Point f a) = V_P (U.Vector (f a))+newtype instance U.MVector s (Point f a) = MV_P (U.MVector s (f a))+instance U.Unbox (f a) => U.Unbox (Point f a)++instance U.Unbox (f a) => M.MVector U.MVector (Point f a) where+ {-# INLINE basicLength #-}+ {-# INLINE basicUnsafeSlice #-}+ {-# INLINE basicOverlaps #-}+ {-# INLINE basicUnsafeNew #-}+ {-# INLINE basicUnsafeRead #-}+ {-# INLINE basicUnsafeWrite #-}+ basicLength (MV_P v) = M.basicLength v+ basicUnsafeSlice m n (MV_P v) = MV_P (M.basicUnsafeSlice m n v)+ basicOverlaps (MV_P v) (MV_P u) = M.basicOverlaps v u+ basicUnsafeNew n = MV_P `liftM` M.basicUnsafeNew n+ basicUnsafeRead (MV_P v) i = P `liftM` M.basicUnsafeRead v i+ basicUnsafeWrite (MV_P v) i (P x) = M.basicUnsafeWrite v i x+ basicInitialize (MV_P v) = M.basicInitialize v+ {-# INLINE basicInitialize #-}++instance U.Unbox (f a) => G.Vector U.Vector (Point f a) where+ {-# INLINE basicUnsafeFreeze #-}+ {-# INLINE basicUnsafeThaw #-}+ {-# INLINE basicLength #-}+ {-# INLINE basicUnsafeSlice #-}+ {-# INLINE basicUnsafeIndexM #-}+ basicUnsafeFreeze (MV_P v) = V_P `liftM` G.basicUnsafeFreeze v+ basicUnsafeThaw ( V_P v) = MV_P `liftM` G.basicUnsafeThaw v+ basicLength ( V_P v) = G.basicLength v+ basicUnsafeSlice m n (V_P v) = V_P (G.basicUnsafeSlice m n v)+ basicUnsafeIndexM (V_P v) i = P `liftM` G.basicUnsafeIndexM v i
+ src/Linear/Algebra.hs view
@@ -0,0 +1,136 @@+{-# LANGUAGE CPP #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE FlexibleInstances #-}+-----------------------------------------------------------------------------+-- |+-- License : BSD-style (see the file LICENSE)+-- Maintainer : Edward Kmett <ekmett@gmail.com>+-- Stability : provisional+-- Portability : portable+--+-----------------------------------------------------------------------------+module Linear.Algebra+ ( Algebra(..)+ , Coalgebra(..)+ , multRep, unitalRep+ , comultRep, counitalRep+ ) where++import Control.Lens hiding (index)+import Data.Functor.Rep+import Data.Complex+import Data.Void+import Linear.Vector+import Linear.Quaternion+import Linear.Conjugate+import Linear.V0+import Linear.V1+import Linear.V2+import Linear.V3+import Linear.V4++-- | An associative unital algebra over a ring+class Num r => Algebra r m where+ mult :: (m -> m -> r) -> m -> r+ unital :: r -> m -> r++multRep :: (Representable f, Algebra r (Rep f)) => f (f r) -> f r+multRep ffr = tabulate $ mult (index . index ffr)++unitalRep :: (Representable f, Algebra r (Rep f)) => r -> f r+unitalRep = tabulate . unital++instance Num r => Algebra r Void where+ mult _ _ = 0+ unital _ _ = 0++instance Num r => Algebra r (E V0) where+ mult _ _ = 0+ unital _ _ = 0++instance Num r => Algebra r (E V1) where+ mult f _ = f ex ex+ unital r _ = r++instance Num r => Algebra r () where+ mult f () = f () ()+ unital r () = r++instance (Algebra r a, Algebra r b) => Algebra r (a, b) where+ mult f (a,b) = mult (\a1 a2 -> mult (\b1 b2 -> f (a1,b1) (a2,b2)) b) a+ unital r (a,b) = unital r a * unital r b++instance Num r => Algebra r (E Complex) where+ mult f = \ i -> c^.el i where+ c = (f ee ee - f ei ei) :+ (f ee ei + f ei ee)+ unital r i = (r :+ 0)^.el i++instance (Num r, TrivialConjugate r) => Algebra r (E Quaternion) where+ mult f = index $ Quaternion+ (f ee ee - (f ei ei + f ej ej + f ek ek))+ (V3 (f ee ei + f ei ee + f ej ek - f ek ej)+ (f ee ej + f ej ee + f ek ei - f ei ek)+ (f ee ek + f ek ee + f ei ej - f ej ei))+ unital r = index (Quaternion r 0)++-- | A coassociative counital coalgebra over a ring+class Num r => Coalgebra r m where+ comult :: (m -> r) -> m -> m -> r+ counital :: (m -> r) -> r++comultRep :: (Representable f, Coalgebra r (Rep f)) => f r -> f (f r)+comultRep fr = tabulate $ \i -> tabulate $ \j -> comult (index fr) i j++counitalRep :: (Representable f, Coalgebra r (Rep f)) => f r -> r+counitalRep = counital . index++instance Num r => Coalgebra r Void where+ comult _ _ _ = 0+ counital _ = 0++instance Num r => Coalgebra r () where+ comult f () () = f ()+ counital f = f ()++instance Num r => Coalgebra r (E V0) where+ comult _ _ _ = 0+ counital _ = 0++instance Num r => Coalgebra r (E V1) where+ comult f _ _ = f ex+ counital f = f ex++instance Num r => Coalgebra r (E V2) where+ comult f = index . index v where+ v = V2 (V2 (f ex) 0) (V2 0 (f ey))+ counital f = f ex + f ey++instance Num r => Coalgebra r (E V3) where+ comult f = index . index q where+ q = V3 (V3 (f ex) 0 0)+ (V3 0 (f ey) 0)+ (V3 0 0 (f ez))+ counital f = f ex + f ey + f ez++instance Num r => Coalgebra r (E V4) where+ comult f = index . index v where+ v = V4 (V4 (f ex) 0 0 0) (V4 0 (f ey) 0 0) (V4 0 0 (f ez) 0) (V4 0 0 0 (f ew))+ counital f = f ex + f ey + f ez + f ew++instance Num r => Coalgebra r (E Complex) where+ comult f = \i j -> c^.el i.el j where+ c = (f ee :+ 0) :+ (0 :+ f ei)+ counital f = f ee + f ei++instance (Num r, TrivialConjugate r) => Coalgebra r (E Quaternion) where+ comult f = index . index+ (Quaternion (Quaternion (f ee) (V3 0 0 0))+ (V3 (Quaternion 0 (V3 (f ei) 0 0))+ (Quaternion 0 (V3 0 (f ej) 0))+ (Quaternion 0 (V3 0 0 (f ek)))))+ counital f = f ee + f ei + f ej + f ek++instance (Coalgebra r m, Coalgebra r n) => Coalgebra r (m, n) where+ comult f (a1, b1) (a2, b2) = comult (\a -> comult (\b -> f (a, b)) b1 b2) a1 a2+ counital k = counital $ \a -> counital $ \b -> k (a,b)
+ src/Linear/Binary.hs view
@@ -0,0 +1,27 @@+-----------------------------------------------------------------------------+-- |+-- Copyright : (C) 2013-2015 Edward Kmett and Anthony Cowley+-- License : BSD-style (see the file LICENSE)+--+-- Maintainer : Edward Kmett <ekmett@gmail.com>+-- Stability : experimental+-- Portability : non-portable+--+-- Serialization of statically-sized types with the "Data.Binary"+-- library.+------------------------------------------------------------------------------+module Linear.Binary+ ( putLinear+ , getLinear+ ) where++import Data.Binary+import Data.Foldable (traverse_)++-- | Serialize a linear type.+putLinear :: (Binary a, Foldable t) => t a -> Put+putLinear = traverse_ put++-- | Deserialize a linear type.+getLinear :: (Binary a, Applicative t, Traversable t) => Get (t a)+getLinear = sequenceA $ pure get
src/Linear/Conjugate.hs view
@@ -1,7 +1,9 @@+{-# LANGUAGE CPP #-}+{-# LANGUAGE DefaultSignatures #-}+ ----------------------------------------------------------------------------- -- |--- Module : Linear.Conjugate--- Copyright : (C) 2012-2013 Edward Kmett,+-- Copyright : (C) 2012-2015 Edward Kmett -- License : BSD-style (see the file LICENSE) -- -- Maintainer : Edward Kmett <ekmett@gmail.com>@@ -12,12 +14,18 @@ ---------------------------------------------------------------------------- module Linear.Conjugate ( Conjugate(..)+ , TrivialConjugate ) where import Data.Complex hiding (conjugate) import Data.Int import Data.Word+import Foreign.C.Types (CFloat, CDouble) +-- $setup+-- >>> import Data.Complex (Complex (..))++ -- | An involutive ring class Num a => Conjugate a where -- | Conjugate a value. This defaults to the trivial involution.@@ -28,8 +36,18 @@ -- >>> conjugate 1 -- 1 conjugate :: a -> a+#ifndef HLINT+ default conjugate :: TrivialConjugate a => a -> a conjugate = id+#endif +-- | Requires and provides a default definition such that+--+-- @+-- 'conjugate' = 'id'+-- @+class Conjugate a => TrivialConjugate a+ instance Conjugate Integer instance Conjugate Int instance Conjugate Int64@@ -43,7 +61,26 @@ instance Conjugate Word8 instance Conjugate Double instance Conjugate Float+instance Conjugate CFloat+instance Conjugate CDouble+ instance (Conjugate a, RealFloat a) => Conjugate (Complex a) where {-# SPECIALIZE instance Conjugate (Complex Float) #-} {-# SPECIALIZE instance Conjugate (Complex Double) #-} conjugate (a :+ b) = conjugate a :+ negate b++instance TrivialConjugate Integer+instance TrivialConjugate Int+instance TrivialConjugate Int64+instance TrivialConjugate Int32+instance TrivialConjugate Int16+instance TrivialConjugate Int8+instance TrivialConjugate Word+instance TrivialConjugate Word64+instance TrivialConjugate Word32+instance TrivialConjugate Word16+instance TrivialConjugate Word8+instance TrivialConjugate Double+instance TrivialConjugate Float+instance TrivialConjugate CFloat+instance TrivialConjugate CDouble
− src/Linear/Core.hs
@@ -1,25 +0,0 @@-{-# LANGUAGE RankNTypes #-}--------------------------------------------------------------------------------- |--- Module : Linear.Core--- Copyright : (C) 2012-2013 Edward Kmett,--- License : BSD-style (see the file LICENSE)------ Maintainer : Edward Kmett <ekmett@gmail.com>--- Stability : experimental--- Portability : non-portable------ Corepresentable functors as vector spaces------------------------------------------------------------------------------module Linear.Core- ( Core(..)- ) where---- |--- A 'Functor' @f@ is corepresentable if it is isomorphic to @(x -> a)@--- for some x. Nearly all such functors can be represented by choosing @x@ to be--- the set of lenses that are polymorphic in the contents of the 'Functor',--- that is to say @x = 'Rep' f@ is a valid choice of 'x' for (nearly) every--- 'Representable' 'Functor'.-class Functor f => Core f where- core :: ((forall g x. Functor g => (x -> g x) -> f x -> g (f x)) -> a) -> f a
+ src/Linear/Covector.hs view
@@ -0,0 +1,73 @@+{-# LANGUAGE CPP, MultiParamTypeClasses, FlexibleInstances, FlexibleContexts #-}+-----------------------------------------------------------------------------+-- |+-- License : BSD-style (see the file LICENSE)+-- Maintainer : Edward Kmett <ekmett@gmail.com>+-- Stability : provisional+-- Portability : portable+--+-- Operations on affine spaces.+-----------------------------------------------------------------------------+module Linear.Covector+ ( Covector(..)+ , ($*)+ ) where++import Control.Applicative+import Control.Monad+import Data.Functor.Plus hiding (zero)+import qualified Data.Functor.Plus as Plus+import Data.Functor.Bind+import Data.Functor.Rep as Rep+import Linear.Algebra++-- | Linear functionals from elements of an (infinite) free module to a scalar++newtype Covector r a = Covector { runCovector :: (a -> r) -> r }++infixr 0 $*++($*) :: Representable f => Covector r (Rep f) -> f r -> r+Covector f $* m = f (Rep.index m)++instance Functor (Covector r) where+ fmap f (Covector m) = Covector $ \k -> m (k . f)++instance Apply (Covector r) where+ Covector mf <.> Covector ma = Covector $ \k -> mf $ \f -> ma (k . f)++instance Applicative (Covector r) where+ pure a = Covector $ \k -> k a+ Covector mf <*> Covector ma = Covector $ \k -> mf $ \f -> ma $ k . f++instance Bind (Covector r) where+ Covector m >>- f = Covector $ \k -> m $ \a -> runCovector (f a) k++instance Monad (Covector r) where+#if !(MIN_VERSION_base(4,11,0))+ return a = Covector $ \k -> k a+#endif+ Covector m >>= f = Covector $ \k -> m $ \a -> runCovector (f a) k++instance Num r => Alt (Covector r) where+ Covector m <!> Covector n = Covector $ \k -> m k + n k++instance Num r => Plus (Covector r) where+ zero = Covector (const 0)++instance Num r => Alternative (Covector r) where+ Covector m <|> Covector n = Covector $ \k -> m k + n k+ empty = Covector (const 0)++instance Num r => MonadPlus (Covector r) where+ Covector m `mplus` Covector n = Covector $ \k -> m k + n k+ mzero = Covector (const 0)++instance Coalgebra r m => Num (Covector r m) where+ Covector f + Covector g = Covector $ \k -> f k + g k+ Covector f - Covector g = Covector $ \k -> f k - g k+ Covector f * Covector g = Covector $ \k -> f $ \m -> g $ comult k m+ negate (Covector f) = Covector $ \k -> negate (f k)+ abs _ = error "Covector.abs: undefined"+ signum _ = error "Covector.signum: undefined"+ fromInteger n = Covector $ \ k -> fromInteger n * counital k
src/Linear/Epsilon.hs view
@@ -1,7 +1,6 @@ ----------------------------------------------------------------------------- -- |--- Module : Linear.Epsilon--- Copyright : (C) 2012 Edward Kmett+-- Copyright : (C) 2012-2015 Edward Kmett -- License : BSD-style (see the file LICENSE) -- Maintainer : Edward Kmett <ekmett@gmail.com> -- Stability : provisional@@ -12,6 +11,8 @@ module Linear.Epsilon ( Epsilon(..) ) where+import Data.Complex (Complex, magnitude)+import Foreign.C.Types (CFloat, CDouble) -- | Provides a fairly subjective test to see if a quantity is near zero. --@@ -37,3 +38,14 @@ -- | @'abs' a '<=' 1e-12@ instance Epsilon Double where nearZero a = abs a <= 1e-12++-- | @'abs' a '<=' 1e-6@+instance Epsilon CFloat where+ nearZero a = abs a <= 1e-6++-- | @'abs' a '<=' 1e-12@+instance Epsilon CDouble where+ nearZero a = abs a <= 1e-12++instance (Epsilon a, RealFloat a) => Epsilon (Complex a) where+ nearZero = nearZero . magnitude
src/Linear/Instances.hs view
@@ -1,77 +1,14 @@-{-# OPTIONS_GHC -fno-warn-orphans #-}+{-# LANGUAGE Safe #-} ----------------------------------------------------------------------------- -- |--- Module : Linear.Instances--- Copyright : (C) 2012 Edward Kmett+-- Copyright : (C) 2012-2015 Edward Kmett -- License : BSD-style (see the file LICENSE) -- Maintainer : Edward Kmett <ekmett@gmail.com> -- Stability : provisional -- Portability : portable ----- Orphans+-- Re-exports orphan instances for @Complex@ from the @base-orphans@ package. ----------------------------------------------------------------------------- module Linear.Instances () where -import Control.Applicative-import Data.Complex-import Data.Foldable-import Data.Functor.Bind-import Data.HashMap.Lazy as HashMap-import Data.Hashable-import Data.Semigroup-import Data.Semigroup.Foldable-import Data.Semigroup.Traversable-import Data.Traversable--instance (Hashable k, Eq k) => Apply (HashMap k) where- (<.>) = HashMap.intersectionWith id--instance (Hashable k, Eq k) => Bind (HashMap k) where- -- this is needlessly painful- m >>- f = HashMap.fromList $ do- (k, a) <- HashMap.toList m- case HashMap.lookup k (f a) of- Just b -> [(k,b)]- Nothing -> []--instance Functor Complex where- fmap f (a :+ b) = f a :+ f b- {-# INLINE fmap #-}--instance Apply Complex where- (a :+ b) <.> (c :+ d) = a c :+ b d--instance Applicative Complex where- pure a = a :+ a- (a :+ b) <*> (c :+ d) = a c :+ b d--instance Bind Complex where- (a :+ b) >>- f = a' :+ b' where- a' :+ _ = f a- _ :+ b' = f b- {-# INLINE (>>-) #-}--instance Monad Complex where- return a = a :+ a- {-# INLINE return #-}-- (a :+ b) >>= f = a' :+ b' where- a' :+ _ = f a- _ :+ b' = f b- {-# INLINE (>>=) #-}--instance Foldable Complex where- foldMap f (a :+ b) = f a `mappend` f b- {-# INLINE foldMap #-}--instance Traversable Complex where- traverse f (a :+ b) = (:+) <$> f a <*> f b- {-# INLINE traverse #-}--instance Foldable1 Complex where- foldMap1 f (a :+ b) = f a <> f b- {-# INLINE foldMap1 #-}--instance Traversable1 Complex where- traverse1 f (a :+ b) = (:+) <$> f a <.> f b- {-# INLINE traverse1 #-}+import Data.Orphans ()
src/Linear/Matrix.hs view
@@ -1,7 +1,14 @@------------------------------------------------------------------------------+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE CPP #-}+{-# LANGUAGE RankNTypes #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE Trustworthy #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE TypeOperators #-}++--------------------------------------------------------------------------- -- |--- Module : Linear.Matrix--- Copyright : (C) 2012-2013 Edward Kmett,+-- Copyright : (C) 2012-2015 Edward Kmett -- License : BSD-style (see the file LICENSE) -- -- Maintainer : Edward Kmett <ekmett@gmail.com>@@ -9,67 +16,128 @@ -- Portability : non-portable -- -- Simple matrix operation for low-dimensional primitives.-----------------------------------------------------------------------------+--------------------------------------------------------------------------- module Linear.Matrix- ( (!*!), (!*) , (*!), (!!*), (*!!)+ ( (!*!), (!+!), (!-!), (!*), (*!), (!!*), (*!!), (!!/)+ , column , adjoint- , M22, M33, M44, M43, m33_to_m44, m43_to_m44- , det22, det33, inv22, inv33- , eye3, eye4- , trace+ , M22, M23, M24, M32, M33, M34, M42, M43, M44+ , m33_to_m44, m43_to_m44+ , det22, det33, det44, inv22, inv33, inv44+ , identity+ , Trace(..) , translation+ , transpose , fromQuaternion , mkTransformation+ , mkTransformationMat+ , _m22, _m23, _m24+ , _m32, _m33, _m34+ , _m42, _m43, _m44+ , lu+ , luFinite+ , forwardSub+ , forwardSubFinite+ , backwardSub+ , backwardSubFinite+ , luSolve+ , luSolveFinite+ , luInv+ , luInvFinite+ , luDet+ , luDetFinite ) where -import Control.Applicative-import Control.Monad (join)+import Control.Lens hiding (index)+import Control.Lens.Internal.Context import Data.Distributive import Data.Foldable as Foldable-import Data.Functor.Apply-import Linear.Epsilon-import Linear.Metric+import Data.Functor.Rep+import GHC.TypeLits import Linear.Quaternion+import Linear.V import Linear.V2 import Linear.V3 import Linear.V4-import Linear.Vector ((*^))+import Linear.Vector import Linear.Conjugate+import Linear.Trace -- $setup--- >>> import Data.Complex--- >>> import Data.IntMap+-- >>> import Control.Lens hiding (index)+-- >>> import Data.Complex (Complex (..))+-- >>> import Linear.V2+-- >>> import Linear.V3+-- >>> import Linear.V+-- >>> import qualified Data.IntMap as IntMap -- >>> import Debug.SimpleReflect.Vars +-- | This is a generalization of 'Control.Lens.inside' to work over any corepresentable 'Functor'.+--+-- @+-- 'column' :: 'Representable' f => 'Lens' s t a b -> 'Lens' (f s) (f t) (f a) (f b)+-- @+--+-- In practice it is used to access a column of a matrix.+--+-- >>> V2 (V3 1 2 3) (V3 4 5 6) ^._x+-- V3 1 2 3+--+-- >>> V2 (V3 1 2 3) (V3 4 5 6) ^.column _x+-- V2 1 4+column :: Representable f => LensLike (Context a b) s t a b -> Lens (f s) (f t) (f a) (f b)+column l f es = o <$> f i where+ go = l (Context id)+ i = tabulate $ \ e -> ipos $ go (index es e)+ o eb = tabulate $ \ e -> ipeek (index eb e) (go (index es e))+ infixl 7 !*!--- | Matrix product. This can compute mixed dense-dense, sparse-dense and sparse-sparse matrix products.+-- | Matrix product. This can compute any combination of sparse and dense multiplication. -- -- >>> V2 (V3 1 2 3) (V3 4 5 6) !*! V3 (V2 1 2) (V2 3 4) (V2 4 5) -- V2 (V2 19 25) (V2 43 58) ----- >>> V2 (fromList [(1,2)]) (fromList [(2,3)]) !*! fromList [(1,V3 0 0 1), (2, V3 0 0 5)]+-- >>> V2 (IntMap.fromList [(1,2)]) (IntMap.fromList [(2,3)]) !*! IntMap.fromList [(1,V3 0 0 1), (2, V3 0 0 5)] -- V2 (V3 0 0 2) (V3 0 0 15)-(!*!) :: (Functor m, Foldable r, Apply r, Distributive n, Num a) => m (r a) -> r (n a) -> m (n a)-f !*! g = fmap (\r -> Foldable.sum . liftF2 (*) r <$> g') f- where g' = distribute g+(!*!) :: (Functor m, Foldable t, Additive t, Additive n, Num a) => m (t a) -> t (n a) -> m (n a)+f !*! g = fmap (\ f' -> Foldable.foldl' (^+^) zero $ liftI2 (*^) f' g) f +infixl 6 !+!+-- | Entry-wise matrix addition.+--+-- >>> V2 (V3 1 2 3) (V3 4 5 6) !+! V2 (V3 7 8 9) (V3 1 2 3)+-- V2 (V3 8 10 12) (V3 5 7 9)+(!+!) :: (Additive m, Additive n, Num a) => m (n a) -> m (n a) -> m (n a)+as !+! bs = liftU2 (^+^) as bs++infixl 6 !-!+-- | Entry-wise matrix subtraction.+--+-- >>> V2 (V3 1 2 3) (V3 4 5 6) !-! V2 (V3 7 8 9) (V3 1 2 3)+-- V2 (V3 (-6) (-6) (-6)) (V3 3 3 3)+(!-!) :: (Additive m, Additive n, Num a) => m (n a) -> m (n a) -> m (n a)+as !-! bs = liftU2 (^-^) as bs++infixl 7 !* -- | Matrix * column vector -- -- >>> V2 (V3 1 2 3) (V3 4 5 6) !* V3 7 8 9 -- V2 50 122-infixl 7 *!-(!*) :: (Functor m, Metric r, Num a) => m (r a) -> r a -> m a-m !* v = dot v <$> m--infixl 7 !*+(!*) :: (Functor m, Foldable r, Additive r, Num a) => m (r a) -> r a -> m a+m !* v = fmap (\r -> Foldable.sum $ liftI2 (*) r v) m +infixl 7 *! -- | Row vector * matrix -- -- >>> V2 1 2 *! V2 (V3 3 4 5) (V3 6 7 8) -- V3 15 18 21-(*!) :: (Metric r, Distributive n, Num a) => r a -> r (n a) -> n a-f *! g = dot f <$> distribute g +-- (*!) :: (Metric r, Additive n, Num a) => r a -> r (n a) -> n a+-- f *! g = dot f <$> distribute g++(*!) :: (Num a, Foldable t, Additive f, Additive t) => t a -> t (f a) -> f a+f *! g = sumV $ liftI2 (*^) f g+ infixl 7 *!! -- | Scalar-matrix product --@@ -88,6 +156,12 @@ (!!*) = flip (*!!) {-# INLINE (!!*) #-} +infixl 7 !!/+-- | Matrix-scalar division+(!!/) :: (Functor m, Functor r, Fractional a) => m (r a) -> a -> m (r a)+m !!/ s = fmap (^/ s) m+{-# INLINE (!!/) #-}+ -- | Hermitian conjugate or conjugate transpose -- -- >>> adjoint (V2 (V2 (1 :+ 2) (3 :+ 4)) (V2 (5 :+ 6) (7 :+ 8)))@@ -96,46 +170,59 @@ adjoint = collect (fmap conjugate) {-# INLINE adjoint #-} --- | Compute the trace of a matrix------ >>> trace (V2 (V2 a b) (V2 c d))--- a + d-trace :: (Monad f, Foldable f, Num a) => f (f a) -> a-trace m = Foldable.sum (join m)-{-# INLINE trace #-}- -- * Matrices -- -- Matrices use a row-major representation. -- | A 2x2 matrix with row-major representation type M22 a = V2 (V2 a)+-- | A 2x3 matrix with row-major representation+type M23 a = V2 (V3 a)+-- | A 2x4 matrix with row-major representation+type M24 a = V2 (V4 a)+-- | A 3x2 matrix with row-major representation+type M32 a = V3 (V2 a) -- | A 3x3 matrix with row-major representation type M33 a = V3 (V3 a)--- | A 4x4 matrix with row-major representation-type M44 a = V4 (V4 a)+-- | A 3x4 matrix with row-major representation+type M34 a = V3 (V4 a)+-- | A 4x2 matrix with row-major representation+type M42 a = V4 (V2 a) -- | A 4x3 matrix with row-major representation type M43 a = V4 (V3 a)+-- | A 4x4 matrix with row-major representation+type M44 a = V4 (V4 a) -- | Build a rotation matrix from a unit 'Quaternion'. fromQuaternion :: Num a => Quaternion a -> M33 a fromQuaternion (Quaternion w (V3 x y z)) =- V3 (V3 (1-2*(y2+z2)) (2*(x*y-z*w)) (2*(x*z+y*w)))- (V3 (2*(x*y+z*w)) (1-2*(x2+z2)) (2*(y*z-x*w)))- (V3 (2*(x*z-y*w)) (2*(y*z+x*w)) (1-2*(x2+y2)))- where x2 = x * x- y2 = y * y- z2 = z * z+ V3 (V3 (1-2*(y2+z2)) (2*(xy-zw)) (2*(xz+yw)))+ (V3 (2*(xy+zw)) (1-2*(x2+z2)) (2*(yz-xw)))+ (V3 (2*(xz-yw)) (2*(yz+xw)) (1-2*(x2+y2)))+ where 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+{-# INLINE fromQuaternion #-} +-- | Build a transformation matrix from a rotation matrix and a+-- translation vector. mkTransformationMat :: Num a => M33 a -> V3 a -> M44 a mkTransformationMat (V3 r1 r2 r3) (V3 tx ty tz) = V4 (snoc3 r1 tx) (snoc3 r2 ty) (snoc3 r3 tz) (V4 0 0 0 1) where snoc3 (V3 x y z) = V4 x y z+{-# INLINE mkTransformationMat #-} -- |Build a transformation matrix from a rotation expressed as a -- 'Quaternion' and a translation vector. mkTransformation :: Num a => Quaternion a -> V3 a -> M44 a mkTransformation = mkTransformationMat . fromQuaternion+{-# INLINE mkTransformation #-} -- | Convert from a 4x3 matrix to a 4x4 matrix, extending it with the @[ 0 0 0 1 ]@ column vector m43_to_m44 :: Num a => M43 a -> M44 a@@ -148,39 +235,78 @@ (V4 d e f 0) (V4 g h i 0) (V4 j k l 1)-{-# ANN m43_to_m44 "HLint: ignore Use camelCase" #-} -- | Convert a 3x3 matrix to a 4x4 matrix extending it with 0's in the new row and column. m33_to_m44 :: Num a => M33 a -> M44 a m33_to_m44 (V3 r1 r2 r3) = V4 (vector r1) (vector r2) (vector r3) (point 0)-{-# ANN m33_to_m44 "HLint: ignore Use camelCase" #-} --- |3x3 identity matrix.+-- |The identity matrix for any dimension vector. ----- >>> eye3+-- >>> identity :: M44 Int+-- V4 (V4 1 0 0 0) (V4 0 1 0 0) (V4 0 0 1 0) (V4 0 0 0 1)+-- >>> identity :: V3 (V3 Int) -- V3 (V3 1 0 0) (V3 0 1 0) (V3 0 0 1)-eye3 :: Num a => M33 a-eye3 = V3 (V3 1 0 0)- (V3 0 1 0)- (V3 0 0 1)+identity :: (Num a, Traversable t, Applicative t) => t (t a)+identity = scaled (pure 1) --- |4x4 identity matrix.------ >>> eye4--- V4 (V4 1 0 0 0) (V4 0 1 0 0) (V4 0 0 1 0) (V4 0 0 0 1)-eye4 :: Num a => M44 a-eye4 = V4 (V4 1 0 0 0)- (V4 0 1 0 0)- (V4 0 0 1 0)- (V4 0 0 0 1)+-- |Extract the translation vector (first three entries of the last+-- column) from a 3x4 or 4x4 matrix.+translation :: (Representable t, R3 t, R4 v) => Lens' (t (v a)) (V3 a)+translation = column _w._xyz+{-+translation f rs = aux <$> f (view _w <$> view _xyz rs)+ where aux (V3 x y z) = (_x._w .~ x) . (_y._w .~ y) . (_z._w .~ z) $ rs +-- translation :: (R3 t, R4 v, Functor f, Functor t) => (V3 a -> f (V3 a)) -> t (v a) -> f (t a)+-- translation = (. fmap (^._w)) . _xyz where+-- x ^. l = getConst (l Const x)+-} --- |Extract the translation vector (first three entries of the last--- column) from a 3x4 or 4x4 matrix-translation :: (R3 t, R4 v, Functor f, Functor t) => (V3 a -> f (V3 a)) -> t (v a) -> f (t a)-translation = (. fmap (^._w)) . _xyz where- x ^. l = getConst (l Const x)+-- |Extract a 2x2 matrix from a matrix of higher dimensions by dropping excess+-- rows and columns.+_m22 :: (Representable t, R2 t, R2 v) => Lens' (t (v a)) (M22 a)+_m22 = column _xy._xy +-- |Extract a 2x3 matrix from a matrix of higher dimensions by dropping excess+-- rows and columns.+_m23 :: (Representable t, R2 t, R3 v) => Lens' (t (v a)) (M23 a)+_m23 = column _xyz._xy++-- |Extract a 2x4 matrix from a matrix of higher dimensions by dropping excess+-- rows and columns.+_m24 :: (Representable t, R2 t, R4 v) => Lens' (t (v a)) (M24 a)+_m24 = column _xyzw._xy++-- |Extract a 3x2 matrix from a matrix of higher dimensions by dropping excess+-- rows and columns.+_m32 :: (Representable t, R3 t, R2 v) => Lens' (t (v a)) (M32 a)+_m32 = column _xy._xyz++-- |Extract a 3x3 matrix from a matrix of higher dimensions by dropping excess+-- rows and columns.+_m33 :: (Representable t, R3 t, R3 v) => Lens' (t (v a)) (M33 a)+_m33 = column _xyz._xyz++-- |Extract a 3x4 matrix from a matrix of higher dimensions by dropping excess+-- rows and columns.+_m34 :: (Representable t, R3 t, R4 v) => Lens' (t (v a)) (M34 a)+_m34 = column _xyzw._xyz++-- |Extract a 4x2 matrix from a matrix of higher dimensions by dropping excess+-- rows and columns.+_m42 :: (Representable t, R4 t, R2 v) => Lens' (t (v a)) (M42 a)+_m42 = column _xy._xyzw++-- |Extract a 4x3 matrix from a matrix of higher dimensions by dropping excess+-- rows and columns.+_m43 :: (Representable t, R4 t, R3 v) => Lens' (t (v a)) (M43 a)+_m43 = column _xyz._xyzw++-- |Extract a 4x4 matrix from a matrix of higher dimensions by dropping excess+-- rows and columns.+_m44 :: (Representable t, R4 t, R4 v) => Lens' (t (v a)) (M44 a)+_m44 = column _xyzw._xyzw+ -- |2x2 matrix determinant. -- -- >>> det22 (V2 (V2 a b) (V2 c d))@@ -199,29 +325,49 @@ (V3 g h i)) = a * (e*i-f*h) - d * (b*i-c*h) + g * (b*f-c*e) {-# INLINE det33 #-} +-- |4x4 matrix determinant.+det44 :: Num a => M44 a -> a+det44 (V4 (V4 i00 i01 i02 i03)+ (V4 i10 i11 i12 i13)+ (V4 i20 i21 i22 i23)+ (V4 i30 i31 i32 i33)) =+ let+ s0 = i00 * i11 - i10 * i01+ s1 = i00 * i12 - i10 * i02+ s2 = i00 * i13 - i10 * i03+ s3 = i01 * i12 - i11 * i02+ s4 = i01 * i13 - i11 * i03+ s5 = i02 * i13 - i12 * i03++ c5 = i22 * i33 - i32 * i23+ c4 = i21 * i33 - i31 * i23+ c3 = i21 * i32 - i31 * i22+ c2 = i20 * i33 - i30 * i23+ c1 = i20 * i32 - i30 * i22+ c0 = i20 * i31 - i30 * i21+ in s0 * c5 - s1 * c4 + s2 * c3 + s3 * c2 - s4 * c1 + s5 * c0+{-# INLINE det44 #-}+ -- |2x2 matrix inverse. -- -- >>> inv22 $ V2 (V2 1 2) (V2 3 4)--- Just (V2 (V2 (-2.0) 1.0) (V2 1.5 (-0.5)))-inv22 :: (Epsilon a, Floating a) => M22 a -> Maybe (M22 a)-inv22 m@(V2 (V2 a b) (V2 c d))- | nearZero det = Nothing- | otherwise = Just $ (1 / det) *!! V2 (V2 d (-b)) (V2 (-c) a)+-- V2 (V2 (-2.0) 1.0) (V2 1.5 (-0.5))+inv22 :: Fractional a => M22 a -> M22 a+inv22 m@(V2 (V2 a b) (V2 c d)) = (1 / det) *!! V2 (V2 d (-b)) (V2 (-c) a) where det = det22 m {-# INLINE inv22 #-} -- |3x3 matrix inverse. -- -- >>> inv33 $ V3 (V3 1 2 4) (V3 4 2 2) (V3 1 1 1)--- Just (V3 (V3 0.0 0.5 (-1.0)) (V3 (-0.5) (-0.75) 3.5) (V3 0.5 0.25 (-1.5)))-inv33 :: (Epsilon a, Floating a) => M33 a -> Maybe (M33 a)+-- V3 (V3 0.0 0.5 (-1.0)) (V3 (-0.5) (-0.75) 3.5) (V3 0.5 0.25 (-1.5))+inv33 :: Fractional a => M33 a -> M33 a inv33 m@(V3 (V3 a b c) (V3 d e f) (V3 g h i))- | nearZero det = Nothing- | otherwise = Just $ (1 / det) *!! V3 (V3 a' b' c')- (V3 d' e' f')- (V3 g' h' i')+ = (1 / det) *!! V3 (V3 a' b' c')+ (V3 d' e' f')+ (V3 g' h' i') where a' = cofactor (e,f,h,i) b' = cofactor (c,b,i,h) c' = cofactor (b,c,e,f)@@ -234,3 +380,352 @@ cofactor (q,r,s,t) = det22 (V2 (V2 q r) (V2 s t)) det = det33 m {-# INLINE inv33 #-}+++-- | 'transpose' is just an alias for 'distribute'+--+-- > transpose (V3 (V2 1 2) (V2 3 4) (V2 5 6))+-- V2 (V3 1 3 5) (V3 2 4 6)+transpose :: (Distributive g, Functor f) => f (g a) -> g (f a)+transpose = distribute+{-# INLINE transpose #-}++-- |4x4 matrix inverse.+inv44 :: Fractional a => M44 a -> M44 a+inv44 (V4 (V4 i00 i01 i02 i03)+ (V4 i10 i11 i12 i13)+ (V4 i20 i21 i22 i23)+ (V4 i30 i31 i32 i33)) =+ let s0 = i00 * i11 - i10 * i01+ s1 = i00 * i12 - i10 * i02+ s2 = i00 * i13 - i10 * i03+ s3 = i01 * i12 - i11 * i02+ s4 = i01 * i13 - i11 * i03+ s5 = i02 * i13 - i12 * i03+ c5 = i22 * i33 - i32 * i23+ c4 = i21 * i33 - i31 * i23+ c3 = i21 * i32 - i31 * i22+ c2 = i20 * i33 - i30 * i23+ c1 = i20 * i32 - i30 * i22+ c0 = i20 * i31 - i30 * i21+ det = s0 * c5 - s1 * c4 + s2 * c3 + s3 * c2 - s4 * c1 + s5 * c0+ invDet = recip det+ in invDet *!! V4 (V4 (i11 * c5 - i12 * c4 + i13 * c3)+ (-i01 * c5 + i02 * c4 - i03 * c3)+ (i31 * s5 - i32 * s4 + i33 * s3)+ (-i21 * s5 + i22 * s4 - i23 * s3))+ (V4 (-i10 * c5 + i12 * c2 - i13 * c1)+ (i00 * c5 - i02 * c2 + i03 * c1)+ (-i30 * s5 + i32 * s2 - i33 * s1)+ (i20 * s5 - i22 * s2 + i23 * s1))+ (V4 (i10 * c4 - i11 * c2 + i13 * c0)+ (-i00 * c4 + i01 * c2 - i03 * c0)+ (i30 * s4 - i31 * s2 + i33 * s0)+ (-i20 * s4 + i21 * s2 - i23 * s0))+ (V4 (-i10 * c3 + i11 * c1 - i12 * c0)+ (i00 * c3 - i01 * c1 + i02 * c0)+ (-i30 * s3 + i31 * s1 - i32 * s0)+ (i20 * s3 - i21 * s1 + i22 * s0))+{-# INLINE inv44 #-}++-- | Compute the (L, U) decomposition of a square matrix using Crout's+-- algorithm. The 'Index' of the vectors must be 'Integral'.+lu :: ( Num a+ , Fractional a+ , Foldable m+ , Traversable m+ , Applicative m+ , Additive m+ , Ixed (m a)+ , Ixed (m (m a))+ , i ~ Index (m a)+ , i ~ Index (m (m a))+ , Eq i+ , Integral i+ , a ~ IxValue (m a)+ , m a ~ IxValue (m (m a))+ , Num (m a)+ )+ => m (m a)+ -> (m (m a), m (m a))+lu a =+ let n = fromIntegral (length a)+ initU = identity+ initL = zero+ buildLVal !i !j !l !u =+ let go !k !s+ | k == j = s+ | otherwise = go (k+1)+ ( s+ + ( (l ^?! ix i ^?! ix k)+ * (u ^?! ix k ^?! ix j)+ )+ )+ s' = go 0 0+ in l & (ix i . ix j) .~ ((a ^?! ix i ^?! ix j) - s')+ buildL !i !j !l !u+ | i == n = l+ | otherwise = buildL (i+1) j (buildLVal i j l u) u+ buildUVal !i !j !l !u =+ let go !k !s+ | k == j = s+ | otherwise = go (k+1)+ ( s+ + ( (l ^?! ix j ^?! ix k)+ * (u ^?! ix k ^?! ix i)+ )+ )+ s' = go 0 0+ in u & (ix j . ix i) .~ ( ((a ^?! ix j ^?! ix i) - s')+ / (l ^?! ix j ^?! ix j)+ )+ buildU !i !j !l !u+ | i == n = u+ | otherwise = buildU (i+1) j l (buildUVal i j l u)+ buildLU !j !l !u+ | j == n = (l, u)+ | otherwise =+ let l' = buildL j j l u+ u' = buildU j j l' u+ in buildLU (j+1) l' u'+ in buildLU 0 initL initU++-- | Compute the (L, U) decomposition of a square matrix using Crout's+-- algorithm, using the vector's 'Finite' instance to provide an index.+luFinite :: ( Num a+ , Fractional a+ , Functor m+ , Finite m+ , n ~ Size m+ , KnownNat n+ , Num (m a)+ )+ => m (m a)+ -> (m (m a), m (m a))+luFinite a =+ bimap (fmap fromV . fromV)+ (fmap fromV . fromV)+ (lu (fmap toV (toV a)))++-- | Solve a linear system with a lower-triangular matrix of coefficients with+-- forwards substitution.+forwardSub :: ( Num a+ , Fractional a+ , Foldable m+ , Additive m+ , Ixed (m a)+ , Ixed (m (m a))+ , i ~ Index (m a)+ , i ~ Index (m (m a))+ , Eq i+ , Ord i+ , Integral i+ , a ~ IxValue (m a)+ , m a ~ IxValue (m (m a))+ )+ => m (m a)+ -> m a+ -> m a+forwardSub a b =+ let n = fromIntegral (length b)+ initX = zero+ coeff !i !j !s !x+ | j == i = s+ | otherwise = coeff i (j+1) (s + ((a ^?! ix i ^?! ix j) * (x ^?! ix j))) x+ go !i !x+ | i == n = x+ | otherwise = go (i + 1) (x & ix i .~ ( ((b ^?! ix i) - coeff i 0 0 x)+ / (a ^?! ix i ^?! ix i)+ ))+ in go 0 initX++-- | Solve a linear system with a lower-triangular matrix of coefficients with+-- forwards substitution, using the vector's 'Finite' instance to provide an+-- index.+forwardSubFinite :: ( Num a+ , Fractional a+ , Foldable m+ , n ~ Size m+ , KnownNat n+ , Additive m+ , Finite m+ )+ => m (m a)+ -> m a+ -> m a+forwardSubFinite a b = fromV (forwardSub (fmap toV (toV a)) (toV b))++-- | Solve a linear system with an upper-triangular matrix of coefficients with+-- backwards substitution.+backwardSub :: ( Num a+ , Fractional a+ , Foldable m+ , Additive m+ , Ixed (m a)+ , Ixed (m (m a))+ , i ~ Index (m a)+ , i ~ Index (m (m a))+ , Eq i+ , Ord i+ , Integral i+ , a ~ IxValue (m a)+ , m a ~ IxValue (m (m a))+ )+ => m (m a)+ -> m a+ -> m a+backwardSub a b =+ let n = fromIntegral (length b)+ initX = zero+ coeff !i !j !s !x+ | j == n = s+ | otherwise = coeff i+ (j+1)+ (s + ((a ^?! ix i ^?! ix j) * (x ^?! ix j)))+ x+ go !i !x+ | i < 0 = x+ | otherwise = go (i-1)+ (x & ix i .~ ( ((b ^?! ix i) - coeff i (i+1) 0 x)+ / (a ^?! ix i ^?! ix i)+ ))+ in go (n-1) initX++-- | Solve a linear system with an upper-triangular matrix of coefficients with+-- backwards substitution, using the vector's 'Finite' instance to provide an+-- index.+backwardSubFinite :: ( Num a+ , Fractional a+ , Foldable m+ , n ~ Size m+ , KnownNat n+ , Additive m+ , Finite m+ )+ => m (m a)+ -> m a+ -> m a+backwardSubFinite a b = fromV (backwardSub (fmap toV (toV a)) (toV b))++-- | Solve a linear system with LU decomposition.+luSolve :: ( Num a+ , Fractional a+ , Foldable m+ , Traversable m+ , Applicative m+ , Additive m+ , Ixed (m a)+ , Ixed (m (m a))+ , i ~ Index (m a)+ , i ~ Index (m (m a))+ , Eq i+ , Integral i+ , a ~ IxValue (m a)+ , m a ~ IxValue (m (m a))+ , Num (m a)+ )+ => m (m a)+ -> m a+ -> m a+luSolve a b =+ let (l, u) = lu a+ in backwardSub u (forwardSub l b)++-- | Solve a linear system with LU decomposition, using the vector's 'Finite'+-- instance to provide an index.+luSolveFinite :: ( Num a+ , Fractional a+ , Functor m+ , Finite m+ , n ~ Size m+ , KnownNat n+ , Num (m a)+ )+ => m (m a)+ -> m a+ -> m a+luSolveFinite a b = fromV (luSolve (fmap toV (toV a)) (toV b))++-- | Invert a matrix with LU decomposition.+luInv :: ( Num a+ , Fractional a+ , Foldable m+ , Traversable m+ , Applicative m+ , Additive m+ , Distributive m+ , Ixed (m a)+ , Ixed (m (m a))+ , i ~ Index (m a)+ , i ~ Index (m (m a))+ , Eq i+ , Integral i+ , a ~ IxValue (m a)+ , m a ~ IxValue (m (m a))+ , Num (m a)+ )+ => m (m a)+ -> m (m a)+luInv a =+ let n = fromIntegral (length a)+ initA' = zero+ (l, u) = lu a+ go !i !a'+ | i == n = a'+ | otherwise = let e = zero & ix i .~ 1+ a'r = backwardSub u (forwardSub l e)+ in go (i+1) (a' & ix i .~ a'r)+ in transpose (go 0 initA')++-- | Invert a matrix with LU decomposition, using the vector's 'Finite' instance+-- to provide an index.+luInvFinite :: ( Num a+ , Fractional a+ , Functor m+ , Finite m+ , n ~ Size m+ , KnownNat n+ , Num (m a)+ )+ => m (m a)+ -> m (m a)+luInvFinite a = fmap fromV (fromV (luInv (fmap toV (toV a))))++-- | Compute the determinant of a matrix using LU decomposition.+luDet :: ( Num a+ , Fractional a+ , Foldable m+ , Traversable m+ , Applicative m+ , Additive m+ , Trace m+ , Ixed (m a)+ , Ixed (m (m a))+ , i ~ Index (m a)+ , i ~ Index (m (m a))+ , Eq i+ , Integral i+ , a ~ IxValue (m a)+ , m a ~ IxValue (m (m a))+ , Num (m a)+ )+ => m (m a)+ -> a+luDet a =+ let (l, u) = lu a+ p = Foldable.foldl (*) 1+ in p (diagonal l) * p (diagonal u)++-- | Compute the determinant of a matrix using LU decomposition, using the+-- vector's 'Finite' instance to provide an index.+luDetFinite :: ( Num a+ , Fractional a+ , Functor m+ , Finite m+ , n ~ Size m+ , KnownNat n+ , Num (m a)+ )+ => m (m a)+ -> a+luDetFinite = luDet . fmap toV . toV
src/Linear/Metric.hs view
@@ -1,8 +1,9 @@+{-# LANGUAGE CPP #-} {-# LANGUAGE DefaultSignatures #-}+{-# LANGUAGE Trustworthy #-} ----------------------------------------------------------------------------- -- |--- Module : Linear.Metric--- Copyright : (C) 2012-2013 Edward Kmett,+-- Copyright : (C) 2012-2015 Edward Kmett -- License : BSD-style (see the file LICENSE) -- -- Maintainer : Edward Kmett <ekmett@gmail.com>@@ -12,16 +13,25 @@ -- Free metric spaces ---------------------------------------------------------------------------- module Linear.Metric- ( Metric(..), normalize+ ( Metric(..), normalize, project ) where +import Control.Applicative import Data.Foldable as Foldable-import Data.Functor.Apply+import Data.Functor.Compose+import Data.Functor.Identity+import Data.Functor.Product+import Data.Vector (Vector)+import Data.IntMap (IntMap)+import Data.Map (Map)+import Data.HashMap.Strict (HashMap)+import Data.Hashable (Hashable) import Linear.Epsilon import Linear.Vector -- $setup -- >>> import Linear+-- -- | Free and sparse inner product/metric spaces. class Additive f => Metric f where@@ -31,8 +41,10 @@ -- >>> V2 1 2 `dot` V2 3 4 -- 11 dot :: Num a => f a -> f a -> a+#ifndef HLINT default dot :: (Foldable f, Num a) => f a -> f a -> a- dot x y = Foldable.sum $ liftF2 (*) x y+ dot x y = Foldable.sum $ liftI2 (*) x y+#endif -- | Compute the squared norm. The name quadrance arises from -- Norman J. Wildberger's rational trigonometry.@@ -56,8 +68,43 @@ signorm v = fmap (/m) v where m = norm v +instance (Metric f, Metric g) => Metric (Product f g) where+ dot (Pair a b) (Pair c d) = dot a c + dot b d+ quadrance (Pair a b) = quadrance a + quadrance b+ qd (Pair a b) (Pair c d) = qd a c + qd b d+ distance p q = sqrt (qd p q)++instance (Metric f, Metric g) => Metric (Compose f g) where+ dot (Compose a) (Compose b) = quadrance (liftI2 dot a b)+ quadrance = quadrance . fmap quadrance . getCompose+ qd (Compose a) (Compose b) = quadrance (liftI2 qd a b)+ distance (Compose a) (Compose b) = norm (liftI2 qd a b)++instance Metric Identity where+ dot (Identity x) (Identity y) = x * y++instance Metric []++instance Metric Maybe++instance Metric ZipList where+ -- ZipList is missing its Foldable instance+ dot (ZipList x) (ZipList y) = dot x y++instance Metric IntMap++instance Ord k => Metric (Map k)++instance (Hashable k, Eq k) => Metric (HashMap k)++instance Metric Vector+ -- | Normalize a 'Metric' functor to have unit 'norm'. This function -- does not change the functor if its 'norm' is 0 or 1. normalize :: (Floating a, Metric f, Epsilon a) => f a -> f a normalize v = if nearZero l || nearZero (1-l) then v else fmap (/sqrt l) v where l = quadrance v++-- | @project u v@ computes the projection of @v@ onto @u@.+project :: (Metric v, Fractional a) => v a -> v a -> v a+project u v = ((v `dot` u) / quadrance u) *^ u
src/Linear/Plucker.hs view
@@ -1,8 +1,17 @@+{-# LANGUAGE CPP #-} {-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE GADTs #-}+{-# LANGUAGE Trustworthy #-}+{-# LANGUAGE DeriveGeneric #-}+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE DeriveLift #-}+ ----------------------------------------------------------------------------- -- |--- Module : Linear.Plucker--- Copyright : (C) 2012-2013 Edward Kmett,+-- Copyright : (C) 2012-2015 Edward Kmett -- License : BSD-style (see the file LICENSE) -- -- Maintainer : Edward Kmett <ekmett@gmail.com>@@ -17,36 +26,109 @@ , isotropic , (><) , plucker+ , plucker3D+ -- * Operations on lines+ , parallel , intersects+ , LinePass(..)+ , passes+ , quadranceToOrigin+ , closestToOrigin+ , isLine+ , coincides+ , coincides' -- * Basis elements , p01, p02, p03 , p10, p12, p13 , p20, p21, p23 , p30, p31, p32++ , e01, e02, e03, e12, e31, e23 ) where +#if !MIN_VERSION_base(4,18,0) import Control.Applicative+#endif+import Control.DeepSeq (NFData(rnf))+import Control.Monad (liftM)+import Control.Monad.Fix+import Control.Monad.Zip+import Control.Lens as Lens hiding (index, (<.>))+import Data.Binary as Binary+import Data.Bytes.Serial import Data.Distributive import Data.Foldable as Foldable+import qualified Data.Foldable.WithIndex as WithIndex import Data.Functor.Bind+import Data.Functor.Classes+import Data.Functor.Rep+import qualified Data.Functor.WithIndex as WithIndex+import Data.Hashable import Data.Semigroup import Data.Semigroup.Foldable-import Data.Semigroup.Traversable-import Data.Traversable+import Data.Serialize as Cereal+import qualified Data.Traversable.WithIndex as WithIndex+import qualified Data.Vector as V+import qualified Data.Vector.Generic.Mutable as M+import qualified Data.Vector.Generic as G+import qualified Data.Vector.Unboxed.Base as U import Foreign.Ptr (castPtr) import Foreign.Storable (Storable(..)) import GHC.Arr (Ix(..))-import Linear.Core+import GHC.Generics (Generic, Generic1)+#if defined(MIN_VERSION_template_haskell)+import Language.Haskell.TH.Syntax (Lift)+#endif import Linear.Epsilon import Linear.Metric+import Linear.V+import Linear.V2+import Linear.V3 import Linear.V4 import Linear.Vector--{-# ANN module "HLint: ignore Reduce duplication" #-}+import System.Random (Random(..), Uniform)+import System.Random.Stateful (UniformRange(..)) -- | Plücker coordinates for lines in a 3-dimensional space.-data Plucker a = Plucker a a a a a a deriving (Eq,Ord,Show,Read)+data Plucker a = Plucker !a !a !a !a !a !a deriving (Eq,Ord,Show,Read+ ,Generic,Generic1+#if defined(MIN_VERSION_template_haskell)+ ,Lift+#endif+ ) +instance Finite Plucker where+ type Size Plucker = 6+ toV (Plucker a b c d e f) = V (V.fromListN 6 [a,b,c,d,e,f])+ fromV (V v) = Plucker (v V.! 0) (v V.! 1) (v V.! 2) (v V.! 3) (v V.! 4) (v V.! 5)++instance Random a => Random (Plucker a) where+ random g = case random g of+ (a, g1) -> case random g1 of+ (b, g2) -> case random g2 of+ (c, g3) -> case random g3 of+ (d, g4) -> case random g4 of+ (e, g5) -> case random g5 of+ (f, g6) -> (Plucker a b c d e f, g6)+ randomR (Plucker a b c d e f, Plucker a' b' c' d' e' f') g = case randomR (a,a') g of+ (a'', g1) -> case randomR (b,b') g1 of+ (b'', g2) -> case randomR (c,c') g2 of+ (c'', g3) -> case randomR (d,d') g3 of+ (d'', g4) -> case randomR (e,e') g4 of+ (e'', g5) -> case randomR (f,f') g5 of+ (f'', g6) -> (Plucker a'' b'' c'' d'' e'' f'', g6)++instance Uniform a => Uniform (Plucker a) where++instance UniformRange a => UniformRange (Plucker a) where+ uniformRM (Plucker a b c d e f, Plucker a' b' c' d' e' f') g = Plucker+ <$> uniformRM (a, a') g+ <*> uniformRM (b, b') g+ <*> uniformRM (c, c') g+ <*> uniformRM (d, d') g+ <*> uniformRM (e, e') g+ <*> uniformRM (f, f') g+ instance Functor Plucker where fmap g (Plucker a b c d e f) = Plucker (g a) (g b) (g c) (g d) (g e) (g f) {-# INLINE fmap #-}@@ -63,7 +145,13 @@ Plucker (a g) (b h) (c i) (d j) (e k) (f l) {-# INLINE (<*>) #-} -instance Additive Plucker+instance Additive Plucker where+ zero = pure 0+ {-# INLINE zero #-}+ liftU2 = liftA2+ {-# INLINE liftU2 #-}+ liftI2 = liftA2+ {-# INLINE liftI2 #-} instance Bind Plucker where Plucker a b c d e f >>- g = Plucker a' b' c' d' e' f' where@@ -76,8 +164,10 @@ {-# INLINE (>>-) #-} instance Monad Plucker where+#if !(MIN_VERSION_base(4,11,0)) return a = Plucker a a a a a a {-# INLINE return #-}+#endif Plucker a b c d e f >>= g = Plucker a' b' c' d' e' f' where Plucker a' _ _ _ _ _ = g a Plucker _ b' _ _ _ _ = g b@@ -96,14 +186,19 @@ (fmap (\(Plucker _ _ _ _ _ x) -> x) f) {-# INLINE distribute #-} -instance Core Plucker where- core f = Plucker (f p01) (f p02) (f p03) (f p23) (f p31) (f p12)- {-# INLINE core #-}+instance Representable Plucker where+ type Rep Plucker = E Plucker+ tabulate f = Plucker (f e01) (f e02) (f e03) (f e23) (f e31) (f e12)+ {-# INLINE tabulate #-}+ index xs (E l) = view l xs+ {-# INLINE index #-} instance Foldable Plucker where foldMap g (Plucker a b c d e f) = g a `mappend` g b `mappend` g c `mappend` g d `mappend` g e `mappend` g f {-# INLINE foldMap #-}+ null _ = False+ length _ = 6 instance Traversable Plucker where traverse g (Plucker a b c d e f) =@@ -170,6 +265,48 @@ fromRational = pure . fromRational {-# INLINE fromRational #-} +instance Floating a => Floating (Plucker a) where+ pi = pure pi+ {-# INLINE pi #-}+ exp = fmap exp+ {-# INLINE exp #-}+ sqrt = fmap sqrt+ {-# INLINE sqrt #-}+ log = fmap log+ {-# INLINE log #-}+ (**) = liftA2 (**)+ {-# INLINE (**) #-}+ logBase = liftA2 logBase+ {-# INLINE logBase #-}+ sin = fmap sin+ {-# INLINE sin #-}+ tan = fmap tan+ {-# INLINE tan #-}+ cos = fmap cos+ {-# INLINE cos #-}+ asin = fmap asin+ {-# INLINE asin #-}+ atan = fmap atan+ {-# INLINE atan #-}+ acos = fmap acos+ {-# INLINE acos #-}+ sinh = fmap sinh+ {-# INLINE sinh #-}+ tanh = fmap tanh+ {-# INLINE tanh #-}+ cosh = fmap cosh+ {-# INLINE cosh #-}+ asinh = fmap asinh+ {-# INLINE asinh #-}+ atanh = fmap atanh+ {-# INLINE atanh #-}+ acosh = fmap acosh+ {-# INLINE acosh #-}++instance Hashable a => Hashable (Plucker a) where+ hashWithSalt s (Plucker a b c d e f) = s `hashWithSalt` a `hashWithSalt` b `hashWithSalt` c `hashWithSalt` d `hashWithSalt` e `hashWithSalt` f+ {-# INLINE hashWithSalt #-}+ instance Storable a => Storable (Plucker a) where sizeOf _ = 6 * sizeOf (undefined::a) {-# INLINE sizeOf #-}@@ -193,30 +330,46 @@ where ptr' = castPtr ptr {-# INLINE peek #-} --- | Given a pair of points represented by homogeneous coordinates generate Plücker coordinates--- for the line through them.+instance Metric Plucker where+ dot (Plucker a b c d e f) (Plucker g h i j k l) = a*g+b*h+c*i+d*j+e*k+f*l+ {-# INLINE dot #-}++instance Epsilon a => Epsilon (Plucker a) where+ nearZero = nearZero . quadrance+ {-# INLINE nearZero #-}++-- | Given a pair of points represented by homogeneous coordinates+-- generate Plücker coordinates for the line through them, directed+-- from the second towards the first. plucker :: Num a => V4 a -> V4 a -> Plucker a plucker (V4 a b c d) (V4 e f g h) = Plucker (a*f-b*e) (a*g-c*e)- (a*d-h*e)- (c*h-d*g)- (d*f-b*h) (b*g-c*f)+ (a*h-d*e)+ (b*h-d*f)+ (c*h-d*g) {-# INLINE plucker #-} +-- | Given a pair of 3D points, generate Plücker coordinates for the+-- line through them, directed from the second towards the first.+plucker3D :: Num a => V3 a -> V3 a -> Plucker a+plucker3D p q = Plucker a b c d e f+ where V3 a b c = p - q+ V3 d e f = p `cross` q+ -- | These elements form a basis for the Plücker space, or the Grassmanian manifold @Gr(2,V4)@. -- -- @--- 'p01' :: Lens' ('Plucker' a) a--- 'p02' :: Lens' ('Plucker' a) a--- 'p03' :: Lens' ('Plucker' a) a--- 'p23' :: Lens' ('Plucker' a) a--- 'p31' :: Lens' ('Plucker' a) a--- 'p12' :: Lens' ('Plucker' a) a+-- 'p01' :: 'Lens'' ('Plucker' a) a+-- 'p02' :: 'Lens'' ('Plucker' a) a+-- 'p03' :: 'Lens'' ('Plucker' a) a+-- 'p23' :: 'Lens'' ('Plucker' a) a+-- 'p31' :: 'Lens'' ('Plucker' a) a+-- 'p12' :: 'Lens'' ('Plucker' a) a -- @-p01, p02, p03, p23, p31, p12 :: Functor f => (a -> f a) -> Plucker a -> f (Plucker a)+p01, p02, p03, p23, p31, p12 :: Lens' (Plucker a) a p01 g (Plucker a b c d e f) = (\a' -> Plucker a' b c d e f) <$> g a p02 g (Plucker a b c d e f) = (\b' -> Plucker a b' c d e f) <$> g b p03 g (Plucker a b c d e f) = (\c' -> Plucker a b c' d e f) <$> g c@@ -233,12 +386,12 @@ -- | These elements form an alternate basis for the Plücker space, or the Grassmanian manifold @Gr(2,V4)@. -- -- @--- 'p10' :: 'Num' a => Lens' ('Plucker' a) a--- 'p20' :: 'Num' a => Lens' ('Plucker' a) a--- 'p30' :: 'Num' a => Lens' ('Plucker' a) a--- 'p32' :: 'Num' a => Lens' ('Plucker' a) a--- 'p13' :: 'Num' a => Lens' ('Plucker' a) a--- 'p21' :: 'Num' a => Lens' ('Plucker' a) a+-- 'p10' :: 'Num' a => 'Lens'' ('Plucker' a) a+-- 'p20' :: 'Num' a => 'Lens'' ('Plucker' a) a+-- 'p30' :: 'Num' a => 'Lens'' ('Plucker' a) a+-- 'p32' :: 'Num' a => 'Lens'' ('Plucker' a) a+-- 'p13' :: 'Num' a => 'Lens'' ('Plucker' a) a+-- 'p21' :: 'Num' a => 'Lens'' ('Plucker' a) a -- @ p10, p20, p30, p32, p13, p21 :: (Functor f, Num a) => (a -> f a) -> Plucker a -> f (Plucker a) p10 = anti p01@@ -257,35 +410,296 @@ anti :: (Functor f, Num a) => ((a -> f a) -> r) -> (a -> f a) -> r anti k f = k (fmap negate . f . negate) +e01, e02, e03, e23, e31, e12 :: E Plucker+e01 = E p01+e02 = E p02+e03 = E p03+e23 = E p23+e31 = E p31+e12 = E p12++instance WithIndex.FunctorWithIndex (E Plucker) Plucker where+ imap f (Plucker a b c d e g) = Plucker (f e01 a) (f e02 b) (f e03 c) (f e23 d) (f e31 e) (f e12 g)+ {-# INLINE imap #-}++instance WithIndex.FoldableWithIndex (E Plucker) Plucker where+ ifoldMap f (Plucker a b c d e g) = f e01 a `mappend` f e02 b `mappend` f e03 c+ `mappend` f e23 d `mappend` f e31 e `mappend` f e12 g+ {-# INLINE ifoldMap #-}++instance WithIndex.TraversableWithIndex (E Plucker) Plucker where+ itraverse f (Plucker a b c d e g) = Plucker <$> f e01 a <*> f e02 b <*> f e03 c+ <*> f e23 d <*> f e31 e <*> f e12 g+ {-# INLINE itraverse #-}++#if !MIN_VERSION_lens(5,0,0)+instance Lens.FunctorWithIndex (E Plucker) Plucker where imap = WithIndex.imap+instance Lens.FoldableWithIndex (E Plucker) Plucker where ifoldMap = WithIndex.ifoldMap+instance Lens.TraversableWithIndex (E Plucker) Plucker where itraverse = WithIndex.itraverse+#endif++type instance Index (Plucker a) = E Plucker+type instance IxValue (Plucker a) = a++instance Ixed (Plucker a) where+ ix i = el i+ {-# INLINE ix #-}++instance Each (Plucker a) (Plucker b) a b where+ each = traverse+ {-# INLINE each #-}++ -- | Valid Plücker coordinates @p@ will have @'squaredError' p '==' 0@ -- -- That said, floating point makes a mockery of this claim, so you may want to use 'nearZero'.-squaredError :: (Eq a, Num a) => Plucker a -> a+squaredError :: Num a => Plucker a -> a squaredError v = v >< v {-# INLINE squaredError #-} -- | This isn't th actual metric because this bilinear form gives rise to an isotropic quadratic space infixl 5 >< (><) :: Num a => Plucker a -> Plucker a -> a-Plucker a b c d e f >< Plucker g h i j k l = a*g+b*h+c*i-d*j-e*k-f*l+Plucker a b c d e f >< Plucker g h i j k l = a*l-b*k+c*j+d*i-e*h+f*g {-# INLINE (><) #-} --- | Checks if the line is near-isotropic (isotropic vectors in this quadratic space represent lines in real 3d space)+-- | Checks if the line is near-isotropic (isotropic vectors in this+-- quadratic space represent lines in real 3d space). isotropic :: Epsilon a => Plucker a -> Bool isotropic a = nearZero (a >< a) {-# INLINE isotropic #-} --- | Checks if the two vectors intersect (or nearly intersect)-intersects :: Epsilon a => Plucker a -> Plucker a -> Bool-intersects a b = nearZero (a >< b)+-- | Checks if two lines intersect (or nearly intersect).+intersects :: (Epsilon a, Ord a) => Plucker a -> Plucker a -> Bool+intersects a b = not (a `parallel` b) && passes a b == Coplanar+-- intersects :: Epsilon a => Plucker a -> Plucker a -> Bool+-- intersects a b = nearZero (a >< b) {-# INLINE intersects #-} -instance Metric Plucker where- dot (Plucker a b c d e f) (Plucker g h i j k l) = a*g+b*h+c*i+d*j+e*k+f*l- {-# INLINE dot #-}+-- | Describe how two lines pass each other.+data LinePass = Coplanar+ -- ^ The lines are coplanar (parallel or intersecting).+ | Clockwise+ -- ^ The lines pass each other clockwise (right-handed+ -- screw)+ | Counterclockwise+ -- ^ The lines pass each other counterclockwise+ -- (left-handed screw).+ deriving (Eq, Show,Generic) -instance Epsilon a => Epsilon (Plucker a) where- nearZero = nearZero . quadrance- {-# INLINE nearZero #-}+-- | Check how two lines pass each other. @passes l1 l2@ describes+-- @l2@ when looking down @l1@.+passes :: (Epsilon a, Ord a) => Plucker a -> Plucker a -> LinePass+passes a b+ | nearZero s = Coplanar+ | s > 0 = Counterclockwise+ | otherwise = Clockwise+ where s = (u1 `dot` v2) + (u2 `dot` v1)+ V2 u1 v1 = toUV a+ V2 u2 v2 = toUV b+{-# INLINE passes #-} +-- | Checks if two lines are parallel.+parallel :: Epsilon a => Plucker a -> Plucker a -> Bool+parallel a b = nearZero $ u1 `cross` u2+ where V2 u1 _ = toUV a+ V2 u2 _ = toUV b+{-# INLINE parallel #-}++-- | Represent a Plücker coordinate as a pair of 3-tuples, typically+-- denoted U and V.+toUV :: Plucker a -> V2 (V3 a)+toUV (Plucker a b c d e f) = V2 (V3 a b c) (V3 d e f)++-- | Checks if two lines coincide in space. In other words, undirected equality.+coincides :: (Epsilon a, Fractional a) => Plucker a -> Plucker a -> Bool+coincides p1 p2 = Foldable.all nearZero $ (s *^ p2) - p1+ where s = maybe 1 getFirst . getOptionCompat . fold $ saveDiv <$> p1 <*> p2+ saveDiv x y | nearZero y = optionCompat Nothing+ | otherwise = optionCompat . Just $ First (x / y)+{-# INLINABLE coincides #-}++-- | Checks if two lines coincide in space, and have the same+-- orientation.+coincides' :: (Epsilon a, Fractional a, Ord a) => Plucker a -> Plucker a -> Bool+coincides' p1 p2 = Foldable.all nearZero ((s *^ p2) - p1) && s > 0+ where s = maybe 1 getFirst . getOptionCompat . fold $ saveDiv <$> p1 <*> p2+ saveDiv x y | nearZero y = optionCompat Nothing+ | otherwise = optionCompat . Just $ First (x / y)+{-# INLINABLE coincides' #-}++-- The coincides and coincides' functions above require the use of a Maybe type+-- with the following Monoid instance:+--+-- instance Semigroup a => Monoid (Maybe a) where ...+--+-- Unfortunately, Maybe has only had such an instance since base-4.11. Prior+-- to that, its Monoid instance had an instance context of Monoid a, which is+-- too strong. To compensate, we use CPP to define an OptionCompat type+-- synonym, which is an alias for Maybe on recent versions of base and an alias+-- for Data.Semigroup.Option on older versions of base. We don't want to use+-- Option on recent versions of base, as it is deprecated.+#if MIN_VERSION_base(4,11,0)+type OptionCompat = Maybe++optionCompat :: Maybe a -> OptionCompat a+optionCompat = id++getOptionCompat :: OptionCompat a -> Maybe a+getOptionCompat = id+#else+type OptionCompat = Option++optionCompat :: Maybe a -> OptionCompat a+optionCompat = Option++getOptionCompat :: OptionCompat a -> Maybe a+getOptionCompat = getOption+#endif++-- | The minimum squared distance of a line from the origin.+quadranceToOrigin :: Fractional a => Plucker a -> a+quadranceToOrigin p = (v `dot` v) / (u `dot` u)+ where V2 u v = toUV p+{-# INLINE quadranceToOrigin #-}++-- | The point where a line is closest to the origin.+closestToOrigin :: Fractional a => Plucker a -> V3 a+closestToOrigin p = normalizePoint $ V4 x y z (u `dot` u)+ where V2 u v = toUV p+ V3 x y z = v `cross` u+{-# INLINE closestToOrigin #-}++-- | Not all 6-dimensional points correspond to a line in 3D. This+-- predicate tests that a Plücker coordinate lies on the Grassmann+-- manifold, and does indeed represent a 3D line.+isLine :: Epsilon a => Plucker a -> Bool+isLine p = nearZero $ u `dot` v+ where V2 u v = toUV p+{-# INLINE isLine #-}+ -- TODO: drag some more stuff out of my thesis++data instance U.Vector (Plucker a) = V_Plucker !Int (U.Vector a)+data instance U.MVector s (Plucker a) = MV_Plucker !Int (U.MVector s a)+instance U.Unbox a => U.Unbox (Plucker a)++instance U.Unbox a => M.MVector U.MVector (Plucker a) where+ basicLength (MV_Plucker n _) = n+ basicUnsafeSlice m n (MV_Plucker _ v) = MV_Plucker n (M.basicUnsafeSlice (6*m) (6*n) v)+ basicOverlaps (MV_Plucker _ v) (MV_Plucker _ u) = M.basicOverlaps v u+ basicUnsafeNew n = liftM (MV_Plucker n) (M.basicUnsafeNew (6*n))+ basicUnsafeRead (MV_Plucker _ a) i =+ do let o = 6*i+ x <- M.basicUnsafeRead a o+ y <- M.basicUnsafeRead a (o+1)+ z <- M.basicUnsafeRead a (o+2)+ w <- M.basicUnsafeRead a (o+3)+ v <- M.basicUnsafeRead a (o+4)+ u <- M.basicUnsafeRead a (o+5)+ return (Plucker x y z w v u)+ basicUnsafeWrite (MV_Plucker _ a) i (Plucker x y z w v u) =+ do let o = 6*i+ M.basicUnsafeWrite a o x+ M.basicUnsafeWrite a (o+1) y+ M.basicUnsafeWrite a (o+2) z+ M.basicUnsafeWrite a (o+3) w+ M.basicUnsafeWrite a (o+4) v+ M.basicUnsafeWrite a (o+5) u+ basicInitialize (MV_Plucker _ v) = M.basicInitialize v++instance U.Unbox a => G.Vector U.Vector (Plucker a) where+ basicUnsafeFreeze (MV_Plucker n v) = liftM ( V_Plucker n) (G.basicUnsafeFreeze v)+ basicUnsafeThaw ( V_Plucker n v) = liftM (MV_Plucker n) (G.basicUnsafeThaw v)+ basicLength ( V_Plucker n _) = n+ basicUnsafeSlice m n (V_Plucker _ v) = V_Plucker n (G.basicUnsafeSlice (6*m) (6*n) v)+ basicUnsafeIndexM (V_Plucker _ a) i =+ do let o = 6*i+ x <- G.basicUnsafeIndexM a o+ y <- G.basicUnsafeIndexM a (o+1)+ z <- G.basicUnsafeIndexM a (o+2)+ w <- G.basicUnsafeIndexM a (o+3)+ v <- G.basicUnsafeIndexM a (o+4)+ u <- G.basicUnsafeIndexM a (o+5)+ return (Plucker x y z w v u)++instance MonadZip Plucker where+ mzipWith = liftA2++instance MonadFix Plucker where+ mfix f = Plucker (let Plucker a _ _ _ _ _ = f a in a)+ (let Plucker _ a _ _ _ _ = f a in a)+ (let Plucker _ _ a _ _ _ = f a in a)+ (let Plucker _ _ _ a _ _ = f a in a)+ (let Plucker _ _ _ _ a _ = f a in a)+ (let Plucker _ _ _ _ _ a = f a in a)++instance NFData a => NFData (Plucker a) where+ rnf (Plucker a b c d e f) = rnf a `seq` rnf b `seq` rnf c+ `seq` rnf d `seq` rnf e `seq` rnf f++instance Serial1 Plucker where+ serializeWith = traverse_+ deserializeWith k = Plucker <$> k <*> k <*> k <*> k <*> k <*> k++instance Serial a => Serial (Plucker a) where+ serialize = serializeWith serialize+ deserialize = deserializeWith deserialize++instance Binary a => Binary (Plucker a) where+ put = serializeWith Binary.put+ get = deserializeWith Binary.get++instance Serialize a => Serialize (Plucker a) where+ put = serializeWith Cereal.put+ get = deserializeWith Cereal.get++instance Eq1 Plucker where+ liftEq k (Plucker a1 b1 c1 d1 e1 f1)+ (Plucker a2 b2 c2 d2 e2 f2)+ = k a1 a2 && k b1 b2 && k c1 c2 && k d1 d2 && k e1 e2 && k f1 f2+instance Ord1 Plucker where+ liftCompare k (Plucker a1 b1 c1 d1 e1 f1)+ (Plucker a2 b2 c2 d2 e2 f2)+ = k a1 a2 `mappend` k b1 b2 `mappend` k c1 c2 `mappend` k d1 d2 `mappend` k e1 e2 `mappend` k f1 f2+instance Read1 Plucker where+ liftReadsPrec k _ z = readParen (z > 10) $ \r ->+ [ (Plucker a b c d e f, r7)+ | ("Plucker",r1) <- lex r+ , (a,r2) <- k 11 r1+ , (b,r3) <- k 11 r2+ , (c,r4) <- k 11 r3+ , (d,r5) <- k 11 r4+ , (e,r6) <- k 11 r5+ , (f,r7) <- k 11 r6+ ]+instance Show1 Plucker where+ liftShowsPrec k _ z (Plucker a b c d e f) = showParen (z > 10) $+ showString "Plucker " . k 11 a . showChar ' ' . k 11 b . showChar ' ' . k 11 c . showChar ' ' . k 11 d . showChar ' ' . k 11 e . showChar ' ' . k 11 f++instance Field1 (Plucker a) (Plucker a) a a where+ _1 f (Plucker x y z u v w) = f x <&> \x' -> Plucker x' y z u v w++instance Field2 (Plucker a) (Plucker a) a a where+ _2 f (Plucker x y z u v w) = f y <&> \y' -> Plucker x y' z u v w++instance Field3 (Plucker a) (Plucker a) a a where+ _3 f (Plucker x y z u v w) = f z <&> \z' -> Plucker x y z' u v w++instance Field4 (Plucker a) (Plucker a) a a where+ _4 f (Plucker x y z u v w) = f u <&> \u' -> Plucker x y z u' v w++instance Field5 (Plucker a) (Plucker a) a a where+ _5 f (Plucker x y z u v w) = f v <&> \v' -> Plucker x y z u v' w++instance Field6 (Plucker a) (Plucker a) a a where+ _6 f (Plucker x y z u v w) = f w <&> \w' -> Plucker x y z u v w'++instance Semigroup a => Semigroup (Plucker a) where+ (<>) = liftA2 (<>)++instance Monoid a => Monoid (Plucker a) where+ mempty = pure mempty+#if !(MIN_VERSION_base(4,11,0))+ mappend = liftA2 mappend+#endif
+ src/Linear/Plucker/Coincides.hs view
@@ -0,0 +1,38 @@+{-# LANGUAGE GADTs #-}+---------------------------------------------------------------------------------+-- |+-- Copyright : (C) 2012-2015 Edward Kmett+-- License : BSD-style (see the file LICENSE)+--+-- Maintainer : Edward Kmett <ekmett@gmail.com>+-- Stability : experimental+-- Portability : non-portable+--+-- Utility for working with Plücker coordinates for lines in 3d homogeneous space.+----------------------------------------------------------------------------------+module Linear.Plucker.Coincides+ ( Coincides(..)+ ) where++import Linear.Epsilon+import Linear.Plucker++-- | When lines are represented as Plücker coordinates, we have the+-- ability to check for both directed and undirected+-- equality. Undirected equality between 'Line's (or a 'Line' and a+-- 'Ray') checks that the two lines coincide in 3D space. Directed+-- equality, between two 'Ray's, checks that two lines coincide in 3D,+-- and have the same direction. To accomodate these two notions of+-- equality, we use an 'Eq' instance on the 'Coincides' data type.+--+-- For example, to check the /directed/ equality between two lines,+-- @p1@ and @p2@, we write, @Ray p1 == Ray p2@.+data Coincides a where+ Line :: (Epsilon a, Fractional a) => Plucker a -> Coincides a+ Ray :: (Epsilon a, Fractional a, Ord a) => Plucker a -> Coincides a++instance Eq (Coincides a) where+ Line a == Line b = coincides a b+ Line a == Ray b = coincides a b+ Ray a == Line b = coincides a b+ Ray a == Ray b = coincides' a b
+ src/Linear/Projection.hs view
@@ -0,0 +1,260 @@+{-# LANGUAGE CPP #-}+---------------------------------------------------------------------------+-- |+-- Copyright : (C) 2015 Edward Kmett+-- License : BSD-style (see the file LICENSE)+--+-- Maintainer : Edward Kmett <ekmett@gmail.com>+-- Stability : experimental+-- Portability : non-portable+--+-- Common projection matrices: e.g. perspective/orthographic transformation+-- matrices.+--+-- Analytically derived inverses are also supplied, because they can be+-- much more accurate in practice than computing them through general+-- purpose means+---------------------------------------------------------------------------+module Linear.Projection+ ( lookAt+ , perspective, inversePerspective+ , infinitePerspective, inverseInfinitePerspective+ , frustum, inverseFrustum+ , ortho, inverseOrtho+ ) where++import Control.Lens hiding (index)+import Linear.V3+import Linear.V4+import Linear.Matrix+import Linear.Epsilon+import Linear.Metric++-- $setup+-- >>> import Linear.Matrix+-- >>> import Linear.V2+-- >>> import Linear.V4++-- | Build a look at view matrix+lookAt+ :: (Epsilon a, Floating a)+ => V3 a -- ^ Eye+ -> V3 a -- ^ Center+ -> V3 a -- ^ Up+ -> M44 a+lookAt eye center up =+ V4 (V4 (xa^._x) (xa^._y) (xa^._z) xd)+ (V4 (ya^._x) (ya^._y) (ya^._z) yd)+ (V4 (-za^._x) (-za^._y) (-za^._z) zd)+ (V4 0 0 0 1)+ where za = normalize $ center - eye+ xa = normalize $ cross za up+ ya = cross xa za+ xd = -dot xa eye+ yd = -dot ya eye+ zd = dot za eye++-- | Build a matrix for a symmetric perspective-view frustum+perspective+ :: Floating a+ => a -- ^ FOV (y direction, in radians)+ -> a -- ^ Aspect ratio+ -> a -- ^ Near plane+ -> a -- ^ Far plane+ -> M44 a+perspective fovy aspect near far =+ V4 (V4 x 0 0 0)+ (V4 0 y 0 0)+ (V4 0 0 z w)+ (V4 0 0 (-1) 0)+ where tanHalfFovy = tan $ fovy / 2+ x = 1 / (aspect * tanHalfFovy)+ y = 1 / tanHalfFovy+ fpn = far + near+ fmn = far - near+ oon = 0.5/near+ oof = 0.5/far+ -- z = 1 / (near/fpn - far/fpn) -- would be better by .5 bits+ z = -fpn/fmn+ w = 1/(oof-oon) -- 13 bits error reduced to 0.17+ -- w = -(2 * far * near) / fmn++#ifdef HERBIE+{-# ANN perspective "NoHerbie" #-}+#endif++-- | Build an inverse perspective matrix+inversePerspective+ :: Floating a+ => a -- ^ FOV (y direction, in radians)+ -> a -- ^ Aspect ratio+ -> a -- ^ Near plane+ -> a -- ^ Far plane+ -> M44 a+inversePerspective fovy aspect near far =+ V4 (V4 a 0 0 0 )+ (V4 0 b 0 0 )+ (V4 0 0 0 (-1))+ (V4 0 0 c d )+ where tanHalfFovy = tan $ fovy / 2+ a = aspect * tanHalfFovy+ b = tanHalfFovy+ c = oon - oof+ d = oon + oof+ oon = 0.5/near+ oof = 0.5/far+++-- | Build a perspective matrix per the classic @glFrustum@ arguments.+frustum+ :: Floating a+ => a -- ^ Left+ -> a -- ^ Right+ -> a -- ^ Bottom+ -> a -- ^ Top+ -> a -- ^ Near+ -> a -- ^ Far+ -> M44 a+frustum l r b t n f =+ V4 (V4 x 0 a 0)+ (V4 0 y e 0)+ (V4 0 0 c d)+ (V4 0 0 (-1) 0)+ where+ rml = r-l+ tmb = t-b+ fmn = f-n+ x = 2*n/rml+ y = 2*n/tmb+ a = (r+l)/rml+ e = (t+b)/tmb+ c = negate (f+n)/fmn+ d = (-2*f*n)/fmn++inverseFrustum+ :: Floating a+ => a -- ^ Left+ -> a -- ^ Right+ -> a -- ^ Bottom+ -> a -- ^ Top+ -> a -- ^ Near+ -> a -- ^ Far+ -> M44 a+inverseFrustum l r b t n f =+ V4 (V4 rx 0 0 ax)+ (V4 0 ry 0 by)+ (V4 0 0 0 (-1))+ (V4 0 0 rd cd)+ where+ hrn = 0.5/n+ hrnf = 0.5/(n*f)+ rx = (r-l)*hrn+ ry = (t-b)*hrn+ ax = (r+l)*hrn+ by = (t+b)*hrn+ cd = (f+n)*hrnf+ rd = (n-f)*hrnf++-- | Build a matrix for a symmetric perspective-view frustum with a far plane at infinite+infinitePerspective+ :: Floating a+ => a -- ^ FOV (y direction, in radians)+ -> a -- ^ Aspect Ratio+ -> a -- ^ Near plane+ -> M44 a+infinitePerspective fovy a n =+ V4 (V4 x 0 0 0)+ (V4 0 y 0 0)+ (V4 0 0 (-1) w)+ (V4 0 0 (-1) 0)+ where+ t = n*tan(fovy/2)+ b = -t+ l = b*a+ r = t*a+ x = (2*n)/(r-l)+ y = (2*n)/(t-b)+ w = -2*n++inverseInfinitePerspective+ :: Floating a+ => a -- ^ FOV (y direction, in radians)+ -> a -- ^ Aspect Ratio+ -> a -- ^ Near plane+ -> M44 a+inverseInfinitePerspective fovy a n =+ V4 (V4 rx 0 0 0)+ (V4 0 ry 0 0)+ (V4 0 0 0 (-1))+ (V4 0 0 rw (-rw))+ where+ t = n*tan(fovy/2)+ b = -t+ l = b*a+ r = t*a+ hrn = 0.5/n+ rx = (r-l)*hrn+ ry = (t-b)*hrn+ rw = -hrn++-- | Build an orthographic perspective matrix from 6 clipping planes.+-- This matrix takes the region delimited by these planes and maps it+-- to normalized device coordinates between [-1,1]+--+-- This call is designed to mimic the parameters to the OpenGL @glOrtho@+-- call, so it has a slightly strange convention: Notably: the near and+-- far planes are negated.+--+-- Consequently:+--+-- @+-- 'ortho' l r b t n f !* 'V4' l b (-n) 1 = 'V4' (-1) (-1) (-1) 1+-- 'ortho' l r b t n f !* 'V4' r t (-f) 1 = 'V4' 1 1 1 1+-- @+--+-- Examples:+--+-- >>> ortho 1 2 3 4 5 6 !* V4 1 3 (-5) 1+-- V4 (-1.0) (-1.0) (-1.0) 1.0+--+-- >>> ortho 1 2 3 4 5 6 !* V4 2 4 (-6) 1+-- V4 1.0 1.0 1.0 1.0+ortho+ :: Fractional a+ => a -- ^ Left+ -> a -- ^ Right+ -> a -- ^ Bottom+ -> a -- ^ Top+ -> a -- ^ Near+ -> a -- ^ Far+ -> M44 a+ortho l r b t n f =+ V4 (V4 (-2*x) 0 0 ((r+l)*x))+ (V4 0 (-2*y) 0 ((t+b)*y))+ (V4 0 0 (2*z) ((f+n)*z))+ (V4 0 0 0 1)+ where x = recip(l-r)+ y = recip(b-t)+ z = recip(n-f)++-- | Build an inverse orthographic perspective matrix from 6 clipping planes+inverseOrtho+ :: Fractional a+ => a -- ^ Left+ -> a -- ^ Right+ -> a -- ^ Bottom+ -> a -- ^ Top+ -> a -- ^ Near+ -> a -- ^ Far+ -> M44 a+inverseOrtho l r b t n f =+ V4 (V4 x 0 0 c)+ (V4 0 y 0 d)+ (V4 0 0 z e)+ (V4 0 0 0 1)+ where x = 0.5*(r-l)+ y = 0.5*(t-b)+ z = 0.5*(n-f)+ c = 0.5*(l+r)+ d = 0.5*(b+t)+ e = -0.5*(n+f)
src/Linear/Quaternion.hs view
@@ -1,8 +1,18 @@-{-# LANGUAGE DeriveDataTypeable, PatternGuards, ScopedTypeVariables #-}+{-# LANGUAGE DeriveDataTypeable #-}+{-# LANGUAGE PatternGuards #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE CPP #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE Trustworthy #-}+{-# LANGUAGE DeriveGeneric #-}+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE DeriveLift #-}+ ----------------------------------------------------------------------------- -- |--- Module : Linear.Quaternion--- Copyright : (C) 2012-2013 Edward Kmett,+-- Copyright : (C) 2012-2015 Edward Kmett -- License : BSD-style (see the file LICENSE) -- -- Maintainer : Edward Kmett <ekmett@gmail.com>@@ -15,6 +25,7 @@ ( Quaternion(..) , Complicated(..) , Hamiltonian(..)+ , ee, ei, ej, ek , slerp , asinq , acosq@@ -29,31 +40,82 @@ ) where import Control.Applicative+import Control.DeepSeq (NFData(rnf))+import Control.Monad (liftM)+import Control.Monad.Fix+import Control.Monad.Zip+import Control.Lens as Lens hiding ((<.>))+import Data.Binary as Binary+import Data.Bytes.Serial import Data.Complex (Complex((:+))) import Data.Data import Data.Distributive-import Data.Traversable import Data.Foldable+import qualified Data.Foldable.WithIndex as WithIndex import Data.Functor.Bind+import Data.Functor.Classes+import Data.Functor.Rep+import qualified Data.Functor.WithIndex as WithIndex+import Data.Hashable+import Data.Hashable.Lifted+#if !(MIN_VERSION_base(4,11,0))+import Data.Semigroup (Semigroup(..))+#endif+import Data.Serialize as Cereal import GHC.Arr (Ix(..)) import qualified Data.Foldable as F-import Data.Monoid+import qualified Data.Traversable.WithIndex as WithIndex+import qualified Data.Vector as V+import qualified Data.Vector.Generic.Mutable as M+import qualified Data.Vector.Generic as G+import qualified Data.Vector.Unboxed.Base as U import Foreign.Ptr (castPtr, plusPtr) import Foreign.Storable (Storable(..))-import Linear.Core+import GHC.Generics (Generic, Generic1)+#if defined(MIN_VERSION_template_haskell)+import Language.Haskell.TH.Syntax (Lift)+#endif import Linear.Epsilon import Linear.Conjugate import Linear.Metric+import Linear.V+import Linear.V2 import Linear.V3+import Linear.V4 import Linear.Vector import Prelude hiding (any)--{-# ANN module "HLint: ignore Reduce duplication" #-}+import System.Random (Random(..), Uniform)+import System.Random.Stateful (UniformRange(..)) -- | Quaternions-data Quaternion a = Quaternion a {-# UNPACK #-}!(V3 a)- deriving (Eq,Ord,Read,Show,Data,Typeable)+data Quaternion a = Quaternion !a {-# UNPACK #-}!(V3 a)+ deriving (Eq,Ord,Read,Show,Data+ ,Generic,Generic1+#if defined(MIN_VERSION_template_haskell)+ ,Lift+#endif+ ) +instance Finite Quaternion where+ type Size Quaternion = 4+ toV (Quaternion a (V3 b c d)) = V (V.fromListN 4 [a, b, c, d])+ fromV (V v) = Quaternion (v V.! 0) (V3 (v V.! 1) (v V.! 2) (v V.! 3))++instance Random a => Random (Quaternion a) where+ random g = case random g of+ (a, g') -> case random g' of+ (b, g'') -> (Quaternion a b, g'')+ randomR (Quaternion a b, Quaternion c d) g = case randomR (a,c) g of+ (e, g') -> case randomR (b,d) g' of+ (f, g'') -> (Quaternion e f, g'')++instance Uniform a => Uniform (Quaternion a) where++instance UniformRange a => UniformRange (Quaternion a) where+ uniformRM (Quaternion a b, Quaternion c d) g = Quaternion+ <$> uniformRM (a, c) g+ <*> uniformRM (b, d) g+ instance Functor Quaternion where fmap f (Quaternion e v) = Quaternion (f e) (fmap f v) {-# INLINE fmap #-}@@ -70,7 +132,13 @@ Quaternion f fv <*> Quaternion a v = Quaternion (f a) (fv <*> v) {-# INLINE (<*>) #-} -instance Additive Quaternion+instance Additive Quaternion where+ zero = pure 0+ {-# INLINE zero #-}+ liftU2 = liftA2+ {-# INLINE liftU2 #-}+ liftI2 = liftA2+ {-# INLINE liftI2 #-} instance Bind Quaternion where Quaternion a (V3 b c d) >>- f = Quaternion a' (V3 b' c' d') where@@ -106,15 +174,49 @@ inRange (l1,u1) i1 && inRange (l2,u2) i2 {-# INLINE inRange #-} -instance Core Quaternion where- core f = Quaternion (f _e) (V3 (f _i) (f _j) (f _k))- {-# INLINE core #-}+instance Representable Quaternion where+ type Rep Quaternion = E Quaternion+ tabulate f = Quaternion (f ee) (V3 (f ei) (f ej) (f ek))+ {-# INLINE tabulate #-}+ index xs (E l) = view l xs+ {-# INLINE index #-} +instance WithIndex.FunctorWithIndex (E Quaternion) Quaternion where+ imap f (Quaternion a (V3 b c d)) = Quaternion (f ee a) $ V3 (f ei b) (f ej c) (f ek d)+ {-# INLINE imap #-}++instance WithIndex.FoldableWithIndex (E Quaternion) Quaternion where+ ifoldMap f (Quaternion a (V3 b c d)) = f ee a `mappend` f ei b `mappend` f ej c `mappend` f ek d+ {-# INLINE ifoldMap #-}++instance WithIndex.TraversableWithIndex (E Quaternion) Quaternion where+ itraverse f (Quaternion a (V3 b c d)) = Quaternion <$> f ee a <*> (V3 <$> f ei b <*> f ej c <*> f ek d)+ {-# INLINE itraverse #-}++#if !MIN_VERSION_lens(5,0,0)+instance Lens.FunctorWithIndex (E Quaternion) Quaternion where imap = WithIndex.imap+instance Lens.FoldableWithIndex (E Quaternion) Quaternion where ifoldMap = WithIndex.ifoldMap+instance Lens.TraversableWithIndex (E Quaternion) Quaternion where itraverse = WithIndex.itraverse+#endif++type instance Index (Quaternion a) = E Quaternion+type instance IxValue (Quaternion a) = a++instance Ixed (Quaternion a) where+ ix i = el i+ {-# INLINE ix #-}++instance Each (Quaternion a) (Quaternion b) a b where+ each = traverse+ {-# INLINE each #-}+ instance Foldable Quaternion where foldMap f (Quaternion e v) = f e `mappend` foldMap f v {-# INLINE foldMap #-} foldr f z (Quaternion e v) = f e (F.foldr f z v) {-# INLINE foldr #-}+ null _ = False+ length _ = 4 instance Traversable Quaternion where traverse f (Quaternion e v) = Quaternion <$> f e <*> traverse f v@@ -167,6 +269,14 @@ ik = isInfinite k {-# INLINE signum #-} +instance Hashable a => Hashable (Quaternion a) where+ hashWithSalt s (Quaternion a b) = s `hashWithSalt` a `hashWithSalt` b+ {-# INLINE hashWithSalt #-}++instance Hashable1 Quaternion where+ liftHashWithSalt h s (Quaternion a b) = liftHashWithSalt h (h s a) b+ {-# INLINE liftHashWithSalt #-}+ qNaN :: RealFloat a => Quaternion a qNaN = Quaternion fNaN (V3 fNaN fNaN fNaN) where fNaN = 0/0 {-# INLINE qNaN #-}@@ -186,7 +296,7 @@ (r0*q3-r1*q2+r2*q1-r3*q0)) ^/ (r0*r0 + r1*r1 + r2*r2 + r3*r3) {-# INLINE (/) #-}- recip q = q ^/ quadrance q+ recip q@(Quaternion e v) = Quaternion e (negate v) ^/ quadrance q {-# INLINE recip #-} fromRational x = Quaternion (fromRational x) 0 {-# INLINE fromRational #-}@@ -197,17 +307,12 @@ -- | A vector space that includes the basis elements '_e' and '_i' class Complicated t where- -- |- -- @- -- '_e' :: Lens' (t a) a- -- @- _e :: Functor f => (a -> f a) -> t a -> f (t a)- -- |- -- @- -- '_i' :: Lens' (t a) a- -- @- _i :: Functor f => (a -> f a) -> t a -> f (t a)+ _e, _i :: Lens' (t a) a +ee, ei :: Complicated t => E t+ee = E _e+ei = E _i+ instance Complicated Complex where _e f (a :+ b) = (:+ b) <$> f a {-# INLINE _e #-}@@ -222,22 +327,13 @@ -- | A vector space that includes the basis elements '_e', '_i', '_j' and '_k' class Complicated t => Hamiltonian t where- -- |- -- @- -- '_j' :: Lens' (t a) a- -- @- _j :: Functor f => (a -> f a) -> t a -> f (t a)- -- |- -- @- -- '_k' :: Lens' (t a) a- -- @- _k :: Functor f => (a -> f a) -> t a -> f (t a)- -- |- -- @- -- '_ijk' :: Lens' (t a) (V3 a)- -- @- _ijk :: Functor f => (V3 a -> f (V3 a)) -> t a -> f (t a)+ _j, _k :: Lens' (t a) a+ _ijk :: Lens' (t a) (V3 a) +ej, ek :: Hamiltonian t => E t+ej = E _j+ek = E _k+ instance Hamiltonian Quaternion where _j f (Quaternion a v) = Quaternion a <$> _y f v {-# INLINE _j #-}@@ -280,6 +376,29 @@ pow q t = exp (t *^ log q) {-# INLINE pow #-} +sqrte2pqiq :: (Floating a, Ord a) => a -> a -> a+sqrte2pqiq e qiq -- = sqrt (e*e + qiq)+ | e < - 1.5097698010472593e153 = -(qiq/e) - e+ | e < 5.582399551122541e57 = sqrt (e*e + qiq) -- direct definition+ | otherwise = (qiq/e) + e+-- {-# SPECIALIZE sqrte2pqiq :: Double -> Double -> Double #-}+-- {-# SPECIALIZE sqrte2pqiq :: Float -> Float -> Float #-}+#ifdef HERBIE+{-# ANN sqrte2pqiq "NoHerbie" #-}+#endif++tanrhs :: (Floating a, Ord a) => a -> a -> a -> a+tanrhs sai ai d -- = cosh ai * (sai / ai) / d -- improved from 6.04 bits of error to 0.19 bits+ | sai < -4.618902267687042e-52 = (sai / d / ai) * cosh ai+ | sai < 1.038530535935153e-39 = (cosh ai * sai) / ai / d+ | otherwise = (sai / d / ai) * cosh ai+-- {-# SPECIALIZE tanrhs :: Double -> Double -> Double -> Double #-}+-- {-# SPECIALIZE tanrhs :: Float -> Float -> Float -> Float #-}+#ifdef HERBIE+{-# ANN tanrhs "NoHerbie" #-}+#endif++ -- ehh.. instance RealFloat a => Floating (Quaternion a) where {-# SPECIALIZE instance Floating (Quaternion Float) #-}@@ -288,18 +407,21 @@ {-# INLINE pi #-} exp q@(Quaternion e v) | qiq == 0 = Quaternion (exp e) v- | ai <- sqrt qiq, ee <- exp e = reimagine (ee * cos ai) (ee * (sin ai / ai)) q+ | ai <- sqrt qiq, exe <- exp e = reimagine (exe * cos ai) (exe * (sin ai / ai)) q where qiq = qi q {-# INLINE exp #-}- log q@(Quaternion e v@(V3 _i j k))+ log q@(Quaternion e v) | qiq == 0 = if e >= 0- then Quaternion (log e) v- else Quaternion (log (negate e)) (V3 pi j k) -- mmm, pi- | ai <- sqrt qiq, m <- sqrt (e*e + qiq) = reimagine (log m) (atan2 m e / ai) q+ then Quaternion (log e) v -- Using v rather than 0 preserves negative zeros+ else Quaternion (negate (log (negate e))) v -- negative scalar: negate quaternion, take log, negate again, preserves negative zeros+ | ai <- sqrt qiq = reimagine (log m) (acos (e / m) / ai) q where qiq = qi q+ m = sqrte2pqiq e qiq {-# INLINE log #-}+ x ** y = exp (y * log x) {-# INLINE (**) #-}+ sqrt q@(Quaternion e v) | m == 0 = q | qiq == 0 = if e > 0@@ -307,38 +429,44 @@ else Quaternion 0 (V3 (sqrt (negate e)) 0 0) | im <- sqrt (0.5*(m-e)) / sqrt qiq = Quaternion (0.5*(m+e)) (v^*im) where qiq = qi q- m = sqrt (e*e + qiq)+ m = sqrte2pqiq e qiq {-# INLINE sqrt #-}+ cos q@(Quaternion e v) | qiq == 0 = Quaternion (cos e) v- | ai <- sqrt qiq = reimagine (cos e * cosh ai) (- sin e * (sinh ai / ai)) q+ | ai <- sqrt qiq = reimagine (cos e * cosh ai) (- sin e / ai / sinh ai) q -- 0.15 bits error where qiq = qi q {-# INLINE cos #-}+ sin q@(Quaternion e v) | qiq == 0 = Quaternion (sin e) v- | ai <- sqrt qiq = reimagine (sin e * cosh ai) (cos e * (sinh ai / ai)) q+ | ai <- sqrt qiq = reimagine (sin e * cosh ai) (cos e * sinh ai / ai) q where qiq = qi q {-# INLINE sin #-}+ tan q@(Quaternion e v) | qiq == 0 = Quaternion (tan e) v | ai <- sqrt qiq, ce <- cos e, sai <- sinh ai, d <- ce*ce + sai*sai =- reimagine (ce * sin e / d) (cosh ai * (sai / ai) / d) q+ reimagine (ce * sin e / d) (tanrhs sai ai d) q where qiq = qi q {-# INLINE tan #-}+ sinh q@(Quaternion e v) | qiq == 0 = Quaternion (sinh e) v- | ai <- sqrt qiq = reimagine (sinh e * cos ai) (cosh e * (sin ai / ai)) q+ | ai <- sqrt qiq = reimagine (sinh e * cos ai) (cosh e * sin ai / ai) q where qiq = qi q {-# INLINE sinh #-}+ cosh q@(Quaternion e v) | qiq == 0 = Quaternion (cosh e) v- | ai <- sqrt qiq = reimagine (cosh e * cos ai) ((sinh e * sin ai) / ai) q+ | ai <- sqrt qiq = reimagine (cosh e * cos ai) (sin ai * (sinh e / ai)) q where qiq = qi q {-# INLINE cosh #-}+ tanh q@(Quaternion e v) | qiq == 0 = Quaternion (tanh e) v | ai <- sqrt qiq, se <- sinh e, cai <- cos ai, d <- se*se + cai*cai =- reimagine ((cosh e * se) / d) ((cai * (sin ai / ai)) / d) q+ reimagine (cosh e * se / d) (tanhrhs cai ai d) q where qiq = qi q {-# INLINE tanh #-} @@ -356,6 +484,15 @@ atanh = cut atanh {-# INLINE atanh #-} +tanhrhs :: (Floating a, Ord a) => a -> a -> a -> a+tanhrhs cai ai d -- = cai * (sin ai / ai) / d+ | d >= -4.2173720203427147e-29 && d < 4.446702369113811e64 = cai / (d * (ai / sin ai))+ | otherwise = cai * (1 / ai / sin ai) / d+-- {-# SPECIALIZE tanhrhs :: Double -> Double -> Double -> Double #-}+-- {-# SPECIALIZE tanhrhs :: Float -> Float -> Float -> Float #-}+#ifdef HERBIE+{-# ANN tanhrhs "NoHerbie" #-}+#endif -- | Helper for calculating with specific branch cuts cut :: RealFloat a => (Complex a -> Complex a) -> Quaternion a -> Quaternion a@@ -424,14 +561,15 @@ {-# INLINE atanhq #-} -- | Spherical linear interpolation between two quaternions.+ slerp :: RealFloat a => Quaternion a -> Quaternion a -> a -> Quaternion a slerp q p t | 1.0 - cosphi < 1e-8 = q- | phi <- acos cosphi, r <- recip (sin phi)- = (sin ((1-t)*phi)*r *^ q ^+^ f (sin (t*phi)*r) *^ p) ^/ sin phi+ | otherwise = ((sin ((1-t)*phi) *^ q) + sin (t*phi) *^ f p) ^/ sin phi where- dqp = dot q p- (cosphi, f) = if dqp < 0 then (-dqp, negate) else (dqp, id)+ dqp = dot q p+ (cosphi, f) = if dqp < 0 then (-dqp, negate) else (dqp, id)+ phi = acos cosphi {-# SPECIALIZE slerp :: Quaternion Float -> Quaternion Float -> Float -> Quaternion Float #-} {-# SPECIALIZE slerp :: Quaternion Double -> Quaternion Double -> Double -> Quaternion Double #-} @@ -449,6 +587,117 @@ -- | @'axisAngle' axis theta@ builds a 'Quaternion' representing a -- rotation of @theta@ radians about @axis@. axisAngle :: (Epsilon a, Floating a) => V3 a -> a -> Quaternion a-axisAngle axis theta = normalize $ Quaternion (cos half) $ sin half *^ axis+axisAngle axis theta = Quaternion (cos half) (sin half *^ normalize axis) where half = theta / 2 {-# INLINE axisAngle #-}++data instance U.Vector (Quaternion a) = V_Quaternion !Int (U.Vector a)+data instance U.MVector s (Quaternion a) = MV_Quaternion !Int (U.MVector s a)+instance U.Unbox a => U.Unbox (Quaternion a)++instance U.Unbox a => M.MVector U.MVector (Quaternion a) where+ basicLength (MV_Quaternion n _) = n+ basicUnsafeSlice m n (MV_Quaternion _ v) = MV_Quaternion n (M.basicUnsafeSlice (4*m) (4*n) v)+ basicOverlaps (MV_Quaternion _ v) (MV_Quaternion _ u) = M.basicOverlaps v u+ basicUnsafeNew n = liftM (MV_Quaternion n) (M.basicUnsafeNew (4*n))+ basicUnsafeRead (MV_Quaternion _ v) i =+ do let o = 4*i+ x <- M.basicUnsafeRead v o+ y <- M.basicUnsafeRead v (o+1)+ z <- M.basicUnsafeRead v (o+2)+ w <- M.basicUnsafeRead v (o+3)+ return (Quaternion x (V3 y z w))+ basicUnsafeWrite (MV_Quaternion _ v) i (Quaternion x (V3 y z w)) =+ do let o = 4*i+ M.basicUnsafeWrite v o x+ M.basicUnsafeWrite v (o+1) y+ M.basicUnsafeWrite v (o+2) z+ M.basicUnsafeWrite v (o+3) w+ basicInitialize (MV_Quaternion _ v) = M.basicInitialize v++instance U.Unbox a => G.Vector U.Vector (Quaternion a) where+ basicUnsafeFreeze (MV_Quaternion n v) = liftM ( V_Quaternion n) (G.basicUnsafeFreeze v)+ basicUnsafeThaw ( V_Quaternion n v) = liftM (MV_Quaternion n) (G.basicUnsafeThaw v)+ basicLength ( V_Quaternion n _) = n+ basicUnsafeSlice m n (V_Quaternion _ v) = V_Quaternion n (G.basicUnsafeSlice (4*m) (4*n) v)+ basicUnsafeIndexM (V_Quaternion _ v) i =+ do let o = 4*i+ x <- G.basicUnsafeIndexM v o+ y <- G.basicUnsafeIndexM v (o+1)+ z <- G.basicUnsafeIndexM v (o+2)+ w <- G.basicUnsafeIndexM v (o+3)+ return (Quaternion x (V3 y z w))++instance MonadZip Quaternion where+ mzipWith = liftA2++instance MonadFix Quaternion where+ mfix f = Quaternion (let Quaternion a _ = f a in a)+ (V3 (let Quaternion _ (V3 a _ _) = f a in a)+ (let Quaternion _ (V3 _ a _) = f a in a)+ (let Quaternion _ (V3 _ _ a) = f a in a))++instance NFData a => NFData (Quaternion a) where+ rnf (Quaternion a b) = rnf a `seq` rnf b++instance Serial1 Quaternion where+ serializeWith f (Quaternion a b) = f a >> serializeWith f b+ deserializeWith f = Quaternion <$> f <*> deserializeWith f++instance Serial a => Serial (Quaternion a) where+ serialize = serializeWith serialize+ deserialize = deserializeWith deserialize++instance Binary a => Binary (Quaternion a) where+ put = serializeWith Binary.put+ get = deserializeWith Binary.get++instance Serialize a => Serialize (Quaternion a) where+ put = serializeWith Cereal.put+ get = deserializeWith Cereal.get++instance Eq1 Quaternion where+ liftEq f (Quaternion a b) (Quaternion c d) = f a c && liftEq f b d+instance Ord1 Quaternion where+ liftCompare f (Quaternion a b) (Quaternion c d) = f a c `mappend` liftCompare f b d+instance Show1 Quaternion where+ liftShowsPrec f g d (Quaternion a b) = showsBinaryWith f (liftShowsPrec f g) "Quaternion" d a b+instance Read1 Quaternion where+ liftReadsPrec f g = readsData $ readsBinaryWith f (liftReadsPrec f g) "Quaternion" Quaternion++instance Field1 (Quaternion a) (Quaternion a) a a where+ _1 f (Quaternion w xyz) = f w <&> \w' -> Quaternion w' xyz++instance Field2 (Quaternion a) (Quaternion a) a a where+ _2 f (Quaternion w (V3 x y z)) = f x <&> \x' -> Quaternion w (V3 x' y z)++instance Field3 (Quaternion a) (Quaternion a) a a where+ _3 f (Quaternion w (V3 x y z)) = f y <&> \y' -> Quaternion w (V3 x y' z)++instance Field4 (Quaternion a) (Quaternion a) a a where+ _4 f (Quaternion w (V3 x y z)) = f z <&> \z' -> Quaternion w (V3 x y z')++instance Semigroup a => Semigroup (Quaternion a) where+ (<>) = liftA2 (<>)++instance Monoid a => Monoid (Quaternion a) where+ mempty = pure mempty+#if !(MIN_VERSION_base(4,11,0))+ mappend = liftA2 mappend+#endif++instance R1 Quaternion where+ _x f (Quaternion w (V3 x y z)) = f x <&> \x' -> Quaternion w (V3 x' y z)++instance R2 Quaternion where+ _y f (Quaternion w (V3 x y z)) = f y <&> \y' -> Quaternion w (V3 x y' z)+ _xy f (Quaternion w (V3 x y z)) = f (V2 x y) <&> \(V2 x' y') -> Quaternion w (V3 x' y' z)++instance R3 Quaternion where+ _z f (Quaternion w (V3 x y z)) = f z <&> \z' -> Quaternion w (V3 x y z')+ _xyz f (Quaternion w xyz) = Quaternion w <$> f xyz++instance R4 Quaternion where+ _w f (Quaternion w xyz) = f w <&> \w' -> Quaternion w' xyz+ _xyzw f (Quaternion w (V3 x y z)) = f (V4 x y z w) <&> \(V4 x' y' z' w') -> Quaternion w' (V3 x' y' z')+
+ src/Linear/Trace.hs view
@@ -0,0 +1,116 @@+{-# LANGUAGE CPP #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE DefaultSignatures #-}+{-# LANGUAGE PolyKinds #-}+{-# LANGUAGE Trustworthy #-}+---------------------------------------------------------------------------+-- |+-- Copyright : (C) 2012-2015 Edward Kmett+-- License : BSD-style (see the file LICENSE)+--+-- Maintainer : Edward Kmett <ekmett@gmail.com>+-- Stability : experimental+-- Portability : non-portable+--+-- Simple matrix operation for low-dimensional primitives.+---------------------------------------------------------------------------+module Linear.Trace+ ( Trace(..)+ , frobenius+ ) where++import Control.Monad as Monad+import Linear.V0+import Linear.V1+import Linear.V2+import Linear.V3+import Linear.V4+import Linear.Plucker+import Linear.Quaternion+import Linear.V+import Linear.Vector+import Data.Complex+import Data.Distributive+import Data.Foldable as Foldable+import Data.Functor.Bind as Bind+import Data.Functor.Compose+import Data.Functor.Product+import Data.Hashable+import Data.HashMap.Lazy+import Data.IntMap (IntMap)+import Data.Map (Map)++-- $setup+-- >>> import Data.Complex+-- >>> import Debug.SimpleReflect.Vars+-- >>> import Linear.V2++class Functor m => Trace m where+ -- | Compute the trace of a matrix+ --+ -- >>> trace (V2 (V2 a b) (V2 c d))+ -- a + d+ trace :: Num a => m (m a) -> a+#ifndef HLINT+ default trace :: (Foldable m, Num a) => m (m a) -> a+ trace = Foldable.sum . diagonal+ {-# INLINE trace #-}+#endif++ -- | Compute the diagonal of a matrix+ --+ -- >>> diagonal (V2 (V2 a b) (V2 c d))+ -- V2 a d+ diagonal :: m (m a) -> m a+#ifndef HLINT+ default diagonal :: Monad m => m (m a) -> m a+ diagonal = Monad.join+ {-# INLINE diagonal #-}+#endif++instance Trace IntMap where+ diagonal = Bind.join+ {-# INLINE diagonal #-}++instance Ord k => Trace (Map k) where+ diagonal = Bind.join+ {-# INLINE diagonal #-}++instance (Eq k, Hashable k) => Trace (HashMap k) where+ diagonal = Bind.join+ {-# INLINE diagonal #-}++instance Dim n => Trace (V n)+instance Trace V0+instance Trace V1+instance Trace V2+instance Trace V3+instance Trace V4+instance Trace Plucker+instance Trace Quaternion++instance Trace Complex where+ trace ((a :+ _) :+ (_ :+ b)) = a + b+ {-# INLINE trace #-}+ diagonal ((a :+ _) :+ (_ :+ b)) = a :+ b+ {-# INLINE diagonal #-}++instance (Trace f, Trace g) => Trace (Product f g) where+ trace (Pair xx yy) = trace (pfst <$> xx) + trace (psnd <$> yy) where+ pfst (Pair x _) = x+ psnd (Pair _ y) = y+ {-# INLINE trace #-}+ diagonal (Pair xx yy) = diagonal (pfst <$> xx) `Pair` diagonal (psnd <$> yy) where+ pfst (Pair x _) = x+ psnd (Pair _ y) = y+ {-# INLINE diagonal #-}++instance (Distributive g, Trace g, Trace f) => Trace (Compose g f) where+ trace = trace . fmap (fmap trace . distribute) . getCompose . fmap getCompose+ {-# INLINE trace #-}+ diagonal = Compose . fmap diagonal . diagonal . fmap distribute . getCompose . fmap getCompose+ {-# INLINE diagonal #-}++-- | Compute the <http://mathworld.wolfram.com/FrobeniusNorm.html Frobenius norm> of a matrix.+frobenius :: (Num a, Foldable f, Additive f, Additive g, Distributive g, Trace g) => f (g a) -> a+frobenius m = trace $ fmap (\ f' -> Foldable.foldl' (^+^) zero $ liftI2 (*^) f' m) (distribute m)
+ src/Linear/V.hs view
@@ -0,0 +1,594 @@+{-# LANGUAGE CPP #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE KindSignatures #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE GeneralizedNewtypeDeriving #-}+{-# LANGUAGE DefaultSignatures #-}+{-# LANGUAGE Rank2Types #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE EmptyDataDecls #-}+{-# LANGUAGE MultiParamTypeClasses, FlexibleContexts, FlexibleInstances, UndecidableInstances #-}+{-# LANGUAGE DeriveDataTypeable #-}+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE PolyKinds #-}+{-# LANGUAGE RoleAnnotations #-}+{-# LANGUAGE Trustworthy #-}+{-# LANGUAGE DeriveGeneric #-}++-----------------------------------------------------------------------------+-- |+-- Copyright : (C) 2012-2015 Edward Kmett+-- License : BSD-style (see the file LICENSE)+--+-- Maintainer : Edward Kmett <ekmett@gmail.com>+-- Stability : experimental+-- Portability : non-portable+--+-- n-D Vectors+----------------------------------------------------------------------------++module Linear.V+ ( V(V,toVector)+#ifdef MIN_VERSION_template_haskell+ , int+#endif+ , dim+ , Dim(..)+ , reifyDim+ , reifyVector+ , reifyDimNat+ , reifyVectorNat+ , fromVector+ , Finite(..)+ , _V, _V'+ ) where++import Control.Applicative+import Control.DeepSeq (NFData)+import Control.Monad+import Control.Monad.Fix+import Control.Monad.Trans.State+import Control.Monad.Zip+import Control.Lens as Lens+import Data.Binary as Binary+import Data.Bytes.Serial+import Data.Complex+import Data.Data+import Data.Distributive+import Data.Foldable as Foldable+import qualified Data.Foldable.WithIndex as WithIndex+import Data.Functor.Bind+import Data.Functor.Classes+import Data.Functor.Rep as Rep+import qualified Data.Functor.WithIndex as WithIndex+import Data.Hashable+import Data.Hashable.Lifted+import Data.Kind+import Data.Reflection as R+import Data.Serialize as Cereal+import qualified Data.Traversable.WithIndex as WithIndex+import qualified Data.Vector as V+import Data.Vector (Vector)+import Data.Vector.Fusion.Util (Box(..))+import qualified Data.Vector.Generic as G+import qualified Data.Vector.Unboxed as U+import qualified Data.Vector.Generic.Mutable as M+import Foreign.Ptr+import Foreign.Storable+import GHC.TypeLits+import GHC.Generics (Generic, Generic1)+#if !(MIN_VERSION_reflection(1,3,0)) && defined(MIN_VERSION_template_haskell)+import Language.Haskell.TH+#endif+import Linear.Epsilon+import Linear.Metric+import Linear.Vector+import Prelude as P+#if !(MIN_VERSION_base(4,11,0))+import Data.Semigroup+#endif+import System.Random (Random(..))+import System.Random.Stateful (Uniform(..), UniformRange(..))++class Dim n where+ reflectDim :: p n -> Int++type role V nominal representational++class Finite v where+ type Size (v :: Type -> Type) :: Nat -- this should allow kind k, for Reifies k Int+ toV :: v a -> V (Size v) a+ default toV :: Foldable v => v a -> V (Size v) a+ toV = V . V.fromList . Foldable.toList+ fromV :: V (Size v) a -> v a++instance Finite Complex where+ type Size Complex = 2+ toV (a :+ b) = V (V.fromListN 2 [a, b])+ fromV (V v) = (v V.! 0) :+ (v V.! 1)++_V :: (Finite u, Finite v) => Iso (V (Size u) a) (V (Size v) b) (u a) (v b)+_V = iso fromV toV++_V' :: Finite v => Iso (V (Size v) a) (V (Size v) b) (v a) (v b)+_V' = iso fromV toV++instance Finite (V (n :: Nat)) where+ type Size (V n) = n+ toV = id+ fromV = id++newtype V n a = V { toVector :: V.Vector a } deriving (Eq,Ord,Show,Read,NFData+ ,Generic,Generic1+ )++dim :: forall n a. Dim n => V n a -> Int+dim _ = reflectDim (Proxy :: Proxy n)+{-# INLINE dim #-}++instance KnownNat n => Dim (n :: Nat) where+ reflectDim = fromInteger . natVal+ {-# INLINE reflectDim #-}++instance (Dim n, Random a) => Random (V n a) where+ random = runState (V <$> V.replicateM (reflectDim (Proxy :: Proxy n)) (state random))+ randomR (V ls,V hs) = runState (V <$> V.zipWithM (\l h -> state $ randomR (l,h)) ls hs)++instance (Dim n, Uniform a) => Uniform (V n a) where+ uniformM g = V <$> V.replicateM (reflectDim (Proxy :: Proxy n)) (uniformM g)++instance (Dim n, UniformRange a) => UniformRange (V n a) where+ uniformRM (V ls, V hs) g = V <$> V.zipWithM (\l h -> uniformRM (l, h) g) ls hs+#if (MIN_VERSION_random(1,3,0))+ isInRange (V ls, V hs) (V xs) = V.and $ V.zipWith3 (\l h x -> isInRange (l, h) x) ls hs xs+#endif++data ReifiedDim (s :: Type)++retagDim :: (Proxy s -> a) -> proxy (ReifiedDim s) -> a+retagDim f _ = f Proxy+{-# INLINE retagDim #-}++instance Reifies s Int => Dim (ReifiedDim s) where+ reflectDim = retagDim reflect+ {-# INLINE reflectDim #-}++reifyDimNat :: Int -> (forall (n :: Nat). KnownNat n => Proxy n -> r) -> r+reifyDimNat i f = R.reifyNat (fromIntegral i) f+{-# INLINE reifyDimNat #-}++reifyVectorNat :: forall a r. Vector a -> (forall (n :: Nat). KnownNat n => V n a -> r) -> r+reifyVectorNat v f = reifyNat (fromIntegral $ V.length v) $ \(Proxy :: Proxy n) -> f (V v :: V n a)+{-# INLINE reifyVectorNat #-}++reifyDim :: Int -> (forall (n :: Type). Dim n => Proxy n -> r) -> r+reifyDim i f = R.reify i (go f) where+ go :: (Proxy (ReifiedDim n) -> a) -> proxy n -> a+ go g _ = g Proxy+{-# INLINE reifyDim #-}++reifyVector :: forall a r. Vector a -> (forall (n :: Type). Dim n => V n a -> r) -> r+reifyVector v f = reifyDim (V.length v) $ \(Proxy :: Proxy n) -> f (V v :: V n a)+{-# INLINE reifyVector #-}++instance Dim n => Dim (V n a) where+ reflectDim _ = reflectDim (Proxy :: Proxy n)+ {-# INLINE reflectDim #-}++instance (Dim n, Semigroup a) => Semigroup (V n a) where+ (<>) = liftA2 (<>)++instance (Dim n, Monoid a) => Monoid (V n a) where+ mempty = pure mempty+#if !(MIN_VERSION_base(4,11,0))+ mappend = liftA2 mappend+#endif++instance Functor (V n) where+ fmap f (V as) = V (fmap f as)+ {-# INLINE fmap #-}++instance WithIndex.FunctorWithIndex Int (V n) where+ imap f (V as) = V (Lens.imap f as)+ {-# INLINE imap #-}++instance Foldable (V n) where+ fold (V as) = fold as+ {-# INLINE fold #-}+ foldMap f (V as) = Foldable.foldMap f as+ {-# INLINE foldMap #-}+ foldr f z (V as) = V.foldr f z as+ {-# INLINE foldr #-}+ foldl f z (V as) = V.foldl f z as+ {-# INLINE foldl #-}+ foldr' f z (V as) = V.foldr' f z as+ {-# INLINE foldr' #-}+ foldl' f z (V as) = V.foldl' f z as+ {-# INLINE foldl' #-}+ foldr1 f (V as) = V.foldr1 f as+ {-# INLINE foldr1 #-}+ foldl1 f (V as) = V.foldl1 f as+ {-# INLINE foldl1 #-}+ length (V as) = V.length as+ {-# INLINE length #-}+ null (V as) = V.null as+ {-# INLINE null #-}+ toList (V as) = V.toList as+ {-# INLINE toList #-}+ elem a (V as) = V.elem a as+ {-# INLINE elem #-}+ maximum (V as) = V.maximum as+ {-# INLINE maximum #-}+ minimum (V as) = V.minimum as+ {-# INLINE minimum #-}+ sum (V as) = V.sum as+ {-# INLINE sum #-}+ product (V as) = V.product as+ {-# INLINE product #-}++instance WithIndex.FoldableWithIndex Int (V n) where+ ifoldMap f (V as) = ifoldMap f as+ {-# INLINE ifoldMap #-}++instance Traversable (V n) where+ traverse f (V as) = V <$> traverse f as+ {-# INLINE traverse #-}++instance WithIndex.TraversableWithIndex Int (V n) where+ itraverse f (V as) = V <$> itraverse f as+ {-# INLINE itraverse #-}++#if !MIN_VERSION_lens(5,0,0)+instance Lens.FunctorWithIndex Int (V n) where imap = WithIndex.imap+instance Lens.FoldableWithIndex Int (V n) where ifoldMap = WithIndex.ifoldMap+instance Lens.TraversableWithIndex Int (V n) where itraverse = WithIndex.itraverse+#endif++instance Apply (V n) where+ V as <.> V bs = V (V.zipWith id as bs)+ {-# INLINE (<.>) #-}++instance Dim n => Applicative (V n) where+ pure = V . V.replicate (reflectDim (Proxy :: Proxy n))+ {-# INLINE pure #-}++ V as <*> V bs = V (V.zipWith id as bs)+ {-# INLINE (<*>) #-}++instance Bind (V n) where+ V as >>- f = V $ V.generate (V.length as) $ \i ->+ toVector (f (as `V.unsafeIndex` i)) `V.unsafeIndex` i+ {-# INLINE (>>-) #-}++instance Dim n => Monad (V n) where+#if !(MIN_VERSION_base(4,11,0))+ return = V . V.replicate (reflectDim (Proxy :: Proxy n))+ {-# INLINE return #-}+#endif+ V as >>= f = V $ V.generate (reflectDim (Proxy :: Proxy n)) $ \i ->+ toVector (f (as `V.unsafeIndex` i)) `V.unsafeIndex` i+ {-# INLINE (>>=) #-}++instance Dim n => Additive (V n) where+ zero = pure 0+ {-# INLINE zero #-}+ liftU2 f (V as) (V bs) = V (V.zipWith f as bs)+ {-# INLINE liftU2 #-}+ liftI2 f (V as) (V bs) = V (V.zipWith f as bs)+ {-# INLINE liftI2 #-}++instance (Dim n, Num a) => Num (V n a) where+ V as + V bs = V $ V.zipWith (+) as bs+ {-# INLINE (+) #-}+ V as - V bs = V $ V.zipWith (-) as bs+ {-# INLINE (-) #-}+ V as * V bs = V $ V.zipWith (*) as bs+ {-# INLINE (*) #-}+ negate = fmap negate+ {-# INLINE negate #-}+ abs = fmap abs+ {-# INLINE abs #-}+ signum = fmap signum+ {-# INLINE signum #-}+ fromInteger = pure . fromInteger+ {-# INLINE fromInteger #-}++instance (Dim n, Fractional a) => Fractional (V n a) where+ recip = fmap recip+ {-# INLINE recip #-}+ V as / V bs = V $ V.zipWith (/) as bs+ {-# INLINE (/) #-}+ fromRational = pure . fromRational+ {-# INLINE fromRational #-}++instance (Dim n, Floating a) => Floating (V n a) where+ pi = pure pi+ {-# INLINE pi #-}+ exp = fmap exp+ {-# INLINE exp #-}+ sqrt = fmap sqrt+ {-# INLINE sqrt #-}+ log = fmap log+ {-# INLINE log #-}+ V as ** V bs = V $ V.zipWith (**) as bs+ {-# INLINE (**) #-}+ logBase (V as) (V bs) = V $ V.zipWith logBase as bs+ {-# INLINE logBase #-}+ sin = fmap sin+ {-# INLINE sin #-}+ tan = fmap tan+ {-# INLINE tan #-}+ cos = fmap cos+ {-# INLINE cos #-}+ asin = fmap asin+ {-# INLINE asin #-}+ atan = fmap atan+ {-# INLINE atan #-}+ acos = fmap acos+ {-# INLINE acos #-}+ sinh = fmap sinh+ {-# INLINE sinh #-}+ tanh = fmap tanh+ {-# INLINE tanh #-}+ cosh = fmap cosh+ {-# INLINE cosh #-}+ asinh = fmap asinh+ {-# INLINE asinh #-}+ atanh = fmap atanh+ {-# INLINE atanh #-}+ acosh = fmap acosh+ {-# INLINE acosh #-}++instance Dim n => Distributive (V n) where+ distribute f = V $ V.generate (reflectDim (Proxy :: Proxy n)) $ \i -> fmap (\(V v) -> V.unsafeIndex v i) f+ {-# INLINE distribute #-}++instance Hashable a => Hashable (V n a) where+ hashWithSalt s0 (V v) =+ V.foldl' (\s a -> s `hashWithSalt` a) s0 v+ `hashWithSalt` V.length v++instance Dim n => Hashable1 (V n) where+ liftHashWithSalt h s0 (V v) =+ V.foldl' (\s a -> h s a) s0 v+ `hashWithSalt` V.length v+ {-# INLINE liftHashWithSalt #-}++instance (Dim n, Storable a) => Storable (V n a) where+ sizeOf _ = reflectDim (Proxy :: Proxy n) * sizeOf (undefined:: a)+ {-# INLINE sizeOf #-}+ alignment _ = alignment (undefined :: a)+ {-# INLINE alignment #-}+ poke ptr (V xs) = Foldable.forM_ [0..reflectDim (Proxy :: Proxy n)-1] $ \i ->+ pokeElemOff ptr' i (V.unsafeIndex xs i)+ where ptr' = castPtr ptr+ {-# INLINE poke #-}+ peek ptr = V <$> V.generateM (reflectDim (Proxy :: Proxy n)) (peekElemOff ptr')+ where ptr' = castPtr ptr+ {-# INLINE peek #-}++instance (Dim n, Epsilon a) => Epsilon (V n a) where+ nearZero = nearZero . quadrance+ {-# INLINE nearZero #-}++instance Dim n => Metric (V n) where+ dot (V a) (V b) = V.sum $ V.zipWith (*) a b+ {-# INLINE dot #-}++-- TODO: instance (Dim n, Ix a) => Ix (V n a)++fromVector :: forall n a. Dim n => Vector a -> Maybe (V n a)+fromVector v+ | V.length v == reflectDim (Proxy :: Proxy n) = Just (V v)+ | otherwise = Nothing++#if !(MIN_VERSION_reflection(1,3,0)) && defined(MIN_VERSION_template_haskell)+data Z -- 0+data D (n :: *) -- 2n+data SD (n :: *) -- 2n+1+data PD (n :: *) -- 2n-1++instance Reifies Z Int where+ reflect _ = 0+ {-# INLINE reflect #-}++retagD :: (Proxy n -> a) -> proxy (D n) -> a+retagD f _ = f Proxy+{-# INLINE retagD #-}++retagSD :: (Proxy n -> a) -> proxy (SD n) -> a+retagSD f _ = f Proxy+{-# INLINE retagSD #-}++retagPD :: (Proxy n -> a) -> proxy (PD n) -> a+retagPD f _ = f Proxy+{-# INLINE retagPD #-}++instance Reifies n Int => Reifies (D n) Int where+ reflect = (\n -> n+n) <$> retagD reflect+ {-# INLINE reflect #-}++instance Reifies n Int => Reifies (SD n) Int where+ reflect = (\n -> n+n+1) <$> retagSD reflect+ {-# INLINE reflect #-}++instance Reifies n Int => Reifies (PD n) Int where+ reflect = (\n -> n+n-1) <$> retagPD reflect+ {-# INLINE reflect #-}++-- | This can be used to generate a template haskell splice for a type level version of a given 'int'.+--+-- This does not use GHC TypeLits, instead it generates a numeric type by hand similar to the ones used+-- in the \"Functional Pearl: Implicit Dimurations\" paper by Oleg Kiselyov and Chung-Chieh Shan.+int :: Int -> TypeQ+int n = case quotRem n 2 of+ (0, 0) -> conT ''Z+ (q,-1) -> conT ''PD `appT` int q+ (q, 0) -> conT ''D `appT` int q+ (q, 1) -> conT ''SD `appT` int q+ _ -> error "ghc is bad at math"+#endif++instance Dim n => Representable (V n) where+ type Rep (V n) = Int+ tabulate = V . V.generate (reflectDim (Proxy :: Proxy n))+ {-# INLINE tabulate #-}+ index (V xs) i = xs V.! i+ {-# INLINE index #-}++type instance Index (V n a) = Int+type instance IxValue (V n a) = a++instance Ixed (V n a) where+ ix i f v@(V as)+ | i < 0 || i >= V.length as = pure v+ | otherwise = vLens i f v+ {-# INLINE ix #-}++instance Dim n => MonadZip (V n) where+ mzip (V as) (V bs) = V $ V.zip as bs+ mzipWith f (V as) (V bs) = V $ V.zipWith f as bs++instance Dim n => MonadFix (V n) where+ mfix f = tabulate $ \r -> let a = Rep.index (f a) r in a++instance Each (V n a) (V n b) a b where+ each = traverse+ {-# INLINE each #-}++instance (Bounded a, Dim n) => Bounded (V n a) where+ minBound = pure minBound+ {-# INLINE minBound #-}+ maxBound = pure maxBound+ {-# INLINE maxBound #-}++vConstr :: Constr+vConstr = mkConstr vDataType "variadic" [] Prefix+{-# NOINLINE vConstr #-}++vDataType :: DataType+vDataType = mkDataType "Linear.V.V" [vConstr]+{-# NOINLINE vDataType #-}++instance (Typeable (V n), Typeable (V n a), Dim n, Data a) => Data (V n a) where+ gfoldl f z (V as) = z (V . V.fromList) `f` V.toList as+ toConstr _ = vConstr+ gunfold k z c = case constrIndex c of+ 1 -> k (z (V . V.fromList))+ _ -> error "gunfold"+ dataTypeOf _ = vDataType+ dataCast1 f = gcast1 f++instance Dim n => Serial1 (V n) where+ serializeWith = traverse_+ deserializeWith f = sequenceA $ pure f++instance (Dim n, Serial a) => Serial (V n a) where+ serialize = traverse_ serialize+ deserialize = sequenceA $ pure deserialize++instance (Dim n, Binary a) => Binary (V n a) where+ put = serializeWith Binary.put+ get = deserializeWith Binary.get++instance (Dim n, Serialize a) => Serialize (V n a) where+ put = serializeWith Cereal.put+ get = deserializeWith Cereal.get++instance Eq1 (V n) where+ liftEq f0 (V as0) (V bs0) = go f0 (V.toList as0) (V.toList bs0) where+ go _ [] [] = True+ go f (a:as) (b:bs) = f a b && go f as bs+ go _ _ _ = False++instance Ord1 (V n) where+ liftCompare f0 (V as0) (V bs0) = go f0 (V.toList as0) (V.toList bs0) where+ go f (a:as) (b:bs) = f a b `mappend` go f as bs+ go _ [] [] = EQ+ go _ _ [] = GT+ go _ [] _ = LT++instance Show1 (V n) where+ liftShowsPrec _ g d (V as) = showParen (d > 10) $ showString "V " . g (V.toList as)++instance Dim n => Read1 (V n) where+ liftReadsPrec _ g d = readParen (d > 10) $ \r ->+ [ (V (V.fromList as), r2)+ | ("V",r1) <- lex r+ , (as, r2) <- g r1+ , P.length as == reflectDim (Proxy :: Proxy n)+ ]++data instance U.Vector (V n a) = V_VN {-# UNPACK #-} !Int !(U.Vector a)+data instance U.MVector s (V n a) = MV_VN {-# UNPACK #-} !Int !(U.MVector s a)+instance (Dim n, U.Unbox a) => U.Unbox (V n a)++instance (Dim n, U.Unbox a) => M.MVector U.MVector (V n a) where+ {-# INLINE basicLength #-}+ {-# INLINE basicUnsafeSlice #-}+ {-# INLINE basicOverlaps #-}+ {-# INLINE basicUnsafeNew #-}+ {-# INLINE basicUnsafeRead #-}+ {-# INLINE basicUnsafeWrite #-}+ basicLength (MV_VN n _) = n+ basicUnsafeSlice m n (MV_VN _ v) = MV_VN n (M.basicUnsafeSlice (d*m) (d*n) v)+ where d = reflectDim (Proxy :: Proxy n)+ basicOverlaps (MV_VN _ v) (MV_VN _ u) = M.basicOverlaps v u+ basicUnsafeNew n = liftM (MV_VN n) (M.basicUnsafeNew (d*n))+ where d = reflectDim (Proxy :: Proxy n)+ basicUnsafeRead (MV_VN _ v) i =+ liftM V $ V.generateM d (\j -> M.basicUnsafeRead v (d*i+j))+ where d = reflectDim (Proxy :: Proxy n)+ basicUnsafeWrite (MV_VN _ v0) i (V vn0) = let d0 = V.length vn0 in go v0 vn0 d0 (d0*i) 0+ where+ go v vn d o j+ | j >= d = return ()+ | otherwise = do+ a <- liftBox $ G.basicUnsafeIndexM vn j+ M.basicUnsafeWrite v o a+ go v vn d (o+1) (j+1)+ basicInitialize (MV_VN _ v) = M.basicInitialize v+ {-# INLINE basicInitialize #-}++liftBox :: Monad m => Box a -> m a+liftBox (Box a) = return a+{-# INLINE liftBox #-}++instance (Dim n, U.Unbox a) => G.Vector U.Vector (V n a) where+ {-# INLINE basicUnsafeFreeze #-}+ {-# INLINE basicUnsafeThaw #-}+ {-# INLINE basicLength #-}+ {-# INLINE basicUnsafeSlice #-}+ {-# INLINE basicUnsafeIndexM #-}+ basicUnsafeFreeze (MV_VN n v) = liftM ( V_VN n) (G.basicUnsafeFreeze v)+ basicUnsafeThaw ( V_VN n v) = liftM (MV_VN n) (G.basicUnsafeThaw v)+ basicLength ( V_VN n _) = n+ basicUnsafeSlice m n (V_VN _ v) = V_VN n (G.basicUnsafeSlice (d*m) (d*n) v)+ where d = reflectDim (Proxy :: Proxy n)+ basicUnsafeIndexM (V_VN _ v) i =+ liftM V $ V.generateM d (\j -> G.basicUnsafeIndexM v (d*i+j))+ where d = reflectDim (Proxy :: Proxy n)++vLens :: Int -> Lens' (V n a) a+vLens i = \f (V v) -> f (v V.! i) <&> \a -> V (v V.// [(i, a)])+{-# INLINE vLens #-}++instance ( 1 <= n) => Field1 (V n a) (V n a) a a where _1 = vLens 0+instance ( 2 <= n) => Field2 (V n a) (V n a) a a where _2 = vLens 1+instance ( 3 <= n) => Field3 (V n a) (V n a) a a where _3 = vLens 2+instance ( 4 <= n) => Field4 (V n a) (V n a) a a where _4 = vLens 3+instance ( 5 <= n) => Field5 (V n a) (V n a) a a where _5 = vLens 4+instance ( 6 <= n) => Field6 (V n a) (V n a) a a where _6 = vLens 5+instance ( 7 <= n) => Field7 (V n a) (V n a) a a where _7 = vLens 6+instance ( 8 <= n) => Field8 (V n a) (V n a) a a where _8 = vLens 7+instance ( 9 <= n) => Field9 (V n a) (V n a) a a where _9 = vLens 8+instance (10 <= n) => Field10 (V n a) (V n a) a a where _10 = vLens 9+instance (11 <= n) => Field11 (V n a) (V n a) a a where _11 = vLens 10+instance (12 <= n) => Field12 (V n a) (V n a) a a where _12 = vLens 11+instance (13 <= n) => Field13 (V n a) (V n a) a a where _13 = vLens 12+instance (14 <= n) => Field14 (V n a) (V n a) a a where _14 = vLens 13+instance (15 <= n) => Field15 (V n a) (V n a) a a where _15 = vLens 14+instance (16 <= n) => Field16 (V n a) (V n a) a a where _16 = vLens 15+instance (17 <= n) => Field17 (V n a) (V n a) a a where _17 = vLens 16+instance (18 <= n) => Field18 (V n a) (V n a) a a where _18 = vLens 17+instance (19 <= n) => Field19 (V n a) (V n a) a a where _19 = vLens 18
src/Linear/V0.hs view
@@ -1,10 +1,17 @@+{-# LANGUAGE CPP #-} {-# LANGUAGE DeriveDataTypeable #-} {-# LANGUAGE TypeFamilies #-} {-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE DeriveGeneric #-}+{-# LANGUAGE Trustworthy #-}+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE DeriveLift #-}+ ----------------------------------------------------------------------------- -- |--- Module : Linear.V0--- Copyright : (C) 2012-2013 Edward Kmett,+-- Copyright : (C) 2012-2015 Edward Kmett -- License : BSD-style (see the file LICENSE) -- -- Maintainer : Edward Kmett <ekmett@gmail.com>@@ -18,22 +25,50 @@ ) where import Control.Applicative+import Control.DeepSeq (NFData(rnf))+import Control.Lens as Lens+import Control.Monad.Fix+import Control.Monad.Zip+import Data.Binary -- binary+import Data.Bytes.Serial -- bytes import Data.Data import Data.Distributive import Data.Foldable+import qualified Data.Foldable.WithIndex as WithIndex+import Data.Functor.Bind+import Data.Functor.Classes+import Data.Functor.Rep+import qualified Data.Functor.WithIndex as WithIndex+import Data.Hashable+import Data.Hashable.Lifted import Data.Ix-import Data.Traversable+#if !(MIN_VERSION_base(4,11,0)) import Data.Semigroup-import Data.Functor.Bind+#endif+import Data.Serialize -- cereal+import qualified Data.Traversable.WithIndex as WithIndex+import qualified Data.Vector as V import Foreign.Storable (Storable(..))-import Linear.Core+import GHC.Generics (Generic, Generic1)+#if defined(MIN_VERSION_template_haskell)+import Language.Haskell.TH.Syntax (Lift)+#endif+import qualified Data.Vector.Generic.Mutable as M+import qualified Data.Vector.Generic as G+import qualified Data.Vector.Unboxed.Base as U import Linear.Metric import Linear.Epsilon import Linear.Vector+import Linear.V+import System.Random (Random(..), Uniform)+import System.Random.Stateful (UniformRange(..)) import Prelude hiding (sum) -- $setup+-- >>> import Control.Applicative -- >>> import Control.Lens+-- >>> import qualified Data.Foldable as F+-- >>> let sum xs = F.sum xs -- | A 0-dimensional vector --@@ -43,8 +78,45 @@ -- >>> V0 + V0 -- V0 ---data V0 a = V0 deriving (Eq,Ord,Show,Read,Ix,Enum,Data,Typeable)+data V0 a = V0 deriving (Eq,Ord,Show,Read,Ix,Enum,Data+ ,Generic,Generic1+#if defined(MIN_VERSION_template_haskell)+ ,Lift+#endif+ ) +instance Finite V0 where+ type Size V0 = 0+ toV _ = V V.empty+ fromV _ = V0++instance Random (V0 a) where+ random g = (V0, g)+ randomR _ g = (V0, g)+ randomRs _ _ = repeat V0+ randoms _ = repeat V0++instance Uniform (V0 a) where++instance UniformRange (V0 a) where+ uniformRM (_, _) _ = pure V0++instance Serial1 V0 where+ serializeWith _ = serialize+ deserializeWith _ = deserialize++instance Serial (V0 a) where+ serialize V0 = return ()+ deserialize = return V0++instance Binary (V0 a) where+ put V0 = return ()+ get = return V0++instance Serialize (V0 a) where+ put V0 = return ()+ get = return V0+ instance Functor V0 where fmap _ V0 = V0 {-# INLINE fmap #-}@@ -54,6 +126,8 @@ instance Foldable V0 where foldMap _ V0 = mempty {-# INLINE foldMap #-}+ null _ = True+ length _ = 0 instance Traversable V0 where traverse _ V0 = pure V0@@ -61,23 +135,40 @@ instance Apply V0 where V0 <.> V0 = V0- {-@ INLINE (<.>) #-}+ {-# INLINE (<.>) #-} instance Applicative V0 where pure _ = V0 {-# INLINE pure #-} V0 <*> V0 = V0- {-@ INLINE (<*>) #-}+ {-# INLINE (<*>) #-} -instance Additive V0+instance Semigroup (V0 a) where+ _ <> _ = V0 +instance Monoid (V0 a) where+ mempty = V0+#if !(MIN_VERSION_base(4,11,0))+ mappend _ _ = V0+#endif++instance Additive V0 where+ zero = V0+ {-# INLINE zero #-}+ liftU2 _ V0 V0 = V0+ {-# INLINE liftU2 #-}+ liftI2 _ V0 V0 = V0+ {-# INLINE liftI2 #-}+ instance Bind V0 where V0 >>- _ = V0 {-# INLINE (>>-) #-} instance Monad V0 where+#if !(MIN_VERSION_base(4,11,0)) return _ = V0 {-# INLINE return #-}+#endif V0 >>= _ = V0 {-# INLINE (>>=) #-} @@ -105,23 +196,67 @@ fromRational _ = V0 {-# INLINE fromRational #-} +instance Floating (V0 a) where+ pi = V0+ {-# INLINE pi #-}+ exp V0 = V0+ {-# INLINE exp #-}+ sqrt V0 = V0+ {-# INLINE sqrt #-}+ log V0 = V0+ {-# INLINE log #-}+ V0 ** V0 = V0+ {-# INLINE (**) #-}+ logBase V0 V0 = V0+ {-# INLINE logBase #-}+ sin V0 = V0+ {-# INLINE sin #-}+ tan V0 = V0+ {-# INLINE tan #-}+ cos V0 = V0+ {-# INLINE cos #-}+ asin V0 = V0+ {-# INLINE asin #-}+ atan V0 = V0+ {-# INLINE atan #-}+ acos V0 = V0+ {-# INLINE acos #-}+ sinh V0 = V0+ {-# INLINE sinh #-}+ tanh V0 = V0+ {-# INLINE tanh #-}+ cosh V0 = V0+ {-# INLINE cosh #-}+ asinh V0 = V0+ {-# INLINE asinh #-}+ atanh V0 = V0+ {-# INLINE atanh #-}+ acosh V0 = V0+ {-# INLINE acosh #-}+ instance Metric V0 where dot V0 V0 = 0 {-# INLINE dot #-} -instance Core V0 where- core _ = V0- {-# INLINE core #-}- instance Distributive V0 where distribute _ = V0 {-# INLINE distribute #-} -instance Epsilon a => Epsilon (V0 a) where+instance Hashable (V0 a) where+ hash V0 = 0+ {-# INLINE hash #-}+ hashWithSalt s V0 = s+ {-# INLINE hashWithSalt #-}++instance Hashable1 V0 where+ liftHashWithSalt _ s V0 = s+ {-# INLINE liftHashWithSalt #-}++instance Epsilon (V0 a) where nearZero _ = True {-# INLINE nearZero #-} -instance Storable a => Storable (V0 a) where+instance Storable (V0 a) where sizeOf _ = 0 {-# INLINE sizeOf #-} alignment _ = 1@@ -130,3 +265,97 @@ {-# INLINE poke #-} peek _ = return V0 {-# INLINE peek #-}++instance WithIndex.FunctorWithIndex (E V0) V0 where+ imap _ V0 = V0+ {-# INLINE imap #-}++instance WithIndex.FoldableWithIndex (E V0) V0 where+ ifoldMap _ V0 = mempty+ {-# INLINE ifoldMap #-}++instance WithIndex.TraversableWithIndex (E V0) V0 where+ itraverse _ V0 = pure V0+ {-# INLINE itraverse #-}++#if !MIN_VERSION_lens(5,0,0)+instance Lens.FunctorWithIndex (E V0) V0 where imap = WithIndex.imap+instance Lens.FoldableWithIndex (E V0) V0 where ifoldMap = WithIndex.ifoldMap+instance Lens.TraversableWithIndex (E V0) V0 where itraverse = WithIndex.itraverse+#endif++instance Representable V0 where+ type Rep V0 = E V0+ tabulate _ = V0+ {-# INLINE tabulate #-}+ index xs (E l) = view l xs+ {-# INLINE index #-}++type instance Index (V0 a) = E V0+type instance IxValue (V0 a) = a++instance Ixed (V0 a) where+ ix i = el i+ {-# INLINE ix #-}++instance Each (V0 a) (V0 b) a b where+ each = traverse+ {-# INLINE each #-}++newtype instance U.Vector (V0 a) = V_V0 Int+newtype instance U.MVector s (V0 a) = MV_V0 Int+instance U.Unbox (V0 a)++instance M.MVector U.MVector (V0 a) where+ {-# INLINE basicLength #-}+ {-# INLINE basicUnsafeSlice #-}+ {-# INLINE basicOverlaps #-}+ {-# INLINE basicUnsafeNew #-}+ {-# INLINE basicUnsafeRead #-}+ {-# INLINE basicUnsafeWrite #-}+ basicLength (MV_V0 n) = n+ basicUnsafeSlice _ n _ = MV_V0 n+ basicOverlaps _ _ = False+ basicUnsafeNew n = return (MV_V0 n)+ basicUnsafeRead _ _ = return V0+ basicUnsafeWrite _ _ _ = return ()+ basicInitialize _ = return ()+ {-# INLINE basicInitialize #-}++instance G.Vector U.Vector (V0 a) where+ {-# INLINE basicUnsafeFreeze #-}+ {-# INLINE basicUnsafeThaw #-}+ {-# INLINE basicLength #-}+ {-# INLINE basicUnsafeSlice #-}+ {-# INLINE basicUnsafeIndexM #-}+ basicUnsafeFreeze (MV_V0 n) = return (V_V0 n)+ basicUnsafeThaw (V_V0 n) = return (MV_V0 n)+ basicLength (V_V0 n) = n+ basicUnsafeSlice _ n _ = V_V0 n+ basicUnsafeIndexM _ _ = return V0++instance MonadZip V0 where+ mzip V0 V0 = V0+ mzipWith _ V0 V0 = V0+ munzip V0 = (V0, V0)++instance MonadFix V0 where+ mfix _ = V0++instance Bounded (V0 a) where+ minBound = V0+ {-# INLINE minBound #-}+ maxBound = V0+ {-# INLINE maxBound #-}++instance NFData (V0 a) where+ rnf V0 = ()++instance Eq1 V0 where+ liftEq _ _ _ = True+instance Ord1 V0 where+ liftCompare _ _ _ = EQ+instance Show1 V0 where+ liftShowsPrec _ _ = showsPrec+instance Read1 V0 where+ liftReadsPrec _ _ = readsPrec
+ src/Linear/V1.hs view
@@ -0,0 +1,400 @@+{-# LANGUAGE CPP #-}+{-# LANGUAGE DeriveDataTypeable #-}+{-# LANGUAGE DeriveFunctor #-}+{-# LANGUAGE DeriveFoldable #-}+{-# LANGUAGE DeriveTraversable #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE GeneralizedNewtypeDeriving #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE Trustworthy #-}+{-# LANGUAGE DeriveGeneric #-}+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE DeriveLift #-}++-----------------------------------------------------------------------------+-- |+-- Copyright : (C) 2012-2015 Edward Kmett+-- License : BSD-style (see the file LICENSE)+--+-- Maintainer : Edward Kmett <ekmett@gmail.com>+-- Stability : experimental+-- Portability : non-portable+--+-- 1-D Vectors+----------------------------------------------------------------------------+module Linear.V1+ ( V1(..)+ , R1(..)+ , ex+ ) where++import Control.Applicative+import Control.DeepSeq (NFData)+import Control.Monad (liftM)+import Control.Monad.Fix+import Control.Monad.Zip+import Control.Lens as Lens+import Data.Binary as Binary+import Data.Bytes.Serial+import Data.Serialize as Cereal+import Data.Data+import Data.Distributive+import Data.Foldable+import qualified Data.Foldable.WithIndex as WithIndex+import Data.Functor.Bind+import Data.Functor.Classes+import Data.Functor.Rep+import qualified Data.Functor.WithIndex as WithIndex+import Data.Hashable+import Data.Hashable.Lifted+import Data.Semigroup.Foldable+import qualified Data.Traversable.WithIndex as WithIndex+import qualified Data.Vector as V+import Linear.V+import Foreign.Storable (Storable)+import GHC.Arr (Ix(..))+import GHC.Generics (Generic, Generic1)+#if defined(MIN_VERSION_template_haskell)+import Language.Haskell.TH.Syntax (Lift)+#endif+import Linear.Metric+import Linear.Epsilon+import Linear.Vector+import Prelude hiding (sum)+import System.Random (Random(..), Uniform)+import System.Random.Stateful (UniformRange(..))+#if !(MIN_VERSION_base(4,11,0))+import Data.Semigroup+#endif++import qualified Data.Vector.Generic.Mutable as M+import qualified Data.Vector.Generic as G+import qualified Data.Vector.Unboxed.Base as U++-- $setup+-- >>> import Control.Applicative+-- >>> import Control.Lens+-- >>> import qualified Data.Foldable as F+-- >>> let sum xs = F.sum xs++-- | A 1-dimensional vector+--+-- >>> pure 1 :: V1 Int+-- V1 1+--+-- >>> V1 2 + V1 3+-- V1 5+--+-- >>> V1 2 * V1 3+-- V1 6+--+-- >>> sum (V1 2)+-- 2++--data V2 a = V2 !a !a deriving (Eq,Ord,Show,Read,Data)+newtype V1 a = V1 a+ deriving (Eq,Ord,Show,Read,Data,+ Functor,Traversable,+ Epsilon,Storable,NFData+ ,Generic,Generic1+#if defined(MIN_VERSION_template_haskell)+ ,Lift+#endif+ )++instance Foldable V1 where+ foldMap f (V1 a) = f a+#if MIN_VERSION_base(4,13,0)+ foldMap' f (V1 a) = f a+#endif+ null _ = False+ length _ = 1++instance Finite V1 where+ type Size V1 = 1+ toV (V1 a) = V (V.singleton a)+ fromV (V v) = V1 (v V.! 0)++instance Foldable1 V1 where+ foldMap1 f (V1 a) = f a+ {-# INLINE foldMap1 #-}++instance Traversable1 V1 where+ traverse1 f (V1 a) = V1 <$> f a+ {-# INLINE traverse1 #-}++instance Apply V1 where+ V1 f <.> V1 x = V1 (f x)+ {-# INLINE (<.>) #-}++instance Applicative V1 where+ pure = V1+ {-# INLINE pure #-}+ V1 f <*> V1 x = V1 (f x)+ {-# INLINE (<*>) #-}++instance Additive V1 where+ zero = pure 0+ {-# INLINE zero #-}+ liftU2 = liftA2+ {-# INLINE liftU2 #-}+ liftI2 = liftA2+ {-# INLINE liftI2 #-}++instance Bind V1 where+ V1 a >>- f = f a+ {-# INLINE (>>-) #-}++instance Monad V1 where+#if !(MIN_VERSION_base(4,11,0))+ return = V1+ {-# INLINE return #-}+#endif+ V1 a >>= f = f a+ {-# INLINE (>>=) #-}++instance Num a => Num (V1 a) where+ (+) = liftA2 (+)+ {-# INLINE (+) #-}+ (-) = liftA2 (-)+ {-# INLINE (-) #-}+ (*) = liftA2 (*)+ {-# INLINE (*) #-}+ negate = fmap negate+ {-# INLINE negate #-}+ abs = fmap abs+ {-# INLINE abs #-}+ signum = fmap signum+ {-# INLINE signum #-}+ fromInteger = pure . fromInteger+ {-# INLINE fromInteger #-}++instance Fractional a => Fractional (V1 a) where+ recip = fmap recip+ {-# INLINE recip #-}+ (/) = liftA2 (/)+ {-# INLINE (/) #-}+ fromRational = pure . fromRational+ {-# INLINE fromRational #-}++instance Floating a => Floating (V1 a) where+ pi = pure pi+ {-# INLINE pi #-}+ exp = fmap exp+ {-# INLINE exp #-}+ sqrt = fmap sqrt+ {-# INLINE sqrt #-}+ log = fmap log+ {-# INLINE log #-}+ (**) = liftA2 (**)+ {-# INLINE (**) #-}+ logBase = liftA2 logBase+ {-# INLINE logBase #-}+ sin = fmap sin+ {-# INLINE sin #-}+ tan = fmap tan+ {-# INLINE tan #-}+ cos = fmap cos+ {-# INLINE cos #-}+ asin = fmap asin+ {-# INLINE asin #-}+ atan = fmap atan+ {-# INLINE atan #-}+ acos = fmap acos+ {-# INLINE acos #-}+ sinh = fmap sinh+ {-# INLINE sinh #-}+ tanh = fmap tanh+ {-# INLINE tanh #-}+ cosh = fmap cosh+ {-# INLINE cosh #-}+ asinh = fmap asinh+ {-# INLINE asinh #-}+ atanh = fmap atanh+ {-# INLINE atanh #-}+ acosh = fmap acosh+ {-# INLINE acosh #-}++instance Hashable a => Hashable (V1 a) where+ hash (V1 a) = hash a+ hashWithSalt s (V1 a) = s `hashWithSalt` a++instance Hashable1 V1 where+ liftHashWithSalt h s (V1 a) = h s a+ {-# INLINE liftHashWithSalt #-}++instance Metric V1 where+ dot (V1 a) (V1 b) = a * b+ {-# INLINE dot #-}++-- | A space that has at least 1 basis vector '_x'.+class R1 t where+ -- |+ -- >>> V1 2 ^._x+ -- 2+ --+ -- >>> V1 2 & _x .~ 3+ -- V1 3+ --+ _x :: Lens' (t a) a++ex :: R1 t => E t+ex = E _x++instance R1 V1 where+ _x f (V1 a) = V1 <$> f a+ {-# INLINE _x #-}++instance R1 Identity where+ _x f (Identity a) = Identity <$> f a+ {-# INLINE _x #-}++instance Distributive V1 where+ distribute f = V1 (fmap (\(V1 x) -> x) f)+ {-# INLINE distribute #-}++instance Ix a => Ix (V1 a) where+ {-# SPECIALISE instance Ix (V1 Int) #-}++ range (V1 l1, V1 u1) =+ [ V1 i1 | i1 <- range (l1,u1) ]+ {-# INLINE range #-}++ unsafeIndex (V1 l1,V1 u1) (V1 i1) = unsafeIndex (l1,u1) i1+ {-# INLINE unsafeIndex #-}++ inRange (V1 l1,V1 u1) (V1 i1) = inRange (l1,u1) i1+ {-# INLINE inRange #-}++instance Representable V1 where+ type Rep V1 = E V1+ tabulate f = V1 (f ex)+ {-# INLINE tabulate #-}+ index xs (E l) = view l xs+ {-# INLINE index #-}++instance WithIndex.FunctorWithIndex (E V1) V1 where+ imap f (V1 a) = V1 (f ex a)+ {-# INLINE imap #-}++instance WithIndex.FoldableWithIndex (E V1) V1 where+ ifoldMap f (V1 a) = f ex a+ {-# INLINE ifoldMap #-}++instance WithIndex.TraversableWithIndex (E V1) V1 where+ itraverse f (V1 a) = V1 <$> f ex a+ {-# INLINE itraverse #-}++#if !MIN_VERSION_lens(5,0,0)+instance Lens.FunctorWithIndex (E V1) V1 where imap = WithIndex.imap+instance Lens.FoldableWithIndex (E V1) V1 where ifoldMap = WithIndex.ifoldMap+instance Lens.TraversableWithIndex (E V1) V1 where itraverse = WithIndex.itraverse+#endif++type instance Index (V1 a) = E V1+type instance IxValue (V1 a) = a++instance Ixed (V1 a) where+ ix i = el i+ {-# INLINE ix #-}++instance Each (V1 a) (V1 b) a b where+ each f (V1 x) = V1 <$> f x+ {-# INLINE each #-}++newtype instance U.Vector (V1 a) = V_V1 (U.Vector a)+newtype instance U.MVector s (V1 a) = MV_V1 (U.MVector s a)+instance U.Unbox a => U.Unbox (V1 a)++instance U.Unbox a => M.MVector U.MVector (V1 a) where+ {-# INLINE basicLength #-}+ {-# INLINE basicUnsafeSlice #-}+ {-# INLINE basicOverlaps #-}+ {-# INLINE basicUnsafeNew #-}+ {-# INLINE basicUnsafeRead #-}+ {-# INLINE basicUnsafeWrite #-}+ basicLength (MV_V1 v) = M.basicLength v+ basicUnsafeSlice m n (MV_V1 v) = MV_V1 (M.basicUnsafeSlice m n v)+ basicOverlaps (MV_V1 v) (MV_V1 u) = M.basicOverlaps v u+ basicUnsafeNew n = liftM MV_V1 (M.basicUnsafeNew n)+ basicUnsafeRead (MV_V1 v) i = liftM V1 (M.basicUnsafeRead v i)+ basicUnsafeWrite (MV_V1 v) i (V1 x) = M.basicUnsafeWrite v i x+ basicInitialize (MV_V1 v) = M.basicInitialize v+ {-# INLINE basicInitialize #-}++instance U.Unbox a => G.Vector U.Vector (V1 a) where+ {-# INLINE basicUnsafeFreeze #-}+ {-# INLINE basicUnsafeThaw #-}+ {-# INLINE basicLength #-}+ {-# INLINE basicUnsafeSlice #-}+ {-# INLINE basicUnsafeIndexM #-}+ basicUnsafeFreeze (MV_V1 v) = liftM V_V1 (G.basicUnsafeFreeze v)+ basicUnsafeThaw (V_V1 v) = liftM MV_V1 (G.basicUnsafeThaw v)+ basicLength (V_V1 v) = G.basicLength v+ basicUnsafeSlice m n (V_V1 v) = V_V1 (G.basicUnsafeSlice m n v)+ basicUnsafeIndexM (V_V1 v) i = liftM V1 (G.basicUnsafeIndexM v i)++instance MonadZip V1 where+ mzip (V1 a) (V1 b) = V1 (a, b)+ mzipWith f (V1 a) (V1 b) = V1 (f a b)+ munzip (V1 (a,b)) = (V1 a, V1 b)++instance MonadFix V1 where+ mfix f = V1 (let V1 a = f a in a)++instance Bounded a => Bounded (V1 a) where+ minBound = pure minBound+ {-# INLINE minBound #-}+ maxBound = pure maxBound+ {-# INLINE maxBound #-}++instance Serial1 V1 where+ serializeWith f (V1 a) = f a+ deserializeWith m = V1 `liftM` m++instance Serial a => Serial (V1 a) where+ serialize (V1 a) = serialize a+ deserialize = V1 `liftM` deserialize++instance Binary a => Binary (V1 a) where+ put = serializeWith Binary.put+ get = deserializeWith Binary.get++instance Serialize a => Serialize (V1 a) where+ put = serializeWith Cereal.put+ get = deserializeWith Cereal.get++instance Random a => Random (V1 a) where+ random g = case random g of (a, g') -> (V1 a, g')+ randoms g = V1 <$> randoms g+ randomR (V1 a, V1 b) g = case randomR (a, b) g of (a', g') -> (V1 a', g')+ randomRs (V1 a, V1 b) g = V1 <$> randomRs (a, b) g++instance Uniform a => Uniform (V1 a) where++instance UniformRange a => UniformRange (V1 a) where+ uniformRM (V1 a, V1 b) g = V1 <$> uniformRM (a, b) g++instance Eq1 V1 where+ liftEq f (V1 a) (V1 b) = f a b+instance Ord1 V1 where+ liftCompare f (V1 a) (V1 b) = f a b+instance Show1 V1 where+ liftShowsPrec f _ d (V1 a) = showParen (d >= 10) $ showString "V1 " . f d a+instance Read1 V1 where+ liftReadsPrec f _ = readsData $ readsUnaryWith f "V1" V1++instance Field1 (V1 a) (V1 b) a b where+ _1 f (V1 x) = V1 <$> f x++instance Semigroup a => Semigroup (V1 a) where+ (<>) = liftA2 (<>)++instance Monoid a => Monoid (V1 a) where+ mempty = pure mempty+#if !(MIN_VERSION_base(4,11,0))+ mappend = liftA2 mappend+#endif+
src/Linear/V2.hs view
@@ -1,11 +1,18 @@ {-# LANGUAGE DeriveDataTypeable #-} {-# LANGUAGE TypeFamilies #-} {-# LANGUAGE ScopedTypeVariables #-}--- {-# OPTIONS_GHC -fno-warn-name-shadowing #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE RankNTypes #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE CPP #-}+{-# LANGUAGE Trustworthy #-}+{-# LANGUAGE DeriveGeneric #-}+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE DeriveLift #-}+ ----------------------------------------------------------------------------- -- |--- Module : Linear.V2--- Copyright : (C) 2012-2013 Edward Kmett,+-- Copyright : (C) 2012-2015 Edward Kmett -- License : BSD-style (see the file LICENSE) -- -- Maintainer : Edward Kmett <ekmett@gmail.com>@@ -16,30 +23,63 @@ ---------------------------------------------------------------------------- module Linear.V2 ( V2(..)+ , R1(..) , R2(..)+ , _yx+ , ex, ey , perp+ , angle+ , unangle+ , crossZ ) where import Control.Applicative+import Control.DeepSeq (NFData(rnf))+import Control.Monad (liftM)+import Control.Monad.Fix+import Control.Monad.Zip+import Control.Lens as Lens hiding ((<.>))+import Data.Binary as Binary+import Data.Bytes.Serial import Data.Data import Data.Distributive import Data.Foldable-import Data.Traversable+import qualified Data.Foldable.WithIndex as WithIndex+import Data.Functor.Bind+import Data.Functor.Classes+import Data.Functor.Rep+import qualified Data.Functor.WithIndex as WithIndex+import Data.Hashable+import Data.Hashable.Lifted import Data.Semigroup import Data.Semigroup.Foldable-import Data.Semigroup.Traversable-import Data.Functor.Bind+import Data.Serialize as Cereal+import qualified Data.Traversable.WithIndex as WithIndex+import qualified Data.Vector as V import Foreign.Ptr (castPtr) import Foreign.Storable (Storable(..)) import GHC.Arr (Ix(..))-import Linear.Core+import GHC.Generics (Generic, Generic1)+#if defined(MIN_VERSION_template_haskell)+import Language.Haskell.TH.Syntax (Lift)+#endif+import qualified Data.Vector.Generic.Mutable as M+import qualified Data.Vector.Generic as G+import qualified Data.Vector.Unboxed.Base as U import Linear.Metric import Linear.Epsilon+import Linear.V import Linear.Vector+import Linear.V1 (R1(..),ex) import Prelude hiding (sum)+import System.Random (Random(..), Uniform)+import System.Random.Stateful (UniformRange(..)) -- $setup+-- >>> import Control.Applicative -- >>> import Control.Lens+-- >>> import qualified Data.Foldable as F+-- >>> let sum xs = F.sum xs -- | A 2-dimensional vector --@@ -55,8 +95,34 @@ -- >>> sum (V2 1 2) -- 3 -data V2 a = V2 a a deriving (Eq,Ord,Show,Read,Data,Typeable)+data V2 a = V2 !a !a deriving+ (Eq,Ord,Show,Read,Data+ ,Generic,Generic1+#if defined(MIN_VERSION_template_haskell)+ ,Lift+#endif+ ) +instance Finite V2 where+ type Size V2 = 2+ toV (V2 a b) = V (V.fromListN 2 [a,b])+ fromV (V v) = V2 (v V.! 0) (v V.! 1)++instance Random a => Random (V2 a) where+ random g = case random g of+ (a, g') -> case random g' of+ (b, g'') -> (V2 a b, g'')+ {-# inline random #-}+ randomR (V2 a b, V2 c d) g = case randomR (a, c) g of+ (x, g') -> case randomR (b, d) g' of+ (y, g'') -> (V2 x y, g'')+ {-# inline randomR #-}++instance Uniform a => Uniform (V2 a) where++instance UniformRange a => UniformRange (V2 a) where+ uniformRM (V2 a b, V2 c d) g = V2 <$> uniformRM (a, c) g <*> uniformRM (b, d) g+ instance Functor V2 where fmap f (V2 a b) = V2 (f a) (f b) {-# INLINE fmap #-}@@ -66,6 +132,12 @@ instance Foldable V2 where foldMap f (V2 a b) = f a `mappend` f b {-# INLINE foldMap #-}+#if MIN_VERSION_base(4,13,0)+ foldMap' f (V2 a b) = f a `mappend` f b+ {-# INLINE foldMap' #-}+#endif+ null _ = False+ length _ = 2 instance Traversable V2 where traverse f (V2 a b) = V2 <$> f a <*> f b@@ -81,16 +153,30 @@ instance Apply V2 where V2 a b <.> V2 d e = V2 (a d) (b e)- {-@ INLINE (<.>) #-}+ {-# INLINE (<.>) #-} instance Applicative V2 where pure a = V2 a a {-# INLINE pure #-} V2 a b <*> V2 d e = V2 (a d) (b e)- {-@ INLINE (<*>) #-}+ {-# INLINE (<*>) #-} -instance Additive V2+instance Hashable a => Hashable (V2 a) where+ hashWithSalt s (V2 a b) = s `hashWithSalt` a `hashWithSalt` b+ {-# INLINE hashWithSalt #-} +instance Hashable1 V2 where+ liftHashWithSalt h s (V2 a b) = s `h` a `h` b+ {-# INLINE liftHashWithSalt #-}++instance Additive V2 where+ zero = pure 0+ {-# INLINE zero #-}+ liftU2 = liftA2+ {-# INLINE liftU2 #-}+ liftI2 = liftA2+ {-# INLINE liftI2 #-}+ instance Bind V2 where V2 a b >>- f = V2 a' b' where V2 a' _ = f a@@ -98,8 +184,10 @@ {-# INLINE (>>-) #-} instance Monad V2 where+#if !(MIN_VERSION_base(4,11,0)) return a = V2 a a {-# INLINE return #-}+#endif V2 a b >>= f = V2 a' b' where V2 a' _ = f a V2 _ b' = f b@@ -129,26 +217,50 @@ fromRational = pure . fromRational {-# INLINE fromRational #-} +instance Floating a => Floating (V2 a) where+ pi = pure pi+ {-# INLINE pi #-}+ exp = fmap exp+ {-# INLINE exp #-}+ sqrt = fmap sqrt+ {-# INLINE sqrt #-}+ log = fmap log+ {-# INLINE log #-}+ (**) = liftA2 (**)+ {-# INLINE (**) #-}+ logBase = liftA2 logBase+ {-# INLINE logBase #-}+ sin = fmap sin+ {-# INLINE sin #-}+ tan = fmap tan+ {-# INLINE tan #-}+ cos = fmap cos+ {-# INLINE cos #-}+ asin = fmap asin+ {-# INLINE asin #-}+ atan = fmap atan+ {-# INLINE atan #-}+ acos = fmap acos+ {-# INLINE acos #-}+ sinh = fmap sinh+ {-# INLINE sinh #-}+ tanh = fmap tanh+ {-# INLINE tanh #-}+ cosh = fmap cosh+ {-# INLINE cosh #-}+ asinh = fmap asinh+ {-# INLINE asinh #-}+ atanh = fmap atanh+ {-# INLINE atanh #-}+ acosh = fmap acosh+ {-# INLINE acosh #-}+ instance Metric V2 where dot (V2 a b) (V2 c d) = a * c + b * d {-# INLINE dot #-} -- | A space that distinguishes 2 orthogonal basis vectors '_x' and '_y', but may have more.-class R2 t where- -- |- -- >>> V2 1 2 ^._x- -- 1- --- -- >>> V2 1 2 & _x .~ 3- -- V2 3 2- --- -- @- -- '_x' :: Lens' (t a) a- -- @- _x :: Functor f => (a -> f a) -> t a -> f (t a)- _x = _xy._x- {-# INLINE _x #-}-+class R1 t => R2 t where -- | -- >>> V2 1 2 ^._y -- 2@@ -156,31 +268,32 @@ -- >>> V2 1 2 & _y .~ 3 -- V2 1 3 --- -- @- -- '_y' :: Lens' (t a) a- -- @- _y :: Functor f => (a -> f a) -> t a -> f (t a)+ _y :: Lens' (t a) a _y = _xy._y {-# INLINE _y #-} - -- |- -- @- -- '_xy' :: Lens' (t a) ('V2' a)- -- @- _xy :: Functor f => (V2 a -> f (V2 a)) -> t a -> f (t a)+ _xy :: Lens' (t a) (V2 a) -instance R2 V2 where+-- |+-- >>> V2 1 2 ^. _yx+-- V2 2 1+_yx :: R2 t => Lens' (t a) (V2 a)+_yx f = _xy $ \(V2 a b) -> f (V2 b a) <&> \(V2 b' a') -> V2 a' b'+{-# INLINE _yx #-}++ey :: R2 t => E t+ey = E _y++instance R1 V2 where _x f (V2 a b) = (`V2` b) <$> f a {-# INLINE _x #-}++instance R2 V2 where _y f (V2 a b) = V2 a <$> f b {-# INLINE _y #-} _xy = id {-# INLINE _xy #-} -instance Core V2 where- core f = V2 (f _x) (f _y)- {-# INLINE core #-}- instance Distributive V2 where distribute f = V2 (fmap (\(V2 x _) -> x) f) (fmap (\(V2 _ y) -> y) f) {-# INLINE distribute #-}@@ -223,3 +336,156 @@ inRange (V2 l1 l2,V2 u1 u2) (V2 i1 i2) = inRange (l1,u1) i1 && inRange (l2,u2) i2 {-# INLINE inRange #-}++instance Representable V2 where+ type Rep V2 = E V2+ tabulate f = V2 (f ex) (f ey)+ {-# INLINE tabulate #-}+ index xs (E l) = view l xs+ {-# INLINE index #-}++instance WithIndex.FunctorWithIndex (E V2) V2 where+ imap f (V2 a b) = V2 (f ex a) (f ey b)+ {-# INLINE imap #-}++instance WithIndex.FoldableWithIndex (E V2) V2 where+ ifoldMap f (V2 a b) = f ex a `mappend` f ey b+ {-# INLINE ifoldMap #-}++instance WithIndex.TraversableWithIndex (E V2) V2 where+ itraverse f (V2 a b) = V2 <$> f ex a <*> f ey b+ {-# INLINE itraverse #-}++#if !MIN_VERSION_lens(5,0,0)+instance Lens.FunctorWithIndex (E V2) V2 where imap = WithIndex.imap+instance Lens.FoldableWithIndex (E V2) V2 where ifoldMap = WithIndex.ifoldMap+instance Lens.TraversableWithIndex (E V2) V2 where itraverse = WithIndex.itraverse+#endif++type instance Index (V2 a) = E V2+type instance IxValue (V2 a) = a++instance Ixed (V2 a) where+ ix i = el i+ {-# INLINE ix #-}++instance Each (V2 a) (V2 b) a b where+ each = traverse+ {-# INLINE each #-}++data instance U.Vector (V2 a) = V_V2 {-# UNPACK #-} !Int !(U.Vector a)+data instance U.MVector s (V2 a) = MV_V2 {-# UNPACK #-} !Int !(U.MVector s a)+instance U.Unbox a => U.Unbox (V2 a)++instance U.Unbox a => M.MVector U.MVector (V2 a) where+ {-# INLINE basicLength #-}+ {-# INLINE basicUnsafeSlice #-}+ {-# INLINE basicOverlaps #-}+ {-# INLINE basicUnsafeNew #-}+ {-# INLINE basicUnsafeRead #-}+ {-# INLINE basicUnsafeWrite #-}+ basicLength (MV_V2 n _) = n+ basicUnsafeSlice m n (MV_V2 _ v) = MV_V2 n (M.basicUnsafeSlice (2*m) (2*n) v)+ basicOverlaps (MV_V2 _ v) (MV_V2 _ u) = M.basicOverlaps v u+ basicUnsafeNew n = liftM (MV_V2 n) (M.basicUnsafeNew (2*n))+ basicUnsafeRead (MV_V2 _ v) i =+ do let o = 2*i+ x <- M.basicUnsafeRead v o+ y <- M.basicUnsafeRead v (o+1)+ return (V2 x y)+ basicUnsafeWrite (MV_V2 _ v) i (V2 x y) =+ do let o = 2*i+ M.basicUnsafeWrite v o x+ M.basicUnsafeWrite v (o+1) y+ basicInitialize (MV_V2 _ v) = M.basicInitialize v+ {-# INLINE basicInitialize #-}++instance U.Unbox a => G.Vector U.Vector (V2 a) where+ {-# INLINE basicUnsafeFreeze #-}+ {-# INLINE basicUnsafeThaw #-}+ {-# INLINE basicLength #-}+ {-# INLINE basicUnsafeSlice #-}+ {-# INLINE basicUnsafeIndexM #-}+ basicUnsafeFreeze (MV_V2 n v) = liftM ( V_V2 n) (G.basicUnsafeFreeze v)+ basicUnsafeThaw ( V_V2 n v) = liftM (MV_V2 n) (G.basicUnsafeThaw v)+ basicLength ( V_V2 n _) = n+ basicUnsafeSlice m n (V_V2 _ v) = V_V2 n (G.basicUnsafeSlice (2*m) (2*n) v)+ basicUnsafeIndexM (V_V2 _ v) i =+ do let o = 2*i+ x <- G.basicUnsafeIndexM v o+ y <- G.basicUnsafeIndexM v (o+1)+ return (V2 x y)++instance MonadZip V2 where+ mzipWith = liftA2++instance MonadFix V2 where+ mfix f = V2 (let V2 a _ = f a in a)+ (let V2 _ a = f a in a)++angle :: Floating a => a -> V2 a+angle a = V2 (cos a) (sin a)++unangle :: (Floating a, Ord a) => V2 a -> a+unangle a@(V2 ax ay) =+ let alpha = asin $ ay / norm a+ in if ax < 0+ then pi - alpha+ else alpha++-- | The Z-component of the cross product of two vectors in the XY-plane.+--+-- >>> crossZ (V2 1 0) (V2 0 1)+-- 1+crossZ :: Num a => V2 a -> V2 a -> a+crossZ (V2 x1 y1) (V2 x2 y2) = x1*y2 - y1*x2+{-# INLINE crossZ #-}++instance Bounded a => Bounded (V2 a) where+ minBound = pure minBound+ {-# INLINE minBound #-}+ maxBound = pure maxBound+ {-# INLINE maxBound #-}++instance NFData a => NFData (V2 a) where+ rnf (V2 a b) = rnf a `seq` rnf b++instance Serial1 V2 where+ serializeWith = traverse_+ deserializeWith k = V2 <$> k <*> k++instance Serial a => Serial (V2 a) where+ serialize = serializeWith serialize+ deserialize = deserializeWith deserialize++instance Binary a => Binary (V2 a) where+ put = serializeWith Binary.put+ get = deserializeWith Binary.get++instance Serialize a => Serialize (V2 a) where+ put = serializeWith Cereal.put+ get = deserializeWith Cereal.get++instance Eq1 V2 where+ liftEq f (V2 a b) (V2 c d) = f a c && f b d+instance Ord1 V2 where+ liftCompare f (V2 a b) (V2 c d) = f a c `mappend` f b d+instance Read1 V2 where+ liftReadsPrec f _ = readsData $ readsBinaryWith f f "V2" V2+instance Show1 V2 where+ liftShowsPrec f _ d (V2 a b) = showsBinaryWith f f "V2" d a b++instance Field1 (V2 a) (V2 a) a a where+ _1 f (V2 x y) = f x <&> \x' -> V2 x' y++instance Field2 (V2 a) (V2 a) a a where+ _2 f (V2 x y) = f y <&> \y' -> V2 x y'++instance Semigroup a => Semigroup (V2 a) where+ (<>) = liftA2 (<>)++instance Monoid a => Monoid (V2 a) where+ mempty = pure mempty+#if !(MIN_VERSION_base(4,11,0))+ mappend = liftA2 mappend+#endif
src/Linear/V3.hs view
@@ -1,8 +1,18 @@-{-# LANGUAGE DeriveDataTypeable, ScopedTypeVariables #-}+{-# LANGUAGE DeriveDataTypeable #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE RankNTypes #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE CPP #-}+{-# LANGUAGE Trustworthy #-}+{-# LANGUAGE DeriveGeneric #-}+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE DeriveLift #-}+ ----------------------------------------------------------------------------- -- |--- Module : Linear.V3--- Copyright : (C) 2012-2013 Edward Kmett,+-- Copyright : (C) 2012-2015 Edward Kmett -- License : BSD-style (see the file LICENSE) -- -- Maintainer : Edward Kmett <ekmett@gmail.com>@@ -14,33 +24,76 @@ module Linear.V3 ( V3(..) , cross, triple+ , R1(..) , R2(..)+ , _yx , R3(..)+ , _xz, _yz, _zx, _zy+ , _xzy, _yxz, _yzx, _zxy, _zyx+ , ex, ey, ez ) where +#if !MIN_VERSION_base(4,18,0) import Control.Applicative+#endif+import Control.DeepSeq (NFData(rnf))+import Control.Monad (liftM)+import Control.Monad.Fix+import Control.Monad.Zip+import Control.Lens as Lens hiding ((<.>))+import Data.Binary as Binary -- binary+import Data.Bytes.Serial -- bytes import Data.Data import Data.Distributive import Data.Foldable+import qualified Data.Foldable.WithIndex as WithIndex import Data.Functor.Bind-import Data.Traversable+import Data.Functor.Classes+import Data.Functor.Rep+import qualified Data.Functor.WithIndex as WithIndex+import Data.Hashable+import Data.Hashable.Lifted+#if !(MIN_VERSION_base(4,11,0)) import Data.Semigroup+#endif import Data.Semigroup.Foldable-import Data.Semigroup.Traversable+import Data.Serialize as Cereal -- cereal+import qualified Data.Traversable.WithIndex as WithIndex+import qualified Data.Vector as V+import qualified Data.Vector.Generic.Mutable as M+import qualified Data.Vector.Generic as G+import qualified Data.Vector.Unboxed.Base as U import Foreign.Ptr (castPtr) import Foreign.Storable (Storable(..)) import GHC.Arr (Ix(..))-import Linear.Core+import GHC.Generics (Generic, Generic1)+#if defined(MIN_VERSION_template_haskell)+import Language.Haskell.TH.Syntax (Lift)+#endif import Linear.Epsilon import Linear.Metric+import Linear.V import Linear.V2 import Linear.Vector+import System.Random (Random(..), Uniform)+import System.Random.Stateful (UniformRange(..)) -{-# ANN module "HLint: ignore Reduce duplication" #-}+-- $setup+-- >>> import Control.Lens hiding (index) -- | A 3-dimensional vector-data V3 a = V3 a a a deriving (Eq,Ord,Show,Read,Data,Typeable)+data V3 a = V3 !a !a !a deriving (Eq,Ord,Show,Read,Data+ ,Generic,Generic1+#if defined(MIN_VERSION_template_haskell)+ ,Lift+#endif+ ) +instance Finite V3 where+ type Size V3 = 3+ toV (V3 a b c) = V (V.fromListN 3 [a,b,c])+ fromV (V v) = V3 (v V.! 0) (v V.! 1) (v V.! 2)+ instance Functor V3 where fmap f (V3 a b c) = V3 (f a) (f b) (f c) {-# INLINE fmap #-}@@ -50,7 +103,31 @@ instance Foldable V3 where foldMap f (V3 a b c) = f a `mappend` f b `mappend` f c {-# INLINE foldMap #-}+#if MIN_VERSION_base(4,13,0)+ foldMap' f (V3 a b c) = (f a `mappend` f b) `mappend` f c+ {-# INLINE foldMap' #-}+#endif+ null _ = False+ length _ = 3 +instance Random a => Random (V3 a) where+ random g = case random g of+ (a, g') -> case random g' of+ (b, g'') -> case random g'' of+ (c, g''') -> (V3 a b c, g''')+ randomR (V3 a b c, V3 a' b' c') g = case randomR (a,a') g of+ (a'', g') -> case randomR (b,b') g' of+ (b'', g'') -> case randomR (c,c') g'' of+ (c'', g''') -> (V3 a'' b'' c'', g''')++instance Uniform a => Uniform (V3 a) where++instance UniformRange a => UniformRange (V3 a) where+ uniformRM (V3 a b c, V3 a' b' c') g = V3+ <$> uniformRM (a, a') g+ <*> uniformRM (b, b') g+ <*> uniformRM (c, c') g+ instance Traversable V3 where traverse f (V3 a b c) = V3 <$> f a <*> f b <*> f c {-# INLINE traverse #-}@@ -73,7 +150,13 @@ V3 a b c <*> V3 d e f = V3 (a d) (b e) (c f) {-# INLINE (<*>) #-} -instance Additive V3+instance Additive V3 where+ zero = pure 0+ {-# INLINE zero #-}+ liftU2 = liftA2+ {-# INLINE liftU2 #-}+ liftI2 = liftA2+ {-# INLINE liftI2 #-} instance Bind V3 where V3 a b c >>- f = V3 a' b' c' where@@ -83,8 +166,10 @@ {-# INLINE (>>-) #-} instance Monad V3 where+#if !(MIN_VERSION_base(4,11,0)) return a = V3 a a a {-# INLINE return #-}+#endif V3 a b c >>= f = V3 a' b' c' where V3 a' _ _ = f a V3 _ b' _ = f b@@ -115,6 +200,52 @@ fromRational = pure . fromRational {-# INLINE fromRational #-} +instance Floating a => Floating (V3 a) where+ pi = pure pi+ {-# INLINE pi #-}+ exp = fmap exp+ {-# INLINE exp #-}+ sqrt = fmap sqrt+ {-# INLINE sqrt #-}+ log = fmap log+ {-# INLINE log #-}+ (**) = liftA2 (**)+ {-# INLINE (**) #-}+ logBase = liftA2 logBase+ {-# INLINE logBase #-}+ sin = fmap sin+ {-# INLINE sin #-}+ tan = fmap tan+ {-# INLINE tan #-}+ cos = fmap cos+ {-# INLINE cos #-}+ asin = fmap asin+ {-# INLINE asin #-}+ atan = fmap atan+ {-# INLINE atan #-}+ acos = fmap acos+ {-# INLINE acos #-}+ sinh = fmap sinh+ {-# INLINE sinh #-}+ tanh = fmap tanh+ {-# INLINE tanh #-}+ cosh = fmap cosh+ {-# INLINE cosh #-}+ asinh = fmap asinh+ {-# INLINE asinh #-}+ atanh = fmap atanh+ {-# INLINE atanh #-}+ acosh = fmap acosh+ {-# INLINE acosh #-}++instance Hashable a => Hashable (V3 a) where+ hashWithSalt s (V3 a b c) = s `hashWithSalt` a `hashWithSalt` b `hashWithSalt` c+ {-# INLINE hashWithSalt #-}++instance Hashable1 V3 where+ liftHashWithSalt h s (V3 a b c) = s `h` a `h` b `h` c+ {-# INLINE liftHashWithSalt #-}+ instance Metric V3 where dot (V3 a b c) (V3 d e f) = a * d + b * e + c * f {-# INLINABLE dot #-}@@ -126,19 +257,51 @@ -- | A space that distinguishes 3 orthogonal basis vectors: '_x', '_y', and '_z'. (It may have more) class R2 t => R3 t where -- |- -- @- -- '_z' :: Lens' (t a) a- -- @- _z :: Functor f => (a -> f a) -> t a -> f (t a)- -- |- -- @- -- '_xyz' :: Lens' (t a) ('V3' a)- -- @- _xyz :: Functor f => (V3 a -> f (V3 a)) -> t a -> f (t a)+ -- >>> V3 1 2 3 ^. _z+ -- 3+ _z :: Lens' (t a) a -instance R2 V3 where+ _xyz :: Lens' (t a) (V3 a)++_xz, _yz, _zx, _zy :: R3 t => Lens' (t a) (V2 a)++_xz f = _xyz $ \(V3 a b c) -> f (V2 a c) <&> \(V2 a' c') -> V3 a' b c'+{-# INLINE _xz #-}++_yz f = _xyz $ \(V3 a b c) -> f (V2 b c) <&> \(V2 b' c') -> V3 a b' c'+{-# INLINE _yz #-}++_zx f = _xyz $ \(V3 a b c) -> f (V2 c a) <&> \(V2 c' a') -> V3 a' b c'+{-# INLINE _zx #-}++_zy f = _xyz $ \(V3 a b c) -> f (V2 c b) <&> \(V2 c' b') -> V3 a b' c'+{-# INLINE _zy #-}++_xzy, _yxz, _yzx, _zxy, _zyx :: R3 t => Lens' (t a) (V3 a)++_xzy f = _xyz $ \(V3 a b c) -> f (V3 a c b) <&> \(V3 a' c' b') -> V3 a' b' c'+{-# INLINE _xzy #-}++_yxz f = _xyz $ \(V3 a b c) -> f (V3 b a c) <&> \(V3 b' a' c') -> V3 a' b' c'+{-# INLINE _yxz #-}++_yzx f = _xyz $ \(V3 a b c) -> f (V3 b c a) <&> \(V3 b' c' a') -> V3 a' b' c'+{-# INLINE _yzx #-}++_zxy f = _xyz $ \(V3 a b c) -> f (V3 c a b) <&> \(V3 c' a' b') -> V3 a' b' c'+{-# INLINE _zxy #-}++_zyx f = _xyz $ \(V3 a b c) -> f (V3 c b a) <&> \(V3 c' b' a') -> V3 a' b' c'+{-# INLINE _zyx #-}++ez :: R3 t => E t+ez = E _z++instance R1 V3 where _x f (V3 a b c) = (\a' -> V3 a' b c) <$> f a {-# INLINE _x #-}++instance R2 V3 where _y f (V3 a b c) = (\b' -> V3 a b' c) <$> f b {-# INLINE _y #-} _xy f (V3 a b c) = (\(V2 a' b') -> V3 a' b' c) <$> f (V2 a b)@@ -150,10 +313,6 @@ _xyz = id {-# INLINE _xyz #-} -instance Core V3 where- core f = V3 (f _x) (f _y) (f _z)- {-# INLINE core #-}- instance Storable a => Storable (V3 a) where sizeOf _ = 3 * sizeOf (undefined::a) {-# INLINE sizeOf #-}@@ -202,3 +361,153 @@ inRange (l1,u1) i1 && inRange (l2,u2) i2 && inRange (l3,u3) i3 {-# INLINE inRange #-}++instance Representable V3 where+ type Rep V3 = E V3+ tabulate f = V3 (f ex) (f ey) (f ez)+ {-# INLINE tabulate #-}+ index xs (E l) = view l xs+ {-# INLINE index #-}++instance WithIndex.FunctorWithIndex (E V3) V3 where+ imap f (V3 a b c) = V3 (f ex a) (f ey b) (f ez c)+ {-# INLINE imap #-}++instance WithIndex.FoldableWithIndex (E V3) V3 where+ ifoldMap f (V3 a b c) = f ex a `mappend` f ey b `mappend` f ez c+ {-# INLINE ifoldMap #-}++instance WithIndex.TraversableWithIndex (E V3) V3 where+ itraverse f (V3 a b c) = V3 <$> f ex a <*> f ey b <*> f ez c+ {-# INLINE itraverse #-}++#if !MIN_VERSION_lens(5,0,0)+instance Lens.FunctorWithIndex (E V3) V3 where imap = WithIndex.imap+instance Lens.FoldableWithIndex (E V3) V3 where ifoldMap = WithIndex.ifoldMap+instance Lens.TraversableWithIndex (E V3) V3 where itraverse = WithIndex.itraverse+#endif++type instance Index (V3 a) = E V3+type instance IxValue (V3 a) = a++instance Ixed (V3 a) where+ ix i = el i+ {-# INLINE ix #-}++instance Each (V3 a) (V3 b) a b where+ each = traverse+ {-# INLINE each #-}++data instance U.Vector (V3 a) = V_V3 {-# UNPACK #-} !Int !(U.Vector a)+data instance U.MVector s (V3 a) = MV_V3 {-# UNPACK #-} !Int !(U.MVector s a)+instance U.Unbox a => U.Unbox (V3 a)++instance U.Unbox a => M.MVector U.MVector (V3 a) where+ {-# INLINE basicLength #-}+ {-# INLINE basicUnsafeSlice #-}+ {-# INLINE basicOverlaps #-}+ {-# INLINE basicUnsafeNew #-}+ {-# INLINE basicUnsafeRead #-}+ {-# INLINE basicUnsafeWrite #-}+ basicLength (MV_V3 n _) = n+ basicUnsafeSlice m n (MV_V3 _ v) = MV_V3 n (M.basicUnsafeSlice (3*m) (3*n) v)+ basicOverlaps (MV_V3 _ v) (MV_V3 _ u) = M.basicOverlaps v u+ basicUnsafeNew n = liftM (MV_V3 n) (M.basicUnsafeNew (3*n))+ basicUnsafeRead (MV_V3 _ v) i =+ do let o = 3*i+ x <- M.basicUnsafeRead v o+ y <- M.basicUnsafeRead v (o+1)+ z <- M.basicUnsafeRead v (o+2)+ return (V3 x y z)+ basicUnsafeWrite (MV_V3 _ v) i (V3 x y z) =+ do let o = 3*i+ M.basicUnsafeWrite v o x+ M.basicUnsafeWrite v (o+1) y+ M.basicUnsafeWrite v (o+2) z+ basicInitialize (MV_V3 _ v) = M.basicInitialize v+ {-# INLINE basicInitialize #-}++instance U.Unbox a => G.Vector U.Vector (V3 a) where+ {-# INLINE basicUnsafeFreeze #-}+ {-# INLINE basicUnsafeThaw #-}+ {-# INLINE basicLength #-}+ {-# INLINE basicUnsafeSlice #-}+ {-# INLINE basicUnsafeIndexM #-}+ basicUnsafeFreeze (MV_V3 n v) = liftM ( V_V3 n) (G.basicUnsafeFreeze v)+ basicUnsafeThaw ( V_V3 n v) = liftM (MV_V3 n) (G.basicUnsafeThaw v)+ basicLength ( V_V3 n _) = n+ basicUnsafeSlice m n (V_V3 _ v) = V_V3 n (G.basicUnsafeSlice (3*m) (3*n) v)+ basicUnsafeIndexM (V_V3 _ v) i =+ do let o = 3*i+ x <- G.basicUnsafeIndexM v o+ y <- G.basicUnsafeIndexM v (o+1)+ z <- G.basicUnsafeIndexM v (o+2)+ return (V3 x y z)++instance MonadZip V3 where+ mzipWith = liftA2++instance MonadFix V3 where+ mfix f = V3 (let V3 a _ _ = f a in a)+ (let V3 _ a _ = f a in a)+ (let V3 _ _ a = f a in a)++instance Bounded a => Bounded (V3 a) where+ minBound = pure minBound+ {-# INLINE minBound #-}+ maxBound = pure maxBound+ {-# INLINE maxBound #-}++instance NFData a => NFData (V3 a) where+ rnf (V3 a b c) = rnf a `seq` rnf b `seq` rnf c++instance Serial1 V3 where+ serializeWith = traverse_+ deserializeWith k = V3 <$> k <*> k <*> k++instance Serial a => Serial (V3 a) where+ serialize = serializeWith serialize+ deserialize = deserializeWith deserialize++instance Binary a => Binary (V3 a) where+ put = serializeWith Binary.put+ get = deserializeWith Binary.get++instance Serialize a => Serialize (V3 a) where+ put = serializeWith Cereal.put+ get = deserializeWith Cereal.get++instance Eq1 V3 where+ liftEq k (V3 a b c) (V3 d e f) = k a d && k b e && k c f+instance Ord1 V3 where+ liftCompare k (V3 a b c) (V3 d e f) = k a d `mappend` k b e `mappend` k c f+instance Read1 V3 where+ liftReadsPrec k _ d = readParen (d > 10) $ \r ->+ [ (V3 a b c, r4)+ | ("V3",r1) <- lex r+ , (a,r2) <- k 11 r1+ , (b,r3) <- k 11 r2+ , (c,r4) <- k 11 r3+ ]+instance Show1 V3 where+ liftShowsPrec f _ d (V3 a b c) = showParen (d > 10) $+ showString "V3 " . f 11 a . showChar ' ' . f 11 b . showChar ' ' . f 11 c++instance Field1 (V3 a) (V3 a) a a where+ _1 f (V3 x y z) = f x <&> \x' -> V3 x' y z++instance Field2 (V3 a) (V3 a) a a where+ _2 f (V3 x y z) = f y <&> \y' -> V3 x y' z++instance Field3 (V3 a) (V3 a) a a where+ _3 f (V3 x y z) = f z <&> \z' -> V3 x y z'++instance Semigroup a => Semigroup (V3 a) where+ (<>) = liftA2 (<>)++instance Monoid a => Monoid (V3 a) where+ mempty = pure mempty+#if !(MIN_VERSION_base(4,11,0))+ mappend = liftA2 mappend+#endif+
src/Linear/V4.hs view
@@ -1,8 +1,18 @@-{-# LANGUAGE DeriveDataTypeable, ScopedTypeVariables #-}+{-# LANGUAGE DeriveDataTypeable #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE RankNTypes #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE CPP #-}+{-# LANGUAGE Trustworthy #-}+{-# LANGUAGE DeriveGeneric #-}+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE DeriveLift #-}+ ----------------------------------------------------------------------------- -- |--- Module : Linear.V4--- Copyright : (C) 2012-2013 Edward Kmett,+-- Copyright : (C) 2012-2015 Edward Kmett -- License : BSD-style (see the file LICENSE) -- -- Maintainer : Edward Kmett <ekmett@gmail.com>@@ -13,36 +23,85 @@ ---------------------------------------------------------------------------- module Linear.V4 ( V4(..)- , vector, point+ , vector, point, normalizePoint+ , R1(..) , R2(..)+ , _yx , R3(..)+ , _xz, _yz, _zx, _zy+ , _xzy, _yxz, _yzx, _zxy, _zyx , R4(..)+ , _xw, _yw, _zw, _wx, _wy, _wz+ , _xyw, _xzw, _xwy, _xwz, _yxw, _yzw, _ywx, _ywz, _zxw, _zyw, _zwx, _zwy+ , _wxy, _wxz, _wyx, _wyz, _wzx, _wzy+ , _xywz, _xzyw, _xzwy, _xwyz, _xwzy, _yxzw , _yxwz, _yzxw, _yzwx, _ywxz+ , _ywzx, _zxyw, _zxwy, _zyxw, _zywx, _zwxy, _zwyx, _wxyz, _wxzy, _wyxz+ , _wyzx, _wzxy, _wzyx+ , ex, ey, ez, ew ) where +#if !MIN_VERSION_base(4,18,0) import Control.Applicative+#endif+import Control.DeepSeq (NFData(rnf))+import Control.Monad (liftM)+import Control.Monad.Fix+import Control.Monad.Zip+import Control.Lens as Lens hiding ((<.>))+import Data.Binary as Binary+import Data.Bytes.Serial import Data.Data import Data.Distributive import Data.Foldable+import qualified Data.Foldable.WithIndex as WithIndex import Data.Functor.Bind+import Data.Functor.Classes+import Data.Functor.Rep+import qualified Data.Functor.WithIndex as WithIndex+import Data.Hashable+import Data.Hashable.Lifted+#if !(MIN_VERSION_base(4,11,0)) import Data.Semigroup+#endif import Data.Semigroup.Foldable-import Data.Semigroup.Traversable-import Data.Traversable+import Data.Serialize as Cereal+import qualified Data.Traversable.WithIndex as WithIndex+import qualified Data.Vector as V+import qualified Data.Vector.Generic.Mutable as M+import qualified Data.Vector.Generic as G+import qualified Data.Vector.Unboxed.Base as U import Foreign.Ptr (castPtr) import Foreign.Storable (Storable(..)) import GHC.Arr (Ix(..))-import Linear.Core+import GHC.Generics (Generic, Generic1)+#if defined(MIN_VERSION_template_haskell)+import Language.Haskell.TH.Syntax (Lift)+#endif import Linear.Epsilon import Linear.Metric+import Linear.V import Linear.V2 import Linear.V3 import Linear.Vector+import System.Random (Random(..), Uniform)+import System.Random.Stateful (UniformRange(..)) -{-# ANN module "HLint: ignore Reduce duplication" #-}+-- $setup+-- >>> import Control.Lens hiding (index) -- | A 4-dimensional vector.-data V4 a = V4 a a a a deriving (Eq,Ord,Show,Read,Data,Typeable)+data V4 a = V4 !a !a !a !a deriving (Eq,Ord,Show,Read,Data+ ,Generic,Generic1+#if defined(MIN_VERSION_template_haskell)+ ,Lift+#endif+ ) +instance Finite V4 where+ type Size V4 = 4+ toV (V4 a b c d) = V (V.fromListN 4 [a,b,c,d])+ fromV (V v) = V4 (v V.! 0) (v V.! 1) (v V.! 2) (v V.! 3)+ instance Functor V4 where fmap f (V4 a b c d) = V4 (f a) (f b) (f c) (f d) {-# INLINE fmap #-}@@ -52,7 +111,34 @@ instance Foldable V4 where foldMap f (V4 a b c d) = f a `mappend` f b `mappend` f c `mappend` f d {-# INLINE foldMap #-}+#if MIN_VERSION_base(4,13,0)+ foldMap' f (V4 a b c d) = ((f a `mappend` f b) `mappend` f c) `mappend` f d+ {-# INLINE foldMap' #-}+#endif+ null _ = False+ length _ = 4 +instance Random a => Random (V4 a) where+ random g = case random g of+ (a, g') -> case random g' of+ (b, g'') -> case random g'' of+ (c, g''') -> case random g''' of+ (d, g'''') -> (V4 a b c d, g'''')+ randomR (V4 a b c d, V4 a' b' c' d') g = case randomR (a,a') g of+ (a'', g') -> case randomR (b,b') g' of+ (b'', g'') -> case randomR (c,c') g'' of+ (c'', g''') -> case randomR (d,d') g''' of+ (d'', g'''') -> (V4 a'' b'' c'' d'', g'''')++instance Uniform a => Uniform (V4 a) where++instance UniformRange a => UniformRange (V4 a) where+ uniformRM (V4 a b c d, V4 a' b' c' d') g = V4+ <$> uniformRM (a, a') g+ <*> uniformRM (b, b') g+ <*> uniformRM (c, c') g+ <*> uniformRM (d, d') g+ instance Traversable V4 where traverse f (V4 a b c d) = V4 <$> f a <*> f b <*> f c <*> f d {-# INLINE traverse #-}@@ -75,7 +161,13 @@ V4 a b c d <.> V4 e f g h = V4 (a e) (b f) (c g) (d h) {-# INLINE (<.>) #-} -instance Additive V4+instance Additive V4 where+ zero = pure 0+ {-# INLINE zero #-}+ liftU2 = liftA2+ {-# INLINE liftU2 #-}+ liftI2 = liftA2+ {-# INLINE liftI2 #-} instance Bind V4 where V4 a b c d >>- f = V4 a' b' c' d' where@@ -86,8 +178,10 @@ {-# INLINE (>>-) #-} instance Monad V4 where+#if !(MIN_VERSION_base(4,11,0)) return a = V4 a a a a {-# INLINE return #-}+#endif V4 a b c d >>= f = V4 a' b' c' d' where V4 a' _ _ _ = f a V4 _ b' _ _ = f b@@ -119,6 +213,44 @@ fromRational = pure . fromRational {-# INLINE fromRational #-} +instance Floating a => Floating (V4 a) where+ pi = pure pi+ {-# INLINE pi #-}+ exp = fmap exp+ {-# INLINE exp #-}+ sqrt = fmap sqrt+ {-# INLINE sqrt #-}+ log = fmap log+ {-# INLINE log #-}+ (**) = liftA2 (**)+ {-# INLINE (**) #-}+ logBase = liftA2 logBase+ {-# INLINE logBase #-}+ sin = fmap sin+ {-# INLINE sin #-}+ tan = fmap tan+ {-# INLINE tan #-}+ cos = fmap cos+ {-# INLINE cos #-}+ asin = fmap asin+ {-# INLINE asin #-}+ atan = fmap atan+ {-# INLINE atan #-}+ acos = fmap acos+ {-# INLINE acos #-}+ sinh = fmap sinh+ {-# INLINE sinh #-}+ tanh = fmap tanh+ {-# INLINE tanh #-}+ cosh = fmap cosh+ {-# INLINE cosh #-}+ asinh = fmap asinh+ {-# INLINE asinh #-}+ atanh = fmap atanh+ {-# INLINE atanh #-}+ acosh = fmap acosh+ {-# INLINE acosh #-}+ instance Metric V4 where dot (V4 a b c d) (V4 e f g h) = a * e + b * f + c * g + d * h {-# INLINE dot #-}@@ -130,22 +262,176 @@ (fmap (\(V4 _ _ _ w) -> w) f) {-# INLINE distribute #-} +instance Hashable a => Hashable (V4 a) where+ hashWithSalt s (V4 a b c d) = s `hashWithSalt` a `hashWithSalt` b `hashWithSalt` c `hashWithSalt` d+ {-# INLINE hashWithSalt #-}++instance Hashable1 V4 where+ liftHashWithSalt h s (V4 a b c d) = s `h` a `h` b `h` c `h` d+ {-# INLINE liftHashWithSalt #-}+ -- | A space that distinguishes orthogonal basis vectors '_x', '_y', '_z', '_w'. (It may have more.) class R3 t => R4 t where -- |- -- @- -- '_w' :: Lens' (t a) a- -- @- _w :: Functor f => (a -> f a) -> t a -> f (t a)- -- |- -- @- -- '_xyzw' :: Lens' (t a) ('V4' a)- -- @- _xyzw :: Functor f => (V4 a -> f (V4 a)) -> t a -> f (t a)+ -- >>> V4 1 2 3 4 ^._w+ -- 4+ _w :: Lens' (t a) a+ _xyzw :: Lens' (t a) (V4 a) -instance R2 V4 where+_xw, _yw, _zw, _wx, _wy, _wz :: R4 t => Lens' (t a) (V2 a)+_xw f = _xyzw $ \(V4 a b c d) -> f (V2 a d) <&> \(V2 a' d') -> V4 a' b c d'+{-# INLINE _xw #-}++_yw f = _xyzw $ \(V4 a b c d) -> f (V2 b d) <&> \(V2 b' d') -> V4 a b' c d'+{-# INLINE _yw #-}++_zw f = _xyzw $ \(V4 a b c d) -> f (V2 c d) <&> \(V2 c' d') -> V4 a b c' d'+{-# INLINE _zw #-}++_wx f = _xyzw $ \(V4 a b c d) -> f (V2 d a) <&> \(V2 d' a') -> V4 a' b c d'+{-# INLINE _wx #-}++_wy f = _xyzw $ \(V4 a b c d) -> f (V2 d b) <&> \(V2 d' b') -> V4 a b' c d'+{-# INLINE _wy #-}++_wz f = _xyzw $ \(V4 a b c d) -> f (V2 d c) <&> \(V2 d' c') -> V4 a b c' d'+{-# INLINE _wz #-}++_xyw, _xzw, _xwy, _xwz, _yxw, _yzw, _ywx, _ywz, _zxw, _zyw, _zwx, _zwy, _wxy, _wxz, _wyx, _wyz, _wzx, _wzy :: R4 t => Lens' (t a) (V3 a)+_xyw f = _xyzw $ \(V4 a b c d) -> f (V3 a b d) <&> \(V3 a' b' d') -> V4 a' b' c d'+{-# INLINE _xyw #-}++_xzw f = _xyzw $ \(V4 a b c d) -> f (V3 a c d) <&> \(V3 a' c' d') -> V4 a' b c' d'+{-# INLINE _xzw #-}++_xwy f = _xyzw $ \(V4 a b c d) -> f (V3 a d b) <&> \(V3 a' d' b') -> V4 a' b' c d'+{-# INLINE _xwy #-}++_xwz f = _xyzw $ \(V4 a b c d) -> f (V3 a d c) <&> \(V3 a' d' c') -> V4 a' b c' d'+{-# INLINE _xwz #-}++_yxw f = _xyzw $ \(V4 a b c d) -> f (V3 b a d) <&> \(V3 b' a' d') -> V4 a' b' c d'+{-# INLINE _yxw #-}++_yzw f = _xyzw $ \(V4 a b c d) -> f (V3 b c d) <&> \(V3 b' c' d') -> V4 a b' c' d'+{-# INLINE _yzw #-}++_ywx f = _xyzw $ \(V4 a b c d) -> f (V3 b d a) <&> \(V3 b' d' a') -> V4 a' b' c d'+{-# INLINE _ywx #-}++_ywz f = _xyzw $ \(V4 a b c d) -> f (V3 b d c) <&> \(V3 b' d' c') -> V4 a b' c' d'+{-# INLINE _ywz #-}++_zxw f = _xyzw $ \(V4 a b c d) -> f (V3 c a d) <&> \(V3 c' a' d') -> V4 a' b c' d'+{-# INLINE _zxw #-}++_zyw f = _xyzw $ \(V4 a b c d) -> f (V3 c b d) <&> \(V3 c' b' d') -> V4 a b' c' d'+{-# INLINE _zyw #-}++_zwx f = _xyzw $ \(V4 a b c d) -> f (V3 c d a) <&> \(V3 c' d' a') -> V4 a' b c' d'+{-# INLINE _zwx #-}++_zwy f = _xyzw $ \(V4 a b c d) -> f (V3 c d b) <&> \(V3 c' d' b') -> V4 a b' c' d'+{-# INLINE _zwy #-}++_wxy f = _xyzw $ \(V4 a b c d) -> f (V3 d a b) <&> \(V3 d' a' b') -> V4 a' b' c d'+{-# INLINE _wxy #-}++_wxz f = _xyzw $ \(V4 a b c d) -> f (V3 d a c) <&> \(V3 d' a' c') -> V4 a' b c' d'+{-# INLINE _wxz #-}++_wyx f = _xyzw $ \(V4 a b c d) -> f (V3 d b a) <&> \(V3 d' b' a') -> V4 a' b' c d'+{-# INLINE _wyx #-}++_wyz f = _xyzw $ \(V4 a b c d) -> f (V3 d b c) <&> \(V3 d' b' c') -> V4 a b' c' d'+{-# INLINE _wyz #-}++_wzx f = _xyzw $ \(V4 a b c d) -> f (V3 d c a) <&> \(V3 d' c' a') -> V4 a' b c' d'+{-# INLINE _wzx #-}++_wzy f = _xyzw $ \(V4 a b c d) -> f (V3 d c b) <&> \(V3 d' c' b') -> V4 a b' c' d'+{-# INLINE _wzy #-}++_xywz, _xzyw, _xzwy, _xwyz, _xwzy, _yxzw , _yxwz, _yzxw, _yzwx, _ywxz+ , _ywzx, _zxyw, _zxwy, _zyxw, _zywx, _zwxy, _zwyx, _wxyz, _wxzy, _wyxz+ , _wyzx, _wzxy, _wzyx :: R4 t => Lens' (t a) (V4 a)+_xywz f = _xyzw $ \(V4 a b c d) -> f (V4 a b d c) <&> \(V4 a' b' d' c') -> V4 a' b' c' d'+{-# INLINE _xywz #-}++_xzyw f = _xyzw $ \(V4 a b c d) -> f (V4 a c b d) <&> \(V4 a' c' b' d') -> V4 a' b' c' d'+{-# INLINE _xzyw #-}++_xzwy f = _xyzw $ \(V4 a b c d) -> f (V4 a c d b) <&> \(V4 a' c' d' b') -> V4 a' b' c' d'+{-# INLINE _xzwy #-}++_xwyz f = _xyzw $ \(V4 a b c d) -> f (V4 a d b c) <&> \(V4 a' d' b' c') -> V4 a' b' c' d'+{-# INLINE _xwyz #-}++_xwzy f = _xyzw $ \(V4 a b c d) -> f (V4 a d c b) <&> \(V4 a' d' c' b') -> V4 a' b' c' d'+{-# INLINE _xwzy #-}++_yxzw f = _xyzw $ \(V4 a b c d) -> f (V4 b a c d) <&> \(V4 b' a' c' d') -> V4 a' b' c' d'+{-# INLINE _yxzw #-}++_yxwz f = _xyzw $ \(V4 a b c d) -> f (V4 b a d c) <&> \(V4 b' a' d' c') -> V4 a' b' c' d'+{-# INLINE _yxwz #-}++_yzxw f = _xyzw $ \(V4 a b c d) -> f (V4 b c a d) <&> \(V4 b' c' a' d') -> V4 a' b' c' d'+{-# INLINE _yzxw #-}++_yzwx f = _xyzw $ \(V4 a b c d) -> f (V4 b c d a) <&> \(V4 b' c' d' a') -> V4 a' b' c' d'+{-# INLINE _yzwx #-}++_ywxz f = _xyzw $ \(V4 a b c d) -> f (V4 b d a c) <&> \(V4 b' d' a' c') -> V4 a' b' c' d'+{-# INLINE _ywxz #-}++_ywzx f = _xyzw $ \(V4 a b c d) -> f (V4 b d c a) <&> \(V4 b' d' c' a') -> V4 a' b' c' d'+{-# INLINE _ywzx #-}++_zxyw f = _xyzw $ \(V4 a b c d) -> f (V4 c a b d) <&> \(V4 c' a' b' d') -> V4 a' b' c' d'+{-# INLINE _zxyw #-}++_zxwy f = _xyzw $ \(V4 a b c d) -> f (V4 c a d b) <&> \(V4 c' a' d' b') -> V4 a' b' c' d'+{-# INLINE _zxwy #-}++_zyxw f = _xyzw $ \(V4 a b c d) -> f (V4 c b a d) <&> \(V4 c' b' a' d') -> V4 a' b' c' d'+{-# INLINE _zyxw #-}++_zywx f = _xyzw $ \(V4 a b c d) -> f (V4 c b d a) <&> \(V4 c' b' d' a') -> V4 a' b' c' d'+{-# INLINE _zywx #-}++_zwxy f = _xyzw $ \(V4 a b c d) -> f (V4 c d a b) <&> \(V4 c' d' a' b') -> V4 a' b' c' d'+{-# INLINE _zwxy #-}++_zwyx f = _xyzw $ \(V4 a b c d) -> f (V4 c d b a) <&> \(V4 c' d' b' a') -> V4 a' b' c' d'+{-# INLINE _zwyx #-}++_wxyz f = _xyzw $ \(V4 a b c d) -> f (V4 d a b c) <&> \(V4 d' a' b' c') -> V4 a' b' c' d'+{-# INLINE _wxyz #-}++_wxzy f = _xyzw $ \(V4 a b c d) -> f (V4 d a c b) <&> \(V4 d' a' c' b') -> V4 a' b' c' d'+{-# INLINE _wxzy #-}++_wyxz f = _xyzw $ \(V4 a b c d) -> f (V4 d b a c) <&> \(V4 d' b' a' c') -> V4 a' b' c' d'+{-# INLINE _wyxz #-}++_wyzx f = _xyzw $ \(V4 a b c d) -> f (V4 d b c a) <&> \(V4 d' b' c' a') -> V4 a' b' c' d'+{-# INLINE _wyzx #-}++_wzxy f = _xyzw $ \(V4 a b c d) -> f (V4 d c a b) <&> \(V4 d' c' a' b') -> V4 a' b' c' d'+{-# INLINE _wzxy #-}++_wzyx f = _xyzw $ \(V4 a b c d) -> f (V4 d c b a) <&> \(V4 d' c' b' a') -> V4 a' b' c' d'+{-# INLINE _wzyx #-}++ew :: R4 t => E t+ew = E _w++instance R1 V4 where _x f (V4 a b c d) = (\a' -> V4 a' b c d) <$> f a {-# INLINE _x #-}++instance R2 V4 where _y f (V4 a b c d) = (\b' -> V4 a b' c d) <$> f b {-# INLINE _y #-} _xy f (V4 a b c d) = (\(V2 a' b') -> V4 a' b' c d) <$> f (V2 a b)@@ -163,10 +449,6 @@ _xyzw = id {-# INLINE _xyzw #-} -instance Core V4 where- core f = V4 (f _x) (f _y) (f _z) (f _w)- {-# INLINE core #-}- instance Storable a => Storable (V4 a) where sizeOf _ = 4 * sizeOf (undefined::a) {-# INLINE sizeOf #-}@@ -183,16 +465,27 @@ where ptr' = castPtr ptr {-# INLINE peek #-} --- | Convert a 3-dimensional affine vector into a 4-dimensional homogeneous vector.+-- | Convert a 3-dimensional affine vector into a 4-dimensional homogeneous vector,+-- i.e. sets the @w@ coordinate to 0. vector :: Num a => V3 a -> V4 a vector (V3 a b c) = V4 a b c 0 {-# INLINE vector #-} --- | Convert a 3-dimensional affine point into a 4-dimensional homogeneous vector.+-- | Convert a 3-dimensional affine point into a 4-dimensional homogeneous vector,+-- i.e. sets the @w@ coordinate to 1. point :: Num a => V3 a -> V4 a point (V3 a b c) = V4 a b c 1 {-# INLINE point #-} +-- | Convert 4-dimensional projective coordinates to a 3-dimensional+-- point. This operation may be denoted, @euclidean [x:y:z:w] = (x\/w,+-- y\/w, z\/w)@ where the projective, homogenous, coordinate+-- @[x:y:z:w]@ is one of many associated with a single point @(x\/w,+-- y\/w, z\/w)@.+normalizePoint :: Fractional a => V4 a -> V3 a+normalizePoint (V4 a b c w) = (1/w) *^ V3 a b c+{-# INLINE normalizePoint #-}+ instance Epsilon a => Epsilon (V4 a) where nearZero = nearZero . quadrance {-# INLINE nearZero #-}@@ -219,3 +512,147 @@ inRange (l1,u1) i1 && inRange (l2,u2) i2 && inRange (l3,u3) i3 && inRange (l4,u4) i4 {-# INLINE inRange #-}++instance Representable V4 where+ type Rep V4 = E V4+ tabulate f = V4 (f ex) (f ey) (f ez) (f ew)+ {-# INLINE tabulate #-}+ index xs (E l) = view l xs+ {-# INLINE index #-}++instance WithIndex.FunctorWithIndex (E V4) V4 where+ imap f (V4 a b c d) = V4 (f ex a) (f ey b) (f ez c) (f ew d)+ {-# INLINE imap #-}++instance WithIndex.FoldableWithIndex (E V4) V4 where+ ifoldMap f (V4 a b c d) = f ex a `mappend` f ey b `mappend` f ez c `mappend` f ew d+ {-# INLINE ifoldMap #-}++instance WithIndex.TraversableWithIndex (E V4) V4 where+ itraverse f (V4 a b c d) = V4 <$> f ex a <*> f ey b <*> f ez c <*> f ew d+ {-# INLINE itraverse #-}++#if !MIN_VERSION_lens(5,0,0)+instance Lens.FunctorWithIndex (E V4) V4 where imap = WithIndex.imap+instance Lens.FoldableWithIndex (E V4) V4 where ifoldMap = WithIndex.ifoldMap+instance Lens.TraversableWithIndex (E V4) V4 where itraverse = WithIndex.itraverse+#endif++type instance Index (V4 a) = E V4+type instance IxValue (V4 a) = a++instance Ixed (V4 a) where+ ix i = el i++instance Each (V4 a) (V4 b) a b where+ each = traverse++data instance U.Vector (V4 a) = V_V4 {-# UNPACK #-} !Int !(U.Vector a)+data instance U.MVector s (V4 a) = MV_V4 {-# UNPACK #-} !Int !(U.MVector s a)+instance U.Unbox a => U.Unbox (V4 a)++instance U.Unbox a => M.MVector U.MVector (V4 a) where+ basicLength (MV_V4 n _) = n+ basicUnsafeSlice m n (MV_V4 _ v) = MV_V4 n (M.basicUnsafeSlice (4*m) (4*n) v)+ basicOverlaps (MV_V4 _ v) (MV_V4 _ u) = M.basicOverlaps v u+ basicUnsafeNew n = liftM (MV_V4 n) (M.basicUnsafeNew (4*n))+ basicUnsafeRead (MV_V4 _ v) i =+ do let o = 4*i+ x <- M.basicUnsafeRead v o+ y <- M.basicUnsafeRead v (o+1)+ z <- M.basicUnsafeRead v (o+2)+ w <- M.basicUnsafeRead v (o+3)+ return (V4 x y z w)+ basicUnsafeWrite (MV_V4 _ v) i (V4 x y z w) =+ do let o = 4*i+ M.basicUnsafeWrite v o x+ M.basicUnsafeWrite v (o+1) y+ M.basicUnsafeWrite v (o+2) z+ M.basicUnsafeWrite v (o+3) w+ basicInitialize (MV_V4 _ v) = M.basicInitialize v++instance U.Unbox a => G.Vector U.Vector (V4 a) where+ basicUnsafeFreeze (MV_V4 n v) = liftM ( V_V4 n) (G.basicUnsafeFreeze v)+ basicUnsafeThaw ( V_V4 n v) = liftM (MV_V4 n) (G.basicUnsafeThaw v)+ basicLength ( V_V4 n _) = n+ basicUnsafeSlice m n (V_V4 _ v) = V_V4 n (G.basicUnsafeSlice (4*m) (4*n) v)+ basicUnsafeIndexM (V_V4 _ v) i =+ do let o = 4*i+ x <- G.basicUnsafeIndexM v o+ y <- G.basicUnsafeIndexM v (o+1)+ z <- G.basicUnsafeIndexM v (o+2)+ w <- G.basicUnsafeIndexM v (o+3)+ return (V4 x y z w)++instance MonadZip V4 where+ mzipWith = liftA2++instance MonadFix V4 where+ mfix f = V4 (let V4 a _ _ _ = f a in a)+ (let V4 _ a _ _ = f a in a)+ (let V4 _ _ a _ = f a in a)+ (let V4 _ _ _ a = f a in a)++instance Bounded a => Bounded (V4 a) where+ minBound = pure minBound+ {-# INLINE minBound #-}+ maxBound = pure maxBound+ {-# INLINE maxBound #-}++instance NFData a => NFData (V4 a) where+ rnf (V4 a b c d) = rnf a `seq` rnf b `seq` rnf c `seq` rnf d++instance Serial1 V4 where+ serializeWith = traverse_+ deserializeWith k = V4 <$> k <*> k <*> k <*> k++instance Serial a => Serial (V4 a) where+ serialize = serializeWith serialize+ deserialize = deserializeWith deserialize++instance Binary a => Binary (V4 a) where+ put = serializeWith Binary.put+ get = deserializeWith Binary.get++instance Serialize a => Serialize (V4 a) where+ put = serializeWith Cereal.put+ get = deserializeWith Cereal.get++instance Eq1 V4 where+ liftEq k (V4 a b c d) (V4 e f g h) = k a e && k b f && k c g && k d h+instance Ord1 V4 where+ liftCompare k (V4 a b c d) (V4 e f g h) = k a e `mappend` k b f `mappend` k c g `mappend` k d h+instance Read1 V4 where+ liftReadsPrec k _ z = readParen (z > 10) $ \r ->+ [ (V4 a b c d, r5)+ | ("V4",r1) <- lex r+ , (a,r2) <- k 11 r1+ , (b,r3) <- k 11 r2+ , (c,r4) <- k 11 r3+ , (d,r5) <- k 11 r4+ ]+instance Show1 V4 where+ liftShowsPrec f _ z (V4 a b c d) = showParen (z > 10) $+ showString "V4 " . f 11 a . showChar ' ' . f 11 b . showChar ' ' . f 11 c . showChar ' ' . f 11 d++instance Field1 (V4 a) (V4 a) a a where+ _1 f (V4 x y z w) = f x <&> \x' -> V4 x' y z w++instance Field2 (V4 a) (V4 a) a a where+ _2 f (V4 x y z w) = f y <&> \y' -> V4 x y' z w++instance Field3 (V4 a) (V4 a) a a where+ _3 f (V4 x y z w) = f z <&> \z' -> V4 x y z' w++instance Field4 (V4 a) (V4 a) a a where+ _4 f (V4 x y z w) = f w <&> \w' -> V4 x y z w'++instance Semigroup a => Semigroup (V4 a) where+ (<>) = liftA2 (<>)++instance Monoid a => Monoid (V4 a) where+ mempty = pure mempty+#if !(MIN_VERSION_base(4,11,0))+ mappend = liftA2 mappend+#endif+
src/Linear/Vector.hs view
@@ -1,9 +1,14 @@ {-# LANGUAGE CPP #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE RankNTypes #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE Trustworthy #-} {-# LANGUAGE DefaultSignatures #-} ----------------------------------------------------------------------------- -- |--- Module : Linear.Vector--- Copyright : (C) 2012 Edward Kmett+-- Copyright : (C) 2012-2015 Edward Kmett -- License : BSD-style (see the file LICENSE) -- Maintainer : Edward Kmett <ekmett@gmail.com> -- Stability : provisional@@ -13,41 +18,103 @@ ----------------------------------------------------------------------------- module Linear.Vector ( Additive(..)+ , E(..) , negated , (^*) , (*^) , (^/)+ , sumV , basis , basisFor+ , scaled+ , outer+ , unit ) where import Control.Applicative+import Control.Lens import Data.Complex-import Data.Foldable (foldMap)-import Data.Functor.Bind-import Data.Functor.Identity+import Data.Foldable as Foldable (forM_, foldl')+import Data.Functor.Compose+import Data.Functor.Product import Data.HashMap.Lazy as HashMap import Data.Hashable import Data.IntMap as IntMap import Data.Map as Map-import Data.Monoid (Sum(..))-import Data.Traversable (Traversable, mapAccumL)+import qualified Data.Vector as Vector+import Data.Vector (Vector)+import qualified Data.Vector.Mutable as Mutable+import GHC.Generics import Linear.Instances () -- $setup--- >>> import Control.Lens -- >>> import Linear.V2 +-- | Basis element+newtype E t = E { el :: forall x. Lens' (t x) x }+ infixl 6 ^+^, ^-^ infixl 7 ^*, *^, ^/ +class GAdditive f where+ gzero :: Num a => f a+ gliftU2 :: (a -> a -> a) -> f a -> f a -> f a+ gliftI2 :: (a -> b -> c) -> f a -> f b -> f c++instance GAdditive U1 where+ gzero = U1+ {-# INLINE gzero #-}+ gliftU2 _ U1 U1 = U1+ {-# INLINE gliftU2 #-}+ gliftI2 _ U1 U1 = U1+ {-# INLINE gliftI2 #-}++instance (GAdditive f, GAdditive g) => GAdditive (f :*: g) where+ gzero = gzero :*: gzero+ {-# INLINE gzero #-}+ gliftU2 f (a :*: b) (c :*: d) = gliftU2 f a c :*: gliftU2 f b d+ {-# INLINE gliftU2 #-}+ gliftI2 f (a :*: b) (c :*: d) = gliftI2 f a c :*: gliftI2 f b d+ {-# INLINE gliftI2 #-}++instance (Additive f, GAdditive g) => GAdditive (f :.: g) where+ gzero = Comp1 $ gzero <$ (zero :: f Int)+ {-# INLINE gzero #-}+ gliftU2 f (Comp1 a) (Comp1 b) = Comp1 $ liftU2 (gliftU2 f) a b+ {-# INLINE gliftU2 #-}+ gliftI2 f (Comp1 a) (Comp1 b) = Comp1 $ liftI2 (gliftI2 f) a b+ {-# INLINE gliftI2 #-}++instance Additive f => GAdditive (Rec1 f) where+ gzero = Rec1 zero+ {-# INLINE gzero #-}+ gliftU2 f (Rec1 g) (Rec1 h) = Rec1 (liftU2 f g h)+ {-# INLINE gliftU2 #-}+ gliftI2 f (Rec1 g) (Rec1 h) = Rec1 (liftI2 f g h)+ {-# INLINE gliftI2 #-}++instance GAdditive f => GAdditive (M1 i c f) where+ gzero = M1 gzero+ {-# INLINE gzero #-}+ gliftU2 f (M1 g) (M1 h) = M1 (gliftU2 f g h)+ {-# INLINE gliftU2 #-}+ gliftI2 f (M1 g) (M1 h) = M1 (gliftI2 f g h)+ {-# INLINE gliftI2 #-}++instance GAdditive Par1 where+ gzero = Par1 0+ gliftU2 f (Par1 a) (Par1 b) = Par1 (f a b)+ {-# INLINE gliftU2 #-}+ gliftI2 f (Par1 a) (Par1 b) = Par1 (f a b)+ {-# INLINE gliftI2 #-}+ -- | A vector is an additive group with additional structure.-class Bind f => Additive f where+class Functor f => Additive f where -- | The zero vector zero :: Num a => f a #ifndef HLINT- default zero :: (Applicative f, Num a) => f a- zero = pure 0+ default zero :: (GAdditive (Rep1 f), Generic1 f, Num a) => f a+ zero = to1 gzero #endif -- | Compute the sum of two vectors@@ -55,48 +122,162 @@ -- >>> V2 1 2 ^+^ V2 3 4 -- V2 4 6 (^+^) :: Num a => f a -> f a -> f a-#ifndef HLINT- default (^+^) :: (Applicative f, Num a) => f a -> f a -> f a- (^+^) = liftA2 (+)+ (^+^) = liftU2 (+) {-# INLINE (^+^) #-}-#endif -- | Compute the difference between two vectors --- -- >>> V2 4 5 - V2 3 1+ -- >>> V2 4 5 ^-^ V2 3 1 -- V2 1 4 (^-^) :: Num a => f a -> f a -> f a-#ifndef HLINT- default (^-^) :: (Applicative f, Num a) => f a -> f a -> f a- (^-^) = liftA2 (-)- {-# INLINE (^-^) #-}-#endif+ x ^-^ y = x ^+^ negated y -- | Linearly interpolate between two vectors.+ --+ -- /Since linear version 1.23, interpolation direction has been reversed; now/+ --+ -- > lerp 0 a b == a+ -- > lerp 1 a b == b lerp :: Num a => a -> f a -> f a -> f a- lerp alpha u v = alpha *^ u ^+^ (1 - alpha) *^ v+ lerp alpha u v = (1 - alpha) *^ u ^+^ alpha *^ v {-# INLINE lerp #-} + -- | Apply a function to merge the 'non-zero' components of two vectors, unioning the rest of the values.+ --+ -- * For a dense vector this is equivalent to 'liftA2'.+ --+ -- * For a sparse vector this is equivalent to 'unionWith'.+ liftU2 :: (a -> a -> a) -> f a -> f a -> f a+#ifndef HLINT+ default liftU2 :: Applicative f => (a -> a -> a) -> f a -> f a -> f a+ liftU2 = liftA2+ {-# INLINE liftU2 #-}+#endif++ -- | Apply a function to the components of two vectors.+ --+ -- * For a dense vector this is equivalent to 'liftA2'.+ --+ -- * For a sparse vector this is equivalent to 'intersectionWith'.+ liftI2 :: (a -> b -> c) -> f a -> f b -> f c+#ifndef HLINT+ default liftI2 :: Applicative f => (a -> b -> c) -> f a -> f b -> f c+ liftI2 = liftA2+ {-# INLINE liftI2 #-}+#endif++instance (Additive f, Additive g) => Additive (Product f g) where+ zero = Pair zero zero+ liftU2 f (Pair a b) (Pair c d) = Pair (liftU2 f a c) (liftU2 f b d)+ liftI2 f (Pair a b) (Pair c d) = Pair (liftI2 f a c) (liftI2 f b d)+ Pair a b ^+^ Pair c d = Pair (a ^+^ c) (b ^+^ d)+ Pair a b ^-^ Pair c d = Pair (a ^-^ c) (b ^-^ d)+ lerp alpha (Pair a b) (Pair c d) = Pair (lerp alpha a c) (lerp alpha b d)++instance (Additive f, Additive g) => Additive (Compose f g) where+ zero = Compose $ zero <$ (zero :: f Int)+ {-# INLINE zero #-}+ Compose a ^+^ Compose b = Compose $ liftU2 (^+^) a b+ {-# INLINE (^+^) #-}+ Compose a ^-^ Compose b = Compose $ liftU2 (^-^) a b+ {-# INLINE (^-^) #-}+ liftU2 f (Compose a) (Compose b) = Compose $ liftU2 (liftU2 f) a b+ {-# INLINE liftU2 #-}+ liftI2 f (Compose a) (Compose b) = Compose $ liftI2 (liftI2 f) a b+ {-# INLINE liftI2 #-}++instance Additive ZipList where+ zero = ZipList []+ {-# INLINE zero #-}+ liftU2 f (ZipList xs) (ZipList ys) = ZipList (liftU2 f xs ys)+ {-# INLINE liftU2 #-}+ liftI2 = liftA2+ {-# INLINE liftI2 #-}++instance Additive Vector where+ zero = mempty+ {-# INLINE zero #-}+ liftU2 f u v = case compare lu lv of+ LT | lu == 0 -> v+ | otherwise -> Vector.modify (\ w -> Foldable.forM_ [0..lu-1] $ \i -> Mutable.unsafeWrite w i $ f (Vector.unsafeIndex u i) (Vector.unsafeIndex v i)) v+ EQ -> Vector.zipWith f u v+ GT | lv == 0 -> u+ | otherwise -> Vector.modify (\ w -> Foldable.forM_ [0..lv-1] $ \i -> Mutable.unsafeWrite w i $ f (Vector.unsafeIndex u i) (Vector.unsafeIndex v i)) u+ where+ lu = Vector.length u+ lv = Vector.length v+ {-# INLINE liftU2 #-}+ liftI2 = Vector.zipWith+ {-# INLINE liftI2 #-}++instance Additive Maybe where+ zero = Nothing+ {-# INLINE zero #-}+ liftU2 f (Just a) (Just b) = Just (f a b)+ liftU2 _ Nothing ys = ys+ liftU2 _ xs Nothing = xs+ {-# INLINE liftU2 #-}+ liftI2 = liftA2+ {-# INLINE liftI2 #-}++instance Additive [] where+ zero = []+ {-# INLINE zero #-}+ liftU2 f = go where+ go (x:xs) (y:ys) = f x y : go xs ys+ go [] ys = ys+ go xs [] = xs+ {-# INLINE liftU2 #-}+ liftI2 = Prelude.zipWith+ {-# INLINE liftI2 #-}+ instance Additive IntMap where zero = IntMap.empty- (^+^) = IntMap.unionWith (+)- xs ^-^ ys = IntMap.unionWith (+) xs (negated ys)+ {-# INLINE zero #-}+ liftU2 = IntMap.unionWith+ {-# INLINE liftU2 #-}+ liftI2 = IntMap.intersectionWith+ {-# INLINE liftI2 #-} instance Ord k => Additive (Map k) where zero = Map.empty- (^+^) = Map.unionWith (+)- xs ^-^ ys = Map.unionWith (+) xs (negated ys)+ {-# INLINE zero #-}+ liftU2 = Map.unionWith+ {-# INLINE liftU2 #-}+ liftI2 = Map.intersectionWith+ {-# INLINE liftI2 #-} instance (Eq k, Hashable k) => Additive (HashMap k) where zero = HashMap.empty- (^+^) = HashMap.unionWith (+)- xs ^-^ ys = HashMap.unionWith (+) xs (negated ys)+ {-# INLINE zero #-}+ liftU2 = HashMap.unionWith+ {-# INLINE liftU2 #-}+ liftI2 = HashMap.intersectionWith+ {-# INLINE liftI2 #-} -instance Additive ((->) b)+instance Additive ((->) b) where+ zero = const 0+ {-# INLINE zero #-}+ liftU2 = liftA2+ {-# INLINE liftU2 #-}+ liftI2 = liftA2+ {-# INLINE liftI2 #-} -instance Additive Complex+instance Additive Complex where+ zero = 0 :+ 0+ {-# INLINE zero #-}+ liftU2 f (a :+ b) (c :+ d) = f a c :+ f b d+ {-# INLINE liftU2 #-}+ liftI2 f (a :+ b) (c :+ d) = f a c :+ f b d+ {-# INLINE liftI2 #-} -instance Additive Identity+instance Additive Identity where+ zero = Identity 0+ {-# INLINE zero #-}+ liftU2 = liftA2+ {-# INLINE liftU2 #-}+ liftI2 = liftA2+ {-# INLINE liftI2 #-} -- | Compute the negation of a vector --@@ -106,6 +287,14 @@ negated = fmap negate {-# INLINE negated #-} +-- | Sum over multiple vectors+--+-- >>> sumV [V2 1 1, V2 3 4]+-- V2 4 5+sumV :: (Foldable f, Additive v, Num a) => f (v a) -> v a+sumV = Foldable.foldl' (^+^) zero+{-# INLINE sumV #-}+ -- | Compute the left scalar product -- -- >>> 2 *^ V2 3 4@@ -127,23 +316,39 @@ f ^/ a = fmap (/a) f {-# INLINE (^/) #-} --- @setElement i x v@ sets the @i@'th element of @v@ to @x@.-setElement :: Traversable t => Int -> a -> t a -> t a-setElement i x = snd . mapAccumL aux 0- where aux j y = let j' = j + 1- y' = if i == j then x else y- in j' `seq` (j', y')- -- | Produce a default basis for a vector space. If the dimensionality -- of the vector space is not statically known, see 'basisFor'.-basis :: (Applicative t, Traversable t, Num a) => [t a]-basis = [ setElement k 1 z | k <- [0..n - 1] ]- where z = pure 0- n = getSum $ foldMap (const (Sum 1)) z+basis :: (Additive t, Traversable t, Num a) => [t a]+basis = basisFor (zero :: Additive v => v Int) -- | Produce a default basis for a vector space from which the -- argument is drawn.-basisFor :: (Traversable t, Enum a, Num a) => t a -> [t a]-basisFor v = [ setElement k 1 z | k <- [0..n-1] ]- where z = 0 <$ v- n = getSum $ foldMap (const (Sum 1)) v+basisFor :: (Traversable t, Num a) => t b -> [t a]+basisFor = \t ->+ ifoldMapOf traversed ?? t $ \i _ ->+ return $+ iover traversed ?? t $ \j _ ->+ if i == j then 1 else 0+{-# INLINABLE basisFor #-}++-- | Produce a diagonal (scale) matrix from a vector.+--+-- >>> scaled (V2 2 3)+-- V2 (V2 2 0) (V2 0 3)+scaled :: (Traversable t, Num a) => t a -> t (t a)+scaled = \t -> iter t (\i x -> iter t (\j _ -> if i == j then x else 0))+ where+ iter :: Traversable t => t a -> (Int -> a -> b) -> t b+ iter x f = iover traversed f x+{-# INLINE scaled #-}++-- | Create a unit vector.+--+-- >>> unit _x :: V2 Int+-- V2 1 0+unit :: (Additive t, Num a) => ASetter' (t a) a -> t a+unit l = set' l 1 zero++-- | Outer (tensor) product of two vectors+outer :: (Functor f, Functor g, Num a) => f a -> g a -> f (g a)+outer a b = fmap (\x->fmap (*x) b) a
+ tests/Prop/Quaternion.hs view
@@ -0,0 +1,28 @@+{-# OPTIONS_GHC -Wno-orphans #-}+module Prop.Quaternion (tests) where++import Linear.Quaternion (Quaternion(..))+import Linear.Epsilon (nearZero)+import Linear.Vector (lerp)+import Test.QuickCheck (Arbitrary(..))+import Test.Tasty (TestTree, testGroup)+import Test.Tasty.QuickCheck (testProperty)++import Prop.V3 ()++instance Arbitrary a => Arbitrary (Quaternion a) where+ arbitrary = Quaternion <$> arbitrary <*> arbitrary++prop_lerp0 :: Quaternion Double -> Quaternion Double -> Bool+prop_lerp0 a b = nearZero (lerp 0 a b - a)++prop_lerp1 :: Quaternion Double -> Quaternion Double -> Bool+prop_lerp1 a b = nearZero (lerp 1 a b - b)++tests :: [TestTree]+tests =+ [ testGroup "lerp"+ [ testProperty "lerp 0 a b == a" prop_lerp0+ , testProperty "lerp 1 a b == b" prop_lerp1+ ]+ ]
+ tests/Prop/V3.hs view
@@ -0,0 +1,8 @@+{-# OPTIONS_GHC -Wno-orphans #-}+module Prop.V3 () where++import Linear.V3 (V3(..))+import Test.QuickCheck (Arbitrary(..))++instance Arbitrary a => Arbitrary (V3 a) where+ arbitrary = V3 <$> arbitrary <*> arbitrary <*> arbitrary
+ tests/Test.hs view
@@ -0,0 +1,24 @@+{-# LANGUAGE CPP #-}+module Main (main) where++import Test.Tasty (defaultMain, testGroup, TestTree)++import qualified Prop.Quaternion+import qualified Unit.Binary+import qualified Unit.Plucker+import qualified Unit.V++tests :: [TestTree]+tests =+ [ testGroup "Property tests"+ [ testGroup "Quaternion" Prop.Quaternion.tests+ ]+ , testGroup "Unit tests"+ [ testGroup "Binary" Unit.Binary.tests+ , testGroup "Plucker" Unit.Plucker.tests+ , testGroup "V" Unit.V.tests+ ]+ ]++main :: IO ()+main = defaultMain $ testGroup "linear" tests
+ tests/Unit/Binary.hs view
@@ -0,0 +1,20 @@+module Unit.Binary (tests) where++import Data.Binary.Put+import Data.Binary.Get+import Linear+import qualified Data.ByteString.Lazy as BS+import Test.Tasty (TestTree)+import Test.Tasty.HUnit ((@?=), testCase)++originalVecs :: (V3 Float, V2 Char)+originalVecs = (V3 1 2 3, V2 'a' 'b')++bytes :: BS.ByteString+bytes = runPut $ do putLinear $ fst originalVecs+ putLinear $ snd originalVecs++tests :: [TestTree]+tests = [ testCase "Serialized length" $ BS.length bytes @?= 3*13+2+ , testCase "Deserialization" $ deserialized @?= originalVecs ]+ where deserialized = runGet ((,) <$> getLinear <*> getLinear) bytes
+ tests/Unit/Plucker.hs view
@@ -0,0 +1,36 @@+module Unit.Plucker (tests) where+import Linear+import Linear.Plucker+import Linear.Plucker.Coincides+import Test.Tasty (TestTree)+import Test.Tasty.HUnit ((@?=), testCase)++ln2,ln3,ln4,ln5,ln6,ln7,ln8,ln9 :: Plucker Float+ln2 = plucker3D (V3 1 3 0) (V3 1 3 (-2)) -- starting line+ln3 = plucker3D (V3 2 3 0) (V3 2 3 (-2)) -- parallel+ln4 = plucker3D (V3 2 4 0) (V3 1 4 (-2)) -- ccw+ln5 = plucker3D (V3 (-2) 4 0) (V3 2 4 (-2)) -- cw+ln6 = plucker3D (V3 2 3 0) (V3 1 3 (-2)) -- intersect+ln7 = plucker3D (V3 1 3 0) (V3 1 3 2) -- reversed+ln8 = plucker3D (V3 0 4 4) (V3 0 (-4) (-4)) -- through origin+ln9 = Plucker 1 2 3 4 5 6 -- not a 3D line++tests :: [TestTree]+tests = [ testCase "parallel" $ parallel ln2 ln3 @?= True+ , testCase "CCW" $ passes ln2 ln4 @?= Counterclockwise+ , testCase "CW" $ passes ln2 ln5 @?= Clockwise+ , testCase "intersect1" $ intersects ln2 ln6 @?= True+ , testCase "intersect2" $ intersects ln2 ln3 @?= False+ , testCase "line equality 1" $ Line ln2 == Line ln2 @?= True+ , testCase "line equality 2" $ Line ln2 == Line ln7 @?= True+ , testCase "line equality 3" $ Line ln2 == Ray ln7 @?= True+ , testCase "line equality 4" $ Ray ln2 == Line ln7 @?= True+ , testCase "ray equality 1" $ Ray ln2 == Ray ln7 @?= False+ , testCase "ray equality 2" $ Ray ln2 == Ray (3 *^ ln2) @?= True+ , testCase "ray equality 3" $ Ray ln2 == Ray (negate ln7) @?= True+ , testCase "quadrance" $ nearZero (quadranceToOrigin ln2 - 10) @?= True+ , testCase "closest 1" $+ nearZero (qd (V3 1 3 0) $ closestToOrigin ln2) @?= True+ , testCase "closest 2" $ nearZero (qd 0 $ closestToOrigin ln8) @?= True+ , testCase "isLine 1" $ isLine ln2 @?= True+ , testCase "isLine 2" $ isLine ln9 @?= False ]
+ tests/Unit/V.hs view
@@ -0,0 +1,14 @@+{-# LANGUAGE DataKinds #-}+module Unit.V (tests) where++import Control.DeepSeq (rnf)+import qualified Data.Vector.Unboxed as U (fromList)+import Linear.V (V)+import Test.Tasty (TestTree)+import Test.Tasty.HUnit ((@?=), testCase)++v10 :: V 10 Int+v10 = return 5++tests :: [TestTree]+tests = [ testCase "GH124" $ rnf (U.fromList [v10]) @?= () ]
tests/doctests.hs view
@@ -1,30 +1,19 @@+-----------------------------------------------------------------------------+-- |+-- Module : Main (doctests)+-- Copyright : (C) 2012-14 Edward Kmett+-- License : BSD-style (see the file LICENSE)+-- Maintainer : Edward Kmett <ekmett@gmail.com>+-- Stability : provisional+-- Portability : portable+--+-- This module exists to add dependencies+----------------------------------------------------------------------------- module Main where -import Build_doctests (deps)-import Control.Applicative-import Control.Monad-import Data.List-import System.Directory-import System.FilePath-import Test.DocTest- main :: IO ()-main = getSources >>= \sources -> doctest $- "-isrc"- : "-idist/build/autogen"- : "-optP-include"- : "-optPdist/build/autogen/cabal_macros.h"- : "-hide-all-packages"- : map ("-package="++) deps ++ sources--getSources :: IO [FilePath]-getSources = filter (isSuffixOf ".hs") <$> go "src"- where- go dir = do- (dirs, files) <- getFilesAndDirectories dir- (files ++) . concat <$> mapM go dirs--getFilesAndDirectories :: FilePath -> IO ([FilePath], [FilePath])-getFilesAndDirectories dir = do- c <- map (dir </>) . filter (`notElem` ["..", "."]) <$> getDirectoryContents dir- (,) <$> filterM doesDirectoryExist c <*> filterM doesFileExist c+main = do+ putStrLn "This test-suite exists only to add dependencies"+ putStrLn "To run doctests: "+ putStrLn " cabal build all --enable-tests"+ putStrLn " cabal-docspec"
− travis/cabal-apt-install
@@ -1,16 +0,0 @@-#!/bin/sh-set -eu--sudo apt-get -q update-sudo apt-get -q -y install dctrl-tools--# Try installing some of the build-deps with apt-get for speed.-eval "$(- printf '%s' "grep-aptavail -n -sPackage '(' -FFALSE -X FALSE ')'"- 2>/dev/null cabal install "$@" --dry-run -v | \- sed -nre "s/^([^ ]+)-[0-9.]+ \(.*$/ -o '(' -FPackage -X libghc-\1-dev ')'/p" | \- xargs -d'\n'-)" | sort -u | xargs -d'\n' sudo apt-get -q -y install -- libghc-quickcheck2-dev--# Install whatever is still needed with cabal.-cabal install "$@"
− travis/config
@@ -1,16 +0,0 @@--- This provides a custom ~/.cabal/config file for use when hackage is down that should work on unix------ This is particularly useful for travis-ci to get it to stop complaining--- about a broken build when everything is still correct on our end.------ This uses Luite Stegeman's mirror of hackage provided by his 'hdiff' site instead------ To enable this, uncomment the before_script in .travis.yml--remote-repo: hdiff.luite.com:http://hdiff.luite.com/packages/archive-remote-repo-cache: ~/.cabal/packages-world-file: ~/.cabal/world-build-summary: ~/.cabal/logs/build.log-remote-build-reporting: anonymous-install-dirs user-install-dirs global