1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
| module Main where
import Control.Monad
import Control.Parallel.Strategies
import Data.Array.ST
import Data.Array.Unboxed
import Data.Bits
import Math.NumberTheory.Primes
import System.Environment
import System.IO
import Text.Printf
safePrime :: Int -> Integer
safePrime bits | bits < 100 =
head $ do
let p0 = 2^bits 1
p <- [p0, p0 2 ..]
guard (isPrime p && isPrime (div p 2))
return p
safePrime bits =
head $ do
let stepSize = 24 * fromIntegral bits
step <- [0, stepSize ..]
let p0 = 2^bits step 1
p1 = 2^bits step stepSize 1
let s = sieve p1 p0
p <- [p0, p0 2 .. p1]
guard (s ! p && bailliePSW p && bailliePSW (div p 2))
return p
where
sieve a b =
runSTUArray $ do
s <- newArray (a, b) True
let primes' = take 1000 primes
mapM_ (\i -> writeArray s i False) $ do
p <- primes'
let i0 = a + mod (a) p
a' = div a 2
j0 = a' + mod (a') p
[i0, i0 + p .. b] ++ [2*j0 + 1, 2*(j0 + p) + 1 .. b]
return s
main :: IO ()
main = do
hSetBuffering stdout LineBuffering
bs <- fmap (map read) getArgs
let ps = parMap rdeepseq (\b -> (b, safePrime b)) bs
forM_ ps $ \(b, p) ->
printf "2^%d - %d\n" b (2^b p)
|