1-- |
2-- Module      : Crypto.Number.Prime
3-- License     : BSD-style
4-- Maintainer  : Vincent Hanquez <vincent@snarc.org>
5-- Stability   : experimental
6-- Portability : Good
7
8{-# LANGUAGE BangPatterns #-}
9module Crypto.Number.Prime
10    (
11      generatePrime
12    , generateSafePrime
13    , isProbablyPrime
14    , findPrimeFrom
15    , findPrimeFromWith
16    , primalityTestMillerRabin
17    , primalityTestNaive
18    , primalityTestFermat
19    , isCoprime
20    ) where
21
22import Crypto.Number.Compat
23import Crypto.Number.Generate
24import Crypto.Number.Basic (sqrti, gcde)
25import Crypto.Number.ModArithmetic (expSafe)
26import Crypto.Random.Types
27import Crypto.Random.Probabilistic
28import Crypto.Error
29
30import Data.Bits
31
32-- | Returns if the number is probably prime.
33-- First a list of small primes are implicitely tested for divisibility,
34-- then a fermat primality test is used with arbitrary numbers and
35-- then the Miller Rabin algorithm is used with an accuracy of 30 recursions.
36isProbablyPrime :: Integer -> Bool
37isProbablyPrime !n
38    | any (\p -> p `divides` n) (filter (< n) firstPrimes) = False
39    | n >= 2 && n <= 2903                                  = True
40    | primalityTestFermat 50 (n `div` 2) n                 = primalityTestMillerRabin 30 n
41    | otherwise                                            = False
42
43-- | Generate a prime number of the required bitsize (i.e. in the range
44-- [2^(b-1)+2^(b-2), 2^b)).
45--
46-- May throw a 'CryptoError_PrimeSizeInvalid' if the requested size is less
47-- than 5 bits, as the smallest prime meeting these conditions is 29.
48-- This function requires that the two highest bits are set, so that when
49-- multiplied with another prime to create a key, it is guaranteed to be of
50-- the proper size.
51generatePrime :: MonadRandom m => Int -> m Integer
52generatePrime bits = do
53    if bits < 5 then
54        throwCryptoError $ CryptoFailed $ CryptoError_PrimeSizeInvalid
55    else do
56        sp <- generateParams bits (Just SetTwoHighest) True
57        let prime = findPrimeFrom sp
58        if prime < 1 `shiftL` bits then
59            return $ prime
60        else generatePrime bits
61
62-- | Generate a prime number of the form 2p+1 where p is also prime.
63-- it is also knowed as a Sophie Germaine prime or safe prime.
64--
65-- The number of safe prime is significantly smaller to the number of prime,
66-- as such it shouldn't be used if this number is supposed to be kept safe.
67--
68-- May throw a 'CryptoError_PrimeSizeInvalid' if the requested size is less than
69-- 6 bits, as the smallest safe prime with the two highest bits set is 59.
70generateSafePrime :: MonadRandom m => Int -> m Integer
71generateSafePrime bits = do
72    if bits < 6 then
73        throwCryptoError $ CryptoFailed $ CryptoError_PrimeSizeInvalid
74    else do
75        sp <- generateParams bits (Just SetTwoHighest) True
76        let p = findPrimeFromWith (\i -> isProbablyPrime (2*i+1)) (sp `div` 2)
77        let val = 2 * p + 1
78        if val < 1 `shiftL` bits then
79            return $ val
80        else generateSafePrime bits
81
82-- | Find a prime from a starting point where the property hold.
83findPrimeFromWith :: (Integer -> Bool) -> Integer -> Integer
84findPrimeFromWith prop !n
85    | even n        = findPrimeFromWith prop (n+1)
86    | otherwise     =
87        if not (isProbablyPrime n)
88            then findPrimeFromWith prop (n+2)
89            else
90                if prop n
91                    then n
92                    else findPrimeFromWith prop (n+2)
93
94-- | Find a prime from a starting point with no specific property.
95findPrimeFrom :: Integer -> Integer
96findPrimeFrom n =
97    case gmpNextPrime n of
98        GmpSupported p -> p
99        GmpUnsupported -> findPrimeFromWith (\_ -> True) n
100
101-- | Miller Rabin algorithm return if the number is probably prime or composite.
102-- the tries parameter is the number of recursion, that determines the accuracy of the test.
103primalityTestMillerRabin :: Int -> Integer -> Bool
104primalityTestMillerRabin tries !n =
105    case gmpTestPrimeMillerRabin tries n of
106        GmpSupported b -> b
107        GmpUnsupported -> probabilistic run
108  where
109    run
110        | n <= 3     = error "Miller-Rabin requires tested value to be > 3"
111        | even n     = return False
112        | tries <= 0 = error "Miller-Rabin tries need to be > 0"
113        | otherwise  = loop <$> generateTries tries
114
115    !nm1 = n-1
116    !nm2 = n-2
117
118    (!s,!d) = (factorise 0 nm1)
119
120    generateTries 0 = return []
121    generateTries t = do
122        v  <- generateBetween 2 nm2
123        vs <- generateTries (t-1)
124        return (v:vs)
125
126    -- factorise n-1 into the form 2^s*d
127    factorise :: Integer -> Integer -> (Integer, Integer)
128    factorise !si !vi
129        | vi `testBit` 0 = (si, vi)
130        | otherwise     = factorise (si+1) (vi `shiftR` 1) -- probably faster to not shift v continuously, but just once.
131    expmod = expSafe
132
133    -- when iteration reach zero, we have a probable prime
134    loop []     = True
135    loop (w:ws) = let x = expmod w d n
136                   in if x == (1 :: Integer) || x == nm1
137                          then loop ws
138                          else loop' ws ((x*x) `mod` n) 1
139
140    -- loop from 1 to s-1. if we reach the end then it's composite
141    loop' ws !x2 !r
142        | r == s    = False
143        | x2 == 1   = False
144        | x2 /= nm1 = loop' ws ((x2*x2) `mod` n) (r+1)
145        | otherwise = loop ws
146
147{-
148    n < z -> witness to test
149              1373653 [2,3]
150              9080191 [31,73]
151              4759123141 [2,7,61]
152              2152302898747 [2,3,5,7,11]
153              3474749660383 [2,3,5,7,11,13]
154              341550071728321 [2,3,5,7,11,13,17]
155-}
156
157-- | Probabilitic Test using Fermat primility test.
158-- Beware of Carmichael numbers that are Fermat liars, i.e. this test
159-- is useless for them. always combines with some other test.
160primalityTestFermat :: Int -- ^ number of iterations of the algorithm
161                    -> Integer -- ^ starting a
162                    -> Integer -- ^ number to test for primality
163                    -> Bool
164primalityTestFermat n a p = and $ map expTest [a..(a+fromIntegral n)]
165    where !pm1 = p-1
166          expTest i = expSafe i pm1 p == 1
167
168-- | Test naively is integer is prime.
169-- while naive, we skip even number and stop iteration at i > sqrt(n)
170primalityTestNaive :: Integer -> Bool
171primalityTestNaive n
172    | n <= 1    = False
173    | n == 2    = True
174    | even n    = False
175    | otherwise = search 3
176        where !ubound = snd $ sqrti n
177              search !i
178                  | i > ubound    = True
179                  | i `divides` n = False
180                  | otherwise     = search (i+2)
181
182-- | Test is two integer are coprime to each other
183isCoprime :: Integer -> Integer -> Bool
184isCoprime m n = case gcde m n of (_,_,d) -> d == 1
185
186-- | List of the first primes till 2903.
187firstPrimes :: [Integer]
188firstPrimes =
189    [ 2    , 3    , 5    , 7    , 11   , 13   , 17   , 19   , 23   , 29
190    , 31   , 37   , 41   , 43   , 47   , 53   , 59   , 61   , 67   , 71
191    , 73   , 79   , 83   , 89   , 97   , 101  , 103  , 107  , 109  , 113
192    , 127  , 131  , 137  , 139  , 149  , 151  , 157  , 163  , 167  , 173
193    , 179  , 181  , 191  , 193  , 197  , 199  , 211  , 223  , 227  , 229
194    , 233  , 239  , 241  , 251  , 257  , 263  , 269  , 271  , 277  , 281
195    , 283  , 293  , 307  , 311  , 313  , 317  , 331  , 337  , 347  , 349
196    , 353  , 359  , 367  , 373  , 379  , 383  , 389  , 397  , 401  , 409
197    , 419  , 421  , 431  , 433  , 439  , 443  , 449  , 457  , 461  , 463
198    , 467  , 479  , 487  , 491  , 499  , 503  , 509  , 521  , 523  , 541
199    , 547  , 557  , 563  , 569  , 571  , 577  , 587  , 593  , 599  , 601
200    , 607  , 613  , 617  , 619  , 631  , 641  , 643  , 647  , 653  , 659
201    , 661  , 673  , 677  , 683  , 691  , 701  , 709  , 719  , 727  , 733
202    , 739  , 743  , 751  , 757  , 761  , 769  , 773  , 787  , 797  , 809
203    , 811  , 821  , 823  , 827  , 829  , 839  , 853  , 857  , 859  , 863
204    , 877  , 881  , 883  , 887  , 907  , 911  , 919  , 929  , 937  , 941
205    , 947  , 953  , 967  , 971  , 977  , 983  , 991  , 997  , 1009 , 1013
206    , 1019 , 1021 , 1031 , 1033 , 1039 , 1049 , 1051 , 1061 , 1063 , 1069
207    , 1087 , 1091 , 1093 , 1097 , 1103 , 1109 , 1117 , 1123 , 1129 , 1151
208    , 1153 , 1163 , 1171 , 1181 , 1187 , 1193 , 1201 , 1213 , 1217 , 1223
209    , 1229 , 1231 , 1237 , 1249 , 1259 , 1277 , 1279 , 1283 , 1289 , 1291
210    , 1297 , 1301 , 1303 , 1307 , 1319 , 1321 , 1327 , 1361 , 1367 , 1373
211    , 1381 , 1399 , 1409 , 1423 , 1427 , 1429 , 1433 , 1439 , 1447 , 1451
212    , 1453 , 1459 , 1471 , 1481 , 1483 , 1487 , 1489 , 1493 , 1499 , 1511
213    , 1523 , 1531 , 1543 , 1549 , 1553 , 1559 , 1567 , 1571 , 1579 , 1583
214    , 1597 , 1601 , 1607 , 1609 , 1613 , 1619 , 1621 , 1627 , 1637 , 1657
215    , 1663 , 1667 , 1669 , 1693 , 1697 , 1699 , 1709 , 1721 , 1723 , 1733
216    , 1741 , 1747 , 1753 , 1759 , 1777 , 1783 , 1787 , 1789 , 1801 , 1811
217    , 1823 , 1831 , 1847 , 1861 , 1867 , 1871 , 1873 , 1877 , 1879 , 1889
218    , 1901 , 1907 , 1913 , 1931 , 1933 , 1949 , 1951 , 1973 , 1979 , 1987
219    , 1993 , 1997 , 1999 , 2003 , 2011 , 2017 , 2027 , 2029 , 2039 , 2053
220    , 2063 , 2069 , 2081 , 2083 , 2087 , 2089 , 2099 , 2111 , 2113 , 2129
221    , 2131 , 2137 , 2141 , 2143 , 2153 , 2161 , 2179 , 2203 , 2207 , 2213
222    , 2221 , 2237 , 2239 , 2243 , 2251 , 2267 , 2269 , 2273 , 2281 , 2287
223    , 2293 , 2297 , 2309 , 2311 , 2333 , 2339 , 2341 , 2347 , 2351 , 2357
224    , 2371 , 2377 , 2381 , 2383 , 2389 , 2393 , 2399 , 2411 , 2417 , 2423
225    , 2437 , 2441 , 2447 , 2459 , 2467 , 2473 , 2477 , 2503 , 2521 , 2531
226    , 2539 , 2543 , 2549 , 2551 , 2557 , 2579 , 2591 , 2593 , 2609 , 2617
227    , 2621 , 2633 , 2647 , 2657 , 2659 , 2663 , 2671 , 2677 , 2683 , 2687
228    , 2689 , 2693 , 2699 , 2707 , 2711 , 2713 , 2719 , 2729 , 2731 , 2741
229    , 2749 , 2753 , 2767 , 2777 , 2789 , 2791 , 2797 , 2801 , 2803 , 2819
230    , 2833 , 2837 , 2843 , 2851 , 2857 , 2861 , 2879 , 2887 , 2897 , 2903
231    ]
232
233{-# INLINE divides #-}
234divides :: Integer -> Integer -> Bool
235divides i n = n `mod` i == 0
236