-- Power series

data Power a = Series [a]

instance Add a => Add (Power a) where
  (Series (x:xs)) + (Series (y:ys)) = Series ((x+y):zs) where
       Series zs = (Series xs) + (Series ys)
  (Series (x:xs)) - (Series (y:ys)) = Series ((x-y):zs) where
       Series zs = (Series xs) - (Series ys)
  negate (Series (x:xs)) = Series ((-x):ys) where 
       Series ys = - (Series xs)
  zero = Series zeros 

instance LeftMul a b => LeftMul a (Power b) where
   a * (Series bs) = Series [ a*b | b <- bs ]

instance Div a b => Div (Power a) b where
   (Series as) / b = Series [ a/b | a <- as ]

instance (LeftMul a b, Add b) => LeftMul (Power a) (Power b) where
   (Series as) * (Series bs) = Series cs where
      cs = [ sum [ a*b | (a,b) <- zip (rtake n as) bs ] | n <- [1..] ]

instance ( Div (Power a) a, LeftMul (Power a) (Power a),
           Mult (Power a), Add (Power a), 
           Add a) => Div (Power a) (Power a) where
   xs / ys = xs * (invert ys)

invert :: (Mult (Power a), Add (Power a),
           Div (Power a) a) => (Power a) -> (Power a)
invert x = (Series (unit:ys))/x0
    where  Series (xs@(x0:_)) = x
           ys = [ f(i) | i <- [1..] ]
           f(n) = -(xs!!n)/x0 - sum [f(n-i)*(xs!!i)/x0 | i <- [1..n-1]] 

instance (Mult a, Add a) => Mult (Power a) where
    unit = const unit

zeros :: Add a => [a]
zeros = zero:zeros

const :: Add a => a -> Power a
const a = Series (a:zeros)

t :: (Mult a, Add a) => Power a
t = Series (zero:unit:zeros)

timesX :: Add a => (Power a) -> (Power a)
timesX (Series as) = Series (zero:as)

divX :: (Power a) -> (Power a)
divX   (Series as) = Series as' where _:as' = as

ptake :: Int -> (Power a) -> [a]
ptake n (Series xs) = take n xs

-- take and reverse
rtake :: Int -> [a] -> [a]
rtake n xs = f n xs [] where
   f 0 _ as          = as
   f (n+1) (x:xs) as = f n xs (x:as)

zip :: [a] -> [b] -> [(a,b)]
zip    []    _      = []
zip     _    []     = []
zip (a:as) (b:bs)   = (a,b):zip as bs

-- repeated differences
diffs :: (Div (Power a) Int,
          Add (Power a),
          Mult ((Power a,Int) -> (Power a,Int))) => [a] -> [a]
diffs xs = [ x | ( Series (x:_),_) <-
                 [(next^n) (Series xs,0) | n <- [0..]]] 
      where next (as,m) = let m' = m+1 in (((divX as) - as) / m',m')

rtake0 :: Add a => Int -> [a] -> [a]
rtake0 n xs = (rtake n xs) ++ zeros

undiff :: (LeftMul Int (Power a),
           Add (Power a)) => Int -> [a] -> [a]
undiff n xs = ys where 
      Series ys = z - n*(timesX z)
      z         = Series xs

-- rebuild coefficients
undiffs :: (LeftMul Int (Power a),
           Add (Power a)) => Int -> Int -> [a] -> [a]
undiffs d 0 xs     = rtake0 (d+1) xs
undiffs d (n+1) xs = (take (n+2) (undiff (d-n-1) (undiffs d n xs)))
                     ++ rtake0 (d-n-1) xs

--  interpolation for integral polynomials
coeffs :: Int -> (Int -> Int) -> [Int]
coeffs d f = rtake (d+1) cs where
         cs = undiffs d d (diffs [ f(n) | n <- [0..] ])

print :: (Text a, Add a, Mult a, Ord a) => Char -> a -> [a] -> String
print c x [] = show x
print c x xs | x == zero = pr c 1 xs
             | otherwise = (show x)++(spr c 1 xs)
    where
      pr c _ []                 = "0"
      pr c n (x:xs) | x == zero = pr c (n+1) xs
                    | otherwise = (t x c n)++(spr c (n+1) xs)
      spr c _ []                = ""
      spr c n (x:xs) | x == zero = spr c (n+1) xs
                     | otherwise = (if x < zero then "" else "+")++
                                   (t x c n)++(spr c (n+1) xs)
      t x c n = (if x == unit then "" else
                if x == -unit then "-" else show x)++[c]++
                (if n == 1 then "" else "^"++show(n))

-- pretty print 
display d f = print 'x' c0 cs where
       c0:cs = coeffs d f
