The sieve of Sundaram: Difference between revisions

m
m (syntax highlighting fixup automation)
 
(15 intermediate revisions by 8 users not shown)
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>
 
Line 139 ⟶ 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>
 
Line 357 ⟶ 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#}}==
Line 453 ⟶ 709:
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}}==
Line 623 ⟶ 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}}==
Line 1,270 ⟶ 1,672:
@sieve[$k] = 0;
 
hyperrace for (1 .. $k).batch(1000) -> \@i {
for @i my-> int $j = i; {
while (my int $lj = i + $j + 2 * i * $j) < $k {;
while (my int $l = $i + $j + 2 * $i * $j++) < $k {
@sieve[$l] = 1;
$j = $j + 1;
}
}
}
 
Line 1,390 ⟶ 1,793:
=={{header|Wren}}==
{{libheader|Wren-fmt}}
{{libheader|Wren-seq}}
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.
<syntaxhighlight lang="ecmascriptwren">import "./fmt" for Fmt
import "/seq" for Lst
 
var sos = Fn.new { |n|
Line 1,448 ⟶ 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:")
for Fmt.tprint(chunk"$3d", in Lst.chunks(primes[0..99], 10)) Fmt.print("$3d", chunk)
Fmt.print("\nThe $,d Sundaram prime is $,d", 1e6, primes[1e6-1])
 
1,983

edits