f408d2f8400eeee95f31eb672d82c05ae5de40ab
[nofib.git] / imaginary / bernouilli / Main.hs
1 -- There was a lot of discussion about various ways of computing
2 -- Bernouilli numbers (whatever they are) on haskell-cafe in March 2003
3 -- Here's one of the programs.
4
5 -- It's not a very good test, I suspect, because it manipulates big integers,
6 -- and so probably spends most of its time in GMP.
7
8 import Data.Ratio
9 import System.Environment
10
11 -- powers = [[r^n | r<-[2..]] | n<-1..]
12 -- type signature required for compilers lacking the monomorphism restriction
13 powers :: [[Integer]]
14 powers = [2..] : map (zipWith (*) (head powers)) powers
15
16 -- powers = [[(-1)^r * r^n | r<-[2..]] | n<-1..]
17 -- type signature required for compilers lacking the monomorphism restriction
18 neg_powers :: [[Integer]]
19 neg_powers =
20 map (zipWith (\n x -> if n then x else -x) (iterate not True)) powers
21
22 pascal:: [[Integer]]
23 pascal = [1,2,1] : map (\line -> zipWith (+) (line++[0]) (0:line)) pascal
24
25 bernoulli 0 = 1
26 bernoulli 1 = -(1%2)
27 bernoulli n | odd n = 0
28 bernoulli n =
29 (-1)%2
30 + sum [ fromIntegral ((sum $ zipWith (*) powers (tail $ tail combs)) -
31 fromIntegral k) %
32 fromIntegral (k+1)
33 | (k,combs)<- zip [2..n] pascal]
34 where powers = (neg_powers!!(n-1))
35
36 main = do
37 [arg] <- getArgs
38 let n = (read arg)::Int
39 putStr $ "Bernoulli of " ++ (show n) ++ " is "
40 print (bernoulli n)