1a697827d6cef278bc5b7a4f444723b0104cb09c
[nofib.git] / imaginary / digits-of-e1 / Main.lhs
1 Compute the digits of "e" using continued fractions.
2 Original program due to Dale Thurston, Aug 2001
3
4 > module Main where
5 > import System.Environment (getArgs)
6
7 > type ContFrac = [Integer]
8
9 Compute the decimal representation of e progressively.
10
11 A continued fraction expansion for e is
12
13 [2,1,2,1,1,4,1,1,6,1,...]
14
15 > eContFrac :: ContFrac
16 > eContFrac = 2:aux 2 where aux n = 1:n:1:aux (n+2)
17
18 We need a general function that applies an arbitrary linear fractional
19 transformation to a legal continued fraction, represented as a list of
20 positive integers.  The complicated guard is to see if we can output a
21 digit regardless of what the input is; i.e., to see if the interval
22 [1,infinity) is mapped into [k,k+1) for some k.
23
24 > -- ratTrans (a,b,c,d) x: compute (a + bx)/(c+dx) as a continued fraction
25 > ratTrans :: (Integer,Integer,Integer,Integer) -> ContFrac -> ContFrac
26 > -- Output a digit if we can
27 > ratTrans (a,b,c,d) xs |
28 >   ((signum c == signum d) || (abs c < abs d)) && -- No pole in range
29 >   (c+d)*q <= a+b && (c+d)*q + (c+d) > a+b       -- Next digit is determined
30 >      = q:ratTrans (c,d,a-q*c,b-q*d) xs
31 >   where q = b `div` d
32 > ratTrans (a,b,c,d) (x:xs) = ratTrans (b,a+x*b,d,c+x*d) xs
33
34 Finally, we convert a continued fraction to digits by repeatedly multiplying by 10.
35
36 > toDigits :: ContFrac -> [Integer]
37 > toDigits (x:xs) = x:toDigits (ratTrans (10,0,0,1) xs)
38
39 > e :: [Integer]
40 > e = toDigits eContFrac
41
42 > main = do
43 >       [digits] <- getArgs
44 >       print (take (read digits) e)
45
46