{-# LANGUAGE ScopedTypeVariables, FlexibleContexts #-}
-- | Univariate polynomials parametrised by the variable name.
module Algebra.UPoly
( UPoly(..)
, deg
, Qx, x
, toUPoly, monomial
, lt, deriv
, cont, isPrimitive
, toPrimitive, propToPrimitive, gaussLemma
, gcdUPolyWitness
, sqfr, sqfrDec
) where
import Data.List
import Test.QuickCheck
import Control.Monad (liftM)
import Algebra.TypeChar.Char hiding (Z,Q)
import Algebra.Structures.Field
import Algebra.Structures.FieldOfFractions
import Algebra.Structures.EuclideanDomain
import Algebra.Structures.ExplicitUnits
import Algebra.Structures.BezoutDomain
import Algebra.Structures.GCDDomain
import Algebra.Structures.PruferDomain
-- import Algebra.Structures.StronglyDiscrete
import Algebra.Ideal
import Algebra.Z
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
isZeroPoly :: CommutativeRing r => UPoly r x -> Bool
isZeroPoly (UP []) = True
isZeroPoly _ = False
-- | 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 (\x -> sumRing . replicate x) [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 . norm
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
norm (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 | norm g <= norm r =
qr (q <+> monomial (lt r </> lt g) (norm r - norm g))
(r <-> monomial (lt r </> lt g) (norm r - norm g) <*> g)
| otherwise = (q,r)
instance (Field k, Eq k) => PruferDomain (UPoly k x) where
calcUVW = calcUVW_B
instance (ExplicitUnits a, Eq a) => ExplicitUnits (UPoly a x) where
unit (UP [a]) = case unit a of
Just a' -> Just (UP [a'])
Nothing -> Nothing
unit _ = Nothing
-- 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
-------------------------------------------------------------------------------
-- Proof that if A is a GCD domain then A[x] is a GCD domain following
-- section 4.4 in A course in constructive algebra.
-- Some test polynomials:
test1 :: UPoly Z X_
test1 = toUPoly [1,2,3]
test2 :: UPoly Z X_
test2 = toUPoly [2,4,6,8,10]
test3 :: UPoly Q X_
test3 = toUPoly [inv 2, inv 3, inv 4]
-- | Compute the content of a polynomial, i.e. the gcd of the coefficients.
cont :: (GCDDomain a, Eq a) => UPoly a x -> a
cont (UP xs) = case filter (/= zero) xs of
[] -> error "cont: Can't compute the content of the zero polynomial"
xs' -> ggcd xs'
-- *Algebra.UPoly> cont test1
-- 1
-- *Algebra.UPoly> cont test2
-- 2
-- | If all coefficients are relatively prime then the polynomial is primitive.
isPrimitive :: (ExplicitUnits a, GCDDomain a, Eq a) => UPoly a x -> Bool
isPrimitive = isUnit . cont
-- *Algebra.UPoly> isPrimitive test1
-- True
-- *Algebra.UPoly> isPrimitive test2
-- False
-- | Lemma 4.2: Given a polynomial p in K[x] where K=Quot(A) we can find c in K
-- and q primitive in A[x] such that p = cq.
toPrimitive :: (GCDDomain a, Eq a)
=> UPoly (FieldOfFractions a) x
-> (FieldOfFractions a, UPoly a x)
toPrimitive p@(UP xs) = (c,q)
where
c0' = toFieldOfFractions $ productRing $ map denominator xs
c0 = inv c0'
g0 = map (fromFieldOfFractions . (c0' <*>)) xs
cg0 = toFieldOfFractions $ cont $ UP g0
c = c0 <*> cg0
q = toUPoly $ map (\x -> fromFieldOfFractions (toFieldOfFractions x </> cg0)) g0
-- *Algebra.UPoly> toPrimitive test3
-- (1/12,6+4x+3x^2)
-- Specification of toPrimitive.
propToPrimitive :: (ExplicitUnits a, GCDDomain a, Eq a) => UPoly (FieldOfFractions a) x -> Property
propToPrimitive p =
not (isZeroPoly p) ==>
p == toUPoly (map (\x -> c <*> toFieldOfFractions x) q) && isPrimitive (UP q)
where (c,UP q) = toPrimitive p
-- *Algebra.UPoly> quickCheck (propToPrimitive :: UPoly Q X_ -> Property)
-- +++ OK, passed 100 tests.
-- | Gauss lemma says that if p and q are polynomials over a GCD domain then
-- cont(pq) = cont(p) * cont(q).
gaussLemma :: (ExplicitUnits a, GCDDomain a, Eq a) => UPoly a x -> UPoly a x -> Property
gaussLemma p q =
not (isZeroPoly p) && not (isZeroPoly q) ==>
cont (p <*> q) ~= cont p <*> cont q
-- *Algebra.UPoly> quickCheck (gaussLemma :: UPoly Z X_ -> UPoly Z X_ -> Property)
-- +++ OK, passed 100 tests.
liftUPoly :: (GCDDomain a, Eq a) => UPoly a x -> UPoly (FieldOfFractions a) x
liftUPoly (UP xs) = toUPoly $ map toFieldOfFractions xs
-- | Proof that if A is a GCD domain then A[x] also is a GCD domain. This also
-- computes witnesses that the computed GCD divides the given polynomials.
gcdUPolyWitness :: (GCDDomain a, Eq a)
=> UPoly a x
-> UPoly a x
-> (UPoly a x, UPoly a x, UPoly a x)
gcdUPolyWitness p q = (constUPoly d <*> h, constUPoly x <*> a, constUPoly y <*> b)
where
(h',a',b') = gcd' (liftUPoly p) (liftUPoly q)
(_,h) = toPrimitive h'
(_,a) = toPrimitive a'
(_,b) = toPrimitive b'
(d,x,y) = gcd' (cont p) (cont q)
test4, test5 :: UPoly Z X_
test4 = toUPoly [6,7,1]
test5 = toUPoly [-6,-5,1]
-- *Algebra.UPoly> gcdUPolyWitness test4 test5
-- (1+x,6+x,-6+x)
-- *Algebra.UPoly> gcdUPolyWitness test4 test4
-- (6+7x+x^2,1,1)
-- *Algebra.UPoly> gcdUPolyWitness test1 test2
-- (1,1+2x+3x^2,2+4x+6x^2+8x^3+10x^4)
-- This does not work:
-- instance (GCDDomain a, Eq a) => GCDDomain (UPoly a x) where
-- gcd' = gcdUPolyWitness
-------------------------------------------------------------------------------
-- 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 | norm 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 :: (CommutativeRing a, Eq a) => UPoly a x -> UPoly a x -> Property
propPD s p = deg s > 1 && deg p > 1 ==> constUPoly c <*> s == p <*> q <+> r
where (c,q,r) = pseudoDivide s p