{-# LANGUAGE CPP #-}
{-
	Copyright (C) 2011-2015 Dr. Alistair Ward

	This program is free software: you can redistribute it and/or modify
	it under the terms of the GNU General Public License as published by
	the Free Software Foundation, either version 3 of the License, or
	(at your option) any later version.

	This program is distributed in the hope that it will be useful,
	but WITHOUT ANY WARRANTY; without even the implied warranty of
	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
	GNU General Public License for more details.

	You should have received a copy of the GNU General Public License
	along with this program.  If not, see <http://www.gnu.org/licenses/>.
-}
{- |
 [@AUTHOR@]	Dr. Alistair Ward

 [@DESCRIPTION@]

	* Describes a <https://en.wikipedia.org/wiki/Univariate> polynomial and operations on it.

	* <https://en.wikipedia.org/wiki/Polynomial>.

	* <http://mathworld.wolfram.com/Polynomial.html>.
-}

module Factory.Data.Polynomial(
-- * Types
-- ** Type-synonyms
--	MonomialList,
-- ** Data-types
	Polynomial,
-- * Constants
	zero,
	one,
-- * Functions
	evaluate,
	getDegree,
	getLeadingTerm,
	lift,
	mod',
	normalise,
--	pruneCoefficients,
	raiseModulo,
	realCoefficientsToFrac,
	terms,
-- ** Constructors
	mkConstant,
	mkLinear,
	mkPolynomial,
-- ** Operators
	(*=),
-- ** Predicates
	areCongruentModulo,
	inAscendingOrder,
	inDescendingOrder,
--	inOrder,
	isMonic,
	isMonomial,
	isNormalised,
	isPolynomial,
--	isReduced,
	isZero
) where

import			Control.Arrow((&&&))
import qualified	Control.Arrow
import qualified	Data.List
import			Factory.Data.Monomial((<*>), (</>), (<=>), (=~))
import qualified	Factory.Data.Monomial		as Data.Monomial
import qualified	Factory.Data.QuotientRing	as Data.QuotientRing
import			Factory.Data.Ring((=*=), (=+=), (=-=))
import qualified	Factory.Data.Ring		as Data.Ring

#if MIN_VERSION_base(4,8,0)
import Prelude hiding ((<*>))	-- The "Prelude" from 'base-4.8' exports this symbol.
#endif

infixl 7 *=	-- Same as (*).

-- | The guts of a 'Polynomial'.
type MonomialList coefficient exponent	= [Data.Monomial.Monomial coefficient exponent]

{- |
	* The type of an arbitrary /univariate/ polynomial;
	actually it's more general, since it permits negative powers (<https://en.wikipedia.org/wiki/Laurent_polynomial>s).
	It can't describe /multivariate/ polynomials, which would require a list of /exponents/.
	Rather than requiring the /exponent/ to implement the /type-class/ 'Integral', this is implemented at the function-level, as required.

	* The structure permits gaps between /exponents/,
	in which /coefficients/ are inferred to be zero, thus enabling efficient representation of sparse polynomials.

	* CAVEAT: the 'MonomialList' is required to;
	be ordered by /descending/ exponent (ie. reverse <https://en.wikipedia.org/wiki/Monomial_order>);
	have had zero coefficients removed;
	and to have had /like/ terms merged;
	so the raw data-constructor isn't exported.
-}
newtype {- Integral exponent => -} Polynomial coefficient exponent	= MkPolynomial {
	getMonomialList	:: MonomialList coefficient exponent	-- ^ Accessor.
} deriving (Eq, Show)

-- | Makes /Polynomial/ a 'Data.Ring.Ring', over the /field/ composed from all possible /coefficients/; <https://en.wikipedia.org/wiki/Polynomial_ring>.
instance (
	Eq	c,
	Num	c,
	Num	e,
	Ord	e
 ) => Data.Ring.Ring (Polynomial c e) where
	MkPolynomial [] =*= _	= zero
	_ =*= MkPolynomial []	= zero
	polynomialL =*= polynomialR
--		| polynomialL == one			= polynomialR	-- Counterproductive.
--		| polynomialR == one			= polynomialL	-- Counterproductive.
		| terms polynomialL > terms polynomialR	= polynomialL `times` polynomialR
		| otherwise				= polynomialR `times` polynomialL
		where
			l `times` r	= {-# SCC "times" #-} Data.Ring.sum' (recip 2) {-TODO-} 10 {-empirical-} . map (l *=) $ getMonomialList r

	MkPolynomial [] =+= p				= p
	p =+= MkPolynomial []				= p
	MkPolynomial listL =+= MkPolynomial listR	= {-# SCC "merge" #-} MkPolynomial $ merge listL listR	where
		merge [] r			= r
		merge l []			= l
		merge l@(lh : ls) r@(rh : rs)	= case lh <=> rh of
			GT	-> lh : merge ls r
			LT	-> rh : merge l rs
			_	-> case lh `Data.Monomial.shiftCoefficient` Data.Monomial.getCoefficient rh of
				(0, _)		-> merge ls rs
				monomial	-> monomial : merge ls rs

	additiveInverse		= lift (Data.Monomial.negateCoefficient `map`)
	multiplicativeIdentity	= one
	additiveIdentity	= zero

{-
	Override the default implementation,
	in order to take advantage of the symmetry under reflection about the main diagonal,
	in the square matrix of products formed from the multiplication of each term by each term.
	Eg:
		(ax^3 + bx^2 + cx + d)^2 = [
			(a^2x^6 + abx^5 + acx^4 + adx^3) +
			(bax^5 + b^2x^4 + bcx^3 + bdx^2) +
			(cax^4 + cbx^3 + c^2x^2 + cdx) +
			(dax^3 + dbx^2 + dcx + d^2)
		]

		= (a^2x^6 + b^2x^4 + c^2x^2 + d^2) + 2 * [ax^3 * (bx^2 + cx + d) + bx^2 * (cx + d) + cx * (d)]
-}
	square (MkPolynomial [])	= zero
	square p			= Data.Ring.sum' (recip 2) {-TODO-} 10 {-empirical-} $ diagonal : corners	where
		diagonal	= {-# SCC "diagonal" #-} map Data.Monomial.square `lift` p
		corners		= {-# SCC "corners" #-} uncurry (
			zipWith (*=)
		 ) $ map MkPolynomial . init {-remove terminal null-} . Data.List.tails . tail &&& map Data.Monomial.double $ getMonomialList p

-- | Defines the ability to divide /polynomials/.
instance (
	Eq		c,
	Fractional	c,
	Num		e,
	Ord		e
 ) => Data.QuotientRing.QuotientRing (Polynomial c e)	where
{-
	Uses /Euclidian division/.
	<https://en.wikipedia.org/wiki/Polynomial_long_division>.
	<http://demonstrations.wolfram.com/PolynomialLongDivision/>.
-}
	_ `quotRem'` MkPolynomial []		= error "Factory.Data.Polynomial.quotRem':\tzero denominator."
	polynomialN `quotRem'` polynomialD	= longDivide polynomialN	where
--		longDivide :: (Fractional c, Num e, Ord e) => Polynomial c e -> (Polynomial c e, Polynomial c e)
		longDivide (MkPolynomial [])	= (zero, zero)	-- Exactly divides.
		longDivide numerator
			| Data.Monomial.getExponent quotient < 0	= (zero, numerator)	-- Indivisible remainder.
			| otherwise					= Control.Arrow.first (lift (quotient :)) $ longDivide (numerator =-= polynomialD *= quotient )
			where
--				quotient :: (Fractional c, Num e) => Data.Monomial.Monomial c e
				quotient	= getLeadingTerm numerator </> getLeadingTerm polynomialD

{- |
	* Transforms the data behind the constructor.

	* CAVEAT: similar to 'Data.Functor.fmap', but 'Polynomial' isn't an instance of 'Data.Functor.Functor' since we may want to operate on both /type-parameters/.

	* CAVEAT: the caller is required to re-'normalise' the resulting polynomial depending on the nature of the transformation of the data.
-}
lift :: (MonomialList c1 e1 -> MonomialList c2 e2) -> Polynomial c1 e1 -> Polynomial c2 e2
lift transform	= MkPolynomial . transform . getMonomialList

-- | Returns the number of non-zero terms in the polynomial.
terms :: Polynomial c e -> Int
terms (MkPolynomial l)	= length l

-- | Return the highest-degree monomial.
getLeadingTerm :: Polynomial c e -> Data.Monomial.Monomial c e
getLeadingTerm (MkPolynomial [])	= error "Factory.Data.Polynomial.getLeadingTerm:\tzero polynomial has no leading term."
getLeadingTerm (MkPolynomial (m : _))	= m

-- | Removes terms with a /coefficient/ of zero.
pruneCoefficients :: (Eq c, Num c) => Polynomial c e -> Polynomial c e
pruneCoefficients (MkPolynomial [])	= zero
pruneCoefficients p			= filter ((/= 0) . Data.Monomial.getCoefficient) `lift` p

-- | Sorts into /descending order/ of exponents, groups /like/ exponents, and calls 'pruneCoefficients'.
normalise :: (Eq c, Num c, Ord e) => Polynomial c e -> Polynomial c e
normalise	= pruneCoefficients . lift (
	map (
		foldr ((+) . Data.Monomial.getCoefficient) 0 &&& Data.Monomial.getExponent . head
	) . Data.List.groupBy (=~) . Data.List.sortBy (flip (<=>))
 )

-- | Constructs an arbitrary /zeroeth-degree polynomial/, ie. independent of the /indeterminate/.
mkConstant :: (Eq c, Num c, Num e) => c -> Polynomial c e
mkConstant 0	= zero
mkConstant c	= MkPolynomial [(c, 0)]

-- | Constructs an arbitrary /first-degree polynomial/.
mkLinear :: (Eq c, Num c, Num e)
	=> c	-- ^ Gradient.
	-> c	-- ^ Constant.
	-> Polynomial c e
mkLinear m c	= pruneCoefficients $ MkPolynomial [(m, 1), (c, 0)]

-- | Smart constructor. Constructs an arbitrary /polynomial/.
mkPolynomial :: (Eq c, Num c, Ord e) => MonomialList c e -> Polynomial c e
mkPolynomial []	= zero
mkPolynomial l	= normalise $ MkPolynomial l

-- | Constructs a /polynomial/ with zero terms.
zero :: Polynomial c e
zero	= MkPolynomial []

-- | Constructs a constant /monomial/, independent of the /indeterminate/.
one :: (Eq c, Num c, Num e) => Polynomial c e
one	= mkConstant 1

-- | True if all /exponents/ are in the order defined by the specified comparator.
inOrder :: (e -> e -> Bool) -> Polynomial c e -> Bool
inOrder comparator p
	| any ($ p) [isZero, isMonomial]	= True
	| otherwise				= and . uncurry (zipWith comparator) . (init &&& tail) . map Data.Monomial.getExponent $ getMonomialList p

-- | True if the /exponents/ of successive terms are in /ascending/ order.
inAscendingOrder :: Ord e => Polynomial c e -> Bool
inAscendingOrder	= inOrder (<=)

-- | True if the /exponents/ of successive terms are in /descending/ order.
inDescendingOrder :: Ord e => Polynomial c e -> Bool
inDescendingOrder	= inOrder (>=)

-- | True if no term has a /coefficient/ of zero.
isReduced :: (Eq c, Num c) => Polynomial c e -> Bool
isReduced	= all ((/= 0) . Data.Monomial.getCoefficient) . getMonomialList

-- | True if no term has a /coefficient/ of zero and the /exponents/ of successive terms are in /descending/ order.
isNormalised :: (Eq c, Num c, Ord e) => Polynomial c e -> Bool
isNormalised polynomial	= all ($ polynomial) [isReduced, inDescendingOrder]

{- |
	* 'True' if the /leading coefficient/ is one.

	* <https://en.wikipedia.org/wiki/Monic_polynomial#Classifications>.
-}
isMonic :: (Eq c, Num c) => Polynomial c e -> Bool
isMonic (MkPolynomial [])	= False	-- All coefficients are zero, and have therefore been removed.
isMonic p			= (== 1) . Data.Monomial.getCoefficient $ getLeadingTerm p

-- | True if there are zero terms.
isZero :: Polynomial c e -> Bool
isZero (MkPolynomial [])	= True
isZero _			= False

-- | True if there's exactly one term.
isMonomial :: Polynomial c e -> Bool
isMonomial (MkPolynomial [])	= True
isMonomial _			= False

-- | True if all /exponents/ are /positive/ integers as required.
isPolynomial :: Integral e => Polynomial c e -> Bool
isPolynomial	= all Data.Monomial.isMonomial . getMonomialList

{- |
	* 'True' if the two specified /polynomials/ are /congruent/ in /modulo/-arithmetic.

	* <http://planetmath.org/encyclopedia/PolynomialCongruence.html>.
-}
areCongruentModulo :: (Integral c, Num e, Ord e)
	=> Polynomial c e	-- ^ LHS.
	-> Polynomial c e	-- ^ RHS.
	-> c			-- ^ Modulus.
	-> Bool
areCongruentModulo _ _ 0	= error "Factory.Data.Polynomial.areCongruentModulo:\tzero modulus."
areCongruentModulo _ _ 1	= True
areCongruentModulo l r	modulus
	| l == r	= True
	| otherwise	= all ((== 0) . (`mod` modulus) . Data.Monomial.getCoefficient) . getMonomialList $ l =-= r

{- |
	* Return the /degree/ (AKA /order/) of the /polynomial/.

	* <https://en.wikipedia.org/wiki/Degree_of_a_polynomial>.

	* <http://mathworld.wolfram.com/PolynomialDegree.html>.
-}
getDegree :: Num e => Polynomial c e -> e
getDegree (MkPolynomial [])	= -1	-- CAVEAT: debatable, but makes some operations more robust and consistent.
getDegree p			= Data.Monomial.getExponent $ getLeadingTerm p

{- |
	* Scale-up the specified /polynomial/ by a constant /monomial/ factor.

	* <https://en.wikipedia.org/wiki/Scalar_multiplication>.
-}
(*=) :: (Eq c, Num c, Num e) => Polynomial c e -> Data.Monomial.Monomial c e -> Polynomial c e
polynomial *= monomial
	| Data.Monomial.getCoefficient monomial == 1	= map (`Data.Monomial.shiftExponent` Data.Monomial.getExponent monomial) `lift` polynomial
	| otherwise					= map (monomial <*>) `lift` polynomial

{- |
	* Raise a /polynomial/ to the specified positive integral power, but using /modulo/-arithmetic.

	* Whilst one could naively implement this as @(x Data.Ring.=^ n) `mod` m@, this will result in arithmetic operatons on unnecessarily big integers.
-}
raiseModulo :: (Integral c, Integral power, Num e, Ord e, Show power)
	=> Polynomial c e	-- ^ The base.
	-> power		-- ^ The exponent to which the base should be raised.
	-> c			-- ^ The modulus.
	-> Polynomial c e	-- ^ The result.
raiseModulo _ _ 0			= error "Factory.Data.Polynomial.raiseModulo:\tzero modulus."
raiseModulo _ _ 1			= zero
raiseModulo _ 0 modulus			= mkConstant $ 1 `mod` modulus
raiseModulo polynomial power modulus
	| power < 0			= error $ "Factory.Data.Polynomial.raiseModulo:\tthe result isn't guaranteed to be a polynomial, for power=" ++ show power
	| first `elem` [zero, one]	= first	-- Eg 'raiseModulo (mkPolynomial [(3,1)]) 100 3' or 'raiseModulo (mkPolynomial [(3,1),(1,0)]) 100 3'.
	| otherwise			= slave power
	where
--		first :: Integral c => Polynomial c e
		first	= polynomial `mod'` modulus

--		slave :: (Integral c, Integral power, Num e, Ord e) => power -> Polynomial c e
		slave 1	= first
		slave n	= (`mod'` modulus) . (if r == 0 {-even-} then id else (polynomial =*=)) . Data.Ring.square $ slave q {-recurse-}	where
			(q, r)	= n `quotRem` 2

-- | Reduces all the coefficients using /modular/ arithmetic.
mod' :: Integral c
	=> Polynomial c e
	-> c	-- ^ Modulus.
	-> Polynomial c e
mod' p modulus	= pruneCoefficients $ map (`Data.Monomial.mod'` modulus) `lift` p

{- |
	* Evaluate the /polynomial/ at a specific /indeterminate/.

	* CAVEAT: requires positive exponents; but it wouldn't really be a /polynomial/ otherwise.

	* If the /polynomial/ is very sparse, this may be inefficient,
	since it /memoizes/ the complete sequence of powers up to the polynomial's /degree/.
-}
evaluate :: (Num n, Integral e, Show e)
	=> n	-- ^ The /indeterminate/.
	-> Polynomial n e
	-> n	-- ^ The Result.
evaluate x	= foldr ((+) . raise) 0 . getMonomialList	where
	powers	= iterate (* x) 1

	raise monomial
		| exponent' < 0	= error $ "Factory.Data.Polynomial.evaluate.raise:\tnegative exponent; " ++ show exponent'
		| otherwise	= Data.Monomial.getCoefficient monomial * Data.List.genericIndex powers exponent'
		where
			exponent'	= Data.Monomial.getExponent monomial

-- | Convert the type of the /coefficient/s.
realCoefficientsToFrac :: (Real r, Fractional f) => Polynomial r e -> Polynomial f e
realCoefficientsToFrac	= lift (Data.Monomial.realCoefficientToFrac `map`)