Sequence of primes by trial division

From Rosetta Code
Task
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).

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

ATS

<lang ATS>(* // Lazy-evaluation: // sieve for primes

  • )

(* ****** ****** *) // // How to compile: // with no GC: // patscc -D_GNU_SOURCE -DATS_MEMALLOC_LIBC -o sieve sieve.dats // with Boehm-GC: // patscc -D_GNU_SOURCE -DATS_MEMALLOC_GCBDW -o sieve sieve.dats -lgc // (* ****** ****** *) //

  1. include

"share/atspre_staload.hats" // (* ****** ****** *)

  1. define :: stream_cons
  2. define cons stream_cons
  3. define nil stream_nil

(* ****** ****** *) // fun from{n:int} (n: int n) 

 :<!laz> stream (intGte(n)) = $delay (cons{intGte(n)}(n, from (n+1)))

// (* ****** ****** *)

typedef N2 = intGte(2)

(* ****** ****** *)

fun sieve ( 

 ns: stream N2

) :<!laz> 

 stream (N2) = $delay

( let 

 val-cons(n, ns) = !ns

in 

 cons{N2}(n, sieve (stream_filter_cloref<N2> (ns, lam x => g1int_nmod(x, n) > 0)))

end : stream_con (N2) ) // end of [sieve]

//

val primes: stream (N2) = sieve (from(2))

//

macdef prime_get (n) = stream_nth_exn (primes, ,(n))

//

implement main0 () = begin // println! ("prime 1000 = ", prime_get (1000)) ; // = 7927 (* println! ("prime 5000 = ", prime_get (5000)) ; // = 48619 println! ("prime 10000 = ", prime_get (10000)) ; // = 104743

  • )

// end // end of [main0]

(* ****** ****** *)</lang>

Befunge

Based on the test in the Primality by trial division task, this list all primes between 2 and 1,000,000.

<lang befunge>2>:::"}"8*:*>`#@_48*:**2v v_v#`\*:%*:*84\/*:*84::+< v#>::48*:*/\48*:*%%!#v_1^ <^+1$_.#<5#<5#<+#<,#<<0:\</lang>

Output:
2
3
5
7
11
13
.
.
.
999931
999953
999959
999961
999979
999983

C++

<lang cpp>

  1. include <math.h>
  2. include <iostream>
  3. include <iomanip>

bool isPrime( unsigned u ) { 

   if( u < 4 ) return u > 1;
   if( /*!( u % 2 ) ||*/ !( u % 3 ) ) return false;

   unsigned q = static_cast<unsigned>( sqrt( static_cast<long double>( u ) ) ),
            c = 5;
   while( c <= q ) {
       if( !( u % c ) || !( u % ( c + 2 ) ) ) return false;
       c += 6;
   }
   return true;

} int main( int argc, char* argv[] ) { 

   unsigned mx = 100000000,
            wid = static_cast<unsigned>( log10( static_cast<long double>( mx ) ) ) + 1;

   std::cout << "[" << std::setw( wid ) << 2 << " ";
   unsigned u = 3, p = 1; // <- start computing from 3
   while( u < mx ) {
       if( isPrime( u ) ) { std::cout << std::setw( wid ) << u << " "; p++; }
       u += 2;
   }
   std::cout << "]\n\n Found " << p << " primes.\n\n";
   return 0;

} </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 ...

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]

Eiffel

<lang Eiffel> class APPLICATION

create make

feature

make do sequence (1, 27) end

sequence (lower, upper: INTEGER) -- Sequence of primes from 'lower' to 'upper'. require lower_positive: lower > 0 upper_positive: upper > 0 lower_smaller: lower < upper local i: INTEGER do io.put_string ("Sequence of primes from " + lower.out + " up to " + upper.out + ".%N") i := lower if i \\ 2 = 0 then i := i + 1 end from until i > upper loop if is_prime (i) then io.put_integer (i) io.put_new_line end i := i + 2 end end

feature {NONE}

is_prime (n: INTEGER): BOOLEAN -- Is 'n' a prime number? require positiv_input: n > 0 local i: INTEGER max: REAL_64 math: DOUBLE_MATH do create math if n = 2 then Result := True elseif n <= 1 or n \\ 2 = 0 then Result := False else Result := True max := math.sqrt (n) from i := 3 until i > max loop if n \\ i = 0 then Result := False end i := i + 2 end end end

end </lang>

Output:
Sequence of primes from 1 to 27.)
3
5
7
11
13
17
19
23

ERRE

<lang ERRE> PROGRAM PRIME_GENERATOR

!$DOUBLE

BEGIN 

  PRINT(CHR$(12);) !CLS
  N=1
  LOOP
    N+=1
    FOR F=2 TO N DO
      IF F=N THEN PRINT(N;) EXIT END IF
      EXIT IF N=F*INT(N/F)
    END FOR
  END LOOP

END PROGRAM </lang> You must press Ctrl+Break to stop the program.

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  547
 557  563  569  571  577  587  593  599  601  607  613  617  619  631  641  643
 647  653  659  661  673  677  683  691  701  709  719  727  733  739  743  751
 757  761  769  773  787  797  809  811  821  823  827  829  839  853  857  859
 863  877  881  883  887  907  911  919  929  937  941  947  953  967  971  977
 983  991  997  1009  1013  1019  1021  1031  1033  1039  1049  1051  1061
 1063  1069  1087  1091  1093  1097  1103  1109  1117  1123  1129  1151  1153
 1163  1171  1181  1187  1193  1201  1213  1217  1223  1229  1231  1237  1249
 1259  1277  1279  1283  1289  1291  1297  1301  1303  1307  1319  1321  1327
 1361  1367  1373  1381  1399  1409  1423  1427  1429  1433  1439  1447  1451
 1453  1459  1471  1481  1483  1487  1489  1493  1499  1511  1523  1531  1543
 1549  1553  1559  1567  1571  1579  1583  1597  1601 ^C

Fortran

This version was written for an IBM1130, which offered 16-bit integers. Storage size was 32K words. It was written before the revised dogma that one is not a prime number became common. With punched-card input, using only capital letters was normal. This system offered Fortran IV which lacks the many later developments such as if ... then style statements, thus the profusion of statement labels and arithmetic if-statements. In if (expression) negative,zero,positive the sign of the expression is examined to choose the appropriate label to jump to. Labels can only be numbers, alas. There was also no MOD function for determining the remainder.

This routine does not attempt a calculation of sqrt(n), a time-consuming and also potentially inacurrate scheme. For instance, using Trunc(Log10(n)) + 1 to determine how many digits are needed to print n seems an obvious ad-hoc ploy, and gave 1 for n = 1 to 9, but it also gave 1 for n = 10, because Log10(10) was calculated as 0.9999964 or so (single precision on an IBM 390 system), which truncates to zero.

Instead, since many successive numbers are to be tested for primality, advantage can be taken of the fact that prime numbers only up to PRIME(LP) need be tried as factors all the way up to N = PRIME(LP + 1)**2 = XP2. This is similar to starting a pass through a Sieve of Eratoshenes at P*P rather than at 2*P. Thus, 5 is the largest factor to try, even beyond 5*5, all the way up to 49, because, if the number were divisible by 7, 7*2 would already have been checked because of 2, 7*3 because of 3, and so on. Only when 7*7 is needed will a further possible factor have to be tried. Likewise, although the last possible factor to try for N up to the integer limit of 32767 is 181 because the square of the next prime (191) exceeds 32767, in order for the method to be able to know this, the PRIME array must have space for this surplus prime. However, it does not know this properly because the square of 191 does exceed 32767 and so its value in XP2 will be incorrect, but this doesn't matter because only equality to XP2 is checked for and there will never be call to try 191 as a factor because 181 suffices up to the integer limit and the iteration will stop by then. Fortunately, 32767 is not divisible by three so that value will not be excluded as a possible candidate for N, and so the search can correctly end after inspecting the largest possible integer - finding it divisible by seven.

This method avoids considering multiples of two and three, leading to the need to pre-load array PRIME and print the first few values explicitly rather than flounder about with special startup tricks. Even so, in order not to pre-load with 7, and to correctly start the factor testing with 5, the first few primes are found with some wasted effort because 5 is not needed at the start. Storing the primes as found has the obvious advantage of enabling divisions only by prime numbers, but care with the startup is needed to ensure that primes have indeed been stored before they are called for.

<lang Fortran> CONCOCTED BY R.N.MCLEAN, APPLIED MATHS COURSE, AUCKLAND UNIVERSITY, MCMLXXI. 

     INTEGER ENUFF,PRIME(44)

CALCULATION SHOWS PRIME(43) = 181, AND PRIME(44) = 191. 

     INTEGER N,F,Q,XP2
     INTEGER INC,IP,LP,PP
     INTEGER ALINE(20),LL,I
     DATA ENUFF/44/
     DATA PP/4/
     DATA PRIME(1),PRIME(2),PRIME(3),PRIME(4)/1,2,3,5/

COPY THE KNOWN PRIMES TO THE OUTPUT LINE. 

     DO 1 I = 1,PP
   1   ALINE(I) = PRIME(I)
     LL = PP
     LP = 3
     XP2 = PRIME(LP + 1)**2
     N = 5
     INC = 4

CONSIDER ANOTHER CANDIDATE. VIA INC, DODGE MULTIPLES OF 2 AND 3. 

  10 INC = 6 - INC
     N = N + INC
     IF (N - XP2) 20,11,20
  11 LP = LP + 1
     XP2 = PRIME(LP + 1)**2
     GO TO 40

CHECK SUCCESSIVE PRIMES AS FACTORS, STARTING WITH PRIME(4) = 5. 

  20 IP = 4
  21 F = PRIME(IP)
     Q = N/F
     IF (Q*F - N) 22,40,22
  22 IP = IP + 1
     IF (IP - LP) 21,21,30

CAUGHT ANOTHER PRIME. 

  30 IF (PP - ENUFF) 31,32,32
  31 PP = PP + 1
     PRIME(PP) = N
  32 IF (LL - 20) 35,33,33
  33 WRITE (6,34) (ALINE(I), I = 1,LL)
  34 FORMAT (20I6)
     LL = 0
  35 LL = LL + 1
     ALINE(LL) = N

COMPLETED? 

  40 IF (N - 32767) 10,41,41
  41 WRITE (6,34) (ALINE(I), I = 1,LL)
     END</lang>

Start of output: 

    1     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   547   557   563   569   571   577   587   593   599   601   607   613   617   619   631   641   643   647   653
  659   661   673   677   683   691   701   709   719   727   733   739   743   751   757   761   769   773   787   797
  809   811   821   823   827   829   839   853   857   859   863   877   881   883   887   907   911   919   929   937
  941   947   953   967   971   977   983   991   997  1009  1013  1019  1021  1031  1033  1039  1049  1051  1061  1063
 1069  1087  1091  1093  1097  1103  1109  1117  1123  1129  1151  1153  1163  1171  1181  1187  1193  1201  1213  1217

etc. and ends 

32423 32429 32441 32443 32467 32479 32491 32497 32503 32507 32531 32533 32537 32561 32563 32569 32573 32579 32587 32603
32609 32611 32621 32633 32647 32653 32687 32693 32707 32713 32717 32719 32749

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 <- 3
 q := 9
 ps := make(chan int)
 go Sieve(ps)                   // separate primes supply
 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)    // n == q == p*p
       go Filter(nums, ch1, p)  // creation of a filter by p, at p*p
       nums = ch1
   	p = <-ps                 // next prime
   	q = p*p                  //   and its square
   }
 }

}

func primes (c chan<- int) { 

 c <- 2
 go Sieve(c)

}

func main() { 

 ch := make(chan int)
 go primes(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 classic version has isPrime in the above inlined:

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

                                   True primes]</lang>

It is easy to amend it to test only odd numbers, by only odd primes:

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

                                          True (drop 1 primes)]</lang>

or automatically skip the multiples of 3 as well (a wheel factorization technique), with

<lang haskell>primes3 = 2 : 3 : 5 : [n | n <- scanl (+) 7 (cycle [4,2]), 

                                    foldr (\p r-> p*p > n || rem n p > 0 && r)
                                          True (drop 2 primes)]</lang>

Sieve by trial division

We can repeatedly sieve a stream of candidate numbers for non-divisibility by prime at a time, even if the stream is unbounded, thanks to lazy evaluation: <lang haskell>primesT = sieve [2..] 

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

-- map head -- . iterate (\(p:xs)-> filter ((>0).(`rem`p)) xs) $ [2..]</lang>

This is similar to David Turner's 1983 SASL code.

As shown in Melissa O'Neill's paper "The Genuine Sieve of Eratosthenes", its complexity is quadratic in number of primes produced whereas that of optimal trial division is , and of true SoE it is , in n primes produced.

Indeed as Eratosthenes sieve works by counting, its removal step could be prototyped as (\(p:xs)-> minus xs [p,p+p..]), where minus xs ys == xs \\ ys for any finite and increasing xs and ys.

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

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

-- (\(a,b:_) -> map head a ++ b) . span ((< m).(^2).head) -- $ iterate (\(p:xs) -> filter ((>0).(`rem`p)) xs) [2..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

-- fix $ (2:) . concatMap (fst.fst) -- . iterate (\((h,xs),p:t) -> (span (< head t^2) $ filter ((>0).(`rem`p)) xs, t)) -- . (,) ([3],[4..])</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 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.

Note: if you instead want an iterator, 4&p: y returns the next prime after y.

Julia

I've chosen to solve this task by creating a new iterator type, TDPrimes. TDPrimes contains the upper limit of the sequence. The iteration state is the list of computed primes, and the item returned with each iteration is the current prime. The core of the solution is the next method for TDPrimes, which computes the next prime by trial division of the previously determined primes contained in the iteration state. <lang Julia> type TDPrimes{T<:Integer} 

   plim::T

end

function Base.start{T<:Integer}(pl::TDPrimes{T}) 

   2ones(T, 1)

end

function Base.done{T<:Integer}(pl::TDPrimes{T}, p::Array{T,1}) 

   p[end] > pl.plim

end

function Base.next{T<:Integer}(pl::TDPrimes{T}, p::Array{T,1}) 

   pr = p[end]
   for i in (pr+1):(pl.plim)
       ispr = true
       for j in p
           if i%j == 0
               ispr = false
               break
           end
       end
       if ispr
           push!(p, i)
           return (pr, p)
       end
   end
   push!(p, typemax(T))
   return (pr, p)

end

n = 100 print("The primes <= ", n, " are:\n ")

for i in TDPrimes(n) 

   print(i, " ")

end println() </lang>

Output:
The primes <= 100 are:
    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

jq

Works with: jq version 1.4

This entry uses is_prime/0 as defined at Primality_by_trial_division#jq. <lang jq># Produce a (possibly empty) stream of primes in the range [m,n], i.e. m <= p <= n def primes(m; n): 

 ([m,2] | max) as $m
 | if $m > n then empty
   elif $m == 2 then 2, primes(3;n)
   else (1 + (2 * range($m/2 | floor; (n + 1) /2 | floor))) | select( is_prime )
   end;</lang>

Examples: <lang jq>primes(0;10)</lang> <lang sh>2 3 5 7</lang> Produce an array of primes, p, satisfying 50 <= p <= 99: <lang jq>[primes(50;99)]</lang> 

[53,59,61,67,71,73,79,83,89,97]

Lua

<lang Lua> function isPrime (n) -- Function to test primality by trial division if n < 2 then return false end if n < 4 then return true end if n % 2 == 0 then return false end for d = 3, math.sqrt(n), 2 do if n % d == 0 then return false end end return true end

local start = arg[1] or 1 -- Default limits may be overridden with local bound = arg[2] or 100 -- command line arguments, for example: -- lua primes.lua 1000 2000 for i = start, bound do if isPrime(i) then io.write(i .. " ") end end </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

MATLAB

<lang MATLAB>function primeList = sieveOfEratosthenes(lastNumber) 

   list = (2:lastNumber); %Construct list of numbers
   primeList = []; %Preallocate prime list
   
   while( list(1)^2 <lastNumber )
       
       primeList = [primeList list(1)]; %add prime to the prime list
       list( mod(list,list(1))==0 ) = []; %filter out all multiples of the current prime
       
   end
   
   primeList = [primeList list]; %The rest of the numbers in the list are primes

end</lang>

Sample Output:
sieveOfEratosthenes(30)

ans =

    2     3     5     7    11    13    17    19    23    29

Oforth

isPrime function is from Primality by trial division page

<lang Oforth>: primeSeq(n) { n seq filter(#isPrime) }</lang>

PARI/GP

<lang parigp>trial(n)={ 

 if(n < 4, return(n > 1)); /* Handle negatives */
 forprime(p=2,sqrt(n),
   if(n%p == 0, return(0))
 );
 1

};

select(trial, [1..100])</lang>


Pascal

Library: primTrial
Works with: Free Pascal

Hiding the work in a existing unit. <lang Pascal> program PrimeRng; uses 

 primTrial;

var 

 Range : ptPrimeList;
 i : integer;

Begin 

 Range := PrimeRange(1000*1000*1000,1000*1000*1000+100);
 For i := Low(Range) to High(Range) do
   write(Range[i]:12);
 writeln;

end.</lang>

output
  1000000007  1000000009  1000000021  1000000033  1000000087  1000000093  1000000097

Perl

<lang perl>sub isprime { 

 my $n = shift;
 return ($n >= 2) if $n < 4;
 return unless $n % 2 && $n % 3;
 my $sqrtn = int(sqrt($n));
 for (my $i = 5; $i <= $sqrtn; $i += 6) {
   return unless $n % $i && $n % ($i+2);
 }
 1;

}

print join(" ", grep { isprime($_) } 0 .. 100 ), "\n"; print join(" ", grep { isprime($_) } 12345678 .. 12345678+100 ), "\n";</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
12345701 12345709 12345713 12345727 12345731 12345743 12345769

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

PowerShell

Works with: PowerShell version 4.0

<lang PowerShell> function eratosthenes ($n) { 

   if($n -ge 1){
       $prime = @(1..($n+1) | foreach{$true})
       $prime[1] = $false
       $m = [Math]::Floor([Math]::Sqrt($n))
       for($i = 2; $i -le $m; $i++) {
           if($prime[$i]) {
               for($j = $i*$i; $j -le $n; $j += $i) {
                   $prime[$j] = $false
               }
           }
       }
       1..$n | where{$prime[$_]}
   } else {
       "$n must be equal or greater than 1"
   }

} function sieve-start-end ($start,$end) { 

   (eratosthenes $end) | where{$_ -ge $start}

} "$(sieve-start-end 100 200)" </lang> Output: 

101 103 107 109 113 127 131 137 139 149 151 157 163 167 173 179 181 191 193 197 199


Python

Using the basic prime() function from: "Primality by trial division" <lang Python> def prime(a): 

   return not (a < 2 or any(a % x == 0 for x in xrange(2, int(a**0.5) + 1)))

def primes_below(n): 

   return [i for i in range(n) if prime(i)]

</lang>

Output:
>>> primes_below(100)
[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]

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 division.*/ parse arg n . /*get optional number of primes to find*/ if n== then n=26 /*Not specified? Then assume default.*/ tell=n>0; n=abs(n) /*N is negative? Then don't display. */ @.1=2; @.2=3; @.3=5; @.4=7; @.5=11; @.6=13; #=5; s=@.#+2 

                                      /*   [↑]  is the number of low primes. */
     do p=1  for #  while  p<=n       /* [↓]  find primes, but don't show 'em*/
     if tell  then say right(@.p,9)   /*display some pre-defined low primes. */
     !.p=@.p**2                       /*also compute the squared value of P. */
     end   /*p*/                      /* [↑]  allows faster loop (below).    */
                                      /* [↓]  N:  default lists up to 101 #s.*/
  do j=s  by 2  while  #<n            /*start with the next odd prime.       */
  if j//  3    ==0   then iterate     /*is this number a triple?             */
  parse var j  -1 _                 /*obtain the last digit of the  J  var.*/
  if _         ==5   then iterate     /*is this number a 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 prev. prime?  ¬prime     ___    */
         end   /*j*/                  /* [↑]   only divide up to     √ J     */
  #=#+1                               /*bump the counter of number of primes.*/
  @.#=j;   !.#=j*j                    /*define this prime; define its square.*/
  if tell  then say right(j, 9)       /*maybe display this prime ──► terminal*/
  end     /*j*/                       /* [↑]  only display N number of primes*/
                                      /* [↓]  display number of primes found.*/

say # ' primes found.' /*stick a fork in it, we're all done. */</lang> output is the same as the 1st 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.

Sidef

Using the is_prime() function from: "Primality by trial division" <lang ruby>func prime_seq(amount, callback) { 

   var (counter, number) = (0, 0);
   while (counter < amount) {
       if (is_prime(number)) {
           callback(number);
           ++counter;
       }
       ++number;
   }

}

prime_seq(100, {|p| say p}); # prints the first 100 primes</lang>

Swift

<lang swift>import Foundation

extension SequenceType { 

 func takeWhile(include: Generator.Element -> Bool) -> AnyGenerator<Generator.Element> {
   var g = self.generate()
   return anyGenerator { g.next().flatMap{include($0) ? $0 : nil }}
 }

}

var pastPrimes = [2]

var primes = anyGenerator { 

 _ -> Int? in
 defer {
   pastPrimes.append(pastPrimes.last!)
   let c = pastPrimes.count - 1
   for p in anyGenerator({++pastPrimes[c]}) {
     let lim = Int(sqrt(Double(p)))
     if (!pastPrimes.takeWhile{$0 <= lim}.contains{p % $0 == 0}) { break }
   }
 }
 return pastPrimes.last

}</lang>

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

lang Eiffel> class SEQUENCE_OF_PRIMES_BY_TRIAL_DIVISIION feature sequence(lower, upper: INTEGER) require lower_positive: lower >0 upper_positive: upper >0 lower_smaller: lower

Retrieved from "https://rosettacode.org/wiki/Sequence_of_primes_by_trial_division?oldid=250143"