diff --git a/arithmetic.cabal b/arithmetic.cabal
--- a/arithmetic.cabal
+++ b/arithmetic.cabal
@@ -1,5 +1,5 @@
 name: arithmetic
-version: 1.3
+version: 1.4
 category: Number Theory
 synopsis: Natural number arithmetic
 license: MIT
@@ -12,8 +12,9 @@
   This package implements a library of natural number arithmetic functions,
   including Montgomery multiplication, the Miller-Rabin primality test,
   Lucas sequences, the Williams p+1 factorization method, continued fraction
-  representations of natural number square roots, the Jacobi symbol and the
-  Tonelli-Shanks algorithm for finding square roots modulo a prime.
+  representations of natural number square roots, the Jacobi symbol, the
+  Tonelli-Shanks algorithm for finding square roots modulo a prime, and
+  the Chakravala method for solving Pell's equation.
 
 Library
   build-depends:
@@ -32,6 +33,7 @@
     Arithmetic.Lucas,
     Arithmetic.Modular,
     Arithmetic.Montgomery,
+    Arithmetic.Pell,
     Arithmetic.Polynomial,
     Arithmetic.Prime,
     Arithmetic.Prime.Factor,
diff --git a/src/Arithmetic/Pell.hs b/src/Arithmetic/Pell.hs
new file mode 100644
--- /dev/null
+++ b/src/Arithmetic/Pell.hs
@@ -0,0 +1,71 @@
+{- |
+module: Arithmetic.Pell
+description: Pell's equation (a^2 = n*b^2 + 1)
+license: MIT
+
+maintainer: Joe Leslie-Hurd <joe@gilith.com>
+stability: provisional
+portability: portable
+-}
+module Arithmetic.Pell
+where
+
+import OpenTheory.Primitive.Natural
+
+import Arithmetic.Utility
+import qualified Arithmetic.Modular as Modular
+import qualified Arithmetic.Quadratic as Quadratic
+
+-------------------------------------------------------------------------------
+-- Using the Chakravala method to find the fundamental solution of
+-- Pell's equation
+--
+--   a^2 = n*b^2 + 1
+--
+-- (where n is not a square).
+-------------------------------------------------------------------------------
+
+chakravala :: Natural -> [(Natural,Natural,Natural)]
+chakravala n =
+    if sqrtN * sqrtN == n then []
+    else let a = minM 0 1 in reduce a 1
+  where
+    reduce a b = (a,b,k) : (if k == 1 && a2 > nb2 then [] else reduce a' b')
+      where
+        a2 = a * a
+        nb2 = n * b * b
+        k =  distance a2 nb2
+        a' = (a * m + n * b) `div` k
+        b' = (a + b * m) `div` k
+        j = case Modular.divide k (Modular.negate k a) b of
+              Just i -> i
+              Nothing -> error "pell: couldn't divide"
+        m = minM j k
+
+    sqrtN = Quadratic.rootFloor n
+
+    minM j k =
+        if n - m_0 * m_0 <= m_1 * m_1 - n then m_0 else m_1
+      where
+        m_0 = sqrtN - Modular.subtract k sqrtN j
+        m_1 = m_0 + k
+
+-------------------------------------------------------------------------------
+-- Finding all integer solutions of Pell's equation
+-------------------------------------------------------------------------------
+
+solutions :: Natural -> [(Natural,Natural)]
+solutions n =
+    (a0,b0) : (if null l then [] else go a1 b1)
+  where
+    a0 = 1
+    b0 = 0
+    l = chakravala n
+    (a1,b1,_) = last l
+    go a b = (a,b) : go (a1 * a + n * b1 * b) (a1 * b + b1 * a)
+
+solution :: Natural -> (Natural,Natural)
+solution n =
+    case solutions n of
+      _ : ab : _ -> ab
+      _ -> error "Pell's equation a^2 = n*b^2 + 1 has no nontrivial integer solution when n is square"
diff --git a/src/Arithmetic/Quadratic.hs b/src/Arithmetic/Quadratic.hs
--- a/src/Arithmetic/Quadratic.hs
+++ b/src/Arithmetic/Quadratic.hs
@@ -12,6 +12,7 @@
 
 import OpenTheory.Primitive.Natural
 import qualified Data.List as List
+import qualified Data.Maybe as Maybe
 
 import Arithmetic.Utility
 import qualified Arithmetic.ContinuedFraction as ContinuedFraction
@@ -33,6 +34,15 @@
     if sqrtn * sqrtn == n then sqrtn else sqrtn + 1
   where
     sqrtn = rootFloor n
+
+destSquare :: Natural -> Maybe Natural
+destSquare n =
+    if sqrtn * sqrtn == n then Just sqrtn else Nothing
+  where
+    sqrtn = rootFloor n
+
+isSquare :: Natural -> Bool
+isSquare = Maybe.isJust . destSquare
 
 rootContinuedFraction :: Natural -> ContinuedFraction.ContinuedFraction
 rootContinuedFraction n =
diff --git a/src/Arithmetic/Utility.hs b/src/Arithmetic/Utility.hs
--- a/src/Arithmetic/Utility.hs
+++ b/src/Arithmetic/Utility.hs
@@ -14,6 +14,9 @@
 import OpenTheory.Natural.Divides
 import qualified OpenTheory.Natural.Bits as Bits
 
+distance :: Natural -> Natural -> Natural
+distance m n = if m <= n then n - m else m - n
+
 functionPower :: (a -> a) -> Natural -> a -> a
 functionPower f =
     loop
diff --git a/src/Test.hs b/src/Test.hs
--- a/src/Test.hs
+++ b/src/Test.hs
@@ -27,6 +27,7 @@
 import qualified Arithmetic.Prime.Factor as Factor
 import qualified Arithmetic.Modular as Modular
 import qualified Arithmetic.Montgomery as Montgomery
+import qualified Arithmetic.Pell as Pell
 import qualified Arithmetic.Polynomial as Polynomial
 import qualified Arithmetic.Quadratic as Quadratic
 import qualified Arithmetic.Ring as Ring
@@ -511,6 +512,12 @@
     p = Polynomial.fromCoefficients r (map (Ring.fromNatural r) ps)
     q = Polynomial.fromCoefficients r (map (Ring.fromNatural r) (qs ++ [1]))
 
+propPellEquation :: Natural -> Bool
+propPellEquation n =
+    Quadratic.isSquare n || a * a == n * b * b + 1
+  where
+    (a,b) = Pell.solution n
+
 {-
 np = (0 :: Natural)
 ps = ([] :: [Natural])
@@ -585,4 +592,5 @@
        check "Polynomial quotient remainder" propPolynomialQuotientRemainder
        check "Polynomial quotient remainder monic"
          propPolynomialQuotientRemainderMonic
+       check "Pell equation solution" propPellEquation
        return ()