-- |-- 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 #-}moduleMath.NumberTheory.Primes.Factorisation.Certified(certifiedFactorisation,certificateFactorisation,provenFactorisation)whereimportSystem.RandomimportControl.Monad.State.StrictimportControl.ApplicativeimportData.MaybeimportData.BitsimportMath.NumberTheory.Primes.Factorisation.MontgomeryimportMath.NumberTheory.Primes.Testing.Certificates.InternalimportMath.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=mapfst.certificateFactorisation-- | @'certificateFactorisation' n@ produces a 'provenFactorisation'-- with a default bound of @100000@.certificateFactorisation::Integer->[((Integer,Int),PrimalityProof)]certificateFactorisationn=provenFactorisation100000n-- | @'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=[]provenFactorisationbdn|n<2=error"provenFactorisation: argument not positive"|bd<2000=provenFactorisation2000n|otherwise=test$casesmallFactorsbdnof(sfs,mb)->map(\t@(p,_)->(t,smallCertp))sfs++casembofNothing->[]Justk->certiFactorisation(Just$bd*(bd+2))primeCheck(randomR.(,)6)(mkStdGen$fromIntegraln`xor`0xdeadbeef)Nothingk-- | verify that we indeed have a correct primality prooftest::[((Integer,Int),PrimalityProof)]->[((Integer,Int),PrimalityProof)]test(t@((p,_),prf):more)|p==cprimeprf&&checkPrimalityProofprf=t:testmore|otherwise=error(invalidpprf)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->MaybePrimalityProofprimeCheckn|bailliePSWn=casecertifyBPSWnofproof@Pocklington{}->Justproof_->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::MaybeInteger-- ^ Lower bound for composite divisors->(Integer->MaybePrimalityProof)-- ^ A primality test->(Integer->g->(Integer,g))-- ^ A PRNG->g-- ^ Initial PRNG state->MaybeInt-- ^ Estimated number of digits of the smallest prime factor->Integer-- ^ The number to factorise->[((Integer,Int),PrimalityProof)]-- ^ List of prime factors, exponents and primality proofscertiFactorisationprimeBoundprimeTestprngseedmbdigsn=caseptestnofJustproof->[((n,1),proof)]Nothing->evalState(factndigits)seedwheredigits=fromMaybe8mbdigsmult1xs=xsmultjxs=[((p,j*k),c)|((p,k),c)<-xs]vdbxs=[(p,2*e)|(p,e)<-xs]dbl(u,v)=(mult2u,vdbv)ptest=caseprimeBoundofJustbd->\k->ifk<=bdthen(Just$smallCertk)elseprimeTestkNothing->primeTestrndRk=state(\gen->prngkgen)factmdigs=dolet(b1,b2,ct)=findParmsdigs(pfs,cfs)<-repFactmb1b2ctifnullcfsthenreturnpfselsedonfs<-forMcfs$\(k,j)->multj<$>factk(ifnullpfsthendigs+4elsedigs)return(mergeAll$pfs:nfs)repFactmb1b2count|count<0=return([],[(m,1)])|otherwise=dos<-rndRmcasemontgomeryFactorisationmb1b2sofNothing->repFactmb1b2(count-1)Justd->dolet!cof=m`quot`dcasegcdcofdof1->do(dp,dc)<-caseptestdofJustproof->return([((d,1),proof)],[])Nothing->repFactdb1b2(count-1)(cp,cc)<-caseptestcofofJustproof->return([((cof,1),proof)],[])Nothing->repFactcofb1b2(count-1)return(mergedpcp,dc++cc)g->doletd'=d`quot`gc'=cof`quot`g(dp,dc)<-caseptestd'ofJustproof->return([((d',1),proof)],[])Nothing->repFactd'b1b2(count-1)(cp,cc)<-caseptestc'ofJustproof->return([((c',1),proof)],[])Nothing->repFactc'b1b2(count-1)(gp,gc)<-caseptestgofJustproof->return([((g,2),proof)],[])Nothing->dbl<$>repFactgb1b2(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)]mergexxs@(x@((p,e),c):xs)yys@(y@((q,d),_):ys)=casecomparepqofLT->x:mergexsyysEQ->((p,e+d),c):mergexsysGT->y:mergexxsysmerge[]ys=ysmergexs_=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]=xsmergeAll(xs:ys:zss)=merge(mergexsys)(mergeAllzss)-- | message for an invalid proof, should never be usedinvalid::Integer->PrimalityProof->Stringinvalidpprf=unlines["\nInvalid primality proof constructed, please report this to the package maintainer!","The supposed prime was:\n",showp,"\nThe presumed proof was:\n",showprf]