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
- The article on Wikipedia.
ALGOL 68
To run this with Algol 68G, you will need to increase the heap size by specifying e.g. -heap=64M
on the command line.
<lang algol68>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</lang>
- 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
AppleScript
The "nth prime" calculation here's gleaned from the Python and Julia solutions and the limitations to marking partly from the Phix.
<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 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()</lang>
- Output:
<lang applescript>"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"</lang>
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.
<lang csharp>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; }
}</lang>
- 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.
F#
<lang fsharp> // 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())) </lang>
- 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
Go
<lang go>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)
}</lang>
- 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
<lang 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</lang>
- 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
JavaScript
<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();
})();</lang>
- 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
Julia
<lang julia> """ 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))
</lang>
- 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)
Nim
<lang 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()</lang>
- 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
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
<lang python>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)
</lang>
- 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. <lang python>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()</lang>
- 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
Raku
<lang perl6>my $nth = 1_000_000;
my $k = Int.new: 2.4 * $nth * log($nth) / 2;
my int @sieve;
@sieve[$k] = 0;
hyper for 1 .. $k -> \i {
my int $j = i; while (my int $l = i + $j + 2 * i * $j) < $k { @sieve[$l] = 1; $j = $j + 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;
}</lang>
- 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. <lang 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 ?</lang>
- 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
Wren
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. <lang ecmascript>import "/fmt" for Fmt import "/seq" for Lst
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:") for (chunk in Lst.chunks(primes[0..99], 10)) Fmt.print("$3d", chunk) 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])</lang>
- 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