{-# LANGUAGE NoImplicitPrelude #-}
module MathObj.Polynomial.Core (
horner, hornerCoeffVector, hornerArgVector,
normalize,
shift, unShift,
equal,
add, sub, negate,
scale, collinear,
tensorProduct, tensorProductAlt,
mul, mulShear, mulShearTranspose,
divMod, divModRev,
stdUnit,
progression, differentiate, integrate, integrateInt,
mulLinearFactor,
alternate, dilate, shrink,
) where
import qualified Algebra.Module as Module
import qualified Algebra.Field as Field
import qualified Algebra.IntegralDomain as Integral
import qualified Algebra.Ring as Ring
import qualified Algebra.Additive as Additive
import qualified Algebra.ZeroTestable as ZeroTestable
import qualified Data.List as List
import NumericPrelude.List (zipWithOverlap, )
import Data.Tuple.HT (mapPair, mapFst, forcePair, )
import Data.List.HT
(dropWhileRev, switchL, shear, shearTranspose, outerProduct, )
import qualified NumericPrelude.Base as P
import qualified NumericPrelude.Numeric as NP
import NumericPrelude.Base hiding (const, reverse, )
import NumericPrelude.Numeric hiding (divMod, negate, stdUnit, )
{-# INLINE horner #-}
horner :: Ring.C a => a -> [a] -> a
horner x = foldr (\c val -> c+x*val) zero
{-# INLINE hornerCoeffVector #-}
hornerCoeffVector :: Module.C a v => a -> [v] -> v
hornerCoeffVector x = foldr (\c val -> c+x*>val) zero
{-# INLINE hornerArgVector #-}
hornerArgVector :: (Module.C a v, Ring.C v) => v -> [a] -> v
hornerArgVector x = foldr (\c val -> c*>one+val*x) zero
{-# INLINE normalize #-}
normalize :: (ZeroTestable.C a) => [a] -> [a]
normalize = dropWhileRev isZero
{-# INLINE shift #-}
shift :: (Additive.C a) => [a] -> [a]
shift [] = []
shift l = zero : l
{-# INLINE unShift #-}
unShift :: [a] -> [a]
unShift [] = []
unShift (_:xs) = xs
{-# INLINE equal #-}
equal :: (Eq a, ZeroTestable.C a) => [a] -> [a] -> Bool
equal x y = and (zipWithOverlap isZero isZero (==) x y)
add, sub :: (Additive.C a) => [a] -> [a] -> [a]
add = (+)
sub = (-)
{-# INLINE negate #-}
negate :: (Additive.C a) => [a] -> [a]
negate = map NP.negate
{-# INLINE scale #-}
scale :: Ring.C a => a -> [a] -> [a]
scale s = map (s*)
collinear :: (Eq a, Ring.C a) => [a] -> [a] -> Bool
collinear (x:xs) (y:ys) =
if x==zero && y==zero
then collinear xs ys
else scale x ys == scale y xs
collinear xs ys =
all (==zero) xs && all (==zero) ys
{-# INLINE tensorProduct #-}
tensorProduct :: Ring.C a => [a] -> [a] -> [[a]]
tensorProduct = outerProduct (*)
tensorProductAlt :: Ring.C a => [a] -> [a] -> [[a]]
tensorProductAlt xs ys = map (flip scale ys) xs
{-# INLINE mul #-}
mul :: Ring.C a => [a] -> [a] -> [a]
mul [] = P.const []
mul xs = foldr (\y zs -> let (v:vs) = scale y xs in v : add vs zs) []
{-# INLINE mulShear #-}
mulShear :: Ring.C a => [a] -> [a] -> [a]
mulShear xs ys = map sum (shear (tensorProduct xs ys))
{-# INLINE mulShearTranspose #-}
mulShearTranspose :: Ring.C a => [a] -> [a] -> [a]
mulShearTranspose xs ys = map sum (shearTranspose (tensorProduct xs ys))
divMod :: (ZeroTestable.C a, Field.C a) => [a] -> [a] -> ([a], [a])
divMod x y =
mapPair (List.reverse, List.reverse) $
divModRev (List.reverse x) (List.reverse y)
divModRev :: (ZeroTestable.C a, Field.C a) => [a] -> [a] -> ([a], [a])
divModRev x y =
case dropWhile isZero y of
[] -> error "MathObj.Polynomial: division by zero"
y0:ys ->
let
aux xs' =
forcePair .
switchL
([], xs')
(P.const $
let (x0:xs) = xs'
q0 = x0/y0
in mapFst (q0:) . aux (sub xs (scale q0 ys)))
in aux x (drop (length ys) x)
{-# INLINE stdUnit #-}
stdUnit :: (ZeroTestable.C a, Ring.C a) => [a] -> a
stdUnit x = case normalize x of
[] -> one
l -> last l
{-# INLINE progression #-}
progression :: Ring.C a => [a]
progression = iterate (one+) one
{-# INLINE differentiate #-}
differentiate :: (Ring.C a) => [a] -> [a]
differentiate = zipWith (*) progression . drop 1
{-# INLINE integrate #-}
integrate :: (Field.C a) => a -> [a] -> [a]
integrate c x = c : zipWith (/) x progression
{-# INLINE integrateInt #-}
integrateInt :: (ZeroTestable.C a, Integral.C a) => a -> [a] -> [a]
integrateInt c x =
c : zipWith Integral.divChecked x progression
{-# INLINE mulLinearFactor #-}
mulLinearFactor :: Ring.C a => a -> [a] -> [a]
mulLinearFactor x yt@(y:ys) = Additive.negate (x*y) : yt - scale x ys
mulLinearFactor _ [] = []
{-# INLINE alternate #-}
alternate :: Additive.C a => [a] -> [a]
alternate = zipWith ($) (cycle [id, Additive.negate])
{-# INLINE shrink #-}
shrink :: Ring.C a => a -> [a] -> [a]
shrink k = zipWith (*) (iterate (k*) one)
{-# INLINE dilate #-}
dilate :: Field.C a => a -> [a] -> [a]
dilate = shrink . Field.recip