diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,17 @@
+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/Setup.hs b/Setup.hs
new file mode 100644
--- /dev/null
+++ b/Setup.hs
@@ -0,0 +1,6 @@
+module Main ( main ) where
+
+import Distribution.Simple
+
+main :: IO ()
+main = defaultMain
diff --git a/arithmetic.cabal b/arithmetic.cabal
new file mode 100644
--- /dev/null
+++ b/arithmetic.cabal
@@ -0,0 +1,60 @@
+name: arithmetic
+version: 1.0
+category: Number Theory
+synopsis: Natural number arithmetic
+license: MIT
+license-file: LICENSE
+cabal-version: >= 1.8.0.2
+build-type: Simple
+author: Joe Leslie-Hurd <joe@gilith.com>
+maintainer: Joe Leslie-Hurd <joe@gilith.com>
+description:
+  This package implements a library of natural number arithmetic functions,
+  including Montgomery multiplication.
+
+Library
+  build-depends:
+    base >= 4.0 && < 5.0,
+    random >= 1.0.1.1 && < 2.0,
+    QuickCheck >= 2.4.0.1 && < 3.0,
+    opentheory-primitive >= 1.0 && < 2.0,
+    opentheory >= 1.0 && < 2.0,
+    opentheory-bits >= 1.0 && < 2.0,
+    opentheory-divides >= 1.0 && < 2.0
+  hs-source-dirs: src
+  ghc-options: -Wall
+  exposed-modules:
+    Arithmetic.Modular,
+    Arithmetic.Montgomery,
+    Arithmetic.Prime,
+    Arithmetic.Random
+
+executable arithmetic
+  build-depends:
+    base >= 4.0 && < 5.0,
+    random >= 1.0.1.1 && < 2.0,
+    QuickCheck >= 2.4.0.1 && < 3.0,
+    opentheory-primitive >= 1.0 && < 2.0,
+    opentheory >= 1.0 && < 2.0,
+    opentheory-bits >= 1.0 && < 2.0,
+    opentheory-divides >= 1.0 && < 2.0
+  hs-source-dirs: src
+  ghc-options: -Wall
+  main-is: Main.hs
+
+test-suite arithmetic-test
+  type: exitcode-stdio-1.0
+  build-depends:
+    base >= 4.0 && < 5.0,
+    random >= 1.0.1.1 && < 2.0,
+    QuickCheck >= 2.4.0.1 && < 3.0,
+    opentheory-primitive >= 1.0 && < 2.0,
+    opentheory >= 1.0 && < 2.0,
+    opentheory-bits >= 1.0 && < 2.0,
+    opentheory-divides >= 1.0 && < 2.0
+  hs-source-dirs: src
+  ghc-options: -Wall
+  main-is: Test.hs
+  other-modules:
+    IntegerDivides,
+    NaturalDivides
diff --git a/src/Arithmetic/Modular.hs b/src/Arithmetic/Modular.hs
new file mode 100644
--- /dev/null
+++ b/src/Arithmetic/Modular.hs
@@ -0,0 +1,62 @@
+{- |
+module: Arithmetic.Modular
+description: Modular arithmetic
+license: MIT
+
+maintainer: Joe Leslie-Hurd <joe@gilith.com>
+stability: provisional
+portability: portable
+-}
+module Arithmetic.Modular
+where
+
+import OpenTheory.Primitive.Natural
+import qualified OpenTheory.Natural.Bits as Bits
+
+multiplyExponential :: (a -> a -> a) -> a -> a -> Natural -> a
+multiplyExponential mult =
+    multExp
+  where
+    multExp z x k =
+        if k == 0 then z else multExp z' x' k'
+      where
+        z' = if Bits.headBits k then mult z x else z
+        x' = mult x x
+        k' = Bits.tailBits k
+
+functionPower :: (a -> a) -> Natural -> a -> a
+functionPower f =
+    loop
+  where
+    loop n x =
+       if n == 0 then x
+       else let x' = f x in x' `seq` loop (n - 1) x'
+
+normalize :: Natural -> Natural -> Natural
+normalize n x = x `mod` n
+
+add :: Natural -> Natural -> Natural -> Natural
+add n x y = normalize n (x + y)
+
+negate :: Natural -> Natural -> Natural
+negate n x =
+    if y == 0 then y else n - y
+  where
+    y = normalize n x
+
+subtract :: Natural -> Natural -> Natural -> Natural
+subtract n x y =
+    if y <= x then normalize n (x - y)
+    else Arithmetic.Modular.negate n (y - x)
+
+multiply :: Natural -> Natural -> Natural -> Natural
+multiply n x y = normalize n (x * y)
+
+square :: Natural -> Natural -> Natural
+square n x = multiply n x x
+
+exp :: Natural -> Natural -> Natural -> Natural
+exp n = multiplyExponential (multiply n) 1
+
+exp2 :: Natural -> Natural -> Natural -> Natural
+exp2 n x k = functionPower (square n) k x
diff --git a/src/Arithmetic/Montgomery.hs b/src/Arithmetic/Montgomery.hs
new file mode 100644
--- /dev/null
+++ b/src/Arithmetic/Montgomery.hs
@@ -0,0 +1,168 @@
+{- |
+module: Arithmetic.Montgomery
+description: Modular arithmetic using Montgomery multiplication
+license: MIT
+
+maintainer: Joe Leslie-Hurd <joe@gilith.com>
+stability: provisional
+portability: portable
+-}
+module Arithmetic.Montgomery
+where
+
+import OpenTheory.Primitive.Natural
+import qualified OpenTheory.Natural.Bits as Bits
+import OpenTheory.Natural.Divides
+
+import qualified Arithmetic.Modular as Modular
+
+data Parameters = Parameters
+    {nParameters :: Natural,
+     wParameters :: Natural,
+     sParameters :: Natural,
+     kParameters :: Natural,
+     rParameters :: Natural,
+     r2Parameters :: Natural,
+     zParameters :: Natural}
+  deriving Show
+
+data Montgomery = Montgomery
+    {pMontgomery :: Parameters,
+     nMontgomery :: Natural}
+  deriving Show
+
+align :: Natural -> Natural -> Natural
+align b n = if n == 0 then 0 else (((n - 1) `div` b) + 1) * b
+
+customParameters :: Natural -> Natural -> Parameters
+customParameters n w =
+    Parameters
+      {nParameters = n,
+       wParameters = w,
+       sParameters = s,
+       kParameters = k,
+       rParameters = r,
+       r2Parameters = r2,
+       zParameters = z}
+  where
+    w2 = shiftLeft 1 w
+    (_,(s,k)) = egcd w2 n
+    r = w2 `mod` n
+    r2 = (r * r) `mod` n
+    z = w2 + n - r
+
+alignedParameters :: Natural -> Natural -> Parameters
+alignedParameters b n = customParameters n (align b (Bits.width n))
+
+standardParameters :: Natural -> Parameters
+standardParameters = alignedParameters 64
+
+-- normalize p a `mod` n = a `mod` n
+-- normalize p a < 2 ^ w
+normalize :: Parameters -> Natural -> Montgomery
+normalize p =
+    loop
+  where
+    w = wParameters p
+    r = rParameters p
+
+    loop a =
+        if x == 0 then
+          Montgomery
+            {pMontgomery = p,
+             nMontgomery = a}
+        else
+          loop ((a - shiftLeft x w) + x * r)
+      where
+        x = shiftRight a w
+
+-- normalize1 p a `mod` n = a `mod` n
+-- a < 2 ^ w + n ==> normalize1 p a < 2 ^ w
+normalize1 :: Parameters -> Natural -> Montgomery
+normalize1 p a =
+    Montgomery {pMontgomery = p, nMontgomery = b}
+  where
+    n = nParameters p
+    w = wParameters p
+    b = if Bits.bit a w then a - n else a
+
+-- reduce p a `mod` n = (a * s) `mod` n
+-- a <= r * x ==> reduce p a < x + n
+reduce :: Parameters -> Natural -> Natural
+reduce p a =
+    shiftRight (a + Bits.bound (a * k) w * n) w
+  where
+    n = nParameters p
+    w = wParameters p
+    k = kParameters p
+
+toNatural :: Montgomery -> Natural
+toNatural a =
+    if b < n then b else 0
+  where
+    p = pMontgomery a
+    n = nParameters p
+    b = reduce p (nMontgomery a)
+
+fromNatural :: Parameters -> Natural -> Montgomery
+fromNatural p =
+    multiply r2 . normalize p
+  where
+    r2 = Montgomery {pMontgomery = p, nMontgomery = r2Parameters p}
+
+zero :: Parameters -> Montgomery
+zero p = Montgomery {pMontgomery = p, nMontgomery = 0}
+
+one :: Parameters -> Montgomery
+one p = Montgomery {pMontgomery = p, nMontgomery = rParameters p}
+
+two :: Parameters -> Montgomery
+two p = double (one p)
+
+add :: Montgomery -> Montgomery -> Montgomery
+add a b = normalize (pMontgomery a) (nMontgomery a + nMontgomery b)
+
+double :: Montgomery -> Montgomery
+double a = add a a
+
+negate :: Montgomery -> Montgomery
+negate a =
+    normalize1 p (z - nMontgomery a)
+  where
+    p = pMontgomery a
+    z = zParameters p
+
+subtract :: Montgomery -> Montgomery -> Montgomery
+subtract a b = add a (Arithmetic.Montgomery.negate b)
+
+multiply :: Montgomery -> Montgomery -> Montgomery
+multiply a b =
+    normalize1 p (reduce p (nMontgomery a * nMontgomery b))
+  where
+    p = pMontgomery a
+
+square :: Montgomery -> Montgomery
+square a = multiply a a
+
+exp :: Montgomery -> Natural -> Montgomery
+exp a =
+    Modular.multiplyExponential multiply (one p) a
+  where
+    p = pMontgomery a
+
+exp2 :: Montgomery -> Natural -> Montgomery
+exp2 a k = Modular.functionPower square k a
+
+modexp :: Natural -> Natural -> Natural -> Natural
+modexp n a k =
+    toNatural m
+  where
+    p = standardParameters n
+    m = Arithmetic.Montgomery.exp (fromNatural p a) k
+
+modexp2 :: Natural -> Natural -> Natural -> Natural
+modexp2 n a k =
+    toNatural m
+  where
+    p = standardParameters n
+    m = exp2 (fromNatural p a) k
diff --git a/src/Arithmetic/Prime.hs b/src/Arithmetic/Prime.hs
new file mode 100644
--- /dev/null
+++ b/src/Arithmetic/Prime.hs
@@ -0,0 +1,80 @@
+{- |
+module: Arithmetic.Prime
+description: Generating random primes
+license: MIT
+
+maintainer: Joe Leslie-Hurd <joe@gilith.com>
+stability: provisional
+portability: portable
+-}
+module Arithmetic.Prime
+where
+
+import OpenTheory.Primitive.Natural
+import OpenTheory.Primitive.Random as Random
+import OpenTheory.Natural
+import qualified OpenTheory.Natural.Bits as Bits
+import qualified OpenTheory.Natural.Uniform as Uniform
+
+import Arithmetic.Random
+import qualified Arithmetic.Modular as Modular
+
+factorTwos :: Natural -> (Int,Natural)
+factorTwos n =
+   if Bits.headBits n then (0,n) else (r + 1, s)
+  where
+    (r,s) = factorTwos (Bits.tailBits n)
+
+millerRabinWitness :: Natural -> Natural -> Bool
+millerRabinWitness n =
+    \a -> witness (Modular.exp n a s) r
+  where
+    witness x i =
+        if i == 0 then x /= 1
+        else if x2 == 1 then not (x == 1 || x == n1)
+        else witness x2 (i - 1)
+      where
+        x2 = Modular.square n x
+
+    (r,s) = factorTwos n1
+
+    n1 = n - 1
+
+millerRabin :: Int -> Natural -> Random.Random -> Bool
+millerRabin t n =
+    \r -> n == 2 || (n /= 1 && naturalOdd n && trials t r)
+  where
+    trials i r =
+        i == 0 || (trial r1 && trials (i - 1) r2)
+      where
+        (r1,r2) = Random.split r
+
+    trial = not . millerRabinWitness n . range
+
+    range r = Uniform.random (n - 3) r + 2
+
+isPrime :: Natural -> Random.Random -> Bool
+isPrime = millerRabin 100
+
+previousPrime :: Natural -> Random.Random -> Natural
+previousPrime n r =
+    if isPrime n r1 then n else previousPrime (n - 2) r2
+  where
+    (r1,r2) = Random.split r
+
+randomPrime :: Int -> Random.Random -> Natural
+randomPrime w =
+    loop
+  where
+    loop r =
+        case oddPrime r1 of
+          Nothing -> loop r2
+          Just n -> n
+      where
+        (r1,r2) = Random.split r
+
+    oddPrime r =
+        if isPrime n r2 then Just n else Nothing
+      where
+        n = randomOdd w r1
+        (r1,r2) = Random.split r
diff --git a/src/Arithmetic/Random.hs b/src/Arithmetic/Random.hs
new file mode 100644
--- /dev/null
+++ b/src/Arithmetic/Random.hs
@@ -0,0 +1,55 @@
+{- |
+module: Arithmetic.Random
+description: Generating random natural numbers of a given width
+license: MIT
+
+maintainer: Joe Leslie-Hurd <joe@gilith.com>
+stability: provisional
+portability: portable
+-}
+module Arithmetic.Random
+where
+
+import Data.Bits
+import OpenTheory.Primitive.Natural
+import OpenTheory.Primitive.Random as Random
+import qualified OpenTheory.Natural.Bits as Bits
+import OpenTheory.Natural.Divides
+import qualified OpenTheory.Natural.Uniform as Uniform
+
+randomWidth :: Int -> Random.Random -> Natural
+randomWidth w r =
+    n + Uniform.random n r
+  where
+    n = shiftL 1 (w - 1)
+
+randomOdd :: Int -> Random.Random -> Natural
+randomOdd w r = Bits.cons True (randomWidth (w - 1) r)
+
+randomCoprime :: Int -> Random.Random -> (Natural,Natural)
+randomCoprime w =
+    loop
+  where
+    loop r =
+        case gen r1 of
+          Just ab -> ab
+          Nothing -> loop r2
+      where
+        (r1,r2) = Random.split r
+
+    gen r =
+        if g == 1 then Just (a,b) else Nothing
+      where
+        a = randomWidth w r1
+        b = randomWidth w r2
+        (g,_) = egcd a b
+        (r1,r2) = Random.split r
+
+uniformInteger :: Integer -> Random.Random -> Integer
+uniformInteger n r = fromIntegral (Uniform.random (fromIntegral n) r)
+
+randomCoprimeInteger :: Int -> Random.Random -> (Integer,Integer)
+randomCoprimeInteger w r =
+    (fromIntegral a, fromIntegral b)
+  where
+    (a,b) = randomCoprime w r
diff --git a/src/IntegerDivides.hs b/src/IntegerDivides.hs
new file mode 100644
--- /dev/null
+++ b/src/IntegerDivides.hs
@@ -0,0 +1,31 @@
+{- |
+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/Main.hs b/src/Main.hs
new file mode 100644
--- /dev/null
+++ b/src/Main.hs
@@ -0,0 +1,232 @@
+{- |
+module: Main
+description: Computing natural number arithmetic operations
+license: MIT
+
+maintainer: Joe Leslie-Hurd <joe@gilith.com>
+stability: provisional
+portability: portable
+-}
+module Main
+  ( main )
+where
+
+import qualified Data.List as List
+import System.Console.GetOpt
+import qualified System.Environment as Environment
+import qualified System.Random
+import OpenTheory.Primitive.Natural
+import qualified OpenTheory.Primitive.Random as Random
+import qualified OpenTheory.Natural.Uniform as Uniform
+
+import Arithmetic.Random
+import qualified Arithmetic.Modular as Modular
+import qualified Arithmetic.Montgomery as Montgomery
+
+--------------------------------------------------------------------------------
+-- Helper functions
+--------------------------------------------------------------------------------
+
+getPrefixString :: String -> (a -> String) -> [a] -> String -> a
+getPrefixString k p xs s =
+    case filter (List.isPrefixOf s . p) xs of
+      [] -> usage $ "bad " ++ k ++ " name: " ++ s
+      [x] -> x
+      _ : _ : _ -> usage $ "ambiguous " ++ k ++ " name: " ++ s
+
+setToString :: (a -> String) -> [a] -> String
+setToString p xs = "{" ++ List.intercalate "," (map p xs) ++ "}"
+
+--------------------------------------------------------------------------------
+-- Operations
+--------------------------------------------------------------------------------
+
+data Operation =
+    Modexp
+  | Timelock
+  deriving Show
+
+operations :: [Operation]
+operations = [Modexp,Timelock]
+
+operationToString :: Operation -> String
+operationToString oper =
+   case oper of
+     Modexp -> "modexp"
+     Timelock -> "timelock"
+
+stringToOperation :: String -> Operation
+stringToOperation = getPrefixString "operation" operationToString operations
+
+--------------------------------------------------------------------------------
+-- Algorithms
+--------------------------------------------------------------------------------
+
+data Algorithm =
+    Modular
+  | Montgomery
+  deriving Show
+
+algorithms :: [Algorithm]
+algorithms = [Modular,Montgomery]
+
+algorithmToString :: Algorithm -> String
+algorithmToString oper =
+   case oper of
+     Modular -> "modular"
+     Montgomery -> "montgomery"
+
+stringToAlgorithm :: String -> Algorithm
+stringToAlgorithm = getPrefixString "algorithm" algorithmToString algorithms
+
+--------------------------------------------------------------------------------
+-- Natural number inputs
+--------------------------------------------------------------------------------
+
+data InputNatural =
+    Fixed Natural
+  | Width Int
+  deriving Show
+
+stringToInputNatural :: String -> InputNatural
+stringToInputNatural s =
+    case s of
+      '[' : s' -> case reads s' of
+                    [(w,"]")] -> Width w
+                    _ -> usage "bad N argument"
+      _ -> case reads s of
+            [(n,"")] -> Fixed n
+            _ -> usage "bad N argument"
+
+uniformInputNatural :: InputNatural -> Random.Random -> Natural
+uniformInputNatural (Fixed n) _ = n
+uniformInputNatural (Width w) r = Uniform.random (2 ^ w) r
+
+oddInputNatural :: InputNatural -> Random.Random -> Natural
+oddInputNatural (Fixed n) _ = n
+oddInputNatural (Width w) r = randomOdd w r
+
+getInputs ::
+    Operation -> InputNatural -> Maybe InputNatural -> Maybe InputNatural ->
+    Random.Random -> (Natural,Natural,Natural)
+getInputs oper wn wx wk r =
+    (n,x,k)
+  where
+    n = oddInputNatural wn rn
+
+    x = case wx of
+          Nothing -> Uniform.random n rx
+          Just w -> uniformInputNatural w rx
+
+    k = case wk of
+          Nothing -> case oper of
+                       Modexp -> Uniform.random n rk
+                       Timelock -> 1000000
+          Just w -> uniformInputNatural w rk
+
+    (rn,r') = Random.split r
+    (rx,rk) = Random.split r'
+
+--------------------------------------------------------------------------------
+-- Options
+--------------------------------------------------------------------------------
+
+data Options = Options
+    {optOperation :: Operation,
+     optAlgorithm :: Algorithm,
+     optModulus :: InputNatural,
+     optBase :: Maybe InputNatural,
+     optExponent :: Maybe InputNatural}
+  deriving Show
+
+defaultOptions :: Options
+defaultOptions =
+  Options
+    {optOperation = Modexp,
+     optAlgorithm = Montgomery,
+     optModulus = Width 50,
+     optBase = Nothing,
+     optExponent = Nothing}
+
+options :: [OptDescr (Options -> Options)]
+options =
+    [Option [] ["operation"]
+       (ReqArg (\s opts -> opts {optOperation = stringToOperation s}) "OPERATION")
+       "select operation",
+     Option [] ["algorithm"]
+       (ReqArg (\s opts -> opts {optAlgorithm = stringToAlgorithm s}) "ALGORITHM")
+       "select algorithm",
+     Option [] ["modulus"]
+       (ReqArg (\s opts -> opts {optModulus = stringToInputNatural s}) "N")
+       "select modulus",
+     Option [] ["base"]
+       (ReqArg (\s opts -> opts {optBase = Just (stringToInputNatural s)}) "N")
+       "select base",
+     Option [] ["exponent"]
+       (ReqArg (\s opts -> opts {optExponent = Just (stringToInputNatural s)}) "N")
+       "select exponent"]
+
+processOptions :: [String] -> Either [String] (Options,[String])
+processOptions args =
+    case getOpt Permute options args of
+      (opts,work,[]) -> Right (foldl (flip id) defaultOptions opts, work)
+      (_,_,errs) -> Left errs
+
+processArguments :: [String] -> Options
+processArguments args =
+    case processOptions args of
+      Left errs -> usage (concat errs)
+      Right (opts,work) ->
+        case work of
+          [] -> opts
+          _ : _ -> usage "too many arguments"
+
+usage :: String -> a
+usage err =
+    error $ err ++ "\n" ++ usageInfo header options ++ footer
+  where
+    header = "Usage: modexp [OPTION...]"
+
+    footer =
+      "where OPERATION is one of " ++
+      setToString operationToString operations ++ ",\n" ++
+      "ALGORITHM is one of " ++
+      setToString algorithmToString algorithms ++ ",\n" ++
+      "and N is either a natural number or has the form [bitwidth]."
+
+--------------------------------------------------------------------------------
+-- Computation
+--------------------------------------------------------------------------------
+
+type Computation = Natural -> Natural -> Natural -> Natural
+
+computation :: Operation -> Algorithm -> Computation
+computation Modexp Modular = Modular.exp
+computation Modexp Montgomery = Montgomery.modexp
+computation Timelock Modular = Modular.exp2
+computation Timelock Montgomery = Montgomery.modexp2
+
+computationToString ::
+    Operation -> Natural -> Natural -> Natural -> Natural -> String
+computationToString Modexp n x k y =
+    "( " ++ show x ++ " ^ " ++ show k ++ " ) `mod` " ++
+    show n ++ " == " ++ show y
+computationToString Timelock n x k y =
+    "( " ++ show x ++ " ^ 2 ^ " ++ show k ++ " ) `mod` " ++
+    show n ++ " == " ++ show y
+
+--------------------------------------------------------------------------------
+-- Main program
+--------------------------------------------------------------------------------
+
+main :: IO ()
+main =
+    do args <- Environment.getArgs
+       r <- fmap Random.fromInt System.Random.randomIO
+       let opts = processArguments args
+       let oper = optOperation opts
+       let (n,x,k) = getInputs oper (optModulus opts) (optBase opts)
+                       (optExponent opts) r
+       let y = computation oper (optAlgorithm opts) n x k
+       putStrLn $ computationToString oper n x k y
+       return ()
diff --git a/src/NaturalDivides.hs b/src/NaturalDivides.hs
new file mode 100644
--- /dev/null
+++ b/src/NaturalDivides.hs
@@ -0,0 +1,37 @@
+{- |
+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
new file mode 100644
--- /dev/null
+++ b/src/Test.hs
@@ -0,0 +1,316 @@
+{- |
+module: Main
+description: Testing the modular exponentiation computation
+license: MIT
+
+maintainer: Joe Leslie-Hurd <joe@gilith.com>
+stability: provisional
+portability: portable
+-}
+module Main
+  ( main )
+where
+
+import qualified Test.QuickCheck as QuickCheck
+import OpenTheory.Primitive.Natural
+import OpenTheory.Natural
+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.Modular as Modular
+import qualified Arithmetic.Montgomery as Montgomery
+
+propIntegerEgcdDivides :: Integer -> Integer -> Bool
+propIntegerEgcdDivides a b =
+    let (g,_) = IntegerDivides.egcd a b in
+    IntegerDivides.divides g a && IntegerDivides.divides g b
+
+propIntegerEgcdEquation :: Integer -> Integer -> Bool
+propIntegerEgcdEquation a b =
+    let (g,(s,t)) = IntegerDivides.egcd a b in
+    s * a + t * b == g
+
+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
+
+propNaturalEgcdDivides :: Natural -> Natural -> Bool
+propNaturalEgcdDivides a b =
+    let (g,_) = NaturalDivides.egcd a b in
+    NaturalDivides.divides g a && NaturalDivides.divides g b
+
+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
+
+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
+
+propIntegerChineseRemainder :: Int -> Random.Random -> Bool
+propIntegerChineseRemainder w r =
+    n `mod` a == x && n `mod` b == y && n < a * b
+  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
+
+propNaturalChineseRemainder :: Int -> Random.Random -> Bool
+propNaturalChineseRemainder 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
+    (r1,r23) = Random.split r
+    (r2,r3) = Random.split r23
+
+randomMontgomeryParameters :: Int -> Random.Random -> Montgomery.Parameters
+randomMontgomeryParameters w r = Montgomery.standardParameters (randomOdd w r)
+
+propMontgomeryInvariant :: Int -> Random.Random -> Bool
+propMontgomeryInvariant nw rnd =
+    naturalOdd n &&
+    n < w2 &&
+    s * w2 == k * n + 1 &&
+    s < n &&
+    k < w2 &&
+    r == w2 `mod` n &&
+    r2 == (r * r) `mod` n &&
+    z `mod` n == 0 &&
+    w2 <= z &&
+    z < w2 + n
+  where
+    Montgomery.Parameters
+      {Montgomery.nParameters = n,
+       Montgomery.wParameters = w,
+       Montgomery.sParameters = s,
+       Montgomery.kParameters = k,
+       Montgomery.rParameters = r,
+       Montgomery.r2Parameters = r2,
+       Montgomery.zParameters = z} = randomMontgomeryParameters nw rnd
+
+    w2 = shiftLeft 1 w
+
+propMontgomeryNormalize :: Int -> Random.Random -> Bool
+propMontgomeryNormalize nw rnd =
+    b `mod` n == a `mod` n &&
+    b < w2
+  where
+    p = randomMontgomeryParameters nw r1
+    a = Uniform.random (w2 * w2) r2
+    b = Montgomery.nMontgomery (Montgomery.normalize p a)
+
+    n = Montgomery.nParameters p
+    w = Montgomery.wParameters p
+    w2 = shiftLeft 1 w
+    (r1,r2) = Random.split rnd
+
+propMontgomeryReduce :: Int -> Random.Random -> Bool
+propMontgomeryReduce nw rnd =
+    b `mod` n == (a * s) `mod` n &&
+    b < w2 + n
+  where
+    p = randomMontgomeryParameters nw r1
+    a = Uniform.random (w2 * w2) r2
+    b = Montgomery.reduce p a
+
+    n = Montgomery.nParameters p
+    w = Montgomery.wParameters p
+    s = Montgomery.sParameters p
+    w2 = shiftLeft 1 w
+    (r1,r2) = Random.split rnd
+
+propMontgomeryReduceSmall :: Int -> Random.Random -> Bool
+propMontgomeryReduceSmall nw rnd =
+    b `mod` n == (a * s) `mod` n &&
+    b <= n
+  where
+    p = randomMontgomeryParameters nw r1
+    a = Uniform.random w2 r2
+    b = Montgomery.reduce p a
+
+    n = Montgomery.nParameters p
+    w = Montgomery.wParameters p
+    s = Montgomery.sParameters p
+    w2 = shiftLeft 1 w
+    (r1,r2) = Random.split rnd
+
+propMontgomeryToNatural :: Int -> Random.Random -> Bool
+propMontgomeryToNatural nw rnd =
+    b == (a * s) `mod` n
+  where
+    p = randomMontgomeryParameters nw r1
+    a = Uniform.random w2 r2
+    b = Montgomery.toNatural (Montgomery.normalize p a)
+
+    n = Montgomery.nParameters p
+    w = Montgomery.wParameters p
+    s = Montgomery.sParameters p
+    w2 = shiftLeft 1 w
+    (r1,r2) = Random.split rnd
+
+propMontgomeryFromNatural :: Int -> Random.Random -> Bool
+propMontgomeryFromNatural nw rnd =
+    b == a `mod` n
+  where
+    p = randomMontgomeryParameters nw r1
+    a = Uniform.random (w2 * w2) r2
+    b = Montgomery.toNatural (Montgomery.fromNatural p a)
+
+    n = Montgomery.nParameters p
+    w = Montgomery.wParameters p
+    w2 = shiftLeft 1 w
+    (r1,r2) = Random.split rnd
+
+propMontgomeryZero :: Int -> Random.Random -> Bool
+propMontgomeryZero nw rnd =
+    Montgomery.toNatural (Montgomery.zero p) == 0
+  where
+    p = randomMontgomeryParameters nw rnd
+
+propMontgomeryOne :: Int -> Random.Random -> Bool
+propMontgomeryOne nw rnd =
+    Montgomery.toNatural (Montgomery.one p) == 1
+  where
+    p = randomMontgomeryParameters nw rnd
+
+propMontgomeryTwo :: Int -> Random.Random -> Bool
+propMontgomeryTwo nw rnd =
+    Montgomery.toNatural (Montgomery.two p) == 2
+  where
+    p = randomMontgomeryParameters nw rnd
+
+propMontgomeryAdd :: Int -> Random.Random -> Bool
+propMontgomeryAdd nw rnd =
+    Montgomery.toNatural c ==
+      Modular.add n (Montgomery.toNatural a) (Montgomery.toNatural b) &&
+    Montgomery.nMontgomery c < w2
+  where
+    p = randomMontgomeryParameters nw r1
+    a = Montgomery.normalize p (Uniform.random w2 r2)
+    b = Montgomery.normalize p (Uniform.random w2 r3)
+    c = Montgomery.add a b
+
+    n = Montgomery.nParameters p
+    w = Montgomery.wParameters p
+    w2 = shiftLeft 1 w
+    (r1,r23) = Random.split rnd
+    (r2,r3) = Random.split r23
+
+propMontgomeryNegate :: Int -> Random.Random -> Bool
+propMontgomeryNegate nw rnd =
+    Montgomery.toNatural b == Modular.negate n (Montgomery.toNatural a) &&
+    Montgomery.nMontgomery b < w2
+  where
+    p = randomMontgomeryParameters nw r1
+    a = Montgomery.normalize p (Uniform.random w2 r2)
+    b = Montgomery.negate a
+
+    n = Montgomery.nParameters p
+    w = Montgomery.wParameters p
+    w2 = shiftLeft 1 w
+    (r1,r2) = Random.split rnd
+
+propMontgomeryMultiply :: Int -> Random.Random -> Bool
+propMontgomeryMultiply nw rnd =
+    Montgomery.toNatural c ==
+      Modular.multiply n (Montgomery.toNatural a) (Montgomery.toNatural b) &&
+    Montgomery.nMontgomery c < w2
+  where
+    p = randomMontgomeryParameters nw r1
+    a = Montgomery.normalize p (Uniform.random w2 r2)
+    b = Montgomery.normalize p (Uniform.random w2 r3)
+    c = Montgomery.multiply a b
+
+    n = Montgomery.nParameters p
+    w = Montgomery.wParameters p
+    w2 = shiftLeft 1 w
+    (r1,r23) = Random.split rnd
+    (r2,r3) = Random.split r23
+
+propMontgomeryModexp :: Int -> Random.Random -> Bool
+propMontgomeryModexp w r =
+    Montgomery.modexp n x k == Modular.exp n x k
+  where
+    n = randomOdd w r1
+    x = Uniform.random n r2
+    k = Uniform.random n r3
+
+    (r1,r23) = Random.split r
+    (r2,r3) = Random.split r23
+
+propMontgomeryModexp2 :: Int -> Random.Random -> Bool
+propMontgomeryModexp2 w r =
+    Montgomery.modexp2 n x k == Modular.exp2 n x k
+  where
+    n = randomOdd w r1
+    x = Uniform.random n r2
+    k = Uniform.random (fromIntegral w) r3
+
+    (r1,r23) = Random.split r
+    (r2,r3) = Random.split r23
+
+propFermat :: Int -> Random.Random -> Bool
+propFermat w r =
+    Montgomery.modexp n a n == a
+  where
+    n = randomPrime w r1
+    a = Uniform.random n r2
+    (r1,r2) = Random.split r
+
+checkWidthProp ::
+    QuickCheck.Testable prop => Int -> String -> (Int -> prop) -> IO ()
+checkWidthProp w s p =
+    check (s ++ " (" ++ show w ++ " bit)\n  ") (p w)
+
+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
+      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
+       mapM_ checkWidthProps ws
+       return ()
+  where
+    ws = takeWhile (\n -> n <= 128) (iterate ((*) 2) 4)