arithmoi 0.1.0.2 → 0.2.0.0

raw patch · 13 files changed

+684/−13 lines, 13 filesPVP ok

version bump matches the API change (PVP)

API changes (from Hackage documentation)

Files

• Changes +2/−0
• Math/NumberTheory/Moduli.hs +1/−1
• Math/NumberTheory/Primes/Factorisation.hs +3/−0
• Math/NumberTheory/Primes/Factorisation/Certified.hs +170/−0
• Math/NumberTheory/Primes/Factorisation/Montgomery.hs +2/−1
• Math/NumberTheory/Primes/Factorisation/TrialDivision.hs +73/−0
• Math/NumberTheory/Primes/Testing.hs +6/−8
• Math/NumberTheory/Primes/Testing/Certificates.hs +36/−0
• Math/NumberTheory/Primes/Testing/Certificates/Internal.hs +349/−0
• Math/NumberTheory/Primes/Testing/Certified.hs +28/−0
• Math/NumberTheory/Utils.hs +4/−0
• TODO +3/−1
• arithmoi.cabal +7/−2

@@ -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:

@@ -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

@@ -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

@@ -0,0 +1,170 @@+-- |+-- Module:      Math.NumberTheory.Primes.Factorisation.Certified+-- Copyright:   (c) 2011 Daniel Fischer+-- Licence:     MIT+-- Maintainer:  Daniel Fischer <daniel.is.fischer@googlemail.com>+-- 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+                    ]

@@ -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)

@@ -0,0 +1,73 @@+-- |+-- Module:      Math.NumberTheory.Primes.Factorisation.TrialDivision+-- Copyright:   (c) 2011 Daniel Fischer+-- Licence:     MIT+-- Maintainer:  Daniel Fischer <daniel.is.fischer@googlemail.com>+-- 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))

@@ -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.

@@ -0,0 +1,36 @@+-- |+-- Module:      Math.NumberTheory.Primes.Testing.Certificates+-- Copyright:   (c) 2011 Daniel Fischer+-- Licence:     MIT+-- Maintainer:  Daniel Fischer <daniel.is.fischer@googlemail.com>+-- 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+

@@ -0,0 +1,349 @@+-- |+-- Module:      Math.NumberTheory.Primes.Testing.Certificates.Internal+-- Copyright:   (c) 2011 Daniel Fischer+-- Licence:     MIT+-- Maintainer:  Daniel Fischer <daniel.is.fischer@googlemail.com>+-- 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,"+                    , "<daniel.is.fischer@googlemail.com>"+                    , "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

@@ -0,0 +1,28 @@+-- |+-- Module:      Math.NumberTheory.Primes.Testing.Certified+-- Copyright:   (c) 2011 Daniel Fischer+-- Licence:     MIT+-- Maintainer:  Daniel Fischer <daniel.is.fischer@googlemail.com>+-- 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.

@@ -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

@@ -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?

@@ -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