Sequence of primes by trial division

Sequence of primes by trial division
You are encouraged to solve this task according to the task description, using any language you may know.

Generate a sequence of primes by means of trial division.

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

Trial division can be used to produce short ranges of primes: <lang haskell>primesFromTo n m = filter isPrime [n..m]</lang>

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.

The filtering function is equivalent to noDivsBy defined as part of isPrime function, except that the comparisons testing for the square root are no longer needed and so are spared.

Perl 6

Here is a straightforward implementation of the naive algorithm. <lang perl6>constant @primes = 2, 3, { first * %% none(@_), (@_[* - 1], * + 2 ... *) } ... *;

say @primes[^100];</lang>

Output:

2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97 101 103 107 109 113 127 131 137 139 149 151 157 163 167 173 179 181 191 193 197 199 211 223 227 229 233 239 241 251 257 263 269 271 277 281 283 293 307 311 313 317 331 337 347 349 353 359 367 373 379 383 389 397 401 409 419 421 431 433 439 443 449 457 461 463 467 479 487 491 499 503 509 521 523 541

REXX

This is an open-ended approach and it's a simple implementation and could be optimized more with some easy programming.
The method used is to divided all odd numbers by all previous odd primes up to and including the   √   of the odd number. <lang rexx>/*REXX pgm lists a sequence of primes by testing primality by trial div.*/ parse arg n . /*let user choose how many, maybe*/ if n== then n=26 /*if not, then assume the default*/ if n<1 then exit /*handle special case of 0 primes*/ @.1=2; say right(@.1, 9) /*display 2 as a special case. */

1. =1 /*# is the number of primes found*/

                                      /* [↑]  N default lists up to 101*/
  do j=      3  by 2  while  #<n      /*start with the first odd prime.*/
                                      /* [↓]  divide by the primes. ___*/
         do k=2  to #  while  !.k<=j  /*divide J with all primes ≤ √ J */
         if j//@.k==0  then iterate j /*÷ by a previous prime?  ¬prime.*/
         end   /*j*/                  /* [↑]   only divide up to  √J.  */
  #=#+1                               /*bump the prime number counter. */
  @.#=j;   !.#=j*j                    /*define this prime & its square.*/
  say right(j, 9)                     /*display this prime to the term.*/
  end     /*j*/                       /* [↑]  only display  N  primes. */
                                      /* [↓]  display the prime count. */

say # ' primes found.' /*stick a fork in it, we're done.*/</lang> output using the default input of:   26

        2
        3
        5
        7
       11
       13
       17
       19
       23
       29
       31
       37
       41
       43
       47
       53
       59
       61
       67
       71
       73
       79
       83
       89
       97
      101
26  primes found.

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.

<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>

Optimized with postponed processing

Since a prime's multiples that count start from its square, we should only add them when we reach that square.

<lang Racket>#lang lazy (define nats (cons 1 (map add1 nats))) (define (sift n l) (filter (λ(x) (not (zero? (modulo x n)))) l)) (define (when-bigger n l f)

 (if (< (car l) n) (cons (car l) (when-bigger n (cdr l) f)) (f l)))

(define (sieve l ps)

 (cons (car l) (when-bigger (* (car ps) (car ps)) (cdr l)
                            (λ(t) (sieve (sift (car ps) t) (cdr ps))))))

(define primes (sieve (cdr nats) primes)) (!! (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>