{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE ViewPatterns #-}
module Math.NumberTheory.Primes.Testing.Certified
( isCertifiedPrime
) where
import Data.List (foldl')
import Data.Bits ((.&.))
import Data.Mod
import Data.Proxy
import GHC.Integer.GMP.Internals (powModInteger)
import GHC.TypeNats (SomeNat(..), someNatVal)
import Math.NumberTheory.Roots (integerSquareRoot)
import Math.NumberTheory.Primes (unPrime)
import Math.NumberTheory.Primes.Factorisation.TrialDivision (trialDivisionPrimeTo, trialDivisionTo, trialDivisionWith)
import Math.NumberTheory.Primes.Factorisation.Montgomery (montgomeryFactorisation, smallFactors, findParms)
import Math.NumberTheory.Primes.Testing.Probabilistic (bailliePSW, isPrime, isStrongFermatPP, lucasTest)
import Math.NumberTheory.Primes.Sieve.Eratosthenes (primeList, primeSieve)
import Math.NumberTheory.Utils (splitOff)
isCertifiedPrime :: Integer -> Bool
isCertifiedPrime n
| n < 0 = isCertifiedPrime (-n)
| otherwise = isPrime n && ((n < bpbd) || checkPrimalityProof (certifyBPSW n))
where
bpbd = 100000000000000000
data PrimalityProof
= Pocklington { cprime :: !Integer
, _factorisedPart, _cofactor :: !Integer
, _knownFactors :: ![(Integer, Word, Integer, PrimalityProof)]
}
| TrialDivision { cprime :: !Integer
, _tdLimit :: !Integer }
| Trivial { cprime :: !Integer
}
deriving Show
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 = powModInteger base (pm1 `quot` pf) p
y = powModInteger x pf p
isTrivialPrime :: Integer -> Bool
isTrivialPrime n = n `elem` trivialPrimes
trivialPrimes :: [Integer]
trivialPrimes = [2,3,5,7,11,13,17,19,23,29]
smallCert :: Integer -> PrimalityProof
smallCert n
| n < 30 = Trivial n
| otherwise = TrialDivision n (integerSquareRoot n + 1)
certify :: Integer -> Maybe PrimalityProof
certify n
| n < 2 = error "Only numbers larger than 1 can be certified"
| n < 31 = case trialDivisionWith trivialPrimes n of
((p,_):_) | p < n -> Nothing
| otherwise -> Just (Trivial n)
_ -> error "Impossible"
| n < billi = let r2 = integerSquareRoot n + 2 in
case trialDivisionTo r2 n of
((p,_):_) | p < n -> Nothing
| otherwise -> Just (TrialDivision n r2)
_ -> error "Impossible"
| otherwise = case smallFactors (fromInteger (abs n)) of
([], Just _) | not (isStrongFermatPP n 2) -> Nothing
| not (lucasTest n) -> Nothing
| otherwise -> Just (certifyBPSW n)
((toInteger -> p,_):_, _)
| p == n -> Just (TrialDivision n (min 100000 n))
| otherwise -> Nothing
_ -> error ("***Error factorising " ++ show n ++ "! Please report this to maintainer of arithmoi.")
where
billi = 1000000000000
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
Nothing -> error ("***Error in factorisation code: " ++ show p
++ " was supposed to be prime but isn't.\n"
++ "Please report this to the maintainer.\n\n")
Just ppr ->(p,e,bs,ppr)
where
q = nm1 `quot` p
x = powModInteger bs q n
y = powModInteger x p n
g = gcd n (x-1)
findDecomposition :: Integer -> (Integer, [(Integer, Word, Bool)], Integer)
findDecomposition n = go 1 n [] prms
where
sr = integerSquareRoot n
pbd = min 1000000 (sr+20)
prms = map unPrime $ 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)
| 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
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
findLoop :: Integer -> Word -> Word -> Word -> Integer -> Either Integer Integer
findLoop _ _ _ 0 s = Left s
findLoop n lo hi ct s
| n <= s+2 = Left 6
| otherwise = case someNatVal (fromInteger n) of
SomeNat (_ :: Proxy t) -> case montgomeryFactorisation lo hi (fromInteger s :: Mod t) of
Nothing -> findLoop n lo hi (ct-1) (s+1)
Just fct
| bailliePSW fct -> Right fct
| otherwise -> Right (findFactor fct 8 (s+1))
bpswMessage :: Integer -> String
bpswMessage n = unlines
[ "\n***Congratulations! You found a Baillie PSW pseudoprime!"
, "Please report this finding to the maintainers:"
, "<daniel.is.fischer@googlemail.com>,"
, "<andrew.lelechenko@gmail.com>"
, "The number in question is:\n"
, show n
, "\nOther parties like wikipedia might also be interested."
, "\nSorry for aborting your programm, but this is a major discovery."
]
found :: Integer -> String
found g = "\nA nontrivial divisor is:\n" ++ show g
fermat :: Integer -> String
fermat b = "\nThe Fermat test fails for base\n" ++ show b