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.Ideal
import Algebra.Z
import Algebra.Q
newtype CommutativeRing r => UPoly r x = UP [r]
deriving (Eq,Ord)
deg :: CommutativeRing r => UPoly r x -> Integer
deg (UP xs) | length xs < 2 = 0
| otherwise = toInteger (length xs) 1
type Qx = UPoly Q X_
x :: Qx
x = UP [zero,one]
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
monomial :: CommutativeRing r => r -> Integer -> UPoly r x
monomial a i = UP $ replicate (fromInteger i) zero ++ [a]
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]
deriv :: CommutativeRing r => UPoly r x -> UPoly r x
deriv (UP ps) = UP $ zipWith (\x -> sumRing . replicate x) [1..] (tail ps)
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
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
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 ps1) (length qs1) : m ps qs (r+1)
c (1) _ _ _ = zero
c r k m n | r > m || k > n = c (r1) (k+1) m n
| otherwise = ps !! r <*> qs !! k <+> c (r1) (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
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
instance (Field k, Eq k) => EuclideanDomain (UPoly k x) where
norm (UP ps) = fromIntegral (length ps) 1
quotientRemainder f g = qr zero f
where
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
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]
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'
isPrimitive :: (ExplicitUnits a, GCDDomain a, Eq a) => UPoly a x -> Bool
isPrimitive = isUnit . cont
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
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
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
liftUPoly :: (GCDDomain a, Eq a) => UPoly a x -> UPoly (FieldOfFractions a) x
liftUPoly (UP xs) = toUPoly $ map toFieldOfFractions xs
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]
sqfr :: (Num k, Field k) => UPoly k x -> UPoly k x
sqfr f = f `quotient` euclidAlg f f'
where f' = deriv f
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
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 (mn)
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 (m2n)
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