diff --git a/.gitignore b/.gitignore
new file mode 100644
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,3 @@
+.cabal-sandbox
+cabal.sandbox.config
+dist
diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,21 @@
+The MIT License (MIT)
+
+Copyright © 2014 Johan Kiviniemi
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
diff --git a/README.md b/README.md
new file mode 100644
--- /dev/null
+++ b/README.md
@@ -0,0 +1,7 @@
+# `quadratic-irrational`
+
+An implementation of [quadratic irrationals][qi] with support for conversion
+from and to [periodic continued fractions][pcf].
+
+[qi]:  http://en.wikipedia.org/wiki/Quadratic_irrational
+[pcf]: http://en.wikipedia.org/wiki/Periodic_continued_fraction
diff --git a/Setup.hs b/Setup.hs
new file mode 100644
--- /dev/null
+++ b/Setup.hs
@@ -0,0 +1,2 @@
+import Distribution.Simple
+main = defaultMain
diff --git a/quadratic-irrational.cabal b/quadratic-irrational.cabal
new file mode 100644
--- /dev/null
+++ b/quadratic-irrational.cabal
@@ -0,0 +1,55 @@
+name: quadratic-irrational
+category: Math, Algorithms, Data
+version: 0.0.1
+license: MIT
+license-file: LICENSE
+author: Johan Kiviniemi <devel@johan.kiviniemi.name>
+maintainer: Johan Kiviniemi <devel@johan.kiviniemi.name>
+stability: provisional
+homepage: https://github.com/ion1/quadratic-irrational
+bug-reports: https://github.com/ion1/quadratic-irrational/issues
+copyright: Copyright © 2014 Johan Kiviniemi
+synopsis: An implementation of quadratic irrationals
+description:
+  An implementation of
+  <http://en.wikipedia.org/wiki/Quadratic_irrational quadratic irrationals>
+  with support for conversion from and to
+  <http://en.wikipedia.org/wiki/Periodic_continued_fraction periodic continued fractions>.
+tested-with: GHC == 7.6.3
+
+build-type: Simple
+cabal-version: >= 1.10
+extra-source-files:
+  .gitignore
+  README.md
+
+source-repository head
+  type: git
+  location: https://github.com/ion1/quadratic-irrational.git
+
+library
+  exposed-modules: Numeric.QuadraticIrrational
+                 , Numeric.QuadraticIrrational.CyclicList
+  hs-source-dirs: src
+  build-depends: base >= 4.6 && < 4.8
+               , arithmoi == 0.4.*
+               , containers == 0.5.*
+               , mtl == 2.1.*
+               , transformers == 0.3.*
+  default-language: Haskell2010
+  ghc-options: -Wall -O2 -funbox-strict-fields
+
+test-suite test-quadratic-irrational
+  type: exitcode-stdio-1.0
+  main-is: Main.hs
+  other-modules: QuadraticIrrational
+               , CyclicList
+  hs-source-dirs: tests
+  build-depends: base
+               , quadratic-irrational
+               , numbers == 3000.*
+               , QuickCheck >= 2.7 && < 3
+               , tasty == 0.8.*
+               , tasty-quickcheck == 0.8.*
+  default-language: Haskell2010
+  ghc-options: -Wall -O2 -funbox-strict-fields
diff --git a/src/Numeric/QuadraticIrrational.hs b/src/Numeric/QuadraticIrrational.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/QuadraticIrrational.hs
@@ -0,0 +1,412 @@
+{-# LANGUAGE ViewPatterns #-}
+
+-- |
+-- Module      : Numeric.QuadraticIrrational
+-- Description : An implementation of quadratic irrationals
+-- Copyright   : © 2014 Johan Kiviniemi
+-- License     : MIT
+-- Maintainer  : Johan Kiviniemi <devel@johan.kiviniemi.name>
+-- Stability   : provisional
+-- Portability : ViewPatterns
+--
+-- An implementation of
+-- <http://en.wikipedia.org/wiki/Quadratic_irrational quadratic irrationals>
+-- with support for conversion from and to
+-- <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
+  , qiToFloat
+  , qiAddI, qiSubI, qiMulI, qiDivI
+  , qiAddR, qiSubR, qiMulR, qiDivR
+  , qiNegate, qiRecip, qiAdd, qiSub, qiMul, qiDiv, qiPow
+  , qiFloor, continuedFractionToQI, qiToContinuedFraction
+  , module Numeric.QuadraticIrrational.CyclicList
+  ) where
+
+import Control.Applicative
+import Control.Monad.State
+import Data.List
+import Data.Maybe
+import Data.Ratio
+import qualified Data.Set as Set
+import Math.NumberTheory.Powers.Squares
+import Math.NumberTheory.Primes.Factorisation
+import Text.Read
+
+import Numeric.QuadraticIrrational.CyclicList
+
+-- | @(a + b √c) / d@
+data QI = QI !Integer
+             !Integer
+             !Integer
+             !Integer
+  deriving (Eq)
+
+instance Show QI where
+  showsPrec p (QI a b c d) = showParen (p > 10)
+                           $ showString "qi " . showsPrec 11 a
+                           . showChar   ' '   . showsPrec 11 b
+                           . showChar   ' '   . showsPrec 11 c
+                           . showChar   ' '   . showsPrec 11 d
+
+instance Read QI where
+  readPrec = parens rQI
+    where
+      rQI = prec 10 $ do
+        Ident "qi" <- lexP
+        qi <$> step readPrec <*> step readPrec <*> step readPrec
+           <*> step readPrec
+
+  readListPrec = readListPrecDefault
+
+type QITuple = (Integer, Integer, Integer, Integer)
+
+-- | Given @a@, @b@, @c@ and @d@ such that @n = (a + b √c)/d@, constuct a 'QI'
+-- corresponding to @n@.
+qi :: Integer  -- ^ a
+   -> Integer  -- ^ b
+   -> Integer  -- ^ c
+   -> Integer  -- ^ d
+   -> QI
+qi a b (nonNegative "qi" -> c) (nonZero "qi" -> d)
+  | b == 0    = reduceCons a 0 0 d
+  | c == 0    = reduceCons a 0 0 d
+  | c == 1    = reduceCons (a + b) 0 0 d
+  | otherwise = simplifyReduceCons a b c d
+{-# INLINE qi #-}
+
+-- Simplify @b √c@ before constructing a 'QI'.
+simplifyReduceCons :: Integer -> Integer -> Integer -> Integer -> QI
+simplifyReduceCons a b (nonZero "simplifyReduceCons" -> c) d
+  | c' == 1   = reduceCons (a + b') 0 0 d
+  | otherwise = reduceCons a b' c' d
+  where ~(b', c') = separateSquareFactors b c
+{-# INLINE simplifyReduceCons #-}
+
+-- | Given @b@ and @c@ such that @n = b √c@, return a potentially simplified
+-- @(b, c)@.
+separateSquareFactors :: Integer -> Integer -> (Integer, Integer)
+separateSquareFactors b (nonNegative "separateSquareFactors" -> c) =
+  case foldl' go (1,1) (factorise c) of
+    ~(bMul, c') -> (b*bMul, c')
+  where
+    go :: (Integer, Integer) -> (Integer, Int) -> (Integer, Integer)
+    go ~(i, j) ~(fac, pow) =
+      i `seq` j `seq` fac `seq` pow `seq`
+        if even pow
+          then (i*fac^(pow     `div` 2), j)
+          else (i*fac^((pow-1) `div` 2), j*fac)
+
+-- Reduce the @a@, @b@, @d@ factors before constructing a 'QI'.
+reduceCons :: Integer -> Integer -> Integer -> Integer -> QI
+reduceCons a b c (nonZero "reduceCons" -> d) =
+  QI (a `quot` q) (b `quot` q) c (d `quot` q)
+  where q = signum d * (a `gcd` b `gcd` d)
+{-# INLINE reduceCons #-}
+
+-- | Given @a@, @b@ and @c@ such that @n = a + b √c@, constuct a 'QI'
+-- corresponding to @n@.
+qi' :: Rational  -- ^ a
+    -> Rational  -- ^ b
+    -> Integer   -- ^ c
+    -> QI
+qi' a b (nonNegative "qi'" -> c) = n
+  where
+    -- (aN/aD) + (bN/bD) √c = ((aN bD) + (bN aD) √c) / (aD bD)
+    n = qi (aN*bD) (bN*aD) c (aD*bD)
+    (aN, aD) = (numerator a, denominator a)
+    (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 -> (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 -> (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 -> (Integer, Integer, Integer, Integer)
+unQI n = runQI n (,,,)
+{-# INLINE unQI #-}
+
+-- | Given @n@ such that @n = a + b √c@, return @(a, b, c)@.
+unQI' :: QI -> (Rational, Rational, Integer)
+unQI' n = runQI' n (,,)
+{-# INLINE unQI' #-}
+
+-- | The constant zero. @qi 0 0 0 1@
+qiZero :: QI
+qiZero = qi 0 0 0 1
+{-# INLINE qiZero #-}
+
+-- | The constant one. @qi 1 0 0 1@
+qiOne :: QI
+qiOne  = qi 1 0 0 1
+{-# INLINE qiOne #-}
+
+-- | Check if the value is zero.
+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 :: 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'.
+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)
+{-# INLINE qiAddI #-}
+
+-- | Add a 'Rational' to a 'QI'.
+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)
+{-# INLINE qiAddR #-}
+
+-- | Subtract an 'Integer' from a 'QI'.
+qiSubI :: QI -> Integer -> QI
+qiSubI n x = qiAddI n (negate x)
+{-# INLINE qiSubI #-}
+
+-- | Subtract a 'Rational' from a 'QI'.
+qiSubR :: QI -> Rational -> QI
+qiSubR n x = qiAddR n (negate x)
+{-# INLINE qiSubR #-}
+
+-- | Multiply a 'QI' by an 'Integer'.
+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)
+{-# INLINE qiMulI #-}
+
+-- | Multiply a 'QI' by a 'Rational'.
+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)
+{-# INLINE qiMulR #-}
+
+-- | Divice a 'QI' by an 'Integer'.
+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)
+{-# INLINE qiDivI #-}
+
+-- | Divice a 'QI' by a 'Rational'.
+qiDivR :: QI -> Rational -> QI
+qiDivR n (nonZero "qiDivR" -> x) = qiMulR n (recip x)
+{-# INLINE qiDivR #-}
+
+-- | Negate a 'QI'.
+qiNegate :: QI -> QI
+qiNegate n = qiModify n $ \a b d ->
+  a `seq` b `seq` d `seq` (negate a, negate b, d)
+{-# INLINE qiNegate #-}
+
+-- | Compute the reciprocal of a 'QI'.
+qiRecip :: QI -> Maybe QI
+qiRecip n@(unQI -> ~(a,b,c,d))
+  -- 1/((a + b √c)/d)                       =
+  -- d/(a + b √c)                           =
+  -- d (a − b √c) / ((a + b √c) (a − b √c)) =
+  -- d (a − b √c) / (a² − b² c)             =
+  -- (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)))
+  where denom = (a*a - b*b * c)
+
+-- | Add two 'QI's if the square root terms are the same or zeros.
+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')))
+  | otherwise = Nothing
+
+-- | Subtract two 'QI's if the square root terms are the same or zeros.
+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.
+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'
+  -- n = (a a' + a b' √c' + a' b √c + b b' √c √c')/(d d')
+  -- 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')))
+  | otherwise = Nothing
+
+-- | Divide two 'QI's if the square root terms are the same or zeros.
+qiDiv :: QI -> QI -> Maybe QI
+qiDiv n n' = qiMul n =<< qiRecip n'
+
+-- | Exponentiate a 'QI' to an 'Integer' power.
+qiPow :: QI -> Integer -> QI
+qiPow num (nonNegative "qiPow" -> pow) = go num pow
+  where
+    go _ 0 = qiOne
+    go n 1 = n
+    go n p
+      | even p    = go  (sudoQIMul n n) (p     `div` 2)
+      | otherwise = go' (sudoQIMul n n) ((p-1) `div` 2) n
+
+    -- Like go but multiplied with n'.
+    go' _ 0 n' = n'
+    go' n 1 n' = sudoQIMul n n'
+    go' n p n'
+      | even p    = go' (sudoQIMul n n) (p     `div` 2) n'
+      | otherwise = go' (sudoQIMul n n) ((p-1) `div` 2) (sudoQIMul n n')
+
+    -- Multiplying a QI with its own power will always succeed.
+    sudoQIMul n n' = case qiMul n n' of ~(Just m) -> m
+
+-- | Compute the floor of a 'QI'.
+qiFloor :: QI -> Integer
+qiFloor (unQI -> ~(a,b,c,d)) =
+  -- n = (a + b √c)/d
+  -- n d = a + b √c
+  -- n d = a + signum b · √(b² c)
+  n_d `div` d
+  where
+    n_d = a + min (signum b * b2cLow) (signum b * b2cHigh)
+
+    ~(b2cLow, b2cHigh) = iSqrtBounds (b*b * c)
+
+-- | Convert a (possibly periodic) continued fraction to a 'QI'.
+continuedFractionToQI :: (Integer, CycList Integer) -> QI
+continuedFractionToQI (i0_, is_) = qiAddI (go is_) i0_
+  where
+    go (NonCyc as)   = goNonCyc as qiZero
+    go (Cyc as b bs) = goNonCyc as (goCyc (b:bs))
+
+    goNonCyc ((pos -> i):is) final = sudoQIRecip (qiAddI (goNonCyc is final) i)
+    goNonCyc []              final = final
+
+    goCyc is = sudoQIRecip (solvePeriodic is)
+
+    -- x = (a x + b) / (c x + d)
+    -- x (c x + d) = a x + b
+    -- c x² + d x = a x + b
+    -- c x² + (d − a) x − b = 0
+    -- Apply quadratic formula, positive solution only.
+    solvePeriodic is =
+      case solvePeriodic' is of
+        ~(a,b,c,d) ->
+          a `seq` b `seq` c `seq` d `seq`
+            qfPos c (d - a) (negate b)
+      where
+        qfPos i j k = qi (negate j) 1 (j*j - 4*i*k) (2*i)
+
+    -- i + 1/((a x + b) / (c x + d))      =
+    -- i + (c x + d)/(a x + b)            =
+    -- ((a i x + b i + c x + d)/(a x + b) =
+    -- ((a i + c) x + (b i + d))/(a x + b)
+    solvePeriodic' ((pos -> i):is) =
+      case solvePeriodic' is of
+        ~(a,b,c,d) ->
+          a `seq` b `seq` c `seq` d `seq` i `seq`
+            (a*i+c, b*i+d, a, b)
+
+    -- x = (1 x + 0) / (0 x + 1)
+    solvePeriodic' [] = (1,0,0,1)
+
+    sudoQIRecip n =
+      fromMaybe (error "continuedFractionToQI: Divide by zero") (qiRecip n)
+
+    pos = positive "continuedFractionToQI"
+
+-- | Convert a 'QI' into a (possibly periodic) continued fraction.
+qiToContinuedFraction :: QI
+                      -> (Integer, CycList Integer)
+qiToContinuedFraction num
+  | Just isLoopQI <- loopQI =
+      case break isLoopQI cfs of
+        (preLoop, ~(i:postLoop)) ->
+          let is = takeWhile (not . isLoopQI) postLoop
+          in  (i0, Cyc (map snd preLoop) (snd i) (map snd is))
+  | otherwise =
+      (i0, NonCyc (map snd cfs))
+  where
+    (i0, cfs) = qiToContinuedFractionList num
+
+    loopQI :: Maybe ((QITuple,a) -> Bool)
+    loopQI = evalState (go cfs) Set.empty
+      where
+        go ((n,_) : xs) = do
+          haveSeen <- gets (Set.member n)
+          modify (Set.insert n)
+          if haveSeen
+            then return (Just ((== n) . fst))
+            else go xs
+        go [] = return Nothing
+
+qiToContinuedFractionList :: QI -> (Integer, [(QITuple, Integer)])
+qiToContinuedFractionList num =
+  case go (Just num) of
+    -- There is always a first number.
+    ~((_,i) : xs) -> (i, xs)
+  where
+    go (Just n) = (unQI n, i) : go (qiRecip (qiSubI n i))
+      where i = qiFloor n
+    go Nothing  = []
+
+iSqrtBounds :: Integer -> (Integer, Integer)
+iSqrtBounds n = (low, high)
+  where
+    low = integerSquareRoot n
+    high | low*low == n = low
+         | otherwise    = low + 1
+
+nonNegative :: (Num a, Ord a, Show a) => String -> a -> a
+nonNegative name = validate name "non-negative" (>= 0)
+{-# INLINE nonNegative #-}
+
+positive :: (Num a, Ord a, Show a) => String -> a -> a
+positive name = validate name "positive" (> 0)
+{-# INLINE positive #-}
+
+nonZero :: (Num a, Eq a, Show a) => String -> a -> a
+nonZero name = validate name "non-zero" (/= 0)
+{-# INLINE nonZero #-}
+
+validate :: Show a => String -> String -> (a -> Bool) -> a -> a
+validate name expected f a
+  | f a = a
+  | otherwise =
+      error ("Numeric.QuadraticIrrational." ++ name ++ ": Got " ++ show a
+              ++ ", expected " ++ expected)
+{-# INLINE validate #-}
diff --git a/src/Numeric/QuadraticIrrational/CyclicList.hs b/src/Numeric/QuadraticIrrational/CyclicList.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/QuadraticIrrational/CyclicList.hs
@@ -0,0 +1,29 @@
+-- |
+-- Module      : Numeric.QuadraticIrrational.CyclicList
+-- Description : A container for a possibly cyclic list.
+-- Copyright   : © 2014 Johan Kiviniemi
+-- License     : MIT
+-- Maintainer  : Johan Kiviniemi <devel@johan.kiviniemi.name>
+-- Stability   : provisional
+-- Portability : portable
+
+module Numeric.QuadraticIrrational.CyclicList
+  ( CycList (..)
+  ) where
+
+import Data.Foldable
+import Data.Monoid
+
+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
+                 -- followed by the tail of the cyclic list.
+  deriving (Eq, Ord, Read, Show)
+
+instance Functor CycList where
+  fmap f (NonCyc as) = NonCyc (fmap f as)
+  fmap f (Cyc as b bs) = Cyc (fmap f as) (f b) (fmap f bs)
+
+instance Foldable CycList where
+  foldMap f (NonCyc as)   = foldMap f as
+  foldMap f (Cyc as b bs) = foldMap f as <> foldMap f (cycle (b:bs))
diff --git a/tests/CyclicList.hs b/tests/CyclicList.hs
new file mode 100644
--- /dev/null
+++ b/tests/CyclicList.hs
@@ -0,0 +1,41 @@
+{-# OPTIONS_GHC -fno-warn-orphans #-}
+
+module CyclicList (tests) where
+
+import Control.Applicative
+import qualified Data.Foldable as F
+import Test.Tasty
+import Test.Tasty.QuickCheck
+
+import Numeric.QuadraticIrrational.CyclicList
+
+instance Arbitrary a => Arbitrary (CycList a) where
+  arbitrary = oneof [ NonCyc <$> arbitrary
+                    , Cyc <$> arbitrary <*> arbitrary <*> arbitrary
+                    ]
+
+  shrink (NonCyc as)   = [ NonCyc as'   | as' <- shrink as ]
+  shrink (Cyc as b bs) = [ Cyc as' b bs | as' <- shrink as ]
+                      ++ [ Cyc as b' bs | b'  <- shrink b  ]
+                      ++ [ Cyc as b bs' | bs' <- shrink bs ]
+
+tests :: TestTree
+tests =
+  testGroup "CyclicList"
+    [ testProperty "fmap" . withListEquiv $ \asC asL ->
+        initEq' (fmap (*10) asC) (fmap (*10) asL)
+    , testProperty "toList" . withListEquiv $ \asC asL ->
+        take 1000 (F.toList asC) === take 1000 asL
+    ]
+
+withListEquiv :: (CycList Integer -> [Integer] -> b) -> CycList Integer -> b
+withListEquiv f cl@(NonCyc as)   = f cl as
+withListEquiv f cl@(Cyc as b bs) = f cl (as ++ cycle (b:bs))
+
+initEq :: Eq a => CycList a -> [a] -> Bool
+initEq (NonCyc as)   cs = take 1000 cs == take 1000 as
+initEq (Cyc as b bs) cs = take 1000 cs == take 1000 (as ++ cycle (b:bs))
+
+initEq' :: (Eq a, Show a) => CycList a -> [a] -> Property
+initEq' (NonCyc as)   cs = take 1000 cs === take 1000 as
+initEq' (Cyc as b bs) cs = take 1000 cs === take 1000 (as ++ cycle (b:bs))
diff --git a/tests/Main.hs b/tests/Main.hs
new file mode 100644
--- /dev/null
+++ b/tests/Main.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
+    ]
diff --git a/tests/QuadraticIrrational.hs b/tests/QuadraticIrrational.hs
new file mode 100644
--- /dev/null
+++ b/tests/QuadraticIrrational.hs
@@ -0,0 +1,158 @@
+{-# LANGUAGE ViewPatterns #-}
+
+{-# OPTIONS_GHC -fno-warn-orphans #-}
+
+module QuadraticIrrational (tests) where
+
+import Control.Applicative
+import Data.Number.CReal
+import Test.Tasty
+import Test.Tasty.QuickCheck
+
+import Numeric.QuadraticIrrational
+
+-- Slow but precise.
+type RefFloat = CReal
+
+instance Arbitrary QI where
+  arbitrary = consQI <$> arbitrary <*> arbitrary <*> arbitrary <*> arbitrary
+    where
+     consQI a b (NonNegative c) (NonZero d) = qi a b c d
+
+  shrink (unQI -> ~(a,b,c,d)) =
+    [ qi a' b  c  d  | a' <- shrink a ] ++
+    [ qi a  b' c  d  | b' <- shrink b ] ++
+    [ qi a  b  c' d  | NonNegative c' <- shrink (NonNegative c) ] ++
+    [ qi a  b  c  d' | NonZero     d' <- shrink (NonZero     d) ]
+
+tests :: TestTree
+tests =
+  testGroup "QuadraticIrrational"
+    [ testGroup "Construction/destruction/conversion"
+      [ testProperty "qi/runQI" $ \a b (NonNegative c) (NonZero d) ->
+          runQI (qi a b c d) $ \a' b' c' d' ->
+            approxEq' (approxQI a b c d) (approxQI a' b' c' d')
+
+      , testProperty "qi/runQI'" $ \a b (NonNegative c) (NonZero d) ->
+          runQI' (qi a b c d) $ \a' b' c' ->
+            approxEq' (approxQI a b c d) (approxQI' a' b' c')
+
+      , testProperty "qi'/runQI" $ \a b (NonNegative c) ->
+          runQI (qi' a b c) $ \a' b' c' d' ->
+            approxEq' (approxQI' a b c) (approxQI a' b' c' d')
+
+      , testProperty "qi'/runQI'" $ \a b (NonNegative c) ->
+          runQI' (qi' a b c) $ \a' b' c' ->
+            approxEq' (approxQI' a b c) (approxQI' a' b' c')
+
+      , 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 "qiToFloat" $ \a b (NonNegative c) (NonZero d) ->
+          approxEq' (qiToFloat (qi a b c d)) (approxQI a b c d)
+      ]
+
+    , testGroup "Numerical operations"
+      [ testProperty "qiAddI" $ \n x ->
+          approxEq' (qiToFloat (qiAddI n x)) (qiToFloat n + fromInteger x)
+
+      , testProperty "qiSubI" $ \n x ->
+          approxEq' (qiToFloat (qiSubI n x)) (qiToFloat n - fromInteger x)
+
+      , testProperty "qiMulI" $ \n x ->
+          approxEq' (qiToFloat (qiMulI n x)) (qiToFloat n * fromInteger x)
+
+      , testProperty "qiDivI" $ \n x ->
+          x /= 0 ==>
+            approxEq' (qiToFloat (qiDivI n x)) (qiToFloat n / fromInteger x)
+
+      , testProperty "qiAddR" $ \n x ->
+          approxEq' (qiToFloat (qiAddR n x)) (qiToFloat n + fromRational x)
+
+      , testProperty "qiSubR" $ \n x ->
+          approxEq' (qiToFloat (qiSubR n x)) (qiToFloat n - fromRational x)
+
+      , testProperty "qiMulR" $ \n x ->
+          approxEq' (qiToFloat (qiMulR n x)) (qiToFloat n * fromRational x)
+
+      , testProperty "qiDivR" $ \n x ->
+          x /= 0 ==>
+            approxEq' (qiToFloat (qiDivR n x)) (qiToFloat n / fromRational x)
+
+      , testProperty "qiNegate" $ \n ->
+          approxEq' (qiToFloat (qiNegate n)) (negate (qiToFloat n))
+
+      , testProperty "qiRecip" $ \n ->
+          not (approxEq (qiToFloat n) 0)
+            ==> let ~(Just nr) = qiRecip n
+                in  approxEq' (qiToFloat nr) (recip (qiToFloat n))
+
+      , testProperty "qiAdd" . withCompatiblePair $ \n n' ->
+          let ~(Just r) = qiAdd n n'
+          in  approxEq' (qiToFloat r) (qiToFloat n + qiToFloat n')
+
+      , testProperty "qiSub" . withCompatiblePair $ \n n' ->
+          let ~(Just r) = qiSub n n'
+          in  approxEq' (qiToFloat r) (qiToFloat n - qiToFloat n')
+
+      , testProperty "qiMul" . withCompatiblePair $ \n n' ->
+          let ~(Just r) = qiMul n n'
+          in  approxEq' (qiToFloat r) (qiToFloat n * qiToFloat n')
+
+      , testProperty "qiDiv" . withCompatiblePair $ \n n' ->
+          let ~(Just r) = qiDiv n n'
+          in  not (approxEq (qiToFloat n') 0)
+                ==> approxEq' (qiToFloat r) (qiToFloat n / qiToFloat n')
+
+      , testProperty "qiPow" $ \n (NonNegative p) ->
+          -- Limit the power for speed.
+          (p <= 10) ==>
+            approxEq' (qiToFloat (qiPow n p))
+                      -- CReal seems to diverge in 0 ** 1, use (^).
+                      (qiToFloat n ^ p)
+
+      , testProperty "qiFloor" $ \n ->
+          qiFloor n === floor (qiToFloat n :: RefFloat)
+
+      , testProperty "qiToContinuedFraction/continuedFractionToQI" $ \n ->
+          let cf  = qiToContinuedFraction n
+              len = case cf of
+                      (_, NonCyc _)   -> 0
+                      (_, Cyc _ _ xs) -> length xs
+          -- Limit the length of the periodic part for speed.
+          in (len <= 100) ==>
+               approxEq' (qiToFloat n) (qiToFloat (continuedFractionToQI cf))
+      ]
+    ]
+
+withCompatiblePair :: Testable p
+                   => (QI -> QI -> p) -> QI -> QI -> Bool -> Bool -> Property
+withCompatiblePair f n0_ n1_ c0Zero c1Zero =
+  counterexample ("n0 = " ++ show n0) . counterexample ("n1 = " ++ show n1) $
+    f n0 n1
+  where
+    n0 = runQI n0_ $ \a b c d ->
+           qi a b (if c0Zero then 0 else c) d
+
+    n1 = runQI n0_ $ \_ _ c _ -> runQI n1_ $ \a b _ d ->
+           qi a b (if c1Zero then 0 else c) d
+
+approxQI :: Integer -> Integer -> Integer -> Integer -> RefFloat
+approxQI a b c d =
+  (fromInteger a + fromInteger b * sqrt (fromInteger c)) / fromInteger d
+
+approxQI' :: Rational -> Rational -> Integer -> RefFloat
+approxQI' a b c =
+  fromRational a + fromRational b * sqrt (fromInteger c)
+
+approxEq :: RefFloat -> RefFloat -> Bool
+approxEq a b = abs (b - a) < 1e-6 * maximum [ 1, abs a, abs b ]
+
+approxEq' :: RefFloat -> RefFloat -> Property
+approxEq' a b =
+  counterexample (show a ++ " is not approximately " ++ show b ++ " (diff = "
+                   ++ show (abs (b - a)) ++ ")")
+                 (approxEq a b)