arithmetic 1.1 → 1.2

raw patch · 16 files changed

+1280/−522 lines, 16 filesdep +containersdep ~opentheory-primitivePVP ok

version bump matches the API change (PVP)

Dependencies added: containers

Dependency ranges changed: opentheory-primitive

API changes (from Hackage documentation)

Files

• arithmetic.cabal +20/−11
• src/Arithmetic/Lucas.hs +28/−0
• src/Arithmetic/Modular.hs +8/−21
• src/Arithmetic/Montgomery.hs +3/−3
• src/Arithmetic/Prime.hs +44/−17
• src/Arithmetic/Prime/Factor.hs +220/−0
• src/Arithmetic/Prime/Sieve.hs +50/−0
• src/Arithmetic/Quadratic.hs +180/−0
• src/Arithmetic/Random.hs +35/−23
• src/Arithmetic/Smooth.hs +0/−86
• src/Arithmetic/SquareRoot.hs +0/−72
• src/Arithmetic/Utility.hs +46/−0
• src/Arithmetic/Utility/Heap.hs +84/−0
• src/Arithmetic/Williams.hs +117/−0
• src/Main.hs +183/−88
• src/Test.hs +262/−201

@@ -1,5 +1,5 @@ name: arithmetic-version: 1.1+version: 1.2 category: Number Theory synopsis: Natural number arithmetic license: MIT@@ -10,39 +10,47 @@ maintainer: Joe Leslie-Hurd <joe@gilith.com> description:   This package implements a library of natural number arithmetic functions,-  including Montgomery multiplication and continued fractions.+  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.  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,+    containers >= 0.4.2.1 && < 1.0,+    opentheory-primitive >= 1.8 && < 2.0,     opentheory >= 1.0 && < 2.0,     opentheory-bits >= 1.0 && < 2.0,-    opentheory-divides >= 1.0 && < 2.0,-    opentheory-prime >= 1.0 && < 2.0+    opentheory-divides >= 1.0 && < 2.0   hs-source-dirs: src   ghc-options: -Wall   exposed-modules:     Arithmetic.ContinuedFraction,+    Arithmetic.Lucas,     Arithmetic.Modular,     Arithmetic.Montgomery,     Arithmetic.Prime,+    Arithmetic.Prime.Factor,+    Arithmetic.Prime.Sieve,+    Arithmetic.Quadratic,     Arithmetic.Random,-    Arithmetic.Smooth,-    Arithmetic.SquareRoot+    Arithmetic.Utility,+    Arithmetic.Utility.Heap,+    Arithmetic.Williams  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,+    containers >= 0.4.2.1 && < 1.0,+    opentheory-primitive >= 1.8 && < 2.0,     opentheory >= 1.0 && < 2.0,     opentheory-bits >= 1.0 && < 2.0,-    opentheory-divides >= 1.0 && < 2.0,-    opentheory-prime >= 1.0 && < 2.0+    opentheory-divides >= 1.0 && < 2.0   hs-source-dirs: src   ghc-options: -Wall   main-is: Main.hs@@ -53,7 +61,8 @@     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,+    containers >= 0.4.2.1 && < 1.0,+    opentheory-primitive >= 1.8 && < 2.0,     opentheory >= 1.0 && < 2.0,     opentheory-bits >= 1.0 && < 2.0,     opentheory-divides >= 1.0 && < 2.0,

@@ -0,0 +1,28 @@+{- |+module: Arithmetic.Lucas+description: Lucas sequences+license: MIT++maintainer: Joe Leslie-Hurd <joe@gilith.com>+stability: provisional+portability: portable+-}+module Arithmetic.Lucas+where++advance :: (a -> a -> a) -> (a -> a -> a) -> a -> a -> a -> a -> a+advance sub mult p q x y = sub (mult p y) (mult q x)++sequence :: (a -> a -> a) -> (a -> a -> a) -> a -> a -> a -> a -> [a]+sequence sub mult p q =+    go+  where+    go x y = x : go y (advance sub mult p q x y)++uSequence :: a -> a -> (a -> a -> a) -> (a -> a -> a) -> a -> a -> [a]+uSequence zero one sub mult p q =+    Arithmetic.Lucas.sequence sub mult p q zero one++vSequence :: a -> (a -> a -> a) -> (a -> a -> a) -> a -> a -> [a]+vSequence two sub mult p q =+    Arithmetic.Lucas.sequence sub mult p q two p

@@ -12,26 +12,8 @@  import OpenTheory.Primitive.Natural import OpenTheory.Natural.Divides-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'+import Arithmetic.Utility  normalize :: Natural -> Natural -> Natural normalize n x = x `mod` n@@ -39,6 +21,9 @@ add :: Natural -> Natural -> Natural -> Natural add n x y = normalize n (x + y) +double :: Natural -> Natural -> Natural+double n x = add n x x+ negate :: Natural -> Natural -> Natural negate n x =     if y == 0 then y else n - y@@ -60,11 +45,13 @@ exp n = multiplyExponential (multiply n) 1  exp2 :: Natural -> Natural -> Natural -> Natural-exp2 n x k = functionPower (square n) k x+exp2 n x k = if k == 0 then normalize n x else functionPower (square n) k x  invert :: Natural -> Natural -> Maybe Natural invert n x =-    if g == 1 then Just s else Nothing+    if n == 1 then Just 0+    else if g == 1 then Just s+    else Nothing   where     (g,(s,_)) = egcd x n 

@@ -14,7 +14,7 @@ import qualified OpenTheory.Natural.Bits as Bits import OpenTheory.Natural.Divides -import qualified Arithmetic.Modular as Modular+import Arithmetic.Utility  data Parameters = Parameters     {nParameters :: Natural,@@ -146,12 +146,12 @@  exp :: Montgomery -> Natural -> Montgomery exp a =-    Modular.multiplyExponential multiply (one p) a+    multiplyExponential multiply (one p) a   where     p = pMontgomery a  exp2 :: Montgomery -> Natural -> Montgomery-exp2 a k = Modular.functionPower square k a+exp2 a k = functionPower square k a  modexp :: Natural -> Natural -> Natural -> Natural modexp n a k =

@@ -13,17 +13,15 @@ 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 Arithmetic.Utility import qualified Arithmetic.Modular as Modular+import qualified Arithmetic.Prime.Sieve as Sieve -factorTwos :: Natural -> (Int,Natural)-factorTwos n =-   if Bits.headBits n then (0,n) else (r + 1, s)-  where-    (r,s) = factorTwos (Bits.tailBits n)+primes :: [Natural]+primes = 2 : Sieve.advance 1 4 Sieve.initial  millerRabinWitness :: Natural -> Natural -> Bool millerRabinWitness n =@@ -40,9 +38,9 @@      n1 = n - 1 -millerRabin :: Int -> Natural -> Random.Random -> Bool+millerRabin :: Natural -> Natural -> Random.Random -> Bool millerRabin t n =-    \r -> n == 2 || (n /= 1 && naturalOdd n && trials t r)+    \r -> n == 2 || n == 3 || (n /= 1 && naturalOdd n && trials t r)   where     trials i r =         i == 0 || (trial r1 && trials (i - 1) r2)@@ -58,23 +56,52 @@  previousPrime :: Natural -> Random.Random -> Natural previousPrime n r =-    if isPrime n r1 then n else previousPrime (n - 2) r2+    if isPrime n r1 then n else previousPrime (n - 1) r2   where     (r1,r2) = Random.split r -randomPrime :: Int -> Random.Random -> Natural-randomPrime w =-    loop+nextPrime :: Natural -> Random.Random -> Natural+nextPrime n r =+    if isPrime n r1 then n else nextPrime (n + 1) r2   where-    loop r =-        case oddPrime r1 of-          Nothing -> loop r2-          Just n -> n+    (r1,r2) = Random.split r++nextPrime3Mod4 :: Natural -> Random.Random -> Natural+nextPrime3Mod4 =+    \n -> go ((4 * (n `div` 4)) + 3)+  where+    go n r =+        if isPrime n r1 then n else go (n + 4) r2       where         (r1,r2) = Random.split r -    oddPrime r =+nextPrime5Mod8 :: Natural -> Random.Random -> Natural+nextPrime5Mod8 =+    \n -> go ((8 * ((n + 2) `div` 8)) + 5)+  where+    go n r =+        if isPrime n r1 then n else go (n + 8) r2+      where+        (r1,r2) = Random.split r++randomPrime :: Natural -> Random.Random -> Natural+randomPrime w =+    randomMaybe gen+  where+    gen r =         if isPrime n r2 then Just n else Nothing       where         n = randomOdd w r1         (r1,r2) = Random.split r++randomPrime3Mod4 :: Natural -> Random.Random -> Natural+randomPrime3Mod4 w =+    randomFilter check (randomPrime w)+  where+    check p = p `mod` 4 == 3++randomPrime5Mod8 :: Natural -> Random.Random -> Natural+randomPrime5Mod8 w =+    randomFilter check (randomPrime w)+  where+    check p = p `mod` 8 == 5

@@ -0,0 +1,220 @@+{- |+module: Arithmetic.Prime.Factor+description: Factorized natural numbers+license: MIT++maintainer: Joe Leslie-Hurd <joe@gilith.com>+stability: provisional+portability: portable+-}+module Arithmetic.Prime.Factor+where++import OpenTheory.Primitive.Natural+import qualified OpenTheory.Natural.Bits as Bits+import qualified Data.Map as Map+import qualified Data.Maybe as Maybe+import qualified OpenTheory.Primitive.Random as Random++import Arithmetic.Prime+import Arithmetic.Random+import Arithmetic.Utility++newtype Factor = Factor {unFactor :: Map.Map Natural Natural}++one :: Factor+one = Factor {unFactor = Map.empty}++isOne :: Factor -> Bool+isOne = Map.null . unFactor++primePower :: Natural -> Natural -> Factor+primePower p k = if k == 0 then one else Factor {unFactor = Map.singleton p k}++destPrimePower :: Factor -> Maybe (Natural,Natural)+destPrimePower f =+    if Map.size m == 1 then Maybe.listToMaybe (Map.toList m) else Nothing+  where+    m = unFactor f++isPrimePower :: Factor -> Bool+isPrimePower = Maybe.isJust . destPrimePower++prime :: Natural -> Factor+prime p = primePower p 1++destPrime :: Factor -> Maybe Natural+destPrime f =+    case destPrimePower f of+      Just (p,1) -> Just p+      _ -> Nothing++isPrime :: Factor -> Bool+isPrime = Maybe.isJust . destPrime++multiply :: Factor -> Factor -> Factor+multiply f1 f2 =+    Factor {unFactor = Map.unionWith (+) m1 m2}+  where+    m1 = unFactor f1+    m2 = unFactor f2++exp :: Factor -> Natural -> Factor+exp f n =+    if n == 0 then one+    else if n == 1 then f+    else Factor {unFactor = Map.map ((*) n) (unFactor f)}++root :: Natural -> Factor -> (Factor,Factor)+root n f =+    if n == 0 then error "Arithmetic.Prime.Factor.root: n == 0"+    else if n == 1 then (f,one)+    else (fq,fr)+  where+    m = unFactor f+    fq = Factor {unFactor = Map.mapMaybe nq m}+    fr = Factor {unFactor = Map.mapMaybe nr m}+    nq k = mz (k `div` n)+    nr k = mz (k `mod` n)+    mz k = if k == 0 then Nothing else Just k++destRoot :: Natural -> Factor -> Maybe Factor+destRoot n f =+    if isOne fr then Just fq else Nothing+  where+    (fq,fr) = root n f++isRoot :: Natural -> Factor -> Bool+isRoot n = Maybe.isJust . destRoot n++gcd :: Factor -> Factor -> Factor+gcd f1 f2 =+    Factor {unFactor = Map.intersectionWith min m1 m2}+  where+    m1 = unFactor f1+    m2 = unFactor f2++trialDivision :: [Natural] -> Natural -> (Factor,Natural)+trialDivision =+    go+  where+    go [] n = (one,n)+    go (p : ps) n =+        if n <= 1 then (one,n)+        else (multiply f (primePower p r), m)+      where+        (r,s) = factorOut p n+        (f,m) = go ps s++destSmooth :: [Natural] -> Natural -> Maybe Factor+destSmooth ps n =+    if m == 1 then Just f else Nothing+  where+    (f,m) = trialDivision ps n++isSmooth :: [Natural] -> Natural -> Bool+isSmooth ps n = Maybe.isJust (destSmooth ps n)++nextSmooth :: [Natural] -> Natural -> Factor+nextSmooth ps =+    go+  where+    go n =+        case destSmooth ps n of+          Nothing -> go (n + 1)+          Just f -> f++multiplicative :: (Natural -> Natural -> a) -> (a -> a -> a) -> a -> Factor -> a+multiplicative pkA multA oneA f =+    case Map.foldrWithKey inc Nothing (unFactor f) of+      Nothing -> oneA+      Just x -> x+  where+    inc p k acc = mult (pkA p k) acc+    mult x Nothing = Just x+    mult x (Just y) = Just (multA x y)++toNatural :: Factor -> Natural+toNatural = multiplicative (^) (*) 1++totient :: Factor -> Natural+totient =+    multiplicative tot (*) 1+  where+    tot p k = (p ^ (k - 1)) * (p - 1)++instance Show Factor where+  show =+      multiplicative showPK (\s t -> s ++ " * " ++ t) "1"+    where+      showPK p k = show p ++ showExp k+      showExp k = if k == 1 then "" else "^" ++ show k++factorPower :: Natural -> Natural -> Maybe (Natural,Natural)+factorPower pmin n =+    go n primes+  where+    go _ [] = error "out of primes!"+    go s (p : ps) =+        if t < pmin then Nothing+        else if t ^ p == n then Just (t,p)+        else go t ps+      where+        t = bisect 1 s++        bisect l u =+            if m == l then l+            else if m ^ p <= n then bisect m u+            else bisect l m+          where+            m = (l + u) `div` 2++factor :: Natural -> (Natural -> Random.Random -> Maybe Natural) ->+          Natural -> Random.Random -> Maybe Factor+factor k ff =+    trial+  where+    (ptrials,pmin) = (init ps, last ps)+      where+        ps = take (fromIntegral (k + 1)) primes++    trial n rnd =+        if m == 1 then Just f+        else mmult (Just f) (go m rnd)+      where+        (f,m) = trialDivision ptrials n++    go n rnd =+        if Arithmetic.Prime.isPrime n r1 then Just (prime n)+        else+          case factorPower pmin n of+            Just (m,i) -> mexp (go m r2) i+            Nothing -> case ff n r2 of+                         Nothing -> Nothing+                         Just m -> mmult (go m r3) (go (n `div` m) r4)+      where+        (r1,r24) = Random.split rnd+        (r2,r34) = Random.split r24+        (r3,r4) = Random.split r34++    mmult (Just f1) (Just f2) = Just (multiply f1 f2)+    mmult _ _ = Nothing++    mexp (Just f) i = Just (Arithmetic.Prime.Factor.exp f i)+    mexp Nothing _ = Nothing++randomRSA :: Natural -> Random.Random -> Factor+randomRSA w =+    randomFilter check gen+  where+    check f = not (isPrimePower f) && Bits.width (toNatural f) == w++    gen rnd =+        multiply (prime p1) (prime p2)+      where+        p1 = randomPrime w1 r1+        p2 = randomPrime w2 r2+        (r1,r2) = Random.split rnd++    w1 = w `div` 2+    w2 = w - w1

@@ -0,0 +1,50 @@+{- |+module: Arithmetic.Prime.Sieve+description: The genuine sieve of Eratosphenes+license: MIT++maintainer: Joe Leslie-Hurd <joe@gilith.com>+stability: provisional+portability: portable+-}+module Arithmetic.Prime.Sieve+where++import OpenTheory.Primitive.Natural++import qualified Arithmetic.Utility.Heap as Heap++newtype Sieve = Sieve { unSieve :: Heap.Heap (Natural,Natural) }++instance Show Sieve where+  show s = show (unSieve s)++initial :: Sieve+initial =+    Sieve (Heap.empty lep)+  where+    lep (kp1,_) (kp2,_) = kp1 <= kp2++-- let p = 2 * m + 1+-- 2m' + 1 = p * p = (2m + 1) * (2m + 1) = 2(((2m + 1) + 1) * m) + 1+-- Therefore, m' = ((2m + 1) + 1) * m = (p + 1) * m+add :: Natural -> Sieve -> (Natural,Sieve)+add m (Sieve ps) =+    (p, Sieve (Heap.add (m',p) ps))+  where+    p = 2 * m + 1+    m' = (p + 1) * m++bump :: Sieve -> (Natural,Sieve)+bump (Sieve ps) =+    case Heap.remove ps of+      Nothing -> error "GenuineSieve.bump"+      Just ((kp,p),ps') -> (kp, Sieve (Heap.add (kp + p, p) ps'))++advance :: Natural -> Natural -> Sieve -> [Natural]+advance m n s =+    if m < n+      then let (p,s') = add m s in p : advance m' n s'+      else let (n',s') = bump s in advance (if m == n then m' else m) n' s'+  where+    m' = m + 1

@@ -0,0 +1,180 @@+{- |+module: Arithmetic.Quadratic+description: Natural number square root+license: MIT++maintainer: Joe Leslie-Hurd <joe@gilith.com>+stability: provisional+portability: portable+-}+module Arithmetic.Quadratic+where++import OpenTheory.Primitive.Natural+import qualified Data.List as List++import Arithmetic.Utility+import qualified Arithmetic.ContinuedFraction as ContinuedFraction+import qualified Arithmetic.Modular as Modular++rootFloor :: Natural -> Natural+rootFloor 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++rootCeiling :: Natural -> Natural+rootCeiling n =+    if sqrtn * sqrtn == n then sqrtn else sqrtn + 1+  where+    sqrtn = rootFloor n++rootContinuedFraction :: Natural -> ContinuedFraction.ContinuedFraction+rootContinuedFraction n =+    ContinuedFraction.ContinuedFraction (sqrtn,qs)+  where+    sqrtn = rootFloor n++    ps = rootContinuedFractionPeriodicTail n sqrtn++    qs = if null ps then [] else cycle ps++rootContinuedFractionPeriodic :: Natural -> [Natural]+rootContinuedFractionPeriodic n =+    rootContinuedFractionPeriodicTail n sqrtn+  where+    sqrtn = rootFloor n++rootContinuedFractionPeriodicTail :: Natural -> Natural -> [Natural]+rootContinuedFractionPeriodicTail 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))++data Residue = Residue | NonResidue | Zero+    deriving (Eq,Ord,Show)++-- The first argument (the modulus) must be an odd natural+jacobiSymbol :: Natural -> Natural -> Residue+jacobiSymbol =+    \n -> if n == 1 then const Residue else go False n+  where+    go f n m =+        if p == 0 then Zero+        else if s == 1 then if g then NonResidue else Residue+        else go h s n+      where+        p = m `mod` n+        (r,s) = factorTwos p+        n8 = n `mod` 8+        n8_17 = n8 == 1 || n8 == 7+        n4_1 = n8 == 1 || n8 == 5+        s4_1 = s `mod` 4 == 1+        g = if even r || n8_17 then f else not f+        h = if n4_1 || s4_1 then g else not g++-- The first argument (the modulus) must be an odd natural greater than 1+isResidue :: Natural -> Natural -> Bool+isResidue n m =+    case jacobiSymbol n m of+      Residue -> True+      _ -> False++-- The first argument (the modulus) must be an odd natural greater than 1+isNonResidue :: Natural -> Natural -> Bool+isNonResidue n m =+    case jacobiSymbol n m of+      NonResidue -> True+      _ -> False++-- The first argument (the modulus) must be an odd natural greater than 1+nextResidue :: Natural -> Natural -> Natural+nextResidue n =+    loop+  where+    loop m = if isResidue n m then m else loop (m + 1)++-- The first argument (the modulus) must be an odd natural greater than 1+nextNonResidue :: Natural -> Natural -> Natural+nextNonResidue n =+    loop+  where+    loop m = if isNonResidue n m then m else loop (m + 1)++-- The first argument (the modulus) must be a prime congruent to 3 mod 4+-- The second argument must be a residue modulo the prime+rootModuloPrime3Mod4 :: Natural -> Natural -> Natural+rootModuloPrime3Mod4 p =+    \n -> Modular.exp p n k+  where+    k = (p + 1) `div` 4++-- The first argument (the modulus) must be a prime congruent to 5 mod 8+-- The second argument must be a residue modulo the prime+rootModuloPrime5Mod8 :: Natural -> Natural -> Natural+rootModuloPrime5Mod8 p =+    go+  where+    go n =+        Modular.multiply p n (Modular.multiply p v (i - 1))+      where+        m = Modular.double p n+        v = Modular.exp p m k+        i = Modular.multiply p m (Modular.square p v)++    k = (p - 5) `div` 8++-- The Tonelli-Shanks algorithm+-- The first argument (the modulus) must be a prime+-- The second argument must be a residue modulo the prime+rootModuloPrime :: Natural -> Natural -> Natural+rootModuloPrime p =+    if p == 2 then Modular.normalize p+    else if r == 1 then rootModuloPrime3Mod4 p+    else if r == 2 then rootModuloPrime5Mod8 p+    else tonelliShanks+  where+    (r,s) = factorTwos (p - 1)+    z = Modular.exp p (nextNonResidue p 2) s++    tonelliShanks n =+        tonelliShanksLoop z d t r+      where+        d = Modular.exp p n ((s + 1) `div` 2)+        t = Modular.exp p n s++    tonelliShanksLoop c d t m =+        if t == 1 then d else tonelliShanksLoop b2 db tb2 i+      where+        i = tonelliShanksMin t 1+        b = Modular.exp2 p c (m - (i + 1))+        b2 = Modular.square p b+        db = Modular.multiply p d b+        tb2 = Modular.multiply p t b2++    tonelliShanksMin t i =+        if t2 == 1 then i else tonelliShanksMin t2 (i + 1)+      where+        t2 = Modular.square p t

@@ -10,46 +10,58 @@ 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+randomPairWith :: (a -> b -> c) -> (Random.Random -> a) ->+                  (Random.Random -> b) -> Random.Random -> c+randomPairWith f ra rb r =+    f (ra r1) (rb r2)   where-    n = shiftL 1 (w - 1)+    (r1,r2) = Random.split r -randomOdd :: Int -> Random.Random -> Natural-randomOdd w r = Bits.cons True (randomWidth (w - 1) r)+randomPair ::+    (Random.Random -> a) -> (Random.Random -> b) -> Random.Random -> (a,b)+randomPair = randomPairWith (,) -randomCoprime :: Int -> Random.Random -> (Natural,Natural)-randomCoprime w =+randomMaybe :: (Random.Random -> Maybe a) -> Random.Random -> a+randomMaybe g =     loop   where     loop r =-        case gen r1 of-          Just ab -> ab+        case g r1 of+          Just a -> a           Nothing -> loop r2       where         (r1,r2) = Random.split r -    gen r =-        if g == 1 then Just (a,b) else Nothing+randomFilter :: (a -> Bool) -> (Random.Random -> a) -> Random.Random -> a+randomFilter p g =+    randomMaybe gp+  where+    gp r =+        if p x then Just x else Nothing       where-        a = randomWidth w r1-        b = randomWidth w r2-        (g,_) = egcd a b-        (r1,r2) = Random.split r+        x = g r -uniformInteger :: Integer -> Random.Random -> Integer-uniformInteger n r = fromIntegral (Uniform.random (fromIntegral n) r)+randomWidth :: Natural -> Random.Random -> Natural+randomWidth w r =+    n + Uniform.random n r+  where+    n = shiftLeft 1 (w - 1) -randomCoprimeInteger :: Int -> Random.Random -> (Integer,Integer)-randomCoprimeInteger w r =-    (fromIntegral a, fromIntegral b)+randomOdd :: Natural -> Random.Random -> Natural+randomOdd w r = Bits.cons True (randomWidth (w - 1) r)++randomCoprime :: Natural -> Random.Random -> (Natural,Natural)+randomCoprime w =+    randomMaybe gen   where-    (a,b) = randomCoprime w r+    gen r =+        if g == 1 then Just (a,b) else Nothing+      where+        (a,b) = randomPair (randomWidth w) (randomWidth w) r+        (g,_) = egcd a b

@@ -1,86 +0,0 @@-{- |-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

@@ -1,72 +0,0 @@-{- |-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))

@@ -0,0 +1,46 @@+{- |+module: Arithmetic.Utility+description: Utility functions+license: MIT++maintainer: Joe Leslie-Hurd <joe@gilith.com>+stability: provisional+portability: portable+-}+module Arithmetic.Utility+where++import OpenTheory.Primitive.Natural+import OpenTheory.Natural.Divides+import qualified OpenTheory.Natural.Bits as Bits++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'++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++factorTwos :: Natural -> (Natural,Natural)+factorTwos n =+   if Bits.headBits n then (0,n) else (r + 1, s)+  where+    (r,s) = factorTwos (Bits.tailBits n)++factorOut :: Natural -> Natural -> (Natural,Natural)+factorOut p =+    go 0+  where+    go k n = if divides p n then go (k + 1) (n `div` p) else (k,n)

@@ -0,0 +1,84 @@+{- |+module: Arithmetic.Utility.Heap+description: Leftist heaps+license: MIT++maintainer: Joe Leslie-Hurd <joe@gilith.com>+stability: provisional+portability: portable+-}+module Arithmetic.Utility.Heap+  ( Heap,+    size,+    isEmpty,+    empty,+    add,+    remove,+    toList )+where++data Node a =+    E+  | T Int a (Node a) (Node a)+  deriving Show++data Heap a =+    Heap (a -> a -> Bool) Int (Node a)++singleton :: a -> Node a+singleton a = T 1 a E E++rank :: Node a -> Int+rank E = 0+rank (T r _ _ _) = r++mkT :: a -> Node a -> Node a -> Node a+mkT a x y =+    if rx <= ry+      then T (rx + 1) a y x+      else T (ry + 1) a x y+  where+    rx = rank x+    ry = rank y++merge :: (a -> a -> Bool) -> Node a -> Node a -> Node a+merge le =+    mrg+  where+    mrg n1 n2 =+      case n1 of+        E -> n2+        T _ a1 x1 y1 ->+          case n2 of+            E -> n1+            T _ a2 x2 y2 ->+              if le a1 a2+                then mkT a1 x1 (mrg y1 n2)+                else mkT a2 x2 (mrg n1 y2)++size :: Heap a -> Int+size (Heap _ k _) = k++isEmpty :: Heap a -> Bool+isEmpty h = size h == 0++empty :: (a -> a -> Bool) -> Heap a+empty le = Heap le 0 E++add :: a -> Heap a -> Heap a+add a (Heap le k n) = Heap le (k + 1) (merge le (singleton a) n)++remove :: Heap a -> Maybe (a, Heap a)+remove (Heap le k n) =+    case n of+      E -> Nothing+      T _ a x y -> Just (a, Heap le (k - 1) (merge le x y))++toList :: Heap a -> [a]+toList h =+    case remove h of+      Nothing -> []+      Just (a,h') -> a : toList h'++instance Show a => Show (Heap a) where+  show = show . toList

@@ -0,0 +1,117 @@+{- |+module: Arithmetic.Williams+description: Williams p+1 factorization method+license: MIT++maintainer: Joe Leslie-Hurd <joe@gilith.com>+stability: provisional+portability: portable+-}+module Arithmetic.Williams+where++--import Debug.Trace(trace)+import OpenTheory.Primitive.Natural+import qualified OpenTheory.Natural.Bits as Bits+import qualified OpenTheory.Primitive.Random as Random+import qualified OpenTheory.Natural.Uniform as Uniform++import Arithmetic.Prime+import Arithmetic.Utility+import qualified Arithmetic.Lucas as Lucas+import qualified Arithmetic.Modular as Modular++sequence :: a -> a -> (a -> a -> a) -> (a -> a -> a) -> a -> [a]+sequence one two sub mult p = Lucas.vSequence two sub mult p one++nthExp :: a -> (a -> a -> a) -> (a -> a -> a) -> a -> Natural -> Natural -> a+nthExp two sub mult p n k =+    if k == 0 then p+    else if n == 0 then two+    else functionPower nthSeq k p+  where+    l = init (Bits.toList n)+    sq z = sub (mult z z) two+    nthSeq v =+        w+      where+        (w,_) = foldr inc (v, sq v) l+        inc b (x,y) =+           if b then (z, sq y) else (sq x, z)+         where+           z = sub (mult x y) v++nth :: a -> (a -> a -> a) -> (a -> a -> a) -> a -> Natural -> a+nth two sub mult p n = nthExp two sub mult p n 1++base :: Natural -> Natural -> Random.Random -> Either Natural [Natural]+base n =+    go+  where+    go x rnd =+        if x == 0 then Right []+        else mcons (gen r1) (go (x - 1) r2)+      where+        (r1,r2) = Random.split rnd++    mcons (Right v) (Right vs) = Right (v : vs)+    mcons _ vs = vs++    gen rnd =+        if 1 < g then Left g else Right v+      where+        v = Uniform.random (n - 3) rnd + 2+        g = gcd n v++method :: Natural -> [Natural] -> [Natural] -> Maybe Natural+method n =+    loop+  where+    w = Bits.width n++    loop [] _ = Nothing+    loop _ [] = Nothing+    loop vs (p : ps) =+        case fltr vs p k of+          Left g -> Just g+          Right vs' -> loop vs' ps+      where+        -- log_p n = log_2 n / log_2 p <= |n| / (|p| - 1)+        k = w `div` (Bits.width p - 1)++    fltr [] _ _ = Right []+    fltr (v : vs) p k =+        case check v p k of+          Left g -> Left g+          Right v' -> mcons v' (fltr vs p k)++    mcons (Just v) (Right vs) = Right (v : vs)+    mcons _ vs = vs++    check v p k =+        if g == n then Right Nothing+        else if 1 < g then+          --trace ("Williams p+1 method succeeded with prime " ++ show p) $+          Left g+        else Right (Just (pow v p k))+      where+        g = gcd n (v - 2)++    pow =+        nthExp two sub mult+      where+        two = Modular.normalize n 2+        sub = Modular.subtract n+        mult = Modular.multiply n++-- Works for odd numbers at least 5+factor :: Natural -> Maybe Natural ->+          Natural -> Random.Random -> Maybe Natural+factor x k n rnd =+    case base n x rnd of+      Left g -> Just g+      Right vs -> method n vs ps+  where+    ps = case k of+           Just m -> take (fromIntegral m) primes+           Nothing -> primes

@@ -20,8 +20,10 @@ import qualified OpenTheory.Natural.Uniform as Uniform  import Arithmetic.Random+import qualified Arithmetic.Prime.Factor as Factor import qualified Arithmetic.Modular as Modular import qualified Arithmetic.Montgomery as Montgomery+import qualified Arithmetic.Williams as Williams  -------------------------------------------------------------------------------- -- Helper functions@@ -42,22 +44,33 @@ --------------------------------------------------------------------------------  data Operation =-    Modexp+    Factor+  | Modexp   | Timelock   deriving Show  operations :: [Operation]-operations = [Modexp,Timelock]+operations = [Factor,Modexp,Timelock]  operationToString :: Operation -> String operationToString oper =    case oper of+     Factor -> "factor"      Modexp -> "modexp"      Timelock -> "timelock" +operationsToString :: [Operation] -> String+operationsToString = setToString operationToString+ stringToOperation :: String -> Operation stringToOperation = getPrefixString "operation" operationToString operations +getOperation :: [String] -> (Operation,[String])+getOperation args =+    case args of+      [] -> usage "no operation specified"+      h : t -> (stringToOperation h, t)+ -------------------------------------------------------------------------------- -- Algorithms --------------------------------------------------------------------------------@@ -65,17 +78,30 @@ data Algorithm =     Modular   | Montgomery+  | Williams   deriving Show  algorithms :: [Algorithm]-algorithms = [Modular,Montgomery]+algorithms = [Modular,Montgomery,Williams] +possibleAlgorithms :: Operation -> [Algorithm]+possibleAlgorithms Factor = [Williams]+possibleAlgorithms Modexp = [Modular,Montgomery]+possibleAlgorithms Timelock = [Modular,Montgomery]++defaultAlgorithm :: Operation -> Algorithm+defaultAlgorithm = last . possibleAlgorithms+ algorithmToString :: Algorithm -> String algorithmToString oper =    case oper of      Modular -> "modular"      Montgomery -> "montgomery"+     Williams -> "williams" +algorithmsToString :: [Algorithm] -> String+algorithmsToString = setToString algorithmToString+ stringToAlgorithm :: String -> Algorithm stringToAlgorithm = getPrefixString "algorithm" algorithmToString algorithms @@ -85,7 +111,7 @@  data InputNatural =     Fixed Natural-  | Width Int+  | Width Natural   deriving Show  stringToInputNatural :: String -> InputNatural@@ -98,34 +124,24 @@             [(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+widthInputNatural :: InputNatural -> Random.Random -> Natural+widthInputNatural (Fixed n) _ = n+widthInputNatural (Width w) r = randomWidth 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+rsaInputNatural :: InputNatural -> Random.Random -> Natural+rsaInputNatural (Fixed n) _ = n+rsaInputNatural (Width w) rnd = Factor.toNatural (Factor.randomRSA w rnd) -    (rn,r') = Random.split r-    (rx,rk) = Random.split r'+getInput :: Operation -> String -> Maybe InputNatural -> InputNatural+getInput oper s m =+    case m of+      Just n -> n+      Nothing -> usage $ "specify " ++ s ++ " parameter for " +++                         operationToString oper ++ " operation"  -------------------------------------------------------------------------------- -- Options@@ -133,100 +149,179 @@  data Options = Options     {optOperation :: Operation,-     optAlgorithm :: Algorithm,-     optModulus :: InputNatural,-     optBase :: Maybe InputNatural,-     optExponent :: Maybe InputNatural}+     optA :: Algorithm,+     optN :: Maybe InputNatural,+     optX :: Maybe InputNatural,+     optK :: Maybe InputNatural}   deriving Show -defaultOptions :: Options-defaultOptions =+nullOptions :: Options+nullOptions =   Options-    {optOperation = Modexp,-     optAlgorithm = Montgomery,-     optModulus = Width 50,-     optBase = Nothing,-     optExponent = Nothing}+    {optOperation = Factor,+     optA = Williams,+     optN = Nothing,+     optX = Nothing,+     optK = 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")+    [Option ['a'] []+       (algorithmArg (\alg opts -> opts {optA = alg}))        "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"]+     Option ['n'] []+       (inputNaturalArg (\n opts -> opts {optN = n}))+       "select n parameter",+     Option ['x'] []+       (inputNaturalArg (\x opts -> opts {optX = x}))+       "select x parameter",+     Option ['k'] []+       (inputNaturalArg (\k opts -> opts {optK = k}))+       "select k parameter"]+  where+    algorithmArg f = ReqArg (\s -> f (stringToAlgorithm s)) "ALGORITHM"+    inputNaturalArg f =+        ReqArg (\s -> f (Just (stringToInputNatural s))) "NATURAL" -processOptions :: [String] -> Either [String] (Options,[String])-processOptions args =+processOptions :: Options -> [String] -> Either [String] (Options,[String])+processOptions opts args =     case getOpt Permute options args of-      (opts,work,[]) -> Right (foldl (flip id) defaultOptions opts, work)+      (opts',args',[]) -> Right (foldl (flip id) opts opts', args')       (_,_,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"+processOperation :: Options -> Operation -> Options+processOperation opts oper =+    opts {optOperation = oper, optA = defaultAlgorithm oper}  usage :: String -> a usage err =     error $ err ++ "\n" ++ usageInfo header options ++ footer   where-    header = "Usage: modexp [OPTION...]"+    header = "Usage: arithmetic OPERATION [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]."+      "where OPERATION is one of " ++ operationsToString operations ++ ",\n" +++      "  ( factor.........factorize n                       )\n" +++      "  ( modexp.........compute (x ^ k) `mod` n           )\n" +++      "  ( timelock.......compute (x ^ 2 ^ k) `mod` n       )\n" +++      "ALGORITHM is one of " ++ algorithmsToString algorithms ++ ",\n" +++      "  ( modular........naive modular arithmetic          )\n" +++      "  ( montgomery.....Montgomery multiplication         )\n" +++      "  ( williams.......Williams p+1 factorization method )\n" +++      "and NATURAL is either a natural number or has the form [bitwidth]." +usageOperation :: Operation -> a+usageOperation oper =+    error $ err ++ "\n" ++ usageInfo header options ++ footer+  where+    err = "bad algorithm"++    algs = possibleAlgorithms oper++    header = "Usage: arithmetic " ++ operationToString oper ++ " [OPTION...]"++    footer =+      "where ALGORITHM is one of " ++ algorithmsToString algs ++ ",\n" +++      "and NATURAL is either a natural number or has the form [bitwidth]."+ -------------------------------------------------------------------------------- -- Computation -------------------------------------------------------------------------------- -type Computation = Natural -> Natural -> Natural -> Natural+computeFactorWilliams :: Options ->+                         Natural -> Random.Random -> Maybe Factor.Factor+computeFactorWilliams opts n rnd =+    Factor.factor 1000 (Williams.factor x k) n r3+  where+    x = case optX opts of+          Nothing -> 5+          Just w -> widthInputNatural w r1+    k = case optK opts of+          Nothing -> Nothing+          Just w -> Just (widthInputNatural w r2)+    (r1,r23) = Random.split rnd+    (r2,r3) = Random.split r23 -computation :: Operation -> Algorithm -> Computation-computation Modexp Modular = Modular.exp-computation Modexp Montgomery = Montgomery.modexp-computation Timelock Modular = Modular.exp2-computation Timelock Montgomery = Montgomery.modexp2+computeFactor :: Operation -> Options -> Random.Random -> String+computeFactor oper opts rnd =+    case x of+      Nothing -> error $ "factorization failed for " ++ show n+      Just f -> show n ++ (if Factor.isPrime f then " is prime"+                           else " == " ++ show f)+  where+    n = rsaInputNatural (getInput oper "n" (optN opts)) r1+    x = case optA opts of+          Williams -> computeFactorWilliams opts n r2+          _ -> usageOperation oper+    (r1,r2) = Random.split rnd -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+computeModexp :: Operation -> Options -> Random.Random -> String+computeModexp oper opts rnd =+    "( " ++ show x ++ " ^ " ++ show k ++ " ) `mod` " ++ show n +++    " == " ++ show y+  where+    n = oddInputNatural (getInput oper "n" (optN opts)) r1+    x = case optX opts of+          Nothing -> Uniform.random n r2+          Just w -> widthInputNatural w r2+    k = case optK opts of+          Nothing -> Uniform.random n r3+          Just w -> widthInputNatural w r3+    f = case optA opts of+          Modular -> Modular.exp+          Montgomery -> Montgomery.modexp+          _ -> usageOperation oper+    y = f n x k+    (r1,r23) = Random.split rnd+    (r2,r3) = Random.split r23 +computeTimelock :: Operation -> Options -> Random.Random -> String+computeTimelock oper opts rnd =+    "( " ++ show x ++ " ^ 2 ^ " ++ show k ++ " ) `mod` " ++ show n +++     " == " ++ show y+  where+    n = oddInputNatural (getInput oper "n" (optN opts)) r1+    x = case optX opts of+          Nothing -> Uniform.random n r2+          Just w -> widthInputNatural w r2+    k = widthInputNatural (getInput oper "k" (optK opts)) r3+    f = case optA opts of+          Modular -> Modular.exp2+          Montgomery -> Montgomery.modexp2+          _ -> usageOperation oper+    y = f n x k+    (r1,r23) = Random.split rnd+    (r2,r3) = Random.split r23++compute :: Options -> Random.Random -> String+compute opts =+    case oper of+      Factor -> computeFactor oper opts+      Modexp -> computeModexp oper opts+      Timelock -> computeTimelock oper opts+  where+    oper = optOperation opts+ -------------------------------------------------------------------------------- -- Main program -------------------------------------------------------------------------------- +processArguments :: [String] -> Options+processArguments cmd =+    case processOptions opts args of+      Left errs -> usage (concat errs)+      Right (opts',work) ->+        case work of+          [] -> opts'+          _ : _ -> usage "too many arguments"+  where+    (oper,args) = getOperation cmd+    opts = processOperation nullOptions oper+ main :: IO () main =     do args <- Environment.getArgs-       r <- fmap Random.fromInt System.Random.randomIO+       rnd <- 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+       putStrLn $ compute opts rnd        return ()

@@ -1,6 +1,6 @@ {- | module: Main-description: Testing the modular exponentiation computation+description: Testing the natural number arithmetic library license: MIT  maintainer: Joe Leslie-Hurd <joe@gilith.com>@@ -15,102 +15,81 @@ import OpenTheory.Primitive.Natural import OpenTheory.Natural import OpenTheory.Natural.Divides+import qualified OpenTheory.Natural.Bits as Bits+import qualified OpenTheory.Natural.Prime as Prime import qualified OpenTheory.Primitive.Random as Random import qualified OpenTheory.Natural.Uniform as Uniform-import OpenTheory.Primitive.Test  import Arithmetic.Random import Arithmetic.Prime import qualified Arithmetic.ContinuedFraction as ContinuedFraction+import qualified Arithmetic.Prime.Factor as Factor import qualified Arithmetic.Modular as Modular import qualified Arithmetic.Montgomery as Montgomery-import qualified Arithmetic.Smooth as Smooth-import qualified Arithmetic.SquareRoot as SquareRoot+import qualified Arithmetic.Quadratic as Quadratic+import qualified Arithmetic.Williams as Williams -propEgcdDivides :: Natural -> Natural -> Bool-propEgcdDivides a b =-    divides g a && divides g b-  where-    (g,_) = egcd a b+propPrimes :: Natural -> Bool+propPrimes k =+    primes !! (fromIntegral k) ==+    Prime.primes !! (fromIntegral k) -propEgcdEquation :: Natural -> Natural -> Bool-propEgcdEquation ap b =-    s * a == t * b + g+propRandomPrime :: Natural -> Random.Random -> Bool+propRandomPrime wp rnd =+    Bits.width p == w &&+    isPrime p r2   where-    a = ap + 1-    (g,(s,t)) = egcd a b+    w = wp + 2+    p = randomPrime w r1+    (r1,r2) = Random.split rnd -propEgcdBound :: Natural -> Natural -> Bool-propEgcdBound ap b =-    s < max b 2 && t < a+propRandomRSA :: Natural -> Random.Random -> Bool+propRandomRSA wp rnd =+    Bits.width n == w &&+    not (isPrime n r2)   where-    a = ap + 1-    (_,(s,t)) = egcd a b+    w = wp + 5+    n = Factor.toNatural (Factor.randomRSA w r1)+    (r1,r2) = Random.split rnd -propSmoothInjective :: Natural -> Natural -> Bool-propSmoothInjective k np =-    Smooth.toNatural (Smooth.fromNatural k n) == n+propTrialDivision :: Natural -> Natural -> Bool+propTrialDivision k np =+    Factor.toNatural f * m == n &&+    all (\p -> not (divides p m)) ps   where     n = np + 1--propFloorSqrt :: Natural -> Bool-propFloorSqrt n =-    sq s <= n && n < sq (s + 1)-  where-    s = SquareRoot.floor n-    sq i = i * i--propCeilingSqrt :: Natural -> Bool-propCeilingSqrt n =-    (s == 0 || sq (s - 1) < n) && n <= sq s-  where-    s = SquareRoot.ceiling n-    sq i = i * i--propContinuedFractionSqrt :: Natural -> Bool-propContinuedFractionSqrt n =-    cf == spec-  where-    cf = ContinuedFraction.toDouble (SquareRoot.continuedFraction n)-    spec = sqrt (fromIntegral n)--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 = chineseRemainder a b x y-    (r1,r23) = Random.split r-    (r2,r3) = Random.split r23+    ps = take (fromIntegral k) primes+    (f,m) = Factor.trialDivision ps n -propModularNegate :: Int -> Random.Random -> Bool-propModularNegate nw rnd =+propModularNegate :: Natural -> Random.Random -> Bool+propModularNegate np rnd =     Modular.add n a b == 0 &&     b < n   where-    n = randomWidth nw r1-    a = Uniform.random n r2+    n = np + 1+    a = Uniform.random n rnd     b = Modular.negate n a-    (r1,r2) = Random.split rnd -propModularInvert :: Int -> Random.Random -> Bool-propModularInvert nw rnd =+propModularInvert :: Natural -> Natural -> Bool+propModularInvert np a =     case Modular.invert n a of-      Nothing -> True-      Just b -> Modular.multiply n a b == 1 && b < n+      Nothing -> gcd n a /= 1+      Just b -> Modular.multiply n a b == Modular.normalize n 1 && b < n   where-    n = randomWidth nw r1-    a = Uniform.random n r2-    (r1,r2) = Random.split rnd+    n = np + 1 -randomMontgomeryParameters :: Int -> Random.Random -> Montgomery.Parameters-randomMontgomeryParameters w r = Montgomery.standardParameters (randomOdd w r)+propFermat :: Natural -> Random.Random -> Bool+propFermat pp rnd =+    Modular.exp p a p == a+  where+    p = nextPrime (pp + 3) r1+    a = Uniform.random p r2+    (r1,r2) = Random.split rnd -propMontgomeryInvariant :: Int -> Random.Random -> Bool-propMontgomeryInvariant nw rnd =+propMontgomeryInvariant :: Natural -> Natural -> Bool+propMontgomeryInvariant np b =     naturalOdd n &&+    1 < n &&     n < w2 &&     s * w2 == k * n + 1 &&     s < n &&@@ -128,214 +107,296 @@        Montgomery.kParameters = k,        Montgomery.rParameters = r,        Montgomery.r2Parameters = r2,-       Montgomery.zParameters = z} = randomMontgomeryParameters nw rnd-+       Montgomery.zParameters = z} =+      Montgomery.customParameters (2 * np + 3) (Bits.width n + b)     w2 = shiftLeft 1 w -propMontgomeryNormalize :: Int -> Random.Random -> Bool-propMontgomeryNormalize nw rnd =+propMontgomeryNormalize :: Natural -> Natural -> Bool+propMontgomeryNormalize np a =     b `mod` n == a `mod` n &&     b < w2   where-    p = randomMontgomeryParameters nw r1-    a = Uniform.random (w2 * w2) r2+    n = 2 * np + 3+    p = Montgomery.standardParameters n     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 =+propMontgomeryReduce :: Natural -> Natural -> Bool+propMontgomeryReduce np a =     b `mod` n == (a * s) `mod` n &&     b < w2 + n   where-    p = randomMontgomeryParameters nw r1-    a = Uniform.random (w2 * w2) r2+    n = 2 * np + 3+    p = Montgomery.standardParameters n     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 =+propMontgomeryReduceSmall :: Natural -> Natural -> Bool+propMontgomeryReduceSmall np ap =     b `mod` n == (a * s) `mod` n &&     b <= n   where-    p = randomMontgomeryParameters nw r1-    a = Uniform.random w2 r2+    n = 2 * np + 3+    p = Montgomery.standardParameters n+    a = ap `mod` w2     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 =+propMontgomeryToNatural :: Natural -> Natural -> Bool+propMontgomeryToNatural np a =     b == (a * s) `mod` n   where-    p = randomMontgomeryParameters nw r1-    a = Uniform.random w2 r2+    n = 2 * np + 3+    p = Montgomery.standardParameters n     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 =+propMontgomeryFromNatural :: Natural -> Natural -> Bool+propMontgomeryFromNatural np a =     b == a `mod` n   where-    p = randomMontgomeryParameters nw r1-    a = Uniform.random (w2 * w2) r2+    n = 2 * np + 3+    p = Montgomery.standardParameters n     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 =+propMontgomeryZero :: Natural -> Bool+propMontgomeryZero np =     Montgomery.toNatural (Montgomery.zero p) == 0   where-    p = randomMontgomeryParameters nw rnd+    n = 2 * np + 3+    p = Montgomery.standardParameters n -propMontgomeryOne :: Int -> Random.Random -> Bool-propMontgomeryOne nw rnd =+propMontgomeryOne :: Natural -> Bool+propMontgomeryOne np =     Montgomery.toNatural (Montgomery.one p) == 1   where-    p = randomMontgomeryParameters nw rnd+    n = 2 * np + 3+    p = Montgomery.standardParameters n -propMontgomeryTwo :: Int -> Random.Random -> Bool-propMontgomeryTwo nw rnd =+propMontgomeryTwo :: Natural -> Bool+propMontgomeryTwo np =     Montgomery.toNatural (Montgomery.two p) == 2   where-    p = randomMontgomeryParameters nw rnd+    n = 2 * np + 3+    p = Montgomery.standardParameters n -propMontgomeryAdd :: Int -> Random.Random -> Bool-propMontgomeryAdd nw rnd =+propMontgomeryAdd :: Natural -> Natural -> Natural -> Bool+propMontgomeryAdd np ap bp =     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)+    n = 2 * np + 3+    p = Montgomery.standardParameters n+    a = Montgomery.normalize p ap+    b = Montgomery.normalize p bp     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 =+propMontgomeryNegate :: Natural -> Natural -> Bool+propMontgomeryNegate np ap =     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)+    n = 2 * np + 3+    p = Montgomery.standardParameters n+    a = Montgomery.normalize p ap     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 =+propMontgomeryMultiply :: Natural -> Natural -> Natural -> Bool+propMontgomeryMultiply np ap bp =     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)+    n = 2 * np + 3+    p = Montgomery.standardParameters n+    a = Montgomery.normalize p ap+    b = Montgomery.normalize p bp     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 =+propMontgomeryModexp :: Natural -> Natural -> Natural -> Bool+propMontgomeryModexp np x k =     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+    n = 2 * np + 3 -propMontgomeryModexp2 :: Int -> Random.Random -> Bool-propMontgomeryModexp2 w r =+propMontgomeryModexp2 :: Natural -> Natural -> Natural -> Bool+propMontgomeryModexp2 np x k =     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+    n = 2 * np + 3 -    (r1,r23) = Random.split r-    (r2,r3) = Random.split r23+propRootFloor :: Natural -> Bool+propRootFloor n =+    sq s <= n && n < sq (s + 1)+  where+    s = Quadratic.rootFloor n+    sq i = i * i -propFermat :: Int -> Random.Random -> Bool-propFermat w r =-    Montgomery.modexp n a n == a+propRootCeiling :: Natural -> Bool+propRootCeiling n =+    (s == 0 || sq (s - 1) < n) && n <= sq s   where-    n = randomPrime w r1-    a = Uniform.random n r2-    (r1,r2) = Random.split r+    s = Quadratic.rootCeiling n+    sq i = i * i -checkWidthProp ::-    QuickCheck.Testable prop => Int -> String -> (Int -> prop) -> IO ()-checkWidthProp w s p =-    check (s ++ " (" ++ show w ++ " bit)\n  ") (p w)+propRootContinuedFraction :: Natural -> Bool+propRootContinuedFraction n =+    cf == spec+  where+    cf = ContinuedFraction.toDouble (Quadratic.rootContinuedFraction n)+    spec = sqrt (fromIntegral n) -checkWidthProps :: Int -> IO ()-checkWidthProps w =-   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 ()+propJacobiSymbol :: Natural -> Natural -> Random.Random -> Bool+propJacobiSymbol np m rnd =+    case Quadratic.jacobiSymbol n m of+      Quadratic.Zero -> not coprime+      Quadratic.Residue -> coprime && (mr || not (isPrime n rnd))+      Quadratic.NonResidue -> coprime && not mr+  where+    coprime = gcd m n == 1+    n = 2 * np + 1+    mn = Modular.normalize n m+    mr = any (\k -> Modular.square n k == mn) [1..np] +propRootModuloPrime3Mod4 :: Natural -> Random.Random -> Bool+propRootModuloPrime3Mod4 pp rnd =+    Modular.square p r == a+  where+    p = nextPrime3Mod4 pp r1+    a = randomFilter (Quadratic.isResidue p) (Uniform.random p) r2+    r = Quadratic.rootModuloPrime3Mod4 p a+    (r1,r2) = Random.split rnd++propRootModuloPrime5Mod8 :: Natural -> Random.Random -> Bool+propRootModuloPrime5Mod8 pp rnd =+    Modular.square p r == a+  where+    p = nextPrime5Mod8 pp r1+    a = randomFilter (Quadratic.isResidue p) (Uniform.random p) r2+    r = Quadratic.rootModuloPrime5Mod8 p a+    (r1,r2) = Random.split rnd++propRootModuloPrime :: Natural -> Random.Random -> Bool+propRootModuloPrime pp rnd =+    Modular.square p r == a+  where+    p = nextPrime pp r1+    a = randomFilter (Quadratic.isResidue p) (Uniform.random p) r2+    r = Quadratic.rootModuloPrime p a+    (r1,r2) = Random.split rnd++propWilliamsNth :: Natural -> Natural -> Natural -> Bool+propWilliamsNth np p k =+    Williams.sequence one two sub mult p !! (fromIntegral k) ==+    Williams.nth two sub mult p k+  where+    n = np + 1+    one = 1+    two = 2+    sub = Modular.subtract n+    mult = Modular.multiply n++propWilliamsNthProduct :: Natural -> Natural -> Natural -> Natural -> Bool+propWilliamsNthProduct np pp i j =+    Williams.nth two sub mult p (i * j) ==+    Williams.nth two sub mult (Williams.nth two sub mult p i) j+  where+    n = np + 1+    p = pp + 1+    two = Modular.normalize n 2+    sub = Modular.subtract n+    mult = Modular.multiply n++propWilliamsNthExp :: Natural -> Natural -> Natural -> Natural -> Bool+propWilliamsNthExp np p m k =+    Williams.nthExp two sub mult p m k ==+    Williams.nth two sub mult p (m ^ k)+  where+    n = np + 1+    two = Modular.normalize n 2+    sub = Modular.subtract n+    mult = Modular.multiply n++propWilliamsNthEqTwo :: Natural -> Natural -> Natural -> Random.Random -> Bool+propWilliamsNthEqTwo pp a mp rnd =+    Williams.nth two sub mult a m == two+  where+    p = nextPrime (pp + 3) rnd+    d = sub (mult a a) 4+    t = case Quadratic.jacobiSymbol p d of+          Quadratic.Zero -> 2+          Quadratic.Residue -> p - 1+          Quadratic.NonResidue -> p + 1+    m = mp * t+    two = Modular.normalize p 2+    sub = Modular.subtract p+    mult = Modular.multiply p++propWilliamsFactor :: Natural -> Natural -> Natural -> Random.Random -> Bool+propWilliamsFactor np x k rnd =+    case Williams.factor x (Just k) n rnd of+      Nothing -> True+      Just p -> 1 < p && p < n && divides p n+  where+    n = 2 * np + 5++check :: QuickCheck.Testable prop => String -> prop -> IO ()+check desc prop =+    do putStr (desc ++ "\n  ")+       res <- QuickCheck.quickCheckWithResult args prop+       case res of+         QuickCheck.Failure {} -> error "Proposition failed"+         _ -> return ()+  where+    args = QuickCheck.stdArgs {QuickCheck.maxSuccess = 1000}+ main :: IO () main =-    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+    do check "Sieve of Eratosphenes" propPrimes+       check "Generating random primes" propRandomPrime+       check "Generating random RSA moduli" propRandomRSA+       check "Trial division" propTrialDivision+       check "Modular negate" propModularNegate+       check "Modular invert" propModularInvert+       check "Fermat's little theorem" propFermat+       check "Montgomery invariant" propMontgomeryInvariant+       check "Montgomery normalize" propMontgomeryNormalize+       check "Montgomery reduce" propMontgomeryReduce+       check "Montgomery reduce small" propMontgomeryReduceSmall+       check "Montgomery toNatural" propMontgomeryToNatural+       check "Montgomery fromNatural" propMontgomeryFromNatural+       check "Montgomery zero" propMontgomeryZero+       check "Montgomery one" propMontgomeryOne+       check "Montgomery two" propMontgomeryTwo+       check "Montgomery add" propMontgomeryAdd+       check "Montgomery negate" propMontgomeryNegate+       check "Montgomery multiply" propMontgomeryMultiply+       check "Montgomery modexp" propMontgomeryModexp+       check "Montgomery modexp2" propMontgomeryModexp2+       check "Floor square root" propRootFloor+       check "Ceiling square root" propRootCeiling+       check "Continued fraction square root" propRootContinuedFraction+       check "Jacobi symbol" propJacobiSymbol+       check "Square root modulo prime congruent to 3 mod 4"+         propRootModuloPrime3Mod4+       check "Square root modulo prime congruent to 5 mod 8"+         propRootModuloPrime5Mod8+       check "Square root modulo prime" propRootModuloPrime+       check "Williams sequence" propWilliamsNth+       check "Williams sequence product" propWilliamsNthProduct+       check "Williams sequence exponential" propWilliamsNthExp+       check "Williams sequence equals two" propWilliamsNthEqTwo+       check "Williams factorization works" propWilliamsFactor        return ()-  where-    ws = takeWhile (\n -> n <= 128) (iterate ((*) 2) 4)