diff --git a/AUTHORS b/AUTHORS new file mode 100644 --- /dev/null +++ b/AUTHORS @@ -0,0 +1,2 @@ +Melissa O'Neill +Leon P Smith diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,10 @@ +Copyright (c) 2009, Melding Monads +All rights reserved. + +Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. + * Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. + * Neither the name of Melding Monads nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/NumberSieves.cabal b/NumberSieves.cabal new file mode 100644 --- /dev/null +++ b/NumberSieves.cabal @@ -0,0 +1,18 @@ +Name: NumberSieves +Version: 0.0 +Description: This package includes the Sieve of O'Neill and two generalizations of the Sieve of Eratosthenes. The Sieve of O'Neill is a fully incremental primality sieve based on priority queues. The other two are array based, and are not incremental. One sieves the smallest prime factor, and is useful if you want to factor a large quantity of small numbers. The other sieves Euler's Totient, which is the number of positive integers relatively prime and less than a given number. +License: BSD3 +License-file: LICENSE +Author: Melissa O'Neill, Leon P Smith +Maintainer: Leon P Smith +Build-Type: Simple +Category: Number Theory, Algorithms +Synopsis: Number Theoretic Sieves: primes, factorization, and Euler's Totient +Cabal-Version: >=1.2 + +Library + Build-Depends: base >= 3 && < 5, + array >= 0 && < 1 + Exposed-Modules: NumberTheory.Sieve.ONeill + NumberTheory.Sieve.Factor + NumberTheory.Sieve.Phi diff --git a/NumberTheory/Sieve/Factor.hs b/NumberTheory/Sieve/Factor.hs new file mode 100644 --- /dev/null +++ b/NumberTheory/Sieve/Factor.hs @@ -0,0 +1,152 @@ +{-# LANGUAGE BangPatterns #-} + +----------------------------------------------------------------------------- +-- | +-- Module : NumberTheory.Sieve.Factor +-- Copyright : (c) Leon P Smith 2009 +-- License : BSD3 +-- +-- Maintainer : leon@melding-monads.com +-- Stability : experimental +-- Portability : portable +-- +-- This is an array-based generalization of the Sieve of Eratosthenes that +-- associates a prime divisor to each number in the sieve. This is useful +-- if you want to factor a large quantity of small numbers +-- +-- This code contains two simple optimizations: even numbers are not +-- represented in the array, reducing both time and space by 50%. +-- Secondly, the smallest prime divisor is sieved, and the prime numbers +-- are represented by "0" in the array instead of themselves. This allows +-- the divisors to be stored in half the number of bits, reducing space +-- consumption by another 50%. +-- +-- Currently, this sieve is limited to numbers less than 2^32, and consumes +-- one byte per number in the sieve on average. Thus if you want to find +-- the smallest divisor of every number up to 2^32, you will need 4 GB +-- of storage. +-- +----------------------------------------------------------------------------- + + +module NumberTheory.Sieve.Factor + ( FactorSieve() + , findFactor + , sieveBound + , isPrime + , factor + , sieve + ) where + +import Control.Monad +import Control.Monad.ST +import Data.Array.ST +import Data.Array.MArray +import Data.Array.Unboxed + +import Data.Bits +import Data.Word + +-- Note that if you want to sieve numbers beyond 2^32, you probably do not +-- want to simply change this type to (UArray Word64 Word32), as this would +-- result in 4 GB of wasted space. Increasing the limit in a sensible +-- fashion would require the use of multiple arrays with heterogenous types. + +newtype FactorSieve = FactorSieve (UArray Word32 Word16) + +instance Show FactorSieve where + show fs = "<>" + +instance Eq FactorSieve where + a == b = sieveBound' a == sieveBound' b + +instance Ord FactorSieve where + compare a b = compare (sieveBound' a) (sieveBound' b) + +-- | Returns the upper limit of a sieve. +sieveBound :: Integral a => FactorSieve -> a +sieveBound (FactorSieve arr) = fromIntegral (shiftL (snd (bounds arr)) 1) + +sieveBound' (FactorSieve arr) = (snd (bounds arr)) + +-- | Is a number prime? +isPrime :: Integral a => FactorSieve -> a -> Bool +isPrime (FactorSieve arr) n + | even n = n == 2 + | n <= 1 = False + | otherwise = arr ! ix (fromIntegral n) == 0 + +-- | Returns the smallest prime divisor of a given number in the sieve. +findFactor :: Integral a => FactorSieve -> a -> a +findFactor (FactorSieve arr) n + | even n = 2 + | d == 0 = n + | otherwise = fromIntegral d + where + d = arr ! ix (fromIntegral n) + +ix n = shiftR n 1 + +{-# INLINE for #-} +for min max step f = loop min + where + loop !x + | x <= max = f x >> loop (x + step) + | otherwise = return () + +flsqrt x = floor (sqrt (fromIntegral x)) + +-- | Finds the smallest prime divisor of every number up to a given limit. +sieve :: Integral a => a -> FactorSieve +sieve limit + | not (0 <= intlim && intlim < 2^32) + = error ("NumberTheory.Sieve.Factor.sieve: out of bounds: " ++ show limit) + | otherwise + = FactorSieve $ runSTUArray $ do + arr <- newArray (ix 3, lim) 0 + for 3 (flsqrt limit) 2 + (\i -> do + i' <- readArray arr (ix i) + when (i' == 0) + (for (ix (i*i)) lim i + (\i' -> do + i'' <- readArray arr i' + when (i'' == 0) + (writeArray arr i' (fromIntegral i))))) + return arr + where + intlim = fromIntegral limit :: Integer + lim = ix (fromIntegral limit - 1) :: Word32 + +-- FIXME: clean up the definition of 'factor', but maintain speed + +-- | Factors a number completely using a sieve, e.g. +-- +-- > factor (sieve 1000) 360 == [(2,3),(3,2),(5,1)] +factor :: Integral a => FactorSieve -> a -> [(a,a)] +factor (FactorSieve arr) n = start (fromIntegral n) + where + p x y = ((,) $! fromIntegral x) $! fromIntegral y + + start n + | even n = loop0 (shiftR n 1) 1 + | otherwise = loop1 n + + loop0 !n !i + | even n = loop0 (shiftR n 1) (i + 1) + | otherwise = p 2 i : loop1 n + + loop1 !n + | n == 1 = [] + | d == 0 = [p n 1] + | otherwise = loop1' (n `div` d) d 1 + where + d = fromIntegral (arr ! ix n) + + loop1' !n !d !i + | d' == 0 = if n == d then [p d (i+1)] else [p d i, p n 1] + | d' == d = loop1' (n `div` d ) d (i+1) + | otherwise = p d i : loop1' (n `div` d') d' 1 + where + d' = fromIntegral (arr ! ix n) + diff --git a/NumberTheory/Sieve/ONeill.hs b/NumberTheory/Sieve/ONeill.hs new file mode 100644 --- /dev/null +++ b/NumberTheory/Sieve/ONeill.hs @@ -0,0 +1,210 @@ +----------------------------------------------------------------------------- +-- | +-- Module : NumberTheory.Sieve.ONeill +-- Copyright : (c) 2006-2009 Melissa O'Neill +-- License : BSD3 +-- +-- Maintainer : leon@melding-monads.com +-- Stability : experimental +-- Portability : portable +-- +-- Incremental primality sieve based on priority queues, described +-- in the paper /The Genuine Sieve of Eratosthenes/ by Melissa O'Neill, +-- /Journal of Functional Programming/, 19(1), pp95-106, Jan 2009. +-- +-- +-- +-- Code is unchanged, other than packaging, from that written by +-- Melissa O'Neill. This version contains optimizations not described +-- in the paper, primarily improving memory consumption. +-- +-- +-- +----------------------------------------------------------------------------- + +-- Generate Primes using ideas from The Sieve of Eratosthenes +-- +-- This code is intended to be a faithful reproduction of the +-- Sieve of Eratosthenes, with the following change from the original +-- - The list of primes is infinite +-- (This change does have consequences for representing the number table +-- from which composites are "crossed out".) +-- +-- (c) 2006-2007 Melissa O'Neill. Code may be used freely so long as +-- this copyright message is retained and changed versions of the file +-- are clearly marked. + +module NumberTheory.Sieve.ONeill + ( primes + , sieve + , calcPrimes + , primesToNth + , primesToLimit + ) where + + +-- Priority Queues; this is essentially a copy-and-paste-job of +-- PriorityQ.hs. By putting the code here, we allow the Haskell +-- compiler to do some whole-program optimization. (Based on ML +-- code by L. Paulson in _ML for the Working Programmer_.) + +data PriorityQ k v = Lf + | Br {-# UNPACK #-} !k v !(PriorityQ k v) !(PriorityQ k v) + deriving (Eq, Ord, Read, Show) + +emptyPQ :: PriorityQ k v +emptyPQ = Lf + +isEmptyPQ :: PriorityQ k v -> Bool +isEmptyPQ Lf = True +isEmptyPQ _ = False + +minKeyValuePQ :: PriorityQ k v -> (k, v) +minKeyValuePQ (Br k v _ _) = (k,v) +minKeyValuePQ _ = error "Empty heap!" + +minKeyPQ :: PriorityQ k v -> k +minKeyPQ (Br k v _ _) = k +minKeyPQ _ = error "Empty heap!" + +minValuePQ :: PriorityQ k v -> v +minValuePQ (Br k v _ _) = v +minValuePQ _ = error "Empty heap!" + +insertPQ :: Ord k => k -> v -> PriorityQ k v -> PriorityQ k v +insertPQ wk wv (Br vk vv t1 t2) + | wk <= vk = Br wk wv (insertPQ vk vv t2) t1 + | otherwise = Br vk vv (insertPQ wk wv t2) t1 +insertPQ wk wv Lf = Br wk wv Lf Lf + +siftdown :: Ord k => k -> v -> PriorityQ k v -> PriorityQ k v -> PriorityQ k v +siftdown wk wv Lf _ = Br wk wv Lf Lf +siftdown wk wv (t @ (Br vk vv _ _)) Lf + | wk <= vk = Br wk wv t Lf + | otherwise = Br vk vv (Br wk wv Lf Lf) Lf +siftdown wk wv (t1 @ (Br vk1 vv1 p1 q1)) (t2 @ (Br vk2 vv2 p2 q2)) + | wk <= vk1 && wk <= vk2 = Br wk wv t1 t2 + | vk1 <= vk2 = Br vk1 vv1 (siftdown wk wv p1 q1) t2 + | otherwise = Br vk2 vv2 t1 (siftdown wk wv p2 q2) + +deleteMinAndInsertPQ :: Ord k => k -> v -> PriorityQ k v -> PriorityQ k v +deleteMinAndInsertPQ wk wv Lf = error "Empty PriorityQ" +deleteMinAndInsertPQ wk wv (Br _ _ t1 t2) = siftdown wk wv t1 t2 + +leftrem :: PriorityQ k v -> (k, v, PriorityQ k v) +leftrem (Br vk vv Lf Lf) = (vk, vv, Lf) +leftrem (Br vk vv t1 t2) = (wk, wv, Br vk vv t t2) where + (wk, wv, t) = leftrem t1 +leftrem _ = error "Empty heap!" + +deleteMinPQ :: Ord k => PriorityQ k v -> PriorityQ k v +deleteMinPQ (Br vk vv Lf _) = Lf +deleteMinPQ (Br vk vv t1 t2) = siftdown wk wv t2 t where + (wk,wv,t) = leftrem t1 +deleteMinPQ _ = error "Empty heap!" + +-- A hybrid of Priority Queues and regular queues. It allows a priority +-- queue to have a feeder queue, filled with items that come in an +-- increasing order. By keeping the feed for the queue separate, we +-- avoid needlessly filling an O(log n) data structure with data that +-- it won't need for a long time. + +type HybridQ k v = (PriorityQ k v, [(k,v)]) + +initHQ :: PriorityQ k v -> [(k,v)] -> HybridQ k v +initHQ pq feeder = (pq, feeder) + +insertHQ :: (Ord k) => k -> v -> HybridQ k v -> HybridQ k v +insertHQ k v (pq, q) = (insertPQ k v pq, q) + +deleteMinAndInsertHQ :: (Ord k) => k -> v -> HybridQ k v -> HybridQ k v +deleteMinAndInsertHQ k v (pq, q) = postRemoveHQ(deleteMinAndInsertPQ k v pq, q) + where + postRemoveHQ mq@(pq, []) = mq + postRemoveHQ mq@(pq, (qk,qv) : qs) + | qk < minKeyPQ pq = (insertPQ qk qv pq, qs) + | otherwise = mq + +minKeyHQ :: HybridQ k v -> k +minKeyHQ (pq, q) = minKeyPQ pq + +minKeyValueHQ :: HybridQ k v -> (k, v) +minKeyValueHQ (pq, q) = minKeyValuePQ pq + + +-- Finally, we have acceptable queues, now on to finding ourselves some +-- primes. + + +-- Here we use a wheel to generate all the number that are not multiples +-- of 2, 3, 5, and 7. We use some hard-coded data for that. + +{-# SPECIALIZE wheel :: [Int] #-} +{-# SPECIALIZE wheel :: [Integer] #-} +wheel :: Integral a => [a] +wheel = 2:4:2:4:6:2:6:4:2:4:6:6:2:6:4:2:6:4:6:8:4:2:4:2:4:8:6:4:6:2:4:6 + :2:6:6:4:2:4:6:2:6:4:2:4:2:10:2:10:wheel + +-- Now generate the primes using that wheel + +-- Sometimes, memoization isn't your friend. Maybe you don't actually want +-- to remember all the primes for the duration of your program and doing so +-- is just wasted space. For that situation, we provide calcPrimes which +-- calculates the infinite list of primes from scratch. + +{-# SPECIALIZE calcPrimes :: () -> [Int] #-} +{-# SPECIALIZE calcPrimes :: () -> [Integer] #-} +-- | An infinite stream of primes +calcPrimes :: Integral a => () -> [a] +calcPrimes () = 2 : 3 : 5 : 7 : sieve 11 wheel + +{-# SPECIALIZE primes :: [Int] #-} +{-# SPECIALIZE primes :: [Integer] #-} +-- | An infinite stream of primes +primes :: Integral a => [a] +primes = calcPrimes () + +{-# SPECIALIZE primesToNth :: Int -> [Integer] #-} +{-# SPECIALIZE primesToNth :: Int -> [Int] #-} +-- | Returns the first @n@ primes +primesToNth :: Integral a => Int -> [a] +primesToNth n = take n (calcPrimes ()) + +{-# SPECIALIZE primesToLimit :: Integer -> [Integer] #-} +{-# SPECIALIZE primesToLimit :: Int -> [Int] #-} +-- | Returns primes up to some limit +primesToLimit :: Integral a => a -> [a] +primesToLimit limit = takeWhile (< limit) (calcPrimes ()) + +-- This version of the sieve takes a wheel, not a list to be sieved. +-- primes1 and primes2 represent the same infinite list of, but they are +-- consumed at different speeds. By creating separate expressions, we +-- avoid retaining all the material between the two points. Sometimes +-- (when you care about space usage) memoization is not your friend. + +{-# SPECIALIZE sieve :: Int -> [Int] -> [Int] #-} +{-# SPECIALIZE sieve :: Integer -> [Integer] -> [Integer] #-} +-- | The first argument specifies which number to start at, +-- the second argument is a wheel of deltas for skipping composites. +-- For example, @primes@ could be defined as +-- +-- > 2 : 3 : sieve 5 (cycle [2,4]) +sieve :: Integral a => a -> [a] -> [a] +sieve n [] = [] +sieve n wheel@(d:ds) = n : (map (\(p,wheel) -> p) primes1) where + primes1 = sieve' (n+d) ds initialTable + primes2 = sieve' (n+d) ds initialTable + initialTable = initHQ (insertPQ (n*n) (n, wheel) emptyPQ) + (map (\(p,wheel) -> (p*p,(p,wheel))) primes2) + sieve' n [] table = [] + sieve' n wheel@(d:ds) table + | nextComposite <= n = sieve' (n+d) ds (adjust table) + | otherwise = (n,wheel) : sieve' (n+d) ds table + where + nextComposite = minKeyHQ table + adjust table + | m <= n = adjust (deleteMinAndInsertHQ m' (p, ds) table) + | otherwise = table + where + (m, (p, d:ds)) = minKeyValueHQ table + m' = m + p * d diff --git a/NumberTheory/Sieve/Phi.hs b/NumberTheory/Sieve/Phi.hs new file mode 100644 --- /dev/null +++ b/NumberTheory/Sieve/Phi.hs @@ -0,0 +1,155 @@ +{-# LANGUAGE BangPatterns #-} + +----------------------------------------------------------------------------- +-- | +-- Module : NumberTheory.Sieve.Phi +-- Copyright : (c) Leon P Smith 2009 +-- 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, +-- 'NumberTheory.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. +-- +----------------------------------------------------------------------------- + + +module NumberTheory.Sieve.Phi + ( sieve + , phi + , isPrime + , sieveBound + ) where + +import Control.Monad +import Control.Monad.ST +import Data.Array.Unboxed +import Data.Array.ST +import Data.Bits +import Data.Word + +newtype PhiSieve = PhiSieve (UArray Word32 Word32) + +instance Show PhiSieve where + show fs = "<>" + +instance Eq PhiSieve where + a == b = sieveBound' a == sieveBound' b + +instance Ord PhiSieve where + compare a b = compare (sieveBound' a) (sieveBound' b) + +-- | The upper limit of the sieve +sieveBound :: Integral a => PhiSieve -> a +sieveBound (PhiSieve arr) = fromIntegral (shiftL (snd (bounds arr)) 1 + 2) + +sieveBound' :: PhiSieve -> Word32 +sieveBound' (PhiSieve arr) = (snd (bounds arr)) + + +ix :: Word32 -> Word32 +ix n = shiftR n 1 + +-- | Calculate Euler's Totient for every integer up to the given limit +sieve :: Integral a => a -> PhiSieve +sieve limit + | not (0 <= intlim && intlim <= fromIntegral (maxBound :: Word32)) + = error ("Data.Numbers.Sieves.Phi: out of bounds: " ++ show limit) + | otherwise + = PhiSieve $ runSTUArray $ do + a <- newArray (0, lim) 1 :: ST s (STUArray s Word32 Word32) + pows <- newArray (0, 19) 0 :: ST s (STUArray s Word32 Word32) + sieve a 3 pows + where + intlim = (fromIntegral limit) :: Integer + lim = ix (fromIntegral limit - 1) :: Word32 + + mult a p i = do + x <- readArray a i + writeArray a i (p * x) + + sieve a p pows + | ix p > lim = return a + | otherwise = do + phi <- readArray a (ix p) + when (phi == 1) $ do + initArray pows (p-1) + adjustPhi a p (ix p) pows + sieve a (p+2) pows + + adjustPhi !a !p !i !pows + | i > lim = return () + | otherwise = do + k <- zeros pows + if k == 0 + then mult a (p-1) i + else mult a ((p-1)*p^k) i + decpows p pows + decpows p pows + adjustPhi a p (i + p) pows + +initArray :: STUArray s Word32 Word32 -> Word32 -> ST s () +initArray !a !val = loop a val 0 + where + loop a val n + | n <= 19 = writeArray a n val >> loop a val (n+1) + | otherwise = return () + +decpows :: Word32 -> STUArray s Word32 Word32 -> ST s () +decpows !p !pows = loop p pows 0 + where + loop p pows n = do + x <- readArray pows n + if x == 0 + then do + writeArray pows n (p-1) + loop p pows (n+1) + else do + writeArray pows n (x-1) + +zeros :: STUArray s Word32 Word32 -> ST s Word32 +zeros !pows = loop pows 0 + where + loop pows n = do + x <- readArray pows n + if x == 0 + then loop pows (n+1) + else return n + +-- | Retrieves the totient from the sieve +phi :: (Integral a, Integral b) => PhiSieve -> a -> b +phi (PhiSieve a) nn + | even n0 = loop 0 (shiftR n0 1) + | otherwise = fromIntegral (a ! (ix n0)) + where + n0 = fromIntegral nn + loop !t n + | even n = loop (t+1) (shiftR n 1) + | otherwise = fromIntegral (shiftL (a ! (ix n)) t) + +-- | Is the given number prime? +isPrime :: Integral a => PhiSieve -> a -> Bool +isPrime (PhiSieve a) n + | even n = n == 2 + | otherwise = a ! (ix nn) == nn - 1 + where + nn = fromIntegral n diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,3 @@ +#!/usr/bin/env runhaskell +import Distribution.Simple +main = defaultMain