{-
Draw bar graphs of the how many times a prime p divides a central
binomial coefficient.  Output is NetPBM.

Copyright 2018 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 #-}
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 System.IO(stderr,hPutStrLn);
import qualified Data.List as List;
--import qualified Control.Monad as Monad;
--import qualified Data.Maybe as Maybe;
--import qualified Data.Map as Map; import Data.Map(Map);
--import qualified Data.Set as Set; import Data.Set(Set);
--import qualified Data.Bifunctor as Bifunctor;

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

main :: IO();
main = getArgs >>= \case{
["go",n,p] -> calcrange (read n) (read p) & genpbm & putStr;
_ -> undefined;
};

factorial :: Integer -> Integer;
factorial n = product [1..n];

central :: Integer -> Integer;
central n = div_exact (factorial (2*n) ) (sqr (factorial n));

sqr :: Integer -> Integer;
sqr x = x*x;

div_exact :: Integer -> Integer -> Integer;
div_exact x y = case divMod x y of{
(z,0) -> z;
_ -> error "did not divide cleanly";
};

findpow :: Integer -> Integer -> Integer;
findpow n p = case divMod n p of {
(x,0) -> 1+findpow x p;
_ -> 0;
};

allrange :: Int -> [Integer];
allrange n = enumFrom 0 & take n;

calcrange :: Int -> Integer -> [Integer];
calcrange n p = allrange n & map (central >>> flip findpow p);

gen_row :: Integer -> Integer -> String;
gen_row m x = List.genericReplicate x "1" ++ repeat "0" & List.genericTake m & unwords;

genpbm :: [Integer] -> String;
genpbm l = unlines (["P1",maximum l & show, length l & show] ++ map (gen_row (maximum l)) l);

} --end