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.


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