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 )
newtype Poly coeff = Poly { fromPoly :: FreeMod coeff X } deriving (Eq,Num,Show)
newtype RisingPoly coeff = RisingPoly { fromRisingPoly :: FreeMod coeff RisingF } deriving (Eq,Show)
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
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
newtype RisingF = RF Int deriving (Eq,Ord,Show)
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) ]
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..m1] ] 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..m1] , i /= j ]
denom = product [ x j x i | i<-[0..m1] , i /= j ]
m = length xs
x i = xs !! i
term i = ZMod.fromList
[ (X 1 , 1 )
, (X 0 , x i )
]