| comb is only called from here: | sumbn n = sum [ bernoulli i * fromIntegral(comb (n+1) i) | i | <- [0 .. n-1] ] Probably I misunderstand what "bernoulli i" stands for. If it is meant Bernoulli number B_i, http://mathworld.wolfram.com/BernoulliNumber.html then the above expression is quite inefficient. Bernoulli numbers with odd indices are all zero, except B_1, which is -1/2. Therefore, the above expression ought to be written as sumbn n = 1 - (fromIntegral (n+1)%(fromIntegral 2)) + sum [ (b' i) * fromIntegral(comb (n+1) i) | i <- [2,4 .. n-1] ] It appears that you compute the sumbn to obtain B_n from the equality sum [bernoulli i * (comb (n+1) i) | i<-[0..n]] == 0 Have you tried bernoulli n = sum [ (sumib k n) % (k+1) | k<-[1..n]] where sumib k n = sum [ (comb k r) * r^n | r<-[2..k]] - sum [ (comb k r) * r^n | r<-[1,3..k]] the advantage of the latter series is that sumib is an Integer, rather than a rational. The powers of r can be memoized. Here's the code -- powers = [[r^n | r<-[1..]] | n<-[1..]] powers = [1..] : map (zipWith (*) (head powers)) powers -- neg_powers = [[(-1)^r * r^n | r<-[1..]] | n<-[1..]] neg_powers = map (zipWith (\n x -> if n then x else -x) (iterate not False)) powers pascal:: [[Integer]] pascal = [1,1] : map (\line -> zipWith (+) (line++[0]) (0:line)) pascal bernoulli n = sum [ fromIntegral (sumib k n) % fromIntegral (k+1) | k<-[1..n]] sumib:: Int -> Int -> Integer sumib k n = sum $ zipWith (*) (neg_powers!!(n-1)) (tail $ pascal!!(k-1)) This code seems to compute 'bernoulli 82' in less then a second, in Hugs (on a 2GHz Pentium IV).