The sieve of Sundaram

From Rosetta Code
Task
The sieve of Sundaram
You are encouraged to solve this task according to the task description, using any language you may know.

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.

Starting with the ordered set of +ve integers, mark every third starting at 4 (4;7;10...).

Step through the set and if the value is not marked output 2*n+1. So from 1 to 4 output 3 5 7.

4 is marked so skip for 5 and 6 output 11 and 13.

7 is marked, so no output but now also mark every fifth starting at 12 (12;17;22...)

as per to 10 and now mark every seventh starting at 17 (17;24;31....)

as per for every further third element (13;16;19...) mark every (9th;11th;13th;...) element.

The output will be the ordered set of odd primes.

Using your function find and output the first 100 and the millionth Sundaram prime.

The faithless amongst you may compare the results with those generated by The sieve of Eratosthenes.

References

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.

11l

Translation of: Python
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)
Output:
   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.

ALGOL 68

Translation of: Nim

To run this with Algol 68G, you will need to increase the heap size by specifying e.g. -heap=64M on the command line.

BEGIN # sieve of Sundaram #
    INT n = 8 000 000;
    INT none = 0, mark1 = 1, mark2 = 2;
    [ 1 : n ]INT mark;
    FOR i FROM LWB mark TO UPB mark DO mark[ i ] := none  OD;
    FOR i FROM   4 BY 3 TO UPB mark DO mark[ i ] := mark1 OD;

    INT            count := 0; # Count of primes.          #
    [ 1 : 100 ]INT list100;    # First 100 primes.         #
    INT            last  := 0; # Millionth prime.          #
    INT            step  := 5; # Current step for marking. #

    FOR i TO n WHILE last = 0 DO
        IF mark[ i ] = none THEN # Add/count a new odd prime.  #
            count +:= 1;
            IF   count <= 100 THEN
                list100[ count ] := 2 * i + 1
            ELIF count = 1 000 000 THEN
                last := 2 * i + 1
            FI
        ELIF mark[ i ] = mark1 THEN # Mark new numbers using current step. #
            IF i > 4 THEN
                FOR k FROM i + step BY step TO n DO
                    IF mark[ k ] = none THEN mark[ k ] := mark2 FI
                OD;
                step +:= 2
            FI
      # ELSE must be mark2 - Ignore this number. #
        FI
    OD; 

    print( ( "First 100 Sundaram primes:", newline ) );
    FOR i FROM LWB list100 TO UPB list100 DO
        print( ( whole( list100[ i ], -3 ) ) );
        IF i MOD 10 = 0 THEN print( ( newline ) ) ELSE print( ( " " ) ) FI
    OD;
    print( ( newline ) );
    IF last = 0 THEN
        print( ( "Not enough values in sieve. Found only ", whole( count, 0 ), newline ) )
    ELSE
        print( ( "The millionth Sundaram prime is ", whole( last, 0 ), newline ) )
    FI
END
Output:
  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 millionth Sundaram prime is 15485867

Amazing Hopper

Translation of: C
#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
Output:
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... */

Version 2

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

#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
Output:
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... */

AppleScript

The "nth prime" calculation here's gleaned from the Python and Julia solutions and the limitations to marking partly from the Phix.

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 indexRange
        set n2 to indexRange
    end if
    
    script o
        property lst : {}
    end script
    
    set {unmarked, marked} to {true, false}
    -- Build a list of 'true's corresponding to the unmarked start numbers implied by the
    -- 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
        set end of o's lst to unmarked
    end repeat
    repeat while (len < limit)
        set o's lst to o's lst & o's lst
        set len to len + len
    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 step from 5 to ((limit * 2) ^ 0.5 as integer) by 2
        -- Like the Phix solution, mark only from half the square of the step size, but adjusted
        -- to sync the repeat to the second slot in each group of three for marking.
        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
    set AppleScript's text item delimiters to delim
    set txt to lst as text
    set AppleScript's text item delimiters to astid
    return txt
end join

on task()
    --set r1 to sieveOfSundaram({1, 100})
    --set r2 to sieveOfSundaram(1000000)
    set r to sieveOfSundaram({1, 1000000})
    set r1 to items 1 thru 100 of r
    set r2 to item 1000000 of r
    set output to {"1st to 100th Sundaram primes:"}
    repeat with i from 1 to 100 by 10
        set end of output to join(items i thru (i + 9) of r1, "  ")
    end repeat
    set end of output to "1,000,000th: "
    set end of output to r2
    
    return join(output, linefeed)
end task

task()
Output:
"1st to 100th 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
1,000,000th: 
15485867"

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;
               }
           }
}
Output:
  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

C#

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.
Heh, nope. It's faster to do the sieving first, then the generation afterwards.

using System;
using System.Collections.Generic;
using System.Linq;
using static System.Console;

class Program
{
    static string fmt(int[] a)
    {
        var sb = new System.Text.StringBuilder();
        for (int i = 0; i < a.Length; i++)
            sb.Append(string.Format("{0,5}{1}",
              a[i], i % 10 == 9 ? "\n" : " "));
        return sb.ToString();
    }

    static void Main(string[] args)
    {
        var sw = System.Diagnostics.Stopwatch.StartNew();
        var pr = PG.Sundaram(15_500_000).Take(1_000_000).ToArray();
        sw.Stop();
        Write("The first 100 odd prime numbers:\n{0}\n",
          fmt(pr.Take(100).ToArray()));
        Write("The millionth odd prime number: {0}", pr.Last());
        Write("\n{0} ms", sw.Elapsed.TotalMilliseconds);
    }
}

class PG
{
    public static IEnumerable<int> Sundaram(int n)
    {
        // yield return 2;
        int i = 1, k = (n + 1) >> 1, t = 1, v = 1, d = 1, s = 1;
        var comps = new bool[k + 1];
        for (; t < k; t = ((++i + (s += d += 2)) << 1) - d - 2)
            while ((t += d + 2) < k)
                comps[t] = true;
        for (; v < k; v++)
            if (!comps[v])
                yield return (v << 1) + 1;
    }
}
Output:

Under 1/5 1/8 of a second @ Tio.run

The first 100 odd prime numbers:
    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 millionth odd prime number: 15485867
124.8262 ms
P.S. for those (possibly faithless) who wish to have a conventional prime number generator, one can uncomment the yield return 2 line at the top of the function.

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;
}
Output:
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

EasyLang

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[]]
Output:
[ 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

F#

// The sieve of Sundaram. Nigel Galloway: August 7th., 2021
let sPrimes()=
  let sSieve=System.Collections.Generic.Dictionary<int,(unit -> int) list>()
  let rec fN g=match g with h::t->(let n=h() in if sSieve.ContainsKey n then sSieve.[n]<-h::sSieve.[n] else sSieve.Add(n,[h])); fN t|_->()
  let     fI n=if sSieve.ContainsKey n then fN sSieve.[n]; sSieve.Remove n|>ignore; None else Some(2*n+1) 
  let     fG n g=let mutable n=n in (fun()->n<-n+g; n)
  let     fE n g=if not(sSieve.ContainsKey n) then sSieve.Add(n,[fG n g]) else sSieve.[n]<-(fG n g)::sSieve.[g] 
  let     fL    =let mutable n,g=4,3 in (fun()->n<-n+3; g<-g+2; fE (n+g) g; n)
  sSieve.Add(4,[fL]); Seq.initInfinite((+)1)|>Seq.choose fI

sPrimes()|>Seq.take 100|>Seq.iter(printf "%d "); printfn ""
printfn "The millionth Sundaram prime is %d" (Seq.item 999999 (sPrimes()))
Output:
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 millionth Sundaram prime is 15485867

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
Output:
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

FreeBASIC

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
Output:
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

Go

Translation of: Wren
Library: Go-rcu
package main

import (
    "fmt"
    "math"
    "rcu"
    "time"
)

func sos(n int) []int {
    if n < 3 {
        return []int{}
    }
    var primes []int
    k := (n-3)/2 + 1
    marked := make([]bool, k) // all false by default
    limit := (int(math.Sqrt(float64(n)))-3)/2 + 1
    for i := 0; i < limit; i++ {
        p := 2*i + 3
        s := (p*p - 3) / 2
        for j := s; j < k; j += p {
            marked[j] = true
        }
    }
    for i := 0; i < k; i++ {
        if !marked[i] {
            primes = append(primes, 2*i+3)
        }
    }
    return primes
}

// odds only
func soe(n int) []int {
    if n < 3 {
        return []int{}
    }
    var primes []int
    k := (n-3)/2 + 1
    marked := make([]bool, k) // all false by default
    limit := (int(math.Sqrt(float64(n)))-3)/2 + 1
    for i := 0; i < limit; i++ {
        if !marked[i] {
            p := 2*i + 3
            s := (p*p - 3) / 2
            for j := s; j < k; j += p {
                marked[j] = true
            }
        }
    }
    for i := 0; i < k; i++ {
        if !marked[i] {
            primes = append(primes, 2*i+3)
        }
    }
    return primes
}

func main() {
    const limit = int(16e6) // say
    start := time.Now()
    primes := sos(limit)
    elapsed := int(time.Since(start).Milliseconds())
    climit := rcu.Commatize(limit)
    celapsed := rcu.Commatize(elapsed)
    million := rcu.Commatize(1e6)
    millionth := rcu.Commatize(primes[1e6-1])
    fmt.Printf("Using the Sieve of Sundaram generated primes up to %s in %s ms.\n\n", climit, celapsed)
    fmt.Println("First 100 odd primes generated by the Sieve of Sundaram:")
    for i, p := range primes[0:100] {
        fmt.Printf("%3d ", p)
        if (i+1)%10 == 0 {
            fmt.Println()
        }
    }
    fmt.Printf("\nThe %s Sundaram prime is %s\n", million, millionth)

    start = time.Now()
    primes = soe(limit)
    elapsed = int(time.Since(start).Milliseconds())
    celapsed = rcu.Commatize(elapsed)
    millionth = rcu.Commatize(primes[1e6-1])
    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)
}
Output:
Using the Sieve of Sundaram generated primes up to 16,000,000 in 62 ms.

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,000 Sundaram prime is 15,485,867

Using the Sieve of Eratosthenes would have generated them in 33 ms.

As a check, the 1,000,000 Sundaram prime would again have been 15,485,867

Haskell

import Data.List (intercalate, transpose)
import Data.List.Split (chunksOf)
import qualified Data.Set as S
import Text.Printf (printf)

--------------------- SUNDARAM PRIMES --------------------

sundaram :: Integral a => a -> [a]
sundaram n =
  [ succ (2 * x)
    | x <- [1 .. m],
      x `S.notMember` excluded
  ]
  where
    m = div (pred n) 2
    excluded =
      S.fromList
        [ 2 * i * j + i + j
          | let fm = fromIntegral m,
            i <- [1 .. floor (sqrt (fm / 2))],
            let fi = fromIntegral i,
            j <- [i .. floor ((fm - fi) / succ (2 * fi))]
        ]

nSundaramPrimes ::
  (Integral a1, RealFrac a2, Floating a2) => a2 -> [a1]
nSundaramPrimes n =
  sundaram $ floor $ (2.4 * n * log n) / 2



--------------------------- TEST -------------------------
main :: IO ()
main = do
  putStrLn "First 100 Sundaram primes (starting at 3):\n"
  (putStrLn . table " " . chunksOf 10) $
    show <$> nSundaramPrimes 100

table :: String -> [[String]] -> String
table gap rows =
  let ws = maximum . fmap length <$> transpose rows
      pw = printf . flip intercalate ["%", "s"] . show
   in unlines $ intercalate gap . zipWith pw ws <$> rows
Output:
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

J

Loosely based on the perl implementation:

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

Task examples:

   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

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;
	}

}
Output:
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

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();
})();
Output:
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

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.

# `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;

For pretty-printing

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;

The Tasks

def hundred:
  Sundaram_primes(100)
  | nwise(10)
  | map(lpad(3))
  | join(" ");

"First hundred:", hundred,
"\nMillionth is \(Sundaram_primes(1000000)[-1])"
Output:
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


Julia

Translation of: Python
"""
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
Python example Wikipedia entry wiki/Sieve_of_Sundaram, to give primes to the
nth prime rather than the Wikipedia function that gives primes less than n.
"""
function sieve_of_Sundaram(nth, print_all=true)
    @assert nth > 0
    k = Int(round(1.2 * nth * log(nth)))  # nth prime is at about n * log(n)
    integers_list = trues(k)
    for i in 1:k
        j = i
        while i + j + 2 * i * j < k
            integers_list[i + j + 2 * i * j + 1] = false
            j += 1
        end
    end
    pcount = 0
    for i in 1:k + 1
        if integers_list[i + 1]
            pcount += 1
            if print_all
                print(lpad(2 * i + 1, 4), pcount % 10 == 0 ? "\n" : "")
            end
            if pcount == nth
                println("\nSundaram primes start with 3. The $(nth)th Sundaram prime is $(2 * i + 1).")
                break
            end
        end
    end
end

sieve_of_Sundaram(100)
@time sieve_of_Sundaram(1000000, false)

println("\nChecking:")
using Primes; @show count(primesmask(15485867))
@time count(primesmask(15485867))
Output:
  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.
  0.127168 seconds (19 allocations: 1.977 MiB)

Checking:
count(primesmask(15485867)) = 1000001
  0.022684 seconds (6 allocations: 5.785 MiB)

Mathematica/Wolfram Language

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]]
Output:
{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

Nim

import strutils

const N = 8_000_000

type Mark {.pure.} = enum None, Mark1, Mark2

var mark: array[1..N, Mark]
for n in countup(4, N, 3): mark[n] = Mark1


var count = 0           # Count of primes.
var list100: seq[int]   # First 100 primes.
var last = 0            # Millionth prime.
var step = 5            # Current step for marking.

for n in 1..N:
  case mark[n]
  of None:
    # Add/count a new odd prime.
    inc count
    if count <= 100:
      list100.add 2 * n + 1
    elif count == 1_000_000:
      last = 2 * n + 1
      break
  of Mark1:
    # Mark new numbers using current step.
    if n > 4:
      for k in countup(n + step, N, step):
        if mark[k] == None: mark[k] = Mark2
      inc step, 2
  of Mark2:
    # Ignore this number.
    discard


echo "First 100 Sundaram primes:"
for i, n in list100:
  stdout.write ($n).align(3), if (i + 1) mod 10 == 0: '\n' else: ' '
echo()
if last == 0:
  quit "Not enough values in sieve. Found only $#.".format(count), QuitFailure
echo "The millionth Sundaram prime is ", ($last).insertSep()
Output:
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

The millionth Sundaram prime is 15_485_867

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;
}
Output:
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

Phix

with javascript_semantics
function sos(integer n)
    if n<3 then return {} end if
    integer r = floor(sqrt(n)),
            k = floor((n-3)/2)+1,
            l = floor((r-3)/2)+1
    sequence primes = {},
             marked = repeat(false,k)
    for i=1 to l do
        integer p = 2*i+1,
                s = (p*p-1)/2
        for j=s to k by p do
            marked[j] = true
        end for
    end for
    for i=1 to k do
        if not marked[i] then
            primes = append(primes, 2*i+1)
        end if
    end for
    return primes
end function

sequence s = sos(16_000_000)
printf(1,"The first 100 odd prime numbers:\n%s\n",{join_by(apply(true,sprintf,{{"%3d"},s[1..100]}),1,10)})
printf(1,"The millionth odd prime number: %,d\n",{s[1_000_000]})
Output:
The first 100 odd prime numbers:
  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 millionth odd prime number: 15,485,867

Python

Python :: Procedural

from numpy import log

def sieve_of_Sundaram(nth, print_all=True):
    """
    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"
    k = int((2.4 * nth * log(nth)) // 2)  # nth prime is at about n * log(n)
    integers_list = [True] * k
    for i in range(1, k):
        j = i
        while i + j + 2 * i * j < k:
            integers_list[i + j + 2 * i * j] = False
            j += 1
    pcount = 0
    for i in range(1, k + 1):
        if integers_list[i]:
            pcount += 1
            if print_all:
                print(f"{2 * i + 1:4}", end=' ')
                if pcount % 10 == 0:
                    print()

            if pcount == nth:
                print(f"\nSundaram primes start with 3. The {nth}th Sundaram prime is {2 * i + 1}.\n")
                break



sieve_of_Sundaram(100, True)

sieve_of_Sundaram(1000000, False)
Output:
   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.

Python :: Functional

Composing functionally, and obtaining slightly better performance by defining a set (rather than list) of exclusions.

'''Sieve of Sundaram'''

from math import floor, log, sqrt
from itertools import islice


# sundaram :: Int -> [Int]
def sundaram(n):
    '''Sundaram prime numbers up to n'''
    m = (n - 1) // 2
    exclusions = {
        2 * i * j + i + j
        for i in range(1, 1 + floor(sqrt(m / 2)))
        for j in range(
            i, 1 + floor((m - i) / (1 + (2 * i)))
        )
    }
    return [
        1 + (2 * x) for x in range(1, 1 + m)
        if not x in exclusions
    ]


# nPrimesBySundaram :: Int -> [Int]
def nPrimesBySundaram(n):
    '''First n primes, by sieve of Sundaram.
    '''
    return list(islice(
        sundaram(
            # Probable limit
            int((2.4 * n * log(n)) // 2)
        ),
        int(n)
    ))


# ------------------------- TEST -------------------------
# main :: IO ()
def main():
    '''First 100 Sundaram primes, 
       and millionth Sundaram prime.
    '''
    print("First hundred Sundaram primes, starting at 3:\n")
    print(table(10)([
        str(s) for s in nPrimesBySundaram(100)
    ]))
    print("\n\nMillionth Sundaram prime, starting at 3:")
    print(
        f'\n\t{nPrimesBySundaram(1E6)[-1]}'
    )


# ----------------------- GENERIC ------------------------

# chunksOf :: Int -> [a] -> [[a]]
def chunksOf(n):
    '''A series of lists of length n, subdividing the
       contents of xs. Where the length of xs is not evenly
       divisible, the final list will be shorter than n.
    '''
    def go(xs):
        return (
            xs[i:n + i] for i in range(0, len(xs), n)
        ) if 0 < n else None
    return go


# table :: Int -> [String] -> String
def table(n):
    '''A list of strings formatted as
       right-justified rows of n columns.
    '''
    def go(xs):
        w = len(xs[-1])
        return '\n'.join(
            ' '.join(row) for row in chunksOf(n)([
                s.rjust(w, ' ') for s in xs
            ])
        )
    return go


# MAIN ---
if __name__ == '__main__':
    main()
Output:
First hundred 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


Millionth Sundaram prime, starting at 3:

    15485867

Racket

Translation of: xxx
#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))
Output:
(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

Raku

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;
}
Output:
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

REXX

For the calculation of the  1,000,000th  Sundaram prime,   it requires a 64-bit version of REXX.

/*REXX program finds & displays  N  Sundaram primes, or displays the Nth Sundaram prime.*/
parse arg n cols .                               /*get optional number of primes to find*/
if    n=='' |    n==","  then    n= 100          /*Not specified?   Then assume default.*/
if cols=='' | cols==","  then cols=  10          /* "      "          "     "       "   */
@.= .;                             lim= 16 * n   /*default value for array; filter limit*/
       do    j=1  for n;   do k=1  for n  until _>lim;  _= j + k + 2*j*k;  @._=
                           end   /*k*/
       end      /*j*/
w= 10                                            /*width of a number in any column.     */
                                     title= 'a list of '   commas(N)   " Sundaram primes"
if cols>0  then say ' index │'center(title,  1 + cols*(w+1)     )
if cols>0  then say '───────┼'center(""   ,  1 + cols*(w+1), '─')
#= 0;                         idx= 1             /*initialize # of Sundaram primes & IDX*/
$=                                               /*a list of Sundaram primes  (so far). */
       do j=1  until #==n                        /*display the output (if cols > 0).    */
       if @.j\==.  then iterate                  /*Is the number not prime?  Then skip. */
       #= # + 1                                  /*bump number of Sundaram primes found.*/
       a= j                                      /*save J for calculating the Nth prime.*/
       if cols<=0  then iterate                  /*Build the list  (to be shown later)? */
       c= commas(j + j + 1)                      /*maybe add commas to  Sundaram  prime.*/
       $= $  right(c, max(w, length(c) ) )       /*add Sundaram prime──►list, allow big#*/
       if #//cols\==0  then iterate              /*have we populated a line of output?  */
       say center(idx, 7)'│'  substr($, 2);  $=  /*display what we have so far  (cols). */
       idx= idx + cols                           /*bump the  index  count for the output*/
       end   /*j*/

if $\==''  then say center(idx, 7)"│"  substr($, 2)  /*possible display residual output.*/
if cols>0  then say '───────┴'center(""   ,  1 + cols*(w+1), '─')
say
say 'found ' commas(#)  " Sundaram primes, and the last Sundaram prime is "  commas(a+a+1)
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 ?
output   when using the default inputs:
 index │                                        a list of  100  Sundaram primes
───────┼───────────────────────────────────────────────────────────────────────────────────────────────────────────────
   1   │          3          5          7         11         13         17         19         23         29         31
  11   │         37         41         43         47         53         59         61         67         71         73
  21   │         79         83         89         97        101        103        107        109        113        127
  31   │        131        137        139        149        151        157        163        167        173        179
  41   │        181        191        193        197        199        211        223        227        229        233
  51   │        239        241        251        257        263        269        271        277        281        283
  61   │        293        307        311        313        317        331        337        347        349        353
  71   │        359        367        373        379        383        389        397        401        409        419
  81   │        421        431        433        439        443        449        457        461        463        467
  91   │        479        487        491        499        503        509        521        523        541        547
───────┴───────────────────────────────────────────────────────────────────────────────────────────────────────────────

found  100  Sundaram primes, and the last Sundaram prime is  547
output   when using the inputs of:     1000000   0
found  1,000,000  Sundaram primes, and the last Sundaram prime is  15,485,867

Ruby

Based on the Python code from the Wikipedia lemma.

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]}"
Output:
[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

Wren

Library: 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.

import "./fmt" for Fmt

var sos = Fn.new { |n|
    if (n < 3) return []
    var primes = []
    var k = ((n-3)/2).floor + 1
    var marked = List.filled(k, true)
    var limit = ((n.sqrt.floor - 3)/2).floor + 1
    limit = limit.max(0)
    for (i in 0...limit) {
        var p = 2*i + 3
        var s = ((p*p - 3)/2).floor
        var j = s
        while (j < k) {
            marked[j] = false
            j = j + p
        }
    }
    for (i in 0...k) {
        if (marked[i]) primes.add(2*i + 3)
    }
    return primes
}

// odds only
var soe = Fn.new { |n|
    if (n < 3) return []
    var primes = []
    var k = ((n-3)/2).floor + 1
    var marked = List.filled(k, true)
    var limit = ((n.sqrt.floor - 3)/2).floor + 1
    limit = limit.max(0)
    for (i in 0...limit) {
        if (marked[i]) {
            var p = 2*i + 3
            var s = ((p*p - 3)/2).floor
            var j = s
            while (j < k) {
                marked[j] = false
                j = j + p
            }
        }
    }
    for (i in 0...k) {
        if (marked[i]) primes.add(2*i + 3)
    }
    return primes
}

var limit = 16e6 // say
var start = System.clock
var primes = sos.call(limit)
var elapsed = ((System.clock - start) * 1000).round
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.tprint("$3d", primes[0..99], 10)
Fmt.print("\nThe $,d Sundaram prime is $,d", 1e6, primes[1e6-1])

start = System.clock
primes = soe.call(limit)
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])
Output:
Using the Sieve of Sundaram generated primes up to 16,000,000 in 1,232 ms.

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,000 Sundaram prime is 15,485,867

Using the Sieve of Eratosthenes would have generated them in 797 ms.

As a check, the 1,000,000 Sundaram prime would again have been 15,485,867