Efficient parallel safe prime search

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)