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