diff --git a/arithmetic.cabal b/arithmetic.cabal
--- a/arithmetic.cabal
+++ b/arithmetic.cabal
@@ -1,5 +1,5 @@
 name: arithmetic
-version: 1.0
+version: 1.1
 category: Number Theory
 synopsis: Natural number arithmetic
 license: MIT
@@ -10,7 +10,7 @@
 maintainer: Joe Leslie-Hurd <joe@gilith.com>
 description:
   This package implements a library of natural number arithmetic functions,
-  including Montgomery multiplication.
+  including Montgomery multiplication and continued fractions.
 
 Library
   build-depends:
@@ -20,14 +20,18 @@
     opentheory-primitive >= 1.0 && < 2.0,
     opentheory >= 1.0 && < 2.0,
     opentheory-bits >= 1.0 && < 2.0,
-    opentheory-divides >= 1.0 && < 2.0
+    opentheory-divides >= 1.0 && < 2.0,
+    opentheory-prime >= 1.0 && < 2.0
   hs-source-dirs: src
   ghc-options: -Wall
   exposed-modules:
+    Arithmetic.ContinuedFraction,
     Arithmetic.Modular,
     Arithmetic.Montgomery,
     Arithmetic.Prime,
-    Arithmetic.Random
+    Arithmetic.Random,
+    Arithmetic.Smooth,
+    Arithmetic.SquareRoot
 
 executable arithmetic
   build-depends:
@@ -37,7 +41,8 @@
     opentheory-primitive >= 1.0 && < 2.0,
     opentheory >= 1.0 && < 2.0,
     opentheory-bits >= 1.0 && < 2.0,
-    opentheory-divides >= 1.0 && < 2.0
+    opentheory-divides >= 1.0 && < 2.0,
+    opentheory-prime >= 1.0 && < 2.0
   hs-source-dirs: src
   ghc-options: -Wall
   main-is: Main.hs
@@ -51,10 +56,8 @@
     opentheory-primitive >= 1.0 && < 2.0,
     opentheory >= 1.0 && < 2.0,
     opentheory-bits >= 1.0 && < 2.0,
-    opentheory-divides >= 1.0 && < 2.0
+    opentheory-divides >= 1.0 && < 2.0,
+    opentheory-prime >= 1.0 && < 2.0
   hs-source-dirs: src
   ghc-options: -Wall
   main-is: Test.hs
-  other-modules:
-    IntegerDivides,
-    NaturalDivides
diff --git a/src/Arithmetic/ContinuedFraction.hs b/src/Arithmetic/ContinuedFraction.hs
new file mode 100644
--- /dev/null
+++ b/src/Arithmetic/ContinuedFraction.hs
@@ -0,0 +1,102 @@
+{- |
+module: Arithmetic.ContinuedFraction
+description: Continued fractions
+license: MIT
+
+maintainer: Joe Leslie-Hurd <joe@gilith.com>
+stability: provisional
+portability: portable
+-}
+module Arithmetic.ContinuedFraction
+where
+
+import OpenTheory.Primitive.Natural
+
+newtype ContinuedFraction =
+    ContinuedFraction {unContinuedFraction :: (Natural,[Natural])}
+  deriving Eq
+
+fromNatural :: Natural -> ContinuedFraction
+fromNatural n = ContinuedFraction (n,[])
+
+goldenRatio :: ContinuedFraction
+goldenRatio = ContinuedFraction (1, repeat 1)
+
+naturalLogarithmBase :: ContinuedFraction
+naturalLogarithmBase =
+    ContinuedFraction (2, go 2)
+  where
+    go n = 1 : n : 1 : go (n + 2)
+
+convergentsFn :: (Natural -> a) -> (a -> a -> a) -> (a -> a -> a) ->
+                 [Natural] -> a -> a -> [a]
+convergentsFn lift add mult =
+    go
+  where
+    go [] _ _ = []
+    go (q : qs) y x =
+        z : go qs z y
+      where
+        z = add (mult (lift q) y) x
+
+numerators :: (Natural -> a) -> (a -> a -> a) -> (a -> a -> a) ->
+              ContinuedFraction -> [a]
+numerators lift add mult (ContinuedFraction (q0,qs)) =
+    x : convergentsFn lift add mult qs x one
+  where
+    x = lift q0
+    one = lift 1
+
+denominators :: (Natural -> a) -> (a -> a -> a) -> (a -> a -> a) ->
+                ContinuedFraction -> [a]
+denominators lift add mult (ContinuedFraction (_,qs)) =
+    one : convergentsFn lift add mult qs one zero
+  where
+    one = lift 1
+    zero = lift 0
+
+convergents :: (Natural -> a) -> (a -> a -> a) -> (a -> a -> a) ->
+               (a -> a -> a) -> ContinuedFraction -> [a]
+convergents lift add mult divf cf =
+    zipWith divf nums dens
+  where
+    nums = numerators lift add mult cf
+    dens = denominators lift add mult cf
+
+unstableConvergents :: Eq a => [a] -> [a]
+unstableConvergents [] = error "empty convergents"
+unstableConvergents (q0 : qs) =
+    q0 : go q0 qs
+  where
+    go _ [] = []
+    go x (h : t) = if x == h then [] else h : go h t
+
+fractionalConvergents :: Fractional a => ContinuedFraction -> [a]
+fractionalConvergents = convergents fromIntegral (+) (*) (/)
+
+rationalConvergents :: ContinuedFraction -> [Rational]
+rationalConvergents = convergents fromIntegral (+) (*) (/)
+
+toDouble :: ContinuedFraction -> Double
+toDouble = last . unstableConvergents . fractionalConvergents
+
+instance Show ContinuedFraction where
+  show = show . toDouble
+
+fromRealFrac :: RealFrac a => a -> ContinuedFraction
+fromRealFrac x =
+    ContinuedFraction (q0, go y)
+  where
+    go s =
+      if s == 0.0 then []
+      else let (q,t) = properFraction (1.0 / s) in q : go t
+
+    (q0,y) = properFraction x
+
+invert :: ContinuedFraction -> Maybe ContinuedFraction
+invert (ContinuedFraction (q0,qs)) =
+    if q0 /= 0 then Just (ContinuedFraction (0, q0 : qs))
+    else
+      case qs of
+        [] -> Nothing
+        h : t -> Just (ContinuedFraction (h,t))
diff --git a/src/Arithmetic/Modular.hs b/src/Arithmetic/Modular.hs
--- a/src/Arithmetic/Modular.hs
+++ b/src/Arithmetic/Modular.hs
@@ -11,6 +11,7 @@
 where
 
 import OpenTheory.Primitive.Natural
+import OpenTheory.Natural.Divides
 import qualified OpenTheory.Natural.Bits as Bits
 
 multiplyExponential :: (a -> a -> a) -> a -> a -> Natural -> a
@@ -60,3 +61,15 @@
 
 exp2 :: Natural -> Natural -> Natural -> Natural
 exp2 n x k = functionPower (square n) k x
+
+invert :: Natural -> Natural -> Maybe Natural
+invert n x =
+    if g == 1 then Just s else Nothing
+  where
+    (g,(s,_)) = egcd x n
+
+divide :: Natural -> Natural -> Natural -> Maybe Natural
+divide n x y =
+    case invert n y of
+      Nothing -> Nothing
+      Just z -> Just (multiply n x z)
diff --git a/src/Arithmetic/Smooth.hs b/src/Arithmetic/Smooth.hs
new file mode 100644
--- /dev/null
+++ b/src/Arithmetic/Smooth.hs
@@ -0,0 +1,86 @@
+{- |
+module: Arithmetic.Smooth
+description: Smooth numbers
+license: MIT
+
+maintainer: Joe Leslie-Hurd <joe@gilith.com>
+stability: provisional
+portability: portable
+-}
+module Arithmetic.Smooth
+where
+
+import qualified Data.List as List
+import OpenTheory.Primitive.Natural
+import qualified OpenTheory.Natural.Bits as Bits
+import OpenTheory.Natural.Divides
+import qualified OpenTheory.Natural.Prime as Prime
+
+factorOut :: Natural -> Natural -> Maybe (Natural,Natural)
+factorOut p =
+    go 0
+  where
+    go k n =
+      if divides p n then go (k + 1) (n `div` p)
+      else if k == 0 then Nothing
+      else Just (k,n)
+
+factorList :: [Natural] -> Natural -> ([(Natural,Natural)],Natural)
+factorList ps n =
+    case ps of
+      [] -> ([],n)
+      p : pt ->
+        case factorOut p n of
+	  Nothing -> factorList pt n
+	  Just (k,m) ->
+            let (pks,q) = factorList pt m in
+            ((p,k) : pks, q)
+
+factorBase :: Natural -> Natural -> ([(Natural,Natural)],Natural)
+factorBase k = factorList (take (fromIntegral k) Prime.primes)
+
+multiplyBase :: ([(Natural,Natural)],Natural) -> Natural
+multiplyBase =
+    \(pks,m) -> foldr mult m pks
+  where
+    mult (p,k) m = (p ^ k) * m
+
+newtype Smooth =
+    Smooth {unSmooth :: ([(Natural,Natural)],Natural)}
+  deriving (Eq,Ord)
+
+instance Show Smooth where
+  show s =
+      if null factors then "1" else List.intercalate "*" factors
+    where
+      factors = map showPk pks ++ showRest
+      showRest = if n == 1 then [] else [showWidth]
+      showWidth = if w < 20 then show n
+                  else "[" ++ show w ++ "]"
+      showPk (p,k) = show p ++ showExp k
+      showExp k = if k == 1 then "" else "^" ++ show k
+      (pks,n) = unSmooth s
+      w = Bits.width n
+
+fromNatural :: Natural -> Natural -> Smooth
+fromNatural k = Smooth . factorBase k
+
+toNatural :: Smooth -> Natural
+toNatural = multiplyBase . unSmooth
+
+factoring :: Smooth -> Maybe [(Natural,Natural)]
+factoring s =
+    if n == 1 then Just pks else Nothing
+  where
+    (pks,n) = unSmooth s
+
+next :: Natural -> Natural -> Smooth
+next k =
+    go
+  where
+    go n =
+        case factoring s of
+          Nothing -> go (n + 1)
+          Just _ -> s
+      where
+        s = fromNatural k n
diff --git a/src/Arithmetic/SquareRoot.hs b/src/Arithmetic/SquareRoot.hs
new file mode 100644
--- /dev/null
+++ b/src/Arithmetic/SquareRoot.hs
@@ -0,0 +1,72 @@
+{- |
+module: Arithmetic.SquareRoot
+description: Natural number square root
+license: MIT
+
+maintainer: Joe Leslie-Hurd <joe@gilith.com>
+stability: provisional
+portability: portable
+-}
+module Arithmetic.SquareRoot
+where
+
+import OpenTheory.Primitive.Natural
+import qualified Data.List as List
+
+import qualified Arithmetic.ContinuedFraction as ContinuedFraction
+
+floor :: Natural -> Natural
+floor n =
+    if n < 2 then n else bisect 0 n
+  where
+    bisect l u =
+        if m == l then l
+	else if m * m <= n then bisect m u
+	else bisect l m
+      where
+        m = (l + u) `div` 2
+
+ceiling :: Natural -> Natural
+ceiling n =
+    if sqrtn * sqrtn == n then sqrtn else sqrtn + 1
+  where
+    sqrtn = Arithmetic.SquareRoot.floor n
+
+continuedFraction :: Natural -> ContinuedFraction.ContinuedFraction
+continuedFraction n =
+    ContinuedFraction.ContinuedFraction (sqrtn,qs)
+  where
+    sqrtn = Arithmetic.SquareRoot.floor n
+
+    ps = continuedFractionPeriodicTail n sqrtn
+
+    qs = if null ps then [] else cycle ps
+
+continuedFractionPeriodic :: Natural -> [Natural]
+continuedFractionPeriodic n =
+    continuedFractionPeriodicTail n sqrtn
+  where
+    sqrtn = Arithmetic.SquareRoot.floor n
+
+continuedFractionPeriodicTail :: Natural -> Natural -> [Natural]
+continuedFractionPeriodicTail n sqrtn =
+    List.unfoldr go (sqrtn,sqrtd)
+  where
+    sqrtd = n - sqrtn * sqrtn
+
+-- (sqrt(n) + a) / b = c + 1 / x ==>
+-- x = b / (sqrt(n) + a - c * b)
+--   = b / (sqrt(n) - (c * b - a))
+--   = (b * (sqrt(n) + (c * b - a))) / (n - (c * b - a)^2)
+    advance (a,b) =
+        (c,(d,e))
+      where
+        c = (sqrtn + a) `div` b
+        d = c * b - a
+        e = (n - d * d) `div` b
+
+    go (a,b) =
+        case b of
+          0 -> Nothing
+          1 -> Just (2 * a, (0,0))
+          _ -> Just (advance (a,b))
diff --git a/src/IntegerDivides.hs b/src/IntegerDivides.hs
deleted file mode 100644
--- a/src/IntegerDivides.hs
+++ /dev/null
@@ -1,31 +0,0 @@
-{- |
-module: IntegerDivides
-description: Integer division algorithms
-license: MIT
-
-maintainer: Joe Leslie-Hurd <joe@gilith.com>
-stability: provisional
-portability: portable
--}
-module IntegerDivides
-where
-
-divides :: Integer -> Integer -> Bool
-divides 0 b = b == 0
-divides a b = abs b `mod` abs a == 0
-
-egcd :: Integer -> Integer -> (Integer,(Integer,Integer))
-egcd a 0 = (a,(1,0))
-egcd a b =
-    (g, (t, s - (a `div` b) * t))
-  where
-    (g,(s,t)) = egcd b (a `mod` b)
-
-chineseRemainder :: Integer -> Integer -> Integer -> Integer -> Integer
-chineseRemainder a b =
-    \x y -> (x * tb + y * sa) `mod` ab
-  where
-    (_,(s,t)) = egcd a b
-    ab = a * b
-    sa = s * a
-    tb = t * b
diff --git a/src/NaturalDivides.hs b/src/NaturalDivides.hs
deleted file mode 100644
--- a/src/NaturalDivides.hs
+++ /dev/null
@@ -1,37 +0,0 @@
-{- |
-module: NaturalDivides
-description: Natural number division algorithms
-license: MIT
-
-maintainer: Joe Leslie-Hurd <joe@gilith.com>
-stability: provisional
-portability: portable
--}
-module NaturalDivides
-where
-
-import OpenTheory.Primitive.Natural
-
-divides :: Natural -> Natural -> Bool
-divides 0 b = b == 0
-divides a b = b `mod` a == 0
-
-egcd :: Natural -> Natural -> (Natural,(Natural,Natural))
-egcd a 0 = (a,(1,0))
-egcd a b =
-    if c == 0
-    then (b, (1, a `div` b - 1))
-    else (g, (u, t + (a `div` b) * u))
-  where
-    c = a `mod` b
-    (g,(s,t)) = egcd c (b `mod` c)
-    u = s + (b `div` c) * t
-
-chineseRemainder :: Natural -> Natural -> Natural -> Natural -> Natural
-chineseRemainder a b =
-    \x y -> (x * tb + y * sa) `mod` ab
-  where
-    (_,(s,t)) = egcd a b
-    ab = a * b
-    sa = s * a
-    tb = (a - t) * b
diff --git a/src/Test.hs b/src/Test.hs
--- a/src/Test.hs
+++ b/src/Test.hs
@@ -14,72 +14,97 @@
 import qualified Test.QuickCheck as QuickCheck
 import OpenTheory.Primitive.Natural
 import OpenTheory.Natural
+import OpenTheory.Natural.Divides
 import qualified OpenTheory.Primitive.Random as Random
 import qualified OpenTheory.Natural.Uniform as Uniform
 import OpenTheory.Primitive.Test
 
-import qualified IntegerDivides
-import qualified NaturalDivides
 import Arithmetic.Random
 import Arithmetic.Prime
+import qualified Arithmetic.ContinuedFraction as ContinuedFraction
 import qualified Arithmetic.Modular as Modular
 import qualified Arithmetic.Montgomery as Montgomery
+import qualified Arithmetic.Smooth as Smooth
+import qualified Arithmetic.SquareRoot as SquareRoot
 
-propIntegerEgcdDivides :: Integer -> Integer -> Bool
-propIntegerEgcdDivides a b =
-    let (g,_) = IntegerDivides.egcd a b in
-    IntegerDivides.divides g a && IntegerDivides.divides g b
+propEgcdDivides :: Natural -> Natural -> Bool
+propEgcdDivides a b =
+    divides g a && divides g b
+  where
+    (g,_) = egcd a b
 
-propIntegerEgcdEquation :: Integer -> Integer -> Bool
-propIntegerEgcdEquation a b =
-    let (g,(s,t)) = IntegerDivides.egcd a b in
-    s * a + t * b == g
+propEgcdEquation :: Natural -> Natural -> Bool
+propEgcdEquation ap b =
+    s * a == t * b + g
+  where
+    a = ap + 1
+    (g,(s,t)) = egcd a b
 
-propIntegerEgcdBound :: Integer -> Integer -> Bool
-propIntegerEgcdBound a b =
-    let (_,(s,t)) = IntegerDivides.egcd a b in
-    abs s <= max ((abs b + 1) `div` 2) 1 &&
-    abs t <= max ((abs a + 1) `div` 2) 1
+propEgcdBound :: Natural -> Natural -> Bool
+propEgcdBound ap b =
+    s < max b 2 && t < a
+  where
+    a = ap + 1
+    (_,(s,t)) = egcd a b
 
-propNaturalEgcdDivides :: Natural -> Natural -> Bool
-propNaturalEgcdDivides a b =
-    let (g,_) = NaturalDivides.egcd a b in
-    NaturalDivides.divides g a && NaturalDivides.divides g b
+propSmoothInjective :: Natural -> Natural -> Bool
+propSmoothInjective k np =
+    Smooth.toNatural (Smooth.fromNatural k n) == n
+  where
+    n = np + 1
 
-propNaturalEgcdEquation :: Natural -> Natural -> Bool
-propNaturalEgcdEquation ap b =
-    let a = ap + 1 in
-    let (g,(s,t)) = NaturalDivides.egcd a b in
-    s * a == t * b + g
+propFloorSqrt :: Natural -> Bool
+propFloorSqrt n =
+    sq s <= n && n < sq (s + 1)
+  where
+    s = SquareRoot.floor n
+    sq i = i * i
 
-propNaturalEgcdBound :: Natural -> Natural -> Bool
-propNaturalEgcdBound ap b =
-    let a = ap + 1 in
-    let (_,(s,t)) = NaturalDivides.egcd a b in
-    s < max b 2 && t < a
+propCeilingSqrt :: Natural -> Bool
+propCeilingSqrt n =
+    (s == 0 || sq (s - 1) < n) && n <= sq s
+  where
+    s = SquareRoot.ceiling n
+    sq i = i * i
 
-propIntegerChineseRemainder :: Int -> Random.Random -> Bool
-propIntegerChineseRemainder w r =
-    n `mod` a == x && n `mod` b == y && n < a * b
+propContinuedFractionSqrt :: Natural -> Bool
+propContinuedFractionSqrt n =
+    cf == spec
   where
-    (a,b) = randomCoprimeInteger w r1
-    x = uniformInteger a r2
-    y = uniformInteger b r3
-    n = IntegerDivides.chineseRemainder a b x y
-    (r1,r23) = Random.split r
-    (r2,r3) = Random.split r23
+    cf = ContinuedFraction.toDouble (SquareRoot.continuedFraction n)
+    spec = sqrt (fromIntegral n)
 
-propNaturalChineseRemainder :: Int -> Random.Random -> Bool
-propNaturalChineseRemainder w r =
+propChineseRemainder :: Int -> Random.Random -> Bool
+propChineseRemainder w r =
     n `mod` a == x && n `mod` b == y && n < a * b
   where
     (a,b) = randomCoprime w r1
     x = Uniform.random a r2
     y = Uniform.random b r3
-    n = NaturalDivides.chineseRemainder a b x y
+    n = chineseRemainder a b x y
     (r1,r23) = Random.split r
     (r2,r3) = Random.split r23
 
+propModularNegate :: Int -> Random.Random -> Bool
+propModularNegate nw rnd =
+    Modular.add n a b == 0 &&
+    b < n
+  where
+    n = randomWidth nw r1
+    a = Uniform.random n r2
+    b = Modular.negate n a
+    (r1,r2) = Random.split rnd
+
+propModularInvert :: Int -> Random.Random -> Bool
+propModularInvert nw rnd =
+    case Modular.invert n a of
+      Nothing -> True
+      Just b -> Modular.multiply n a b == 1 && b < n
+  where
+    n = randomWidth nw r1
+    a = Uniform.random n r2
+    (r1,r2) = Random.split rnd
+
 randomMontgomeryParameters :: Int -> Random.Random -> Montgomery.Parameters
 randomMontgomeryParameters w r = Montgomery.standardParameters (randomOdd w r)
 
@@ -281,35 +306,35 @@
 
 checkWidthProps :: Int -> IO ()
 checkWidthProps w =
-   do checkWidthProp w "Check integer Chinese remainder properties"
-        propIntegerChineseRemainder
-      checkWidthProp w "Check natural Chinese remainder properties"
-        propNaturalChineseRemainder
-      checkWidthProp w "Check Montgomery invariant" propMontgomeryInvariant
-      checkWidthProp w "Check Montgomery normalize" propMontgomeryNormalize
-      checkWidthProp w "Check Montgomery reduce" propMontgomeryReduce
-      checkWidthProp w "Check Montgomery reduce small" propMontgomeryReduceSmall
-      checkWidthProp w "Check Montgomery toNatural" propMontgomeryToNatural
-      checkWidthProp w "Check Montgomery fromNatural" propMontgomeryFromNatural
-      checkWidthProp w "Check Montgomery zero" propMontgomeryZero
-      checkWidthProp w "Check Montgomery one" propMontgomeryOne
-      checkWidthProp w "Check Montgomery two" propMontgomeryTwo
-      checkWidthProp w "Check Montgomery add" propMontgomeryAdd
-      checkWidthProp w "Check Montgomery negate" propMontgomeryNegate
-      checkWidthProp w "Check Montgomery multiply" propMontgomeryMultiply
-      checkWidthProp w "Check Montgomery modexp" propMontgomeryModexp
-      checkWidthProp w "Check Montgomery modexp2" propMontgomeryModexp2
+   do checkWidthProp w "Chinese remainder" propChineseRemainder
+      checkWidthProp w "Modular negate" propModularNegate
+      checkWidthProp w "Modular invert" propModularInvert
+      checkWidthProp w "Montgomery invariant" propMontgomeryInvariant
+      checkWidthProp w "Montgomery normalize" propMontgomeryNormalize
+      checkWidthProp w "Montgomery reduce" propMontgomeryReduce
+      checkWidthProp w "Montgomery reduce small" propMontgomeryReduceSmall
+      checkWidthProp w "Montgomery toNatural" propMontgomeryToNatural
+      checkWidthProp w "Montgomery fromNatural" propMontgomeryFromNatural
+      checkWidthProp w "Montgomery zero" propMontgomeryZero
+      checkWidthProp w "Montgomery one" propMontgomeryOne
+      checkWidthProp w "Montgomery two" propMontgomeryTwo
+      checkWidthProp w "Montgomery add" propMontgomeryAdd
+      checkWidthProp w "Montgomery negate" propMontgomeryNegate
+      checkWidthProp w "Montgomery multiply" propMontgomeryMultiply
+      checkWidthProp w "Montgomery modexp" propMontgomeryModexp
+      checkWidthProp w "Montgomery modexp2" propMontgomeryModexp2
       checkWidthProp w "Fermat's little theorem" propFermat
       return ()
 
 main :: IO ()
 main =
-    do check "Check integer egcd divides\n  " propIntegerEgcdDivides
-       check "Check integer egcd equation\n  " propIntegerEgcdEquation
-       check "Check integer egcd bound\n  " propIntegerEgcdBound
-       check "Check natural egcd divides\n  " propNaturalEgcdDivides
-       check "Check natural egcd equation\n  " propNaturalEgcdEquation
-       check "Check natural egcd bound\n  " propNaturalEgcdBound
+    do check "Check egcd divides\n  " propEgcdDivides
+       check "Check egcd equation\n  " propEgcdEquation
+       check "Check egcd bound\n  " propEgcdBound
+       check "Check smooth injective\n  " propSmoothInjective
+       check "Check floor square root\n  " propFloorSqrt
+       check "Check ceiling square root\n  " propCeilingSqrt
+       check "Check continued fraction square root\n  " propContinuedFractionSqrt
        mapM_ checkWidthProps ws
        return ()
   where