Definitions

   A   safe prime   is a prime   p   and where   (p-1)/2   is also prime.
   The corresponding prime  (p-1)/2   is known as a   Sophie Germain   prime.
   An   unsafe prime   is a prime   p   and where   (p-1)/2   isn't   a prime.
   An   unsafe prime   is a prime that   isn't   a   safe   prime.


Task

   Find and display (on one line) the first   35   safe primes.
   Find and display the   count   of the safe primes below   1,000,000.
   Find and display the   count   of the safe primes below 10,000,000.
   Find and display (on one line) the first   40   unsafe primes.
   Find and display the   count   of the unsafe primes below   1,000,000.
   Find and display the   count   of the unsafe primes below 10,000,000.
   (Optional)   display the   counts   and   "below numbers"   with commas.
Show all output here.



Related Task

   strong and weak primes.


Also see

   The OEIS article:     safe   primes.
   The OEIS article:   unsafe primes.


isqrt = proc (s: int) returns (int)
    x0: int := s/2
    if x0=0 then return(s) end
    x1: int := (x0 + s/x0)/2
    while x1 < x0 do
        x0 := x1
        x1 := (x0 + s/x0)/2
    end
    return(x0)
end isqrt
 
sieve = proc (n: int) returns (array[bool])
    prime: array[bool] := array[bool]$fill(0,n+1,true)
    prime[0] := false
    prime[1] := false
    for p: int in int$from_to(2, isqrt(n)) do
        if prime[p] then
            for c: int in int$from_to_by(p*p,n,p) do
                prime[c] := false
            end
        end
    end
    return(prime)
end sieve
 
start_up = proc ()
    primeinfo = record [
        name: string,
        ps: array[int],
        maxps, n_1e6, n_1e7: int
    ]
 
    po: stream := stream$primary_output()
    prime: array[bool] := sieve(10000000)
 
    safe: primeinfo := primeinfo${
        name: "safe",
        ps: array[int]$[],
        maxps: 35,
        n_1e6: 0,
        n_1e7: 0
    }
 
    unsafe: primeinfo := primeinfo${
        name: "unsafe",
        ps: array[int]$[],
        maxps: 40,
        n_1e6: 0,
        n_1e7: 0
    }
 
    for p: int in int$from_to(2, 10000000) do
        if ~prime[p] then continue end
        ir: primeinfo 
        if prime[(p-1)/2] 
            then ir := safe 
            else ir := unsafe
        end
 
        if array[int]$size(ir.ps) < ir.maxps then
            array[int]$addh(ir.ps,p)
        end
        if p<1000000 then ir.n_1e6 := ir.n_1e6 + 1 end
        if p<10000000 then ir.n_1e7 := ir.n_1e7 + 1 end
    end
 
    for ir: primeinfo in array[primeinfo]$elements(
                       array[primeinfo]$[safe, unsafe]) do
        stream$putl(po, "First " || int$unparse(ir.maxps)
                   || " " || ir.name || " primes:")
        for i: int in array[int]$elements(ir.ps) do
            stream$puts(po, int$unparse(i) || " ")
        end
        stream$putl(po, "\nThere are " || int$unparse(ir.n_1e6)
                      || " " || ir.name || " primes < 1,000,000.")
        stream$putl(po, "There are " || int$unparse(ir.n_1e7)
                      || " " || ir.name || " primes < 1,000,000.\n")
    end
end start_up