{- 	Taken from Doug McIlroy's paper
	"Power series, power serious"
	JFP 9(3) May 1999
-}

module Main where

import IO
import Ratio

infixl 7 .*
infixr 5 :+: 

default (Integer, Rational, Double)

main = do { putStrLn (show (extract 50 (sinx - sqrt (1-cosx^2)))) ;
	    putStrLn (show (extract 50 (sinx/cosx - revert (integral (1/(1+x^2)))))) ;
	    putStrLn (show (extract 50 ts)) ;
	    putStrLn (show (extract 50 tree)) 
	  }
	

-- From Section 6
tree   = 0 :+: forest
forest = compose list tree
list   = 1 :+: list

ts = 1 :+: ts^2
	

-- The main implementation follows
data Ps a = Pz | a :+: Ps a

extract :: Int -> Ps a -> [a]
extract 0 ps 	     = []
extract n Pz 	     = []
extract n (x :+: ps) = x : extract (n-1) ps

deriv:: Num a => Ps a -> Ps a
integral:: Fractional a => Ps a -> Ps a
compose:: Num a => Ps a -> Ps a -> Ps a
revert:: Fractional a => Ps a -> Ps a
toList:: Num a => Ps a -> [a]
takePs:: Num a => Int -> Ps a -> [a]
(.*):: Num a => a -> Ps a -> Ps a
x:: Num a => Ps a
expx, sinx, cosx:: Fractional a => Ps a

c .* Pz = Pz
c .* (f :+: fs) = c*f :+: c.*fs

x = 0 :+: 1 :+: Pz

toList Pz = []						--(0)
toList (f :+: fs) = f : (toList fs)

takePs n fs = take n (toList fs)

instance Num a => Eq (Ps a) where			--(1)
	Pz == Pz = True
	Pz == (f :+: fs) = f==0 && Pz==fs
	fs == Pz = Pz==fs
	(f :+: fs) == (g :+: gs) = f==g && fs==gs

instance Num a => Show (Ps a) where			--(2)
	showsPrec p Pz = showsPrec p [0]
	showsPrec p fs = showsPrec p (toList fs)

instance Num a => Num (Ps a) where
	negate Pz = Pz
	negate (f :+: fs) = -f :+: -fs

	Pz + gs = gs
	fs + Pz = fs
	(f :+: fs) + (g :+: gs) = f+g :+: fs+gs

	Pz * _ = Pz
	_ * Pz = Pz
	(f :+: fs) * (g :+: gs) =
		f*g :+: f.*gs + g.*fs + x*fs*gs		--(3)

	fromInteger 0 = Pz
	fromInteger c = fromInteger c :+: Pz

instance Fractional a => Fractional (Ps a) where
	recip fs = 1/fs

	Pz/Pz = error "power series 0/0"
	Pz / (0 :+: gs) = Pz / gs
	Pz / _ = Pz
	(0 :+: fs) / (0 :+: gs) = fs / gs
	(f :+: fs) / (g :+: gs) = let q = f/g in
		q :+: (fs - q.*gs)/(g :+: gs)

compose Pz _ = Pz
compose (f :+: _) Pz = f :+: Pz
compose (f :+: fs) (0 :+: gs) = f :+: gs*(compose fs (0 :+: gs))
compose (f :+: fs) gs = (f :+: Pz) + gs*(compose fs gs)	--(4)

revert (0 :+: fs) = rs where
	rs = 0 :+: 1/(compose fs rs)
revert (f0 :+: f1 :+: Pz) = -1/f1 :+: 1/f1 :+: Pz	--(5)

deriv Pz = Pz
deriv (_ :+: fs) = deriv1 fs 1 where
	deriv1 Pz _ = Pz
	deriv1 (f :+: fs) n = n*f :+: (deriv1 fs (n+1))

integral fs = 0 :+: (int1 fs 1) where			--(6)
	int1 Pz _ = Pz
	int1 (f :+: fs) n = f/n :+: (int1 fs (n+1))

instance Fractional a => Floating (Ps a) where
	sqrt Pz = Pz
	sqrt (0 :+: 0 :+: fs) = 0 :+: (sqrt fs)
	sqrt (1 :+: fs) = qs where
		qs = 1 + integral((deriv (1:+:fs))/(2.*qs))

expx = 1 + (integral expx)
sinx = integral cosx
cosx  = 1 - (integral sinx)

--(0) Convert power series to a list; used in printing.
--(1) Equality works on polynomials; diverges otherwise.
--(2) Specifies how to print the new data type.
--(3) x*fs*gs replaces 0:fs*gs to avoid extra zero
--    at end of product of polynomials; it works
--    because x is a finite series.
--(4) This extra production works for the composition
--    of polynomials with non-zero constant term,
--    but not for infinite series.
--(5) Special case for reverting a linear function.
--(6) There is no special case for (integral Pz)
--    because this would defeat the property that
--    (integral) emits one term before evaluating
--    its operand--a property used in solving
--    differential equations.