packages feed

arithmetic (empty) → 1.0

raw patch · 11 files changed

+1064/−0 lines, 11 filesdep +QuickCheckdep +basedep +opentheorysetup-changed

Dependencies added: QuickCheck, base, opentheory, opentheory-bits, opentheory-divides, opentheory-primitive, random

Files

+ LICENSE view
@@ -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.
+ Setup.hs view
@@ -0,0 +1,6 @@+module Main ( main ) where++import Distribution.Simple++main :: IO ()+main = defaultMain
+ arithmetic.cabal view
@@ -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
+ src/Arithmetic/Modular.hs view
@@ -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
+ src/Arithmetic/Montgomery.hs view
@@ -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
+ src/Arithmetic/Prime.hs view
@@ -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
+ src/Arithmetic/Random.hs view
@@ -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
+ src/IntegerDivides.hs view
@@ -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
+ src/Main.hs view
@@ -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 ()
+ src/NaturalDivides.hs view
@@ -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
+ src/Test.hs view
@@ -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)