-- | Univariate polynomials

{-# LANGUAGE GeneralizedNewtypeDeriving #-}
module Math.RootLoci.Algebra.Polynomial where

--------------------------------------------------------------------------------

import Data.Array ( assocs ) 

import Math.Combinat.Numbers

import Math.RootLoci.Misc

import qualified Math.RootLoci.Algebra.FreeMod as ZMod
import Math.RootLoci.Algebra.FreeMod ( FreeMod , ZMod , QMod )

--------------------------------------------------------------------------------
-- * Polynomials

-- | Standard univariate polynomials
newtype Poly coeff = Poly { fromPoly :: FreeMod coeff X } deriving (Eq,Num,Show)

-- | Univariate polynomials using /rising factorials/ as a basis function
newtype RisingPoly  coeff = RisingPoly  { fromRisingPoly  :: FreeMod coeff RisingF }  deriving (Eq,Show)

-- | Univariate polynomials using /falling factorials/ as a basis function
newtype FallingPoly coeff = FallingPoly { fromFallingPoly :: FreeMod coeff FallingF } deriving (Eq,Show)

instance (Num c, Show c, Eq c, IsSigned c) => Pretty (Poly        c) where pretty (Poly        p) = pretty p 
instance (Num c, Show c, Eq c, IsSigned c) => Pretty (RisingPoly  c) where pretty (RisingPoly  p) = pretty p 
instance (Num c, Show c, Eq c, IsSigned c) => Pretty (FallingPoly c) where pretty (FallingPoly p) = pretty p 

--------------------------------------------------------------------------------
-- * Monomials 

-- | A power of @x@ (that is, a monomial of the form @x^i@)
newtype X = X Int deriving (Eq,Ord,Show)

instance Monoid X where
  mempty = X 0
  mappend (X e) (X f) = X (e+f)

instance Pretty X where
  pretty (X e) = case e of
    0 -> "1"
    1 -> "x"
    _ -> "x^" ++ show e

--------------------------------------------------------------------------------
-- * Rising and falling factorials 

-- | Rising factorial @x^(k) = x(x+1)(x+2)...(x+k-1)@
newtype RisingF = RF Int deriving (Eq,Ord,Show)

-- | Falling factorial @x_(k) = x(x-1)(x-2)...(x-k+1)@
newtype FallingF = FF Int deriving (Eq,Ord,Show)

instance Pretty RisingF where
  pretty (RF k) = case k of
    0 -> "1"
    1 -> "x"
    _ -> "x^(" ++ show k ++ ")"

instance Pretty FallingF where
  pretty (FF k) = case k of
    0 -> "1"
    1 -> "x"
    _ -> "x_(" ++ show k ++ ")"

risingPoly :: RisingF -> Poly Integer
risingPoly (RF k) = Poly $ ZMod.fromList
  [ (X p, abs c) | (p,c) <- assocs (signedStirling1stArray k) ]

fallingPoly :: FallingF -> Poly Integer
fallingPoly (FF k) = Poly $ ZMod.fromList
  [ (X p,     c) | (p,c) <- assocs (signedStirling1stArray k) ]

--------------------------------------------------------------------------------
-- * Lagrange interpolation

lagrangeInterp :: [(Rational,Rational)] -> Poly Rational
lagrangeInterp = Poly . lagrangeInterp'

lagrangeInterp' :: [(Rational,Rational)] -> QMod X
lagrangeInterp' xys = final where
  final = ZMod.sum [ ZMod.scale (ys!!j) (lagrangePoly' xs j) | j<-[0..m-1] ]  where
  m = length xys
  (xs,ys) = unzip xys

lagrangePoly' :: [Rational] -> Int -> QMod X
lagrangePoly' xs j = ZMod.scale (1/denom) numer where
  numer  = ZMod.product [ term i    | i<-[0..m-1] , i /= j ]
  denom  = product      [ x j - x i | i<-[0..m-1] , i /= j ]
  m      = length xs
  x i    = xs !! i
  term i = ZMod.fromList 
    [ (X 1 ,     1 )
    , (X 0 , - x i )
    ]

--------------------------------------------------------------------------------