1 {- Taken from Doug McIlroy's paper
2 "Power series, power serious"
3 JFP 9(3) May 1999
4 -}
6 module Main where
8 import System.IO
9 import Data.Ratio
10 import System.Environment (getArgs)
12 infixl 7 .*
13 infixr 5 :+:
15 default (Integer, Rational, Double)
17 main = do { (n:_) <- getArgs ;
18 let { p = read n :: Int } ;
19 putStrLn (show (extract p (sinx - sqrt (1-cosx^2)))) ;
20 putStrLn (show (extract p (sinx/cosx - revert (integral (1/(1+x^2)))))) ;
21 putStrLn (show (extract p ts)) ;
22 putStrLn (show (extract p tree))
23 }
26 -- From Section 6
27 tree = 0 :+: forest
28 forest = compose list tree
29 list = 1 :+: list
31 ts = 1 :+: ts^2
34 -- The main implementation follows
35 data Ps a = Pz | a :+: Ps a
37 extract :: Int -> Ps a -> [a]
38 extract 0 ps = []
39 extract n Pz = []
40 extract n (x :+: ps) = x : extract (n-1) ps
42 deriv:: Num a => Ps a -> Ps a
43 integral:: Fractional a => Ps a -> Ps a
44 compose:: (Eq a, Num a) => Ps a -> Ps a -> Ps a
45 revert:: (Eq a, Fractional a) => Ps a -> Ps a
46 toList:: Num a => Ps a -> [a]
47 takePs:: Num a => Int -> Ps a -> [a]
48 (.*):: Num a => a -> Ps a -> Ps a
49 x:: Num a => Ps a
50 expx, sinx, cosx:: Fractional a => Ps a
52 c .* Pz = Pz
53 c .* (f :+: fs) = c*f :+: c.*fs
55 x = 0 :+: 1 :+: Pz
57 toList Pz = [] --(0)
58 toList (f :+: fs) = f : (toList fs)
60 takePs n fs = take n (toList fs)
62 instance (Eq a, Num a) => Eq (Ps a) where --(1)
63 Pz == Pz = True
64 Pz == (f :+: fs) = f==0 && Pz==fs
65 fs == Pz = Pz==fs
66 (f :+: fs) == (g :+: gs) = f==g && fs==gs
68 instance (Show a, Num a) => Show (Ps a) where --(2)
69 showsPrec p Pz = showsPrec p [0]
70 showsPrec p fs = showsPrec p (toList fs)
72 instance Num a => Num (Ps a) where
73 negate Pz = Pz
74 negate (f :+: fs) = -f :+: -fs
76 Pz + gs = gs
77 fs + Pz = fs
78 (f :+: fs) + (g :+: gs) = f+g :+: fs+gs
80 Pz * _ = Pz
81 _ * Pz = Pz
82 (f :+: fs) * (g :+: gs) =
83 f*g :+: f.*gs + g.*fs + x*fs*gs --(3)
85 fromInteger 0 = Pz
86 fromInteger c = fromInteger c :+: Pz
88 instance (Eq a, Fractional a) => Fractional (Ps a) where
89 recip fs = 1/fs
91 Pz/Pz = error "power series 0/0"
92 Pz / (0 :+: gs) = Pz / gs
93 Pz / _ = Pz
94 (0 :+: fs) / (0 :+: gs) = fs / gs
95 (f :+: fs) / (g :+: gs) = let q = f/g in
96 q :+: (fs - q.*gs)/(g :+: gs)
98 compose Pz _ = Pz
99 compose (f :+: _) Pz = f :+: Pz
100 compose (f :+: fs) (0 :+: gs) = f :+: gs*(compose fs (0 :+: gs))
101 compose (f :+: fs) gs = (f :+: Pz) + gs*(compose fs gs) --(4)
103 revert (0 :+: fs) = rs where
104 rs = 0 :+: 1/(compose fs rs)
105 revert (f0 :+: f1 :+: Pz) = -1/f1 :+: 1/f1 :+: Pz --(5)
107 deriv Pz = Pz
108 deriv (_ :+: fs) = deriv1 fs 1 where
109 deriv1 Pz _ = Pz
110 deriv1 (f :+: fs) n = n*f :+: (deriv1 fs (n+1))
112 integral fs = 0 :+: (int1 fs 1) where --(6)
113 int1 Pz _ = Pz
114 int1 (f :+: fs) n = f/n :+: (int1 fs (n+1))
116 instance (Eq a, Fractional a) => Floating (Ps a) where
117 sqrt Pz = Pz
118 sqrt (0 :+: 0 :+: fs) = 0 :+: (sqrt fs)
119 sqrt (1 :+: fs) = qs where
120 qs = 1 + integral((deriv (1:+:fs))/(2.*qs))
122 expx = 1 + (integral expx)
123 sinx = integral cosx
124 cosx = 1 - (integral sinx)
126 --(0) Convert power series to a list; used in printing.
127 --(1) Equality works on polynomials; diverges otherwise.
128 --(2) Specifies how to print the new data type.
129 --(3) x*fs*gs replaces 0:fs*gs to avoid extra zero
130 -- at end of product of polynomials; it works
131 -- because x is a finite series.
132 --(4) This extra production works for the composition
133 -- of polynomials with non-zero constant term,
134 -- but not for infinite series.
135 --(5) Special case for reverting a linear function.
136 --(6) There is no special case for (integral Pz)
137 -- because this would defeat the property that
138 -- (integral) emits one term before evaluating
139 -- its operand--a property used in solving
140 -- differential equations.