module Polynomial where

import Data.List (unfoldr, )
import Data.List.HT (viewL, )


type T a = [a]

{-| polynomial evaluation -}
horner :: Num a => [a] -> a -> a
horner p x = foldr (\c sum' -> c+x*sum') 0 p

{-| Horner's scheme for arbitrary rings -}
hornerRing :: c -> (a -> b -> c) -> (c -> b) -> [a] -> c
hornerRing zero add' scale' =
   foldr (\c sum' -> add' c (scale' sum')) zero

{-| express horner using hornerRing -}
hornerRingPoly :: Num a => [a] -> a -> a
hornerRingPoly p x = hornerRing 0 (+) (x*) p


{-| Evaluation of polynomial and its derivatives.
    n-th derivative divided by (factorial n). -}
hornerMulti :: Num a => [a] -> a -> [a]
hornerMulti p x =
   unfoldr (viewL . scanr1 (\c sum' -> c+x*sum')) p

fromScalar :: a -> [a]
fromScalar = (:[])

-- | add two polynomials or series
add :: Num a => [a] -> [a] -> [a]
{- zipWith (+) would cut the resulting list
   to the length of the shorter operand -}
add [] ys = ys
add xs [] = xs
add (x:xs) (y:ys) = x+y : add xs ys

-- | subtract two polynomials or series
sub :: Num a => [a] -> [a] -> [a]
sub [] ys = map negate ys
sub xs [] = xs
sub (x:xs) (y:ys) = x-y : sub xs ys

neg :: Num a => [a] -> [a]
neg = map negate

-- | scale a polynomial or series by a factor
scale :: Num a => a -> [a] -> [a]
scale s = map (s*)

-- | multiply a polynomial with a monomial
shift :: (Num a, Eq a) => Int -> [a] -> [a]
shift n x =
   if n>=0
   then replicate n 0 ++ x
   else if 0 == head x
        then shift (n+1) (tail x)
        else error "shift impossible"

-- | multiply a polynomial with linear factor
mulLinearFactor :: Num a => a -> [a] -> [a]
mulLinearFactor _ [] = []
mulLinearFactor x yp@(y:ys) =
   x*y : add (scale x ys) yp

mulRevLinearFactor :: Num a => a -> [a] -> [a]
mulRevLinearFactor _ [] = []
mulRevLinearFactor x yp@(y:ys) =
   y : add ys (scale x yp)


translate :: Num a => a -> [a] -> [a]
translate x =
   foldr (\c p -> sub (c:p) (scale x p)) []

shrink :: Num a => a -> [a] -> [a]
shrink k =
   zipWith (*) (iterate (k*) 1)


-- | multiply two polynomials or series
mul :: Num a => [a] -> [a] -> [a]
{- prevent from generation of many zeros
   if the first operand is the empty list -}
mul [] = const []
mul xs = foldr (\y zs -> add (scale y xs) (0:zs)) []


fromRoots :: Num a => [a] -> [a]
fromRoots =
   foldl (flip mulLinearFactor) [1] . map negate

fromReciprocalRoots :: Num a => [a] -> [a]
fromReciprocalRoots =
   foldl (flip mulRevLinearFactor) [1] . map negate


progression :: Num a => [a]
progression = iterate (1+) 1


differentiate :: (Num a) => [a] -> [a]
differentiate x = zipWith (*) (tail x) progression

integrate :: (Fractional a) => a -> [a] -> [a]
integrate c x = c : zipWith (/) x progression