diff --git a/Math/NumberTheory/Primes/Counting/HowPrimeCountingWorks.md b/Math/NumberTheory/Primes/Counting/HowPrimeCountingWorks.md new file mode 100644 --- /dev/null +++ b/Math/NumberTheory/Primes/Counting/HowPrimeCountingWorks.md @@ -0,0 +1,153 @@ +## Algorithm and Implementation Description + +### History and Formulas + +This prime counting implementation is based on the work by Adrien-Marie Legendre, who in about the year 1800 observed that the count of the primes to a limit, Pi(limit), can be determined by the inclusion-exclusion principle as noted at the following link: + +https://en.wikipedia.org/wiki/Prime-counting_function#Algorithms_for_evaluating_%CF%80(x) + +In analysis of this series, it can be seen that it is also an expression of what the Sieve of Eratosthenes does when extended to become the Wheel Sieve where no composite number is culled more than once since the products of the unique primes must always be unique. Note that the terms with the product of odd numbers of primes on the denominator are subtracted and with even primes are added: thus the inclusion/exclusion principle. This expression can then be expressed recursively using a "Phi" expression as follows: + +Phi(m, n) = Phi(m, n - 1) - Phi(floor(m / pn), n - 1) where + `m` is the current limit for this "branch" of the Phi function and + `n` is the one-based index of the nth prime so p1 is 2, p2 is 3, etc + and with the terminating conditions that + Phi(0, n) is 0 and Phi(m, 0) is m. + +Now the Legendre formula for the prime count, Pi(limit), is as follows: + +Pi(limit) = Phi(limit, Pi(sqrt(limit))) + Pi(sqrt(limit)) - 1 + +where from this it can be seen that one needs only know the primes to the square root of the range to be counted (and thus the number of them). + +The above recursive solution is quite inefficient as to complexity due to the deeply nested recursion, but a further terminating condition can stop the splitting of the "tree" a little sooner as follows: + +Phi(m, n) when m is less than or equal to pn must be exactly one + +This can be seen to be true because Phi(m, n) expresses the number of values up to `m` that cannot be divided evenly by any of the primes up to `pn` (this includes one, which can't be divided evenly by all numbers higher than it), thus excludes the sieving primes in the range p1 to pn. Therefore, the Phi value for any `m` <= `pn` must be just one. + +The computational complexity of this recursive expression is as follows: + +O(limit/((log limit)^2)), which while lower than by the Sieve of Eratosthenes due to the log squared in the denominator, is practically quite slow due to the high cost of the multiple recursions and (long) divisions. + +However, it has been realized quite recently (1985) that one does not have to solve that expression recursively by, instead of starting from the top of the "tree" and "splitting" first by the larger primes and building up the prime products down to the lowest primes, one can start at the bottom of the "tree" and divide by the lowest prime first on up; also, one can do all of the divisions by each of the primes in partial sieving "waves" while retaining the results in an array of intermediate results which replace the recursive method's "stack". One then finds that the results of all of the prime products already exist in the as yet unculled values in the partial sieved passes, greatly reducing the amount of computation (especially divisions) of having to build up these prime products from individual primes. + +"Phi's" as described are closely related to "rough" numbers, which are the values remaining after culling of all multiples of base primes up to some "partial sieving" base prime and including a count for the one value, so the number of rough numbers up to a limit sieved by primes up to pn is exactly Phi(that limit, n). Using an array of rough numbers and a parallel array of the Phi's (or Pi's) for each of those rough numbers is an efficient way to organize the intermediate calculations. + +Note that Phi's and Pi's are also closely related and can be converted one to the other because: + +Phi(x, Pi(sqrt(x))) = Pi(x) - Pi(sqrt(x)) + 1 + +Rough numbers: https://en.wikipedia.org/wiki/Rough_number.Applicative + +### Method + +In general, the method is as follows: + +1. All arrays are of size the square root of the counting range limit, as follows: + + a. A `roughs` array to process the remaining rough numbers in culling waves, initialized to all odds starting from 1 (ie. 1, 3, ...), + + b. A `phis` array which contains counts of remaining values after culling up to some level of base prime but is used for emulating the stack during processing initialized with the values as culled by two as in ((`limit`/roughs[`i`] + 1) `div` 2) where `i` is the zero-base index, and + + c. A `pisndxs` array that also contains a simple count of values remaining after culling up to a given culling base prime level but is in `Pi` form with the representations of the already culled base primes included in the count, and with no count for 1 just as for a normal pi; this is initialized with the index numbers as in (0, 1, ...). + + d. Variables are initialized to track the number of base primes processed initialized to zero as well as the current effective size of the `roughs`/`pis` arrays initialized to their current total size. + + + The `pisndxs` array and the `roughs` array are closely related to each other in that if one takes a given rough value at a given index, divides that rough value by two and uses the quotient as an index to this array will produce an element value that when the number of (odd) base primes already used is subtracted, produces the index in the roughs array of the original roughs element to form a connected loop; also, this array is used to look up "Pi" values that are used to add and subtract the counts in the manner as described for the phi/pi calculation. The implementation tries to reduce the number of offsets added and subtracted in the "hottest" loops. + +2. The main loop takes the following recursive "bottom-up" function and + flattens it by using the `pis` array to replace the stack used + by this function: + + For the following function `baseprms` is an array of odd primes plus non-prime one and values from three up to the square root limit from above. + + ```haskell + phip2 x = (x + 1) / 2 + + phi pmult pilmt = + let looppi pi acc = + if pi >= pilmt then acc else + p = unsafeAt baseprms pi + npmult = pmult * p + -- special termination condition for when p >= limit / npmult... + if p * npmult >= limit then acc + pilmit - pi + else looppi (pi + 1) + (acc + phip2(limit / npmult) - phi npmult pi) -- recursion + ``` + + Phi(limit (Pi(sqrt(limit))) is the above function called with: + phi 1 (lenth paseprms). + + Using the recursive function is not efficient because it doesn't have the asymptotic improvements to performance due to partial sieving. + + The loop is for each "base prime" in succession up to the square root of the square root of the `limit` in three phases: + + a. the given base prime value and all multiples of it are marked to be culled if they haven't already been eliminated in previous passes by something like the Sieve of Eratosthenes culling pass per base prime. + + b. this phase processes each of the remaining `roughs` values in turn for all remaining after being marked and splits the processing depending on whether the base prime value multiplied by the remaining rough element is: + + - less than or equal to the square root limit, in which case it compresses the `roughs` and `pis` array by subtracting the next `pis` element from the current one and writing it into an index that doesn't included the marked/skipped elements, or + + - if greater than it subtracts a value looked up by rough value index in the `pisndxs` array and again writes the result into the compression index as per above. + + - in all cases it moves the `roughs` to synchronizes with the pi's indices to move values to the left to fill the holes left by eliminating the marked rough values. + + c. the last phase is to adjust pisndxs values so they address the correct rough/pis values after the above compression of these arrays has taken place. + +3. The partial sieving loop only runs the base prime values up to the square root of the square root of the limit; after that point all of the roughs are prime (other than one) and therefore, the `limit` divided by each of the unique products of pairs of these roughs are added to the final answer from the loop with the usual procedure of looking up the values to be added from the `pisndx` array by using as an index: toIndex (limit / p1 / p2) but with the limitation that p1 can not be larger than (limit / p1 / p1) which satisfies the terminating condition that `p1` must be less than or equal to the quotient: + + `limit` `div` `pi` `div` `p1`. + +4. There are various offset calculations and compenstation to make sure allowance has be made for the various conversions between Phi and Pi. + +### Detailed Method + +In more detail, the "bottom up" partial sieving method is as follows; this implementeds the odds-only optimization so the arrays represent only the odd numbers: + +1. Keep track of the number of base primes that have been processed yet and the current effective number of rough numbers; INITIALIAZE `y` as the square root of the limit and three arrays each of the size `y`, with the first small counts array keeping track of the current number of possible primes up to the number represented by the index value so initialized with each element containing the index value, the second the roughs array containing the current level zero rough numbers as per the current culling passes (starts with nothing culled, so all odd values could be prime), and finally, the third containing the current Phi's/ (Pi's in this case) for each rough number in the rough numbers array (initially limit divided by the rough value corresponding to the index). + +2. For each base prime in order from the lowest to the limit^(1/4) (same as the square root of `y`) do the following: + + a. Do a partial sieve culling pass of the composites boolean array by the Sieve of Eratosthenes marking all multiples of the current base prime, also marking the value representing the current base prime as "non rough". The marked values will be those that can't be rough number. + + b. Update the Phi's table for each value that is still rough, while + eliminating and removing them from the array of roughs by moving them down the table (to the left, to lower indices). The adjustment is done differently depending whether the product of the base prime and the current rough is greater than the cube root of the limit (where limit divided by this product would be less than the square root of the limit and the counts can be just looked up in the counts table) or less than or equal to the square root in which case they can't be looked up in the counts table; however, the value in the counts table can be looked up in the counts table and that can be used to loop up the adjustment value in the phis table. All adjustments are subtracted from the current value of the phis table and, since all accumulated phi values will be subtracted from the first accumulation in the phi table, has the effect of subtracting for the first partial sieve pass, adding for the seconds, subtracting for the third, and so on, thus implementing the inclusion/exclusion principle. + + c. Adjust the counts table so the counts reflect the rough values culled in the last partial sieve culling pass. + + d. Repeat all above steps of the stop 2x steps for each base prime up to the limit where the rough values have been fully sieved, adjusting the number of base primes and the effective used sizes of the roughs and phis arrays for every pass. + +3. All the accumulated phis from the roughs higher than one can be subtracted from the accumulated phi for one to get an amalgamated answer which is still not complete, as per the last steps. At this point there are no further changes made to the counts, roughs, or phis arrays. + +4. The remaining base primes are now all higher than the the square root of `y` and the remaining roughs higher than these base primes are all prime so the following results can all be added (add for even number of primes in the denominator product - inclusion). Also, these products are all higher than `y` so the square root of the limit, so all counts can now be looked up in the counts table as follows: + + For each base prime from above the `y`^(1/2) to limit^(1/3) value multiplied by each value above the first to just below `y`` such that the product of the first prime squared and second primes is never higher than the limit (the terminal condition that Phi(0, n) is 0), divide the limit by this product and use the quotient to index the counts value to be added to the answer. + +5. The final answer can be obtained by doing any necessary conversions from Phi to Pi. + +There are some additional offset and compensation adjustments to make the array indexing and final reaults work out right, especially adjusting for whether particular count/phi values included the current base primes or not; these are negligible in computing cost. + +### Asymptotic Complexity + +It can easily be seen that the memory requirements are proportional to the square root of the limit; however, the computational complexity isn't quite so obvious as follows: + +1. All of the 3a steps represent a Sieve of Eratosthenes to the square root of the limit, so `y` times log log `y`); this is close enough to `y` and the operations are fast, so this is not a significant portion of the time used. + +2. All of the 3b steps take the number of base primes to the square root of `y` or `y`^(1/2)/(log (`y`^(1/2))) or 2*`y`^(1/2)/(log `y`) times the number of roughs which quickly reduces to about `y`/log`y` so the product is 2*`y`^(3/2)/((log `y`)^2) or 8*`limit`^(3/4)/((log `limit`)^2). For large limits, this is clearly a stronger term than the first. + +3. All of the 3c stops are the product of all the base primes to the square root of `y` as above times almost all of the counts values (of size `y`) so this is 2*`y`^(3/2)/(log `y`) or 4*`limit`^(3/4)/(log `limit`), which makes this a slightly stronger term than that of 2 above due to the loss of a log factor on the dividend; however, the operations of this step are faster than those of step 2 because step 2 has divisions for almost all operations, so these two may be effectively about the same until large\counting ranges when this term will be the strongest. + +4. The operations of step 4 would takes a long time at `limit`/((log `y`)^2) but are saved by the terminating condition when p >= `limit`/`p`/`p` so are about `y`^(4/3)/((log `y`)^2) or 4*`limit`^(2/3)/((log `limit`)^2), meaning that it is considerably less strong than the operations of step 2, although also slightly more complex. + +5. Thus, the time complexity is likely mostly controlled by step 2 for smaller counting ranges and by step 3 for large counting ranges meaning that the asymptotic time complexity is about that of 3, which lacks one log factor on the denominator due to having to process most of the count array values for each base prime pass. Thus the asymptotic complexity for large ranges can be expressed as O(limit^(3/4)/(log limit)), although the epirical complexity won't approach that for the ranges where this algorithm is usable due to memory use constraints. + +This implementation only deals with odd values (1, 3, 5 ...), thereby reducing memory storage requirements by a factor of two. + +### Conclusions + +In short, this is an effective implementation of the Legendre algorithm because any of the "Phi" values for the products of unique primes have already been computed simply and are available in the arrays where this isn't available for the fully recursive expression (for instance, the limit/(p1*p2*P3) value will be immediately available where the recursive expression would have to build it +up from all the combinations of limit/(p1*p2), limit/(p1*p3), and limit/(p2*p3), and so on for longer chains of base prime products). + +The worst problem with this algorithm and implementation is its relatively high memory use of O(limit^(1/2)) which in this case is eight times the square root of the counting range in bytes. This means that to sieve to 1e16 takes about 800 Megabytes and to the theoretical maximum limit of 2^19 - 1 would take about 32 Gigabytes, even if one were willing to wait the hours required to return a result. For such larger counting ranges, there are much better algorithms that aren't that much more difficult to implement, such as a variation on this with the Daniel Friedrich Ernst Meissel modification in about 1870 to change `y` to the cube root of the counting range instead of the square root to use memory as the cube root of the counting range and somewhat improve on the execution asymptotic complexity to proportial to the limit to the two thirds power. diff --git a/Math/NumberTheory/Primes/Counting/Impl.hs b/Math/NumberTheory/Primes/Counting/Impl.hs --- a/Math/NumberTheory/Primes/Counting/Impl.hs +++ b/Math/NumberTheory/Primes/Counting/Impl.hs @@ -4,12 +4,13 @@ -- Licence: MIT -- Maintainer: Daniel Fischer -- --- Number of primes not exceeding @n@, @π(n)@, and @n@-th prime. +-- Number of primes not exceeding @limit@, @π(limit)@, and @n@-th prime. -- {-# LANGUAGE BangPatterns #-} {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE ScopedTypeVariables #-} +{-# OPTIONS_GHC -O2 #-} {-# OPTIONS_GHC -fspec-constr-count=24 #-} {-# OPTIONS_GHC -Wno-incomplete-uni-patterns #-} @@ -20,54 +21,214 @@ ) where import Math.NumberTheory.Primes.Sieve.Eratosthenes - (PrimeSieve(..), primeList, primeSieve, psieveFrom, sieveTo, sieveBits, sieveRange) -import Math.NumberTheory.Primes.Sieve.Indexing (toPrim, idxPr) -import Math.NumberTheory.Primes.Counting.Approximate (nthPrimeApprox, approxPrimeCount) + (PrimeSieve(..), primeSieve, psieveFrom) +import Math.NumberTheory.Primes.Sieve.Indexing (toPrim) +import Math.NumberTheory.Primes.Counting.Approximate (nthPrimeApprox) import Math.NumberTheory.Primes.Types -import Math.NumberTheory.Roots +import Math.NumberTheory.Roots (integerSquareRoot) import Math.NumberTheory.Utils.FromIntegral -import Control.Monad.ST -import Data.Bit -import Data.Bits -import Data.Int +import Data.Word (Word64, Word32) +import Data.Bits (Bits(shiftR, (.&.), (.|.))) +import Control.Monad (forM_, when) +import Control.Monad.ST (ST, runST) +import Data.Array.Base (STUArray, MArray(unsafeNewArray_), + unsafeAt, unsafeFreezeSTUArray, unsafeRead, unsafeWrite) +import Data.Bit (Bit(..), unBit, nthBitIndex, countBits) import qualified Data.Vector.Unboxed as U -import qualified Data.Vector.Unboxed.Mutable as MU -- | Maximal allowed argument of 'primeCount'. Currently 8e18. primeCountMaxArg :: Integer primeCountMaxArg = 8000000000000000000 --- | @'primeCount' n == π(n)@ is the number of (positive) primes not exceeding @n@. +-- | @'primeCount' limit == π(limit)@ is the number of primes not exceeding @limit@. -- --- For efficiency, the calculations are done on 64-bit signed integers, therefore @n@ must --- not exceed 'primeCountMaxArg'. +-- For efficiency, the calculations are done on 64-bit unsigned integers, therefore @limit@ +-- must not exceed 'primeCountMaxArg'. -- --- Requires @/O/(n^0.5)@ space, the time complexity is roughly @/O/(n^0.7)@. --- For small bounds, @'primeCount' n@ simply counts the primes not exceeding @n@, --- for bounds from @30000@ on, Meissel's algorithm is used in the improved form due to --- D.H. Lehmer, cf. --- . +-- Requires @/O/(limit^0.5)@ space, the time complexity is roughly @/O/(limit^0.7)@. +-- @'primeCount' limit@ uses Legendre's algorithm in an improved form using "partial sieving" +-- and processing by "splitting" based on whether the product of the "base prime" and a +-- product of higher co-primes is less than or equal to, or greater than the square root of +-- @limit@ where "base primes" are lower than the square root of the square root of @limit@; +-- above this limit all remaining values are primes so are used in unique prime pairs to +-- calculate the additional amounts to add to the accumulated "Phi" for the answer. + +-- NOTE: This is not related to the later work (about 1870) by Daniel Friedrich Ernst Meissel, +-- nor the extension to Meissel's work by Professor D. H. Lehmer in 1959 to adapt these types +-- of algorithms to use on a mainframe computer of that time; both of whose purpose was to +-- reduce the number of calculations for a given counting range and/or reduce the storage +-- requirements. Neither of these used "partial sieving" to get anywhere near the asymptotic +-- complexity of this implementation. + +-- See the "HowPrimeCountingWorks.md" file in this directory for a more complete explanation +-- of how this implementation works. + primeCount :: Integer -> Integer -primeCount n - | n > primeCountMaxArg = error $ "primeCount: can't handle bound " ++ show n - | n < 2 = 0 - | n < 1000 = intToInteger . length . takeWhile (<= n) . map unPrime . primeList . primeSieve $ max 242 n - | n < 30000 = runST $ do - ba <- sieveTo n - let (s, e) = (0, MU.length ba - 1) - ct <- countFromTo s e ba - return (intToInteger $ ct+3) - | otherwise = - let !ub = cop $ fromInteger n - !sr = integerSquareRoot ub - !cr = nxtEnd $ integerCubeRoot ub + 15 - nxtEnd k = k - (k `rem` 30) + 31 - !phn1 = calc ub cr - !cs = cr+6 - !pdf = sieveCount ub cs sr - in phn1 - pdf +primeCount limit + | limit > primeCountMaxArg = error $ "primeCount: can't handle bound " ++ show limit + | limit < 9 = if limit < 2 then 0 else (limit + 1) `div` 2 + | otherwise = + let + + -- initialize constants... + ilimit = fromIntegral limit + sqrtlmt = fromIntegral $ integerSquareRoot ilimit + sqrtsqrtlmt = integerSquareRoot sqrtlmt + maxndx = toIndex sqrtlmt -- last index of arrays + -- `numbps` will be the number of odd base primes up to `limit`^(1/4); + -- `roughssz` is the current effective length of `roughs` after reduction; + -- `phindxs` is the count of odd primes to limit by index not including bps, + -- the above is an index to the `roughs`/`phis` that repr the index; + -- `roughs` are the values remaining after culling base primes < its index, + -- the above values when divided by two is the index for `phindxs`; + -- `phis` is count of odd primes to a limit set by its index... + (numbps, roughssz, phindxs, roughs, phis) = runST $ do -- run in ST monad... + + -- initialize monadic versions of `phindxs`/`roughs`/`phis`... + mphindxs <- unsafeNewArray_ (0, maxndx) :: ST s (STUArray s Int Word32) + -- initially is odd phi of the index repr 3 shr 1 is 1 -> value 1, etc... + forM_ [ 0 .. maxndx ] $ \ i -> unsafeWrite mphindxs i (fromIntegral i) + unsafeWrite mphindxs 0 1 -- for correctness, never used! + mrs <- unsafeNewArray_ (0, maxndx) :: ST s (STUArray s Int Word32) + -- odd values from 1 as in 1, 3 ... to maximum for roughs... + forM_ [ 0 .. maxndx ] $ \ i -> + unsafeWrite mrs i (fromIntegral i * 2 + 1) + mphis <- unsafeNewArray_ (0, maxndx) :: ST s (STUArray s Int Word64) + -- initialized to limit // roughs; are phis including count for one... + forM_ [ 0 .. maxndx ] $ \ i -> do + r <- unsafeRead mrs i; let d = fromIntegral r + unsafeWrite mphis i (fromIntegral $ phip2 $ divide ilimit d) + + -- all work requiring modifying arrays and values done recursively here; + -- the "partial sieving" loop with one `bp` sieving pass per loop; + -- uses `roi` output and `rii` input roughs processing indices... + let loop !nbps !rsilmt = do -- `rslmt` is maximum current ndx for `roughs` + bpw32 <- unsafeRead mrs 1 + let bp = fromIntegral bpw32 + if bp > sqrtsqrtlmt then do -- means `bp` primes <= limit^(1/4) + fmsops <- unsafeFreezeSTUArray mphindxs + fmrs <- unsafeFreezeSTUArray mrs + fmlops <- unsafeFreezeSTUArray mphis + return (nbps, rsilmt + 1, fmsops, fmrs, fmlops) -- done loop! + else do -- for each base prime `bp`... + let -- mark `mrs` values that are multiples of bp if still there + cullmrs cullpos = + if cullpos > sqrtlmt then return () else do + cnt <- unsafeRead mphindxs (cullpos `shiftR` 1) + let ndx = fromIntegral cnt - nbps + tstr <- unsafeRead mrs ndx + when (tstr == fromIntegral cullpos) + (unsafeWrite mrs ndx 0) + cullmrs (cullpos + bp + bp) + + -- recursive function to process all remaining `mrs` by + -- forming products of unique pairs with `bp`... + split rii !roi = + if rii > rsilmt then return (roi - 1) else do + m <- unsafeRead mrs rii -- multiplier may not be prime! + if m == 0 then split (rii + 1) roi else do -- skip marked + -- only unculled values; may not be prime + olv <- unsafeRead mphis rii -- large odd "pi" to adjust + let mbp = fromIntegral m * fromIntegral bp + adjv <- -- depends on condition... + if mbp <= fromIntegral sqrtlmt then do + -- ilimit `div` mbp too large... + let cnti = fromIntegral mbp `shiftR` 1 + adji <- unsafeRead mphindxs cnti + let adjndx = fromIntegral adji - nbps + unsafeRead mphis adjndx +-- adj <- unsafeRead mphis adjndx +-- return $ adj - fromIntegral nbps + else do + -- ilimit `div` mbp in index range; use directly! + let adjndx = toIndex (divide ilimit mbp) + adj <- unsafeRead mphindxs adjndx -- phi form... + return $ fromIntegral adj - fromIntegral nbps + 1 + -- write adjusted value into `mphis` at new offset... + unsafeWrite mphis roi (olv - adjv) + unsafeWrite mrs roi m -- move rougn values in sync + split (rii + 1) (roi + 1) -- recursively loop + + -- update `mphindxs` array for last cull pass... + adjcnt cm !mxci = -- cull multiple and maximum index + if cm < bp then return () else do + ofstc <- unsafeRead mphindxs (cm `shiftR` 1) + let c = ofstc - fromIntegral nbps + e = (cm * bp) `shiftR` 1 + adjci ci = + if ci < e then adjcnt (cm - 2) ci else do + ov <- unsafeRead mphindxs ci + unsafeWrite mphindxs ci (ov - c) + adjci (ci - 1) + adjci mxci + + -- the code that uses the above "let"'s... + unsafeWrite mrs 1 0 -- mark first non-one rough for deletion + cullmrs (bp * bp) -- mark the other rough multiples of `bp` + maxrsi <- split 0 0 -- adjust `roughs` and `phis` for cull + let topcullpnt = (sqrtlmt `div` bp - 1) .|. 1 -- odd <= sqrtlmt + adjcnt topcullpnt maxndx -- update `phindxs` for culling pass + loop (nbps + 1) maxrsi -- recurse for all base primes + + loop 0 maxndx -- calling recursive "partial sieving" loop! + + -- the offset of the sum of the other `phis`... + othroddpis = sum [ unsafeAt phis bpi | bpi <- [ 1 .. roughssz - 1 ] ] + -- subtracted from the first first element of `phis`... + phi0 = unsafeAt phis 0 - othroddpis -- + pi0crct -- to produce an intermediate phi + + -- recursively calculate the additional odd "phis" for all pairs of + -- unique primes/`roughs` starting above limit^(1/4) to limit^(1/2); + -- these are exactly the remaining values in `roughs` above the "one"; + -- Note that all roughs above the first "one" element are now prime; + -- pre-comp of all additional ones for following "pairs" calculation... + phi0adj = fromIntegral $ (roughssz - 2) * (roughssz - 1) `div` 2 + accum p1i !ans = + if p1i >= roughssz - 1 then ans else -- for all roughs skipping "one"... + let p1 = fromIntegral $ unsafeAt roughs p1i -- `p1` - first of prime pair + qp1 = ilimit `div` p1 -- pre divide for "p1" value + ndx = unsafeAt phindxs (toIndex (fromIntegral (qp1 `div` p1))) + endndx = fromIntegral ndx - numbps -- last `p1` index! + adj = fromIntegral $ (endndx - p1i) * (numbps + p1i - 1) + comp p2i !ac = -- `rii` is index of the second of the `roughs` + if p2i > endndx then accum (p1i + 1) ac else -- exit if reach end index! + let p2 = fromIntegral (unsafeAt roughs p2i) -- second rough + cnti = toIndex (divide qp1 p2) -- ndx for comp "pi" + in comp (p2i + 1) (ac + fromIntegral (unsafeAt phindxs cnti)) + in if endndx <= p1i then ans -- terminate if `p1`^3 >= `ilimit`! + -- adjust for ones added and not used due to `endndx` termination... + else comp (p1i + 1) (ans - adj) + + numsqrtprms = fromIntegral $ numbps + roughssz + + -- finally call addition of the phis for all pairs of unique primes added to + -- the offset of the other odd "phis" already adjusted for the subtracting + -- of odd "phis" from odd "phis" plus the total primes to the square root of + -- the `limit` counting range minus one according to the formula... + in fromIntegral $ accum 1 (phi0 + phi0adj) + numsqrtprms - 1 + +-------------------------------------------------------------------------------- +-- Auxiliaries -- +-------------------------------------------------------------------------------- + +{-# INLINE divide #-} +divide :: Word64 -> Word64 -> Int +divide n d = fromIntegral $ n `div` d + +{-# INLINE phip2 #-} +phip2 :: Int -> Int +phip2 x = (x + 1) `shiftR` 1 + +{-# INLINE toIndex #-} +toIndex :: Int -> Int +toIndex x = (x - 1) `shiftR` 1 + +-------------------------------------------------------------------------------- +-- Nth Prime -- +-------------------------------------------------------------------------------- + -- | @'nthPrime' n@ calculates the @n@-th prime. Numbering of primes is -- @1@-based, so @'nthPrime' 1 == 2@. -- @@ -144,271 +305,6 @@ r1 = r0 `quot` 3 r2 = min 7 (if r1 > 5 then r1-1 else r1) --- highSieve :: Integer -> Integer -> Integer -> Integer --- highSieve a surp gap = error "Oh shit" - -sieveCount :: Int64 -> Int64 -> Int64 -> Integer -sieveCount ub cr sr = runST (sieveCountST ub cr sr) - -sieveCountST :: forall s. Int64 -> Int64 -> Int64 -> ST s Integer -sieveCountST ub cr sr = do - let psieves = psieveFrom (int64ToInteger cr) - pisr = approxPrimeCount sr - picr = approxPrimeCount cr - diff = pisr - picr - size = int64ToInt (diff + diff `quot` 50) + 30 - store <- MU.unsafeNew size :: ST s (MU.MVector s Int64) - let feed :: Int64 -> Int -> Int -> U.Vector Bit -> [PrimeSieve] -> ST s Integer - feed voff !wi !ri uar sves - | ri == sieveBits = case sves of - (PS vO ba : more) -> feed (fromInteger vO) wi 0 ba more - _ -> error "prime stream ended prematurely" - | pval > sr = do - stu <- U.unsafeThaw uar - eat 0 0 voff (wi-1) ri stu sves - | unBit (uar `U.unsafeIndex` ri) = do - MU.unsafeWrite store wi (ub `quot` pval) - feed voff (wi+1) (ri+1) uar sves - | otherwise = feed voff wi (ri+1) uar sves - where - pval = voff + toPrim ri - eat :: Integer -> Integer -> Int64 -> Int -> Int -> MU.MVector s Bit -> [PrimeSieve] -> ST s Integer - eat !acc !btw voff !wi !si stu sves - | si == sieveBits = - case sves of - [] -> error "Premature end of prime stream" - (PS vO ba : more) -> do - nstu <- U.unsafeThaw ba - eat acc btw (fromInteger vO) wi 0 nstu more - | wi < 0 = return acc - | otherwise = do - qb <- MU.unsafeRead store wi - let dist = qb - voff - 7 - if dist < intToInt64 sieveRange - then do - let (b,j) = idxPr (dist+7) - !li = (b `shiftL` 3) .|. j - new <- if li < si then return 0 else countFromTo si li stu - let nbtw = btw + intToInteger new + 1 - eat (acc+nbtw) nbtw voff (wi-1) (li+1) stu sves - else do - let (cpl,fds) = dist `quotRem` intToInt64 sieveRange - (b,j) = idxPr (fds+7) - !li = (b `shiftL` 3) .|. j - ctLoop !lac 0 (PS vO ba : more) = do - nstu <- U.unsafeThaw ba - new <- countFromTo 0 li nstu - let nbtw = btw + lac + 1 + intToInteger new - eat (acc+nbtw) nbtw (integerToInt64 vO) (wi-1) (li+1) nstu more - ctLoop lac s (ps : more) = do - let !new = countAll ps - ctLoop (lac + intToInteger new) (s-1) more - ctLoop _ _ [] = error "Primes ended" - new <- countFromTo si (sieveBits-1) stu - ctLoop (intToInteger new) (cpl-1) sves - case psieves of - (PS vO ba : more) -> feed (fromInteger vO) 0 0 ba more - _ -> error "No primes sieved" - -calc :: Int64 -> Int64 -> Integer -calc lim plim = runST (calcST lim plim) - -calcST :: forall s. Int64 -> Int64 -> ST s Integer -calcST lim plim = do - !parr <- sieveTo (int64ToInteger plim) - let (plo, phi) = (0, MU.length parr - 1) - !pct <- countFromTo plo phi parr - !ar1 <- MU.unsafeNew end - MU.unsafeWrite ar1 0 lim - MU.unsafeWrite ar1 1 1 - !ar2 <- MU.unsafeNew end - let go :: Int -> Int -> MU.MVector s Int64 -> MU.MVector s Int64 -> ST s Integer - go cap pix old new - | pix == 2 = coll cap old - | otherwise = do - Bit isp <- MU.unsafeRead parr pix - if isp - then do - let !n = fromInteger (toPrim pix) - !ncap <- treat cap n old new - go ncap (pix-1) new old - else go cap (pix-1) old new - coll :: Int -> MU.MVector s Int64 -> ST s Integer - coll stop ar = - let cgo !acc i - | i < stop = do - !k <- MU.unsafeRead ar i - !v <- MU.unsafeRead ar (i+1) - cgo (acc + int64ToInteger v*cp6 k) (i+2) - | otherwise = return (acc+intToInteger pct+2) - in cgo 0 0 - go 2 start ar1 ar2 - where - (bt,ri) = idxPr plim - !start = 8*bt + ri - !size = int64ToInt $ integerSquareRoot lim `quot` 4 - !end = 2*size - -treat :: Int -> Int64 -> MU.MVector s Int64 -> MU.MVector s Int64 -> ST s Int -treat end n old new = do - qi0 <- locate n 0 (end `quot` 2 - 1) old - let collect stop !acc ix - | ix < end = do - !k <- MU.unsafeRead old ix - if k < stop - then do - v <- MU.unsafeRead old (ix+1) - collect stop (acc-v) (ix+2) - else return (acc,ix) - | otherwise = return (acc,ix) - goTreat !wi !ci qi - | qi < end = do - !key <- MU.unsafeRead old qi - !val <- MU.unsafeRead old (qi+1) - let !q0 = key `quot` n - !r0 = int64ToInt (q0 `rem` 30030) - !nkey = q0 - int8ToInt64 (cpDfAr `U.unsafeIndex` r0) - nk0 = q0 + int8ToInt64 (cpGpAr `U.unsafeIndex` (r0+1) + 1) - !nlim = n*nk0 - (wi1,ci1) <- copyTo end nkey old ci new wi - ckey <- MU.unsafeRead old ci1 - (!acc, !ci2) <- if ckey == nkey - then do - !ov <- MU.unsafeRead old (ci1+1) - return (ov-val,ci1+2) - else return (-val,ci1) - (!tot, !nqi) <- collect nlim acc (qi+2) - MU.unsafeWrite new wi1 nkey - MU.unsafeWrite new (wi1+1) tot - goTreat (wi1+2) ci2 nqi - | otherwise = copyRem end old ci new wi - goTreat 0 0 qi0 - --------------------------------------------------------------------------------- --- Auxiliaries -- --------------------------------------------------------------------------------- - -locate :: Int64 -> Int -> Int -> MU.MVector s Int64 -> ST s Int -locate p low high arr = do - let go lo hi - | lo < hi = do - let !md = (lo+hi) `quot` 2 - v <- MU.unsafeRead arr (2*md) - case compare p v of - LT -> go lo md - EQ -> return (2*md) - GT -> go (md+1) hi - | otherwise = return (2*lo) - go low high - -{-# INLINE copyTo #-} -copyTo :: Int -> Int64 -> MU.MVector s Int64 -> Int - -> MU.MVector s Int64 -> Int -> ST s (Int,Int) -copyTo end lim old oi new ni = do - let go ri wi - | ri < end = do - ok <- MU.unsafeRead old ri - if ok < lim - then do - !ov <- MU.unsafeRead old (ri+1) - MU.unsafeWrite new wi ok - MU.unsafeWrite new (wi+1) ov - go (ri+2) (wi+2) - else return (wi,ri) - | otherwise = return (wi,ri) - go oi ni - -{-# INLINE copyRem #-} -copyRem :: Int -> MU.MVector s Int64 -> Int -> MU.MVector s Int64 -> Int -> ST s Int -copyRem end old oi new ni = do - let len = end - oi - MU.copy (MU.slice ni len new) (MU.slice oi len old) - pure $ ni + len - -{-# INLINE cp6 #-} -cp6 :: Int64 -> Integer -cp6 k = - case k `quotRem` 30030 of - (q,r) -> 5760*int64ToInteger q + - int16ToInteger (cpCtAr `U.unsafeIndex` int64ToInt r) - -cop :: Int64 -> Int64 -cop m = m - int8ToInt64 (cpDfAr `U.unsafeIndex` int64ToInt (m `rem` 30030)) - - --------------------------------------------------------------------------------- --- Ugly helper arrays -- --------------------------------------------------------------------------------- - -cpCtAr :: U.Vector Int16 -cpCtAr = runST $ do - ar <- MU.replicate 30030 1 - let zilch s i - | i < 30030 = MU.unsafeWrite ar i 0 >> zilch s (i+s) - | otherwise = return () - accumulate ct i - | i < 30030 = do - v <- MU.unsafeRead ar i - let !ct' = ct+v - MU.unsafeWrite ar i ct' - accumulate ct' (i+1) - | otherwise = return ar - zilch 2 0 - zilch 6 3 - zilch 10 5 - zilch 14 7 - zilch 22 11 - zilch 26 13 - vec <- accumulate 1 2 - U.unsafeFreeze vec - -cpDfAr :: U.Vector Int8 -cpDfAr = runST $ do - ar <- MU.replicate 30030 0 - let note s i - | i < 30029 = MU.unsafeWrite ar i 1 >> note s (i+s) - | otherwise = return () - accumulate d i - | i < 30029 = do - v <- MU.unsafeRead ar i - if v == 0 - then accumulate 2 (i+2) - else do MU.unsafeWrite ar i d - accumulate (d+1) (i+1) - | otherwise = return ar - note 2 0 - note 6 3 - note 10 5 - note 14 7 - note 22 11 - note 26 13 - vec <- accumulate 2 3 - U.unsafeFreeze vec - -cpGpAr :: U.Vector Int8 -cpGpAr = runST $ do - ar <- MU.replicate 30031 0 - MU.unsafeWrite ar 30030 1 - let note s i - | i < 30029 = MU.unsafeWrite ar i 1 >> note s (i+s) - | otherwise = return () - accumulate d i - | i < 1 = return ar - | otherwise = do - v <- MU.unsafeRead ar i - if v == 0 - then accumulate 2 (i-2) - else do MU.unsafeWrite ar i d - accumulate (d+1) (i-1) - note 2 0 - note 6 3 - note 10 5 - note 14 7 - note 22 11 - note 26 13 - vec <- accumulate 2 30027 - U.unsafeFreeze vec - ------------------------------------------------------------------------------- -- Prime counting @@ -419,13 +315,3 @@ countToNth !n (PS v0 bs : more) = case nthBitIndex (Bit True) n bs of Just i -> v0 + toPrim i Nothing -> countToNth (n - countBits bs) more - --- count all set bits in a chunk, do it wordwise for speed. -countAll :: PrimeSieve -> Int -countAll (PS _ bs) = countBits bs - --- count set bits between two indices (inclusive) --- start and end must both be valid indices and start <= end -countFromTo :: Int -> Int -> MU.MVector s Bit -> ST s Int -countFromTo start end = - fmap countBits . U.unsafeFreeze . MU.slice start (end - start + 1) diff --git a/arithmoi.cabal b/arithmoi.cabal --- a/arithmoi.cabal +++ b/arithmoi.cabal @@ -1,5 +1,5 @@ name: arithmoi -version: 0.13.2.0 +version: 0.13.3.0 cabal-version: 2.0 build-type: Simple license: MIT @@ -20,6 +20,7 @@ tested-with: GHC ==9.0.2 GHC ==9.2.8 GHC ==9.4.8 GHC ==9.6.7 GHC ==9.8.4 GHC ==9.10.2 GHC ==9.12.2 GHC ==9.14.1 extra-doc-files: changelog.md + Math/NumberTheory/Primes/Counting/HowPrimeCountingWorks.md source-repository head type: git @@ -28,6 +29,7 @@ library build-depends: base >=4.15 && <5, + array >=0.5.4.0 && <0.6, bitvec >=1.1 && <1.2, containers >=0.5.11 && <0.9, chimera >=0.3 && <0.5, @@ -37,7 +39,7 @@ ghc-bignum <1.5, infinite-list <0.2, integer-logarithms >=1.0 && <1.1, - integer-roots >=1.0 && <1.1, + integer-roots >= 1.0.4.0 && <1.1, mod >=0.2.1 && <0.3, random >=1.0 && <1.4, transformers >=0.4 && <0.7, @@ -113,7 +115,7 @@ infinite-list, integer-roots >=1.0, mod, - QuickCheck >=2.10 && <2.17, + QuickCheck >=2.10 && <2.19, quickcheck-classes >=0.6.3 && <0.7, semirings >=0.2, smallcheck >=1.2 && <1.3, @@ -181,7 +183,7 @@ mod, random, semirings, - tasty-bench >= 0.4 && < 0.5, + tasty-bench >= 0.4 && < 0.6, vector other-modules: Math.NumberTheory.ArithmeticFunctionsBench diff --git a/benchmark/Math/NumberTheory/RecurrencesBench.hs b/benchmark/Math/NumberTheory/RecurrencesBench.hs --- a/benchmark/Math/NumberTheory/RecurrencesBench.hs +++ b/benchmark/Math/NumberTheory/RecurrencesBench.hs @@ -34,11 +34,11 @@ benchSuite :: Benchmark benchSuite = bgroup "Recurrences" [ bgroup "Bilinear" - [ benchTriangle "binomial" binomial 100 - , benchTriangle "stirling1" stirling1 100 - , benchTriangle "stirling2" stirling2 100 - , benchTriangle "eulerian1" eulerian1 100 - , benchTriangle "eulerian2" eulerian2 100 + [ benchTriangle "binomial" binomial 80 + , benchTriangle "stirling1" stirling1 75 + , benchTriangle "stirling2" stirling2 75 + , benchTriangle "eulerian1" eulerian1 50 + , benchTriangle "eulerian2" eulerian2 50 ] , benchPartition 1000 , bgroup "factorialFactors" diff --git a/changelog.md b/changelog.md --- a/changelog.md +++ b/changelog.md @@ -1,5 +1,13 @@ # Changelog +## 0.13.3.0 + +### Changed + +* Changed primeCount implementation to something a little faster and better documented. +* Added a documentation file: "HowPrimeCountingWorks.md". +* Tweaked some benchmark settings to make them complete faster. + ## 0.13.2.0 ### Changed diff --git a/test-suite/Math/NumberTheory/TestUtils.hs b/test-suite/Math/NumberTheory/TestUtils.hs --- a/test-suite/Math/NumberTheory/TestUtils.hs +++ b/test-suite/Math/NumberTheory/TestUtils.hs @@ -8,6 +8,7 @@ -- {-# LANGUAGE ConstraintKinds #-} +{-# LANGUAGE CPP #-} {-# LANGUAGE DataKinds #-} {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE FlexibleInstances #-} @@ -68,9 +69,11 @@ import Math.NumberTheory.TestUtils.MyCompose import Math.NumberTheory.TestUtils.Wrappers +#if !MIN_VERSION_QuickCheck(2,17,0) instance Arbitrary Natural where arbitrary = fromInteger <$> (arbitrary `suchThat` (>= 0)) shrink = map fromInteger . filter (>= 0) . shrink . toInteger +#endif instance Arbitrary E.EisensteinInteger where arbitrary = (E.:+) <$> arbitrary <*> arbitrary