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
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
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 (*>) [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 . d
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
d (UP ps) = fromIntegral (length ps) 1
quotientRemainder f g = qr zero f
where
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)
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 | d 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 :: 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