diff --git a/.travis.yml b/.travis.yml
--- a/.travis.yml
+++ b/.travis.yml
@@ -1,114 +1,134 @@
+# This Travis job script has been generated by a script via
+#
+#   haskell-ci '--config' 'cabal.haskell-ci' '--output' '.travis.yml' 'quadratic-irrational.cabal'
+#
+# For more information, see https://github.com/haskell-CI/haskell-ci
+#
+# version: 0.2.1
+#
 language: c
-sudo: false
+dist: xenial
 
+git:
+  submodules: false  # whether to recursively clone submodules
+
 cache:
   directories:
     - $HOME/.cabal/packages
     - $HOME/.cabal/store
 
 before_cache:
-  - rm -fv $HOME/.cabal/packages/hackage.haskell.org/build-reports.log
+  - rm -fv $CABALHOME/packages/hackage.haskell.org/build-reports.log
   # remove files that are regenerated by 'cabal update'
-  - rm -fv $HOME/.cabal/packages/hackage.haskell.org/00-index.*
-  - rm -fv $HOME/.cabal/packages/hackage.haskell.org/*.json
-  - rm -fv $HOME/.cabal/packages/hackage.haskell.org/01-index.cache
-  - rm -fv $HOME/.cabal/packages/hackage.haskell.org/01-index.tar
-  - rm -fv $HOME/.cabal/packages/hackage.haskell.org/01-index.tar.idx
+  - rm -fv $CABALHOME/packages/hackage.haskell.org/00-index.*
+  - rm -fv $CABALHOME/packages/hackage.haskell.org/*.json
+  - rm -fv $CABALHOME/packages/hackage.haskell.org/01-index.cache
+  - rm -fv $CABALHOME/packages/hackage.haskell.org/01-index.tar
+  - rm -fv $CABALHOME/packages/hackage.haskell.org/01-index.tar.idx
 
-  - rm -rfv $HOME/.cabal/packages/head.hackage
+  - rm -rfv $CABALHOME/packages/head.hackage
 
 matrix:
   include:
+    - compiler: "ghc-8.8.1"
+      env: GHCHEAD=true
+      addons: {apt: {packages: [ghc-ppa-tools,cabal-install-2.4,ghc-8.8.1], sources: [hvr-ghc]}}
+    - compiler: "ghc-8.6.5"
+      addons: {apt: {packages: [ghc-ppa-tools,cabal-install-2.4,ghc-8.6.5], sources: [hvr-ghc]}}
+    - compiler: "ghc-8.4.4"
+      addons: {apt: {packages: [ghc-ppa-tools,cabal-install-2.4,ghc-8.4.4], sources: [hvr-ghc]}}
+    - compiler: "ghc-8.2.2"
+      addons: {apt: {packages: [ghc-ppa-tools,cabal-install-2.4,ghc-8.2.2], sources: [hvr-ghc]}}
+    - compiler: "ghc-8.0.2"
+      addons: {apt: {packages: [ghc-ppa-tools,cabal-install-2.4,ghc-8.0.2], sources: [hvr-ghc]}}
     - compiler: "ghc-7.10.3"
       env: TEST=--disable-tests
       addons: {apt: {packages: [ghc-ppa-tools,cabal-install-head,ghc-7.10.3], sources: [hvr-ghc]}}
-    - compiler: "ghc-8.0.2"
-      addons: {apt: {packages: [ghc-ppa-tools,cabal-install-head,ghc-8.0.2], sources: [hvr-ghc]}}
-    - compiler: "ghc-8.2.2"
-      addons: {apt: {packages: [ghc-ppa-tools,cabal-install-head,ghc-8.2.2], sources: [hvr-ghc]}}
-    - compiler: "ghc-8.4.3"
-      addons: {apt: {packages: [ghc-ppa-tools,cabal-install-head,ghc-8.4.3], sources: [hvr-ghc]}}
-    - compiler: "ghc-8.6.1"
-      env: GHCHEAD=true
-      addons: {apt: {packages: [ghc-ppa-tools,cabal-install-head,ghc-8.6.1], sources: [hvr-ghc]}}
     - compiler: "ghc-head"
       env: GHCHEAD=true
       addons: {apt: {packages: [ghc-ppa-tools,cabal-install-head,ghc-head], sources: [hvr-ghc]}}
 
   allow_failures:
+    - compiler: "ghc-8.8.1"
     - compiler: "ghc-head"
 
 before_install:
-  - HC=${CC}
+  - HC=/opt/ghc/bin/${CC}
   - HCPKG=${HC/ghc/ghc-pkg}
   - unset CC
+  - CABAL=/opt/ghc/bin/cabal
+  - CABALHOME=$HOME/.cabal
+  - export PATH="$CABALHOME/bin:$PATH"
   - ROOTDIR=$(pwd)
-  - mkdir -p $HOME/.local/bin
-  - "PATH=/opt/ghc/bin:/opt/ghc-ppa-tools/bin:$HOME/local/bin:$PATH"
   - HCNUMVER=$(( $(${HC} --numeric-version|sed -E 's/([0-9]+)\.([0-9]+)\.([0-9]+).*/\1 * 10000 + \2 * 100 + \3/') ))
   - echo $HCNUMVER
 
 install:
-  - cabal --version
+  - ${CABAL} --version
   - echo "$(${HC} --version) [$(${HC} --print-project-git-commit-id 2> /dev/null || echo '?')]"
-  - BENCH=${BENCH---enable-benchmarks}
   - TEST=${TEST---enable-tests}
-  - HADDOCK=${HADDOCK-true}
-  - INSTALLED=${INSTALLED-true}
+  - BENCH=--enable-benchmarks
   - GHCHEAD=${GHCHEAD-false}
-  - travis_retry cabal update -v
-  - "sed -i.bak 's/^jobs:/-- jobs:/' ${HOME}/.cabal/config"
+  - travis_retry ${CABAL} update -v
+  - sed -i.bak 's/^jobs:/-- jobs:/' $CABALHOME/config
   - rm -fv cabal.project cabal.project.local
-  # Overlay Hackage Package Index for GHC HEAD: https://github.com/hvr/head.hackage
   - |
     if $GHCHEAD; then
-      sed -i.bak 's/-- allow-newer:.*/allow-newer: *:base, *:template-haskell, *:ghc, *:Cabal, *:stm/' ${HOME}/.cabal/config
-
-      echo 'repository head.hackage'                                                        >> ${HOME}/.cabal/config
-      echo '   url: http://head.hackage.haskell.org/'                                       >> ${HOME}/.cabal/config
-      echo '   secure: True'                                                                >> ${HOME}/.cabal/config
-      echo '   root-keys: 07c59cb65787dedfaef5bd5f987ceb5f7e5ebf88b904bbd4c5cbdeb2ff71b740' >> ${HOME}/.cabal/config
-      echo '              2e8555dde16ebd8df076f1a8ef13b8f14c66bad8eafefd7d9e37d0ed711821fb' >> ${HOME}/.cabal/config
-      echo '              8f79fd2389ab2967354407ec852cbe73f2e8635793ac446d09461ffb99527f6e' >> ${HOME}/.cabal/config
-      echo '   key-threshold: 3'                                                            >> ${HOME}/.cabal.config
-
-      cabal new-update head.hackage -v
-    fi
-  - grep -Ev -- '^\s*--' ${HOME}/.cabal/config | grep -Ev '^\s*$'
-  - "printf 'packages: \".\"\\n' > cabal.project"
-  - cat cabal.project
-  - if [ -f "./configure.ac" ]; then
-      (cd "." && autoreconf -i);
+      sed -i.bak 's/-- allow-newer:.*/allow-newer: *:base, *:template-haskell, *:ghc, *:Cabal/' $CABALHOME/config
     fi
+  - grep -Ev -- '^\s*--' $CABALHOME/config | grep -Ev '^\s*$'
+  - rm -f cabal.project
+  - touch cabal.project
+  - "printf 'packages: \".\"\\n' >> cabal.project"
+  - "printf 'write-ghc-environment-files: always\\n' >> cabal.project"
+  - touch cabal.project.local
+  - "for pkg in $($HCPKG list --simple-output); do echo $pkg | sed 's/-[^-]*$//' | grep -vE -- '^(quadratic-irrational)$' | sed 's/^/constraints: /' | sed 's/$/ installed/' >> cabal.project.local; done"
+  - cat cabal.project || true
+  - cat cabal.project.local || true
+  - if [ -f "./configure.ac" ]; then (cd "." && autoreconf -i); fi
   - rm -f cabal.project.freeze
-  - cabal new-build -w ${HC} ${TEST} ${BENCH} --project-file="cabal.project" --dep -j2 all
-  - cabal new-build -w ${HC} --disable-tests --disable-benchmarks --project-file="cabal.project" --dep -j2 all
-  - rm -rf "."/.ghc.environment.* "."/dist
+  - ${CABAL} new-freeze -w ${HC} ${TEST} ${BENCH} --project-file="cabal.project" --dry
+  - "cat \"cabal.project.freeze\" | sed -E 's/^(constraints: *| *)//' | sed 's/any.//'"
+  - rm  "cabal.project.freeze"
+  - ${CABAL} new-build -w ${HC} ${TEST} ${BENCH} --project-file="cabal.project" --dep -j2 all
+  - ${CABAL} new-build -w ${HC} --disable-tests --disable-benchmarks --project-file="cabal.project" --dep -j2 all
+  - rm -rf .ghc.environment.* "."/dist
   - DISTDIR=$(mktemp -d /tmp/dist-test.XXXX)
 
 # Here starts the actual work to be performed for the package under test;
 # any command which exits with a non-zero exit code causes the build to fail.
 script:
   # test that source-distributions can be generated
-  - (cd "." && cabal sdist)
-  - mv "."/dist/quadratic-irrational-*.tar.gz ${DISTDIR}/
+  - ${CABAL} new-sdist all
+  - mv dist-newstyle/sdist/*.tar.gz ${DISTDIR}/
   - cd ${DISTDIR} || false
   - find . -maxdepth 1 -name '*.tar.gz' -exec tar -xvf '{}' \;
-  - "printf 'packages: quadratic-irrational-*/*.cabal\\n' > cabal.project"
-  - cat cabal.project
+  - rm -f cabal.project
+  - touch cabal.project
+  - "printf 'packages: \"quadratic-irrational-*/*.cabal\"\\n' >> cabal.project"
+  - "printf 'write-ghc-environment-files: always\\n' >> cabal.project"
+  - touch cabal.project.local
+  - "for pkg in $($HCPKG list --simple-output); do echo $pkg | sed 's/-[^-]*$//' | grep -vE -- '^(quadratic-irrational)$' | sed 's/^/constraints: /' | sed 's/$/ installed/' >> cabal.project.local; done"
+  - cat cabal.project || true
+  - cat cabal.project.local || true
   # this builds all libraries and executables (without tests/benchmarks)
-  - cabal new-build -w ${HC} --disable-tests --disable-benchmarks -f check-bounds all
-
-  # Build with installed constraints for packages in global-db
-  - if $INSTALLED; then echo cabal new-build -w ${HC} --disable-tests --disable-benchmarks $(${HCPKG} list --global --simple-output --names-only | sed 's/\([a-zA-Z0-9-]\{1,\}\) */--constraint="\1 installed" /g') all | sh; else echo "Not building with installed constraints"; fi
+  - ${CABAL} new-build -w ${HC} --disable-tests --disable-benchmarks all
 
   # build & run tests, build benchmarks
-  - cabal new-build -w ${HC} ${TEST} ${BENCH} all
-  - if [ "x$TEST" = "x--enable-tests" ]; then cabal new-test -w ${HC} ${TEST} ${BENCH} -f check-bounds all; fi
+  - ${CABAL} new-build -w ${HC} ${TEST} ${BENCH} all
+  - if [ "x$TEST" = "x--enable-tests" ]; then ${CABAL} new-test -w ${HC} ${TEST} ${BENCH} all; fi
 
+  # cabal check
+  - (cd quadratic-irrational-* && ${CABAL} check)
+
   # haddock
-  - rm -rf ./dist-newstyle
-  - if $HADDOCK; then cabal new-haddock -w ${HC} ${TEST} ${BENCH} all; else echo "Skipping haddock generation";fi
+  - ${CABAL} new-haddock -w ${HC} ${TEST} ${BENCH} all
 
-# REGENDATA ["quadratic-irrational.cabal"]
+  # Build without installed constraints for packages in global-db
+  - rm -f cabal.project.local; ${CABAL} new-build -w ${HC} --disable-tests --disable-benchmarks all;
+
+  # Constraint sets
+  - rm -rf cabal.project.local
+
+# REGENDATA ["--config","cabal.haskell-ci","--output",".travis.yml","quadratic-irrational.cabal"]
 # EOF
diff --git a/ChangeLog.md b/ChangeLog.md
--- a/ChangeLog.md
+++ b/ChangeLog.md
@@ -1,3 +1,10 @@
+# 0.1.0 (2019-04-26)
+
+* Allow imaginary square roots, e. g., `qi 1 1 (-5) 1`.
+* Remove `Ord QI` instance: complex values cannot be ordered.
+* Roots of 0 are reduced: `qi a b 0 d` becomes `qi a 0 2 d`.
+* Remove `qiZero` and `qiOne`.
+
 # 0.0.6 (2018-08-29)
 
 * Support GHC up to 8.6.1.
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.6
+version: 0.1.0
 license: MIT
 license-file: LICENSE
 author: Johan Kiviniemi <devel@johan.kiviniemi.name>
@@ -51,7 +51,7 @@
   > 3.5      = [3; 2]
   > (1+√5)/2 = [1; 1, 1, 1, …]
   > √2       = [1; 2, 2, 2, …]
-tested-with: GHC==7.10.3, GHC==8.0.2, GHC==8.2.2, GHC==8.4.3, GHC==8.6.1
+tested-with: GHC==7.10.3, GHC==8.0.2, GHC==8.2.2, GHC==8.4.3, GHC==8.6.5, GHC==8.8.1
 
 build-type: Simple
 cabal-version: >= 1.10
@@ -69,12 +69,16 @@
   exposed-modules: Numeric.QuadraticIrrational
                  , Numeric.QuadraticIrrational.CyclicList
                  , Numeric.QuadraticIrrational.Internal.Lens
+  other-modules:   Numeric.QuadraticIrrational.SquareFree
+                 , Numeric.QuadraticIrrational.Ring
   hs-source-dirs: src
   build-depends: base >= 4.8 && < 5
-               , arithmoi >= 0.4
+               , arithmoi >= 0.7
                , containers >= 0.5 && < 0.7
-               , mtl >= 2.1 && < 2.3
                , transformers >= 0.3 && < 0.6
+  if impl(ghc <8.0)
+    build-depends:
+      semigroups >=0.8
   default-language: Haskell2010
   ghc-options: -Wall -O2 -funbox-strict-fields
 
@@ -88,7 +92,7 @@
                , quadratic-irrational
                , numbers >= 3000.0 && < 3000.3
                , QuickCheck >= 2.7 && < 3
-               , tasty >= 0.8 && < 1.2
+               , tasty >= 0.8
                , tasty-quickcheck >= 0.8 && < 0.11
   default-language: Haskell2010
   ghc-options: -Wall -O2 -funbox-strict-fields
@@ -103,5 +107,6 @@
                , filepath
                , mtl
                , numbers
+               , arithmoi
   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
@@ -1,5 +1,8 @@
-{-# LANGUAGE FlexibleContexts #-}
-{-# LANGUAGE ViewPatterns     #-}
+{-# LANGUAGE FlexibleContexts    #-}
+{-# LANGUAGE GADTs               #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE StandaloneDeriving  #-}
+{-# LANGUAGE ViewPatterns        #-}
 
 -- |
 -- Module      : Numeric.QuadraticIrrational
@@ -57,7 +60,7 @@
   , -- * Lenses
     _qi, _qi', _qiABD, _qiA, _qiB, _qiC, _qiD
   , -- * Numerical operations
-    qiZero, qiOne, qiIsZero
+    qiIsZero
   , qiToFloat
   , qiAddI, qiSubI, qiMulI, qiDivI
   , qiAddR, qiSubR, qiMulR, qiDivR
@@ -69,12 +72,15 @@
   , module Numeric.QuadraticIrrational.CyclicList
   ) where
 
-import Control.Monad.State (evalState, gets, modify)
+import Control.Monad.Trans.State (evalState, gets, modify)
 import Data.Foldable (toList)
 import Data.List (foldl')
 import Data.Maybe (fromMaybe)
+import Data.Proxy
 import Data.Ratio ((%), denominator, numerator)
+import Data.Semigroup
 import qualified Data.Set as Set (empty, insert, member)
+import Data.Type.Equality
 import Math.NumberTheory.Powers.Squares (integerSquareRoot)
 import Math.NumberTheory.Primes.Factorisation (factorise)
 import Text.Read (Lexeme (Ident), Read (readListPrec, readPrec),
@@ -82,19 +88,23 @@
 
 import Numeric.QuadraticIrrational.CyclicList
 import Numeric.QuadraticIrrational.Internal.Lens
+import Numeric.QuadraticIrrational.Ring
+import Numeric.QuadraticIrrational.SquareFree
 
 -- $setup
 -- >>> import Data.Number.CReal
 
--- | @(a + b √c) \/ d@
-data QI = QI !Integer
-             !Integer
-             !Integer
-             !Integer
-  deriving (Eq)
+-- | @(a + b √c)@
+data QI where
+  QI :: KnownSquareFree c => QuadExt c Rational -> QI
 
+instance Eq QI where
+  (QI (q1 :: QuadExt c1 Rational)) == (QI (q2 :: QuadExt c2 Rational)) = case sameSquareFree (Proxy :: Proxy c1) (Proxy :: Proxy c2) of
+      Nothing   -> False
+      Just Refl -> q1 == q2
+
 instance Show QI where
-  showsPrec p (QI a b c d) = showParen (p > 10)
+  showsPrec p q = runQI q $ \a b c d -> showParen (p > 10)
                            $ showString "qi " . showsPrec 11 a
                            . showChar   ' '   . showsPrec 11 b
                            . showChar   ' '   . showsPrec 11 c
@@ -110,61 +120,6 @@
 
   readListPrec = readListPrecDefault
 
-instance Ord QI where
-  compare (QI a b c d) (QI a' b' c' d') = res
-    where
-      -- (a + b √c)/d   ⋛ (a' + b' √c')/d'
-      -- (a + b √c) d'  ⋛ (a' + b' √c') d
-      -- a d' + b d' √c ⋛ a' d + b' d √c'
-      -- a d' − a' d    ⋛ b' d √c' − b d' √c
-      --
-      -- let i = a d' − a' d
-      --     j = b' d √c'
-      --     k = b d' √c
-      i = a * d' - a' * d
-      sqJ = sq b' * sq d  * c'
-      sqK = sq b  * sq d' * c
-
-      -- i ⋛ j − k
-      --
-      -- sign (b' d √c') = sign b' because d  ≥ 0 and c' ≥ 0
-      -- sign (b d' √c)  = sign b  because d' ≥ 0 and c  ≥ 0
-      --
-      -- if j − k < 0 then (sign b') j² − (sign b) k² < 0
-      --
-      -- (sign i) |i| ⋛ sign ((sign b') j² − (sign b) k²) |j − k|
-      --
-      -- let snL = sign i
-      --     snR = sign ((sign b') j² − (sign b) k²)
-      snL = signum i
-      snR = signum (signum b' * sqJ - signum b * sqK)
-
-      -- snL |i|                ⋛ snR |j − k|
-      -- snL i²                 ⋛ snR (j − k)²
-      -- snL i²                 ⋛ snR (j² + k² − 2 j k)
-      -- snL i² − snR (j² + k²) ⋛ snR (−2) j k
-      -- snL i² − snR (j² + k²) ⋛ snR (−2) b b' d d' √c √c'
-      --
-      -- let q = snL i² − snR (j² + k²)
-      --     r = snR (−2) b b' d d' √c √c'
-      q = snL * sq i - snR * (sqJ + sqK)
-      sqR = 4 * sq b * sq b' * sq d * sq d' * c * c'
-
-      -- q ⋛ r
-      --
-      -- sign (snR (−2) b b' d d' √c √c') = sign (snR (−2) b b')
-      --
-      -- let snL' = sign q
-      --     snR' = sign (snR (−2) b b')
-      snL' = signum q
-      snR' = signum (snR * (-2) * b * b')
-
-      -- snL' |q| ⋛ snR' |r|
-      -- snL' q²  ⋛ snR' r²
-      res = compare (snL' * sq q) (snR' * sqR)
-
-      sq x = x*x
-
 type QITuple = (Integer, Integer, Integer, Integer)
 
 -- | Given @a@, @b@, @c@ and @d@ such that @n = (a + b √c)\/d@, constuct a 'QI'
@@ -173,18 +128,23 @@
 -- >>> qi 3 4 5 6
 -- qi 3 4 5 6
 --
+-- >>> qi 3 0 42 1
+-- qi 3 0 42 1
+--
+-- >>> qi 1 1 (-1) 1
+-- qi 1 1 (-1) 1
+--
 -- 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 (-30) (-40) (-5) (-60)
+-- qi 3 4 (-5) 6
 --
--- >>> qi 3 0 42 1
--- qi 3 0 0 1
+-- If @c = 0@ then @b@ is zeroed and @c@ is put equal to 2:
 --
 -- >>> qi 3 42 0 1
--- qi 3 0 0 1
+-- qi 3 0 2 1
 --
 -- The @b √c@ term is simplified:
 --
@@ -194,33 +154,30 @@
 -- If @c = 1@ (after simplification) then @b@ is moved to @a@:
 --
 -- >>> qi 1 5 (2*2) 1
--- qi 11 0 0 1
+-- qi 11 0 2 1
 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
+qi a _ 0 (nonZero "qi" -> d) = simplifyReduceCons a 0 2 d
+qi a b c (nonZero "qi" -> d) = 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
+qiNoSimpl a b (nonZero "qiNoSimpl" -> c) (nonZero "qiNoSimpl" -> d)
+  | b == 0    = reduceCons a 0 2 d
+  | c == 0    = reduceCons a 0 2 d
+  | c == 1    = reduceCons (a + b) 0 2 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
-  | c' == 1   = reduceCons (a + b') 0 0 d
+  | c' == 1   = reduceCons (a + b') 0 2 d
   | otherwise = reduceCons a b' c' d
   where ~(b', c') = separateSquareFactors b c
 {-# INLINE simplifyReduceCons #-}
@@ -228,11 +185,11 @@
 -- | 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')
+separateSquareFactors b (nonZero "separateSquareFactors" -> c) =
+  case foldl' go (1,1) (factorise (abs c)) of
+    ~(bMul, c') -> (b*bMul, c' * signum c)
   where
-    go :: (Integer, Integer) -> (Integer, Int) -> (Integer, Integer)
+    go :: Integral a => (Integer, Integer) -> (Integer, a) -> (Integer, Integer)
     go ~(i, j) ~(fac, pow) =
       i `seq` j `seq` fac `seq` pow `seq`
         if even pow
@@ -241,9 +198,9 @@
 
 -- 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)
+reduceCons a b c (nonZero "reduceCons" -> d) = case someSquareFreeVal c of
+  Just (SomeSquareFree (Proxy :: Proxy c)) -> QI (QuadExt (a % d) (b % d) :: QuadExt c Rational)
+  Nothing -> error $ "reduceCons: square root term is not square free and equals to " ++ show c
 {-# INLINE reduceCons #-}
 
 -- | Given @a@, @b@ and @c@ such that @n = a + b √c@, constuct a 'QI'
@@ -251,11 +208,13 @@
 --
 -- >>> qi' 0.5 0.7 2
 -- qi 5 7 2 10
+-- >>> qi' 0.3 0.4 (-3)
+-- qi 3 4 (-3) 10
 qi' :: Rational  -- ^ a
     -> Rational  -- ^ b
     -> Integer   -- ^ c
     -> QI
-qi' a b (nonNegative "qi'" -> c) = n
+qi' a b c = n
   where
     -- (aN/aD) + (bN/bD) √c = ((aN bD) + (bN aD) √c) / (aD bD)
     n = qi (aN*bD) (bN*aD) c (aD*bD)
@@ -267,22 +226,35 @@
 --
 -- >>> runQI (qi 3 4 5 6) (\a b c d -> (a,b,c,d))
 -- (3,4,5,6)
+-- >>> runQI (qi 1 1 (-1) 1) (\a b c d -> (a,b,c,d))
+-- (1,1,-1,1)
 runQI :: QI -> (Integer -> Integer -> Integer -> Integer -> a) -> a
-runQI (QI a b c d) f = f a b c d
+runQI (QI (QuadExt a b :: QuadExt c Rational)) f = f a' b' c d
+  where
+    c = squareFreeVal (Proxy :: Proxy c)
+    d = lcm (denominator a) (denominator b)
+    a' = numerator a * (d `quot` denominator a)
+    b' = numerator b * (d `quot` denominator b)
 {-# 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' 0.3 0.4 (-3)) (\a b c -> (a, b, c))
+-- (3 % 10,2 % 5,-3)
 runQI' :: QI -> (Rational -> Rational -> Integer -> a) -> a
-runQI' (QI a b c d) f = f (a % d) (b % d) c
+runQI' (QI (QuadExt a b :: QuadExt c Rational)) f = f a b c
+  where
+    c = squareFreeVal (Proxy :: Proxy 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 1 1 (-1) 1)
+-- (1,1,-1,1)
 unQI :: QI -> (Integer, Integer, Integer, Integer)
 unQI n = runQI n (,,,)
 {-# INLINE unQI #-}
@@ -291,6 +263,8 @@
 --
 -- >>> unQI' (qi' 0.5 0.7 2)
 -- (1 % 2,7 % 10,2)
+-- >>> unQI' (qi' 0.3 0.4 (-3))
+-- (3 % 10,2 % 5,-3)
 unQI' :: QI -> (Rational, Rational, Integer)
 unQI' n = runQI' n (,,)
 {-# INLINE unQI' #-}
@@ -299,9 +273,13 @@
 --
 -- >>> view _qi (qi 3 4 5 6)
 -- (3,4,5,6)
+-- >>> view _qi (qi 1 1 (-1) 1)
+-- (1,1,-1,1)
 --
 -- >>> over _qi (\(a,b,c,d) -> (a+10, b+10, c+10, d+10)) (qi 3 4 5 6)
 -- qi 13 14 15 16
+-- >>> over _qi (\(a,b,c,d) -> (a+10, b+10, c+10, d+10)) (qi 1 1 (-1) 1)
+-- qi 4 0 2 1
 _qi :: Lens' QI (Integer, Integer, Integer, Integer)
 _qi f n = (\ ~(a',b',c',d') -> qi a' b' c' d') <$> f (unQI n)
 {-# INLINE _qi #-}
@@ -310,9 +288,13 @@
 --
 -- >>> view _qi' (qi' 0.5 0.7 2)
 -- (1 % 2,7 % 10,2)
+-- >>> view _qi' (qi' 0.3 0.4 (-3))
+-- (3 % 10,2 % 5,-3)
 --
 -- >>> over _qi' (\(a,b,c) -> (a/5, b/6, c*3)) (qi 3 4 5 6)
 -- qi 9 10 15 90
+-- >>> over _qi' (\(a,b,c) -> (a/5, b/6, c*3)) (qi 1 1 (-1) 1)
+-- qi 6 5 (-3) 30
 _qi' :: Lens' QI (Rational, Rational, Integer)
 _qi' f n = (\ ~(a',b',c') -> qi' a' b' c') <$> f (unQI' n)
 {-# INLINE _qi' #-}
@@ -322,9 +304,13 @@
 --
 -- >>> view _qiABD (qi 3 4 5 6)
 -- (3,4,6)
+-- >>> view _qiABD (qi 1 1 (-1) 1)
+-- (1,1,1)
 --
 -- >>> over _qiABD (\(a,b,d) -> (a+10, b+10, d+10)) (qi 3 4 5 6)
 -- qi 13 14 5 16
+-- >>> over _qiABD (\(a,b,d) -> (a+10, b+10, d+10)) (qi 1 1 (-1) 1)
+-- qi 1 1 (-1) 1
 _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)
@@ -335,9 +321,13 @@
 --
 -- >>> view _qiA (qi 3 4 5 6)
 -- 3
+-- >>> view _qiA (qi 1 1 (-1) 1)
+-- 1
 --
 -- >>> over _qiA (+ 10) (qi 3 4 5 6)
 -- qi 13 4 5 6
+-- >>> over _qiA (+ 10) (qi 1 1 (-1) 1)
+-- qi 11 1 (-1) 1
 _qiA :: Lens' QI Integer
 _qiA = _qiABD . go
   where go f ~(a,b,d) = (\a' -> (a',b,d)) <$> f a
@@ -347,9 +337,13 @@
 --
 -- >>> view _qiB (qi 3 4 5 6)
 -- 4
+-- >>> view _qiB (qi 1 1 (-1) 1)
+-- 1
 --
 -- >>> over _qiB (+ 10) (qi 3 4 5 6)
 -- qi 3 14 5 6
+-- >>> over _qiB (+ 10) (qi 1 1 (-1) 1)
+-- qi 1 11 (-1) 1
 _qiB :: Lens' QI Integer
 _qiB = _qiABD . go
   where go f ~(a,b,d) = (\b' -> (a,b',d)) <$> f b
@@ -359,9 +353,13 @@
 --
 -- >>> view _qiC (qi 3 4 5 6)
 -- 5
+-- >>> view _qiC (qi 1 1 (-1) 1)
+-- -1
 --
 -- >>> over _qiC (+ 10) (qi 3 4 5 6)
 -- qi 3 4 15 6
+-- >>> over _qiC (+ 10) (qi 1 1 (-1) 1)
+-- qi 4 0 2 1
 _qiC :: Lens' QI Integer
 _qiC = _qi . go
   where go f ~(a,b,c,d) = (\c' -> (a,b,c',d)) <$> f c
@@ -371,42 +369,32 @@
 --
 -- >>> view _qiD (qi 3 4 5 6)
 -- 6
+-- >>> view _qiD (qi 1 1 (-1) 1)
+-- 1
 --
 -- >>> over _qiD (+ 10) (qi 3 4 5 6)
 -- qi 3 4 5 16
+-- >>> over _qiD (+ 10) (qi 1 1 (-1) 1)
+-- qi 1 1 (-1) 11
 _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.
---
--- >>> 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]
+-- >>> map qiIsZero [qi 0 0 0 1, qi 1 0 0 1, qiSubR (qi 7 0 0 2) 3.5, qi 0 9 0 1, qi 0 0 (-1) 5]
+-- [True,False,True,True,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
+qiIsZero (QI q) = q == 0
 {-# INLINE qiIsZero #-}
 
 -- | Convert a 'QI' number into a 'Floating' one.
+-- Returns NaN for imaginary values.
 --
 -- >>> qiToFloat (qi 3 4 5 6) == ((3 + 4 * sqrt 5)/6 :: Double)
 -- True
+-- >>> qiToFloat (qi 2 3 (-5) 7)
+-- NaN
 qiToFloat :: Floating a => QI -> a
 qiToFloat (unQI -> ~(a,b,c,d)) =
   (fromInteger a + fromInteger b * sqrt (fromInteger c)) / fromInteger d
@@ -416,219 +404,213 @@
 --
 -- >>> qi 3 4 5 6 `qiAddI` 1
 -- qi 9 4 5 6
+-- >>> qi 1 1 (-1) 1 `qiAddI` 1
+-- qi 2 1 (-1) 1
 qiAddI :: QI -> Integer -> QI
-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)
+qiAddI (QI q) x = QI (q + fromInteger x)
 {-# INLINE qiAddI #-}
 
 -- | Add a 'Rational' to a 'QI'.
 --
 -- >>> qi 3 4 5 6 `qiAddR` 1.2
 -- qi 51 20 5 30
+-- >>> qi 1 1 (-1) 1 `qiAddR` 1.2
+-- qi 11 5 (-1) 5
 qiAddR :: QI -> Rational -> QI
-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)
+qiAddR (QI q) x = QI (q + fromRational x)
 {-# INLINE qiAddR #-}
 
 -- | Subtract an 'Integer' from a 'QI'.
 --
 -- >>> qi 3 4 5 6 `qiSubI` 1
 -- qi (-3) 4 5 6
+-- >>> qi 1 1 (-1) 1 `qiSubI` 1
+-- qi 0 1 (-1) 1
 qiSubI :: QI -> Integer -> QI
-qiSubI n x = qiAddI n (negate x)
+qiSubI (QI q) x = QI (q - fromInteger x)
 {-# INLINE qiSubI #-}
 
 -- | Subtract a 'Rational' from a 'QI'.
 --
 -- >>> qi 3 4 5 6 `qiSubR` 1.2
 -- qi (-21) 20 5 30
+-- >>> qi 1 1 (-1) 1 `qiSubR` 1.2
+-- qi (-1) 5 (-1) 5
 qiSubR :: QI -> Rational -> QI
-qiSubR n x = qiAddR n (negate x)
+qiSubR (QI q) x = QI (q - fromRational x)
 {-# INLINE qiSubR #-}
 
 -- | Multiply a 'QI' by an 'Integer'.
 --
 -- >>> qi 3 4 5 6 `qiMulI` 2
 -- qi 3 4 5 3
+-- >>> qi 1 1 (-1) 1 `qiMulI` 2
+-- qi 2 2 (-1) 1
 qiMulI :: QI -> Integer -> QI
-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)
+qiMulI (QI q) x = QI (q * fromInteger x)
 {-# INLINE qiMulI #-}
 
 -- | Multiply a 'QI' by a 'Rational'.
 --
 -- >>> qi 3 4 5 6 `qiMulR` 0.5
 -- qi 3 4 5 12
+-- >>> qi 1 1 (-1) 1 `qiMulR` 0.5
+-- qi 1 1 (-1) 2
 qiMulR :: QI -> Rational -> QI
-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)
+qiMulR (QI q) x = QI (q * fromRational x)
 {-# INLINE qiMulR #-}
 
--- | Divice a 'QI' by an 'Integer'.
+-- | Divide a 'QI' by an 'Integer'.
 --
 -- >>> qi 3 4 5 6 `qiDivI` 2
 -- qi 3 4 5 12
+-- >>> qi 1 1 (-1) 1 `qiDivI` 2
+-- qi 1 1 (-1) 2
 qiDivI :: QI -> Integer -> QI
-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)
+qiDivI (QI q) (nonZero "qiDivI" -> x) = QI (q * fromRational (1 % x))
 {-# INLINE qiDivI #-}
 
--- | Divice a 'QI' by a 'Rational'.
+-- | Divide a 'QI' by a 'Rational'.
 --
 -- >>> qi 3 4 5 6 `qiDivR` 0.5
 -- qi 3 4 5 3
+-- >>> qi 1 1 (-1) 1 `qiDivR` 0.5
+-- qi 2 2 (-1) 1
 qiDivR :: QI -> Rational -> QI
-qiDivR n (nonZero "qiDivR" -> x) = qiMulR n (recip x)
+qiDivR (QI q) (nonZero "qiDivR" -> x) = QI (q * fromRational (recip x))
 {-# INLINE qiDivR #-}
 
 -- | Negate a 'QI'.
 --
 -- >>> qiNegate (qi 3 4 5 6)
 -- qi (-3) (-4) 5 6
+-- >>> qiNegate (qi 1 1 (-1) 1)
+-- qi (-1) (-1) (-1) 1
 qiNegate :: QI -> QI
-qiNegate n = over _qiABD go n
-  where go ~(a,b,d) = a `seq` b `seq` d `seq` (negate a, negate b, d)
+qiNegate (QI q) = QI (negate q)
 {-# INLINE qiNegate #-}
 
 -- | Compute the reciprocal of a 'QI'.
 --
 -- >>> qiRecip (qi 5 0 0 2)
--- Just (qi 2 0 0 5)
+-- Just (qi 2 0 2 5)
 --
 -- >>> qiRecip (qi 0 1 5 2)
 -- Just (qi 0 2 5 5)
 --
--- >>> qiRecip qiZero
+-- >>> qiRecip (qi 1 1 (-1) 1)
+-- Just (qi 1 (-1) (-1) 2)
+--
+-- >>> qiRecip (qi 0 0 0 1)
 -- Nothing
 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 (set _qiABD (a * d, negate (b * d), denom) n)
-  where denom = (a*a - b*b * c)
+qiRecip (QI 0) = Nothing
+qiRecip (QI q) = Just (QI (recip q))
 
 -- | Add two 'QI's if the square root terms are the same or zeros.
 --
--- >>> qi 3 4 5 6 `qiAdd` qiOne
+-- >>> qi 3 4 5 6 `qiAdd` qi 1 0 5 1
 -- Just (qi 9 4 5 6)
 --
 -- >>> qi 3 4 5 6 `qiAdd` qi 3 4 5 6
 -- Just (qi 3 4 5 3)
 --
+-- >>> 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 (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
+qiAdd (QI (q1 :: QuadExt c1 Rational)) (QI (q2 :: QuadExt c2 Rational)) =
+  case sameSquareFree (Proxy :: Proxy c1) (Proxy :: Proxy c2) of
+    Nothing   -> Nothing
+    Just Refl -> Just (QI (q1 + q2))
 
 -- | Subtract two 'QI's if the square root terms are the same or zeros.
 --
--- >>> qi 3 4 5 6 `qiSub` qiOne
+-- >>> qi 3 4 5 6 `qiSub` qi 1 0 5 1
 -- Just (qi (-3) 4 5 6)
 --
 -- >>> qi 3 4 5 6 `qiSub` qi 3 4 5 6
--- Just (qi 0 0 0 1)
+-- Just (qi 0 0 5 1)
 --
+-- >>> qi 3 4 (-5) 6 `qiSub` qi 3 4 (-5) 6
+-- Just (qi 0 0 (-5) 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')
+qiSub (QI (q1 :: QuadExt c1 Rational)) (QI (q2 :: QuadExt c2 Rational)) =
+  case sameSquareFree (Proxy :: Proxy c1) (Proxy :: Proxy c2) of
+    Nothing   -> Nothing
+    Just Refl -> Just (QI (q1 - q2))
 
 -- | 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` qi 0 0 5 1
+-- Just (qi 0 0 5 1)
 --
--- >>> qi 3 4 5 6 `qiMul` qiOne
+-- >>> qi 3 4 5 6 `qiMul` qi 1 0 5 1
 -- Just (qi 3 4 5 6)
 --
 -- >>> qi 3 4 5 6 `qiMul` qi 3 4 5 6
 -- Just (qi 89 24 5 36)
 --
+-- >>> qi 3 4 (-5) 6 `qiMul` qi 3 4 (-5) 6
+-- Just (qi (-71) 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'
-  -- 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 (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
+qiMul (QI (q1 :: QuadExt c1 Rational)) (QI (q2 :: QuadExt c2 Rational)) =
+  case sameSquareFree (Proxy :: Proxy c1) (Proxy :: Proxy c2) of
+    Nothing   -> Nothing
+    Just Refl -> Just (QI (q1 * q2))
 
 -- | Divide two 'QI's if the square root terms are the same or zeros.
 --
--- >>> qi 3 4 5 6 `qiDiv` qiZero
+-- >>> qi 3 4 5 6 `qiDiv` qi 0 0 5 1
 -- Nothing
 --
--- >>> qi 3 4 5 6 `qiDiv` qiOne
+-- >>> qi 3 4 5 6 `qiDiv` qi 1 0 5 1
 -- Just (qi 3 4 5 6)
 --
 -- >>> qi 3 4 5 6 `qiDiv` qi 3 4 5 6
--- Just (qi 1 0 0 1)
+-- Just (qi 1 0 5 1)
 --
 -- >>> qi 3 4 5 6 `qiDiv` qi 0 1 5 1
 -- Just (qi 20 3 5 30)
 --
+-- >>> 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'
+qiDiv _ (QI 0) = Nothing
+qiDiv (QI (q1 :: QuadExt c1 Rational)) (QI (q2 :: QuadExt c2 Rational)) =
+  case sameSquareFree (Proxy :: Proxy c1) (Proxy :: Proxy c2) of
+    Nothing   -> Nothing
+    Just Refl -> Just (QI (q1 / q2))
 
 -- | Exponentiate a 'QI' to an 'Integer' power.
 --
 -- >>> qi 3 4 5 6 `qiPow` 0
--- qi 1 0 0 1
+-- qi 1 0 5 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
+--
+-- >>> qi 1 1 (-1) 1 `qiPow` 4
+-- qi (-4) 0 (-1) 1
 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
+qiPow (QI q) pow = QI (getProduct (stimes pow (Product q)))
 
 -- | Compute the floor of a 'QI'.
+-- Throws an error for imaginary values.
 --
 -- >>> qiFloor (qi 10 0 0 2)
 -- 5
@@ -638,6 +620,10 @@
 --
 -- >>> qiFloor (qi 10 2 5 2)
 -- 7
+--
+-- >>> qiFloor (qi 1 1 (-1) 1)
+-- *** Exception: ...
+-- ...
 qiFloor :: QI -> Integer
 qiFloor (unQI -> ~(a,b,c,d)) =
   -- n = (a + b √c)/d
@@ -654,7 +640,7 @@
 -- @[2; 2] = 2 + 1\/2 = 5\/2@.
 --
 -- >>> continuedFractionToQI (2,CycList [2] [])
--- qi 5 0 0 2
+-- qi 5 0 2 2
 --
 -- The golden ratio is @[1; 1, 1, …]@.
 --
@@ -666,7 +652,7 @@
 continuedFractionToQI :: (Integer, CycList Integer) -> QI
 continuedFractionToQI (i0_, is_) = qiAddI (go is_) i0_
   where
-    go (CycList as bs) = goNonCyc as (if null bs then qiZero else goCyc bs)
+    go (CycList as bs) = goNonCyc as (if null bs then qi 0 0 0 1 else goCyc bs)
 
     goNonCyc ((pos -> i):is) final = sudoQIRecip (qiAddI (goNonCyc is final) i)
     goNonCyc []              final = final
@@ -705,6 +691,7 @@
     pos = positive "continuedFractionToQI"
 
 -- | Convert a 'QI' into a (possibly periodic) simple continued fraction.
+-- Throws an error for imaginary values.
 --
 -- @5\/2 = 2 + 1\/2 = [2; 2]@.
 --
@@ -718,6 +705,10 @@
 --
 -- >>> qiToContinuedFraction (qi 987601513930253257378987883 1 14116473325908285531353005 81983584717737887813195873886)
 -- (0,CycList [83,78,65,75,69] [32,66,65,68,71,69,82])
+--
+-- >>> qiToContinuedFraction (qi 1 1 (-1) 1)
+-- *** Exception: ...
+-- ...
 qiToContinuedFraction :: QI
                       -> (Integer, CycList Integer)
 qiToContinuedFraction num
@@ -772,10 +763,6 @@
     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)
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
@@ -13,6 +13,7 @@
   ( CycList(..)
   ) where
 
+import Data.Foldable (toList)
 import Data.Monoid ((<>))
 
 -- | A container for a possibly cyclic list.
@@ -29,6 +30,10 @@
   deriving (Eq, Ord, Read, Show, Functor)
 
 instance Foldable CycList where
-  foldMap f (CycList as bs) = foldMap f as <> if null bs then mempty else cycleAppend (foldMap f bs)
+  toList (CycList as []) = as
+  toList (CycList as bs) = as <> cycle bs
+
+  foldMap f (CycList as []) = foldMap f as
+  foldMap f (CycList as bs) = foldMap f as <> cycleAppend (foldMap f bs)
     where
       cycleAppend x = x <> cycleAppend x
diff --git a/src/Numeric/QuadraticIrrational/Ring.hs b/src/Numeric/QuadraticIrrational/Ring.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/QuadraticIrrational/Ring.hs
@@ -0,0 +1,48 @@
+{-# LANGUAGE DataKinds           #-}
+{-# LANGUAGE KindSignatures      #-}
+{-# LANGUAGE PolyKinds           #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+
+module Numeric.QuadraticIrrational.Ring
+  ( QuadExt(..)
+  ) where
+
+import Data.Proxy
+import Numeric.QuadraticIrrational.SquareFree
+
+data QuadExt (c :: SquareFree) t = QuadExt t t
+  deriving (Eq)
+
+instance (KnownSquareFree c, Show t) => Show (QuadExt c t) where
+  showsPrec p (QuadExt a b) = showParen (p > 10)
+    $ showString "QuadExt " . showsPrec 11 a . showChar '+' . showsPrec 11 b . showChar '√' . showsPrec 11 c
+    where
+      c = squareFreeVal (Proxy :: Proxy c)
+
+instance (KnownSquareFree c, Num t) => Num (QuadExt c t) where
+  QuadExt a1 b1 + QuadExt a2 b2 = QuadExt (a1 + a2) (b1 + b2)
+  QuadExt a1 b1 - QuadExt a2 b2 = QuadExt (a1 - a2) (b1 - b2)
+  negate (QuadExt a b) = QuadExt (negate a) (negate b)
+  abs q = q
+  signum _ = 1
+  fromInteger n = QuadExt (fromInteger n) 0
+
+  QuadExt a1 b1 * QuadExt a2 b2 = QuadExt a b
+    where
+      a = a1 * a2 + fromInteger c * b1 * b2
+      b = a1 * b2 + a2 * b1
+      c = squareFreeVal (Proxy :: Proxy c)
+
+instance (KnownSquareFree c, Fractional t) => Fractional (QuadExt c t) where
+  fromRational r = QuadExt (fromRational r) 0
+
+  recip (QuadExt a b) = QuadExt (a / d) (negate b / d)
+    where
+      d = a * a - fromInteger c * b * b
+      c = squareFreeVal (Proxy :: Proxy c)
+
+  QuadExt a1 b1 / QuadExt a2 b2 = QuadExt (a / d) (b / d)
+    where
+      QuadExt a b = QuadExt a1 b1 * QuadExt a2 (negate b2) :: QuadExt c t
+      d = a2 * a2 - fromInteger c * b2 * b2
+      c = squareFreeVal (Proxy :: Proxy c)
diff --git a/src/Numeric/QuadraticIrrational/SquareFree.hs b/src/Numeric/QuadraticIrrational/SquareFree.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/QuadraticIrrational/SquareFree.hs
@@ -0,0 +1,48 @@
+{-# LANGUAGE DataKinds           #-}
+{-# LANGUAGE GADTs               #-}
+{-# LANGUAGE KindSignatures      #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE TypeOperators       #-}
+
+module Numeric.QuadraticIrrational.SquareFree
+  ( SquareFree()
+  , KnownSquareFree(..)
+  , SomeSquareFree(..)
+  , someSquareFreeVal
+  , sameSquareFree
+  ) where
+
+import Data.Proxy
+import Data.Type.Equality
+import GHC.TypeLits
+import Math.NumberTheory.ArithmeticFunctions
+import Unsafe.Coerce
+
+data SquareFree = Negative Nat | Positive Nat
+
+class KnownSquareFree (c :: SquareFree) where
+  squareFreeVal :: Proxy c -> Integer
+
+instance KnownNat n => KnownSquareFree ('Negative n) where
+  squareFreeVal (Proxy :: Proxy ('Negative n)) = negate $ natVal (Proxy :: Proxy n)
+
+instance KnownNat n => KnownSquareFree ('Positive n) where
+  squareFreeVal (Proxy :: Proxy ('Positive n)) = natVal (Proxy :: Proxy n)
+
+data SomeSquareFree where
+  SomeSquareFree :: KnownSquareFree c => Proxy c -> SomeSquareFree
+
+someSquareFreeVal :: Integer -> Maybe SomeSquareFree
+someSquareFreeVal 0 = Nothing
+someSquareFreeVal n = case moebius n of
+  MoebiusZ -> Nothing
+  _ -> case someNatVal n of
+    Just (SomeNat (Proxy :: Proxy n)) -> Just (SomeSquareFree (Proxy :: Proxy ('Positive n)))
+    Nothing -> case someNatVal (negate n) of
+      Just (SomeNat (Proxy :: Proxy n)) -> Just (SomeSquareFree (Proxy :: Proxy ('Negative n)))
+      Nothing -> Nothing
+
+sameSquareFree :: (KnownSquareFree a, KnownSquareFree b) => Proxy a -> Proxy b -> Maybe (a :~: b)
+sameSquareFree x y
+  | squareFreeVal x == squareFreeVal y = Just (unsafeCoerce Refl)
+  | otherwise = Nothing
diff --git a/tests/QuadraticIrrational.hs b/tests/QuadraticIrrational.hs
--- a/tests/QuadraticIrrational.hs
+++ b/tests/QuadraticIrrational.hs
@@ -5,9 +5,7 @@
 
 import Data.Number.CReal (CReal)
 import Test.Tasty (TestTree, testGroup)
-import Test.Tasty.QuickCheck (Arbitrary (arbitrary, shrink),
-  NonNegative (NonNegative), NonZero (NonZero), Property, Testable, (===),
-  (==>), conjoin, counterexample, testProperty)
+import Test.Tasty.QuickCheck
 
 import Numeric.QuadraticIrrational
 import Numeric.QuadraticIrrational.Internal.Lens
@@ -87,16 +85,9 @@
           approxEq' (qiToFloat (qi a b c d)) (approxQI a b c d)
 
       , testProperty "compare equals" $ \a ->
-          conjoin [ a === a, compare a a === EQ ]
+          conjoin [ a === a ]
             `const` (a :: QI)
 
-      , testProperty "compare" $ \a b ->
-          let a' = qiToFloat a :: RefFloat
-              b' = qiToFloat b :: RefFloat
-          in  conjoin [ (a == b)    === (a' == b')
-                      , compare a b === compare a' b'
-                      ]
-
       , testProperty "qiAddI" $ \n x ->
           approxEq' (qiToFloat (qiAddI n x)) (qiToFloat n + fromInteger x)
 
@@ -163,7 +154,7 @@
           let cf@(_, CycList _ xs)  = qiToContinuedFraction n
           -- Limit the length of the periodic part for speed.
           in (length xs <= 100) ==>
-               n === continuedFractionToQI cf
+              (qiToFloat n :: Double) === qiToFloat (continuedFractionToQI cf)
 
       , testProperty "continuedFractionApproximate" $ \n ->
           let cf = qiToContinuedFraction n
@@ -173,16 +164,16 @@
     ]
 
 withCompatiblePair :: Testable p
-                   => (QI -> QI -> p) -> QI -> QI -> Bool -> Bool -> Property
-withCompatiblePair f n0_ n1_ c0Zero c1Zero =
+                   => (QI -> QI -> p) -> QI -> QI -> Property
+withCompatiblePair f n0_ n1_ =
   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
+           qi a b c d
 
     n1 = runQI n0_ $ \_ _ c _ -> runQI n1_ $ \a b _ d ->
-           qi a b (if c1Zero then 0 else c) d
+           qi a b c d
 
 approxQI :: Integer -> Integer -> Integer -> Integer -> RefFloat
 approxQI a b c d =
diff --git a/tests/doctests.hs b/tests/doctests.hs
--- a/tests/doctests.hs
+++ b/tests/doctests.hs
@@ -2,9 +2,8 @@
 
 module Main (main) where
 
-import Control.Applicative (empty)
 import Control.Monad (guard)
-import Control.Monad.List (ListT (ListT), liftIO,  runListT)
+import Data.Foldable (toList)
 import Data.List (isSuffixOf)
 import System.Directory (doesDirectoryExist, doesFileExist,
   getDirectoryContents)
@@ -17,17 +16,22 @@
   doctest ("-isrc" : sources)
 
 findSources :: FilePath -> IO [FilePath]
-findSources dir = runListT (goDir dir)
+findSources = goDir
   where
-    goItem :: FilePath -> FilePath -> ListT IO FilePath
-    goItem _ ('.':_) = empty
+    goItem :: FilePath -> FilePath -> IO [FilePath]
+    goItem _ ('.':_) = pure []
     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
+      isDir  <- doesDirectoryExist path
+      isFile <- doesFileExist path
+      if | isDir     -> goDir path
+         | isFile    -> pure $ toList $ goFile path
+         | otherwise -> pure []
 
-    goDir  path = goItem path =<< ListT (getDirectoryContents path)
+    goDir :: FilePath -> IO [FilePath]
+    goDir path = do
+      items <- getDirectoryContents path
+      concat <$> traverse (goItem path) items
+
+    goFile :: FilePath -> Maybe FilePath
     goFile path = path <$ guard (".hs" `isSuffixOf` path)