-- |-- Module: Math.NumberTheory.Primes.Sieve.Misc-- Copyright: (c) 2011 Daniel Fischer-- Licence: MIT-- Maintainer: Daniel Fischer <daniel.is.fischer@googlemail.com>-- Stability: Provisional-- Portability: Non-portable (GHC extensions)--{-# LANGUAGE CPP, BangPatterns, ScopedTypeVariables, MonoLocalBinds, FlexibleContexts #-}#if __GLASGOW_HASKELL__ >= 700{-# OPTIONS_GHC -fspec-constr-count=8 #-}#endif{-# OPTIONS_HADDOCK hide #-}moduleMath.NumberTheory.Primes.Sieve.Misc(-- * TypesFactorSieve,TotientSieve,CarmichaelSieve-- * Functions-- ** Smallest prime factors,factorSieve,sieveFactor,fsBound,fsPrimeTest-- ** Totients,totientSieve,sieveTotient-- ** Carmichael,carmichaelSieve,sieveCarmichael)whereimportControl.Monad.STimportData.Array.Base(unsafeRead,unsafeWrite,unsafeAt)importData.Array.STimportData.Array.UnboxedimportControl.Monad(when)importData.BitsimportGHC.WordimportSystem.RandomimportMath.NumberTheory.Powers.Squares(integerSquareRoot')importMath.NumberTheory.Primes.Sieve.IndexingimportMath.NumberTheory.Primes.Factorisation.MontgomeryimportMath.NumberTheory.Primes.Factorisation.UtilsimportMath.NumberTheory.Utils-- | A compact store of smallest prime factors.dataFactorSieve=FS{-# UNPACK #-}!Word{-# UNPACK #-}!(UArrayIntWord16)-- | A compact store of totients.dataTotientSieve=TS{-# UNPACK #-}!Word{-# UNPACK #-}!(UArrayIntWord)-- | A compact store of values of the Carmichael function.dataCarmichaelSieve=CS{-# UNPACK #-}!Word{-# UNPACK #-}!(UArrayIntWord)-- | @'factorSieve' n@ creates a store of smallest prime factors of the numbers not exceeding @n@.-- If you need to factorise many smallish numbers, this can give a big speedup since it avoids-- many superfluous divisions. However, a too large sieve leads to a slowdown due to cache misses.-- The prime factors are stored as 'Word16' for compactness, so @n@ must be-- smaller than @2^32@.factorSieve::Integer->FactorSievefactorSievebound|4294967295<bound=error"factorSieve: overflow"|bound<8=FS7(array(0,2)[(0,0),(1,0),(2,0)])|otherwise=FSbndfSievewherebnd=fromIntegerboundibnd=fromInteger((bound-3)`quot`2)svbd=(fromInteger(integerSquareRoot'bound)-1)`quot`2fSieve=runSTUArray$dosieve<-newArray(0,ibnd)0::STs(STUArraysIntWord16)letsifti|i<svbd=dosp<-unsafeReadsieveiwhen(sp==0)(mark(2*i+3)(2*i*(i+3)+3))sift(i+1)|otherwise=returnsievemarkpj|j>ibnd=return()|otherwise=dosp<-unsafeReadsievejwhen(sp==0)(unsafeWritesievej$fromIntegralp)markp(j+p)sift0-- | @'fsBound' sieve@ tells the limit to which the sieve stores the smallest prime factors.fsBound::FactorSieve->WordfsBound(FSb_)=b-- | @'fsPrimeTest' sieve n@ checks in the sieve whether @n@ is prime. If @n@ is larger-- than the sieve can handle, an error is raised.fsPrimeTest::FactorSieve->Integer->BoolfsPrimeTestfs@(FSbndsve)n|n<0=fsPrimeTestfs(-n)|n<2=False|fromIntegern.&.(1::Int)==0=n==2|n<=fromIntegralbnd=sve`unsafeAt`(fromInteger(n`shiftR`1)-1)==0|otherwise=error"Out of bounds"-- | @'sieveFactor' fs n@ finds the prime factorisation of @n@ using the 'FactorSieve' @fs@.-- For negative @n@, a factor of @-1@ is included with multiplicity @1@.-- After stripping any present factors @2@, the remaining cofactor @c@ (if larger-- than @1@) is factorised with @fs@. This is most efficient of course if @c@ does not-- exceed the bound with which @fs@ was constructed. If it does, trial division is performed-- until either the cofactor falls below the bound or the sieve is exhausted. In the latter-- case, the elliptic curve method is used to finish the factorisation.sieveFactor::FactorSieve->Integer->[(Integer,Int)]sieveFactor(FSbndsve)=checkwherebound=fromIntegralbndcheck0=error"0 has no prime factorisation"check1=[]checkn|n<0=(-1,1):check(-n)|n<=bound=go2w(fromIntegraln)-- avoid expensive Integer ops if possible|fromIntegern.&.(1::Int)==1=sieveLoopn|otherwise=go2ngo2wn|n.&.1==1=intLoop((n-3)`shiftR`1)|otherwise=caseshiftToOddCountnof(k,m)->(2,k):ifm==1then[]elseintLoop((m-3)`shiftR`1)go2n=caseshiftToOddCountnof(k,m)->(2,k):ifm==1then[]elsesieveLoopmsieveLoopn|bound<n=tdLoopn(integerSquareRoot'n)0|otherwise=intLoop(fromIntegral(n`shiftR`1)-1)intLoop::Word->[(Integer,Int)]intLoop!n=caseunsafeAtsve(fromIntegraln)of0->[(2*fromIntegraln+3,1)]p->letp'=fromIntegralpincountLoopp'(n`quot`p'-1)1countLoop!p!i!c=caseunsafeAtsve(fromIntegrali)of0|p-3==2*i->[(fromIntegralp,c+1)]|otherwise->(fromIntegralp,c):(2*fromIntegrali+3,1):[]q|fromIntegralq==p->countLoopp(i`quot`p-1)(c+1)|otherwise->(fromIntegralp,c):intLoopilstIdx=snd(boundssve)tdLoopnsrix|lstIdx<ix=curven|sr<p=[(n,1)]|pix/=0=tdLoopnsr(ix+1)-- not a prime|otherwise=casesplitOffpnof(0,_)->tdLoopnsr(ix+1)(k,m)->(p,k):casemof1->[]j|j<=bound->intLoop(fromIntegral(j`shiftR`1)-1)|otherwise->tdLoopj(integerSquareRoot'j)(ix+1)wherep=toPrimixpix=unsafeAtsveixcurven=stdGenFactorisation(Just(bound*(bound+2)))(mkStdGen$fromIntegraln`xor`0xdecaf00d)Nothingn-- | @'totientSieve' n@ creates a store of the totients of the numbers not exceeding @n@.-- A 'TotientSieve' only stores values for numbers coprime to @30@ to reduce space usage.-- The maximal admissible value for @n@ is @'fromIntegral' ('maxBound' :: 'Word')@.totientSieve::Integer->TotientSievetotientSievebound|fromIntegral(maxBound::Word)<bound=error"totientSieve: overflow"|bound<8=TS7(array(0,0)[(0,6)])|otherwise=TSbnd(totSievebnd)wherebnd=fromIntegerbound-- | @'sieveTotient' ts n@ finds the totient @&#960;(n)@, i.e. the number of integers @k@ with-- @1 <= k <= n@ and @'gcd' n k == 1@, in other words, the order of the group of units-- in @&#8484;/(n)@, using the 'TotientSieve' @ts@.-- First, factors of @2, 3@ and @5@ are handled individually, if the remaining-- cofactor of @n@ is within the sieve range, its totient is looked up, otherwise-- the calculation falls back on factorisation, first by trial division and-- if necessary, elliptic curves.sieveTotient::TotientSieve->Integer->IntegersieveTotient(TSbndsve)=checkwherebound=fromIntegralbndcheckn|n<1=error"Totient only defined for positive numbers"|n==1=1|otherwise=go2ngo2n=caseshiftToOddCountnof(0,_)->go31n(k,m)->lettt=(shiftL1(k-1))inifm==1thenttelsego3ttmgo3!ttn=casesplitOff3nof(0,_)->go5ttn(k,m)->letnt=tt*(2*3^(k-1))inifm==1thenntelsego5ntmgo5!ttn=casesplitOff5nof(0,_)->sieveLoopttn(k,m)->letnt=tt*(4*5^(k-1))inifm==1thenntelsesieveLoopntmsieveLoop!ttn|bound<n=tdLoopttn(integerSquareRoot'n)0|otherwise=caseunsafeAtsve(toIdxn)ofnt->tt*fromIntegralntlstIdx=snd(boundssve)tdLoop!ttnsrix|lstIdx<ix=curvettn|sr<p'=tt*(n-1)-- n is a prime|pix/=p-1=tdLoopttnsr(ix+1)-- not a prime, next|otherwise=casesplitOffp'nof(0,_)->tdLoopttnsr(ix+1)(k,m)->letnt=tt*ppTotient(p',k)incasemof1->ntj|j<=bound->nt*fromIntegral(unsafeAtsve(toIdxj))|otherwise->tdLoopntj(integerSquareRoot'j)(ix+1)wherep=toPrimixp'=fromIntegralppix=unsafeAtsveixcurvettn=tt*totientFromCanonical(stdGenFactorisation(Just(bound*(bound+2)))(mkStdGen$fromIntegraln`xor`0xdecaf00d)Nothingn)-- | @'carmichaelSieve' n@ creates a store of values of the Carmichael function-- for numbers not exceeding @n@.-- Like a 'TotientSieve', a 'CarmichaelSieve' only stores values for numbers coprime to @30@-- to reduce space usage. The maximal admissible value for @n@ is @'fromIntegral' ('maxBound' :: 'Word')@.carmichaelSieve::Integer->CarmichaelSievecarmichaelSievebound|fromIntegral(maxBound::Word)<bound=error"carmichaelSieve: overflow"|bound<8=CS7(array(0,0)[(0,6)])|otherwise=CSbnd(carSievebnd)wherebnd=fromIntegerbound-- | @'sieveCarmichael' cs n@ finds the value of @&#955;(n)@ (or @&#968;(n)@), the smallest positive-- integer @e@ such that for all @a@ with @gcd a n == 1@ the congruence @a^e &#8801; 1 (mod n)@ holds,-- in other words, the (smallest) exponent of the group of units in @&#8484;/(n)@.-- The strategy is analogous to 'sieveTotient'.sieveCarmichael::CarmichaelSieve->Integer->IntegersieveCarmichael(CSbndsve)=checkwherebound=fromIntegralbndcheckn|n<1=error"Carmichael function only defined for positive numbers"|n==1=1|otherwise=go2ngo2n=caseshiftToOddCountnof(0,_)->go31n(k,m)->lettt=casekof1->12->2_->(shiftL1(k-2))inifm==1thenttelsego3ttmgo3!ttn=casesplitOff3nof(0,_)->go5ttn(k,1)|tt==1->2*3^(k-1)|otherwise->tt*3^(k-1)(k,m)|tt==1->go5(2*3^(k-1))m|otherwise->go5(tt*3^(k-1))mgo5!ttn=casesplitOff5nof(0,_)->sieveLoopttn(k,m)->lettt'=casefromIntegertt.&.(3::Int)of0->tt2->2*tt_->4*ttnt=tt'*5^(k-1)inifm==1thenntelsesieveLoopntmsieveLoop!ttn|bound<n=tdLoopttn(integerSquareRoot'n)0|otherwise=caseunsafeAtsve(toIdxn)ofnt->tt`lcm`fromIntegralntlstIdx=snd(boundssve)tdLoop!ttnsrix|lstIdx<ix=curvettn|sr<p'=tt`lcm`(n-1)-- n is a prime|pix/=p-1=tdLoopttnsr(ix+1)-- not a prime, next|otherwise=casesplitOffp'nof(0,_)->tdLoopttnsr(ix+1)(k,m)->letnt=(lcmtt(p'-1))*p'^(k-1)incasemof1->ntj|j<=bound->nt`lcm`fromIntegral(unsafeAtsve(toIdxj))|otherwise->tdLoopntj(integerSquareRoot'j)(ix+1)wherep=toPrimixp'=fromIntegralppix=unsafeAtsveixcurvettn=tt`lcm`carmichaelFromCanonical(stdGenFactorisation(Just(bound*(bound+2)))(mkStdGen$fromIntegraln`xor`0xdecaf00d)Nothingn)spfSieve::Word->STs(STUArraysIntWord)spfSievebound=dolet(octs,lidx)=idxPrbound!mxidx=8*octs+lidxmxval::Wordmxval=30*fromIntegralocts+fromIntegral(rholidx)!mxsve=integerSquareRoot'mxval(kr,r)=idxPrmxsve!svbd=8*kr+rar<-newArray(0,mxidx)0letstartki=8*(k*(30*k+2*rhoi)+bytei)+idxitickpstpoffjix|mxidx<ix=return()|otherwise=dos<-unsafeReadarixwhen(s==0)(unsafeWritearixp)tickpstpoff((j+1).&.7)(ix+stp*deltaj+tau(off+j))siftix|svbd<ix=returnar|otherwise=doe<-unsafeReadarixwhen(e==0)(doleti=ix.&.7k=ix`shiftR`3!off=i`shiftL`3!stp=ix-i!p=toPrimixtickpstpoffi(startki))sift(ix+1)sift0totSieve::Word->UArrayIntWordtotSievebound=runSTUArray$doar<-spfSievebound(_,lst)<-getBoundsarlettotix|lst<ix=returnar|otherwise=dop<-unsafeReadarixifp==0thenunsafeWritearix(toPrimix-1)elsedolet!n=toPrimix(tp,m)=unFactp(n`quot`p)casemof1->unsafeWritearixtp_->dotm<-unsafeReadar(toIdxm)unsafeWritearix(tp*tm)tot(ix+1)tot0carSieve::Word->UArrayIntWordcarSievebound=runSTUArray$doar<-spfSievebound(_,lst)<-getBoundsarletcarix|lst<ix=returnar|otherwise=dop<-unsafeReadarixifp==0thenunsafeWritearix(toPrimix-1)elsedolet!n=toPrimix(tp,m)=unFactp(n`quot`p)casemof1->unsafeWritearixtp_->dotm<-unsafeReadar(toIdxm)unsafeWritearix(lcmtptm)car(ix+1)car0-- Find the p-part of the totient of (p*m) and the cofactor-- of the p-power in m.{-# INLINE unFact #-}unFact::Word->Word->(Word,Word)unFactpm=go(p-1)mwherego!ttk=casek`quotRem`pof(q,0)->go(p*tt)q_->(tt,k)