{-

Outputs a sorted list of prime numbers of the form a^2^n + b^2^n.

Copyright 2020 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 #-}
module Main where {
import System.Environment(getArgs);
import Control.Exception(assert);
import Debug.Trace(trace);
import Control.Category((>>>));

import Data.Function((&));
import Prelude hiding((.),(>>));
import qualified Data.List as List;
import qualified Data.Bifunctor as Bifunctor;
import qualified Data.List.Ordered as Ordered;
import Data.Ord(comparing);
import Primality2(isPrime_x);
import System.IO(hSetBuffering,stdout,BufferMode(LineBuffering));
import qualified Control.Monad as Monad;

-- to avoid the redundancy warning
trace_placeholder :: ();
trace_placeholder = (trace,assert) & (id >>> error) "trace_placeholder";

type Ii = Integer;

main :: IO();
main = do {
hSetBuffering stdout LineBuffering;
getArgs >>= \case{
["p"] -> tenlines 21 isprime;
["2p"] -> tenlines 21 twiceprime;
["expo",n] -> let { e :: Exponent;
e = read n & Exponent
} in
-- infinite output
mapM_ {- (print1 e) -} print  (nums e & filter (isprime e));
--  mapM_ (print1 e) (nums e );

_ -> undefined;
};
};

data N = N [Ii] deriving (Show,Eq);

newtype Exponent = Exponent Ii deriving (Show);

value :: Exponent -> N -> Ii;
value e (N x) = map (twopow e) x & sum;

twopow :: Exponent -> Ii -> Ii;
twopow (Exponent e) base = base^((2::Ii)^e);

numterms :: Int;
numterms = 2;

nums :: Exponent -> [N];
nums e = N (replicate numterms 0) & alldescendents (comparing (value e));
{-
00
01
02 11
03 12 -- is 12 a child of 02 or of 11
04 13 22
05 14 23
06 15 24 33

000
001
002 011
(003^,012)(012,111^)
003 012 111
(004^,013)(013,022^,112)(112)
004 013 022 112
(005^ 014)(014 023 113)(023 122)(113 122)
005 014 023 113 122
(006^ 015)(015 024 114)(024 033^ 123)(114 123)(123 222^)
006 015 024 033 114 123 222
(007^ 016)(016 025 115)(025 034 124)(034 133)(115 124)(124 133 223) (223)
007 016 025 034 115 124 133 223
(008^ 017)(017 026 116)(026 035 125)(035 044^ 134)(116 125)(125 134 224)(134 233)(224 233)
008 017 026 035 044 116 125 134 224 233

general plan: generate all children; for each child, generate all parents, am i the lowest numbered possible parent?

-}


print1 :: Exponent -> N -> IO ();
print1 e n = putStrLn $ (show n) ++ " " ++ (show $ value e n);

canonical :: N -> Bool;
canonical (N x) = Ordered.isSorted x;

children :: N -> [N];
children (N x) = do {
(p,q) <- allsplits x ;
case q of {
(h:t) -> return $ N $ p ++ (succ h:t);
_ -> error "no tail";
}};

parents :: N -> [N];
parents (N x) = do {
(p,q) <- allsplits x ;
case q of {
(h:t) -> do {
Monad.guard $ 0 < h;
return $ N $ p ++ (pred h:t);
};
_ -> error "no tail";
}};


allsplits :: [a] -> [([a],[a])];
allsplits [] = [];
allsplits l@(x:rest) = ([],l):(allsplits rest & map (Bifunctor.first (\z -> x:z)));

cchildren :: N -> [N];
cchildren x = do {
assert (canonical x) $ return ();
y ::N <- children x;
Monad.guard $ canonical y;
Monad.guard $ x == (parents y & nminimum);
return y;
};

nminimum :: [N] -> N;
nminimum = map (\(N x) -> x) >>> minimum >>> N;

nsort :: [N] -> [N];
nsort = map (\(N x) -> x) >>> List.sort >>> map N;

alldescendents :: (N -> N -> Ordering) -> N -> [N];
alldescendents cmp x = x : (cchildren x & map (alldescendents cmp) & foldr (Ordered.mergeBy cmp) []);

-- about right for search for primes of form x^1024
myisprime :: Ii -> Bool;
myisprime = isPrime_x 100000;

isprime :: Exponent -> N -> Bool;
isprime e n = value e n & myisprime;
{-
10 10 100 no opt 197 sec
10 10 100 w opt 188
10 10 100 w small 187
10 10 200 w opt 187
10 10 1000 187
10 10 10000 152
10 10 30000 150
10 10 100000 161
-}
twiceprime :: Exponent -> N -> Bool;
twiceprime e n = case divMod (value e n) 2 of
{ (q,0) -> myisprime q
; _ -> False
};

oneline :: Int -> (Exponent -> N -> Bool) -> Exponent -> IO ();
oneline numitems ptest e = do {
putStr $ show e ++ " :";
nums e & filter (ptest e) & take numitems & mapM_ (\(N x) -> putStr $ " " ++ show (reverse x));
putStrLn "";
};

tenlines :: Int -> (Exponent -> N -> Bool) -> IO ();
tenlines numitems ptest = [0..11] & mapM_ (Exponent >>> oneline numitems ptest);


} --end