diff --git a/arithmetic.cabal b/arithmetic.cabal --- a/arithmetic.cabal +++ b/arithmetic.cabal @@ -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 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, diff --git a/src/Arithmetic/Lucas.hs b/src/Arithmetic/Lucas.hs new file mode 100644 --- /dev/null +++ b/src/Arithmetic/Lucas.hs @@ -0,0 +1,28 @@ +{- | +module: Arithmetic.Lucas +description: Lucas sequences +license: MIT + +maintainer: Joe Leslie-Hurd +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 diff --git a/src/Arithmetic/Modular.hs b/src/Arithmetic/Modular.hs --- a/src/Arithmetic/Modular.hs +++ b/src/Arithmetic/Modular.hs @@ -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 diff --git a/src/Arithmetic/Montgomery.hs b/src/Arithmetic/Montgomery.hs --- a/src/Arithmetic/Montgomery.hs +++ b/src/Arithmetic/Montgomery.hs @@ -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 = diff --git a/src/Arithmetic/Prime.hs b/src/Arithmetic/Prime.hs --- a/src/Arithmetic/Prime.hs +++ b/src/Arithmetic/Prime.hs @@ -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 diff --git a/src/Arithmetic/Prime/Factor.hs b/src/Arithmetic/Prime/Factor.hs new file mode 100644 --- /dev/null +++ b/src/Arithmetic/Prime/Factor.hs @@ -0,0 +1,220 @@ +{- | +module: Arithmetic.Prime.Factor +description: Factorized natural numbers +license: MIT + +maintainer: Joe Leslie-Hurd +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 diff --git a/src/Arithmetic/Prime/Sieve.hs b/src/Arithmetic/Prime/Sieve.hs new file mode 100644 --- /dev/null +++ b/src/Arithmetic/Prime/Sieve.hs @@ -0,0 +1,50 @@ +{- | +module: Arithmetic.Prime.Sieve +description: The genuine sieve of Eratosphenes +license: MIT + +maintainer: Joe Leslie-Hurd +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 diff --git a/src/Arithmetic/Quadratic.hs b/src/Arithmetic/Quadratic.hs new file mode 100644 --- /dev/null +++ b/src/Arithmetic/Quadratic.hs @@ -0,0 +1,180 @@ +{- | +module: Arithmetic.Quadratic +description: Natural number square root +license: MIT + +maintainer: Joe Leslie-Hurd +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 diff --git a/src/Arithmetic/Random.hs b/src/Arithmetic/Random.hs --- a/src/Arithmetic/Random.hs +++ b/src/Arithmetic/Random.hs @@ -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 diff --git a/src/Arithmetic/Smooth.hs b/src/Arithmetic/Smooth.hs deleted file mode 100644 --- a/src/Arithmetic/Smooth.hs +++ /dev/null @@ -1,86 +0,0 @@ -{- | -module: Arithmetic.Smooth -description: Smooth numbers -license: MIT - -maintainer: Joe Leslie-Hurd -stability: provisional -portability: portable --} -module Arithmetic.Smooth -where - -import qualified Data.List as List -import OpenTheory.Primitive.Natural -import qualified OpenTheory.Natural.Bits as Bits -import OpenTheory.Natural.Divides -import qualified OpenTheory.Natural.Prime as Prime - -factorOut :: Natural -> Natural -> Maybe (Natural,Natural) -factorOut p = - go 0 - where - go k n = - if divides p n then go (k + 1) (n `div` p) - else if k == 0 then Nothing - else Just (k,n) - -factorList :: [Natural] -> Natural -> ([(Natural,Natural)],Natural) -factorList ps n = - case ps of - [] -> ([],n) - p : pt -> - case factorOut p n of - Nothing -> factorList pt n - Just (k,m) -> - let (pks,q) = factorList pt m in - ((p,k) : pks, q) - -factorBase :: Natural -> Natural -> ([(Natural,Natural)],Natural) -factorBase k = factorList (take (fromIntegral k) Prime.primes) - -multiplyBase :: ([(Natural,Natural)],Natural) -> Natural -multiplyBase = - \(pks,m) -> foldr mult m pks - where - mult (p,k) m = (p ^ k) * m - -newtype Smooth = - Smooth {unSmooth :: ([(Natural,Natural)],Natural)} - deriving (Eq,Ord) - -instance Show Smooth where - show s = - if null factors then "1" else List.intercalate "*" factors - where - factors = map showPk pks ++ showRest - showRest = if n == 1 then [] else [showWidth] - showWidth = if w < 20 then show n - else "[" ++ show w ++ "]" - showPk (p,k) = show p ++ showExp k - showExp k = if k == 1 then "" else "^" ++ show k - (pks,n) = unSmooth s - w = Bits.width n - -fromNatural :: Natural -> Natural -> Smooth -fromNatural k = Smooth . factorBase k - -toNatural :: Smooth -> Natural -toNatural = multiplyBase . unSmooth - -factoring :: Smooth -> Maybe [(Natural,Natural)] -factoring s = - if n == 1 then Just pks else Nothing - where - (pks,n) = unSmooth s - -next :: Natural -> Natural -> Smooth -next k = - go - where - go n = - case factoring s of - Nothing -> go (n + 1) - Just _ -> s - where - s = fromNatural k n diff --git a/src/Arithmetic/SquareRoot.hs b/src/Arithmetic/SquareRoot.hs deleted file mode 100644 --- a/src/Arithmetic/SquareRoot.hs +++ /dev/null @@ -1,72 +0,0 @@ -{- | -module: Arithmetic.SquareRoot -description: Natural number square root -license: MIT - -maintainer: Joe Leslie-Hurd -stability: provisional -portability: portable --} -module Arithmetic.SquareRoot -where - -import OpenTheory.Primitive.Natural -import qualified Data.List as List - -import qualified Arithmetic.ContinuedFraction as ContinuedFraction - -floor :: Natural -> Natural -floor n = - if n < 2 then n else bisect 0 n - where - bisect l u = - if m == l then l - else if m * m <= n then bisect m u - else bisect l m - where - m = (l + u) `div` 2 - -ceiling :: Natural -> Natural -ceiling n = - if sqrtn * sqrtn == n then sqrtn else sqrtn + 1 - where - sqrtn = Arithmetic.SquareRoot.floor n - -continuedFraction :: Natural -> ContinuedFraction.ContinuedFraction -continuedFraction n = - ContinuedFraction.ContinuedFraction (sqrtn,qs) - where - sqrtn = Arithmetic.SquareRoot.floor n - - ps = continuedFractionPeriodicTail n sqrtn - - qs = if null ps then [] else cycle ps - -continuedFractionPeriodic :: Natural -> [Natural] -continuedFractionPeriodic n = - continuedFractionPeriodicTail n sqrtn - where - sqrtn = Arithmetic.SquareRoot.floor n - -continuedFractionPeriodicTail :: Natural -> Natural -> [Natural] -continuedFractionPeriodicTail n sqrtn = - List.unfoldr go (sqrtn,sqrtd) - where - sqrtd = n - sqrtn * sqrtn - --- (sqrt(n) + a) / b = c + 1 / x ==> --- x = b / (sqrt(n) + a - c * b) --- = b / (sqrt(n) - (c * b - a)) --- = (b * (sqrt(n) + (c * b - a))) / (n - (c * b - a)^2) - advance (a,b) = - (c,(d,e)) - where - c = (sqrtn + a) `div` b - d = c * b - a - e = (n - d * d) `div` b - - go (a,b) = - case b of - 0 -> Nothing - 1 -> Just (2 * a, (0,0)) - _ -> Just (advance (a,b)) diff --git a/src/Arithmetic/Utility.hs b/src/Arithmetic/Utility.hs new file mode 100644 --- /dev/null +++ b/src/Arithmetic/Utility.hs @@ -0,0 +1,46 @@ +{- | +module: Arithmetic.Utility +description: Utility functions +license: MIT + +maintainer: Joe Leslie-Hurd +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) diff --git a/src/Arithmetic/Utility/Heap.hs b/src/Arithmetic/Utility/Heap.hs new file mode 100644 --- /dev/null +++ b/src/Arithmetic/Utility/Heap.hs @@ -0,0 +1,84 @@ +{- | +module: Arithmetic.Utility.Heap +description: Leftist heaps +license: MIT + +maintainer: Joe Leslie-Hurd +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 diff --git a/src/Arithmetic/Williams.hs b/src/Arithmetic/Williams.hs new file mode 100644 --- /dev/null +++ b/src/Arithmetic/Williams.hs @@ -0,0 +1,117 @@ +{- | +module: Arithmetic.Williams +description: Williams p+1 factorization method +license: MIT + +maintainer: Joe Leslie-Hurd +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 diff --git a/src/Main.hs b/src/Main.hs --- a/src/Main.hs +++ b/src/Main.hs @@ -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 () diff --git a/src/Test.hs b/src/Test.hs --- a/src/Test.hs +++ b/src/Test.hs @@ -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 @@ -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)