-- Power series

-- revsegs [x0,..] = [[x0],[x1,x0],....]
revsegs :: [a] -> [[a]]    
revsegs (x0:xs) = [x0]:f xs [x0] where
         f []      _   = []
         f (y0:ys) acc = chunk:f ys chunk where
           chunk = y0:acc

-- segs [x0,..] = [[x0],[x0,x1],....]
segs :: [a] -> [[a]]
segs xs = [take n xs | n<-[1..] ]

pmul :: (Mult a, Add a) => [a] -> [a] -> [a]
pmul xs ys = zipWith innerproduct (revsegs xs) (segs ys) where
         innerproduct u v = sum  (zipWith (*) u v)

data Pseries a = Series [a]

type Binop a = a -> a -> a

constant :: Add a => a -> (Pseries a)
constant x = Series (x:noughts) where noughts = nought:noughts

instance Add a => Add (Pseries a) where
  (Series xs) + (Series ys) = Series (zipWith (+) xs ys)
  (Series xs) - (Series ys) = Series (zipWith (-) xs ys)
  - (Series xs) = Series (map negate xs)
  nought = Series noughts where noughts = nought:noughts

instance (Mult a, Add a) => LeftMul (Pseries a) (Pseries a) where
   (Series xs) * (Series ys) = Series (pmul xs ys)

instance (Mult a, Add a) => Mult (Pseries a) where
   one = Series (one:noughts) where noughts = nought:noughts

subst :: (Eq a, Mult (Pseries a), Add (Pseries a)) => Binop (Pseries a)
subst (Series xs) us = if 
         consterm us then error "constant term must be nought"
                       else
         foldr1 (+)  ss where
           ss = zipWith (*) (map constant xs) [ us^n | n <- [0..] ] 
           consterm (Series (c:_)) = c /= nought

upto :: Int -> (Pseries a) -> [a]
upto n (Series xs) = take (n+1) xs



