The sieve of Sundaram: Difference between revisions
m
→{{header|EasyLang}}
(Added AppleScript.) |
|||
| (41 intermediate revisions by 21 users not shown) | |||
Line 1:
{{
The sieve of Eratosthenes: you've been there; done that; have the T-shirt. The sieve of Eratosthenes was ancient history when Euclid was a schoolboy. You are ready for something less than 3000 years old. You are ready for The sieve of Sundaram.
Line 22:
;References:
* The article on [[wp:Sieve_of_Sundaram|Wikipedia]].
'''Comment on the Sundaram Sieve'''
In case casual readers and programmers read the above blurb and get the impression that something several thousand years newer must needs be better than the "old" Sieve of Eratosthenes (SoE), do note the only difference between the Sieve of Sundaram (SoS) and the odds-only SoE is that the SoS marks as composite/"culls" according to '''all odd "base" numbers''' as is quite clear in the above description of how to implement it and the above linked Wikipedia article (updated), and the SoE marks as composite/"culls" according to only the previously determined unmarked primes (which are all odd except for two, which is not used for the "odds-only" algorithm); the time complexity (which relates to the execution time) is therefore O(n log n) for the SoS and O(n log log n) for the SoE, which difference can make a huge difference to the time it takes to sieve as the ranges get larger. It takes about a billion "culls" to sieve odds-only to a billion for the SoE, whereas it takes about 2.28 billion "culls" to cull to the same range for the SoS, which implies that the SoS must be about this ratio slower for this range with the memory usage identical. Why would one choose the SoS over the SoE to save a single line of code at the cost of this much extra time? The Wren comparison at the bottom of this page makes that clear, as would implementing the same in any language.
<br><br>
=={{header|11l}}==
{{trans|Python}}
<syntaxhighlight lang="11l">F sieve_of_Sundaram(nth, print_all = 1B)
‘
The sieve of Sundaram is a simple deterministic algorithm for finding all the
prime numbers up to a specified integer. This function is modified from the
Wikipedia entry wiki/Sieve_of_Sundaram, to give primes to their nth rather
than the Wikipedia function that gives primes less than n.
’
assert(nth > 0, ‘nth must be a positive integer’)
V k = Int((2.4 * nth * log(nth)) I/ 2)
V integers_list = [1B] * k
L(i) 1 .< k
V j = Int64(i)
L i + j + 2 * i * j < k
integers_list[Int(i + j + 2 * i * j)] = 0B
j++
V pcount = 0
L(i) 1 .. k
I integers_list[i]
pcount++
I print_all
print(f:‘{2 * i + 1:4}’, end' ‘ ’)
I pcount % 10 == 0
print()
I pcount == nth
print("\nSundaram primes start with 3. The "nth‘th Sundaram prime is ’(2 * i + 1)".\n")
L.break
sieve_of_Sundaram(100, 1B)
sieve_of_Sundaram(1000000, 0B)</syntaxhighlight>
{{out}}
<pre>
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
Sundaram primes start with 3. The 100th Sundaram prime is 547.
Sundaram primes start with 3. The 1000000th Sundaram prime is 15485867.
</pre>
=={{header|ALGOL 68}}==
{{Trans|Nim}}
To run this with Algol 68G, you will need to increase the heap size by specifying e.g. <code>-heap=64M</code> on the command line.
<
INT n = 8 000 000;
INT none = 0, mark1 = 1, mark2 = 2;
Line 68 ⟶ 128:
print( ( "The millionth Sundaram prime is ", whole( last, 0 ), newline ) )
FI
END</
{{out}}
<pre>
Line 83 ⟶ 143:
The millionth Sundaram prime is 15485867
</pre>
=={{header|Amazing Hopper}}==
{{trans|C}}
<syntaxhighlight lang="Amazing Hopper">
#include <jambo.h>
Main
Set break
tiempo inicio = 0, tiempo final = 0
nprimes=1000000, nmax=0
Let ( nmax := Ceil( Mul( nprimes, Sub( Add(Log(nprimes), Log(Log(nprimes))), 0.9385) ) ) )
k=0, Let( k := Div( Minus two 'nmax', 2) )
a=0
Set decimal '0'
Seqspaced(3, {k} Mul by '2' Plus '1', {k} Mul by '2' Plus '1' Div into '2', a)
Unset decimal
i=1
pos inicial sumando=2, pos ini factor = 2, factor factor = 2, suma = 6
sumando = 0, factor = 0
end subloop=0
Tic( tiempo inicio )
Loop
/* calculo las secuencias para las posiciones; ocupa la memoria creada para la primera secuencia,
es mucho, pero si lo hago con ciclos, el loop termina dentro de 6 minutos :D */
Let ( end subloop := {k} Minus 'i', {i} Mul by '2' Plus '1', Div it )
Sequence( pos inicial sumando, 1, end subloop, sumando )
Sequence( pos ini factor, factor factor, end subloop, factor )
Let ( sumando := Add( sumando, factor) )
Set range 'sumando', Set '0', Put 'a' // pongo ceros en las posiciones calculadas
Clr range
/* recalculo índices para nuevas posiciones */
pos inicial sumando += 2 // 2,4,6,8...
pos ini factor += suma // 2, 8, 18, 32
suma += 4 // 10, 14, 18
factor factor += 2 // 2,4,6,8
++i
While ( Less equal ( Mul( Mul( Plus one (i),i ),2), k ) )
Toc( tiempo inicio, tiempo final )
/* Visualización de los primeros 100 primos. Esto podría hacerlo con ciclos,
como lo hace la versión de "C", pero me gusta disparar moscas con un rifle */
Cls
ta=0, Compact 'a', Move to 'a' // elimino los ceros = compacto array
[1:100] Get 'a', Move to 'ta', Redim (ta, 10, 10)
Tok sep ("\t"), Print table 'ta'
Clr interval, Clear 'ta'
/* imprimo el primo número "nprimes" */
Print( nprimes, " th Sundaram prime is ", [ nprimes ] Get 'a', "\n" )
Printnl( "Time = ", tiempo final, " segs" )
End
</syntaxhighlight>
{{out}}
<pre>
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
1000000 th Sundaram prime is 15485867
Time = 9.9078 segs
/* Sí, mi lenguaje es lento para algunas cosas... */
</pre>
<p>Version 2</p>
<p>Reescribí el programa para ver si podía reducir el tiempo de ejecución, y logré rducirlo a unos 5 segundos aproximados en una máquina cuántica. :D </p>
<syntaxhighlight lang="Amazing Hopper">
#include <jambo.h>
Main
tiempo inicio = 0, tiempo final = 0
nprimes=1000000,
nmax=0
Let ( nmax := Ceil( Mul( nprimes, Sub( Add(Log(nprimes), Log(Log(nprimes))), 0.9385) ) ) )
k=0
Let( k := Div( Minus two 'nmax', 2) )
a=0
Set decimal '0'
Seqspaced(3, {k} Mul by '2' Plus '1', {k} Mul by '2' Plus '1' Div into '2', a)
Unset decimal
pos inicial sumando=2
pos ini factor = 2, suma = 6
end subloop=0
i=1
Tic( tiempo inicio )
Loop
Let ( end subloop := 'k' Minus 'i'; 'i' Mul by '2' Plus '1', Div it )
Get sequence( pos inicial sumando, 1, end subloop )
Get sequence( pos ini factor, pos inicial sumando, end subloop )
---Add it---
Get range // usa el rango desde el stack. Se espera que el rango sea una variable,
// por lo que no se quitará desde la memoria hasta un kill (Forget en Jambo)
Set '0', Put 'a'
--- Forget --- // para quitar el rango desde el stack.
pos inicial sumando += 2 // 2,4,6,8...
pos ini factor += suma // 2, 8, 18, 32
suma += 4 // 10, 14, 18
++i
While ( Less equal ( Mul( Mul( Plus one (i),i ),2), k ) )
Toc( tiempo inicio, tiempo final )
Clr range
/* Visualización */
Cls
ta=0, [1:100] Move positives from 'a' Into 'ta'
Redim (ta, 10, 10)
Tok sep ("\t"), Print table 'ta'
Clr interval, Clear 'ta'
/* imprimo el primo número "nprimes" */
Setdecimal(0)
Print( nprimes, " th Sundaram prime is ", [ nprimes ] Get positives from 'a' , "\n" )
Printnl( "Time = ", tiempo final, " segs" )
End
</syntaxhighlight>
{{out}}
<pre>
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
1000000 th Sundaram prime is 15485867
Time = 4.9630 segs
/* Sí, mi lenguaje sigue siendo lento para algunas cosas... */
</pre>
=={{header|AppleScript}}==
The "nth prime" calculation here's gleaned from the Python and Julia solutions and the limitations to marking partly from the Phix.
<syntaxhighlight lang="applescript">on sieveOfSundaram(indexRange)
if (indexRange's class is list) then
set n1 to beginning of indexRange
set n2 to end of indexRange
else
set n1 to
set n2 to
end if
Line 99 ⟶ 317:
end script
set {unmarked, marked} to {true, false}
--
-- 1-based indices. The Python and Julia solutions note that the nth prime is approximately
-- n * 1.2 * log(n), but the number from which it'll be derived is only about half that.
-- 15 is added too here to ensure headroom with lower prime counts.
set limit to (do shell script "echo '" & n2 & " * 0.6 * l(" & n2 & ") + 15'| bc -l") as integer
set len to 1500
repeat len times
Line 115 ⟶ 332:
end repeat
-- Since it's a given that every third slot from 4 on will be "marked" (changed to false), there'll be
-- no need to check these and thus no point in actually marking them! Skip the step = 3 marking sweep
-- and the first slot of every three for marking in the subsequent sweeps.
repeat with
repeat with j from (step * step div 2 - (step * 2 mod 3) * step + step) to (limit - step) by (step * 3)
set item j of o's lst to marked
set item (j + step) of o's lst to marked
end repeat
end repeat
-- Calculate the primes from the indices of the unmarked slots
-- and store them in the list from the beginning.
set i to 1
set item i of o's lst to i * 2 + 1
repeat with n from 2 to limit by 3
if (item n of o's lst) then
set i to i + 1
set item i of o's lst to n * 2 + 1
if (i = n2) then exit repeat
end if
if (item (n + 1) of o's lst) then
set i to i + 1
set item i of o's lst to n * 2 + 3 -- ((n + 1) * 2) + 1)
if (i = n2) then exit repeat
end if
end repeat
-- set beginning of o's lst to 2 -- Uncomment if required.
return items n1 thru n2 of o's lst
end sieveOfSundaram
-- Task code:
on join(lst, delim)
set astid to AppleScript's text item delimiters
Line 155 ⟶ 390:
end task
task()</
{{output}}
<
3 5 7 11 13 17 19 23 29 31
37 41 43 47 53 59 61 67 71 73
Line 170 ⟶ 405:
479 487 491 499 503 509 521 523 541 547
1,000,000th:
15485867"</
=={{header|C}}==
<syntaxhighlight lang="c">
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
int main(void) {
int nprimes = 1000000;
int nmax = ceil(nprimes*(log(nprimes)+log(log(nprimes))-0.9385));
// should be larger than the last prime wanted; See
// https://www.maa.org/sites/default/files/jaroma03200545640.pdf
int i, j, m, k; int *a;
k = (nmax-2)/2;
a = (int *)calloc(k + 1, sizeof(int));
for(i = 0; i <= k; i++)a[i] = 2*i+1;
for (i = 1; (i+1)*i*2 <= k; i++)
for (j = i; j <= (k-i)/(2*i+1); j++) {
m = i + j + 2*i*j;
if(a[m]) a[m] = 0;
}
for (i = 1, j = 0; i <= k; i++)
if (a[i]) {
if(j%10 == 0 && j <= 100)printf("\n");
j++;
if(j <= 100)printf("%3d ", a[i]);
else if(j == nprimes){
printf("\n%d th prime is %d\n",j,a[i]);
break;
}
}
}
</syntaxhighlight>
{{out}}
<pre>
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
1000000 th prime is 15485867
</pre>
=={{header|C#|CSharp}}==
<strike>Generating prime numbers during sieve creation gives a performance boost over completing the sieve and then scanning the sieve for output. There are, of course, a few at the end to scan out.</strike><br/>
Heh, nope. It's faster to do the sieving first, then the generation afterwards.
<
using System.Collections.Generic;
using System.Linq;
Line 217 ⟶ 500:
yield return (v << 1) + 1;
}
}</
{{out}}
Under <strike>1/5</strike> <b><font color=darkgreen>1/8</font></b> of a second @ Tio.run
Line 234 ⟶ 517:
The millionth odd prime number: 15485867
124.8262 ms</pre>P.S. for those (possibly faithless) who wish to have a conventional prime number generator, one can uncomment the <code>yield return 2</code> line at the top of the function.
=={{header|C++}}==
<syntaxhighlight lang="c++">
#include <cmath>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <vector>
std::vector<uint32_t> sieve_of_sundaram(const uint32_t& limit) {
std::vector<uint32_t> primes = {};
if ( limit < 3 ) {
return primes;
}
const uint32_t k = ( limit - 3 ) / 2 + 1;
std::vector<bool> marked(k, true);
for ( uint32_t i = 0; i < ( std::sqrt(limit) - 3 ) / 2 + 1; ++i ) {
uint32_t p = 2 * i + 3;
uint32_t s = ( p * p - 3 ) / 2;
for ( uint32_t j = s; j < k; j += p ) {
marked[j] = false;
}
}
for ( uint32_t i = 0; i < k; ++i ) {
if ( marked[i] ) {
primes.emplace_back(2 * i + 3);
}
}
return primes;
}
int main() {
std::vector<uint32_t> primes = sieve_of_sundaram(16'000'000);
std::cout << "The first 100 odd primes generated by the Sieve of Sundaram:" << std::endl;
for ( uint32_t i = 0; i < 100; ++i ) {
std::cout << std::setw(3) << primes[i] << ( i % 10 == 9 ? "\n" :" " );
}
std::cout << "\n" << "The 1_000_000th Sundaram prime is " << primes[1'000'000 - 1] << std::endl;
}
</syntaxhighlight>
{{ out }}
<pre>
The first 100 odd primes generated by the Sieve of Sundaram:
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
The 1_000_000th Sundaram prime is 15485867
</pre>
=={{header|EasyLang}}==
<syntaxhighlight>
func log n .
return log10 n / log10 2.71828182845904523
.
proc sundaram np . primes[] .
nmax = floor (np * (log np + log log np) - 0.9385) + 1
k = (nmax - 2) / 2
len marked[] k
for i to k
h = 2 * i + 2 * i * i
while h <= k
marked[h] = 1
h += 2 * i + 1
.
.
i = 1
primes[] = [ ]
while np > 0
if marked[i] = 0
np -= 1
primes[] &= 2 * i + 1
.
i += 1
.
.
sundaram 100 primes[]
print primes[]
sundaram 1000000 primes[]
print primes[len primes[]]
</syntaxhighlight>
{{out}}
<pre>
[ 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 ]
15485867
</pre>
=={{header|F_Sharp|F#}}==
<
// The sieve of Sundaram. Nigel Galloway: August 7th., 2021
let sPrimes()=
Line 249 ⟶ 628:
sPrimes()|>Seq.take 100|>Seq.iter(printf "%d "); printfn ""
printfn "The millionth Sundaram prime is %d" (Seq.item 999999 (sPrimes()))
</syntaxhighlight>
{{out}}
<pre>
Line 255 ⟶ 634:
The millionth Sundaram prime is 15485867
</pre>
=={{header|Fortran}}==
<syntaxhighlight lang="fortran">
PROGRAM SUNDARAM
IMPLICIT NONE
!
! Local variables
!
INTEGER(8) :: curr_index
INTEGER(8) :: i
INTEGER(8) :: j
INTEGER :: lim
INTEGER(8) :: mid
INTEGER :: primcount
LOGICAL*1 , ALLOCATABLE , DIMENSION(:) :: primes !Array of booleans representing integers
lim = 10000000 ! Not the number of primes but the storage where the prime marker is held for the millionth prime
ALLOCATE(primes(lim))
primes(1:lim) = .TRUE.
!Set all to .True., we will later block out the known non-primes
mid = lim/2
!Generate primes
DO j = 1 , mid
DO i = 1 , j
curr_index = i + j + (2*i*j)
IF( curr_index>lim )EXIT ! Too big already, leave the loop.
primes(curr_index) = .FALSE. !This candidate will not produce a prime
END DO
END DO
!
i = 0
j = 0
WRITE(6 , *)'The first 100 primes:'
DO WHILE ( i < 100 )
j = j + 1
IF( primes(j) )THEN
WRITE(6 , 34 , ADVANCE = 'no')j*2 + 1 !Take the candidate, multiply by 2, add 1, and you have a prime
34 FORMAT(I0 , 1x)
i = i + 1 ! Counter used for printing
IF( MOD(i,10)==0 )WRITE(6 , *)' '
END IF
END DO
! Now print the millionth prime
primcount = 0
DO i = 1 , lim
IF( primes(i) )THEN
primcount = primcount + 1
IF( primcount==1000000 )THEN
WRITE(6 , 35)'1 millionth Prime Found: ' , (i*2) + 1
35 FORMAT(/ , a , i0)
EXIT
END IF
END IF
END DO
DEALLOCATE(primes)
STOP
END PROGRAM SUNDARAM
</syntaxhighlight>
{{out}}
<pre>
The first 100 primes:
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
1 millionth Prime Found: 15485867
</pre>
=={{header|FreeBASIC}}==
<syntaxhighlight lang="vbnet">Function sieve_of_Sundaram(n As Uinteger) As Uinteger Ptr
If n < 3 Then Return 0
Dim As Uinteger r = Cint(Sqr(n))
Dim As Uinteger k = Cint((n - 3) / 2) + 1
Dim As Uinteger l = Cint((r - 3) / 2) + 1
Dim As Uinteger Ptr primes = Callocate(k, Sizeof(Uinteger))
Dim As Boolean Ptr marked = Callocate(k, Sizeof(Boolean))
For i As Uinteger = 1 To l
Dim As Uinteger p = 2 * i + 1
Dim As Uinteger s = Cint((p * p - 1) / 2)
For j As Uinteger = s To k Step p
marked[j] = True
Next j
Next i
Dim As Uinteger count = 0
For i As Uinteger = 1 To k
If Not marked[i] Then
primes[count] = 2 * i + 1
count += 1
End If
Next i
Return primes
End Function
Const limit As Uinteger = 16e6
Dim As Uinteger Ptr s = sieve_of_Sundaram(limit)
Print "First 100 odd primes generated by the Sieve of Sundaram:"
For i As Uinteger = 0 To 99
Print Using "#####"; s[i];
If (i + 1) Mod 10 = 0 Then Print
Next i
Print !"\nSundaram primes start with 3."
Print !"\nThe 100th Sundaram prime is: "; s[99]
Print !"\nThe 1000000th Sundaram prime is: "; s[999999]
Sleep</syntaxhighlight>
{{out}}
<pre>First 100 odd primes generated by the Sieve of Sundaram:
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
Sundaram primes start with 3.
The 100th Sundaram prime is: 547
The 1000000th Sundaram prime is: 15485867</pre>
=={{header|Go}}==
{{trans|Wren}}
{{libheader|Go-rcu}}
<
import (
Line 343 ⟶ 853:
fmt.Printf("\nUsing the Sieve of Eratosthenes would have generated them in %s ms.\n", celapsed)
fmt.Printf("\nAs a check, the %s Sundaram prime would again have been %s\n", million, millionth)
}</
{{out}}
Line 369 ⟶ 879:
=={{header|Haskell}}==
<
import Data.List.Split (chunksOf)
import qualified Data.Set as S
Line 411 ⟶ 921:
let ws = maximum . fmap length <$> transpose rows
pw = printf . flip intercalate ["%", "s"] . show
in unlines $ intercalate gap . zipWith pw ws <$> rows</
{{Out}}
<pre>First 100 Sundaram primes (starting at 3):
Line 425 ⟶ 935:
421 431 433 439 443 449 457 461 463 467
479 487 491 499 503 509 521 523 541 547</pre>
=={{header|J}}==
Loosely based on the perl implementation:
<syntaxhighlight lang=J>sundaram=: {{
sieve=. -.1{.~k=. <.1.2*(*^.) y
for_i. 1+i.y do.
f=. 1+2*i
j=. (#~ k > ]) (i,f) p. i+i.<.k%f
if. 0=#j do. y{.1+2*I. sieve return. end.
sieve=. 0 j} sieve
end.
}}</syntaxhighlight>
Task examples:
<syntaxhighlight lang=J>
P=: sundaram 1e6
10 10$P
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
{:P
15485867
</syntaxhighlight>
=={{header|Java}}==
<syntaxhighlight lang="java>
import java.util.ArrayList;
import java.util.List;
public final class TheSieveOfSundaram {
public static void main(String[] args) {
List<Integer> primes = sieveOfSundaram(16_000_000);
System.out.println("The first 100 odd primes generated by the Sieve of Sundaram:");
for ( int i = 0; i < 100; i++ ) {
System.out.print(String.format("%3d%s", primes.get(i), ( i % 10 == 9 ? "\n" :" " )));
}
System.out.println();
System.out.println("The 1_000_000th Sundaram prime is " + primes.get(1_000_000 - 1));
}
private static List<Integer> sieveOfSundaram(int limit) {
List<Integer> primes = new ArrayList<Integer>();
if ( limit < 3 ) {
return primes;
}
final int k = ( limit - 3 ) / 2 + 1;
boolean[] marked = new boolean[k];
for ( int i = 0; i < ( (int) Math.sqrt(limit) - 3 ) / 2 + 1; i++ ) {
int p = 2 * i + 3;
int s = ( p * p - 3 ) / 2;
for ( int j = s; j < k; j += p ) {
marked[j] = true;
}
}
for ( int i = 0; i < k; i++ ) {
if ( ! marked[i] ) {
primes.add(2 * i + 3);
}
}
return primes;
}
}
</syntaxhighlight>
{{ out }}
<pre>
The first 100 odd primes generated by the Sieve of Sundaram:
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
The 1_000_000th Sundaram prime is 15485867
</pre>
=={{header|JavaScript}}==
<syntaxhighlight lang="javascript">(() => {
"use strict";
// ----------------- SUNDARAM PRIMES -----------------
// sundaramsUpTo :: Int -> [Int]
const sundaramsUpTo = n => {
const
m = Math.floor(n - 1) / 2,
excluded = new Set(
enumFromTo(1)(
Math.floor(Math.sqrt(m / 2))
)
.flatMap(
i => enumFromTo(i)(
Math.floor((m - i) / (1 + (2 * i)))
)
.flatMap(
j => [(2 * i * j) + i + j]
)
)
);
return enumFromTo(1)(m).flatMap(
x => excluded.has(x) ? (
[]
) : [1 + (2 * x)]
);
};
// nSundaramsPrimes :: Int -> [Int]
const nSundaramPrimes = n =>
sundaramsUpTo(
// Probable limit
Math.floor((2.4 * n * Math.log(n)) / 2)
)
.slice(0, n);
// ---------------------- TEST -----------------------
const main = () => [
"First 100 Sundaram primes",
"(starting at 3):\n",
table(10)(" ")(
nSundaramPrimes(100)
.map(n => `${n}`)
)
].join("\n");
// --------------------- GENERIC ---------------------
// enumFromTo :: Int -> Int -> [Int]
const enumFromTo = m =>
n => Array.from({
length: 1 + n - m
}, (_, i) => m + i);
// --------------------- DISPLAY ---------------------
// chunksOf :: Int -> [a] -> [[a]]
const chunksOf = n => {
// xs split into sublists of length n.
// The last sublist will be short if n
// does not evenly divide the length of xs.
const go = xs => {
const chunk = xs.slice(0, n);
return 0 < chunk.length ? (
[chunk].concat(
go(xs.slice(n))
)
) : [];
};
return go;
};
// justifyRight :: Int -> Char -> String -> String
const justifyRight = n =>
// The string s, preceded by enough padding (with
// the character c) to reach the string length n.
c => s => Boolean(s) ? (
s.padStart(n, c)
) : "";
// table :: Int -> String -> [String] -> String
const table = nCols =>
// A tabulation of a list of values into a given
// number of columns, using a specified gap
// between those columns.
gap => xs => {
const w = xs[xs.length - 1].length;
return chunksOf(nCols)(xs)
.map(
row => row.map(
justifyRight(w)(" ")
).join(gap)
)
.join("\n");
};
return main();
})();</syntaxhighlight>
{{Out}}
<pre>First 100 Sundaram primes
(starting at 3):
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</pre>
=={{header|jq}}==
{{works with|jq}}
'''Works with gojq, the Go implementation of jq''' (*)
The hard part is anticipating how large the sieve must be to ensure
the n-th prime will be included.
Julia uses: (1.2 * nth * log(nth)) which is perhaps larger than always necessary.
Here we employ a naive adaptive approach.
(*) For large sieves, gojq will consume a very large amount of memory.
<syntaxhighlight lang="jq"># `sieve_of_Sundaram` as defined here generates the stream of
# consecutive primes from 3 on but less than or equal to the specified
# limit specified by `.`.
# input: an integer, n
# output: stream of consecutive primes from 3 but less than or equal to n
def sieve_of_Sundaram:
def idiv($b): (. - (. % $b))/$b ;
debug |
round as $n
| if $n < 2 then empty
else
((($n-3) | idiv(2)) + 1) as $k
| [range(0; $k + 1) | 1 ] # integers_list
| reduce range (0; (($n|sqrt) - 3) / 2 + 1) as $i (.;
(2*$i + 3) as $p
| ((($p*$p - 3) | idiv(2))) as $s
| reduce range($s; $k; $p) as $j (.;
if .[$j] then .[$j] = false else . end ) )
| range(0; $k) as $i
| if .[$i] then ($i+1)*2+1 else empty end
end ;
# Emit an array of $n Sundaram primes.
# The first Sundaram prime is 3 so we ensure Sundaram_prime(1) is [3].
# An adaptive definition to ensure generality without being excessively conservative.
def Sundaram_primes($n):
def sieve:
. as $in
| [limit($n; sieve_of_Sundaram)]
| if length == $n then .
else ($n + $in) as $m
| ("... nth_Sundaram_prime(\($n)): \($in) => \($m))" | debug) as $debug
| $m | sieve
end;
if $n < 1 then empty
elif $n <= 100 then ($n | 1.2 * . * log) | sieve
else $n | (1.15 * . * log) | sieve # OK
end;</syntaxhighlight>
'''For pretty-printing'''
<syntaxhighlight lang="jq">def lpad($len): tostring | ($len - length) as $l | (" " * $l)[:$l] + .;
def nwise($n):
def n: if length <= $n then . else .[0:$n] , (.[$n:] | n) end;
n;</syntaxhighlight>
'''The Tasks'''
<syntaxhighlight lang="jq">def hundred:
Sundaram_primes(100)
| nwise(10)
| map(lpad(3))
| join(" ");
"First hundred:", hundred,
"\nMillionth is \(Sundaram_primes(1000000)[-1])"</syntaxhighlight>
{{out}}
<pre>
First hundred:
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
Millionth is 15485867
</pre>
=={{header|Julia}}==
{{trans|Python}}
<
"""
The sieve of Sundaram is a simple deterministic algorithm for finding all the
Line 467 ⟶ 1,272:
using Primes; @show count(primesmask(15485867))
@time count(primesmask(15485867))
</
<pre>
3 5 7 11 13 17 19 23 29 31
Line 489 ⟶ 1,294:
0.022684 seconds (6 allocations: 5.785 MiB)
</pre>
=={{header|Mathematica}}/{{header|Wolfram Language}}==
<syntaxhighlight lang="mathematica">ClearAll[SieveOfSundaram]
SieveOfSundaram[n_Integer] := Module[{i, prefac, k, ints},
k = Floor[(n - 2)/2];
ints = ConstantArray[True, k + 1];
Do[
prefac = 2 i + 1;
If[i + i prefac <= k,
ints[[i + i prefac ;; ;; prefac]] = False
];
,
{i, 1, k + 1}
];
2 Flatten[Position[ints, True]] + 1
]
SieveOfSundaram[600][[;; 100]]
SieveOfSundaram[16000000][[10^6]]</syntaxhighlight>
{{out}}
<pre>{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}
15485867</pre>
=={{header|Nim}}==
<
const N = 8_000_000
Line 533 ⟶ 1,359:
if last == 0:
quit "Not enough values in sieve. Found only $#.".format(count), QuitFailure
echo "The millionth Sundaram prime is ", ($last).insertSep()</
{{out}}
Line 550 ⟶ 1,376:
The millionth Sundaram prime is 15_485_867</pre>
=={{header|Perl}}==
<syntaxhighlight lang="perl">use strict;
use warnings;
use feature 'say';
my @sieve;
my $nth = 1_000_000;
my $k = 2.4 * $nth * log($nth) / 2;
$sieve[$k] = 0;
for my $i (1 .. $k) {
my $j = $i;
while ((my $l = $i + $j + 2 * $i * $j) < $k) {
$sieve[$l] = 1;
$j++
}
}
$sieve[0] = 1;
my @S = (grep { $_ } map { ! $sieve[$_] and 1+$_*2 } 0..@sieve)[0..99];
say "First 100 Sundaram primes:\n" .
(sprintf "@{['%5d' x 100]}", @S) =~ s/(.{50})/$1\n/gr;
my ($count, $index);
for (@sieve) {
$count += !$_;
(say "One millionth: " . (1+2*$index)) and last if $count == $nth;
++$index;
}</syntaxhighlight>
{{out}}
<pre>First 100 Sundaram primes:
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
One millionth: 15485867</pre>
=={{header|Phix}}==
<!--<
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">sos</span><span style="color: #0000FF;">(</span><span style="color: #004080;">integer</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">)</span>
Line 579 ⟶ 1,448:
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"The first 100 odd prime numbers:\n%s\n"</span><span style="color: #0000FF;">,{</span><span style="color: #7060A8;">join_by</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">apply</span><span style="color: #0000FF;">(</span><span style="color: #004600;">true</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">sprintf</span><span style="color: #0000FF;">,{{</span><span style="color: #008000;">"%3d"</span><span style="color: #0000FF;">},</span><span style="color: #000000;">s</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">..</span><span style="color: #000000;">100</span><span style="color: #0000FF;">]}),</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">10</span><span style="color: #0000FF;">)})</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"The millionth odd prime number: %,d\n"</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">s</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1_000_000</span><span style="color: #0000FF;">]})</span>
<!--</
{{out}}
<pre>
Line 599 ⟶ 1,468:
=={{header|Python}}==
===Python :: Procedural===
<
def sieve_of_Sundaram(nth, print_all=True):
Line 634 ⟶ 1,503:
sieve_of_Sundaram(1000000, False)
</
<pre>
3 5 7 11 13 17 19 23 29 31
Line 654 ⟶ 1,523:
===Python :: Functional===
Composing functionally, and obtaining slightly better performance by defining a set (rather than list) of exclusions.
<
from math import floor, log, sqrt
Line 738 ⟶ 1,607:
# MAIN ---
if __name__ == '__main__':
main()</
{{Out}}
<pre>First hundred Sundaram primes, starting at 3:
Line 757 ⟶ 1,626:
15485867</pre>
=={{header|Racket}}==
{{trans|xxx}}
<syntaxhighlight lang="racket">#lang racket
(define (make-sieve-as-set limit)
(let ((marked (for/mutable-set ((i limit)) (add1 i))))
(let loop ((start 4) (step 3))
(cond [(>= start limit) marked]
[else (for ((i (in-range start limit step))) (set-remove! marked i))
(loop (+ start 3) (+ step 2))]))
(define (prime? n)
(and (odd? n)
(let ((idx (quotient (sub1 n) 2)))
(unless (<= idx limit) (error 'out-of-bounds))
(set-member? marked idx))))
(values marked prime?)))
(define (Sieve-of-Sundaram)
(define-values (sieve#1 prime?#1) (make-sieve-as-set 1000))
(displayln (for/list ((i 100) (p (sequence-filter prime?#1 (in-naturals)))) p))
;; this will generate primes *twice* as big, which should include 15485867...
(define-values (sieve#2 prime?#2) (make-sieve-as-set 10000000))
(define sorted-sieve#2 (sort (set->list sieve#2) <))
(displayln (add1 (* 2 (list-ref sorted-sieve#2 (sub1 1000000))))))
(module+ main
(Sieve-of-Sundaram))</syntaxhighlight>
{{out}}
<pre>(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)
15485867</pre>
=={{header|Raku}}==
<syntaxhighlight lang="raku" line>my $nth = 1_000_000;
my $k = Int.new: 2.4 * $nth * log($nth) / 2;
my int @sieve;
@sieve[$k] = 0;
race for (1 .. $k).batch(1000) -> @i {
for @i -> int $i {
my int $j = $i;
while (my int $l = $i + $j + 2 * $i * $j++) < $k {
@sieve[$l] = 1;
}
}
}
@sieve[0] = 1;
say "First 100 Sundaram primes:";
say @sieve.kv.map( { next if $^v; $^k * 2 + 1 } )[^100]».fmt("%4d").batch(10).join: "\n";
say "\nOne millionth:";
my ($count, $index);
for @sieve {
$count += !$_;
say $index * 2 + 1 and last if $count == $nth;
++$index;
}</syntaxhighlight>
{{out}}
<pre>First 100 Sundaram primes:
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
One millionth:
15485867
</pre>
=={{header|REXX}}==
For the calculation of the 1,000,000<sup>th</sup> Sundaram prime, it requires a 64-bit version of REXX.
<
parse arg n cols . /*get optional number of primes to find*/
if n=='' | n=="," then n= 100 /*Not specified? Then assume default.*/
Line 792 ⟶ 1,744:
exit 0 /*stick a fork in it, we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
commas: parse arg ?; do jc=length(?)-3 to 1 by -3; ?=insert(',', ?, jc); end; return ?</
{{out|output|text= when using the default inputs:}}
<pre>
Line 816 ⟶ 1,768:
</pre>
=={{header|Ruby}}==
Based on the Python code from the Wikipedia lemma.
<syntaxhighlight lang="ruby">def sieve_of_sundaram(upto)
n = (2.4 * upto * Math.log(upto)) / 2
k = (n - 3) / 2 + 1
bools = [true] * k
(0..(Integer.sqrt(n) - 3) / 2 + 1).each do |i|
p = 2*i + 3
s = (p*p - 3) / 2
(s..k).step(p){|j| bools[j] = false}
end
bools.filter_map.each_with_index {|b, i| (i + 1) * 2 + 1 if b }
end
p sieve_of_sundaram(100)
n = 1_000_000
puts "\nThe #{n}th sundaram prime is #{sieve_of_sundaram(n)[n-1]}"
</syntaxhighlight>
{{out}}
<pre>[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]
The 1000000th sundaram prime is 15485867
</pre>
=={{header|Wren}}==
{{libheader|Wren-fmt}}
I've worked here from the second (optimized) Python example in the Wikipedia article for SOS which allows an easy transition to an 'odds only' SOE for comparison.
<
var sos = Fn.new { |n|
Line 876 ⟶ 1,849:
Fmt.print("Using the Sieve of Sundaram generated primes up to $,d in $,d ms.\n", limit, elapsed)
System.print("First 100 odd primes generated by the Sieve of Sundaram:")
Fmt.print("\nThe $,d Sundaram prime is $,d", 1e6, primes[1e6-1])
Line 883 ⟶ 1,856:
elapsed = ((System.clock - start) * 1000).round
Fmt.print("\nUsing the Sieve of Eratosthenes would have generated them in $,d ms.", elapsed)
Fmt.print("\nAs a check, the $,d Sundaram prime would again have been $,d", 1e6, primes[1e6-1])</
{{out}}
| |||