module Domain.Math.Polynomial.Rules
( sameConFactor, abcFormula, allPowerFactors, bringAToOne, cancelTerms
, commonFactorVar, commonFactorVarNew, defPowerNat
, distributeDivisionT, distributeDivisionMulti
, distributeTimes, distributionSquare, exposeSameFactor, factorLeftAsSquare
, factorVariablePower, flipEquation, higherSubst, merge, moveToLeft, mulZero
, niceFactors, niceFactorsNew, noDivisionConstant, noLinFormula, oneVar
, prepareSplitSquare, quadraticRuleOrder, removeDivision
, ruleApproximate, ruleNormalizeMixedFraction, ruleNormalizeRational
, ruleNormalizePolynomial
, sameFactor, simplerLinearFactor, simplerPolynomial, simplerSquareRootMulti
, squareBothSides, substBackVar, varToLeft, conditionVarsRHS, fractionProduct
) where
import Control.Monad
import Data.List
import Data.Maybe
import Data.Ord
import Data.Ratio
import Domain.Math.Approximation (precision)
import Domain.Math.CleanUp
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.Views
import Domain.Math.Power.OldViews (powerFactorView)
import Domain.Math.Safe
import Domain.Math.Simplification hiding (simplifyWith, distribution)
import Domain.Math.SquareRoot.Views
import Ideas.Common.Library hiding (terms, simplify, merge, (.*.), (./.))
import Ideas.Utils.Prelude (thd3)
import Ideas.Utils.Uniplate (universe, descend)
import Prelude hiding ( (^) )
import qualified Prelude
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"
quadraticNF :: View Expr (String, (Rational, Rational, Rational))
quadraticNF = polyNormalForm rationalView >>> second quadraticPolyView
commonFactorVar :: Rule (Equation Expr)
commonFactorVar = rhsIsZero commonFactorVarNew
commonFactorVarNew :: Rule Expr
commonFactorVarNew = describe "Common factor variable" $
makeRule (quadreq, "common-factor") $ \expr -> do
(x, (a, b, c)) <- match quadraticNF expr
guard (a /= 0 && b /= 0 && c == 0)
let d = signum a * 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
noLinFormula :: Rule (Equation Expr)
noLinFormula = describe "No linear term ('b=0')" $ liftView myView $
ruleMaybe (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
niceFactors :: Rule (Equation Expr)
niceFactors = rhsIsZero niceFactorsNew
niceFactorsNew :: Rule Expr
niceFactorsNew = describe "Find a nice decomposition" $
makeRule (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)) <- sign <$> matchM (polyNormalForm integerView >>> second quadraticPolyView) expr
guard (a==1)
let ok (i, j) = i+j == b
f (i, j)
| i == j =
(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)]
]
simplerPolynomial :: Rule (Equation Expr)
simplerPolynomial = describe "simpler polynomial" $
rhsIsZero $ liftViewIn (quadraticNF >>> toView swapView) $
makeRule (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)
bringAToOne :: Rule (Equation Expr)
bringAToOne = rhsIsZero $ liftViewIn (quadraticNF >>> toView swapView) $
describe "Bring 'a' to one" $
ruleMaybe (quadreq, "scale") $ \(a, b, c) -> do
guard (a `notElem` [0, 1])
return (1, b/a, c/a)
mulZero :: Rule (OrList (Equation Expr))
mulZero = describe "multiplication is zero" $
makeRule (quadreq, "product-zero") $ oneDisjunct bothSides
where
bothSides eq = oneSide eq ++ oneSide (flipSides eq)
oneSide (lhs :==: rhs) = do
guard (rhs == 0)
(_, xs) <- matchM productView lhs
guard (length xs > 1)
return $ toOrList $ flip map xs $ \e ->
case match (polyNormalForm rationalView >>> second linearPolyView) e of
Just (x, (a, b))
| a == 1 ->
Var x :==: fromRational (-b)
| a == -1 ->
Var x :==: fromRational b
_ -> e :==: 0
oneVar :: ConfigCoverUp
oneVar = configCoverUp
{ configName = "onevar"
, predicateCovered = \a -> p1 a || p2 a
, predicateCombined = const True
, coverLHS = True
, coverRHS = True
}
where
p1 = (\x -> x == 1) . length . vars
p2 a = fromMaybe False $ do
(x, y) <- match timesView a
return (hasSomeVar x /= hasSomeVar y)
simplerSquareRootMulti :: IsTerm a => Rule (Context a)
simplerSquareRootMulti = describe "simpler square root" $
makeRule (quadreq, "simpler-sqrt") $ applyAll $
repeat1 (somewhere (use (makeRule () simplerSqrt)))
where
simplerSqrt :: Expr -> Maybe Expr
simplerSqrt e = do
xs <- f e
guard (not (null xs))
new <- canonical (squareRootViewWith rationalView) e
ys <- f new
guard (xs /= ys)
return new
f :: Expr -> Maybe [Rational]
f e = sort <$> sequence [ match rationalView a | Sqrt a <- universe e ]
cancelTerms :: Rule (Equation Expr)
cancelTerms = describe "Cancel terms" $
makeRule (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)
distributionSquare :: Rule Expr
distributionSquare = describe "distribution for special products" $
rewriteRules (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 b -> (a+b)*(a-b) :~> a^2 - b^2
, \a b -> (a-b)*(a+b) :~> a^2 - b^2
]
squareBothSides :: Rule (OrList (Equation Expr))
squareBothSides = describe "square both sides" $
rewriteRule (quadreq, "square-both") $ \a b ->
singleton (a^2 :==: b^2) :~> toOrList [a :==: b, a :==: -b]
prepareSplitSquare :: Rule (Equation Expr)
prepareSplitSquare = describe "prepare split square" $
liftView myView $
ruleMaybe (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
factorLeftAsSquare :: Rule (Equation Expr)
factorLeftAsSquare = describe "factor left as square" $
makeRule (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)
flipEquation :: Rule (Equation Expr)
flipEquation = describe "flip equation" $
rewriteRule (lineq, "flip") $ \a b ->
(a :==: b) :~> (b :==: a)
conditionVarsRHS :: Equation Expr -> Bool
conditionVarsRHS (lhs :==: rhs) = hasSomeVar rhs && hasNoVar lhs
moveToLeft :: Rule (Equation Expr)
moveToLeft = describe "Move to left" $
ruleMaybe (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
sorted = simplifyWith (sortBy (comparing toPF)) sumView
toPF = fmap (negate . thd3) . match powerFactorView
ruleApproximate :: Rule (Relation Expr)
ruleApproximate = describe "Approximate irrational number" $
makeRule (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
ruleNormalizePolynomial :: Rule Expr
ruleNormalizePolynomial =
describe "normalize polynomial" $
ruleFromView (polyeq, "norm-poly") (polyNormalForm rationalView)
allPowerFactors :: Rule (OrList (Equation Expr))
allPowerFactors = describe "all power factors" $
makeRule (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 = filter (/= 0) (fromPolynomial p1 ++ fromPolynomial p2)
f p = build myView (s1, raise (-n) p)
guard ((s1==s2 || p1==0 || p2==0) && n > 0 && length ts > 1)
return $ toOrList [Var s1 :==: 0, f p1 :==: f p2]
factorVariablePower :: Rule Expr
factorVariablePower = describe "factor variable power" $
makeRule (polyeq, "factor-varpower") $ \expr -> do
let myView = polyNormalForm rationalView
(s, p) <- match (polyNormalForm rationalView) expr
let n = lowestDegree p
guard (n > 0 && length (filter (/=0) (fromPolynomial p)) > 1)
new <- p `safeDiv` (var Prelude.^ n)
return $ Var s .^. fromIntegral n * build myView (s, new)
sameFactor :: Rule (OrList (Equation Expr))
sameFactor = describe "same factor" $
makeRule (quadreq, "same-factor") $ oneDisjunct $ \(lhs :==: rhs) -> do
(b1, xs) <- match productView lhs
(b2, ys) <- match productView rhs
(x, y) <- listToMaybe [ (x, y) | x <- xs, y <- ys, x==y, hasSomeVar x ]
return $ toOrList [ x :==: 0, build productView (b1, xs\\[x]) :==: build productView (b2, ys\\[y]) ]
sameConFactor :: Rule (Equation Expr)
sameConFactor =
describe "same constant factor" $
liftView myView $
makeRule (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 (toView sumView >>> listView (toView 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)" $
makeRule (quadreq, "abc") $ \cor -> do
oreq <- currentInContext cor
lhs :==: rhs <- getSingleton oreq
guard (rhs == 0)
(x, (a, b, c)) <- matchM quadraticNF lhs
let discr = b*b - 4 * a * c
sqD = sqrt (fromRational discr)
eqs = case compare discr 0 of
LT -> false
EQ -> singleton $
Var x :==: (-fromRational b) / (2 * fromRational a)
GT -> toOrList
[ Var x :==: (-fromRational b + sqD) / (2 * fromRational a)
, Var x :==: (-fromRational b - sqD) / (2 * fromRational a)
]
return $ addToClipboard "a" (fromRational a)
$ addToClipboard "b" (fromRational b)
$ addToClipboard "c" (fromRational c)
$ addToClipboard "D" (fromRational discr)
$ replaceInContext eqs cor
higherSubst :: Rule (Context (Equation Expr))
higherSubst = describe "Substitute variable" $
ruleMaybe (polyeq, "subst") $ \ceq -> do
lhs :==: rhs <- fromContext ceq
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)))
return $ addToClipboard "subst" (toExpr (Var "p" :==: Var x .^. fromIntegral n2))
$ replaceInContext (new :==: 0) ceq
substBackVar :: Rule (Context Expr)
substBackVar = describe "Substitute back a variable" $
makeRule (polyeq, "back-subst") $ \ca -> do
a <- fromContext ca
expr <- lookupClipboard "subst" ca
case fromExpr expr of
Just (Var p :==: rhs) -> do
guard (hasVar p a)
return (replaceInContext (subst p rhs a) ca)
_ -> 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" $
liftView (bothSidesView (toView productView)) $
makeRule (polyeq, "expose-factor") $ \((bx, xs) :==: (by, ys)) -> do
(nx, ny) <- [ (xs, new) | x <- xs, isOk x, new <- exposeList x ys ] ++
[ (new, ys) | y <- ys, isOk y, new <- exposeList y xs ]
return ((bx, nx) :==: (by, ny))
where
isOk 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 safeDiv p2 p1 of
Just p3 -> return $ map (\p -> build (polyViewWith rationalView) (s1,p)) [p1, p3]
Nothing -> []
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
distribution :: Expr -> [Expr]
distribution expr = do
(b, xs) <- matchM simpleProductView expr
ys <- rec (combine xs)
return $ build simpleProductView (b, ys)
where
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 ]
varToLeft :: Rule (Relation Expr)
varToLeft = doAfter (fmap collectLikeTerms) $
describe "variable to left" $
ruleTrans (lineq, "var-left") $
inputWith arg minusRule
where
arg = transMaybe $ \eq -> do
(x, a, _) <- matchM (linearViewWith rationalView) (rightHandSide eq)
guard (a/=0)
return (fromRational a * Var x)
removeDivision :: Rule (Relation Expr)
removeDivision = doAfter (fmap (collectLikeTerms . distributeAll)) $
describe "remove division" $
ruleTrans (lineq, "remove-div") $
inputWith arg timesRule
where
arg = transMaybe $ \eq -> do
xs <- matchM sumView (leftHandSide eq)
ys <- matchM sumView (rightHandSide eq)
zs <- forM (xs ++ ys) $ \a -> do
(_, list) <- matchM productView a
return [ (hasSomeVar a, e) | e <- list ]
let f (b, e) = do
(_, this) <- match (divView >>> second integerView) e
return (b, this)
(bs, ns) = unzip (mapMaybe f (concat zs))
guard (or bs)
return (fromInteger (foldr1 lcm ns))
distributeTimes :: Rule Expr
distributeTimes = describe "distribution multiplication" $
makeRule (lineq, "distr-times") $
fmap collectLikeTerms . distribution
distributeDivisionMulti :: IsTerm a => Rule (Context a)
distributeDivisionMulti = describe "distribution division" $
makeRule (quadreq, "distr-div") $ apply $ repeat1 $
somewhere (use (makeRule () distributeDivisionT))
distributeDivisionT :: Expr -> Maybe Expr
distributeDivisionT expr = do
(xs, r) <- match (divView >>> (toView 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" $
ruleMaybe (lineq, "merge") $ \old -> do
let norm = cleanUpSimple old
new = collectLikeTerms norm
f = maybe 0 length . match sumView
guard (f norm > f new)
return new
simplerLinearFactor :: Rule Expr
simplerLinearFactor = describe "simpler linear factor" $
makeRule (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 `notElem` [1, -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 = makeRule s $ \a -> do
b <- canonical v a
guard (a /= b)
return b
rhsIsZero :: Rule Expr -> Rule (Equation Expr)
rhsIsZero r = makeRule (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) = (,) <$> 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) = (:==:) <$> 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
noDivisionConstant :: Rule Expr
noDivisionConstant = makeRule (lineq, "no-div-con") f
where
f (a :/: b) | hasNoVar b && hasSomeVar a =
return ((1/b) * a)
f _ = Nothing
fractionProduct :: Rule Expr
fractionProduct = makeRule (polyeq, "fraction-product") $ \expr -> do
((a, b), (c, d)) <- match (timesView >>> divView *** divView) expr
return ((a .*. c) ./. (b .*. d))
defPowerNat :: Rule Expr
defPowerNat = makeRule (polyeq, "def-power-nat") f
where
f (Sym _ [Var _, _]) = Nothing
f (Sym s [a, Nat n]) | isPowerSymbol s =
return (build productView (False, replicate (fromInteger n) a))
f _ = Nothing