1 {- Taken from Doug McIlroy's paper 2 "Power series, power serious" 3 JFP 9(3) May 1999 4 -} 5 6 module Main where 7 8 import IO 9 import Ratio 10 import System.Environment (getArgs) 11 12 infixl 7 .* 13 infixr 5 :+: 14 15 default (Integer, Rational, Double) 16 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 } 24 25 26 -- From Section 6 27 tree = 0 :+: forest 28 forest = compose list tree 29 list = 1 :+: list 30 31 ts = 1 :+: ts^2 32 33 34 -- The main implementation follows 35 data Ps a = Pz | a :+: Ps a 36 37 extract :: Int -> Ps a -> [a] 38 extract 0 ps = [] 39 extract n Pz = [] 40 extract n (x :+: ps) = x : extract (n-1) ps 41 42 deriv:: Num a => Ps a -> Ps a 43 integral:: Fractional a => Ps a -> Ps a 44 compose:: Num a => Ps a -> Ps a -> Ps a 45 revert:: 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 51 52 c .* Pz = Pz 53 c .* (f :+: fs) = c*f :+: c.*fs 54 55 x = 0 :+: 1 :+: Pz 56 57 toList Pz = [] --(0) 58 toList (f :+: fs) = f : (toList fs) 59 60 takePs n fs = take n (toList fs) 61 62 instance 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 67 68 instance Num a => Show (Ps a) where --(2) 69 showsPrec p Pz = showsPrec p [0] 70 showsPrec p fs = showsPrec p (toList fs) 71 72 instance Num a => Num (Ps a) where 73 negate Pz = Pz 74 negate (f :+: fs) = -f :+: -fs 75 76 Pz + gs = gs 77 fs + Pz = fs 78 (f :+: fs) + (g :+: gs) = f+g :+: fs+gs 79 80 Pz * _ = Pz 81 _ * Pz = Pz 82 (f :+: fs) * (g :+: gs) = 83 f*g :+: f.*gs + g.*fs + x*fs*gs --(3) 84 85 fromInteger 0 = Pz 86 fromInteger c = fromInteger c :+: Pz 87 88 instance Fractional a => Fractional (Ps a) where 89 recip fs = 1/fs 90 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) 97 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) 102 103 revert (0 :+: fs) = rs where 104 rs = 0 :+: 1/(compose fs rs) 105 revert (f0 :+: f1 :+: Pz) = -1/f1 :+: 1/f1 :+: Pz --(5) 106 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)) 111 112 integral fs = 0 :+: (int1 fs 1) where --(6) 113 int1 Pz _ = Pz 114 int1 (f :+: fs) n = f/n :+: (int1 fs (n+1)) 115 116 instance 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)) 121 122 expx = 1 + (integral expx) 123 sinx = integral cosx 124 cosx = 1 - (integral sinx) 125 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. 141