diff --git a/.ghci b/.ghci
deleted file mode 100644
--- a/.ghci
+++ /dev/null
@@ -1,1 +0,0 @@
-:set -isrc -idist/build/autogen -optP-include -optPdist/build/autogen/cabal_macros.h
diff --git a/CHANGELOG.markdown b/CHANGELOG.markdown
--- a/CHANGELOG.markdown
+++ b/CHANGELOG.markdown
@@ -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.
diff --git a/Setup.lhs b/Setup.lhs
--- a/Setup.lhs
+++ b/Setup.lhs
@@ -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))
         ]
diff --git a/linear.cabal b/linear.cabal
--- a/linear.cabal
+++ b/linear.cabal
@@ -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
diff --git a/src/Linear/Affine.hs b/src/Linear/Affine.hs
--- a/src/Linear/Affine.hs
+++ b/src/Linear/Affine.hs
@@ -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)
diff --git a/src/Linear/Instances.hs b/src/Linear/Instances.hs
--- a/src/Linear/Instances.hs
+++ b/src/Linear/Instances.hs
@@ -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
diff --git a/src/Linear/Matrix.hs b/src/Linear/Matrix.hs
--- a/src/Linear/Matrix.hs
+++ b/src/Linear/Matrix.hs
@@ -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
diff --git a/src/Linear/Projection.hs b/src/Linear/Projection.hs
--- a/src/Linear/Projection.hs
+++ b/src/Linear/Projection.hs
@@ -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.
diff --git a/src/Linear/Quaternion.hs b/src/Linear/Quaternion.hs
--- a/src/Linear/Quaternion.hs
+++ b/src/Linear/Quaternion.hs
@@ -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
diff --git a/src/Linear/V.hs b/src/Linear/V.hs
--- a/src/Linear/V.hs
+++ b/src/Linear/V.hs
@@ -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)
diff --git a/tests/doctests.hs b/tests/doctests.hs
--- a/tests/doctests.hs
+++ b/tests/doctests.hs
@@ -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