linear 1.20.1 → 1.20.2
raw patch · 11 files changed
+165/−31 lines, 11 filesdep +HerbiePlugindep ~doctestdep ~vectorsetup-changed
Dependencies added: HerbiePlugin
Dependency ranges changed: doctest, vector
Files
- .ghci +0/−1
- CHANGELOG.markdown +6/−0
- Setup.lhs +4/−0
- linear.cabal +12/−3
- src/Linear/Affine.hs +6/−6
- src/Linear/Instances.hs +8/−0
- src/Linear/Matrix.hs +12/−6
- src/Linear/Projection.hs +16/−4
- src/Linear/Quaternion.hs +50/−8
- src/Linear/V.hs +48/−0
- tests/doctests.hs +3/−3
− .ghci
@@ -1,1 +0,0 @@-:set -isrc -idist/build/autogen -optP-include -optPdist/build/autogen/cabal_macros.h
CHANGELOG.markdown view
@@ -1,3 +1,9 @@+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.
Setup.lhs view
@@ -30,6 +30,10 @@ withTestLBI pkg lbi $ \suite suitecfg -> do rewriteFile (dir </> "Build_" ++ testName suite ++ ".hs") $ unlines [ "module Build_" ++ testName suite ++ " where"+ , ""+ , "autogen_dir :: String"+ , "autogen_dir = " ++ show dir+ , "" , "deps :: [String]" , "deps = " ++ (show $ formatdeps (testDeps libcfg suitecfg)) ]
linear.cabal view
@@ -1,6 +1,6 @@ name: linear category: Math, Algebra-version: 1.20.1+version: 1.20.2 license: BSD3 cabal-version: >= 1.8 license-file: LICENSE@@ -15,7 +15,6 @@ build-type: Custom tested-with: GHC == 7.4.2, GHC == 7.6.1, GHC == 7.8.4, GHC == 7.10.1 extra-source-files:- .ghci .gitignore .travis.yml .vim.custom@@ -32,6 +31,11 @@ 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@@ -62,6 +66,11 @@ if flag(template-haskell) && impl(ghc) build-depends: template-haskell >= 2.7 && < 3.0 + if flag(herbie)+ build-depends: HerbiePlugin >= 0.1 && < 0.2+ ghc-options: -fplugin=Herbie+ cpp-options: -DHERBIE+ exposed-modules: Linear Linear.Affine@@ -98,7 +107,7 @@ build-depends: base, directory >= 1.0 && < 1.3,- doctest >= 0.8 && < 0.10,+ doctest >= 0.8 && < 0.11, filepath >= 1.3 && < 1.5, lens, simple-reflect >= 0.3.1
src/Linear/Affine.hs view
@@ -48,10 +48,10 @@ import Data.Serialize as Cereal import Data.Vector (Vector) import Foreign.Storable-#if defined(__GLASGOW_HASKELL__) && __GLASGOW_HASKELL__ >= 702+#if __GLASGOW_HASKELL__ >= 702 import GHC.Generics (Generic) #endif-#if defined(__GLASGOW_HASKELL__) && __GLASGOW_HASKELL__ >= 706+#if __GLASGOW_HASKELL__ >= 706 import GHC.Generics (Generic1) #endif import Linear.Epsilon@@ -135,13 +135,13 @@ , Traversable, Apply, Additive, Metric , Fractional , Num, Ix, Storable, Epsilon , Hashable-#if defined(__GLASGOW_HASKELL__) && __GLASGOW_HASKELL__ >= 702+#if __GLASGOW_HASKELL__ >= 702 , Generic #endif-#if defined(__GLASGOW_HASKELL__) && __GLASGOW_HASKELL__ >= 706+#if __GLASGOW_HASKELL__ >= 706 , Generic1 #endif-#if defined(__GLASGOW_HASKELL__) && __GLASGOW_HASKELL__ >= 708+#if __GLASGOW_HASKELL__ >= 708 , Typeable, Data #endif )@@ -165,7 +165,7 @@ put (P p) = Cereal.put p get = P `liftM` Cereal.get -#if defined(__GLASGOW_HASKELL__) && __GLASGOW_HASKELL__ < 708+#if __GLASGOW_HASKELL__ < 708 instance forall f. Typeable1 f => Typeable1 (Point f) where typeOf1 _ = mkTyConApp (mkTyCon3 "linear" "Linear.Affine" "Point") [] `mkAppTy` typeOf1 (undefined :: f a)
src/Linear/Instances.hs view
@@ -43,16 +43,20 @@ Just b -> [(k,b)] Nothing -> [] +#if __GLASGOW_HASKELL__ < 711 instance Functor Complex where fmap f (a :+ b) = f a :+ f b {-# INLINE fmap #-}+#endif instance Apply Complex where (a :+ b) <.> (c :+ d) = a c :+ b d +#if __GLASGOW_HASKELL__ < 711 instance Applicative Complex where pure a = a :+ a (a :+ b) <*> (c :+ d) = a c :+ b d+#endif instance Bind Complex where (a :+ b) >>- f = a' :+ b' where@@ -60,6 +64,7 @@ _ :+ b' = f b {-# INLINE (>>-) #-} +#if __GLASGOW_HASKELL__ < 711 instance Monad Complex where return a = a :+ a {-# INLINE return #-}@@ -68,6 +73,7 @@ a' :+ _ = f a _ :+ b' = f b {-# INLINE (>>=) #-}+#endif instance MonadZip Complex where mzipWith = liftA2@@ -75,6 +81,7 @@ instance MonadFix Complex where mfix f = (let a :+ _ = f a in a) :+ (let _ :+ a = f a in a) +#if __GLASGOW_HASKELL__ < 711 instance Foldable Complex where foldMap f (a :+ b) = f a `mappend` f b {-# INLINE foldMap #-}@@ -82,6 +89,7 @@ instance Traversable Complex where traverse f (a :+ b) = (:+) <$> f a <*> f b {-# INLINE traverse #-}+#endif instance Foldable1 Complex where foldMap1 f (a :+ b) = f a <> f b
src/Linear/Matrix.hs view
@@ -185,12 +185,18 @@ -- | 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
src/Linear/Projection.hs view
@@ -69,9 +69,19 @@ where tanHalfFovy = tan $ fovy / 2 x = 1 / (aspect * tanHalfFovy) y = 1 / tanHalfFovy- z = -(far + near) / (far - near)- w = -(2 * far * near) / (far - near)+ fpn = far + near+ fmn = far - near+ oon = 1/near+ oof = 1/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@@ -88,8 +98,10 @@ where tanHalfFovy = tan $ fovy / 2 a = aspect * tanHalfFovy b = tanHalfFovy- c = -(far - near) / (2 * far * near)- d = (far + near) / (2 * far * near)+ c = oon - oof+ d = oon + oof+ oon = 1/near+ oof = 1/far -- | Build a perspective matrix per the classic @glFrustum@ arguments.
src/Linear/Quaternion.hs view
@@ -341,6 +341,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) #-}@@ -356,11 +379,14 @@ | 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+ | ai <- sqrt qiq = reimagine (log m) (atan2 m e / 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@@ -368,38 +394,45 @@ 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+ -- reimagine (cos e * cosh ai) (- sin e * sinh ai / ai) q -- 13.5 bits worse 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 #-} @@ -417,6 +450,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
src/Linear/V.hs view
@@ -52,6 +52,7 @@ import Control.Applicative #endif import Control.DeepSeq (NFData)+import Control.Monad import Control.Monad.Fix import Control.Monad.Zip import Control.Lens as Lens@@ -72,6 +73,9 @@ import Data.Traversable (sequenceA) #endif import Data.Vector as V+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 #ifdef USE_TYPE_LITS@@ -462,3 +466,47 @@ instance Dim n => Ord1 (V n) where compare1 = compare instance Dim n => Show1 (V n) where showsPrec1 = showsPrec instance Dim n => Read1 (V n) where readsPrec1 = readsPrec+++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 _ v) i (V vn) =+ V.imapM_ (\j -> M.basicUnsafeWrite v (d*i+j)) vn+ where d = reflectDim (Proxy :: Proxy n)+#if MIN_VERSION_vector(0,11,0)+ basicInitialize (MV_VN _ v) = M.basicInitialize v+ {-# INLINE basicInitialize #-}+#endif++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)
tests/doctests.hs view
@@ -1,6 +1,6 @@ module Main where -import Build_doctests (deps)+import Build_doctests (autogen_dir, deps) import Control.Applicative import Control.Monad import Data.List@@ -11,9 +11,9 @@ main :: IO () main = getSources >>= \sources -> doctest $ "-isrc"- : "-idist/build/autogen"+ : ("-i" ++ autogen_dir) : "-optP-include"- : "-optPdist/build/autogen/cabal_macros.h"+ : ("-optP" ++ autogen_dir ++ "/cabal_macros.h") : "-hide-all-packages" : map ("-package="++) deps ++ sources