import Data.List
import Data.Ratio

{--
 -- I would have liked to implement a faster algorithm like
 -- https://www.williamstein.org/edu/2006/spring/583/notes/2006-04-14/2006-04-14.pdf
 -- but could not get the fixed-point arithmetic to work in a timely manner 
 -- (using only the base package)...
 -- The below implementation is terribly slow in coparison, 
 -- since it is based on the double sum obtained by Worpotzky 
 -- (using the StirlingNumbers of the 2nd kind).
 -- The binominal coefficient runs reasonably fast and the haskell standard 
 -- for exponentiation is also nice.
 -- Atleast it's memory-friendly (and fast for uneven bernoullis ^^).
 -- Also, full disclosure, I have no idea how to properly optimize in haskell.
 -- Still handing it in, because I spent a couple hours
 -- on it and you never know... :D
 --}
bernoulli :: Integer -> Rational
bernoulli 0 = 0
bernoulli 1 = - recip 2
bernoulli n
    | odd n = 0
    | otherwise = let primesL = primes' n in 
                      foldl' (\acc k -> acc + ((1%(k+1)) * 
                             (foldl' (\acc v -> 
                              let s = ((binom k v primesL) * (v^n)) % 1 in
                                  if even v then acc+s else acc-s) (0%1) [0..k]))
                      ) (0%1) [0..n]

{--
-- Implementatiom of the algorithm described in 
-- 'Computing Binomial Coefficients' -  P. Goetgheluck 
-- (The American Mathematical Monthl Vol. 94, No. 4).
-- TLDR, for all primes < n determine the power e such that
-- p^e is elment of the primefactorization of `n choose k`
-- (where e can be 0) and calculate the product.
-- It's not that much faster for small n or k but was fun to
-- implement.
-- Also, give list of primes as argument to avoid generating it
-- every time in the bernoulli calculation.
--}
binom :: Integer -> Integer -> [Integer] -> Integer
-- special cases are algorithmically covered but take up precious cycles
binom _ 0  _ = 1
binom n 1  _ = n
binom n k primeL = foldl' (\acc p->acc*powIn p n k')1(takeWhile (<=n) primeL)
        where k' = if k> n `div` 2 then n-k else k
              powIn p n k
                      | p > (n - k) = p
                      | p > (n `div` 2) = 1
                      | ((p*p) > n) = if ((n `mod` p) < (k `mod` p)) then p else 1
                      | otherwise = powLoop p n k 0 1 
              powLoop p n k r res
                  | n' == 0 = res
                  | a < b = powLoop p n' k' 1 $! (res*p)
                  | otherwise = powLoop p n' k' 0 res
                  where n' = n `div` p
                        k' = k `div` p
                        a  = n `mod` p
                        b  = (k `mod` p) + r

-- reuse decomposition as a prime checker and avoid building a list
decompositionPrime :: Integer -> Bool
decompositionPrime 1 = False
decompositionPrime n 
    | even n = False
    | otherwise = decHelper 3 n
        where decHelper d n
                | d*d > n = True
                | n `mod` d == 0 = False
                | otherwise = decHelper (d+2) n

-- modified to skip even numbers 
primes' :: Integer -> [Integer]
primes' n = 2:(filter (decompositionPrime) [3,5..n])