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.

Ada

Use the generic function Prime_Numbers.Is_Prime, as specified in Prime decomposition#Ada. The program reads two numbers A and B from the command line and prints all primes between A and B (inclusive).

<lang Ada>with Prime_Numbers, Ada.Text_IO, Ada.Command_Line;

procedure Sequence_Of_Primes is

  package Integer_Numbers is new 
    Prime_Numbers (Natural, 0, 1, 2); 
  use Integer_Numbers;
  
  Start: Natural := Natural'Value(Ada.Command_Line.Argument(1));
  Stop:  Natural := Natural'Value(Ada.Command_Line.Argument(2));
  

begin

  for I in Start .. Stop loop
     if Is_Prime(I) then 
        Ada.Text_IO.Put(Natural'Image(I));
     end if;
  end loop;

end Sequence_Of_Primes;</lang>

Output:

>./sequence_of_primes 50 99
 53 59 61 67 71 73 79 83 89 97

D

Translation of: Haskell

This is a quite inefficient prime generator. <lang d>import std.stdio, std.range, std.algorithm, std.traits,

      std.numeric, std.concurrency;

Generator!(ForeachType!R) nubBy(alias pred, R)(R items) {

   return new typeof(return)({
       ForeachType!R[] seen;

       OUTER: foreach (x; items) {
           foreach (y; seen)
               if (pred(x, y))
                   continue OUTER;
           yield(x);
           seen ~= x;
       }
   });

}

void main() /*@safe*/ {

   sequence!q{n + 2}
   .nubBy!((x, y) => gcd(x, y) > 1)
   .take(20)
   .writeln;

}</lang>

Output:

[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71]

Go

An unbounded cascading filtering method using channels, adapted from the classic concurrent prime sieve example in the "Try Go" window at http://golang.org/, improved by postponing the initiation of the filtering by a prime until the prime's square is seen in the input. <lang go>package main

import "fmt"

func NumsFromBy(from int, by int, ch chan<- int) {

 for i := from; ; i+=by {
   ch <- i
 }

}

func Filter(in <-chan int, out chan<- int, prime int) {

 for {
   i := <-in
   if i%prime != 0 {            // here is the trial division
     out <- i
   }
 }

}

func Sieve(out chan<- int) {

 out <- 2                   
 out <- 3
 q := 9
 ps := make(chan int)
 go Sieve(ps)                   // separate primes supply
 p := <-ps        
 p = <-ps         
 nums := make(chan int)
 go NumsFromBy(5,2,nums)        // end of setup
 for i := 0; ; i++ {
   n := <-nums
   if n < q {
   	out <- n                 // n is prime
   } else {
       ch1 := make(chan int)
       go Filter(nums, ch1, p)  // postponed creation of a filter
       nums = ch1
   	p = <-ps                
   	q = p*p                  
   }
 }

}

func main() {

 ch := make(chan int)
 go Sieve(ch)
 fmt.Print("First twenty:")
 for i := 0; i < 20; i++ {
   fmt.Print(" ", <-ch)
 }
 fmt.Println()

}</lang>

Output:

First twenty: 2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71

A simple iterative method, also unbounded and starting with 2. <lang go>package main

import "fmt"

func newP() func() int {

   n := 1
   return func() int {
       for {
           n++
           // Trial division as naïvely as possible.  For a candidate n,
           // numbers between 1 and n are checked to see if they divide n.
           // If no number divides n, n is prime.
           for f := 2; ; f++ {
               if f == n {
                   return n
               }
               if n%f == 0 { // here is the trial division
                   break
               }
           }
       }
   }

}

func main() {

   p := newP()
   fmt.Print("First twenty:")
   for i := 0; i < 20; i++ {
       fmt.Print(" ", p())
   }
   fmt.Println()

}</lang>

Output:

First twenty: 2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71

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.

The classical version has isPrime inlined, and so can ignore the prime 2 while testing only odd numbers:

<lang haskell>primes = 2 : 3 : [n | n <- [5,7..], foldr (\p r-> p*p > n || rem n p > 0 && r)

                                         True (drop 1 primes)]</lang>

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>primesT = sieve [2..]

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

-- map head -- . iterate (\(p:xs)-> [x | x <- xs, rem x p /= 0]) $ [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 via 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, rem x p /= 0]

-- map fst a ++ b:c where (a,(b,c):_) = span ((< m).(^2).fst) -- . iterate (\(_,p:xs)-> (p, [x | x <- xs, rem 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>primesPT = 2 : 3 : sieve [5,7..] 9 (tail primesPT)

  where
  sieve (x:xs) q ps@(p:t) 
       | x < q     = x : sieve xs q ps    -- inlined (span (< q))
       | otherwise = sieve [y | y <- xs, rem 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, rem y p /= 0], t))) $ ([2,3], ([5,7..], tail ps))</lang>

Segmented Generate and Test

Explicating the run-time list of filters (created implicitly 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).rem 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.

J

Generally, this task should be accomplished in J using (#~ 1&p:). Here we take an approach that's more comparable with the other examples on this page.

Implementation: <lang J>primTrial=:3 :0

 try=. i.&.(p:inv) %: >./ y
 candidate=. (y>1)*y=<.y
 y #~ candidate*(y e.try) = +/ 0= try|/ y

)</lang>

Example use:

<lang J> primTrial 1e6+i.100 1000003 1000033 1000037 1000039 1000081 1000099</lang>

Note that this is a filter - it selects values from its argument which are prime. If no suitable values are found the resulting sequence of primes will be empty.

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

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>

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>

REXX

somewhat optimized

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.
Usage note:   by using a negative number (for the program's argument), the list of primes is suppressed, but the prime count is still shown. <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*/ tell=n>0; n=abs(n) /*N is negative? Don't display.*/ @.1=2; if tell then 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.*/
  if tell  then say right(j, 9)       /*maybe display this prime──►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.

more optimized

This version shows how the REXX program may be optimized further by extending the list of low primes and the special low prime divisions   (//   tests). <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*/ tell=n>0; n=abs(n) /*N is negative? Don't display.*/ @.1=2; @.2=3; @.3=5; @.4=7; @.5=11; @.6=13; #=5; s=@.#+2

                                      /*   [↑]  is the # of low primes.*/
     do p=1  for #  while  p<=n       /* [↓]  don't compute, just list.*/
     if tell  then say right(@.p,9)   /*display some low primes.       */
     !.p=@.p**2                       /*also compute the squared value.*/
     end   /*p*/                      /* [↑]  allows faster loop below.*/
                                      /* [↓]  N default lists up to 101*/
  do j=s  by 2  while  #<n            /*start with the next odd prime. */
  if j//  3    ==0   then iterate     /*is this number a triple?       */
  if right(j,1)==5   then iterate     /* "   "     "   " nickel?       */
  if j//  7    ==0   then iterate     /* "   "     "     lucky?        */
  if j// 11    ==0   then iterate     /* "   "     "  a multiple of 11?*/
                                      /* [↓]  divide by the primes. ___*/
         do k=p  to #  while  !.k<=j  /*divide J by other 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.*/
  if tell  then say right(j, 9)       /*maybe display this prime──►term*/
  end     /*j*/                       /* [↑]  only display  N  primes. */
                                      /* [↓]  display the prime count. */

say # ' primes found.' /*stick a fork in it, we're done.*/</lang> output is the same as the 1^st REXX version.

Ruby

The Prime class in the standard library has several Prime generators. In some methods it can be specified which generator will be used. The generator can be used on it's own: <lang ruby>require "prime"

pg = Prime::TrialDivisionGenerator.new p pg.take(10) # => [2, 3, 5, 7, 11, 13, 17, 19, 23, 29] p pg.next # => 31</lang>

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.

Tcl

As we're generating a sequence of primes, we can use that sequence of primes to describe what we're filtering against. <lang tcl>set primes {} proc havePrime n {

   global primes
   foreach p $primes {

# Do the test-by-trial-division if {$n/$p*$p == $n} {return false}

   }
   return true

} for {set n 2} {$n < 100} {incr n} {

   if {[havePrime $n]} {

lappend primes $n puts -nonewline "$n "

   }

} puts ""</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