module MathObj.Polynomial
(T, fromCoeffs, coeffs,
showsExpressionPrec, const,
evaluate, evaluateCoeffVector, evaluateArgVector,
compose, equal, add, sub, negate,
horner, hornerCoeffVector, hornerArgVector,
shift, unShift,
mul, scale, divMod, divModRev,
tensorProduct, tensorProductAlt,
mulShear, mulShearTranspose,
progression, differentiate, integrate, integrateInt,
fromRoots, alternate, reverse,
translate, dilate, shrink, )
where
import qualified Algebra.Differential as Differential
import qualified Algebra.VectorSpace as VectorSpace
import qualified Algebra.Module as Module
import qualified Algebra.Vector as Vector
import qualified Algebra.Field as Field
import qualified Algebra.PrincipalIdealDomain as PID
import qualified Algebra.Units as Units
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 Algebra.Indexable as Indexable
import Algebra.Module((*>))
import Algebra.ZeroTestable(isZero)
import Control.Monad (liftM, )
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 Test.QuickCheck (Arbitrary(arbitrary,coarbitrary))
import qualified Prelude as P98
import qualified PreludeBase as P
import qualified NumericPrelude as NP
import PreludeBase hiding (const, reverse, )
import NumericPrelude hiding (divMod, negate, stdUnit, )
newtype T a = Cons {coeffs :: [a]}
fromCoeffs :: [a] -> T a
fromCoeffs = lift0
lift0 :: [a] -> T a
lift0 = Cons
lift1 :: ([a] -> [a]) -> (T a -> T a)
lift1 f (Cons x0) = Cons (f x0)
lift2 :: ([a] -> [a] -> [a]) -> (T a -> T a -> T a)
lift2 f (Cons x0) (Cons x1) = Cons (f x0 x1)
instance Functor T where
fmap f (Cons xs) = Cons (map f xs)
plusPrec, appPrec :: Int
plusPrec = 6
appPrec = 10
instance (Show a) => Show (T a) where
showsPrec p (Cons xs) =
showParen (p >= appPrec) (showString "Polynomial.fromCoeffs " . shows xs)
showsExpressionPrec :: (Show a, ZeroTestable.C a, Additive.C a) =>
Int -> String -> T a -> String -> String
showsExpressionPrec p var poly =
if isZero poly
then showString "0"
else
let terms = filter (not . isZero . fst)
(zip (coeffs poly) monomials)
monomials = id :
showString "*" . showString var :
map (\k -> showString "*" . showString var
. showString "^" . shows k)
[(2::Int)..]
showsTerm x showsMon = showsPrec (plusPrec+1) x . showsMon
in showParen (p > plusPrec)
(foldl (.) id $ List.intersperse (showString " + ") $
map (uncurry showsTerm) terms)
horner :: Ring.C a => a -> [a] -> a
horner x = foldr (\c val -> c+x*val) zero
hornerCoeffVector :: Module.C a v => a -> [v] -> v
hornerCoeffVector x = foldr (\c val -> c+x*>val) zero
hornerArgVector :: (Module.C a v, Ring.C v) => v -> [a] -> v
hornerArgVector x = foldr (\c val -> c*>one+val*x) zero
evaluate :: Ring.C a => T a -> a -> a
evaluate (Cons y) x = horner x y
evaluateCoeffVector :: Module.C a v => T v -> a -> v
evaluateCoeffVector (Cons y) x = hornerCoeffVector x y
evaluateArgVector :: (Module.C a v, Ring.C v) => T a -> v -> v
evaluateArgVector (Cons y) x = hornerArgVector x y
compose :: (Ring.C a) => T a -> T a -> T a
compose (Cons x) y = horner y (map const x)
normalize :: (ZeroTestable.C a) => [a] -> [a]
normalize = dropWhileRev isZero
shift :: (Additive.C a) => [a] -> [a]
shift [] = []
shift l = zero : l
unShift :: [a] -> [a]
unShift [] = []
unShift (_:xs) = xs
const :: a -> T a
const x = lift0 [x]
equal :: (Eq a, ZeroTestable.C a) => [a] -> [a] -> Bool
equal x y = and (zipWithOverlap isZero isZero (==) x y)
instance (Eq a, ZeroTestable.C a) => Eq (T a) where
(Cons x) == (Cons y) = equal x y
instance (Indexable.C a, ZeroTestable.C a) => Indexable.C (T a) where
compare = Indexable.liftCompare coeffs
instance (ZeroTestable.C a) => ZeroTestable.C (T a) where
isZero (Cons x) = isZero x
add, sub :: (Additive.C a) => [a] -> [a] -> [a]
add = (+)
sub = ()
negate :: (Additive.C a) => [a] -> [a]
negate = map NP.negate
instance (Additive.C a) => Additive.C (T a) where
(+) = lift2 add
() = lift2 sub
zero = lift0 []
negate = lift1 negate
scale :: Ring.C a => a -> [a] -> [a]
scale s = map (s*)
instance Vector.C T where
zero = zero
(<+>) = (+)
(*>) = Vector.functorScale
instance (Module.C a b) => Module.C a (T b) where
(*>) x = lift1 (x *>)
instance (Field.C a, Module.C a b) => VectorSpace.C a (T b)
tensorProduct :: Ring.C a => [a] -> [a] -> [[a]]
tensorProduct = outerProduct (*)
tensorProductAlt :: Ring.C a => [a] -> [a] -> [[a]]
tensorProductAlt xs ys = map (flip scale ys) xs
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) []
mulShear :: Ring.C a => [a] -> [a] -> [a]
mulShear xs ys = map sum (shear (tensorProduct xs ys))
mulShearTranspose :: Ring.C a => [a] -> [a] -> [a]
mulShearTranspose xs ys = map sum (shearTranspose (tensorProduct xs ys))
instance (Ring.C a) => Ring.C (T a) where
one = const one
fromInteger = const . fromInteger
(*) = lift2 mul
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 =
let (y0:ys) = dropWhile isZero y
aux xs' =
forcePair .
switchL
([], xs')
(P.const $
let (x0:xs) = xs'
q0 = x0/y0
in mapFst (q0:) . aux (sub xs (scale q0 ys)))
in if isZero y
then error "MathObj.Polynomial: division by zero"
else aux x (drop (length y 1) x)
instance (ZeroTestable.C a, Field.C a) => Integral.C (T a) where
divMod (Cons x) (Cons y) =
let (d,m) = divMod x y
in (Cons d, Cons m)
stdUnit :: (ZeroTestable.C a, Ring.C a) => [a] -> a
stdUnit x = case normalize x of
[] -> one
l -> last l
instance (ZeroTestable.C a, Field.C a) => Units.C (T a) where
isUnit (Cons []) = False
isUnit (Cons (x0:xs)) = not (isZero x0) && all isZero xs
stdUnit (Cons x) = const (stdUnit x)
stdUnitInv (Cons x) = const (recip (stdUnit x))
instance (ZeroTestable.C a, Field.C a) => PID.C (T a)
progression :: Ring.C a => [a]
progression = iterate (one+) one
differentiate :: (Ring.C a) => [a] -> [a]
differentiate = zipWith (*) progression . drop 1
integrate :: (Field.C a) => a -> [a] -> [a]
integrate c x = c : zipWith (/) x progression
integrateInt :: (ZeroTestable.C a, Integral.C a) => a -> [a] -> [a]
integrateInt c x =
c : zipWith Integral.safeDiv x progression
instance (Ring.C a) => Differential.C (T a) where
differentiate = lift1 differentiate
fromRoots :: (Ring.C a) => [a] -> T a
fromRoots = Cons . foldl (flip mulLinearFactor) [1]
mulLinearFactor :: Ring.C a => a -> [a] -> [a]
mulLinearFactor x yt@(y:ys) = Additive.negate (x*y) : yt scale x ys
mulLinearFactor _ [] = []
alternate :: Additive.C a => [a] -> [a]
alternate = zipWith ($) (cycle [id, Additive.negate])
reverse :: Additive.C a => T a -> T a
reverse = lift1 alternate
translate :: Ring.C a => a -> T a -> T a
translate d =
lift1 $ foldr (\c p -> [c] + mulLinearFactor d p) []
shrink :: Ring.C a => a -> T a -> T a
shrink k =
lift1 $ zipWith (*) (iterate (k*) one)
dilate :: Field.C a => a -> T a -> T a
dilate = shrink . Field.recip
instance (Arbitrary a, ZeroTestable.C a) => Arbitrary (T a) where
arbitrary = liftM (fromCoeffs . normalize) arbitrary
coarbitrary = undefined
legacyInstance :: a
legacyInstance =
error "legacy Ring.C instance for simple input of numeric literals"
instance (Ring.C a, Eq a, Show a, ZeroTestable.C a) => P98.Num (T a) where
fromInteger = const . fromInteger
negate = Additive.negate
(+) = legacyInstance
(*) = legacyInstance
abs = legacyInstance
signum = legacyInstance
instance (Field.C a, Eq a, Show a, ZeroTestable.C a) => P98.Fractional (T a) where
fromRational = const . fromRational
(/) = legacyInstance