{-# LANGUAGE ScopedTypeVariables, FlexibleContexts #-} -- | Univariate polynomials parametrised by the variable name. module Algebra.UPoly ( UPoly(..) , deg , Qx, x , toUPoly, monomial , lt, deriv , sqfr, sqfrDec ) where import Data.List import Test.QuickCheck import Control.Monad (liftM) import Algebra.TypeChar.Char hiding (Q) import Algebra.Structures.Field import Algebra.Structures.BezoutDomain import Algebra.Structures.EuclideanDomain import Algebra.Structures.StronglyDiscrete import Algebra.Ideal import Algebra.Q -- | Polynomials over a commutative ring, indexed by a phantom type x that -- denote the name of the variable that the polynomial is over. For example -- UPoly Q X_ is Q[x] and UPoly Q T_ is Q[t]. newtype CommutativeRing r => UPoly r x = UP [r] deriving (Eq,Ord) -- | The degree of the polynomial. deg :: CommutativeRing r => UPoly r x -> Integer deg (UP xs) | length xs < 2 = 0 | otherwise = (toInteger $ length xs) - 1 -- | Useful shorthand for Q[x]. type Qx = UPoly Q X_ -- | The variable x in Q[x]. x :: Qx x = UP [zero,one] -- | Take a list and construct a polynomial by removing all zeroes in the end. toUPoly :: (CommutativeRing r, Eq r) => [r] -> UPoly r x toUPoly = UP . reverse . dropWhile (==zero) . reverse -- | Take an element of the ring and the degree of the desired monomial, for -- example: monomial 3 7 = 3x^7 monomial :: CommutativeRing r => r -> Integer -> UPoly r x monomial a i = UP $ replicate (fromInteger i) zero ++ [a] -- | Compute the leading term of a polynomial. lt :: CommutativeRing r => UPoly r x -> r lt (UP []) = zero lt (UP xs) = last xs constUPoly :: (CommutativeRing r, Eq r) => r -> UPoly r x constUPoly x = toUPoly [x] -- | Formal derivative of polynomials in k[x]. deriv :: CommutativeRing r => UPoly r x -> UPoly r x deriv (UP ps) = UP $ zipWith (*>) [1..] (tail ps) -- | Funny integration: integrate :: (Enum b, Field b, Integral k, Field k, Fractional b) => UPoly k x -> UPoly b x integrate (UP ps) = UP $ 0.0 : zipWith (/) (map fromIntegral ps) [1..] instance (CommutativeRing r, Eq r, Show r, Show x) => Show (UPoly r x) where show (UP []) = "0" show (UP ps) = init $ fixSign $ concat [ show' (show (undefined :: x)) p n | (p,n) <- zip ps [0..] , p /= zero ] where show' :: (CommutativeRing r, Show r) => String -> r -> Integer -> String show' x p 0 = show p ++ "+" show' x p 1 = if p == one then x ++ "+" else show p ++ x ++ "+" show' x p n = if p == one then x ++ "^" ++ show n ++ "+" else show p ++ x ++ "^" ++ show n ++ "+" fixSign [] = [] fixSign [x] = [x] fixSign ('+':'-':xs) = '-' : fixSign xs fixSign (x:xs) = x : fixSign xs instance (CommutativeRing r, Eq r, Arbitrary r) => Arbitrary (UPoly r x) where arbitrary = liftM (toUPoly . take 5) arbitrary -- Addition of polynomials. addUP :: (CommutativeRing r, Eq r) => UPoly r x -> UPoly r x -> UPoly r x addUP (UP ps) (UP qs) | length ps >= length qs = add' ps qs | otherwise = add' qs ps where add' a b = toUPoly $ zipWith (<+>) a b ++ drop (length b) a -- Multiplication of polynomials. mulUP :: (CommutativeRing r, Eq r) => UPoly r x -> UPoly r x -> UPoly r x mulUP (UP ps) (UP qs) = toUPoly $ m ps qs 0 where m ps qs r | r > length ps + length qs - 2 = [] | otherwise = c r 0 (length ps-1) (length qs-1) : m ps qs (r+1) c (-1) _ _ _ = zero c r k m n | r > m || k > n = c (r-1) (k+1) m n | otherwise = ps !! r <*> qs !! k <+> c (r-1) (k+1) m n instance (CommutativeRing r, Eq r) => Ring (UPoly r x) where (<+>) = addUP zero = UP [] one = UP [one] neg (UP ps) = UP $ map neg ps (<*>) = mulUP instance (Show r, Field r, Num r, Show x) => Num (UPoly r x) where (+) = (<+>) (-) = (<->) (*) = (<*>) abs = fromInteger . d signum = undefined -- Is it possible to define this? fromInteger x = UP [fromInteger x] instance (CommutativeRing r, Eq r) => CommutativeRing (UPoly r x) where instance (CommutativeRing r, Eq r) => IntegralDomain (UPoly r x) where -- Polynomial rings are Euclidean. instance (Field k, Eq k) => EuclideanDomain (UPoly k x) where d (UP ps) = fromIntegral (length ps) - 1 quotientRemainder f g = qr zero f where -- This is the division algorithm in k[x]. Page 39 in Cox. qr q r | d g <= d r = qr (q <+> monomial (lt r </> lt g) (d r - d g)) (r <-> monomial (lt r </> lt g) (d r - d g) <*> g) | otherwise = (q,r) -- Now that we know that the polynomial ring k[x] is a Bezout domain it is -- possible to implement membership in an ideal of k[x]. f is a member of the -- ideal <f1,...,fn> if the rest is zero when dividing f with the principal -- ideal <h>. -- instance (Field k, Eq k, Show x) => StronglyDiscrete (UPoly k x) where -- member p ps = modulo p h == zero -- where Id [h] = (\(a,_,_) -> a) $ toPrincipal ps -- Square free decomposition. -- Teo Mora; Solving Polynomial Equations Systems I. pg 69-70 -- Works only for Char 0 --TODO: Add check for char --square free associate of f -- | Square free decomposition of a polynomial. sqfr :: (Num k, Field k) => UPoly k x -> UPoly k x sqfr f = f `quotient` euclidAlg f f' where f' = deriv f -- | Distinct power factorization, aka square free decomposition sqfrDec :: (Num k, Field k) => UPoly k x -> [UPoly k x] sqfrDec f = help p q where p = euclidAlg f (deriv f) q = f `quotient` p help p q | d q < 1 = [] | otherwise = t : help (p `quotient` s) s where s = euclidAlg p q t = q `quotient` s -- | Pseudo-division of polynomials. -- -- Given s(x) and p(x) compute c, q(x) and r(x) such that: -- -- cs(x) = p(x)q(x)+r(x), deg r < deg p. pseudoDivide :: (CommutativeRing a, Eq a) => UPoly a x -> UPoly a x -> (a, UPoly a x, UPoly a x) pseudoDivide s p | m < n = (one,zero,s) | otherwise = pD (a' <*> s' <-> b' <*> xmn <*> p') 1 (b' <*> xmn) s2 where n = deg p m = deg s a = lt p a' = constUPoly a b = lt s b' = constUPoly b s' = s <-> monomial b m xmn = monomial one (m-n) p' = p <-> monomial a n s2 = a' <*> s' <-> b' <*> xmn <-> p' pD s k out1 out2 | deg s < n = (a <^> k,out1,out2) | otherwise = pD s3 (k+1) (b2xm2n <+> a' <*> out1) s3 where b2 = lt s m2 = deg s s2' = s <-> monomial b2 m2 b2xm2n = monomial b2 (m2-n) s3 = (a' <*> s) <-> (b2xm2n <*> p) propPD :: Qx -> Qx -> Property propPD s p = deg s > 1 && deg p > 1 ==> constUPoly c <*> s == p <*> q <+> r where (c,q,r) = pseudoDivide s p