{-#LANGUAGE DataKinds#-}
{-#LANGUAGE TypeOperators#-}

module Problem650 where

import qualified Data.HashMap.Strict as M
import Data.Modular
import qualified Data.Vector as V
import Data.Maybe

lim = 20000 :: Int

(*!) = M.unionWith (+)
(/!) a b = M.filter (> 0) $ M.unionWith (-) a b

factorizations :: V.Vector (M.HashMap Int Int)
-- starts at n=1
factorizations = V.fromList (M.empty : [ getFact i | i <- [2 .. lim] ])
  where
    getFact n =
        let
            ps = filter (\p -> mod n p == 0)
                $ takeWhile (\p -> p * p <= n) (2 : [3, 5 ..])
        in
            if null ps
                then M.singleton n 1
                else
                    let p = head ps
                    in M.insertWith (+) p 1 (factorizations V.! ((n `div` p) - 1))

factorialFactorizations :: V.Vector (M.HashMap Int Int)
factorialFactorizations = V.scanl1 (*!) factorizations

f :: V.Vector (M.HashMap Int Int)
f = V.scanl1 (*!) factorialFactorizations

d :: Int -> Int / 1000000007
-- nC0*nC1..nCn = n!^(n+1) / (0! 1! ... n!)^2
-- B(5) = 5!^6/(0!*1!...5!)^2 = 120^6/(34560)^2 = 2500
d n =
    product
        $  M.mapWithKey (\p v -> let p' = toMod p in (p' ^ (v + 1) - 1) * (fromJust $ inv (p' - 1)))
        $  fmap (* (n + 1)) (factorialFactorizations V.! (n - 1))
        /! fmap (* 2)       (f V.! (n - 1))

main :: IO ()
main = print $ sum $ map d [1 .. lim]


