ideas-0.7: src/Domain/Math/Polynomial/Rules.hs
-----------------------------------------------------------------------------
-- Copyright 2010, Open Universiteit Nederland. This file is distributed
-- under the terms of the GNU General Public License. For more information,
-- see the file "LICENSE.txt", which is included in the distribution.
-----------------------------------------------------------------------------
-- |
-- Maintainer : bastiaan.heeren@ou.nl
-- Stability : provisional
-- Portability : portable (depends on ghc)
--
-----------------------------------------------------------------------------
module Domain.Math.Polynomial.Rules
( sameConFactor, abcFormula, allPowerFactors, bringAToOne, cancelTerms
, commonFactorVar, commonFactorVarNew, defPowerNat, distributeDivision
, distributeTimes, distributionSquare, exposeSameFactor, factorLeftAsSquare
, factorVariablePower, flipEquation, higherSubst, merge, moveToLeft, mulZero
, niceFactors, niceFactorsNew, noDivisionConstant, noLinFormula, oneVar
, parentNotNegCheck, prepareSplitSquare, quadraticRuleOrder, removeDivision
, ruleApproximate, ruleNormalizeMixedFraction, ruleNormalizeRational
, sameFactor, simplerLinearFactor, simplerPolynomial, simplerSquareRoot
, squareBothSides, substBackVar, varToLeft
) where
import Common.Library hiding (terms, simplify)
import Common.Uniplate (universe, descend)
import Common.Utils
import Control.Monad
import Data.List
import Data.Maybe
import Data.Ord
import Data.Ratio
import Domain.Math.Approximation (precision)
import Domain.Math.Clipboard
import Domain.Math.Data.OrList
import Domain.Math.Data.Polynomial
import Domain.Math.Data.Relation
import Domain.Math.Equation.BalanceRules
import Domain.Math.Equation.CoverUpRules
import Domain.Math.Expr
import Domain.Math.Numeric.Views
import Domain.Math.Polynomial.CleanUp
import Domain.Math.Polynomial.Views
import Domain.Math.Power.OldViews
import Domain.Math.Simplification hiding (simplifyWith)
import Domain.Math.SquareRoot.Views
import Prelude hiding ( (^) )
quadraticRuleOrder :: [Id]
quadraticRuleOrder =
[ getId coverUpTimes, getId (coverUpMinusRightWith oneVar)
, getId (coverUpMinusLeftWith oneVar), getId (coverUpPlusWith oneVar)
, getId coverUpPower
, getId commonFactorVar, getId simplerPolynomial
, getId niceFactors, getId noLinFormula
, getId cancelTerms, getId sameConFactor, getId distributionSquare
, getId allPowerFactors
]
lineq, quadreq, polyeq :: String
lineq = "algebra.equations.linear"
quadreq = "algebra.equations.quadratic"
polyeq = "algebra.equations.polynomial"
------------------------------------------------------------
-- General form rules: ax^2 + bx + c = 0
quadraticNF :: View Expr (String, (Rational, Rational, Rational))
quadraticNF = polyNormalForm rationalView >>> second quadraticPolyView
-- ax^2 + bx = 0
commonFactorVar :: Rule (Equation Expr)
commonFactorVar = rhsIsZero commonFactorVarNew
-- Maybe to be replaced by more general factorVariablePower??
commonFactorVarNew :: Rule Expr
commonFactorVarNew = describe "Common factor variable" $
makeSimpleRule (quadreq, "common-factor") $ \expr -> do
(x, (a, b, c)) <- match quadraticNF expr
guard (b /= 0 && c == 0)
-- also search for constant factor
let d | a<0 && b<0 = -gcdFrac a b
| otherwise = gcdFrac a b
return (fromRational d .*. Var x .*. (fromRational (a/d) .*. Var x .+. fromRational (b/d)))
gcdFrac :: Rational -> Rational -> Rational
gcdFrac r1 r2 =
if denominator r1 == 1 && denominator r2 == 1
then fromInteger (numerator r1 `gcd` numerator r2)
else 1
-- ax^2 + c = 0
noLinFormula :: Rule (Equation Expr)
noLinFormula = describe "No linear term ('b=0')" $ liftRule myView $
makeSimpleRule (quadreq, "no-lin") $ \((x, (a, b, c)), rhs) -> do
guard (rhs == 0 && b == 0 && c /= 0)
return $ if a>0 then ((x, (a, 0, 0)), -c)
else ((x, (-a, 0, 0)), c)
where
myView = constantRight quadraticNF
-- search for (X+A)*(X+B) decomposition
niceFactors :: Rule (Equation Expr)
niceFactors = rhsIsZero niceFactorsNew
-- search for (X+A)*(X+B) decomposition
niceFactorsNew :: Rule Expr
niceFactorsNew = describe "Find a nice decomposition" $
makeSimpleRuleList (quadreq, "nice-factors") $ \expr -> do
let sign t@(x, (a, b, c)) = if a== -1 then (x, (1, -b, -c)) else t
(x, (a, b, c)) <- liftM sign (matchM (polyNormalForm integerView >>> second quadraticPolyView) expr)
guard (a==1)
let ok (i, j) = i+j == b
f (i, j)
| i == j = -- special case
(Var x + fromInteger i) ^ 2
| otherwise =
(Var x + fromInteger i) * (Var x + fromInteger j)
map f (filter ok (factors c))
where
factors :: Integer -> [(Integer, Integer)]
factors n = [ pair
| let h = (floor :: Double -> Integer) (sqrt (abs (fromIntegral n)))
, a <- [1..h], let b = n `div` a, a*b == n
, pair <- [(a, b), (negate a, negate b)]
]
-- Simplify polynomial by multiplying (or dividing) the terms:
-- 1) If a,b,c are ints, then find gcd
-- 2) If any of a,b,c is a fraction, find lcm of denominators
-- 3) If a<0, then also suggest to change sign (return two solutions)
simplerPolynomial :: Rule (Equation Expr)
simplerPolynomial = describe "simpler polynomial" $
rhsIsZero $ liftRuleIn thisView $
makeSimpleRuleList (quadreq, "simpler-poly") $ \(a, b, c) -> do
r <- findFactor (filter (/=0) [a, b, c])
d <- if a >= 0 then [r] else [-r, r]
guard (d `notElem` [0, 1])
return (a*d, b*d, c*d)
where
thisView = polyNormalForm rationalView >>> swapView >>> first quadraticPolyView
-- Simplified variant of simplerPoly: just bring a to 1.
-- Needed for quadratic strategy without square formula
bringAToOne :: Rule (Equation Expr)
bringAToOne = rhsIsZero $ liftRuleIn thisView $
describe "Bring 'a' to one" $
makeSimpleRule (quadreq, "scale") $ \(a, b, c) -> do
guard (a `notElem` [0, 1])
return (1, b/a, c/a)
where
thisView = polyNormalForm rationalView >>> swapView >>> first quadraticPolyView
------------------------------------------------------------
-- General form rules: expr = 0
-- Rule must be symmetric in side of equation
mulZero :: Rule (OrList (Equation Expr))
mulZero = describe "multiplication is zero" $
makeSimpleRuleList (quadreq, "product-zero") $ oneDisjunct bothSides
where
bothSides eq = oneSide eq `mplus` oneSide (flipSides eq)
oneSide (lhs :==: rhs) = do
guard (rhs == 0)
(_, xs) <- matchM productView lhs
guard (length xs > 1)
return $ orList $ flip map xs $ \e ->
case match (polyNormalForm rationalView >>> second linearPolyView) e of
-- special cases (simplify immediately, as in G&R)
Just (x, (a, b))
| a == 1 ->
Var x :==: fromRational (-b)
| a == -1 ->
Var x :==: fromRational b
_ -> e :==: 0
------------------------------------------------------------
-- Constant form rules: expr = constant
-- Use this configuration for covering-up plus and minus symbols!
-- Prevent (x^2+3x)+5 = 0 to be covered up
oneVar :: ConfigCoverUp
oneVar = configCoverUp
{ configName = "onevar"
, predicateCovered = \a -> p1 a || p2 a
, predicateCombined = hasNoVar
, coverLHS = True
, coverRHS = True
}
where
p1 = (==1) . length . vars
-- predicate p2 tests for cases such as 12*(x^2-3*x)+8 == 56
p2 a = fromMaybe False $ do
(x, y) <- match timesView a
return (hasSomeVar x /= hasSomeVar y)
------------------------------------------------------------
-- Top form rules: expr1 = expr2
-- Do not simplify (5+sqrt 53)/2
simplerSquareRoot :: Rule Expr
simplerSquareRoot = describe "simpler square root" $
makeSimpleRule (quadreq, "simpler-sqrt") $ \e -> do
xs <- f e
guard (not (null xs))
new <- canonical (squareRootViewWith rationalView) e
ys <- f new
guard (xs /= ys)
return new
where
-- return numbers under sqrt symbol
f :: Expr -> Maybe [Rational]
f e = liftM sort $ sequence [ match rationalView a | Sqrt a <- universe e ]
cancelTerms :: Rule (Equation Expr)
cancelTerms = describe "Cancel terms" $
makeSimpleRule (quadreq, "cancel") $ \(lhs :==: rhs) -> do
xs <- match sumView lhs
ys <- match sumView rhs
let zs = filter (`elem` ys) (nub xs)
guard (not (null zs))
let without as = build sumView (as \\ zs)
return (without xs :==: without ys)
-- Two out of three "merkwaardige producten"
distributionSquare :: Rule Expr
distributionSquare = describe "distribution square" $
ruleList (quadreq, "distr-square")
[ \a b -> (a+b)^2 :~> a^2 + 2*a*b + b^2
, \a b -> (a-b)^2 :~> a^2 - 2*a*b + b^2
]
-- a^2 == b^2
squareBothSides :: Rule (OrList (Equation Expr))
squareBothSides = describe "square both sides" $
rule (quadreq, "square-both") $ \a b ->
orList [a^2 :==: b^2] :~> orList [a :==: b, a :==: -b]
-- prepare splitting a square; turn lhs into x^2+bx+c such that (b/2)^2 is c
prepareSplitSquare :: Rule (Equation Expr)
prepareSplitSquare = describe "prepare split square" $
liftRule myView $
makeSimpleRule (quadreq, "prepare-split") $ \((x, (a, b, c)), r) -> do
let newC = (b/2)*(b/2)
newRHS = r + newC - c
guard (a==1 && b/=0 && c /= newC)
return ((x, (a, b, newC)), newRHS)
where
myView = constantRight quadraticNF
-- factor left-hand side into (ax + c)^2
factorLeftAsSquare :: Rule (Equation Expr)
factorLeftAsSquare = describe "factor left as square" $
makeSimpleRule (quadreq, "left-square") $ \(lhs :==: rhs) -> do
guard (hasNoVar rhs)
(x, (a, b, c)) <- match quadraticNF lhs
let h = b/2
guard (a==1 && b/=0 && h*h == c)
return ((Var x + build rationalView h)^2 :==: rhs)
-- flip the two sides of an equation
flipEquation :: Rule (Equation Expr)
flipEquation = doBeforeTrans condition $
describe "flip equation" $
rule (lineq, "flip") $ \a b ->
(a :==: b) :~> (b :==: a)
where
condition = makeTrans $ \eq@(lhs :==: rhs) -> do
guard (hasSomeVar rhs && hasNoVar lhs)
return eq
-- Afterwards, merge and sort
moveToLeft :: Rule (Equation Expr)
moveToLeft = describe "Move to left" $
makeSimpleRule (quadreq, "move-left") $ \(lhs :==: rhs) -> do
guard (rhs /= 0 && hasSomeVar lhs && (hasSomeVar rhs || isComplex lhs))
return (collectLikeTerms (sorted (lhs - rhs)) :==: 0)
where
isComplex = maybe False ((>= 2) . length . filter hasSomeVar)
. match sumView . applyD merge
-- high exponents first, non power-factor terms at the end
sorted = simplifyWith (sortBy (comparing toPF)) sumView
toPF = fmap (negate . thd3) . match powerFactorView
ruleApproximate :: Rule (Relation Expr)
ruleApproximate = describe "Approximate irrational number" $
makeSimpleRule (quadreq, "approx") $ \relation -> do
lhs :==: rhs <- match equationView relation
guard (not (simplify rhs `belongsTo` rationalView))
x <- getVariable lhs
d <- match doubleView rhs
let new = fromDouble (precision 4 d)
return (Var x .~=. new)
ruleNormalizeRational :: Rule Expr
ruleNormalizeRational =
describe "normalize rational number" $
ruleFromView (lineq, "norm-rational") rationalView
ruleNormalizeMixedFraction :: Rule Expr
ruleNormalizeMixedFraction =
describe "normalize mixed fraction" $
ruleFromView (lineq, "norm-mixed") mixedFractionView
-----------------------------------------------------------
-------- Rules From HDE
-- X*A + X*B = X*C + X*D
-- New implementation, but slightly different than original
-- This one does not factor constants
allPowerFactors :: Rule (OrList (Equation Expr))
allPowerFactors = describe "all power factors" $
makeSimpleRule (polyeq, "power-factors") $ oneDisjunct $
\(lhs :==: rhs) -> do
let myView = polyNormalForm rationalView
(s1, p1) <- match myView lhs
(s2, p2) <- match myView rhs
let n | p1 == 0 = lowestDegree p2
| p2 == 0 = lowestDegree p1
| otherwise = lowestDegree p1 `min` lowestDegree p2
ts = terms p1 ++ terms p2
f p = build myView (s1, raise (-n) p)
guard ((s1==s2 || p1==0 || p2==0) && n > 0 && length ts > 1)
return $ orList [Var s1 :==: 0, f p1 :==: f p2]
factorVariablePower :: Rule Expr
factorVariablePower = describe "factor variable power" $
makeSimpleRule (polyeq, "factor-varpower") $ \expr -> do
let myView = polyNormalForm rationalView
(s, p) <- match (polyNormalForm rationalView) expr
let n = lowestDegree p
guard (n > 0 && length (terms p) > 1)
return $ Var s .^. fromIntegral n * build myView (s, raise (-n) p)
-- A*B = A*C implies A=0 or B=C
sameFactor :: Rule (OrList (Equation Expr))
sameFactor = describe "same factor" $
makeSimpleRule (quadreq, "same-factor") $ oneDisjunct $ \(lhs :==: rhs) -> do
(b1, xs) <- match productView lhs
(b2, ys) <- match productView rhs
(x, y) <- safeHead [ (x, y) | x <- xs, y <- ys, x==y, hasSomeVar x ] -- equality is too strong?
return $ orList [ x :==: 0, build productView (b1, xs\\[x]) :==: build productView (b2, ys\\[y]) ]
-- N*(A+B) = N*C + N*D recognize a constant factor on both sides
-- Example: 3(x^2+1/2) = 6+6x
sameConFactor :: Rule (Equation Expr)
sameConFactor =
describe "same constant factor" $
liftRule myView $
makeSimpleRule (quadreq, "same-con-factor") $ \(ps1 :==: ps2) -> do
let (bs, zs) = unzip (ps1 ++ ps2)
(rs, es) = unzip (map (f 1 []) zs)
f r acc [] = (r, reverse acc)
f r acc (x:xs) = case match rationalView x of
Just r2 -> f (r*r2) acc xs
Nothing -> f r (x:acc) xs
c <- whichCon rs
guard (c /= 1)
let make b r e = (b, fromRational (r/c):e)
(newLeft, newRight) = splitAt (length ps1) (zipWith3 make bs rs es)
return (newLeft :==: newRight)
where
myView = bothSidesView (sumView >>> listView productView)
whichCon :: [Rational] -> Maybe Rational
whichCon xs
| all (\x -> denominator x == 1 && x /= 0) xs =
Just (fromInteger (foldr1 gcd (map numerator xs)))
| otherwise = Nothing
abcFormula :: Rule (Context (OrList (Equation Expr)))
abcFormula = describe "quadratic formula (abc formule)" $
makeSimpleRule (quadreq, "abc") $ withCM $ oneDisjunct $ \(lhs :==: rhs) -> do
guard (rhs == 0)
(x, (a, b, c)) <- matchM quadraticNF lhs
addListToClipboard ["a", "b", "c"] (map fromRational [a, b, c])
let discr = b*b - 4 * a * c
sqD = sqrt (fromRational discr)
addToClipboard "D" (fromRational discr)
case compare discr 0 of
LT -> return false
EQ -> return $ return $
Var x :==: (-fromRational b) / (2 * fromRational a)
GT -> return $ orList
[ Var x :==: (-fromRational b + sqD) / (2 * fromRational a)
, Var x :==: (-fromRational b - sqD) / (2 * fromRational a)
]
higherSubst :: Rule (Context (Equation Expr))
higherSubst = describe "Substitute variable" $
makeSimpleRule (polyeq, "subst") $ withCM $ \(lhs :==: rhs) -> do
guard (rhs == 0)
let myView = polyView >>> second trinomialPolyView
(x, ((a, n1), (b, n2), (c, n3))) <- matchM myView lhs
guard (n1 == 0 && n2 > 1 && n3 `mod` n2 == 0 && x /= "p")
let new = build myView ("p", ((a, 0), (b, 1), (c, n3 `div` n2)))
addToClipboard "subst" (toExpr (Var "p" :==: Var x .^. fromIntegral n2))
return (new :==: 0)
substBackVar :: Rule (Context Expr)
substBackVar = describe "Substitute back a variable" $
makeSimpleRule (polyeq, "back-subst") $ withCM $ \a -> do
expr <- lookupClipboard "subst"
case fromExpr expr of
Just (Var p :==: rhs) -> do
guard (hasVar p a)
return (subst p rhs a)
_ -> fail "no subst in clipboard"
where
subst a b (Var c) | a==c = b
subst a b expr = descend (subst a b) expr
exposeSameFactor :: Rule (Equation Expr)
exposeSameFactor = describe "expose same factor" $
liftRule (bothSidesView productView) $
makeSimpleRuleList (polyeq, "expose-factor") $ \((bx, xs) :==: (by, ys)) -> do
(nx, ny) <- [ (xs, new) | x <- xs, suitable x, new <- exposeList x ys ] ++
[ (new, ys) | y <- ys, suitable y, new <- exposeList y xs ]
return ((bx, nx) :==: (by, ny))
where
suitable p = fromMaybe False $ do
(_, _, b) <- match (linearViewWith rationalView) p
guard (b /= 0)
return True
exposeList _ [] = []
exposeList a (b:bs) = map (++bs) (expose a b) ++ map (b:) (exposeList a bs)
expose a b = do
(s1, p1) <- matchM (polyViewWith rationalView) a
(s2, p2) <- matchM (polyViewWith rationalView) b
guard (s1==s2 && p1/=p2)
case division p2 p1 of
Just p3 -> return $ map (\p -> build (polyViewWith rationalView) (s1,p)) [p1, p3]
Nothing -> []
---------------------------------------------------------
-- From LinearEquations
-- Only used for cleaning up
distributeAll :: Expr -> Expr
distributeAll expr =
case expr of
e1 :*: e2 -> let as = fromMaybe [e1] (match sumView e1)
bs = fromMaybe [e2] (match sumView e2)
in build sumView [ a .*. b | a <- as, b <- bs ]
_ -> expr
-- This rule should consider the associativity of multiplication
-- Combine bottom-up, for example: 5*(x-5)*(x+5)
-- However, in -2x(2x+10) (-2x) should be seen as "one term"
distributionT :: Transformation Expr
distributionT = makeTransList f
where
f expr = do
(b, xs) <- matchM simpleProductView expr
ys <- rec (combine xs)
return $ build simpleProductView (b, ys)
combine :: [Expr] -> [Expr]
combine (x:y:rest) | p x && p y = combine ((x*y):rest)
where p = maybe False ((==1) . length) . match sumView
combine [] = []
combine (x:xs) = x : combine xs
rec :: [Expr] -> [[Expr]]
rec (a:b:xs) = map (:xs) (g a b) ++ map (a:) (rec (b:xs))
rec _ = []
g :: Expr -> Expr -> [Expr]
g e1 e2 = do
as <- matchM sumView e1
bs <- matchM sumView e2
guard (length as > 1 || length bs > 1)
return $ build sumView [ a .*. b | a <- as, b <- bs ]
-------------------------------------------------------
-- Rewrite Rules
varToLeft :: Rule (Relation Expr)
varToLeft = doAfter (fmap collectLikeTerms) $
describe "variable to left" $
makeRule (lineq, "var-left") $ flip supply1 minusT $ \eq -> do
(x, a, _) <- match (linearViewWith rationalView) (rightHandSide eq)
guard (a/=0)
return (fromRational a * Var x)
-- factor is always positive due to lcm function
removeDivision :: Rule (Relation Expr)
removeDivision = doAfter (fmap (collectLikeTerms . distributeAll)) $
describe "remove division" $
makeRule (lineq, "remove-div") $ flip supply1 timesT $ \eq -> do
xs <- match sumView (leftHandSide eq)
ys <- match sumView (rightHandSide eq)
-- also consider parts without variables
-- (but at least one participant should have a variable)
zs <- forM (xs ++ ys) $ \a -> do
(_, list) <- match productView a
return [ (hasSomeVar a, e) | e <- list ]
let f (b, e) = do
(_, this) <- match (divView >>> second integerView) e
return (b, this)
case mapMaybe f (concat zs) of
[] -> Nothing
ps -> let (bs, ns) = unzip ps
in if or bs then return (fromInteger (foldr1 lcm ns))
else Nothing
distributeTimes :: Rule Expr
distributeTimes = describe "distribution multiplication" $
makeSimpleRuleList (lineq, "distr-times") $
liftM collectLikeTerms . applyAll distributionT
distributeDivision :: Rule Expr
distributeDivision = describe "distribution division" $
makeSimpleRule (quadreq, "distr-div") $ \expr -> do
(xs, r) <- match (divView >>> (sumView *** rationalView)) expr
guard (length xs > 1)
let ys = map (/fromRational r) xs
return $ build sumView ys
merge :: Rule Expr
merge = describe "merge similar terms" $
makeSimpleRule (lineq, "merge") $ \old -> do
let new = collectLikeTerms old
f = maybe 0 length . match sumView
guard (f old > f new)
return new
simplerLinearFactor :: Rule Expr
simplerLinearFactor = describe "simpler linear factor" $
makeSimpleRule (polyeq, "simpler-linfactor") $ \expr -> do
let myView = polyNormalForm rationalView >>> second linearPolyView
(x, (a, b)) <- match myView expr
let d = (if a<0 then negate else id) (gcdFrac a b)
guard (a /= 0 && b /= 0 && d /= 1 && d /= -1)
return $ fromRational d * build myView (x, (a/d, b/d))
ruleFromView :: (IsId n, Eq a) => n -> View a b -> Rule a
ruleFromView s v = makeSimpleRuleList s $ \a -> do
b <- canonicalM v a
guard (a /= b)
return b
rhsIsZero :: Rule Expr -> Rule (Equation Expr)
rhsIsZero r = makeSimpleRuleList (showId r) $ \(lhs :==: rhs) -> do
guard (rhs == 0)
a <- applyAll r lhs
return (a :==: rhs)
constantRight :: View Expr a -> View (Equation Expr) (a, Rational)
constantRight v = makeView f g
where
f (lhs :==: rhs) = liftM2 (,) (match v lhs) (match rationalView rhs)
g (a, r) = build v a :==: build rationalView r
bothSidesView :: View a b -> View (Equation a) (Equation b)
bothSidesView v = makeView f (fmap (build v))
where
f (lhs :==: rhs) = liftM2 (:==:) (match v lhs) (match v rhs)
findFactor :: Monad m => [Rational] -> m Rational
findFactor rs
| null rs =
fail "no factor"
| all ((==1) . denominator) rs =
return $ Prelude.recip $ fromIntegral $ foldr1 gcd $ map numerator rs
| otherwise =
return $ fromIntegral $ foldr1 lcm $ map denominator rs
parentNotNegCheck :: Rule (Context Expr)
parentNotNegCheck = minorRule $ makeSimpleRule "parent not negate check" $ \c ->
case up c >>= current of
Just (Negate _) -> Nothing
_ -> Just c
noDivisionConstant :: Rule Expr
noDivisionConstant = makeSimpleRule (lineq, "no-div-con") f
where
f (a :/: b) | hasNoVar b && hasSomeVar a =
return ((1/b) * a)
f _ = Nothing
defPowerNat :: Rule Expr
defPowerNat = makeSimpleRule (polyeq, "def-power-nat") f
where
f (Sym _ [Var _, _]) = Nothing -- should not work on x^5
f (Sym s [a, Nat n]) | isPowerSymbol s =
return (build productView (False, replicate (fromInteger n) a))
f _ = Nothing