```{-# 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 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

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
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
--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
```