Determining primality by removing select composits

Calculating many remainders for each prime and composite consumes more time.

I use a factoring function to find factors of single values.

When I use a factoring function to find lists of primes, I use a wheel.

A wheel exploits composite values punctuating an otherwise prime sequence.

The primary advantage of a wheel is reducing significantly, the candidate list and so processing.

My wheel eliminates 2, 3, 5 and 7 multiples from the a list. For the longest time, I used only elimination of 2s, 3s and 5s because the wheel list is only 8 values and simple to construct but there are just to many 7s.

Take the list of primes 11,13,17,19,23,29,31,39 and add 30 to each of the previous 8 values for subsequent values.

Take the first 55 from the list and remove the 7 multiples. The deltas will be this wheel.

What should be faster is subtraction of one list from another where the calc is compare and once only for each.

wl = [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]
-- the so-called wheel: the recurring pattern in a composite punctuated list
n7sl = scanl (+) 11 $cycle wl -- infinite -- short list of composits with no 2's, 3's 5's or 7's; 77%+ reduced n7s = take 150$ n7sl

-- diagonalize and limit the factor's composite list
lls i y = drop i.take y $n7sl -- substract the 1st list from the 2nd, infinite list rms [] _ = [] -- stop when first list is exhausted rms n@(x:xs) p@(y:ys) | x==y = rms xs ys | x>y = y:rms n ys | y>x = rms xs p -- generate the composite list comps [] _ _ _ _ =[] comps (n:ns) [] y i b = comps ns (lls (i+1) y) y (i+1) b comps k@(n:ns) (r:rs) y i b | m>b =comps ns (lls (i+1) y) y (i+1) b | m<=b =m:comps k rs y i b where m = n*r -- the result of comps is not just a diagonalization but 2, top & irregular bottom -- comps is maybe necessary to stop when the limit of a list is reached -- otherwise, I was using a list comprehension which is way less complicated -- comps has to many parameters so this to reduce it to 1 -- it will be subtracted the fixed length wheel numbers and the lazy wheel comp1 n =sort$ comps n7s (lls 0 n) n 0 (11*(last.take n $n7sl)) -- put everything together to generate primes  last$ rms (comp1 5000) n7sl

240641

(0.12 secs, 65,219,592 bytes)

last $rms (comp1 10000) n7sl 481249 (0.24 secs, 143,743,944 bytes) last$ rms (comp1 15000) n7sl

721891

(0.37 secs, 220,099,904 bytes)

Question is this better that factoring? Is there a better way to generate multiples without multiplying in comps or list comprehensions?

Any prime list I generate excludes 2,3,5 & 7 because they are misbehaved, irregular and to well known.

• Well, I had residuals. rms was written to remove 7s from 2s,3s,5s list. Remove sevens: rms Then, in comps I had had trouble with the multiple. It was in a where statement b/c it occurs 3 times. Then, after fixing another problem forgot to put it back. Mar 27, 2019 at 20:33

You should use more Data.List functions - explicit recursion is hard to read.

comp1 y = sort [ m
| (i,n) <- zip [0..149] n7sl
, m <- takeWhile ((<= 11 * last (take y n7sl))) $map (n*)$ drop i $take n n7sl ]  Edit: Not sure whether my version fails to be equal to your original one, but here's a further refactoring of the one in your comment to get rid of some number manipulation. The takeWhile reflects the short-circuiting behavior of the original, explicitly recursive implementation, which uses that bs is ascending. t = sort [ m | bs@(a:_) <- tails n7s, m <- takeWhile (<= 11*last n7s)$ map (a*) bs]

• Well, thank you. Your comprehension removes a fixed amount from each factor sublist. Mine removes a varied amount according to value. Yours does not provide near enough of any. There should n multiples of 11 for comp1 n. Finally, it must work as a select list to subtract from the wheel to produce primes. Yours does not. Mar 28, 2019 at 15:02
• A list comprehension to produce the correct composits n7s = take 50 n7sl; t = sort [ m | (a,i) <- zip n7s [0..], b <- drop i n7s, m <- [a*b], m <= (11*(last n7s))]; rms t n7sl Mar 28, 2019 at 15:46
• Well, I'm so impressed. I am also sorry. When I checked your results with mine I found a mismatch but it was because I had forgot to sort mine so I assumed yours was wrong. My bad. It is right. It is impressive and so clean. Mar 29, 2019 at 16:38
• With all the other advantages of this, it is also the fastest, if you move off the limit calculation with something like >nt lim = sort [ m | bs@(a:_) <- tails n7s, m <- takeWhile (<= lim) $map (a*) bs] with ntr = nt (11*last n7s). It should speed up the primes program nicely, Thanks again! Mar 31, 2019 at 1:51 After simplifications, your infinite code (from your answer) turns out to be equivalent to import Data.List (tails) import Data.List.Ordered (minus, unionAll) primes = ([2,3,5,7] ++) . minus n7sl . unionAll$
[ map (x*) xs | xs@(x:_) <- tails n7sl ]


The main difference is that whereas your code builds the matrix

> mapM_ print $take 20$ zip n7sl [take i n7sl | i <- [1..]]
(11,[11])
(13,[11,13])
(17,[11,13,17])
(19,[11,13,17,19])
(23,[11,13,17,19,23])
(29,[11,13,17,19,23,29])
(31,[11,13,17,19,23,29,31])
(37,[11,13,17,19,23,29,31,37])
(41,[11,13,17,19,23,29,31,37,41])
(43,[11,13,17,19,23,29,31,37,41,43])
(47,[11,13,17,19,23,29,31,37,41,43,47])
(53,[11,13,17,19,23,29,31,37,41,43,47,53])
(59,[11,13,17,19,23,29,31,37,41,43,47,53,59])
(61,[11,13,17,19,23,29,31,37,41,43,47,53,59,61])
(67,[11,13,17,19,23,29,31,37,41,43,47,53,59,61,67])
(71,[11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71])
(73,[11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73])
(79,[11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79])
(83,[11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83])
(89,[11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89])


by rows (in each row, the list is multiplied by the first number in the tuple), and processes it by rows -- putting them through the unionAll process -- the rewrite in this answer is as if working on the same matrix by columns (after the multiplication by the first number in the tuple, again).

Because tails is much less computationally demanding compared with your repeated use of take, and because we use minus from Data.List.Ordered package instead of your rms (with the flipped order of arguments), this runs much much faster now. Testing in GHCi:

> primes !! 1000000
15485867
it :: Integral a => a
(2.12 secs, 3734834160 bytes)

> primes !! 500000
7368791
it :: Integral a => a
(0.98 secs, 1765893576 bytes)

> logBase 2 (2.12 / 0.98)
1.113210610447991


Yes, that's a million primes that it now reaches, in just over 2 seconds, at about n1.1 empirical orders of growth (in n primes produced), which is quite good (the coefficient, not the time; arithmoi's code reaches the one millionth prime in 0.1 seconds on the same computer).

This is not a sieve of Eratosthenes though. It builds the multiples not from primes, but from the 2-3-5-7-wheel enumeration.

For proper testing always compile with the -O2 switch and run the resulting standalone executable at the shell, with +RTS -s run-time options (to get the time and space statistics).

• You so often amaze me and you are always right. 'n7sl' is remove 1st 4 prime multiples wheel. I had since changed the infinite function composites to (w/no take) umap = unionAll $map (\n -> (n*) <$> n7sl) n7sl ; take 50000 $rms umap n7sl But the non-infinite above still outperforms it and yours (IDK) and makes me mad Both mine above eschew trial division. May 28, 2019 at 0:39 • which non-infinite above? please give a test expression to run. May 28, 2019 at 6:08 • The infinite is now way faster. My misconception was remarkable, in a bad way. I thought I had to do inits for an infinite application. It put my factor at the end of each list so I had to use two lists and was penalized it time. tails is congruent with Haskell pattern matching and the factor is now at the head of each list by using tails. @Will Ness is the best. inits or tails diagonalizes the Cartesian product, upper or lower. May 29, 2019 at 17:30 • did you notice that using n7sl there are no non-prime numbers until 121? Your diagonal of my function is all primes. In Haskell Data.Numbers.Primes the wheelSieve take a parameter of the number of initial primes used to create the wheel and the larger wheel speeds it up. In my non-infinite list the larger wheel slowed it down. Maybe it won't in the infinite function. May 29, 2019 at 17:37 Breaking up comps results in simpler, more efficient, easier to control functions. Two lists (x,y) for a Cartesian product are usual. The products are half duplicates. If each column starts with a perfect square, the duplicates are eliminated. Then each column is generated. Each column starts with a perfect square but then must also be specifically trimmed in length. Each multiple is compared to the limititing value. So, about 25% or less of the Cartesian product are used, that is, calculated. calcc limits the value of each factor multiple it calculates just like comps but unlike comps it doesn't use an initial-factor list. Each initial-factor squared is first on the multiple list. Each initial-factor is the head of each factor lists input. I tried making comps use one list only instead of two. The function was short, more complicated but not faster. Two remaining problems with it were the calculation for the limit value which only needs to be calculated one time at the beginning and a problem with all such functions a glob of parameters or superflous where calculations. The functions need the parameters but should not be recalculated time-after-time. runr passes the initial-factor to apply to itself and all other factors. The list sent to calcc starts with a perfect square factor value. The diagonal eliminates duplicates. rrunr sarts runr with a length n factor list and the limit based on that list length. The 11 multiples extend to n. Successive columns shorten to none at about n/10 Checking every perfect square against the limit does not save time. I also tried mergeAll from Data.List.Ordered but it seemed to error around 3/4 the way down a prime list. I still have to investigate. sort is what I still use but the functions now let me create a single list or a list of lists for mergeAll. calcc _ _ []= [] calcc f lim (x:xs) | m<= lim= m:calcc f lim xs | True= [] where m= f*x runr _ []= [] runr lim xss@(x:xs) = (calcc x lim xss)++runr lim xs rrunr n = runr lim ds where ds= take n n7sl lim= 11*last ds  last.rms (sort$ rrunr 5000) $n7sl 240641 (0.07 secs, 50,617,744 bytes) last.rms (sort$ rrunr 10000) $n7sl 481249 (0.18 secs, 115,282,560 bytes) last.rms (sort$ rrunr 15000) $n7sl 721891 (0.25 secs, 179,096,928 bytes) 5/13/2019 I tried unionAll instead of mergeAll because each sub-list is already in order. It did result in a minor speed up. runr produces individual sub-lists for unionAll. runr _ []= [] runr lim xss@(x:xs) = [calcc x lim xss]++runr lim xs  last$ rms (unionAll $rrunr 10000) n7sl 481249 (0.12 secs, 89,250,264 bytes) last$ rms (unionAll $rrunr 15000) n7sl 721891 (0.20 secs, 140,110,392 bytes) last$ rms (unionAll $rrunr 20000) n7sl 962509 (0.27 secs, 192,858,656 bytes) I have an infinite version of this which is much simpler coding but it is a little slower. mupcl _ [] = [] mupcl (x:xs) (y:ys) = [(x*) <$> y]++ mupcl xs ys

uat = unionAll.mupcl n7sl $(\i-> (take i n7sl)) <$> [1..]


take 50000 $rms uat n7sl 611993 (0.45 secs, 269,026,264 bytes) I updated this because finding the nth prime is important and this can. It's still slower than the fixed just before this but has some uses not available to the fixed. I toyed with making the sublists sequences but then found that I could generate the sublist end value and multiply it by all previous values and itself. Sequence would access the end value in about O(1) which is awesome but the end values are in one-to-one correspondence with each sublist. The diagonal of these is just each sublist is one longer than the previous.This is the top diagonal necessary for an infinite list. The previous functions use a bottom diagonal and a limit value. If they used the top diagonal would they end as this function? IDK 5/27/2019 The composite list was being traversed twice and it was driving me nuts. I could not come up with a single function that would multiply two identical infinite lists diagonally. I had developed one function to multiply 2s, 3s & 5s together to form the Hamming list. Saturday morning I looked at that function and followed the pattern for this function. It amazed me it worked. unionAll does the limiting. umap = unionAll$ map (\n -> (n*) <$> n7sl) n7sl  But, I am such an impatient idiot the function is going through two lists, exactly what I didn't want. I have an aversion to tails, more like a phobia because of experience with tails crashing my PC in an infinite function. What I discovered the past two or three months is that mergeAll and unionAll do magic with infinite lists. They automatically limit end points. Talk about no coding. I spent so much time coding end points in this and in the Hamming numbers (which are way fast, too) that I really feel like an idiot. @Will Ness knows way more than I, too. Haskell is pure (no pun) magic with laziness. c2 = unionAll$ map (\xs@(x:_)-> (x*) <$> xs)$ tails n7sl


(minus n7sl c2) !! 1000000

15485941

(1.55 secs, 3,545,902,888 bytes)

And now, it's way faster than the non-infinite above.

This does 75803 primes to 962509 in 0.08 and the non-infinite does 75803 to 962509 in 0.27 and uses a little more memory.

This does

(minus n7sl c2) !! 75803

962509

(0.08 secs, 162,030,280 bytes)

5.30.2019

OMG. It worked! I told @Will Ness it might, but that it would not with the non-infinite. What it is is the size of the wheel. I started calling it a wheel when I saw the deltas to create it. Mine is not a wheel because I have to generate the deltas from a column in Excel. I take 300 of n7sl into a column in Excel them I remove all the 11 multiples. I take the deltas of the remainder then find where it repeats.

My no 11s list is

n11s=[4,2,4,6,2,6,4,2,4,6,6,2,6,4,2,6,4,6,8,4,2,4,2,4,14,4,6,2,10,2,6,6,4,2,4,6,2,10,2,4,2,12,10,2,4,2,4,6,2,6,4,6,6,6,2,6,4,2,6,4,6,8,4,2,4,6,8,6,10,2,4,6,2,6,6,4,2,4,6,2,6,4,2,6,10,2,10,2,4,2,4,6,8,4,2,4,12,2,6,4,2,6,4,6,12,2,4,2,4,8,6,4,6,2,4,6,2,6,10,2,4,6,2,6,4,2,4,2,10,2,10,2,4,6,6,2,6,6,4,6,6,2,6,4,2,6,4,6,8,4,2,6,4,8,6,4,6,2,4,6,8,6,4,2,10,2,6,4,2,4,2,10,2,10,2,4,2,4,8,6,4,2,4,6,6,2,6,4,8,4,6,8,4,2,4,2,4,8,6,4,6,6,6,2,6,6,4,2,4,6,2,6,4,2,4,2,10,2,10,2,6,4,6,2,6,4,2,4,6,6,8,4,2,6,10,8,4,2,4,2,4,8,10,6,2,4,8,6,6,4,2,4,6,2,6,4,6,2,10,2,10,2,4,2,4,6,2,6,4,2,4,6,6,2,6,6,6,4,6,8,4,2,4,2,4,8,6,4,8,4,6,2,6,6,4,2,4,6,8,4,2,4,2,10,2,10,2,4,2,4,6,2,10,2,4,6,8,6,4,2,6,4,6,8,4,6,2,4,8,6,4,6,2,4,6,2,6,6,4,6,6,2,6,6,4,2,10,2,10,2,4,2,4,6,2,6,4,2,10,6,2,6,4,2,6,4,6,8,4,2,4,2,12,6,4,6,2,4,6,2,12,4,2,4,8,6,4,2,4,2,10,2,10,6,2,4,6,2,6,4,2,4,6,6,2,6,4,2,10,6,8,6,4,2,4,8,6,4,6,2,4,6,2,6,6,6,4,6,2,6,4,2,4,2,10,12,2,4,2,10,2,6,4,2,4,6,6,2,10,2,6,4,14,4,2,4,2,4,8,6,4,6,2,4,6,2,6,6,4,2,4,6,2,6,4,2,4,12,2,12]
n11sl = scanl (+) 13 $cycle n11s  The improvement is not stark and going any further with this would not be beneficial. c2 modified to use n11sl c2 = unionAll$ map (\xs@(x:_) -> (x*) <$> xs)$ tails n11sl


Executing it would be

(11:(minus n11sl c2)) !! 1000000

15485941

(1.30 secs, 2,968,456,072 bytes)

(11:(minus n11sl c2)) !! 500000

7368841

(0.55 secs, 1,305,594,104 bytes)