Sequence of primes by trial division

Generate a sequence of primes by means of trial division.

Sequence of primes by trial division is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

Trial division is an algorithm where a candidate number is tested for being a prime by trying to divide it by other numbers. You may use primes, or any numbers of your choosing, as long as the result is indeed a sequence of primes.

The sequence may be bounded (i.e. up to some limit), unbounded, starting from the start (i.e. 2) or above some given value. Organize your function as you wish, in particular, it might resemble a filtering operation, or a sieving operation. If you want to use a ready-made is_prime function, use one from the Primality by trial division page (i.e., add yours there if it isn't there already).

• Closely related to: Primality by trial division, Sieve of Eratosthenes.

Haskell

The following code leaves no "duplicates" in its input by keeping every newly encountered distinct number, defining the equality of any two numbers as them having a GCD above 1: <lang haskell>import Data.List (nubBy)

primes = nubBy (((>1).).gcd) [2..]</lang>

It is very slow.

Sieve by trial division

Filtering candidate numbers for non-divisibility, by prime at a time, is a kind of sieving. One often sees a "naive" version presented as an unbounded sieve of Eratosthenes, similar to David Turner's 1975 SASL code, <lang haskell>primes = sieve [2..] where

  sieve (p:xs) = p : sieve [x | x <- xs, mod x p /= 0]

-- (diff xs [p, p+p, ...]) would be the sieve of Eratosthenes -- map fst . tail -- . iterate (\(_,p:xs)-> (p, [x | x <- xs, mod x p /= 0])) $ (1, [2..])</lang> As shown in Melissa O'Neill's paper "The Genuine Sieve of Eratosthenes", this is rather a sub-optimal trial division algorithm. Its complexity is quadratic in number of primes produced whereas that of optimal trial division is $O(n^{1.5}/(\log n)^{0.5})$, and of true SoE it is $O(n\log n\log \log n)$, in n primes produced.

Bounded sieve by trial division

Bounded formulation has normal trial division complexity, because it can stop early through an explicit guard: <lang haskell>primesTo m = 2 : sieve [3,5..m] where

  sieve (p:xs) | p*p > m   = p : xs
               | otherwise = p : sieve [x | x <- xs, mod x p /= 0]

-- map fst a ++ b:c where (a,(b,c):_) = span ((< m).(^2).fst) -- . iterate (\(_,p:xs)-> (p, [x | x <- xs, mod x p /= 0])) $ (2, [3,5..m])</lang>

Postponed sieve by trial division

To make it unbounded, the guard cannot be simply discarded. The firing up of a filter by a prime should be postponed until its square is seen amongst the candidates (so a bigger chunk of input numbers are taken straight away as primes, between each opening up of a new filter, instead of just one head element in the non-postponed algorithm): <lang haskell>primesT = 2 : 3 : sieve [5,7..] 9 (tail primesT) where

  sieve (x:xs) q ps@(p:t) | x < q     = x : sieve xs q ps  -- inlined (span (< q))
                          | otherwise = sieve [y | y <- xs, mod y p /= 0] (head t^2) t

-- ps = concat . map fst -- . iterate (\(_,(ns,p:t))-> let (h,xs)=span (< p*p) ns in -- (h, ([y | y <- xs, mod y p /= 0], t))) $ ([2,3], ([5,7..], tail ps))</lang>

Segmented Generate and Test

Explicating the list of filters (created implicitly at run-time by the sieves above) as a list of factors to test by on each segment between the consecutive squares of primes (so that no testing is done prematurely), and rearranging to avoid recalculations, leads to the following: <lang haskell>import Data.List (inits)

primesST = 2 : 3 : sieve 5 9 (drop 2 primesST) (inits $ tail primesST) where

  sieve x q ps (fs:ft) = filter (\y-> all ((/=0).mod y) fs) [x,x+2..q-2]
                         ++ sieve (q+2) (head ps^2) (tail ps) ft</lang>

inits makes a stream of (progressively growing) prefixes of an input stream, starting with an empty prefix, here making the fs parameter to get a sequence of values [], [3], [3,5], ....

Runs at empirical $O(n^{1.4...})$ time complexity, in n primes produced. Can be used as a framework for unbounded segmented sieves, replacing divisibility testing with proper sieve of Eratosthenes on arrays, etc.

Scala

Odds-Only "infinite" primes generator using Streams and Co-Inductive Streams

Using Streams, the "unfaithful sieve", i.e. sub-optimal trial division sieve. <lang scala>def sieve(nums: Stream[Int]): Stream[Int] =

   Stream.cons(nums.head, sieve((nums.tail).filter(_ % nums.head != 0)))
 val primes = 2 #:: sieve(Stream.from(3, 2))

 println(primes take 10 toList) //         //List(2, 3, 5, 7, 11, 13, 17, 19, 23, 29)
 println(primes takeWhile (_ < 30) toList) //List(2, 3, 5, 7, 11, 13, 17, 19, 23, 29)</lang>

Output:

Both println statements give the same results

List(2, 3, 5, 7, 11, 13, 17, 19, 23, 29)

The above code is extremely inefficient for larger ranges, both because it tests for primality using computationally expensive divide (modulo) operations and because it sets up deferred tests for division by all of the primes up to each prime candidate, meaning that it has approximately a square law computational complexity with range.

Racket

Infinite list of primes:

Using laziness

This example uses infinite lists (streams) to implement a sieve algorithm that produces all prime numbers.

<lang Racket>#lang lazy (define nats (cons 1 (map add1 nats))) (define (sift n l) (filter (λ(x) (not (zero? (modulo x n)))) l)) (define (sieve l) (cons (first l) (sieve (sift (first l) (rest l))))) (define primes (sieve (rest nats))) (!! (take 25 primes))</lang>

Using threads and channels

Same algorithm as above, but now using threads and channels to produce a channel of all prime numbers (similar to newsqueak). The macro at the top is a convenient wrapper around definitions of channels using a thread that feeds them.

<lang Racket>#lang racket (define-syntax (define-thread-loop stx)

 (syntax-case stx ()
   [(_ (name . args) expr ...)
    (with-syntax ([out! (datum->syntax stx 'out!)])
      #'(define (name . args)
          (define out (make-channel))
          (define (out! x) (channel-put out x))
          (thread (λ() (let loop () expr ... (loop))))
          out))]))

(define-thread-loop (nats) (for ([i (in-naturals 1)]) (out! i))) (define-thread-loop (filter pred? c)

 (let ([x (channel-get c)]) (when (pred? x) (out! x))))

(define (sift n c) (filter (λ(x) (not (zero? (modulo x n)))) c)) (define-thread-loop (sieve c)

 (let ([x (channel-get c)]) (out! x) (set! c (sift x c))))

(define primes (let ([ns (nats)]) (channel-get ns) (sieve ns))) (for/list ([i 25] [x (in-producer (λ() (channel-get primes)))]) x)</lang>

Using generators

Yet another variation of the same algorithm as above, this time using generators.

<lang Racket>#lang racket (require racket/generator) (define nats (generator () (for ([i (in-naturals 1)]) (yield i)))) (define (filter pred g)

 (generator () (for ([i (in-producer g #f)] #:when (pred i)) (yield i))))

(define (sift n g) (filter (λ(x) (not (zero? (modulo x n)))) g)) (define (sieve g)

 (generator ()
            (let loop ([g g]) (let ([x (g)]) (yield x) (loop (sift x g))))))

(define primes (begin (nats) (sieve nats))) (for/list ([i 25] [x (in-producer primes)]) x)</lang>