{-
Search for safe primes below doubly exponential decimal powers of two.

Example usage:

./a.out go 1000000 primediffs.gz 14.2 0

1000000 is sieve size.

primediffs.gz is consecutive differences between the first (say) 100 million primes starting at 3, processed with (x/2)-1, written as bytes, then compressed with gzip.

14.2 is exponent, searching for safe primes less than 2^round(2^14.2)-offset = 2^18820

0 is offset, to continue from interrupted previous computations.

This program is has a space leak.

Copyright 2021 Ken Takusagawa

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or (at
your option) any later version.

This program is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program.  If not, see <http://www.gnu.org/licenses/>.

-}

{-# LANGUAGE ScopedTypeVariables, LambdaCase, PackageImports #-}
{-# LANGUAGE FlexibleContexts #-} -- MArray context with Bool
module Main(main) where {
import System.Environment(getArgs);
import Control.Exception(assert);
import Data.Function((&));
import Control.Category((>>>));
import Prelude hiding((.),(>>),gcd);
--import qualified Prelude;
import System.IO(hSetBuffering,stdout,BufferMode(LineBuffering));
import qualified Data.Array.MArray.Safe as Mutable;
import qualified Data.Array.IO.Safe as Ioarray;
import qualified Primality2;
import Math.NumberTheory.Euclidean(extendedGCD); -- arithmoi
import Data.ByteString.Lazy(ByteString);
import qualified Data.ByteString.Lazy as ByteString;
import qualified Codec.Compression.GZip as Gzip;
import Control.Monad.Extra(unfoldM_);
import Text.Printf(printf);
import System.CPUTime(getCPUTime);
import Data.IORef(IORef,newIORef,readIORef,modifyIORef');

newtype Sievesize = Sievesize Ii deriving(Show);

readsievesize :: String -> Sievesize;
readsievesize = read >>> Sievesize;

main :: IO();
main = getArgs >>= \case{
["go",sievesize,gzprimes,e,offset] -> decimalgz (readsievesize sievesize) gzprimes (read e) (Startoffset$ read offset);
_ -> undefined;
};

type Ii = Integer;

newtype Magicmod = Magicmod Ii deriving Show;
newtype Theprime = Theprime Ii deriving Show;
newtype Center = Center Ii deriving (Show);

{-
we want
center-n == magicmod-1 mod magicmod
center-n == 0 mod theprime
-}

getxy :: Magicmod -> Theprime -> (Ii,Ii);
getxy (Magicmod g) (Theprime p) = let {
(gcd,x,y) = extendedGCD g (negate p)
} in assert (gcd == 1) $ (x,y);

rawfirstn :: Center -> Magicmod -> Theprime -> Ii;
rawfirstn (Center c) g@(Magicmod magicmod) p@(Theprime theprime) = let {
x = getxy g p & fst;
k = div (c+1-magicmod*x) (magicmod*theprime);
} in c - magicmod*(x+k*theprime) + 1;

firstnq :: Center -> Magicmod -> Theprime -> Ii;
firstnq c@(Center center) g@(Magicmod magicmod) p@(Theprime theprime) =
  assert (mod center magicmod == magicmod - 1)
$ assert (r == 0)
$ assert (0 <= q)
$ assert (mod (center - answer) theprime == 0)
$ q where {
answer = rawfirstn c g p;
(q,r) = divMod answer magicmod;
};

newtype Stride = Stride Ii deriving (Show);

linemark :: (Mutable.MArray a Bool m) => a Ii Bool -> Center -> Magicmod -> Theprime -> m ();
linemark arr c g p@(Theprime theprime) = markmany arr (Stride theprime) $ firstnq c g p;

markmany :: (Mutable.MArray a Bool m) => a Ii Bool -> Stride -> Ii -> m();
markmany arr s@(Stride stride) x = do {
maxb <- Mutable.getBounds arr >>= (snd >>> return);
if x <= maxb
then do {
Mutable.writeArray arr x False;
markmany arr s $ x+stride
} else return ();
};

safecentercalc :: Startcenter -> (Center, Magicmod);
safecentercalc = centercalc $ Magicmod 12;

centercalc :: Magicmod -> Startcenter -> (Center, Magicmod);
centercalc (Magicmod g) (Startcenter start) = assert (mod answer g == g-1) $ (Center answer, Magicmod g) where {
answer = start + (g-1 - mod start g);
};

sophiecentercalc :: Startcenter -> (Center, Magicmod);
sophiecentercalc = safecentercalc >>> fst >>> \case {
Center safe -> case divMod safe 2 of {
(q,1) -> assert (mod q 6 == 5) (Center q, Magicmod 6);
_ -> error "center of safe primes should be odd";
}};

newtype Startcenter = Startcenter Ii deriving (Show);

-- cannot derive show on IOUArray
data Arrstuff = Arrstuff (Ioarray.IOUArray Ii Bool) (Center,Magicmod);

makeit :: Sievesize -> (Center,Magicmod) -> IO Arrstuff;
makeit (Sievesize sizemax) cg = do {
arr <- Mutable.newArray (0,sizemax) True;
return $ Arrstuff arr cg;
};

mksafe :: Sievesize -> Startcenter -> IO Arrstuff;
mksafe sievesize = safecentercalc >>> makeit sievesize;

mksophie :: Sievesize -> Startcenter -> IO Arrstuff;
mksophie sievesize = sophiecentercalc >>> makeit sievesize;

doaprime :: Theprime -> Arrstuff -> IO();
doaprime p (Arrstuff arr (c,g)) = linemark arr c g p;

bothmark :: [Arrstuff] -> Theprime -> IO();
bothmark arrs p = mapM_ (doaprime p) arrs;

readone :: Ii -> Arrstuff -> IO Bool;
readone i (Arrstuff arr _) = Mutable.readArray arr i;

-- this assumes it is a safe prime, so both p and (p-1)/2 are prime
isgenerator2 :: Ii -> Bool;
isgenerator2 p = let {
generatortotest = 2
} in Primality2.modular_exponentiation p generatortotest (div (p-1) 2) /= 1;

uncompressprimes :: [Ii] -> [Ii];
-- space leak
--uncompressprimes l = 2:scanl (\old diff -> let {answer = old + 2*(diff+1)} in seq answer answer) 3 l;
uncompressprimes l = 2: uncompressprimes1 3 l where {
uncompressprimes1 accum [] = seq accum $ accum : [];
uncompressprimes1 accum (diff:t) = seq accum $ accum : uncompressprimes1 (accum + 2*(diff+1)) t;
};

frombytestring :: ByteString -> [Ii];
frombytestring = ByteString.unpack >>> map fromIntegral;

blocktest2 :: (Primeblob primeblob) => Sievesize -> Distantcenter -> Startcenter -> primeblob -> IO Bool;
blocktest2 sievesize@(Sievesize sizemax) distant@(Distantcenter distant1) center@(Startcenter start1) compressedprimes = do {
putStrLn $ "blocktest2 distant-start = " ++ (show $ distant1 - start1);
cpu "start";
let {myprimes :: [Theprime] = getprimes compressedprimes & drop 2 };
arr :: [Arrstuff] <- sequence [mksafe sievesize center,mksophie sievesize center];
mapM_ (bothmark arr) myprimes;
cpu "primemark";

workcounter :: IORef Ii <- newIORef 0;
answer <- icountmbool 0 sizemax (onetest distant arr workcounter);
readIORef workcounter >>= (\x -> putStrLn $ "workcounter " ++ show x);
return answer;

-- doing just printing still has a space leak
--Monad.forM_ [0..sizemax-1] $ oneprint distant arr; return True;
};

icountmbool :: (Monad m) => Ii -> Ii -> (Ii -> m Bool) -> m Bool;
icountmbool i lastvalue f = f i >>= \case {
False -> if i == lastvalue then return False
else icountmbool (1+i) lastvalue f;
True -> return True;
};

-- not necessarily this block of primes but could be earlier
newtype Distantcenter = Distantcenter Ii deriving (Show);

onetest :: Distantcenter -> [Arrstuff] -> IORef Ii -> Ii -> IO Bool;
onetest (Distantcenter start) arrs workcounter i = -- hPutStrLn stderr ("onetest " ++ show i ++ " " ++ (arrs & map arrstuffgetshowable & show)) Prelude.>>
mapM (readone i) arrs >>= (and >>> \case{
False -> return False;
True -> do {
let {
  center = case head arrs of {
    Arrstuff _ (Center c, Magicmod 12) -> c;
    _ -> error "magicmod not 12";
  };
  n = center - 12*i;
  nm = let {m = div (n-1) 2;}in [m,n];
  mr = map Primality2.miller_rabin_2 nm;
  sq = map Primality2.is_almost_square nm;
  bp = map Primality2.strong_lucas nm;
};
if and mr && (not (or sq)) && and bp
then do {putStrLn $ "found "++show (start-n);
  return $ isgenerator2 n;
} else do {
  --putStrLn $ "sieve ok but not work " ++ show (n-start);
  modifyIORef' workcounter (+1);
  readIORef workcounter >>= (\x -> if mod x 100 /= 0 then return ()
    else do{ cputime <- getCPUTime ; putStrLn $ "snapshot work " ++ show x ++ " offset " ++ show (start-n) ++ " cpu " ++ show cputime;
         -- cpu "asnapshot"
  });
  return False
}}});

class Primeblob a where {
getprimes :: a -> [Theprime];
};

-- optimizer could try to save this, in which case terrible.
instance Primeblob Diffprimes where {
getprimes (Diffprimes compressedprimes) = compressedprimes & frombytestring & uncompressprimes & map Theprime
};

newtype Diffprimes = Diffprimes ByteString;

newtype Gzprimes = Gzprimes {ungzprimes :: ByteString};

instance Primeblob Gzprimes where {
getprimes (Gzprimes compressedprimes) = compressedprimes & Gzip.decompress & frombytestring & uncompressprimes & map Theprime
};

newtype Startoffset = Startoffset Ii deriving (Show);

blocktestmany :: (Primeblob primeblob) => Sievesize -> primeblob -> Distantcenter -> Startoffset -> IO();
blocktestmany sievesize compressedprimes distant@(Distantcenter start) (Startoffset offset) = do{
unfoldM_ (blocktest1 sievesize distant compressedprimes) (Startcenter $ start-offset);};

blocktest1 :: (Primeblob primeblob) => Sievesize -> Distantcenter -> primeblob -> Startcenter -> IO (Maybe Startcenter);
blocktest1 sievesize@(Sievesize sizemax) distant compressedprimes start@(Startcenter start1) =
--putStr "b" Prelude.>>
 blocktest2 sievesize distant start compressedprimes >>= (\case {
True -> Nothing;
False -> start1 - 12*(sizemax +1) & Startcenter &Just} >>> return);

decimalgz :: Sievesize -> String -> Double -> Startoffset -> IO ();
decimalgz sievesize primefilegz e offset = do {
hSetBuffering stdout LineBuffering;
putStrLn "$Id: safeprimesieve.hs,v 1.28 2021/05/12 06:52:58 kenta Exp kenta $";
gzprimes :: Gzprimes <- ByteString.readFile primefilegz >>= (Gzprimes >>> return);
putStrLn $ "gzmemory " ++ (show $ByteString.length $ ungzprimes gzprimes);
print sievesize;
printf "e %.1f\n" e;
blocktestmany sievesize gzprimes (Distantcenter $ double2pow e) offset;
};

cpu :: String -> IO ();
cpu tag = do {
putStr $ "cpu "++tag++" ";
getCPUTime >>= print;
};

pow2 :: Double -> Ii;
pow2 e = 2**e & round;

double2pow :: Double -> Ii;
double2pow e = 2^pow2 e;

} --end