diff --git a/ChangeLog.md b/ChangeLog.md
new file mode 100644
--- /dev/null
+++ b/ChangeLog.md
@@ -0,0 +1,10 @@
+# 0.0.2 (2014-03-25)
+
+* Add doctests.
+* Fix qiModify potentially constructing `qi 1 0 5 1` instead of the equivalent
+  but simpler `qi 1 0 0 1`.
+* Add lenses.
+
+# 0.0.1 (2014-03-24)
+
+* Initial release.
diff --git a/quadratic-irrational.cabal b/quadratic-irrational.cabal
--- a/quadratic-irrational.cabal
+++ b/quadratic-irrational.cabal
@@ -1,6 +1,6 @@
 name: quadratic-irrational
 category: Math, Algorithms, Data
-version: 0.0.1
+version: 0.0.2
 license: MIT
 license-file: LICENSE
 author: Johan Kiviniemi <devel@johan.kiviniemi.name>
@@ -21,6 +21,7 @@
 cabal-version: >= 1.10
 extra-source-files:
   .gitignore
+  ChangeLog.md
   README.md
 
 source-repository head
@@ -30,6 +31,7 @@
 library
   exposed-modules: Numeric.QuadraticIrrational
                  , Numeric.QuadraticIrrational.CyclicList
+                 , Numeric.QuadraticIrrational.Internal.Lens
   hs-source-dirs: src
   build-depends: base >= 4.6 && < 4.8
                , arithmoi == 0.4.*
@@ -39,9 +41,9 @@
   default-language: Haskell2010
   ghc-options: -Wall -O2 -funbox-strict-fields
 
-test-suite test-quadratic-irrational
+test-suite tasty-tests
   type: exitcode-stdio-1.0
-  main-is: Main.hs
+  main-is: tasty.hs
   other-modules: QuadraticIrrational
                , CyclicList
   hs-source-dirs: tests
@@ -53,3 +55,15 @@
                , tasty-quickcheck == 0.8.*
   default-language: Haskell2010
   ghc-options: -Wall -O2 -funbox-strict-fields
+
+test-suite doctests
+  type: exitcode-stdio-1.0
+  main-is: doctests.hs
+  hs-source-dirs: tests
+  build-depends: base
+               , directory
+               , doctest >= 0.9
+               , filepath
+               , mtl
+  default-language: Haskell2010
+  ghc-options: -threaded -Wall
diff --git a/src/Numeric/QuadraticIrrational.hs b/src/Numeric/QuadraticIrrational.hs
--- a/src/Numeric/QuadraticIrrational.hs
+++ b/src/Numeric/QuadraticIrrational.hs
@@ -15,13 +15,19 @@
 -- <http://en.wikipedia.org/wiki/Periodic_continued_fraction periodic continued fractions>.
 
 module Numeric.QuadraticIrrational
-  ( QI, qi, qi', qiModify, runQI, runQI', unQI, unQI'
-  , qiZero, qiOne, qiIsZero
+  ( -- * Constructors and deconstructors
+    QI, qi, qi', runQI, runQI', unQI, unQI'
+  , -- * Lenses
+    _qi, _qi', _qiABD, _qiA, _qiB, _qiC, _qiD
+  , -- * Numerical operations
+    qiZero, qiOne, qiIsZero
   , qiToFloat
   , qiAddI, qiSubI, qiMulI, qiDivI
   , qiAddR, qiSubR, qiMulR, qiDivR
   , qiNegate, qiRecip, qiAdd, qiSub, qiMul, qiDiv, qiPow
-  , qiFloor, continuedFractionToQI, qiToContinuedFraction
+  , qiFloor
+  , -- * Continued fractions
+    continuedFractionToQI, qiToContinuedFraction
   , module Numeric.QuadraticIrrational.CyclicList
   ) where
 
@@ -36,6 +42,7 @@
 import Text.Read
 
 import Numeric.QuadraticIrrational.CyclicList
+import Numeric.QuadraticIrrational.Internal.Lens
 
 -- | @(a + b √c) / d@
 data QI = QI !Integer
@@ -65,6 +72,32 @@
 
 -- | Given @a@, @b@, @c@ and @d@ such that @n = (a + b √c)/d@, constuct a 'QI'
 -- corresponding to @n@.
+--
+-- >>> qi 3 4 5 6
+-- qi 3 4 5 6
+--
+-- The fractions are reduced:
+--
+-- >>> qi 30 40 5 60
+-- qi 3 4 5 6
+--
+-- If @b = 0@ then @c@ is zeroed and vice versa:
+--
+-- >>> qi 3 0 42 1
+-- qi 3 0 0 1
+--
+-- >>> qi 3 42 0 1
+-- qi 3 0 0 1
+--
+-- The @b √c@ term is simplified:
+--
+-- >>> qi 0 1 (5*5*6) 1
+-- qi 0 5 6 1
+--
+-- If @c = 1@ (after simplification) then @b@ is moved to @a@:
+--
+-- >>> qi 1 5 (2*2) 1
+-- qi 11 0 0 1
 qi :: Integer  -- ^ a
    -> Integer  -- ^ b
    -> Integer  -- ^ c
@@ -77,6 +110,16 @@
   | otherwise = simplifyReduceCons a b c d
 {-# INLINE qi #-}
 
+-- Construct a 'QI' without simplifying @b √c@. Make sure it has already been
+-- simplified.
+qiNoSimpl :: Integer -> Integer -> Integer -> Integer -> QI
+qiNoSimpl a b (nonNegative "qiNoSimpl" -> c) (nonZero "qiNoSimpl" -> d)
+  | b == 0    = reduceCons a 0 0 d
+  | c == 0    = reduceCons a 0 0 d
+  | c == 1    = reduceCons (a + b) 0 0 d
+  | otherwise = reduceCons a b c d
+{-# INLINE qiNoSimpl #-}
+
 -- Simplify @b √c@ before constructing a 'QI'.
 simplifyReduceCons :: Integer -> Integer -> Integer -> Integer -> QI
 simplifyReduceCons a b (nonZero "simplifyReduceCons" -> c) d
@@ -108,6 +151,9 @@
 
 -- | Given @a@, @b@ and @c@ such that @n = a + b √c@, constuct a 'QI'
 -- corresponding to @n@.
+--
+-- >>> qi' 0.5 0.7 2
+-- qi 5 7 2 10
 qi' :: Rational  -- ^ a
     -> Rational  -- ^ b
     -> Integer   -- ^ c
@@ -120,116 +166,253 @@
     (bN, bD) = (numerator b, denominator b)
 {-# INLINE qi' #-}
 
--- | Given a 'QI' corresponding to @n = (a + b √c)/d@, modify @(a, b, d)@.
--- Avoids having to simplify @b √c@.
-qiModify :: QI
-         -> (Integer -> Integer -> Integer -> (Integer, Integer, Integer))
-         -> QI
-qiModify (QI a b c d) f = reduceCons a' b' c d'
-  where (a', b', d') = f a b d
-{-# INLINE qiModify #-}
-
 -- | Given @n@ and @f@ such that @n = (a + b √c)/d@, run @f a b c d@.
+--
+-- >>> runQI (qi 3 4 5 6) (\a b c d -> (a,b,c,d))
+-- (3,4,5,6)
 runQI :: QI -> (Integer -> Integer -> Integer -> Integer -> a) -> a
 runQI (QI a b c d) f = f a b c d
 {-# INLINE runQI #-}
 
 -- | Given @n@ and @f@ such that @n = a + b √c@, run @f a b c@.
+--
+-- >>> runQI' (qi' 0.5 0.7 2) (\a b c -> (a, b, c))
+-- (1 % 2,7 % 10,2)
 runQI' :: QI -> (Rational -> Rational -> Integer -> a) -> a
 runQI' (QI a b c d) f = f (a % d) (b % d) c
 {-# INLINE runQI' #-}
 
 -- | Given @n@ such that @n = (a + b √c)/d@, return @(a, b, c, d)@.
+--
+-- >>> unQI (qi 3 4 5 6)
+-- (3,4,5,6)
 unQI :: QI -> (Integer, Integer, Integer, Integer)
 unQI n = runQI n (,,,)
 {-# INLINE unQI #-}
 
 -- | Given @n@ such that @n = a + b √c@, return @(a, b, c)@.
+--
+-- >>> unQI' (qi' 0.5 0.7 2)
+-- (1 % 2,7 % 10,2)
 unQI' :: QI -> (Rational, Rational, Integer)
 unQI' n = runQI' n (,,)
 {-# INLINE unQI' #-}
 
--- | The constant zero. @qi 0 0 0 1@
+-- | Given a 'QI' corresponding to @n = (a + b √c)/d@, access @(a, b, c, d)@.
+--
+-- >>> view _qi (qi 3 4 5 6)
+-- (3,4,5,6)
+--
+-- >>> over _qi (\(a,b,c,d) -> (a+10, b+10, c+10, d+10)) (qi 3 4 5 6)
+-- qi 13 14 15 16
+_qi :: Lens' QI (Integer, Integer, Integer, Integer)
+_qi f n = (\ ~(a',b',c',d') -> qi a' b' c' d') <$> f (unQI n)
+{-# INLINE _qi #-}
+
+-- | Given a 'QI' corresponding to @n = a + b √c@, access @(a, b, c)@.
+--
+-- >>> view _qi' (qi' 0.5 0.7 2)
+-- (1 % 2,7 % 10,2)
+--
+-- >>> over _qi' (\(a,b,c) -> (a/5, b/6, c*3)) (qi 3 4 5 6)
+-- qi 9 10 15 90
+_qi' :: Lens' QI (Rational, Rational, Integer)
+_qi' f n = (\ ~(a',b',c') -> qi' a' b' c') <$> f (unQI' n)
+{-# INLINE _qi' #-}
+
+-- | Given a 'QI' corresponding to @n = (a + b √c)/d@, access @(a, b, d)@.
+-- Avoids having to simplify @b √c@ upon reconstruction.
+--
+-- >>> view _qiABD (qi 3 4 5 6)
+-- (3,4,6)
+--
+-- >>> over _qiABD (\(a,b,d) -> (a+10, b+10, d+10)) (qi 3 4 5 6)
+-- qi 13 14 5 16
+_qiABD :: Lens' QI (Integer, Integer, Integer)
+_qiABD f (unQI -> ~(a,b,c,d)) =
+  (\ ~(a',b',d') -> qiNoSimpl a' b' c d') <$> f (a,b,d)
+{-# INLINE _qiABD #-}
+
+-- | Given a 'QI' corresponding to @n = (a + b √c)/d@, access @a@. It is more
+-- efficient to use '_qi' or '_qiABD' when modifying multiple terms at once.
+--
+-- >>> view _qiA (qi 3 4 5 6)
+-- 3
+--
+-- >>> over _qiA (+ 10) (qi 3 4 5 6)
+-- qi 13 4 5 6
+_qiA :: Lens' QI Integer
+_qiA = _qiABD . go
+  where go f ~(a,b,d) = (\a' -> (a',b,d)) <$> f a
+
+-- | Given a 'QI' corresponding to @n = (a + b √c)/d@, access @b@. It is more
+-- efficient to use '_qi' or '_qiABD' when modifying multiple terms at once.
+--
+-- >>> view _qiB (qi 3 4 5 6)
+-- 4
+--
+-- >>> over _qiB (+ 10) (qi 3 4 5 6)
+-- qi 3 14 5 6
+_qiB :: Lens' QI Integer
+_qiB = _qiABD . go
+  where go f ~(a,b,d) = (\b' -> (a,b',d)) <$> f b
+
+-- | Given a 'QI' corresponding to @n = (a + b √c)/d@, access @c@. It is more
+-- efficient to use '_qi' or '_qiABD' when modifying multiple terms at once.
+--
+-- >>> view _qiC (qi 3 4 5 6)
+-- 5
+--
+-- >>> over _qiC (+ 10) (qi 3 4 5 6)
+-- qi 3 4 15 6
+_qiC :: Lens' QI Integer
+_qiC = _qi . go
+  where go f ~(a,b,c,d) = (\c' -> (a,b,c',d)) <$> f c
+
+-- | Given a 'QI' corresponding to @n = (a + b √c)/d@, access @d@. It is more
+-- efficient to use '_qi' or '_qiABD' when modifying multiple terms at once.
+--
+-- >>> view _qiD (qi 3 4 5 6)
+-- 6
+--
+-- >>> over _qiD (+ 10) (qi 3 4 5 6)
+-- qi 3 4 5 16
+_qiD :: Lens' QI Integer
+_qiD = _qiABD . go
+  where go f ~(a,b,d) = (\d' -> (a,b,d')) <$> f d
+
+-- | The constant zero.
+--
+-- >>> qiZero
+-- qi 0 0 0 1
 qiZero :: QI
 qiZero = qi 0 0 0 1
 {-# INLINE qiZero #-}
 
--- | The constant one. @qi 1 0 0 1@
+-- | The constant one.
+--
+-- >>> qiOne
+-- qi 1 0 0 1
 qiOne :: QI
 qiOne  = qi 1 0 0 1
 {-# INLINE qiOne #-}
 
 -- | Check if the value is zero.
+--
+-- >>> map qiIsZero [qiZero, qiOne, qiSubR (qi 7 0 0 2) 3.5]
+-- [True,False,True]
 qiIsZero :: QI -> Bool
 -- If b = 0 then c = 0 and vice versa, guaranteed by the constructor.
 qiIsZero (unQI -> ~(a,b,_,_)) = a == 0 && b == 0
 {-# INLINE qiIsZero #-}
 
 -- | Convert a 'QI' number into a 'Floating' one.
+--
+-- >>> qiToFloat (qi 3 4 5 6) == ((3 + 4 * sqrt 5)/6 :: Double)
+-- True
 qiToFloat :: Floating a => QI -> a
 qiToFloat (unQI -> ~(a,b,c,d)) =
   (fromInteger a + fromInteger b * sqrt (fromInteger c)) / fromInteger d
 {-# INLINE qiToFloat #-}
 
 -- | Add an 'Integer' to a 'QI'.
+--
+-- >>> qi 3 4 5 6 `qiAddI` 1
+-- qi 9 4 5 6
 qiAddI :: QI -> Integer -> QI
-qiAddI n x = qiModify n $ \a b d ->
-  a `seq` b `seq` d `seq` x `seq` (a + d*x, b, d)
+qiAddI n x = over _qiABD go n
+  where go ~(a,b,d) = a `seq` b `seq` d `seq` x `seq` (a + d*x, b, d)
 {-# INLINE qiAddI #-}
 
 -- | Add a 'Rational' to a 'QI'.
+--
+-- >>> qi 3 4 5 6 `qiAddR` 1.2
+-- qi 51 20 5 30
 qiAddR :: QI -> Rational -> QI
-qiAddR n x = qiModify n $ \a b d ->
-  -- n = (a + b √c)/d + xN/xD
-  -- n = ((a + b √c) xD)/(d xD) + (d xN)/(d xD)
-  -- n = ((a xD + d xN) + b xD √c)/(d xD)
-  a `seq` b `seq` d `seq` xN `seq` xD `seq` (a*xD + d*xN, b*xD, d*xD)
-  where (xN, xD) = (numerator x, denominator x)
+qiAddR n x = over _qiABD go n
+  where
+    -- n = (a + b √c)/d + xN/xD
+    -- n = ((a + b √c) xD)/(d xD) + (d xN)/(d xD)
+    -- n = ((a xD + d xN) + b xD √c)/(d xD)
+    go ~(a,b,d) =
+      a `seq` b `seq` d `seq` xN `seq` xD `seq` (a*xD + d*xN, b*xD, d*xD)
+    (xN, xD) = (numerator x, denominator x)
 {-# INLINE qiAddR #-}
 
 -- | Subtract an 'Integer' from a 'QI'.
+--
+-- >>> qi 3 4 5 6 `qiSubI` 1
+-- qi (-3) 4 5 6
 qiSubI :: QI -> Integer -> QI
 qiSubI n x = qiAddI n (negate x)
 {-# INLINE qiSubI #-}
 
 -- | Subtract a 'Rational' from a 'QI'.
+--
+-- >>> qi 3 4 5 6 `qiSubR` 1.2
+-- qi (-21) 20 5 30
 qiSubR :: QI -> Rational -> QI
 qiSubR n x = qiAddR n (negate x)
 {-# INLINE qiSubR #-}
 
 -- | Multiply a 'QI' by an 'Integer'.
+--
+-- >>> qi 3 4 5 6 `qiMulI` 2
+-- qi 3 4 5 3
 qiMulI :: QI -> Integer -> QI
-qiMulI n x = qiModify n $ \a b d ->
-  a `seq` b `seq` d `seq` x `seq` (a*x, b*x, d)
+qiMulI n x = over _qiABD go n
+  where go ~(a,b,d) = a `seq` b `seq` d `seq` x `seq` (a*x, b*x, d)
 {-# INLINE qiMulI #-}
 
 -- | Multiply a 'QI' by a 'Rational'.
+--
+-- >>> qi 3 4 5 6 `qiMulR` 0.5
+-- qi 3 4 5 12
 qiMulR :: QI -> Rational -> QI
-qiMulR n x = qiModify n $ \a b d ->
-  -- n = (a + b √c)/d xN/xD
-  -- n = (a xN + b xN √c)/(d xD)
-  a `seq` b `seq` d `seq` xN `seq` xD `seq` (a*xN, b*xN, d*xD)
-  where (xN, xD) = (numerator x, denominator x)
+qiMulR n x = over _qiABD go n
+  where
+    -- n = (a + b √c)/d xN/xD
+    -- n = (a xN + b xN √c)/(d xD)
+    go ~(a,b,d) = a `seq` b `seq` d `seq` xN `seq` xD `seq` (a*xN, b*xN, d*xD)
+    (xN, xD) = (numerator x, denominator x)
 {-# INLINE qiMulR #-}
 
 -- | Divice a 'QI' by an 'Integer'.
+--
+-- >>> qi 3 4 5 6 `qiDivI` 2
+-- qi 3 4 5 12
 qiDivI :: QI -> Integer -> QI
-qiDivI n (nonZero "qiDivI" -> x) = qiModify n $ \a b d ->
-  a `seq` b `seq` d `seq` x `seq` (a, b, d*x)
+qiDivI n (nonZero "qiDivI" -> x) = over _qiABD go n
+  where go ~(a,b,d) = a `seq` b `seq` d `seq` x `seq` (a, b, d*x)
 {-# INLINE qiDivI #-}
 
 -- | Divice a 'QI' by a 'Rational'.
+--
+-- >>> qi 3 4 5 6 `qiDivR` 0.5
+-- qi 3 4 5 3
 qiDivR :: QI -> Rational -> QI
 qiDivR n (nonZero "qiDivR" -> x) = qiMulR n (recip x)
 {-# INLINE qiDivR #-}
 
 -- | Negate a 'QI'.
+--
+-- >>> qiNegate (qi 3 4 5 6)
+-- qi (-3) (-4) 5 6
 qiNegate :: QI -> QI
-qiNegate n = qiModify n $ \a b d ->
-  a `seq` b `seq` d `seq` (negate a, negate b, d)
+qiNegate n = over _qiABD go n
+  where go ~(a,b,d) = a `seq` b `seq` d `seq` (negate a, negate b, d)
 {-# INLINE qiNegate #-}
 
 -- | Compute the reciprocal of a 'QI'.
+--
+-- >>> qiRecip (qi 5 0 0 2)
+-- Just (qi 2 0 0 5)
+--
+-- >>> qiRecip (qi 0 1 5 2)
+-- Just (qi 0 2 5 5)
+--
+-- >>> qiRecip qiZero
+-- Nothing
 qiRecip :: QI -> Maybe QI
 qiRecip n@(unQI -> ~(a,b,c,d))
   -- 1/((a + b √c)/d)                       =
@@ -239,25 +422,55 @@
   -- (a d − b d √c) / (a² − b² c)
   | qiIsZero n = Nothing
   | denom == 0 = error ("qiRecip: Failed for " ++ show n)
-  | otherwise  = Just (qiModify n (\_ _ _ -> (a * d, negate (b * d), denom)))
+  | otherwise  = Just (set _qiABD (a * d, negate (b * d), denom) n)
   where denom = (a*a - b*b * c)
 
 -- | Add two 'QI's if the square root terms are the same or zeros.
+--
+-- >>> qi 3 4 5 6 `qiAdd` qiOne
+-- Just (qi 9 4 5 6)
+--
+-- >>> qi 3 4 5 6 `qiAdd` qi 3 4 5 6
+-- Just (qi 3 4 5 3)
+--
+-- >>> qi 0 1 5 1 `qiAdd` qi 0 1 6 1
+-- Nothing
 qiAdd :: QI -> QI -> Maybe QI
 qiAdd n@(unQI -> ~(a,b,c,d)) n'@(unQI -> ~(a',b',c',d'))
   -- n = (a + b √c)/d + (a' + b' √c')/d'
   -- n = ((a + b √c) d' + (a' + b' √c') d)/(d d')
   -- if c = c' then n = ((a d' + a' d) + (b d' + b' d) √c)/(d d')
-  | c  == 0   = Just (qiModify n' (\_ _ _ -> (a*d' + a'*d,        b'*d, d*d')))
-  | c' == 0   = Just (qiModify n  (\_ _ _ -> (a*d' + a'*d, b*d'       , d*d')))
-  | c  == c'  = Just (qiModify n  (\_ _ _ -> (a*d' + a'*d, b*d' + b'*d, d*d')))
+  | c  == 0   = Just (set _qiABD (a*d' + a'*d,        b'*d, d*d') n')
+  | c' == 0   = Just (set _qiABD (a*d' + a'*d, b*d'       , d*d') n)
+  | c  == c'  = Just (set _qiABD (a*d' + a'*d, b*d' + b'*d, d*d') n)
   | otherwise = Nothing
 
 -- | Subtract two 'QI's if the square root terms are the same or zeros.
+--
+-- >>> qi 3 4 5 6 `qiSub` qiOne
+-- Just (qi (-3) 4 5 6)
+--
+-- >>> qi 3 4 5 6 `qiSub` qi 3 4 5 6
+-- Just (qi 0 0 0 1)
+--
+-- >>> qi 0 1 5 1 `qiSub` qi 0 1 6 1
+-- Nothing
 qiSub :: QI -> QI -> Maybe QI
 qiSub n n' = qiAdd n (qiNegate n')
 
 -- | Multiply two 'QI's if the square root terms are the same or zeros.
+--
+-- >>> qi 3 4 5 6 `qiMul` qiZero
+-- Just (qi 0 0 0 1)
+--
+-- >>> qi 3 4 5 6 `qiMul` qiOne
+-- Just (qi 3 4 5 6)
+--
+-- >>> qi 3 4 5 6 `qiMul` qi 3 4 5 6
+-- Just (qi 89 24 5 36)
+--
+-- >>> qi 0 1 5 1 `qiMul` qi 0 1 6 1
+-- Nothing
 qiMul :: QI -> QI -> Maybe QI
 qiMul n@(unQI -> ~(a,b,c,d)) n'@(unQI -> ~(a',b',c',d'))
   -- n = (a + b √c)/d (a' + b' √c')/d'
@@ -265,16 +478,40 @@
   -- if c = 0  then n = (a a' + a b' √c')/(d d')
   -- if c' = 0 then n = (a a' + a' b √c)/(d d')
   -- if c = c' then n = ((a a' + b b' c) + (a b' + a' b) √c)/(d d')
-  | c  == 0   = Just (qiModify n' (\_ _ _ -> (a*a'         , a*b'       , d*d')))
-  | c' == 0   = Just (qiModify n  (\_ _ _ -> (a*a'         ,        a'*b, d*d')))
-  | c  == c'  = Just (qiModify n  (\_ _ _ -> (a*a' + b*b'*c, a*b' + a'*b, d*d')))
+  | c  == 0   = Just (set _qiABD (a*a'         , a*b'       , d*d') n')
+  | c' == 0   = Just (set _qiABD (a*a'         ,        a'*b, d*d') n)
+  | c  == c'  = Just (set _qiABD (a*a' + b*b'*c, a*b' + a'*b, d*d') n)
   | otherwise = Nothing
 
 -- | Divide two 'QI's if the square root terms are the same or zeros.
+--
+-- >>> qi 3 4 5 6 `qiDiv` qiZero
+-- Nothing
+--
+-- >>> qi 3 4 5 6 `qiDiv` qiOne
+-- Just (qi 3 4 5 6)
+--
+-- >>> qi 3 4 5 6 `qiDiv` qi 3 4 5 6
+-- Just (qi 1 0 0 1)
+--
+-- >>> qi 3 4 5 6 `qiDiv` qi 0 1 5 1
+-- Just (qi 20 3 5 30)
+--
+-- >>> qi 0 1 5 1 `qiDiv` qi 0 1 6 1
+-- Nothing
 qiDiv :: QI -> QI -> Maybe QI
 qiDiv n n' = qiMul n =<< qiRecip n'
 
 -- | Exponentiate a 'QI' to an 'Integer' power.
+--
+-- >>> qi 3 4 5 6 `qiPow` 0
+-- qi 1 0 0 1
+--
+-- >>> qi 3 4 5 6 `qiPow` 1
+-- qi 3 4 5 6
+--
+-- >>> qi 3 4 5 6 `qiPow` 2
+-- qi 89 24 5 36
 qiPow :: QI -> Integer -> QI
 qiPow num (nonNegative "qiPow" -> pow) = go num pow
   where
@@ -295,6 +532,15 @@
     sudoQIMul n n' = case qiMul n n' of ~(Just m) -> m
 
 -- | Compute the floor of a 'QI'.
+--
+-- >>> qiFloor (qi 10 0 0 2)
+-- 5
+--
+-- >>> qiFloor (qi 10 2 2 2)
+-- 6
+--
+-- >>> qiFloor (qi 10 2 5 2)
+-- 7
 qiFloor :: QI -> Integer
 qiFloor (unQI -> ~(a,b,c,d)) =
   -- n = (a + b √c)/d
@@ -307,6 +553,19 @@
     ~(b2cLow, b2cHigh) = iSqrtBounds (b*b * c)
 
 -- | Convert a (possibly periodic) continued fraction to a 'QI'.
+--
+-- @[2; 2] = 2 + 1\/2 = 5\/2@.
+--
+-- >>> continuedFractionToQI (2,NonCyc [2])
+-- qi 5 0 0 2
+--
+-- @[2; 1, 1, 1, 4, 1, 1, 1, 4, …] = √7@.
+--
+-- >>> continuedFractionToQI (2,Cyc [] 1 [1,1,4])
+-- qi 0 1 7 1
+--
+-- >>> continuedFractionToQI (0,Cyc [83,78,65,75,69] 32 [66,65,68,71,69,82])
+-- qi 987601513930253257378987883 1 14116473325908285531353005 81983584717737887813195873886
 continuedFractionToQI :: (Integer, CycList Integer) -> QI
 continuedFractionToQI (i0_, is_) = qiAddI (go is_) i0_
   where
@@ -350,6 +609,19 @@
     pos = positive "continuedFractionToQI"
 
 -- | Convert a 'QI' into a (possibly periodic) continued fraction.
+--
+-- @5\/2 = 2 + 1\/2 = [2; 2]@.
+--
+-- >>> qiToContinuedFraction (qi 5 0 0 2)
+-- (2,NonCyc [2])
+--
+-- @√7 = [2; 1, 1, 1, 4, 1, 1, 1, 4, …]@.
+--
+-- >>> qiToContinuedFraction (qi 0 1 7 1)
+-- (2,Cyc [] 1 [1,1,4])
+--
+-- >>> qiToContinuedFraction (qi 987601513930253257378987883 1 14116473325908285531353005 81983584717737887813195873886)
+-- (0,Cyc [83,78,65,75,69] 32 [66,65,68,71,69,82])
 qiToContinuedFraction :: QI
                       -> (Integer, CycList Integer)
 qiToContinuedFraction num
diff --git a/src/Numeric/QuadraticIrrational/CyclicList.hs b/src/Numeric/QuadraticIrrational/CyclicList.hs
--- a/src/Numeric/QuadraticIrrational/CyclicList.hs
+++ b/src/Numeric/QuadraticIrrational/CyclicList.hs
@@ -14,6 +14,16 @@
 import Data.Foldable
 import Data.Monoid
 
+-- $setup
+-- import Data.Foldable (toList)
+
+-- | A container for a possibly cyclic list.
+--
+-- >>> toList (NonCyc "hello")
+-- "hello"
+--
+-- >>> take 70 (toList (Cyc "prefix " 'c' "ycle"))
+-- "prefix cyclecyclecyclecyclecyclecyclecyclecyclecyclecyclecyclecyclecyc"
 data CycList a = NonCyc [a]  -- ^ A non-cyclic list.
                | Cyc [a] a [a]
                  -- ^ A non-cyclic list followed by the head of a cyclic list
diff --git a/src/Numeric/QuadraticIrrational/Internal/Lens.hs b/src/Numeric/QuadraticIrrational/Internal/Lens.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/QuadraticIrrational/Internal/Lens.hs
@@ -0,0 +1,42 @@
+{-# LANGUAGE Rank2Types #-}
+
+-- |
+-- Module      : Numeric.QuadraticIrrational.Internal.Lens
+-- Description : A tiny implementation of some lens primitives
+-- Copyright   : © 2014 Johan Kiviniemi
+-- License     : MIT
+-- Maintainer  : Johan Kiviniemi <devel@johan.kiviniemi.name>
+-- Stability   : provisional
+-- Portability : Rank2Types
+--
+-- A tiny implementation of some lens primitives. Please see
+-- <http://hackage.haskell.org/package/lens> for proper documentation.
+
+module Numeric.QuadraticIrrational.Internal.Lens
+  ( Lens, Traversal, Lens', Traversal', Getting, Setting
+  , view, over, set
+  ) where
+
+import Control.Applicative
+import Data.Functor.Identity
+
+type Lens      s t a b = Functor     f => (a -> f b) -> s -> f t
+type Traversal s t a b = Applicative f => (a -> f b) -> s -> f t
+
+type Lens'      s a = Functor     f => (a -> f a) -> s -> f s
+type Traversal' s a = Applicative f => (a -> f a) -> s -> f s
+
+type Getting r s a   = (a -> Const r a)  -> s -> Const r s
+type Setting s t a b = (a -> Identity b) -> s -> Identity t
+
+view :: Getting a s a -> s -> a
+view l s = getConst (l Const s)
+{-# INLINE view #-}
+
+over :: Setting s t a b -> (a -> b) -> s -> t
+over l f s = runIdentity (l (f `seq` Identity . f) s)
+{-# INLINE over #-}
+
+set :: Setting s t a b -> b -> s -> t
+set l b s = over l (const b) s
+{-# INLINE set #-}
diff --git a/tests/Main.hs b/tests/Main.hs
deleted file mode 100644
--- a/tests/Main.hs
+++ /dev/null
@@ -1,16 +0,0 @@
-module Main (main) where
-
-import Test.Tasty
-
-import qualified CyclicList
-import qualified QuadraticIrrational
-
-main :: IO ()
-main = defaultMain tests
-
-tests :: TestTree
-tests =
-  testGroup "quadratic-irrational"
-    [ CyclicList.tests
-    , QuadraticIrrational.tests
-    ]
diff --git a/tests/QuadraticIrrational.hs b/tests/QuadraticIrrational.hs
--- a/tests/QuadraticIrrational.hs
+++ b/tests/QuadraticIrrational.hs
@@ -10,6 +10,7 @@
 import Test.Tasty.QuickCheck
 
 import Numeric.QuadraticIrrational
+import Numeric.QuadraticIrrational.Internal.Lens
 
 -- Slow but precise.
 type RefFloat = CReal
@@ -44,19 +45,48 @@
       , testProperty "qi'/runQI'" $ \a b (NonNegative c) ->
           runQI' (qi' a b c) $ \a' b' c' ->
             approxEq' (approxQI' a b c) (approxQI' a' b' c')
+      ]
+    , testGroup "Lenses"
+      [ testProperty "_qi" $ \n a' b' (NonNegative c') (NonZero d') ->
+          let n'  = over _qi (\(a,b,c,d) -> (a+a',b-b',c*c',d*d')) n
+              n'' = runQI n $ \a b c d -> qi (a+a') (b-b') (c*c') (d*d')
+          in  approxEq (qiToFloat n') (qiToFloat n'')
 
-      , testProperty "qiModify" $ \n a' b' (NonZero d') ->
-          runQI n $ \a b c d ->
-            approxEq' (qiToFloat (qiModify n (\a_ b_ d_ ->
-                                                (a_+a', b_-b', d_*d'))))
-                      (qiToFloat (qi (a+a') (b-b') c (d*d')))
+      , testProperty "_qi'" $ \n a' b' (NonNegative c') ->
+          let n'  = over _qi' (\(a,b,c) -> (a+a',b-b',c*c')) n
+              n'' = runQI' n $ \a b c -> qi' (a+a') (b-b') (c*c')
+          in  approxEq (qiToFloat n') (qiToFloat n'')
 
-      , testProperty "qiToFloat" $ \a b (NonNegative c) (NonZero d) ->
-          approxEq' (qiToFloat (qi a b c d)) (approxQI a b c d)
-      ]
+      , testProperty "_qiABD" $ \n a' b' (NonZero d') ->
+          let n'  = over _qiABD (\(a,b,d) -> (a+a',b-b',d*d')) n
+              n'' = runQI n $ \a b c d -> qi (a+a') (b-b') c (d*d')
+          in  approxEq (qiToFloat n') (qiToFloat n'')
 
+      , testProperty "_qiA" $ \n a' ->
+          let n'  = over _qiA (+ a') n
+              n'' = runQI n $ \a b c d -> qi (a+a') b c d
+          in  approxEq (qiToFloat n') (qiToFloat n'')
+
+      , testProperty "_qiB" $ \n b' ->
+          let n'  = over _qiB (+ b') n
+              n'' = runQI n $ \a b c d -> qi a (b+b') c d
+          in  approxEq (qiToFloat n') (qiToFloat n'')
+
+      , testProperty "_qiC" $ \n (NonNegative c') ->
+          let n'  = over _qiC (* c') n
+              n'' = runQI n $ \a b c d -> qi a b (c*c') d
+          in  approxEq (qiToFloat n') (qiToFloat n'')
+
+      , testProperty "_qiD" $ \n (NonZero d') ->
+          let n'  = over _qiD (* d') n
+              n'' = runQI n $ \a b c d -> qi a b c (d*d')
+          in  approxEq (qiToFloat n') (qiToFloat n'')
+      ]
     , testGroup "Numerical operations"
-      [ testProperty "qiAddI" $ \n x ->
+      [ testProperty "qiToFloat" $ \a b (NonNegative c) (NonZero d) ->
+          approxEq' (qiToFloat (qi a b c d)) (approxQI a b c d)
+
+      , testProperty "qiAddI" $ \n x ->
           approxEq' (qiToFloat (qiAddI n x)) (qiToFloat n + fromInteger x)
 
       , testProperty "qiSubI" $ \n x ->
@@ -116,8 +146,9 @@
 
       , testProperty "qiFloor" $ \n ->
           qiFloor n === floor (qiToFloat n :: RefFloat)
-
-      , testProperty "qiToContinuedFraction/continuedFractionToQI" $ \n ->
+      ]
+    , testGroup "Continued fractions"
+      [ testProperty "qiToContinuedFraction/continuedFractionToQI" $ \n ->
           let cf  = qiToContinuedFraction n
               len = case cf of
                       (_, NonCyc _)   -> 0
diff --git a/tests/doctests.hs b/tests/doctests.hs
new file mode 100644
--- /dev/null
+++ b/tests/doctests.hs
@@ -0,0 +1,32 @@
+{-# LANGUAGE MultiWayIf #-}
+
+module Main (main) where
+
+import Control.Applicative
+import Control.Monad
+import Control.Monad.List
+import Data.List
+import System.Directory
+import System.FilePath
+import Test.DocTest
+
+main :: IO ()
+main = do
+  sources <- findSources "src"
+  doctest ("-isrc" : sources)
+
+findSources :: FilePath -> IO [FilePath]
+findSources dir = runListT (goDir dir)
+  where
+    goItem :: FilePath -> FilePath -> ListT IO FilePath
+    goItem _ ('.':_) = empty
+    goItem parent name = do
+      let path = parent </> name
+      isDir  <- liftIO (doesDirectoryExist path)
+      isFile <- liftIO (doesFileExist path)
+      if | isDir     -> goDir  path
+         | isFile    -> goFile path
+         | otherwise -> empty
+
+    goDir  path = goItem path =<< ListT (getDirectoryContents path)
+    goFile path = path <$ guard (".hs" `isSuffixOf` path)
diff --git a/tests/tasty.hs b/tests/tasty.hs
new file mode 100644
--- /dev/null
+++ b/tests/tasty.hs
@@ -0,0 +1,16 @@
+module Main (main) where
+
+import Test.Tasty
+
+import qualified CyclicList
+import qualified QuadraticIrrational
+
+main :: IO ()
+main = defaultMain tests
+
+tests :: TestTree
+tests =
+  testGroup "quadratic-irrational"
+    [ CyclicList.tests
+    , QuadraticIrrational.tests
+    ]