{-# LANGUAGE CPP #-}
{-
Copyright (C) 2011 Dr. Alistair Ward

This program is free software: you can redistribute it and/or modify
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.
-}

module Factory.Math.Implementations.Primality(
-- * Types
-- ** Data-types
Algorithm(..)
-- * Functions
-- ** Predicates
--	isPrimeByAKS,
--	isPrimeByMillerRabin,
--	witnessesCompositeness
) where

import			Control.Arrow((&&&))
import qualified	Control.DeepSeq
import qualified	Data.Numbers.Primes
import qualified	Factory.Data.MonicPolynomial		as Data.MonicPolynomial
import qualified	Factory.Data.Polynomial			as Data.Polynomial
import qualified	Factory.Data.QuotientRing		as Data.QuotientRing
import qualified	Factory.Math.MultiplicativeOrder	as Math.MultiplicativeOrder
import qualified	Factory.Math.Power			as Math.Power
import qualified	Factory.Math.Primality			as Math.Primality
import qualified	Factory.Math.PrimeFactorisation		as Math.PrimeFactorisation
import qualified	ToolShed.Defaultable			as Defaultable

#if MIN_VERSION_parallel(3,0,0)
import qualified	Control.Parallel.Strategies
#endif

-- | The algorithms by which /primality/-testing has been implemented.
data Algorithm factorisationAlgorithm	=
AKS factorisationAlgorithm
| MillerRabin

instance Defaultable.Defaultable (Algorithm factorisationAlgorithm)	where
defaultValue	= MillerRabin

instance Math.PrimeFactorisation.Algorithm factorisationAlgorithm => Math.Primality.Algorithm (Algorithm factorisationAlgorithm)	where
isPrime _ 2	= True	--The only even prime.
isPrime algorithm candidate
| 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.
) . take 5 {-arbitrarily-} \$ Data.Numbers.Primes.primes	--Excludes even numbers, provided at least the 1st prime is tested.
)		= False
| otherwise	= (
case algorithm of
AKS factorisationAlgorithm	-> isPrimeByAKS factorisationAlgorithm
MillerRabin			-> 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.Algorithm factorisationAlgorithm, Control.DeepSeq.NFData i, Integral i) => factorisationAlgorithm -> i -> Bool
isPrimeByAKS factorisationAlgorithm n	= and [
not \$ Math.Power.isPerfectPower n,	--Step 1.
Math.Primality.areCoprime n `all` filter (/= n) [2 .. r],	--Step 3.
#if MIN_VERSION_parallel(3,0,0)
and \$ Control.Parallel.Strategies.parMap Control.Parallel.Strategies.rdeepseq	--Benefits from '+RTS -H100M', which reduces garbage-collections.
#else
all
#endif
(
\a	-> let
--			lhs, rhs :: Data.Polynomial.Polynomial i i
lhs	= Data.Polynomial.raiseModulo (Data.Polynomial.mkLinear 1 a) n {-power-} n {-modulus-}
rhs	= Data.Polynomial.mod' (Data.Polynomial.mkPolynomial [(1, n), (a, 0)]) n
in Data.QuotientRing.areCongruentModulo (
Data.MonicPolynomial.mkMonicPolynomial lhs
) (
Data.MonicPolynomial.mkMonicPolynomial rhs
) (
Data.MonicPolynomial.mkMonicPolynomial modulus
) -- Because all these polynomials are /monic/, one can establish /congruence/ using /integer/-division.
) [
1 .. floor . (* lg) . sqrt \$ fromIntegral r
] --Step 4; (x + a)^n ~ x^n + a mod (x^r - 1, n).
] where
lg :: Double
lg	= logBase 2 \$ fromIntegral n

--	r :: i
r	= fst . head . dropWhile (
(<= floor (Math.Power.square lg)) . snd
) . map (
id &&& Math.MultiplicativeOrder.multiplicativeOrder factorisationAlgorithm n
) \$ Math.Primality.areCoprime n `filter` [2 ..]	--Step 2.

--	modulus :: Data.Polynomial.Polynomial i i
modulus	= Data.Polynomial.mkPolynomial [(1, r), (negate 1, 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 :: Integral i
=> i	-- ^ Candidate integer.
-> i
-> Int
-> i	-- ^ Base.
-> Bool
witnessesCompositeness candidate oddRemainder nPowersOfTwo base	= all (
\$ ((`mod` candidate) . Math.Power.square) `iterate` Math.Power.raiseModulo base oddRemainder candidate	--Repeatedly modulo-square.
) [
(/= 1) . head,					--Check whether the zeroeth modulo-power is incongruent to one.
all (/= pred candidate) . take nPowersOfTwo	--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 :: Integral i => i -> Bool
isPrimeByMillerRabin primeCandidate	= not \$ witnessesCompositeness primeCandidate (
fst \$ last binaryFactors	--Odd-remainder.
) (
length binaryFactors	--The number of times that 'two' can be factored-out from 'predecessor'.
) `any` testBases	where
predecessor	= primeCandidate - 1
binaryFactors	= takeWhile ((== 0) . snd) . tail {-drop the original-} \$ iterate ((`quotRem` 2) . fst) (predecessor, 0)	--Factor-out powers of two.
testBases
| null fewestPrimeBases	= let
millersTestSet	= floor . (* 2 {-Eric Bach-}) . Math.Power.square . toRational {-avoid premature rounding-} \$ log (fromIntegral primeCandidate :: Double {-overflows at 10^851-})
in [2 .. predecessor `min` millersTestSet]
| otherwise		= head fewestPrimeBases `take` Data.Numbers.Primes.primes
where
fewestPrimeBases	= map fst \$ 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.
]