diff --git a/Changes b/Changes --- a/Changes +++ b/Changes @@ -1,2 +1,4 @@ -0.1.0: -First release +0.1.0.1: + Elaborate on overflow, work more on native Ints in Eratosthenes +0.1.0.0: + First release diff --git a/Math/NumberTheory/Primes/Sieve.hs b/Math/NumberTheory/Primes/Sieve.hs --- a/Math/NumberTheory/Primes/Sieve.hs +++ b/Math/NumberTheory/Primes/Sieve.hs @@ -16,7 +16,10 @@ -- where sieving is done, thus sieving primes up to @n@ requires -- @/O/(sqrt n/log n)@ space. module Math.NumberTheory.Primes.Sieve - ( primes + ( -- * Limitations + -- $limits + -- * Sieves and lists + primes , sieveFrom , PrimeSieve , primeSieve @@ -26,3 +29,33 @@ ) where import Math.NumberTheory.Primes.Sieve.Eratosthenes + +-- $limits +-- +-- There are three factors limiting the range of these sieves. +-- +-- (1) Memory +-- +-- (2) Overflow +-- +-- (3) The internal representation of the state +-- +-- An Eratosthenes type sieve needs to store the primes up to the square root of +-- the currently sieved region, thus requires @/O/(n\/log n)@ space.We store @16@ bytes +-- of information per prime, thus a Gigabyte of memory takes you to about @1.6*10^18@. +-- The @log@ doesn't change much in that range, so as a first approximation, doubling +-- the storage increases the sieve range by a factor of four. +-- +-- On a 64-bit system, this is (currently) the only limitation to be concerned with, but +-- with more than four Terabyte of memory, the fact that the internal representation +-- currently limits the sieve range to about @6.8*10^25@ could become relevant. +-- Overflow in array indexing doesn't become a concern before memory and internal +-- representation would allow to sieve past @10^37@. +-- +-- On a 32-bit system, the internal representation imposes no additional limits, +-- but overflow has to be reckoned with. On the one hand, the fact that arrays are +-- 'Int'-indexed restricts the size of the prime store, on the other hand, overflow +-- in calculating the indices to cross off multiples is possible before running out +-- of memory. The former limits the upper bound of the monolithic 'primeSieve' to +-- shortly above @8*10^9@, the latter limits the range of the segmented sieves to +-- about @1.7*10^18@. diff --git a/Math/NumberTheory/Primes/Sieve/Eratosthenes.hs b/Math/NumberTheory/Primes/Sieve/Eratosthenes.hs --- a/Math/NumberTheory/Primes/Sieve/Eratosthenes.hs +++ b/Math/NumberTheory/Primes/Sieve/Eratosthenes.hs @@ -111,8 +111,9 @@ ] -- | List of primes. --- Since the sieve uses unboxed arrays, overflow occurs at some point, --- but not before @10^6*'fromIntegral' ('maxBound' :: 'Int')@ (I forgot where exactly). +-- Since the sieve uses unboxed arrays, overflow occurs at some point. +-- On 64-bit systems, that point is beyond the memory limits, on +-- 32-bit systems, it is at about @1.7*10^18@. primes :: [Integer] primes = 2:3:5:concat [[vO + toPrim i | i <- [0 .. li], unsafeAt bs i] | PS vO bs <- psieveList, let (_,li) = bounds bs] @@ -123,7 +124,7 @@ psieveList :: [PrimeSieve] psieveList = makeSieves plim sqlim 0 0 cache where - plim = 4801 -- prime #647 + plim = 4801 -- prime #647, 644 of them to use sqlim = plim*plim cache = runSTUArray $ do sieve <- sieveTo 4801 @@ -135,12 +136,9 @@ if p then do let !i = indx .&. 7 - k :: Integer - k = fromIntegral (indx `shiftR` 3) - strt1 = (k*(30*k + fromIntegral (2*rho i)) - + fromIntegral (byte i)) `shiftL` 3 - + fromIntegral (idx i) - !strt = fromIntegral strt1 .&. 0xFFFFF + k = indx `shiftR` 3 + strt1 = (k*(30*k + 2*rho i) + byte i) `shiftL` 3 + fromIntegral (idx i) + !strt = fromIntegral (strt1 .&. 0xFFFFF) !skip = fromIntegral (strt1 `shiftR` 20) !ixes = fromIntegral indx `shiftL` 23 + strt `shiftL` 3 + fromIntegral i unsafeWrite new j skip @@ -185,18 +183,18 @@ then unsafeWrite cache pr (w-1) else do ixes <- unsafeRead cache (pr+1) - let !stj = ixes .&. 0x7FFFFF - !ixw = ixes `shiftR` 23 - !i = fromIntegral (ixw .&. 7) - !k = fromIntegral ixw - i + let !stj = fromIntegral ixes .&. 0x7FFFFF -- position of multiple and index of cofactor + !ixw = fromIntegral (ixes `shiftR` 23) -- prime data, up to 41 bits + !i = ixw .&. 7 + !k = ixw - i -- On 32-bits, k > 44717396 means overflow is possible in tick !o = i `shiftL` 3 - !j = fromIntegral stj .&. 7 - !s = fromIntegral stj `shiftR` 3 + !j = stj .&. 7 -- index of cofactor + !s = stj `shiftR` 3 -- index of first multiple to tick off (n, u) <- tick k o j s - let !skip = fromIntegral n `shiftR` 20 - !strt = fromIntegral n .&. 0xFFFFF + let !skip = fromIntegral (n `shiftR` 20) + !strt = fromIntegral (n .&. 0xFFFFF) unsafeWrite cache pr skip - unsafeWrite cache (pr+1) (ixes - stj + strt `shiftL` 3 + fromIntegral u) + unsafeWrite cache (pr+1) ((ixes .&. complement 0x7FFFFF) .|. strt `shiftL` 3 .|. fromIntegral u) treat (pr+2) tick stp off j ix | lastIndex < ix = return (ix - sieveBits, j) @@ -244,8 +242,8 @@ (bt,ix) = idxPr plim !start = 8*bt+ix+1 !nlim = plim+4800 - sieve <- sieveTo nlim - (_,hi) <- getBounds sieve + sieve <- sieveTo nlim -- Implement SieveFromTo for this, it's pretty wasteful when nlim isn't + (_,hi) <- getBounds sieve -- very small anymore more <- countFromToWd start hi sieve new <- unsafeNewArray_ (0,num+2*more) :: ST s (STUArray s Int CacheWord) let copy i @@ -269,7 +267,7 @@ strt1 = strt0 - offset !strt = fromIntegral strt1 .&. 0xFFFFF !skip = fromIntegral (strt1 `shiftR` 20) - !ixes = fromIntegral indx `shiftL` 23 + strt `shiftL` 3 + fromIntegral i + !ixes = fromIntegral indx `shiftL` 23 .|. strt `shiftL` 3 .|. fromIntegral i unsafeWrite new j skip unsafeWrite new (j+1) ixes fill (j+2) (indx+1) @@ -373,7 +371,7 @@ | otherwise = (strt1 - bitOff, i) !strt = fromIntegral strt2 .&. 0xFFFFF !skip = fromIntegral (strt2 `shiftR` 20) - !ixes = fromIntegral indx `shiftL` 23 + strt `shiftL` 3 + fromIntegral r2 + !ixes = fromIntegral indx `shiftL` 23 .|. strt `shiftL` 3 .|. fromIntegral r2 unsafeWrite new j skip unsafeWrite new (j+1) ixes fill (j+2) (indx+1) diff --git a/arithmoi.cabal b/arithmoi.cabal --- a/arithmoi.cabal +++ b/arithmoi.cabal @@ -1,5 +1,5 @@ name : arithmoi -version : 0.1.0.0 +version : 0.1.0.1 cabal-version : >= 1.6 author : Daniel Fischer copyright : (c) 2011 Daniel Fischer