{-
  Copyright 2012 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,GeneralizedNewtypeDeriving,RankNTypes #-}
module Main (main) where{
import System.Environment(getArgs);
import Data.List(groupBy,sortBy);
import Primality;
import IO;
import Data.Array.IArray;
import Data.Array.Unboxed;
import Control.Exception(assert);
import Control.Monad(when);
import Debug.Trace(trace);

main :: IO(());
main = (getArgs >>= (\lambda_case_var ->case lambda_case_var of {
  ["nothing"]-> (return ());
  ["sortPs"]-> ((mapM_ print)(sortPs));
  ["numbers"]-> ((mapM_ print)((map number)(sortPs)));
  ["large", n]-> ((mapM_ print)(getPs(read(n))));
  ["era1"]-> ((mapM_ print)((groupBy same_era)(sortPs)));
  ["era-limit", n, m]-> ((mapM_ single_print2)((take (read m))((drop (read n))((groupBy same_era)(sortPs)))));
  ["era", n]-> (do{
  linebuffering;
  ((mapM_ single_print2)((drop (read n))((groupBy same_era)(sortPs))));
});
  ["era10", n, part]-> (do{
  linebuffering;
  ((mapM_ single_print2)((filter (cpu_divide (read part)))((drop (read n))((groupBy same_era)(sortPs)))));
});
  ["div-io"]-> ((mapM_ div_io)(sortPs));
  ["test-lucas", d, length]-> ((mapM_ test_lucas_pierpont)((take (read length))((drop (read d))(sortPs))));
  ["time", d, amount]-> (print(length((filter time_pierpont)((take (read amount))((drop (read d))(sortPs))))));
  ["nth", n]-> (print(getlog(((!!) sortPs (read n)))))
}));
data P  = P(Integer)(Integer) deriving (Show, Eq);
getlog :: P -> Double;
getlog (P (two :: Integer ) (three :: Integer )) = (((+) (fromInteger two))(((*) log2_3)(fromInteger(three))));
log2_3 :: Double;
log2_3 = ((/) (log 3) (log 2));
ordP :: P -> P -> Ordering;
ordP x y = (compare (getlog x) (getlog y));
getPs1 :: Integer -> [](P);
getPs1 max = (do{
  x :: Integer <- (enumFromTo 1 max);
  y :: Integer <- (enumFromTo 0 (floor ((/) (fromInteger max) log2_3)));
  (return (P x y));
});
getPs :: Integer -> [](P);
getPs max = ((takeWhile ((/=) (P max 0)))((sortBy ordP)(getPs1(max))));
five_isPrimeP :: P -> Bool;
five_isPrimeP x = (miller_rabin_1 (((+) 1)(number(x))) 5);
isPrimeP :: P -> Bool;
isPrimeP x = (((&&) (all_check2 x))(isPrimeP_miller_rabin(x)));
isPrimeP_miller_rabin :: P -> Bool;
isPrimeP_miller_rabin x = (miller_rabin_isPrime(((+) 1)(number(x))));
linebuffering :: IO(());
linebuffering = (hSetBuffering stdout LineBuffering);
run1 :: Integer -> IO(());
run1 n = (do{
  linebuffering;
  (mapM_ (single_print (getPs n)) (enumFromTo 2 n));
});
find_lt :: [](P) -> P -> P;
find_lt ps target = (head((dropWhile (not . isPrimeP))(reverse((takeWhile ((/=) target))(ps)))));
single_print :: [](P) -> Integer -> IO(());
single_print ps i = (do{
  (putStr ((show i) ++ ": "));
  (print((find_lt ps (P i 0))));
});
inc2 :: P -> P;
inc2 (P (two :: Integer ) (three :: Integer )) = (P ((+) 1 two) three);
inc3 :: P -> P;
inc3 (P (two :: Integer ) (three :: Integer )) = (P two ((+) 1 three));
sortPs :: [](P);
sortPs = ((:) (P 1 0) (merge (map inc2 sortPs) (map inc3 sortPs)));
merge :: [](P) -> [](P) -> [](P);
merge x y = (case (ordP (head x) (head y)) of {
  LT-> ((:) (head x) (merge (tail x) y));
  GT-> ((:) (head y) (merge x (tail y)));
  EQ-> ((:) (head x) (merge (tail x) (tail y)))
});
same_era :: P -> P -> Bool;
same_era x y = (let {
era :: P -> Int;
era i = (floor(getlog(i)))
}
 in ((==) (era x) (era y)));
maybe_head :: [](a) -> Maybe(a);
maybe_head l = (case l of {
  ([] )-> Nothing;
  _ -> (Just(head(l)))
});
single_print2 :: [](P) -> IO(());
single_print2 g = (print ((head(g)), (maybe_head((dropWhile (not . isPrimeP))(reverse(g))))));
number :: P -> Integer;
number (P (two :: Integer ) (three :: Integer )) = (gen_number two three);
gen_number :: Integer -> Integer -> Integer;
gen_number two three = ((*) ((^) 2 two) ((^) 3 three));
type Aii  = Array((Int, Int))(Bool);
div_n_array :: Integer -> Aii;
div_n_array n = (let {
minus :: Int;
minus = (pred(pred(fromInteger(n))))
}
 in (array ((0, 0), (minus, minus)) (do{
  a :: Int <- (enumFromTo 0 minus);
  b :: Int <- (enumFromTo 0 minus);
  (return ((a, b), (((/=) 0)(((flip mod) n)(((+) 1)((gen_number (toInteger a) (toInteger b))))))));
})));
type Div_array  = (Aii, Int);
make_div_array :: Integer -> Div_array;
make_div_array n = ((div_n_array n), (fromInteger n));
all_div :: [](Div_array);
all_div = ((map make_div_array)((take 40)((drop 2)(primes(())))));
div_check :: P -> Div_array -> Bool;
div_check (P (two :: Integer ) (three :: Integer )) ((darray :: Aii ), (n :: Int )) = (let {
m :: Integer -> Int;
m x = (mod (fromInteger x) (pred n))
}
 in ((!) darray ((m two), (m three))));
div5 :: P -> Bool;
div5 p = ((div_check p)(head(all_div)));
all_check :: P -> Bool;
all_check p = (and((map (div_check p))(all_div)));
div_io :: P -> IO(());
div_io n = (case ((==) ((/=) 0 (mod ((+) 1 (number n)) 5)) (div5 n)) of {
  False-> (putStrLn ("differs " ++ (show n)));
  _ -> (return ())
});
lucas_primality_test_1factor :: Integer -> Integer -> Integer -> Bool;
lucas_primality_test_1factor test base factor = (case (divMod (pred test) factor) of {
  ((q :: Integer ), (remainder :: Integer ))-> ((assert ((==) 0 remainder)) (not_mod_1 base q test))
});
not_mod_1 :: Integer -> Integer -> Integer -> Bool;
not_mod_1 base exponent modulus = ((/=) 1 (modular_exponentiation base exponent modulus));
data Lucas_result  = Prime_lucas_result | Composite_lucas_result | Unknown_lucas_result deriving (Show, Eq);
lucas_primality_test1 :: Integer -> [](Integer) -> Integer -> Lucas_result;
lucas_primality_test1 test factors base = (case (miller_rabin_1 test base) of {
  False-> Composite_lucas_result;
  _ -> (case (and (map (lucas_primality_test_1factor test base) factors)) of {
  True-> Prime_lucas_result;
  False-> Unknown_lucas_result
})
});
lucas_primality_with_witness :: Integer -> [](Integer) -> (Integer, Bool);
lucas_primality_with_witness test factors = (let {
loop :: Integer -> (Integer, Bool);
loop base = (case (lucas_primality_test1 test factors base) of {
  Prime_lucas_result-> (base, True);
  Composite_lucas_result-> (base, False);
  Unknown_lucas_result-> (loop ((+) 1 base))
})
}
 in (loop 2));
lucas_pierpont :: P -> (Integer, Bool);
lucas_pierpont p = (lucas_primality_with_witness ((+) 1 (number p)) (pierpont_get_factors p));
optimized_lucas_pierpont :: P -> (Integer, Bool);
optimized_lucas_pierpont p = (optimized_lucas_primality_with_witness ((+) 1 (number p)) (pierpont_get_factors p));
pierpont_get_factors :: P -> [](Integer);
pierpont_get_factors (P (two :: Integer ) (three :: Integer )) = ((++) (case two of {
  0-> [];
  _ -> [2]
}) (case three of {
  0-> [];
  _ -> [3]
}));
test_lucas_pierpont :: P -> IO(());
test_lucas_pierpont p = (when ((/=) (isPrimeP_miller_rabin p) (snd(optimized_lucas_pierpont(p)))) (putStrLn ("fail " ++ (show p))));
opt_lucas :: Integer -> [](Integer) -> Integer -> Lucas_result;
opt_lucas test factors base = (let {
all_factors :: Integer;
all_factors = (product factors);
remaining :: [](Integer);
remaining = (do{
  f :: Integer <- factors;
  (return (div all_factors f));
})
}
 in (case (divMod (pred test) all_factors) of {
  ((pre_multiplied :: Integer ), (r :: Integer ))-> ((assert ((==) 0 r)) ((case (miller_rabin_1 test base) of {
  False-> Composite_lucas_result;
  _ -> (let {
almost_there :: Integer;
almost_there = (modular_exponentiation base pre_multiplied test)
}
 in (case (and (do{
  leftover :: Integer <- remaining;
  (return (not_mod_1 almost_there leftover test));
})) of {
  True-> Prime_lucas_result;
  False-> Unknown_lucas_result
}))
}) ))
}));
base_sequence :: [](Integer);
base_sequence = (drop 2 (primes ()));
optimized_lucas_primality_with_witness :: Integer -> [](Integer) -> (Integer, Bool);
optimized_lucas_primality_with_witness test factors = (let {
loop :: [](Integer) -> (Integer, Bool);
loop ((:) (base :: Integer ) (more_bases :: [](Integer) )) = (case (opt_lucas test factors base) of {
  Prime_lucas_result-> (base, True);
  Composite_lucas_result-> (base, False);
  Unknown_lucas_result-> (loop more_bases)
})
}
 in (loop base_sequence));
time_pierpont :: P -> Bool;
time_pierpont p = (snd(lucas_pierpont(p)));
powers_mod :: Int -> Int -> [](Int);
powers_mod base modulus = (let {
f :: Int -> Int;
f x = (mod ((*) x base) modulus)
}
 in ((:) 1 (map f (powers_mod base modulus))));
type Adiv  = UArray(Int)(Int);
single_divisibility_array :: Int -> Int -> Adiv;
single_divisibility_array base modulus = (let {
minus :: Int;
minus = (pred(pred(modulus)))
}
 in (listArray (0, minus) (powers_mod base modulus)));
pair_divisibility_array :: Int -> (Adiv, Adiv);
pair_divisibility_array modulus = ((single_divisibility_array 2 modulus), (single_divisibility_array 3 modulus));
all_pair_divisibilities :: []((Adiv, Adiv));
all_pair_divisibilities = ((map pair_divisibility_array)((map fromInteger)((take 300)((drop 2)(primes(()))))));
div_check_pair :: P -> (Adiv, Adiv) -> Bool;
div_check_pair (P (two :: Integer ) (three :: Integer )) ((atwo :: Adiv ), (athree :: Adiv )) = (let {
modulus :: Int;
modulus = (((+) 2)(snd(bounds(atwo))));
m :: Integer -> Int;
m i = (mod (fromInteger i) (pred modulus))
}
 in ((/=) 0 (mod ((+) 1 ((*) ((!) atwo (m two)) ((!) athree (m three)))) modulus)));
all_check2 :: P -> Bool;
all_check2 p = (and((map (div_check_pair p))(all_pair_divisibilities)));
cpu_divide :: Integer -> [](P) -> Bool;
cpu_divide part g = (let {
answer :: Bool;
answer = ((==) part (mod (case (head g) of {
  (P (two :: Integer ) 0)-> two
}) 10))
}
 in (answer))
}