{-
Compute Goodstein sequences and evaluate the Goodstein function.

Copyright 2016 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, GeneralizedNewtypeDeriving #-}
-- GeneralizedNewtypeDeriving needed to derive Enum
module Main where {
import System.Environment(getArgs);
--import Control.Exception(assert);
--import Debug.Trace(trace);
import Data.Function((&));
import Control.Category((>>>));
import Prelude hiding((.),(>>));
--import System.IO(hSetBuffering,stdout,BufferMode(LineBuffering));
import Data.List;
import Control.Monad;
import Data.Maybe;
--import qualified Data.Map as Map;
--import Data.Map(Map);
--import qualified Data.Set as Set;
--import Data.Set(Set);
import Data.Tuple(swap);

main :: IO();
main = getArgs >>= \case{
["demo19"] -> goodstein_hb_sequence 19 & take 6 & mapM_ print;
["test", b, i] -> test_range (read b) (read i) & print; -- test 100 100000
[n] -> read n & goodstein & print;  -- 4 and larger will take "forever";
-- prints the expressions in hereditary base notation
[n,"expr"] -> read n & goodstein_hb_sequence & mapM_ print;
_ -> error "need argument";
};

newtype Base = Base Integer deriving (Show, Enum);  -- Enum needed for succ

-- output is little-endian
int_to_hb :: Base -> Integer -> HerditaryBase;
int_to_hb base x = zipWith Term (radix_convert base x) (map (int_to_hb base) [0..]);

simplify_term :: Term -> Maybe Term;
simplify_term (Term 0 _) = Nothing;
simplify_term (Term k x) = Just $ Term k $ simplify x;

simplify :: HerditaryBase -> HerditaryBase;
simplify = map simplify_term >>> catMaybes;

-- must be zero or positive.  output is little-endian.
radix_convert :: Base -> Integer -> [Integer];
radix_convert (Base base) = unfoldr $ \n -> if n==0 then Nothing else Just $ swap $ divMod n base;

-- curiously, these Terms may be reordered, i.e., commutative, which is different from normal base N notation.
type HerditaryBase = [Term];

hb_to_int :: Base -> HerditaryBase -> Integer;
hb_to_int base = map (term_to_int base) >>> sum;

-- value is Integer * base ^ HerditaryBase.  base is given exogenously.
data Term = Term Integer HerditaryBase deriving Show;

term_to_int :: Base -> Term -> Integer;
term_to_int _ (Term 0 _) = 0; -- optimization
term_to_int (Base base) (Term k x) = k * base ^ (hb_to_int (Base base) x);

goodstein_tuple :: (Base, Integer) -> (Base, Integer);
goodstein_tuple (base, x) = (succ base
, int_to_hb base x & hb_to_int (succ base) & subtract 1);

goodstein_sequence :: Integer -> [Integer];
goodstein_sequence x = (Base 2,x) & goodstein_sequence_from;

goodstein_sequence_from :: (Base, Integer) -> [Integer];
goodstein_sequence_from = iterate goodstein_tuple >>> map snd >>> takeWhile (>=0);

goodstein :: Integer -> Integer;
goodstein = goodstein_sequence >>> genericLength;

goodstein_hb :: MonadPlus m => m (Base, HerditaryBase) -> m (Base, HerditaryBase);
goodstein_hb input = do {
 (base, x) <- input;
 let { i = hb_to_int (succ base) x & subtract 1 };
 guard $ i>0;
 return (succ base, i & int_to_hb (succ base) & simplify);
};

goodstein_hb_sequence :: Integer -> [(Base, HerditaryBase)];
goodstein_hb_sequence start = (Base 2, int_to_hb (Base 2) start) & goodstein_hb_sequence_from;

goodstein_hb_sequence_from :: (Base, HerditaryBase) -> [(Base, HerditaryBase)];
goodstein_hb_sequence_from = Just >>> iterate goodstein_hb >>> takeWhile isJust >>> map fromJust;

hb_size :: HerditaryBase -> Integer;
hb_size = map (\(Term _ x) -> 1 + hb_size x) >>> sum;

test :: Base -> Integer -> Bool;
test base i = i == (int_to_hb base i & simplify & hb_to_int base);

test_range :: Integer -> Integer -> Bool;
test_range maxbase max_i = and $ do {
base <- [2..maxbase];
i <- [0..max_i];
return $ test (Base base) i;
};

} --end