6893. Power Sums TLE

Hello all , I am trying to solve this problem [ en.wikipedia.org/wiki/Bernoulli_number1 ] . Could some one please tell me how to make this code more faster or this problem is also not solvable in Haskell .
Thank you

import Data.List
import Data.Maybe
import Ratio
import qualified Data.ByteString.Char8 as BS
binomM :: Integer -> Integer -> Integer -> Integer
binomM n k p =  pas !! ( fromIntegral n) !! ( fromIntegral k )   {--mod ( c * modInv d p ) p   where  
	a = max k ( n - k ) 
	b = min k ( n - k )
	c = foldl ( \ x y -> mod ( x * y ) p ) 1 [ n , pred n .. succ a ] 
	d = foldl ( \ x y -> mod ( x * y ) p ) 1 [ 1 .. b ]
--}
pas :: [[Integer]]
pas = iterate ( \x -> map ( flip mod 10007 ) $  1 : zipWith (+) x ( tail x )   ++ [1] ) [1]   
extendedGcd :: Integer -> Integer -> ( Integer , Integer ) 
extendedGcd a b 
   | b == 0 = ( 1 , 0 ) 
   | otherwise = ( t , s - q * t ) where 
	( q , r ) = quotRem a b 
	( s , t ) = extendedGcd b r 
modInv :: Integer -> Integer -> Integer 
modInv  a b 
   | gcd a b /= 1 = error " gcd is not 1 "
   | otherwise = d where 
	d = until ( > 0 ) ( + b  ) . fst.extendedGcd a $ b   
-- fast calculation of bernoulli numbers from http://okmij.org/ftp/Haskell/Bernoulli.hs 
powers = [2..] : map (zipWith (*) (head powers)) powers
neg_powers = 
  map (zipWith (\n x -> if n then x else -x) (iterate not True)) powers
pascal:: [[Integer]]
pascal = [1,2,1] : map (\line -> zipWith (+) (line++[0]) (0:line)) pascal
b' 0 = fromIntegral 1
b' 2 = 1%6
b' n = -(sumbn n)/(fromIntegral (n+1))
sumbn:: Int -> Rational
sumbn n = 1 - (fromIntegral (n+1)%(fromIntegral 2)) + 
    sum [ (b' i) * fromIntegral(comb (n+1) i) | i  <- [2,4 .. n-1] ]
  where
  comb n i = pascal!!(n-2)!!i  -- Binomial coefficient
bernoulli :: Integer -> Ratio Integer  
bernoulli 0 = 1
bernoulli 1 = (1%2)	
bernoulli n | odd n = 0
bernoulli n = 
   (-1)%2 
     + sum [ fromIntegral ((sum $ zipWith (*) powers (tail $ tail combs)) - 
                            fromIntegral k) %
             fromIntegral (k+1)
     | (k,combs)<- zip [2..n] pascal]
  where powers = (neg_powers!!( fromIntegral (n-1)) ) 
primes       = 2:map head (iterate sieve [3,5..])
sieve (p:xs) = [ x | x<-xs, x `rem` p /= 0 ]
b_denom twok 
  = product [ p | p <- takeWhile (<= twok1) primes, 
                  twok `rem` (p-1) == 0]
  where twok1 = twok + 1
--upto here take from http://okmij.org/ftp/Haskell/Bernoulli.hs
solve :: Integer -> [Ratio Integer]   
solve m = [ ( (%) (  binomM ( m + 1 ) k 10007 ) 1  ) * ( bernoulli k)   | k <-[0..m ] ]
helpBernoulli :: Integer -> [Integer] 
helpBernoulli k = ys where 
	xs = solve k 
	ys = map ( flip mod 10007 ) . map (  ( modInv ( k + 1 ) 10007 )* ) . map fun $ xs  
	fun  x  = if a < 0 then (   mod ( (10007 + a ) * modInv b  10007 ) 10007  ) else mod ( a * modInv b 10007 ) 10007  where 
	   a = numerator x 
	   b = denominator x
format ::Integer -> [Integer] -> String
format k []  = show 0 ++ "x^" ++ show 0 
format k t@(x:xs)  =  show x ++ "x^" ++ show k ++ " + " ++ format  ( k - 1 ) xs  
final :: Integer -> String
final k = format ( k + 1) ( helpBernoulli k )
main = interact $ unlines . map ( final  . read ) . lines