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
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 // (* ****** ****** *) //
- include
"share/atspre_staload.hats" // (* ****** ****** *)
- define :: stream_cons
- define cons stream_cons
- 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>
D
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 SEQUENCE_OF_PRIMES_BY_TRIAL_DIVISIION feature sequence(lower, upper: INTEGER) 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 + " to " + upper.out + ".%N") i:= lower if i\\2=0 then i:=i+1 end from until i>upper loop if prime(i)=TRUE then io.put_integer (i) io.put_new_line end i:= i+2 end end feature{NONE} prime(n:INTEGER): BOOLEAN 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 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> Test: <lang Eiffel> class APPLICATION
inherit ARGUMENTS
create make
feature {NONE} -- Initialization
make local i: INTEGER do create prime_sequence prime_sequence.sequence (7, 99999) end prime_sequence: SEQUENCE_OF_PRIMES_BY_TRIAL_DIVISIION end </lang>
- Output:
Sequence of primes from 7 to 99999.) 7 11 13 . . . 99971 99989 99991
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
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 , and of true SoE it is , 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
-- primes = (2:) . concatMap (fst.snd) -- . iterate (\(p:t,(h,xs)) -> (t,span (< head t^2) [y | y <- xs, rem y p /= 0])) -- $ (primes, ([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
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.
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 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
Python
Straightforward naive implementation checking first if the number is even, then all the odd numbers below sqrt(n). The range is then filteres by the is_prime check
<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 is_prime(i)]
</lang>
Some testing:
- Output:
>>> primes_below(100) [1, 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 /*# 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 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.
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