{-# LANGUAGE CPP #-}{-
Copyright (C) 2011 Dr. Alistair Ward
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
-}{- |
[@AUTHOR@] Dr. Alistair Ward
[@DESCRIPTION@]
* Determines whether an integer is prime.
* <http://en.wikipedia.org/wiki/Primality_test>.
* <http://primes.utm.edu/index.html>
* CAVEAT: it doesn't determine the prime-factors of composite numbers, just that they exist.
-}moduleFactory.Math.Implementations.Primality(-- * Types-- ** Data-typesAlgorithm(..)-- * Functions-- ** Predicates-- isPrimeByAKS,-- isPrimeByMillerRabin,-- witnessesCompositeness)whereimportControl.Arrow((&&&))importqualifiedControl.DeepSeqimportqualifiedData.Numbers.PrimesimportqualifiedFactory.Data.MonicPolynomialasData.MonicPolynomialimportqualifiedFactory.Data.PolynomialasData.PolynomialimportqualifiedFactory.Data.QuotientRingasData.QuotientRingimportqualifiedFactory.Math.MultiplicativeOrderasMath.MultiplicativeOrderimportqualifiedFactory.Math.PowerasMath.PowerimportqualifiedFactory.Math.PrimalityasMath.PrimalityimportqualifiedFactory.Math.PrimeFactorisationasMath.PrimeFactorisationimportqualifiedToolShed.DefaultableasDefaultable#if MIN_VERSION_parallel(3,0,0)importqualifiedControl.Parallel.Strategies#endif-- | The algorithms by which /primality/-testing has been implemented.dataAlgorithmfactorisationAlgorithm=AKSfactorisationAlgorithm-- ^ <http://en.wikipedia.org/wiki/AKS_primality_test>.|MillerRabin-- ^ <http://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test>.deriving(Eq,Read,Show)instanceDefaultable.Defaultable(AlgorithmfactorisationAlgorithm)wheredefaultValue=MillerRabininstanceMath.PrimeFactorisation.AlgorithmfactorisationAlgorithm=>Math.Primality.Algorithm(AlgorithmfactorisationAlgorithm)whereisPrime_2=True--The only even prime.isPrimealgorithmcandidate|candidate<2||(any((==0).(candidate`rem`)--The candidate has a small prime-factor, and is therefore composite.).filter((candidate>=).(*2)--The candidate must be at least double the small prime, for it to be a potential factor.).take5{-arbitrarily-}$Data.Numbers.Primes.primes--Excludes even numbers, provided at least the 1st prime is tested.)=False|otherwise=(casealgorithmofAKSfactorisationAlgorithm->isPrimeByAKSfactorisationAlgorithmMillerRabin->isPrimeByMillerRabin)candidate{- |
* An implementation of the /Agrawal-Kayal-Saxena/ primality-test; <http://en.wikipedia.org/wiki/AKS_primality_test>,
using the /Lenstra/ and /Pomerance/ algorithm.
* CAVEAT: this deterministic algorithm has a theoretical time-complexity of @O(log^6)@,
and therefore can't compete with the performance of probabilistic ones.
* The /formal polynomials/ used in this algorithm, are conceptually different from /polynomial functions/;
the /indeterminate/ and its powers, are merely used to name a sequence of pigeon-holes in which /coefficients/ are stored,
and is never substituted for a specific value.
This mind-shift, allows one to introduce concepts like /modular/ arithmetic on polynomials,
which merely represent an operation on their coefficients and the pigeon-hole in which they're placed.
[@Manindra Agrawal, Neeraj Kayal and Nitin Saxena@] <http://www.cse.iitk.ac.in/users/manindra/algebra/primality_v6.pdf>.
[@H. W. Lenstra, Jr. and Carl Pomerance@] <http://www.math.dartmouth.edu/~carlp/PDF/complexity12.pdf>.
[@Salembier and Southerington@] <http://ece.gmu.edu/courses/ECE746/project/F06_Project_resources/Salembier_Southerington_AKS.pdf>,
[@R. Crandall and J. Papadopoulos@] <http://images.apple.com/acg/pdf/aks3.pdf>,
[@Andreas Klappenecker@] <http://faculty.cs.tamu.edu/klappi/629/aks.ps>,
[@Vibhor Bhatt and G. K. Patra@] <http://www.cmmacs.ernet.in/cmmacs/Publications/resch_rep/rrcm0307.pdf>,
-}isPrimeByAKS::(Math.PrimeFactorisation.AlgorithmfactorisationAlgorithm,Control.DeepSeq.NFDatai,Integrali)=>factorisationAlgorithm->i->BoolisPrimeByAKSfactorisationAlgorithmn=and[not$Math.Power.isPerfectPowern,--Step 1.Math.Primality.areCoprimen`all`filter(/=n)[2..r],--Step 3.#if MIN_VERSION_parallel(3,0,0)and$Control.Parallel.Strategies.parMapControl.Parallel.Strategies.rdeepseq--Benefits from '+RTS -H100M', which reduces garbage-collections.#elseall#endif(\a->let-- lhs, rhs :: Data.Polynomial.Polynomial i ilhs=Data.Polynomial.raiseModulo(Data.Polynomial.mkLinear1a)n{-power-}n{-modulus-}rhs=Data.Polynomial.mod'(Data.Polynomial.mkPolynomial[(1,n),(a,0)])ninData.QuotientRing.areCongruentModulo(Data.MonicPolynomial.mkMonicPolynomiallhs)(Data.MonicPolynomial.mkMonicPolynomialrhs)(Data.MonicPolynomial.mkMonicPolynomialmodulus)-- Because all these polynomials are /monic/, one can establish /congruence/ using /integer/-division.)[1..floor.(*lg).sqrt$fromIntegralr]--Step 4; (x + a)^n ~ x^n + a mod (x^r - 1, n).]wherelg::Doublelg=logBase2$fromIntegraln-- r :: ir=fst.head.dropWhile((<=floor(Math.Power.squarelg)).snd).map(id&&&Math.MultiplicativeOrder.multiplicativeOrderfactorisationAlgorithmn)$Math.Primality.areCoprimen`filter`[2..]--Step 2.-- modulus :: Data.Polynomial.Polynomial i imodulus=Data.Polynomial.mkPolynomial[(1,r),(negate1,0)]{- |
* Uses the specified 'base' in an attempt to prove the /compositeness/ of an integer.
* This is the opposite of the /Miller Test/; <http://mathworld.wolfram.com/MillersPrimalityTest.html>.
* If the result is 'True', then the candidate is /composite/; regrettably the converse isn't true.
Amongst the set of possible bases, over three-quarters are /witnesses/ to the compositeness of a /composite/ candidate,
the remainder belong to the subset of /liars/.
In consequence, many false results must be accumulated for different bases, to convincingly identify a prime.
-}witnessesCompositeness::Integrali=>i-- ^ Candidate integer.->i->Int->i-- ^ Base.->BoolwitnessesCompositenesscandidateoddRemaindernPowersOfTwobase=all($((`mod`candidate).Math.Power.square)`iterate`Math.Power.raiseModulobaseoddRemaindercandidate--Repeatedly modulo-square.)[(/=1).head,--Check whether the zeroeth modulo-power is incongruent to one.all(/=predcandidate).takenPowersOfTwo--Check whether any modulo-power is incongruent to -1.]{- |
* Repeatedly calls 'witnessesCompositeness', to progressively increase the probability of detecting a /composite/ number,
until ultimately the candidate integer is proven to be prime.
* Should all bases be tested, then the test is deterministic, but at an efficiency /lower/ than performing prime-factorisation.
* The test becomes deterministic, for any candidate integer, when the number of tests reaches the limit defined by /Eric Bach/.
* A testing of smaller set of bases, is sufficient for candidates smaller than various thresholds; <http://primes.utm.edu/prove/prove2_3.html>.
* <http://en.wikipedia.org/wiki/Miller-Rabin_primality_test>.
* <http://mathworld.wolfram.com/Rabin-MillerStrongPseudoprimeTest.html>
* <http://mathworld.wolfram.com/StrongPseudoprime.html>.
* <http://oeis.org/A014233>, <http://oeis.org/A006945>.
-}isPrimeByMillerRabin::Integrali=>i->BoolisPrimeByMillerRabinprimeCandidate=not$witnessesCompositenessprimeCandidate(fst$lastbinaryFactors--Odd-remainder.)(lengthbinaryFactors--The number of times that 'two' can be factored-out from 'predecessor'.)`any`testBaseswherepredecessor=primeCandidate-1binaryFactors=takeWhile((==0).snd).tail{-drop the original-}$iterate((`quotRem`2).fst)(predecessor,0)--Factor-out powers of two.testBases|nullfewestPrimeBases=letmillersTestSet=floor.(*2{-Eric Bach-}).Math.Power.square.toRational{-avoid premature rounding-}$log(fromIntegralprimeCandidate::Double{-overflows at 10^851-})in[2..predecessor`min`millersTestSet]|otherwise=headfewestPrimeBases`take`Data.Numbers.Primes.primeswherefewestPrimeBases=mapfst$dropWhile((primeCandidate>=).snd)[(0,9),--All odd integers less this, are prime, and require no further verification.(1,2047),(2,1373653),(3,25326001),(4,3215031751),(5,2152302898747),--Jaeschke ...(6,3474749660383),(8,341550071728321),(11,3825123056546413051),--Zhang ...(12,318665857834031151167461),(13,3317044064679887385961981),(14,6003094289670105800312596501),(15,59276361075595573263446330101),(17,564132928021909221014087501701),(19,1543267864443420616877677640751301),(20,10^(36::Int))--At least.]