diff --git a/Changes b/Changes --- a/Changes +++ b/Changes @@ -1,3 +1,5 @@ +0.2.0.0: + Added certificates and certified testing/factorisation 0.1.0.2: Fixed doc bugs 0.1.0.1: diff --git a/Math/NumberTheory/Moduli.hs b/Math/NumberTheory/Moduli.hs --- a/Math/NumberTheory/Moduli.hs +++ b/Math/NumberTheory/Moduli.hs @@ -189,7 +189,7 @@ bse' = if base < 0 || mdl' <= base then base `mod` mdl' else base -- | Specialised worker without input checks. Makes the same assumptions --- as the general version. +-- as the general version 'powerMod''. powerModInteger' :: Integer -> Integer -> Integer -> Integer powerModInteger' base expo md = go e1 w1 1 base where diff --git a/Math/NumberTheory/Primes/Factorisation.hs b/Math/NumberTheory/Primes/Factorisation.hs --- a/Math/NumberTheory/Primes/Factorisation.hs +++ b/Math/NumberTheory/Primes/Factorisation.hs @@ -26,6 +26,8 @@ , FactorSieve , factorSieve , sieveFactor + -- *** Trial division + , trialDivisionTo -- ** Partial factorisation , smallFactors , stdGenFactorisation @@ -65,6 +67,7 @@ import Math.NumberTheory.Primes.Factorisation.Utils import Math.NumberTheory.Primes.Factorisation.Montgomery +import Math.NumberTheory.Primes.Factorisation.TrialDivision import Math.NumberTheory.Primes.Sieve.Misc -- $algorithm diff --git a/Math/NumberTheory/Primes/Factorisation/Certified.hs b/Math/NumberTheory/Primes/Factorisation/Certified.hs new file mode 100644 --- /dev/null +++ b/Math/NumberTheory/Primes/Factorisation/Certified.hs @@ -0,0 +1,170 @@ +-- | +-- Module: Math.NumberTheory.Primes.Factorisation.Certified +-- Copyright: (c) 2011 Daniel Fischer +-- Licence: MIT +-- Maintainer: Daniel Fischer +-- Stability: Provisional +-- Portability: Non-portable (GHC extensions) +-- +-- Factorisation proving the primality of the found factors. +-- +-- For large numbers, this will be very slow in general. +-- Use only if you're paranoid or must be /really/ sure. +{-# LANGUAGE BangPatterns #-} +module Math.NumberTheory.Primes.Factorisation.Certified + ( certifiedFactorisation + , certificateFactorisation + , provenFactorisation + ) where + +import System.Random +import Control.Monad.State.Strict +import Control.Applicative +import Data.Maybe +import Data.Bits + +import Math.NumberTheory.Primes.Factorisation.Montgomery +import Math.NumberTheory.Primes.Testing.Certificates.Internal +import Math.NumberTheory.Primes.Testing.Probabilistic + +-- | @'certifiedFactorisation' n@ produces the prime factorisation +-- of @n@, proving the primality of the factors, but doesn't report the proofs. +certifiedFactorisation :: Integer -> [(Integer,Int)] +certifiedFactorisation = map fst . certificateFactorisation + +-- | @'certificateFactorisation' n@ produces a 'provenFactorisation' +-- with a default bound of @100000@. +certificateFactorisation :: Integer -> [((Integer,Int),PrimalityProof)] +certificateFactorisation n = provenFactorisation 100000 n + +-- | @'provenFactorisation' bound n@ constructs a the prime factorisation of @n@ +-- (which must be positive) together with proofs of primality of the factors, +-- using trial division up to @bound@ (which is arbitrarily replaced by @2000@ +-- if the supplied value is smaller) and elliptic curve factorisation for the +-- remaining factors if necessary. +-- +-- Construction of primality proofs can take a /very/ long time, so this +-- will usually be slow (but should be faster than using 'factorise' and +-- proving the primality of the factors from scratch). +provenFactorisation :: Integer -> Integer -> [((Integer,Int),PrimalityProof)] +provenFactorisation _ 1 = [] +provenFactorisation bd n + | n < 2 = error "provenFactorisation: argument not positive" + | bd < 2000 = provenFactorisation 2000 n + | otherwise = test $ + case smallFactors bd n of + (sfs,mb) -> map (\t@(p,_) -> (t, smallCert p)) sfs + ++ case mb of + Nothing -> [] + Just k -> certiFactorisation (Just $ bd*(bd+2)) primeCheck (randomR . (,) 6) + (mkStdGen $ fromIntegral n `xor` 0xdeadbeef) Nothing k + +-- | verify that we indeed have a correct primality proof +test :: [((Integer,Int),PrimalityProof)] -> [((Integer,Int),PrimalityProof)] +test (t@((p,_),prf):more) + | p == cprime prf && checkPrimalityProof prf = t : test more + | otherwise = error (invalid p prf) +test [] = [] + +-- | produce a proof of primality for primes +-- Only called for (not too small) numbers known to have no small prime factors, +-- so we can directly use BPSW without trial division. +primeCheck :: Integer -> Maybe PrimalityProof +primeCheck n + | bailliePSW n = case certifyBPSW n of + proof@Pocklington{} -> Just proof + _ -> Nothing + | otherwise = Nothing + +-- | produce a certified factorisation +-- Assumes all small prime factors have been stripped before. +-- Since it is not exported, that is known to hold. +-- This is a near duplicate of 'curveFactorisation', I should some time +-- clean this up. +certiFactorisation :: Maybe Integer -- ^ Lower bound for composite divisors + -> (Integer -> Maybe PrimalityProof) + -- ^ A primality test + -> (Integer -> g -> (Integer,g)) -- ^ A PRNG + -> g -- ^ Initial PRNG state + -> Maybe Int -- ^ Estimated number of digits of the smallest prime factor + -> Integer -- ^ The number to factorise + -> [((Integer,Int),PrimalityProof)] + -- ^ List of prime factors, exponents and primality proofs +certiFactorisation primeBound primeTest prng seed mbdigs n + = case ptest n of + Just proof -> [((n,1),proof)] + Nothing -> evalState (fact n digits) seed + where + digits = fromMaybe 8 mbdigs + mult 1 xs = xs + mult j xs = [((p,j*k),c) | ((p,k),c) <- xs] + vdb xs = [(p,2*e) | (p,e) <- xs] + dbl (u,v) = (mult 2 u, vdb v) + ptest = case primeBound of + Just bd -> \k -> if k <= bd then (Just $ smallCert k) else primeTest k + Nothing -> primeTest + rndR k = state (\gen -> prng k gen) + fact m digs = do let (b1,b2,ct) = findParms digs + (pfs,cfs) <- repFact m b1 b2 ct + if null cfs + then return pfs + else do + nfs <- forM cfs $ \(k,j) -> + mult j <$> fact k (if null pfs then digs+4 else digs) + return (mergeAll $ pfs:nfs) + repFact m b1 b2 count + | count < 0 = return ([],[(m,1)]) + | otherwise = do + s <- rndR m + case montgomeryFactorisation m b1 b2 s of + Nothing -> repFact m b1 b2 (count-1) + Just d -> do + let !cof = m `quot` d + case gcd cof d of + 1 -> do + (dp,dc) <- case ptest d of + Just proof -> return ([((d,1),proof)],[]) + Nothing -> repFact d b1 b2 (count-1) + (cp,cc) <- case ptest cof of + Just proof -> return ([((cof,1),proof)],[]) + Nothing -> repFact cof b1 b2 (count-1) + return (merge dp cp, dc ++ cc) + g -> do + let d' = d `quot` g + c' = cof `quot` g + (dp,dc) <- case ptest d' of + Just proof -> return ([((d',1),proof)],[]) + Nothing -> repFact d' b1 b2 (count-1) + (cp,cc) <- case ptest c' of + Just proof -> return ([((c',1),proof)],[]) + Nothing -> repFact c' b1 b2 (count-1) + (gp,gc) <- case ptest g of + Just proof -> return ([((g,2),proof)],[]) + Nothing -> dbl <$> repFact g b1 b2 (count-1) + return (mergeAll [dp,cp,gp], dc ++ cc ++ gc) + +-- | merge two lists of factors, so that the result is strictly increasing (wrt the primes) +merge :: [((Integer,Int),PrimalityProof)] -> [((Integer,Int),PrimalityProof)] -> [((Integer,Int),PrimalityProof)] +merge xxs@(x@((p,e),c):xs) yys@(y@((q,d),_):ys) + = case compare p q of + LT -> x : merge xs yys + EQ -> ((p,e+d),c) : merge xs ys + GT -> y : merge xxs ys +merge [] ys = ys +merge xs _ = xs + +-- | merge a list of lists of factors so that the result is strictly increasing (wrt the primes) +mergeAll :: [[((Integer,Int),PrimalityProof)]] -> [((Integer,Int),PrimalityProof)] +mergeAll [] = [] +mergeAll [xs] = xs +mergeAll (xs:ys:zss) = merge (merge xs ys) (mergeAll zss) + +-- | message for an invalid proof, should never be used +invalid :: Integer -> PrimalityProof -> String +invalid p prf = unlines + [ "\nInvalid primality proof constructed, please report this to the package maintainer!" + , "The supposed prime was:\n" + , show p + , "\nThe presumed proof was:\n" + , show prf + ] diff --git a/Math/NumberTheory/Primes/Factorisation/Montgomery.hs b/Math/NumberTheory/Primes/Factorisation/Montgomery.hs --- a/Math/NumberTheory/Primes/Factorisation/Montgomery.hs +++ b/Math/NumberTheory/Primes/Factorisation/Montgomery.hs @@ -37,6 +37,7 @@ , curveFactorisation -- ** Single curve worker , montgomeryFactorisation + , findParms ) where #include "MachDeps.h" @@ -147,7 +148,7 @@ | ptest n = [(n,1)] | otherwise = evalState (fact n digits) seed where - digits = fromMaybe 6 mbdigs + digits = fromMaybe 8 mbdigs mult 1 xs = xs mult j xs = [(p,j*k) | (p,k) <- xs] dbl (u,v) = (mult 2 u, mult 2 v) diff --git a/Math/NumberTheory/Primes/Factorisation/TrialDivision.hs b/Math/NumberTheory/Primes/Factorisation/TrialDivision.hs new file mode 100644 --- /dev/null +++ b/Math/NumberTheory/Primes/Factorisation/TrialDivision.hs @@ -0,0 +1,73 @@ +-- | +-- Module: Math.NumberTheory.Primes.Factorisation.TrialDivision +-- Copyright: (c) 2011 Daniel Fischer +-- Licence: MIT +-- Maintainer: Daniel Fischer +-- Stability: Provisional +-- Portability: Non-portable (GHC extensions) +-- +-- Factorisation and primality testing using trial division. +-- +-- Used to create and check certificates. +-- Currently not exposed because it's not that useful, is it? +-- But the trial...To functions are exported from other modules. +{-# LANGUAGE BangPatterns #-} +module Math.NumberTheory.Primes.Factorisation.TrialDivision + ( trialDivisionWith + , trialDivisionTo + , trialDivisionPrimeWith + , trialDivisionPrimeTo + ) where + +import Math.NumberTheory.Primes.Sieve.Eratosthenes +import Math.NumberTheory.Powers.Squares +import Math.NumberTheory.Utils + +-- | Factorise an 'Integer' using a given list of numbers considered prime. +-- If the list is not a list of primes containing all relevant primes, the +-- result could be surprising. +trialDivisionWith :: [Integer] -> Integer -> [(Integer,Int)] +trialDivisionWith prs n + | n < 0 = trialDivisionWith prs (-n) + | n == 0 = error "trialDivision of 0" + | n == 1 = [] + | otherwise = go n (integerSquareRoot' n) prs + where + go !m !r (p:ps) + | r < p = [(m,1)] + | otherwise = + case splitOff p m of + (0,_) -> go m r ps + (k,q) -> (p,k) : if q == 1 + then [] + else go q (integerSquareRoot' q) ps + go m _ _ = [(m,1)] + +-- | @'trialDivisionTo' bound n@ produces a factorisation of @n@ using the +-- primes @<= bound@. If @n@ has prime divisors @> bound@, the last entry +-- in the list is the product of all these. If @n <= bound^2@, this is a +-- full factorisation, but very slow if @n@ has large prime divisors. +trialDivisionTo :: Integer -> Integer -> [(Integer,Int)] +trialDivisionTo bd + | bd < 100 = trialDivisionTo 100 + | bd < 10000000 = trialDivisionWith (primeList $ primeSieve bd) + | otherwise = trialDivisionWith (takeWhile (<= bd) $ (psieveList >>= primeList)) + +-- | Check whether a number is coprime to all of the numbers in the list +-- (assuming that list contains only numbers > 1 and is ascending). +trialDivisionPrimeWith :: [Integer] -> Integer -> Bool +trialDivisionPrimeWith prs n + | n < 0 = trialDivisionPrimeWith prs (-n) + | n < 2 = False + | otherwise = go n (integerSquareRoot' n) prs + where + go !m !r (p:ps) = r < p || m `rem` p /= 0 && go m r ps + go _ _ _ = True + +-- | @'trialDivisionPrimeTo' bound n@ tests whether @n@ is coprime to all primes @<= bound@. +-- If @n <= bound^2@, this is a full prime test, but very slow if @n@ has no small prime divisors. +trialDivisionPrimeTo :: Integer -> Integer -> Bool +trialDivisionPrimeTo bd + | bd < 100 = trialDivisionPrimeTo 100 + | bd < 10000000 = trialDivisionPrimeWith (primeList $ primeSieve bd) + | otherwise = trialDivisionPrimeWith (takeWhile (<= bd) $ (psieveList >>= primeList)) diff --git a/Math/NumberTheory/Primes/Testing.hs b/Math/NumberTheory/Primes/Testing.hs --- a/Math/NumberTheory/Primes/Testing.hs +++ b/Math/NumberTheory/Primes/Testing.hs @@ -10,7 +10,8 @@ module Math.NumberTheory.Primes.Testing ( -- * Standard tests isPrime - -- $certificates + , isCertifiedPrime + -- * Partial tests , bailliePSW , millerRabinV , isStrongFermatPP @@ -18,9 +19,13 @@ -- * Using a sieve , FactorSieve , fsIsPrime + -- * Trial division + , trialDivisionPrimeTo ) where import Math.NumberTheory.Primes.Testing.Probabilistic +import Math.NumberTheory.Primes.Testing.Certified +import Math.NumberTheory.Primes.Factorisation.TrialDivision import Math.NumberTheory.Primes.Sieve.Misc -- | Test primality using a 'FactorSieve'. If @n@ is out of bounds @@ -31,10 +36,3 @@ | n <= fromIntegral (fsBound fs) = fsPrimeTest fs n | otherwise = isPrime n --- $certificates --- --- The tests in this module may wrongly consider some composite numbers as prime. --- For the Baillie-PSW test, no pseudoprimes are known, and it is known that none --- exist below @2^64@, so for most practical purposes it can be regarded as conclusive. --- Nevertheless, it is desirable to certify numbers passing it as primes (or find that --- they are composite). The addition of prime certificates is planned for the next release. diff --git a/Math/NumberTheory/Primes/Testing/Certificates.hs b/Math/NumberTheory/Primes/Testing/Certificates.hs new file mode 100644 --- /dev/null +++ b/Math/NumberTheory/Primes/Testing/Certificates.hs @@ -0,0 +1,36 @@ +-- | +-- Module: Math.NumberTheory.Primes.Testing.Certificates +-- Copyright: (c) 2011 Daniel Fischer +-- Licence: MIT +-- Maintainer: Daniel Fischer +-- Stability: Provisional +-- Portability: Non-portable (GHC extensions) +-- +-- Certificates for primality or compositeness. +module Math.NumberTheory.Primes.Testing.Certificates + ( -- * Certificates + Certificate(..) + , argueCertificate + , CompositenessProof + , composite + , PrimalityProof + , cprime + -- * Arguments + , CompositenessArgument(..) + , PrimalityArgument(..) + -- ** Weaken proofs to arguments + , arguePrimality + , argueCompositeness + -- ** Prove valid arguments + , verifyPrimalityArgument + , verifyCompositenessArgument + -- * Determine and prove whether a number is prime or composite + , certify + -- ** Checks for the paranoid + , checkCertificate + , checkCompositenessProof + , checkPrimalityProof + ) where + +import Math.NumberTheory.Primes.Testing.Certificates.Internal + diff --git a/Math/NumberTheory/Primes/Testing/Certificates/Internal.hs b/Math/NumberTheory/Primes/Testing/Certificates/Internal.hs new file mode 100644 --- /dev/null +++ b/Math/NumberTheory/Primes/Testing/Certificates/Internal.hs @@ -0,0 +1,349 @@ +-- | +-- Module: Math.NumberTheory.Primes.Testing.Certificates.Internal +-- Copyright: (c) 2011 Daniel Fischer +-- Licence: MIT +-- Maintainer: Daniel Fischer +-- Stability: Provisional +-- Portability: Non-portable (GHC extensions) +-- +-- Certificates for primality or compositeness. +{-# OPTIONS_HADDOCK hide #-} +module Math.NumberTheory.Primes.Testing.Certificates.Internal + ( Certificate(..) + , CompositenessProof(..) + , PrimalityProof(..) + , CompositenessArgument(..) + , PrimalityArgument(..) + , checkCertificate + , checkCompositenessProof + , checkPrimalityProof + , certify + , trivial + , smallCert + , certifyBPSW + , argueCertificate + , arguePrimality + , argueCompositeness + , verifyPrimalityArgument + , verifyCompositenessArgument + ) where + +import Data.List +import Data.Word +import Data.Bits +import Data.Maybe + +import Math.NumberTheory.Moduli +import Math.NumberTheory.Utils +import Math.NumberTheory.Primes.Factorisation.TrialDivision +import Math.NumberTheory.Primes.Factorisation.Montgomery +import Math.NumberTheory.Primes.Testing.Probabilistic +import Math.NumberTheory.Primes.Sieve.Eratosthenes +import Math.NumberTheory.Powers.Squares + +-- | A certificate of either compositeness or primality of an +-- 'Integer'. Only numbers @> 1@ can be certified, trying to +-- create a certificate for other numbers raises an error. +data Certificate + = Composite !CompositenessProof + | Prime !PrimalityProof + deriving Show + +-- | A proof of compositeness of a positive number. The type is +-- abstract to ensure the validity of proofs. +data CompositenessProof + = Factors { composite, firstFactor, secondFactor :: !Integer } + | StrongFermat { composite, witness :: !Integer } + | LucasSelfridge { composite :: !Integer } + deriving Show + +-- | An argument for compositeness of a number (which must be @> 1@). +-- 'CompositenessProof's translate directly to 'CompositenessArguments', +-- correct arguments can be transformed into proofs. This type allows the +-- manipulation of proofs while maintaining their correctness. +-- The only way to access components of a 'CompositenessProof' except +-- the composite is through this type. +data CompositenessArgument + = Divisors { compo, firstDivisor, secondDivisor :: Integer } + -- ^ @compo == firstDiv*secondDiv@, where all are @> 1@ + | Fermat { compo, fermatBase :: Integer } -- ^ @compo@ fails the strong Fermat test for @fermatBase@ + | Lucas { compo :: Integer } -- ^ @compo@ fails the Lucas-Selfridge test + | Belief { compo :: Integer } -- ^ No particular reason given + deriving (Show, Read, Eq, Ord) + +-- | A proof of primality of a positive number. The type is +-- abstract to ensure the validity of proofs. +data PrimalityProof + = Pocklington { cprime :: !Integer + , factorisedPart, cofactor :: !Integer + , knownFactors :: ![(Integer,Int,Integer,PrimalityProof)] + } + | TrialDivision { cprime, tdLimit :: !Integer } + | Trivial { cprime :: !Integer } + deriving Show + +-- | An argument for primality of a number (which must be @> 1@). +-- 'PrimalityProof's translate directly to 'PrimalityArguments', +-- correct arguments can be transformed into proofs. This type allows the +-- manipulation of proofs while maintaining their correctness. +-- The only way to access components of a 'PrimalityProof' except +-- the prime is through this type. +data PrimalityArgument + = Pock { aprime :: Integer + , largeFactor, smallFactor :: Integer + , factorList :: [(Integer,Int,Integer,PrimalityArgument)] + } -- ^ A suggested Pocklington certificate + | Division { aprime, alimit :: Integer } -- ^ Primality should be provable by trial division to @alimit@ + | Obvious { aprime :: Integer } -- ^ @aprime@ is said to be obviously prime, that holds for primes @< 30@ + | Assumption { aprime :: Integer } -- ^ Primality assumed + deriving (Show, Read, Eq, Ord) + +argueCertificate :: Certificate -> Either CompositenessArgument PrimalityArgument +argueCertificate (Composite proof) = Left (argueCompositeness proof) +argueCertificate (Prime proof) = Right (arguePrimality proof) + +-- | @'arguePrimality'@ transforms a proof of primality into an argument for primality. +arguePrimality :: PrimalityProof -> PrimalityArgument +arguePrimality (TrialDivision p l) = Division p l +arguePrimality (Trivial p) = Obvious p +arguePrimality (Pocklington p a b fcts) = Pock p a b (map argue fcts) + where + argue (x,y,z,prf) = (x,y,z,arguePrimality prf) + +-- | @'verifyPrimalityArgument'@ checks the given argument and constructs a proof from +-- it, if it is valid. For the explicit arguments, this is simple and resonably fast, +-- for an 'Assumption', the verification uses 'certify' and hence may take a long time. +verifyPrimalityArgument :: PrimalityArgument -> Maybe PrimalityProof +verifyPrimalityArgument (Assumption p) + = case certify p of + Composite _ -> Nothing + Prime proof -> Just proof +verifyPrimalityArgument arg + | checkPrimalityProof prf = Just prf + | otherwise = Nothing + where + prf = primProof arg + +-- | not exported, this is the one place where invalid proofs can be constructed +primProof :: PrimalityArgument -> PrimalityProof +primProof (Division p l) = TrialDivision p l +primProof (Obvious p) = Trivial p +primProof (Assumption p) = case certify p of + Composite _ -> Trivial p -- we're faking to not raise an error + Prime proof -> proof +primProof (Pock p a b fcts) = Pocklington p a b (map prove fcts) + where + prove (x,y,z,arg) = (x,y,z,primProof arg) + +-- | @'argueCompositeness'@ transforms a proof of compositeness into an argument +-- for compositeness. +argueCompositeness :: CompositenessProof -> CompositenessArgument +argueCompositeness (Factors c f s) = Divisors c f s +argueCompositeness (StrongFermat c b) = Fermat c b +argueCompositeness (LucasSelfridge c) = Lucas c + +-- | @'verifyCompositenessArgument'@ checks the given argument and constructs a proof from +-- it, if it is valid. For the explicit arguments, this is simple and resonably fast, +-- for a 'Belief', the verification uses 'certify' and hence may take a long time. +verifyCompositenessArgument :: CompositenessArgument -> Maybe CompositenessProof +verifyCompositenessArgument (Belief c) + = case certify c of + Composite proof -> Just proof + Prime _ -> Nothing +verifyCompositenessArgument arg + | checkCompositenessProof prf = Just prf + | otherwise = Nothing + where + prf = compProof arg + +-- | not exported, here is where invalid proofs can be constructed, +-- they must not leak +compProof :: CompositenessArgument -> CompositenessProof +compProof (Divisors c f s) = Factors c f s +compProof (Fermat c b) = StrongFermat c b +compProof (Lucas c) = LucasSelfridge c +compProof (Belief _) = error "Trying to prove by belief" + +-- | Check the validity of a 'Certificate'. Since it should be impossible +-- to construct invalid certificates by the public interface, this should +-- never return 'False'. +checkCertificate :: Certificate -> Bool +checkCertificate (Composite cp) = checkCompositenessProof cp +checkCertificate (Prime pp) = checkPrimalityProof pp + +-- | Check the validity of a 'CompositenessProof'. Since it should be +-- impossible to create invalid proofs by the public interface, this +-- should never return 'False'. +checkCompositenessProof :: CompositenessProof -> Bool +checkCompositenessProof (Factors c a b) = a > 1 && b > 1 && a*b == c +checkCompositenessProof (StrongFermat c w) = w > 1 && c > w && not (isStrongFermatPP c w) +checkCompositenessProof (LucasSelfridge c) = c > 3 && fromIntegral c .&. (1 :: Int) == 1 && lucasTest c + +-- | Check the validity of a 'PrimalityProof'. Since it should be +-- impossible to create invalid proofs by the public interface, this +-- should never return 'False'. +checkPrimalityProof :: PrimalityProof -> Bool +checkPrimalityProof (Trivial n) = isTrivialPrime n +checkPrimalityProof (TrialDivision p b) = p <= b*b && trialDivisionPrimeTo b p +checkPrimalityProof (Pocklington p a b fcts) = b > 0 && a > b && a*b == pm1 && a == ppProd fcts && all verify fcts + where + pm1 = p-1 + ppProd pps = product [pf^e | (pf,e,_,_) <- pps] + verify (pf,_,base,proof) = pf == cprime proof && crit pf base && checkPrimalityProof proof + crit pf base = gcd p (x-1) == 1 && y == 1 + where + x = powerModInteger' base (pm1 `quot` pf) p + y = powerModInteger' x pf p + +-- | @'trivial'@ records a trivially known prime. +-- If the argument is not one of them, an error is raised. +trivial :: Integer -> PrimalityProof +trivial n = fromMaybe oops $ maybeTrivial n + where + oops = error ("trivial: " ++ show n ++ " isn't a trivially known prime.") + +-- | @'maybeTrivial'@ finds out if its argument is a trivially known +-- prime or not and returns the appropriate. +maybeTrivial :: Integer -> Maybe PrimalityProof +maybeTrivial n + | isTrivialPrime n = Just (Trivial n) + | otherwise = Nothing + +-- | @'isTrivialPrime'@ checks whether its argument is a trivially +-- known prime. +isTrivialPrime :: Integer -> Bool +isTrivialPrime n = n `elem` trivialPrimes + +-- | List of trivially known primes. +trivialPrimes :: [Integer] +trivialPrimes = [2,3,5,7,11,13,17,19,23,29] + +-- | Certify a small number. This is not exposed and should only +-- be used where correct. It is always checked after use, though, +-- so it shouldn't be able to lie. +smallCert :: Integer -> PrimalityProof +smallCert n + | n < 30 = Trivial n + | otherwise = TrialDivision n (integerSquareRoot' n + 1) + +-- | @'certify' n@ constructs, for @n > 1@, a proof of either +-- primality or compositeness of @n@. This may take a very long +-- time if the number has no small(ish) prime divisors +certify :: Integer -> Certificate +certify n + | n < 2 = error "Only numbers larger than 1 can be certified" + | n < 31 = case trialDivisionWith trivialPrimes n of + ((p,_):_) | p < n -> Composite (Factors n p (n `quot` p)) + | otherwise -> Prime (Trivial n) + _ -> error "Impossible" + | n < billi = let r2 = integerSquareRoot' n + 2 in + case trialDivisionTo r2 n of + ((p,_):_) | p < n -> Composite (Factors n p (n `quot` p)) + | otherwise -> Prime (TrialDivision n r2) + _ -> error "Impossible" + | otherwise = case smallFactors 100000 n of + ([], Just _) | not (isStrongFermatPP n 2) -> Composite (StrongFermat n 2) + | not (lucasTest n) -> Composite (LucasSelfridge n) + | otherwise -> Prime (certifyBPSW n) -- if it isn't we error and ask for a report. + ((p,_):_, _) | p == n -> Prime (TrialDivision n (min 100000 n)) + | otherwise -> Composite (Factors n p (n `quot` p)) + _ -> error ("***Error factorising " ++ show n ++ "! Please report this to maintainer of arithmoi.") + where + billi = 1000000000000 + +-- | Certify a number known to be not too small, having no small prime divisors and having +-- passed the Baillie PSW test. So we assume it's prime, erroring if not. +-- Since it's presumably a large number, we don't bother with trial division and +-- construct a Pocklington certificate. +certifyBPSW :: Integer -> PrimalityProof +certifyBPSW n = Pocklington n a b kfcts + where + nm1 = n-1 + h = nm1 `quot` 2 + m3 = fromInteger n .&. (3 :: Int) == 3 + (a,pp,b) = findDecomposition nm1 + kfcts0 = map check pp + kfcts = foldl' force [] kfcts0 + force xs t@(_,_,_,prf) = prf `seq` (t:xs) + check (p,e,byTD) = go 2 + where + go bs + | bs > h = error (bpswMessage n) + | x == 1 = if m3 && (p == 2) then (p,e,n-bs,Trivial 2) else go (bs+1) + | g /= 1 = error (bpswMessage n ++ found g) + | y /= 1 = error (bpswMessage n ++ fermat bs) + | byTD = (p,e,bs, smallCert p) + | otherwise = case certify p of + Composite cpr -> error ("***Error in factorisation code: " ++ show p + ++ " was supposed to be prime but isn't.\n" + ++ "Please report this to the maintainer.\n\n" + ++ show cpr) + Prime ppr ->(p,e,bs,ppr) + where + q = nm1 `quot` p + x = powerModInteger' bs q n + y = powerModInteger' x p n + g = gcd n (x-1) + +-- | Find a decomposition of p-1 for the pocklington certificate. +-- Usually bloody slow if p-1 has two (or more) /large/ prime divisors. +findDecomposition :: Integer -> (Integer, [(Integer,Int,Bool)], Integer) +findDecomposition n = go 1 n [] prms + where + sr = integerSquareRoot' n + pbd = min 1000000 (sr+20) + prms = primeList (primeSieve $ pbd) + go a b afs (p:ps) + | a > b = (a,afs,b) + | otherwise = case splitOff p b of + (0,_) -> go a b afs ps + (e,q) -> go (a*p^e) q ((p,e,True):afs) ps + go a b afs [] + | a > b = (a,afs,b) + | bailliePSW b = (b,[(b,1,False)],a) -- Until a Baillie PSW pseudoprime is found, I'm going with this + | e == 0 = error ("Error in factorisation, " ++ show p ++ " was found as a factor of " ++ show b ++ " but isn't.") + | otherwise = go (a*p^e) q ((p,e,False):afs) [] + where + p = findFactor b 8 6 + (e,q) = splitOff p b + +-- | Find a factor of a known composite with approximately digits digits, +-- starting with curve s. Actually, this may loop infinitely, but the +-- loop should not be entered before the heat death of the universe. +findFactor :: Integer -> Int -> Integer -> Integer +findFactor n digits s = case findLoop n lo hi count s of + Left t -> findFactor n (digits+5) t + Right f -> f + where + (lo,hi,count) = findParms digits + +-- | Find a factor or say with which curve to continue. +findLoop :: Integer -> Word -> Word -> Int -> Integer -> Either Integer Integer +findLoop _ _ _ 0 s = Left s +findLoop n lo hi ct s + | n <= s+2 = Left 6 + | otherwise = case montgomeryFactorisation n lo hi s of + Nothing -> findLoop n lo hi (ct-1) (s+1) + Just fct + | bailliePSW fct -> Right fct + | otherwise -> Right (findFactor fct 8 (s+1)) + +-- | Message in the unlikely case a Baillie PSW pseudoprime is found. +bpswMessage :: Integer -> String +bpswMessage n = unlines + [ "\n***Congratulations! You found a Baillie PSW pseudoprime!" + , "Please report this finding to the package maintainer," + , "" + , "The number in question is:\n" + , show n + , "\nOther parties like wikipedia might also be interested." + , "\nSorry for aborting your programme, but this is a major discovery." + ] + +-- | Found a factor +found :: Integer -> String +found g = "\nA nontrivial divisor is:\n" ++ show g + +-- | Fermat failure +fermat :: Integer -> String +fermat b = "\nThe Fermat test fails for base\n" ++ show b diff --git a/Math/NumberTheory/Primes/Testing/Certified.hs b/Math/NumberTheory/Primes/Testing/Certified.hs new file mode 100644 --- /dev/null +++ b/Math/NumberTheory/Primes/Testing/Certified.hs @@ -0,0 +1,28 @@ +-- | +-- Module: Math.NumberTheory.Primes.Testing.Certified +-- Copyright: (c) 2011 Daniel Fischer +-- Licence: MIT +-- Maintainer: Daniel Fischer +-- Stability: Provisional +-- Portability: Non-portable (GHC extensions) +-- +-- Deterministic primality testing. +module Math.NumberTheory.Primes.Testing.Certified (isCertifiedPrime) where + +import Math.NumberTheory.Primes.Testing.Probabilistic +import Math.NumberTheory.Primes.Testing.Certificates.Internal + +-- | @'isCertifiedPrime' n@ tests primality of @n@, first trial division +-- by small primes is performed, then a Baillie PSW test and finally a +-- prime certificate is constructed and verified, provided no step before +-- found @n@ to be composite. Constructing prime certificates can take +-- a /very/ long time, so use this with care. +isCertifiedPrime :: Integer -> Bool +isCertifiedPrime n + | n < 0 = isCertifiedPrime (-n) + | otherwise = isPrime n && ((n < bpbd) || checkPrimalityProof (certifyBPSW n)) + where + bpbd = 100000000000000000 +-- Although it is known that there are no Baillie PSW pseudoprimes below 2^64, +-- use the verified bound 10^17, I don't know whether Gilchrist's result has been +-- verified yet. diff --git a/Math/NumberTheory/Utils.hs b/Math/NumberTheory/Utils.hs --- a/Math/NumberTheory/Utils.hs +++ b/Math/NumberTheory/Utils.hs @@ -178,7 +178,11 @@ -- Int -> Int -> (Int, Int), -- Word -> Word -> (Int, Word) -- #-} +#if __GLASGOW_HASKELL__ >= 700 {-# INLINABLE splitOff #-} +#else +{-# INLINE splitOff #-} +#endif splitOff :: Integral a => a -> a -> (Int, a) splitOff p n = go 0 n where diff --git a/TODO b/TODO --- a/TODO +++ b/TODO @@ -1,5 +1,7 @@ -- Prime Certificates - Atkin sieve - General number field sieve - Portability - Check whether bit twiddling can be as fast as the lookup table for leading and trailing zeros + Using bit twiddling already, faster on my x86_64, not benchmarked on x86 recently, + but it used to be only a marginal difference anyway. +- More Certificates? diff --git a/arithmoi.cabal b/arithmoi.cabal --- a/arithmoi.cabal +++ b/arithmoi.cabal @@ -1,5 +1,5 @@ name : arithmoi -version : 0.1.0.2 +version : 0.2.0.0 cabal-version : >= 1.6 author : Daniel Fischer copyright : (c) 2011 Daniel Fischer @@ -27,7 +27,7 @@ category : Math, Algorithms, Number Theory -tested-with : GHC == 6.12.3, GHC == 7.0.2, GHC == 7.0.3, GHC == 7.2.1 +tested-with : GHC == 6.12.3, GHC == 7.0.2, GHC == 7.0.4, GHC == 7.2.1 extra-source-files : Changes, TODO @@ -48,8 +48,10 @@ Math.NumberTheory.Primes Math.NumberTheory.Primes.Sieve Math.NumberTheory.Primes.Factorisation + Math.NumberTheory.Primes.Factorisation.Certified Math.NumberTheory.Primes.Counting Math.NumberTheory.Primes.Testing + Math.NumberTheory.Primes.Testing.Certificates Math.NumberTheory.Primes.Heap other-modules : Math.NumberTheory.Utils Math.NumberTheory.Logarithms.Internal @@ -58,10 +60,13 @@ Math.NumberTheory.Primes.Counting.Approximate Math.NumberTheory.Primes.Factorisation.Montgomery Math.NumberTheory.Primes.Factorisation.Utils + Math.NumberTheory.Primes.Factorisation.TrialDivision Math.NumberTheory.Primes.Sieve.Eratosthenes Math.NumberTheory.Primes.Sieve.Indexing Math.NumberTheory.Primes.Sieve.Misc Math.NumberTheory.Primes.Testing.Probabilistic + Math.NumberTheory.Primes.Testing.Certified + Math.NumberTheory.Primes.Testing.Certificates.Internal extensions : BangPatterns ghc-options : -O2 -Wall ghc-prof-options : -O2 -auto