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