{-# LANGUAGE BangPatterns #-}------------------------------------------------------------------------------- |-- Module : Math.Sieve.Phi-- Copyright : (c) 2009-2012 Leon P Smith-- License : BSD3---- Maintainer : leon@melding-monads.com-- Stability : experimental-- Portability : portable---- This is a sieve that calculates Euler's Totient, also know as Euler's-- Phi function, for every number up to a given limit. @phi(n)@ is defined-- as the number of positive integers less than @n@ that are relatively-- prime to @n@, i.e. @gcd n x == 1@. Given the prime factorization of-- a number, it can be calculated efficiently from the formulas:---- > phi (p ^ k) = (p-1)*p^(k-1) if p is prime-- > phi (x * y) = phi x * phi y if x and y are relatively prime---- The second case says that @phi@ is a /multiplicative/ function. Since-- any multiplicative function can be calculated from the prime factorization,-- 'Math.Sieve.Factor' can also be applied, however this variant-- avoids a great deal of integer division, and so is significantly faster-- for calculating @phi@ for large quantities of different values.---- This sieve does not represent even numbers directly, and the maximum-- number that can currently be sieved is 2^32. This means that the sieve-- requires two bytes per number sieved on average, and thus sieving up to-- 2^32 requires 8 GB of storage.-------------------------------------------------------------------------------moduleMath.Sieve.Phi(sieve,phi,isPrime,sieveBound)whereimportControl.MonadimportControl.Monad.STimportData.Array.UnboxedimportData.Array.STimportData.BitsimportData.WordnewtypePhiSieve=PhiSieve(UArrayWord32Word32)instanceShowPhiSievewhereshowfs="<<PhiSieve "++show(shiftL(sieveBound'fs)1+2)++">>"instanceEqPhiSievewherea==b=sieveBound'a==sieveBound'binstanceOrdPhiSievewherecompareab=compare(sieveBound'a)(sieveBound'b)-- | The upper limit of the sievesieveBound::Integrala=>PhiSieve->asieveBound(PhiSievearr)=fromIntegral(shiftL(snd(boundsarr))1+2)sieveBound'::PhiSieve->Word32sieveBound'(PhiSievearr)=(snd(boundsarr))ix::Word32->Word32ixn=shiftRn1-- | Calculate Euler's Totient for every integer up to the given limitsieve::(Showa,Integrala)=>a->PhiSievesievelimit|not(0<=intlim&&intlim<=fromIntegral(maxBound::Word32))=error("Math.Sieve.Phi: out of bounds: "++showlimit)|otherwise=PhiSieve$runSTUArray$doa<-newArray(0,lim)1::STs(STUArraysWord32Word32)pows<-newArray(0,19)0::STs(STUArraysWord32Word32)sievea3powswhereintlim=(fromIntegrallimit)::Integerlim=ix(fromIntegrallimit-1)::Word32multapi=dox<-readArrayaiwriteArrayai(p*x)sieveappows|ixp>lim=returna|otherwise=dophi<-readArraya(ixp)when(phi==1)$doinitArraypows(p-1)adjustPhiap(ixp)powssievea(p+2)powsadjustPhi!a!p!i!pows|i>lim=return()|otherwise=dok<-zerospowsifk==0thenmulta(p-1)ielsemulta((p-1)*p^k)idecpowsppowsdecpowsppowsadjustPhiap(i+p)powsinitArray::STUArraysWord32Word32->Word32->STs()initArray!a!val=loopaval0whereloopavaln|n<=19=writeArrayanval>>loopaval(n+1)|otherwise=return()decpows::Word32->STUArraysWord32Word32->STs()decpows!p!pows=loopppows0whereloopppowsn=dox<-readArraypowsnifx==0thendowriteArraypowsn(p-1)loopppows(n+1)elsedowriteArraypowsn(x-1)zeros::STUArraysWord32Word32->STsWord32zeros!pows=looppows0wherelooppowsn=dox<-readArraypowsnifx==0thenlooppows(n+1)elsereturnn-- | Retrieves the totient from the sievephi::(Integrala,Integralb)=>PhiSieve->a->bphi(PhiSievea)nn|evenn0=loop0(shiftRn01)|otherwise=fromIntegral(a!(ixn0))wheren0=fromIntegralnnloop!tn|evenn=loop(t+1)(shiftRn1)|otherwise=fromIntegral(shiftL(a!(ixn))t)-- | Is the given number prime?isPrime::Integrala=>PhiSieve->a->BoolisPrime(PhiSievea)n|evenn=n==2|otherwise=a!(ixnn)==nn-1wherenn=fromIntegraln