Eq-1.0: EqManips/Polynome.hs
{-# OPTIONS_GHC -fno-warn-orphans #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE Rank2Types #-}
module EqManips.Polynome( convertToPolynome
, convertToFormula
, polynomizeFormula
, polyMap
, polyCoeffMap
, scalarToCoeff
, coefToFormula
, isCoeffNull
, prepareFormula
, syntheticDiv
, polyAsFormula
-- | Pack/simplify polynome with only one coefficient
-- and/or null coef.
, simplifyPolynome
) where
import Data.Maybe( fromMaybe )
import Data.Ord( comparing )
import Control.Applicative( (<$>), (<*>) )
import Control.Arrow( (***), second )
import Control.Monad( join )
import Data.Either( partitionEithers )
import Data.List( sortBy, groupBy, foldl' )
import Data.Ratio
import EqManips.Types
import EqManips.Algorithm.Utils
import EqManips.FormulaIterator
import qualified EqManips.ErrorMessages as Err
-- | will pack/simplify internal representation of a polynome.
-- If there is only one null coefficient only subPoly will be present
simplifyPolynome :: Polynome -> Polynome
simplifyPolynome (Polynome v p@[(lastCoeff, PolyRest constant)])
| isCoeffNull lastCoeff = PolyRest constant
| otherwise = Polynome v p
simplifyPolynome (Polynome v p@[(lastCoeff, subPoly)])
| isCoeffNull lastCoeff = subPoly
| otherwise = Polynome v p
simplifyPolynome a = a
polyAsFormula :: Polynome -> FormulaPrim
polyAsFormula (PolyRest coeff) = coefToFormula coeff
polyAsFormula (Polynome _ [(0, a)]) = polyAsFormula a
polyAsFormula p = poly p
-- | Given a formula, it'll try to convert it to a polynome.
-- Formula should be expanded and in list form to get this
-- function to work (nested shit shouldn't work)
convertToPolynome :: Formula ListForm -> Maybe Polynome
convertToPolynome (Formula f) = polynomize
$ prepareFormula f
convertToFormula :: Polynome -> Formula ListForm
convertToFormula = Formula . convertToFormulaPrim
-- | Run across the whole formula and replace what
-- can polynomized by a polynome
polynomizeFormula :: Formula ListForm -> Formula ListForm
polynomizeFormula (Formula f) = Formula $ topDownTraversal converter f
where converter f' = poly <$> convertToPolynome (Formula f')
-- | Convert a polynome into a simpler formula using only
-- basic operators.
convertToFormulaPrim :: Polynome -> FormulaPrim
convertToFormulaPrim (PolyRest coeff) = coefToFormula coeff
convertToFormulaPrim (Polynome var lst) =
foldl' constructor realFirst rest
where constructor a (Left b) = a + b
constructor a (Right b) = a - b
realFirst = either id id felem
(felem : rest) = map elemConverter lst
fvar = Variable var
elemConverter (degree,def) =
degreeOf (convertToFormulaPrim def)
(coefToFormula degree)
degreeOf fdef (CInteger 0)
| isConstantNegative fdef = Right $ negateConstant fdef
| otherwise = Left $ fdef
degreeOf (CInteger 1 ) (CInteger 1) = Left fvar
degreeOf (CInteger (-1)) (CInteger 1) = Right fvar
degreeOf fdef (CInteger 1)
| isConstantNegative fdef = Right $ negateConstant fdef * fvar
| otherwise = Left $ fdef * fvar
degreeOf (CInteger 1) deg = Left $ fvar ** deg
degreeOf (CInteger (-1)) deg = Right $ fvar ** deg
degreeOf fdef deg
| isConstantNegative fdef =
Right $ negateConstant fdef * (fvar ** deg)
| otherwise = Left $ fdef * (fvar ** deg)
-- | Conversion from coef to basic formula. ratio
-- are converted to (a/b), like a division.
coefToFormula :: PolyCoeff -> FormulaPrim
coefToFormula (CoeffFloat f) = CFloat f
coefToFormula (CoeffInt i) = CInteger i
coefToFormula (CoeffRatio r) = if denominator r == 1
then CInteger $ numerator r
else Fraction r
-- | Flatten the formula, remove all the OpSub and replace them
-- by OpAdd. Also bring lowest variables to the front, regardless of
-- their order. Ordering is very important in this function. All
-- the polynome construction is built uppon the ordering created here.
prepareFormula :: FormulaPrim -> FormulaPrim
prepareFormula = polySort . formulaFlatter
polySort :: FormulaPrim -> FormulaPrim
polySort = depthFormulaPrimTraversal `asAMonad` sortBinOp sorter
where lexicalOrder EQ b = b
lexicalOrder a _ = a
invert LT = GT
invert EQ = EQ
invert GT = LT
-- Special sort which bring x in front, followed by others. Lexical
-- order first.
sorter (Poly _ p1) (Poly _ p2) = compare p1 p2
sorter (Poly _ _) _ = LT
sorter _ (Poly _ _) = GT
-- Rules to fine-sort '*' elements
-- (x before y), no regard for formula degree
sorter (Variable v1) (Variable v2) = compare v1 v2
-- x ^ n * y ^ n (n can be one, not shown)
sorter (BinOp _ OpPow [Variable v1, p1])
(BinOp _ OpPow [Variable v2, p2]) =
compare v1 v2 `lexicalOrder` compare p1 p2
-- x * y ^ n
sorter (Variable v1)
(BinOp _ OpPow (Variable v2:_)) =
compare v1 v2 `lexicalOrder` LT
-- x ^ n * y
sorter (BinOp _ OpPow (Variable v1:_))
(Variable v2) = compare v1 v2 `lexicalOrder` GT
-- (x * ...) + y ^ n
sorter (BinOp _ OpMul (Variable v1:_))
(BinOp _ OpPow [Variable v2, _]) = compare v1 v2 `lexicalOrder` LT
-- x ^ n + (y * ...)
sorter (BinOp _ OpPow [Variable v1, _])
(BinOp _ OpMul (Variable v2:_)) = compare v1 v2 `lexicalOrder` GT
-- (x ^ m * ...) + y ^ n
sorter (BinOp _ OpMul (BinOp _ OpPow [Variable v1,p1]:_))
(BinOp _ OpPow [Variable v2, p2]) =
compare v1 v2 `lexicalOrder` compare p1 p2
-- x ^ n + (y ^ m * ...)
sorter (BinOp _ OpPow [Variable v1, p1])
(BinOp _ OpMul (BinOp _ OpPow [Variable v2,p2]:_)) =
compare v1 v2 `lexicalOrder` compare p1 p2
-- Rules to fine sort the '+' elements, lowest variable
-- first (x before y), smallest order first (x before x ^ 15)
-- (x^n * ....) + (y^n * ...)
sorter (BinOp _ OpMul (BinOp _ OpPow (Variable v1: power1):_))
(BinOp _ OpMul (BinOp _ OpPow (Variable v2: power2):_)) =
compare v1 v2 `lexicalOrder` compare power1 power2
-- (x * ...) + (y^n * ...)
sorter (BinOp _ OpMul (Variable v1:_))
(BinOp _ OpMul (BinOp _ OpPow (Variable v2:_):_)) =
compare v1 v2 `lexicalOrder` LT
-- (x^n * ...) + (y * ...)
sorter (BinOp _ OpMul (BinOp _ OpPow (Variable v1:_):_))
(BinOp _ OpMul (Variable v2:_)) = compare v1 v2 `lexicalOrder` GT
-- (x * ...) + (y * ...)
sorter (BinOp _ OpMul (Variable v1:_))
(BinOp _ OpMul (Variable v2:_)) = compare v1 v2
-- x + (y * ...)
sorter (Variable v1)
(BinOp _ OpMul (Variable v2:_)) = compare v1 v2
-- (x * ...) + y
sorter (BinOp _ OpMul (Variable v1:_))
(Variable v2) = compare v1 v2
sorter (BinOp _ OpPow a) (BinOp _ OpPow b) =
case comparing length a b of
LT -> LT
GT -> GT
EQ -> foldl' (\acc (a', b') -> if acc == EQ
then acc
else compare a' b') EQ $ zip a b
-- x ^ n * ?
sorter _ (BinOp _ OpPow (Variable _:_)) = GT
sorter (BinOp _ OpPow (Variable _:_)) _ = LT
-- make sure weird things go at the end.
sorter (Variable _) _ = LT
sorter _ (Variable _) = GT
-- Just reverse the general readable order.
sorter a b = invert $ compare a b
-- | Called when we found an OpSub operator within the
-- formula. -- We assume that the formula as been previously sorted
resign :: FormulaPrim -> [FormulaPrim] -> [FormulaPrim]
resign = globalResign
where globalResign (BinOp _ OpMul (a:xs)) acc
| isFormulaInteger a = case atomicResign a of
Nothing -> binOp OpMul (CInteger (-1):a:xs) : acc
Just a' -> binOp OpMul (a':xs) : acc
globalResign (BinOp _ OpAdd lst) acc = foldr resign acc lst
globalResign a acc = fromMaybe (CInteger (-1) * a) (atomicResign a) : acc
atomicResign (CInteger i) = Just $ CInteger (-i)
atomicResign (CFloat i) = Just $ CFloat (-i)
atomicResign (UnOp _ OpNegate a) = Just a
atomicResign (BinOp _ OpDiv [a,b]) = (\a' -> binOp OpDiv [a', b]) <$> atomicResign a
atomicResign _ = Nothing
-- | Flatten a whole formula, by flattening from the leafs.
formulaFlatter :: FormulaPrim -> FormulaPrim
formulaFlatter = depthFormulaPrimTraversal `asAMonad` listFlatter
-- | Given a formula in LIST form, provide a version
-- with only Pluses.
listFlatter :: FormulaPrim -> FormulaPrim
listFlatter (BinOp _ OpAdd lst) = binOp OpAdd $ foldr flatter [] lst
where flatter (BinOp _ OpSub (x:xs)) acc = x : foldr resign acc xs
flatter (BinOp _ OpAdd lst') acc = lst' ++ acc
flatter x acc = x:acc
listFlatter (BinOp _ OpSub ((BinOp _ OpAdd lst'):xs)) =
binOp OpAdd $ lst' ++ foldr resign [] xs
listFlatter (BinOp _ OpSub (x:xs)) =
binOp OpAdd $ x : foldr resign [] xs
-- Remove the maximum of negation in the multiplication.
-- In the end, keep the needed negation into the first term
listFlatter (BinOp _ OpMul lst) = if foldr countInversion False lst
then let (x:xs) = map cleanSign lst
in binOp OpMul $ resign x xs
else binOp OpMul $ map cleanSign lst
where iodd :: Int -> Bool
iodd = odd
countInversion whole@(UnOp _ OpNegate _) acc =
if iodd . fst $ getUnsignedRoot 0 whole
then not acc
else acc
countInversion _ acc = acc
getUnsignedRoot n (UnOp _ OpNegate something) = getUnsignedRoot (n+1) something
getUnsignedRoot n (something) = (n :: Int, something)
cleanSign whole@(UnOp _ OpNegate _) = snd $ getUnsignedRoot 0 whole
cleanSign a = a
listFlatter a = a
-- | Verify if the coefficient is valid in the context
-- of polynomial. might add a reduction rule here.
evalCoeff :: [FormulaPrim] -> Maybe PolyCoeff
evalCoeff [CInteger i] = Just $ CoeffInt i
evalCoeff [CFloat f] = Just $ CoeffFloat f
evalCoeff [UnOp _ OpNegate (CInteger i)] = Just $ CoeffInt (-i)
evalCoeff [UnOp _ OpNegate (CFloat f)] = Just $ CoeffFloat (-f)
evalCoeff [BinOp _ OpDiv [CInteger a, CInteger b]] = Just . CoeffRatio $ a % b
evalCoeff [UnOp _ OpNegate (BinOp _ OpDiv [CInteger a, CInteger b])] = Just . CoeffRatio $ (-a) % b
evalCoeff _ = Nothing
-- | Given a rest (a leading +c, where c is a constant) and
-- a group of variable and coefficients, try to build a full
-- blown polynomial out of it.
translator :: [FormulaPrim] -- Unnammed rest (var ^ 0)
-> [(String, [(FormulaPrim, FormulaPrim)])] -- Named things x ^ n or y ^ n, n > 0
-> Maybe (Maybe Polynome) -- ^ First maybe: error, nested maybe: empty
translator [] [(var, coefs)] = do
result <- mapM (\(rank, polyn) -> (,) <$> evalCoeff [rank] <*> polynomize polyn) coefs
return . Just $ Polynome var result
translator pow0 [(var, coefs)] = do
result <- mapM (\(rank,polyn) -> (,) <$> evalCoeff [rank] <*> polynomize polyn) coefs
rest <- evalCoeff pow0
return . Just . Polynome var $ (CoeffInt 0, PolyRest rest):result
translator pow0 ((var,coefs):rest) = do
result <- mapM (\ (rank,polyn) -> (,) <$> evalCoeff [rank] <*> polynomize polyn) coefs
subPolynome <- translator pow0 rest
let finalList = case subPolynome of
Nothing -> result
Just p -> (CoeffInt 0, p) : result
return . Just $ Polynome var finalList
translator pow0 [] = return $ PolyRest <$> evalCoeff pow0
-- | Try to transform a formula in polynome.
polynomize :: FormulaPrim -> Maybe Polynome
polynomize wholeFormula@(BinOp _ OpMul _) = polynomize (binOp OpAdd [wholeFormula])
-- HMmm?
polynomize (BinOp _ OpAdd lst) = join -- flatten a maybe level, we don't distingate
. translator pow0 -- cases at the upper level.
. packCoefs
$ varGroup polys
where (polys, pow0) = partitionEithers $ map extractFirstTerm lst
varGroup = groupBy (\(var,_,_) (var',_,_) -> var == var')
coeffGroup = groupBy (\(_,coeff1,_) (_,coeff2,_) -> coeff1 == coeff2)
packCoefs :: [[(String,FormulaPrim,FormulaPrim)]] -> [(String, [(FormulaPrim,FormulaPrim)])]
packCoefs varGrouped = map grouper varGrouped
where nameOfGroup ((varName, _,_):_) = varName
nameOfGroup [] = error Err.polynom_emptyCoeffPack
grouper :: [(String,FormulaPrim,FormulaPrim)] -> (String, [(FormulaPrim,FormulaPrim)])
grouper lst' = (nameOfGroup lst'
, [(coef group, polySort $ binOp OpAdd $ defs group)
| group <- coeffGroup lst'])
defs = map (\(_,_,def) -> def)
coef ((_,c1,_):_) = c1
coef [] = error Err.polynom_emptyCoeffPack
polynomize (BinOp _ OpPow [Variable v, CInteger c]) =
Just $ Polynome v [(CoeffInt c, PolyRest 1)]
polynomize _ = Nothing
-- | Function in charge of extracting variable name (if any), and
-- return the coeff function.
extractFirstTerm :: FormulaPrim
-> Either (String, FormulaPrim, FormulaPrim) FormulaPrim
extractFirstTerm fullFormula@(BinOp _ OpMul lst) = varCoef lst
where varCoef ((BinOp _ OpPow [(Variable v), f]):xs)
| isFormulaConstant f = Left (v, f, multify xs)
varCoef ((Variable v):xs) = Left (v, CInteger 1, multify xs)
varCoef _ = Right fullFormula
multify [] = error $ Err.empty_binop "Polynome.OpMul"
multify [x] = x
multify alist = binOp OpMul alist
extractFirstTerm (BinOp _ OpPow [Variable v, order])
| isFormulaConstant order = Left (v, order, CInteger 1)
extractFirstTerm (Variable v) = Left (v, CInteger 1, CInteger 1)
extractFirstTerm a = Right a
--------------------------------------------------
---- Polynome instances
--------------------------------------------------
-- | Only to map on the polynome coefficients (not the degree
-- of it).
polyCoeffMap :: (PolyCoeff -> PolyCoeff) -> Polynome -> Polynome
polyCoeffMap f = polyMap mapper
where mapper (deg, PolyRest c) = (deg, PolyRest $ f c)
mapper otherCoeff = otherCoeff
-- | polynome mapping
polyMap :: ((PolyCoeff, Polynome) -> (PolyCoeff, Polynome)) -> Polynome -> Polynome
polyMap f (Polynome s lst) = Polynome s $ map (second $ polyMap f) lst
polyMap f rest@(PolyRest _) = snd $ f (CoeffInt 0, rest)
-- | Transform a scalar formula component to
-- a polynome coefficient. If formula is not
-- a scalar, error is called.
scalarToCoeff :: FormulaPrim -> PolyCoeff
scalarToCoeff (UnOp _ OpNegate f) = negate $ scalarToCoeff f
scalarToCoeff (CFloat f) = CoeffFloat f
scalarToCoeff (CInteger i) = CoeffInt i
scalarToCoeff (BinOp _ OpDiv [CInteger a, CInteger b]) = CoeffRatio $ a % b
scalarToCoeff _ = error Err.polynom_coeff_notascalar
-- | Operation on polynome coefficients. Put there
-- to provide automatic Equality derivation for polynome
-- and in the end... Formula
coeffOp :: (forall a. (Num a) => a -> a -> a)
-> PolyCoeff -> PolyCoeff -> PolyCoeff
coeffOp op c1 c2 = eval $ polyCoeffCast c1 c2
where eval (CoeffInt i1, CoeffInt i2) = CoeffInt $ i1 `op` i2
eval (CoeffFloat f1, CoeffFloat f2) = CoeffFloat $ f1 `op` f2
eval (CoeffRatio r1, CoeffRatio r2) = CoeffRatio $ r1 `op` r2
eval _ = error Err.polynom_bad_casting
inf :: PolyCoeff -> PolyCoeff -> Bool
inf = coeffPredicate ((<) :: forall a. (Ord a) => a -> a -> Bool)
-- | Implement the same idea that the one used by the
-- mergesort, only this time it's only used to perform
-- addition or substraction on polynomial.
lockStep :: (Polynome -> Polynome -> Polynome)
-> [(PolyCoeff, Polynome)] -> [(PolyCoeff, Polynome)]
-> [(PolyCoeff, Polynome)]
lockStep op xs [] = map (\(c,v) -> (c, v `op` PolyRest 0)) xs
lockStep op [] ys = map (\(c,v) -> (c, PolyRest 0 `op` v)) ys
lockStep op whole1@((c1, def1):xs) whole2@((c2, def2):ys)
| c1 `inf` c2 =
(c1, def1 `op` PolyRest (CoeffInt 0)) : lockStep op xs whole2
| c1 == c2 =
(c1, def1 `op` def2) : lockStep op xs ys
| otherwise =
(c2, PolyRest (CoeffInt 0) `op` def2) : lockStep op whole1 ys
-- | Tell if a coefficient can be treated as Null
isCoeffNull :: PolyCoeff -> Bool
isCoeffNull (CoeffInt 0) = True
isCoeffNull (CoeffFloat 0.0) = True
isCoeffNull (CoeffRatio r) = numerator r == 0
isCoeffNull _ = False
coeffPropagator :: (forall a. (Num a) => a -> a -> a) -> (PolyCoeff, Polynome) -> (PolyCoeff, Polynome)
coeffPropagator op (degree, PolyRest a) = (degree, PolyRest $ coeffOp op (CoeffInt 0) a)
coeffPropagator op (degree, Polynome v lst) = (degree, Polynome v $ map (coeffPropagator op) lst)
polySimpleOp :: (forall a. (Num a) => a -> a -> a) -> Polynome -> Polynome -> Polynome
polySimpleOp _ (Polynome _ []) _ = error Err.ill_formed_polynomial
polySimpleOp _ _ (Polynome _ []) = error Err.ill_formed_polynomial
polySimpleOp op (PolyRest c1) (PolyRest c2) = PolyRest $ coeffOp op c1 c2
polySimpleOp op left@(PolyRest c1) (Polynome v1 as@((coeff, def):xs))
| isCoeffNull coeff = case def of
PolyRest a -> Polynome v1 $ (CoeffInt 0, PolyRest $ coeffOp op c1 a) : map (coeffPropagator op) xs
_ -> Polynome v1 $ (coeff,polySimpleOp op left def) : map (coeffPropagator op) xs
| otherwise =
Polynome v1 $ (CoeffInt 0, PolyRest $ coeffOp op c1 (CoeffInt 0)) : map (coeffPropagator op) as
polySimpleOp op (Polynome v1 as@((coeff, def):xs)) right@(PolyRest c1)
| isCoeffNull coeff = case def of
PolyRest a -> Polynome v1 $ (CoeffInt 0, PolyRest $ coeffOp op a c1)
: map (coeffPropagator $ flip op) xs
_ -> Polynome v1 $ (coeff,polySimpleOp op def right)
: map (coeffPropagator $ flip op) xs
| otherwise =
Polynome v1 $ (CoeffInt 0, PolyRest $ coeffOp op (CoeffInt 0) c1)
: as
polySimpleOp op (Polynome v1 as@((c, d1):rest)) right@(Polynome v2 bs)
| v1 > v2 = polySimpleOp (flip op) (Polynome v2 bs) (Polynome v1 as)
| v1 == v2 =
let computedCoefs = lockStep op as bs
in if null computedCoefs then PolyRest 0
else Polynome v1 computedCoefs
| isCoeffNull c =
Polynome v1 $ (c, polySimpleOp op d1 right) : map (coeffPropagator $ flip op) rest
| otherwise =
Polynome v1 $ (CoeffInt 0, polySimpleOp op (PolyRest $ CoeffInt 0) right)
: map (coeffPropagator $ flip op) as
-- | Multiply two polynomials between them using the brute force
-- way, algorithm in O(n²)
polyMul :: Polynome -> Polynome -> Polynome
polyMul p@(Polynome _ _) (PolyRest c) = polyCoeffMap (* c) p
polyMul (PolyRest c) p@(Polynome _ _) = polyCoeffMap (c *) p
polyMul (PolyRest c) (PolyRest c2) = PolyRest $ coeffOp (*) c c2
polyMul p1@(Polynome v1 _) p2@(Polynome v2 _) | v1 > v2 = polyMul p2 p1
polyMul (Polynome v1 coefs1) p2@(Polynome v2 coefs2)
| v1 /= v2 {- v1 < v2 by previous line -} =
Polynome v1 $ map (\(order, c) -> (order, polyMul c p2)) coefs1
| otherwise {- v1 == v2 -} =
Polynome v1
{-. map (\lst@((o,_):_) -> (o, foldr1 (+) $ map snd lst))-}
. map (\lst@((o,_):_) -> (o, sum $ map snd lst))
. groupBy (\(o1,_) (o2,_) -> o1 == o2) -- Regroup same order together
$ sortBy (\(c1,_) (c2,_) -> compare c1 c2)
[ (degree1 + degree2, c1 * c2) | (degree1, c1) <- coefs1, (degree2, c2) <- coefs2]
--------------------------------------------------
---- Division
--------------------------------------------------
-- | Expand coefficients of an _UNIVARIATE_ polynomial
-- in an descending way, each integer power given a
-- coefficient (0 if none).
expandCoeff :: Polynome -> Maybe [PolyCoeff]
expandCoeff (PolyRest _) = error ""
expandCoeff (Polynome _ coefs) = snd <$> foldl' sparser (Just (-1, [])) coefs
where sparser (Just (lastNum, lst)) (CoeffInt n, PolyRest r) =
Just (fromInteger n, r : replicate (fromInteger n - lastNum - 1) (CoeffInt 0)
++ lst)
sparser _ _ = Nothing
-- | Tell if a polynomial has only one var
isPolyMonovariate :: Polynome -> Bool
isPolyMonovariate (PolyRest _) = False
isPolyMonovariate (Polynome _ coefs) = all isCoeff coefs
where isCoeff (_,PolyRest _) = True
isCoeff _ = False
-- | Given a power descending list of coefficient, rearrange
-- them to make it normal polynomial
packCoeffs :: [PolyCoeff] -> [(PolyCoeff, Polynome)]
packCoeffs = reverse . snd . foldr packer (0, [])
where packer coeff (n, lst)
| isCoeffNull coeff = (n + 1, lst)
| otherwise = (n + 1, (CoeffInt n, PolyRest coeff) : lst)
-- | Apply an operation on an head of a list given an other list.
-- return Nothing if first list finish after "applied" list.
headApply :: (a -> b -> a) -> [a] -> [b] -> Maybe [a]
headApply _ [] [] = Just []
headApply _ rest [] = Just rest
headApply _ [] _ = Nothing
headApply f (x:xs) (y:ys) = (f x y :) <$> headApply f xs ys
-- | Try to perform a polynomial synthetic division on
-- monovariate polynomial.
syntheticDiv :: Polynome -> Polynome -> (Maybe Polynome, Maybe Polynome)
syntheticDiv polyn@(Polynome var lst1) divisor@(Polynome var' lst2)
| var == var'
&& isPolyMonovariate polyn && isPolyMonovariate divisor
&& fst (last lst1) > fst (last lst2) =
(finalize . packCoeffs . map (/ normalizingCoeff)
*** finalize . packCoeffs)
. splitAt (length coefList + 1 - length divCoeff)
$ firstCoeff : syntheticInnerDiv divCoeff firstCoeff coefList
where Just (firstCoeff: coefList) = expandCoeff polyn
Just (firstDivCoeff:divCoeff) = map negate <$> expandCoeff divisor
normalizingCoeff = negate firstDivCoeff
finalize [] = Nothing
finalize lst = Just $ Polynome var lst
syntheticInnerDiv :: [PolyCoeff]
-> PolyCoeff -> [PolyCoeff] -> [PolyCoeff]
syntheticInnerDiv _ _ [] = []
syntheticInnerDiv diviCoeff prevCoeff polyCoeff =
case endCoeffs of
Just [] -> error "syntheticDiv - empty rest, impossible"
Just (x:xs) -> x : syntheticInnerDiv diviCoeff x xs
Nothing -> polyCoeff
where normalizedCoeff = prevCoeff / normalizingCoeff
endCoeffs = headApply (+) polyCoeff
$ map (normalizedCoeff *) diviCoeff
syntheticDiv _ _ = (Nothing, Nothing)
instance Num PolyCoeff where
fromInteger = CoeffInt
(+) = coeffOp (+)
(-) = coeffOp (-)
(*) = coeffOp (*)
abs (CoeffInt i) = CoeffInt $ abs i
abs (CoeffFloat f) = CoeffFloat $ abs f
abs (CoeffRatio r) = CoeffRatio $ abs r
signum (CoeffInt i) = CoeffInt $ signum i
signum (CoeffFloat f) = CoeffFloat $ signum f
signum (CoeffRatio r) = CoeffRatio $ signum r
instance Fractional PolyCoeff where
a / b = case polyCoeffCast a b of
(CoeffInt i1, CoeffInt i2) -> if i1 `mod` i2 == 0
then CoeffInt $ i1 `div` i2
else CoeffRatio $ i1 % i2
(CoeffFloat f1, CoeffFloat f2) -> CoeffFloat $ f1 / f2
(CoeffRatio r1, CoeffRatio r2) -> CoeffRatio $ r1 / r2
_ -> error Err.polynom_bad_casting
recip (CoeffFloat f) = CoeffFloat $ recip f
recip (CoeffInt i) = CoeffRatio $ 1 % i
recip (CoeffRatio r) = if denominator r' == 1
then CoeffInt $ numerator r'
else CoeffRatio r'
where r' = recip r
fromRational = CoeffRatio
instance Num Polynome where
(+) = polySimpleOp (+)
(-) = polySimpleOp (-)
(*) = polyMul
fromInteger = PolyRest . fromInteger
abs = error "Unimplemented-Abs"
signum = error "Unimplemented-signum"