Merge pull request #351 from JuliaSymbolics/s/flatten_frac

`flatten_fraction`

Author: shashi

src/polyform.jl

@@ -1,4 +1,4 @@
-export PolyForm, simplify_fractions, quick_cancel
+export PolyForm, simplify_fractions, quick_cancel, flatten_fractions
 using Bijections
 using DynamicPolynomials: PolyVar
 
@@ -7,7 +7,8 @@ using DynamicPolynomials: PolyVar
 
 Abstracts a [MultivariatePolynomials.jl](https://juliaalgebra.github.io/MultivariatePolynomials.jl/stable/) as a SymbolicUtils expression and vice-versa.
 
-The SymbolicUtils term interface (`istree`, `operation, and `arguments`) works on PolyForm lazily: the `operation` and `arguments` are created by converting one level of arguments into SymbolicUtils expressions. They may further contain PolyForm within them.
+The SymbolicUtils term interface (`istree`, `operation, and `arguments`) works on PolyForm lazily:
+the `operation` and `arguments` are created by converting one level of arguments into SymbolicUtils expressions. They may further contain PolyForm within them.
 We use this to hold polynomials in memory while doing `simplify_fractions`.
 
     PolyForm{T}(x; Fs=Union{typeof(*),typeof(+),typeof(^)}, recurse=false)
@@ -190,15 +191,16 @@ function TermInterface.arguments(x::PolyForm{T}) where {T}
 
     if MP.nterms(x.p) == 1
         MP.isconstant(x.p) && return [convert(Number, x.p)]
-        c = MP.coefficient(x.p)
-        t = MP.monomial(x.p)
+        t = MP.term(x.p)
+        c = MP.coefficient(t)
+        m = MP.monomial(t)
 
         if !isone(c)
             [c, (unstable_pow(resolve(v), pow)
-                        for (v, pow) in MP.powers(t) if !iszero(pow))...]
+                        for (v, pow) in MP.powers(m) if !iszero(pow))...]
         else
             [unstable_pow(resolve(v), pow)
-                    for (v, pow) in MP.powers(t) if !iszero(pow)]
+                    for (v, pow) in MP.powers(m) if !iszero(pow)]
         end
     else
         ts = MP.terms(x.p)
@@ -223,12 +225,12 @@ Expand expressions by distributing multiplication over addition, e.g.,
 multivariate polynomials implementation.
 `variable_type` can be any subtype of `MultivariatePolynomials.AbstractVariable`.
 """
-expand(expr) = PolyForm(expr, Fs=Union{typeof(+), typeof(*), typeof(^)}, recurse=true)
+expand(expr) = Postwalk(identity)(PolyForm(expr, Fs=Union{typeof(+), typeof(*), typeof(^)}, recurse=true))
 
 
 ## Rational Polynomial form with Div
 
-function polyform_factors(d::Div, pvar2sym, sym2term)
+function polyform_factors(d, pvar2sym, sym2term)
     make(xs) = map(xs) do x
         if x isa Pow && x.base isa Integer && x.exp > 0
             # here we do want to recurse one level, that's why it's wrong to just
@@ -255,11 +257,11 @@ function simplify_div(d::Div)
     end
 end
 
-function add_divs(x::Div, y::Div)
+function add_divs(x, y)
     x_num, x_den = polyform_factors(x, get_pvar2sym(), get_sym2term())
     y_num, y_den = polyform_factors(y, get_pvar2sym(), get_sym2term())
 
-    Div(_mul(x_num, y_den) + _mul(x_den, y_num), _mul(x_den, y_den))
+    (_mul(x_num, y_den) + _mul(x_den, y_num)) / (_mul(x_den, y_den))
 end
 
 """
@@ -278,9 +280,38 @@ function simplify_fractions(x)
     rules = [@rule ~x::isdiv => simplify_div(~x)
              @acrule ~a::isdiv + ~b::isdiv => add_divs(~a,~b)]
 
+    Fixpoint(Postwalk(RestartedChain(rules)))(x)
+end
+
+"""
+    flatten_fractions(x)
+
+Flatten nested fractions that are added together.
+
+```julia
+julia> flatten_fractions((1+(1+1/a)/a)/a)
+(1 + a + a^2) / (a^3)
+```
+"""
+function flatten_fractions(x)
+    rules = [@acrule ~a::(x->x isa Div) + ~b => add_divs(~a,~b)
+             @rule *(~~x, ~a / ~b, ~~y) / ~c => *((~~x)..., ~a, (~~y)...) / (~b * ~c)
+             @rule ~c / *(~~x, ~a / ~b, ~~y)  => (~b * ~c) / *((~~x)..., ~a, (~~y)...)]
     Fixpoint(Postwalk(Chain(rules)))(x)
 end
 
+function fraction_iszero(x)
+    !istree(x) && return _iszero(x)
+    # fast path and then slow path
+    any(_iszero, numerators(flatten_fractions(x))) ||
+    any(_iszero∘expand, numerators(flatten_fractions(x)))
+end
+
+function fraction_isone(x)
+    !istree(x) && return _isone(x)
+    _isone(simplify_fractions(flatten_fractions(x)))
+end
+
 function needs_div_rules(x)
     (x isa Div && !(x.num isa Number) && !(x.den isa Number)) ||
     (istree(x) && operation(x) === (+) && count(has_div, unsorted_arguments(x)) > 1) ||

src/types.jl

@@ -857,6 +857,15 @@ maybe_intcoeff(x) = x
 function (::Type{Div{T}})(n, d, simplified=false; metadata=nothing) where {T}
     _iszero(n) && return zero(typeof(n))
     _isone(d) && return n
+
+    if n isa Div && d isa Div
+        return Div{T}(n.num * d.den, n.den * d.num)
+    elseif n isa Div
+        return Div{T}(n.num, n.den * d)
+    elseif d isa Div
+        return Div{T}(n * d.den, d.num)
+    end
+
     d isa Number && _isone(-d) && return -1 * n
     n isa Rat && d isa Rat && return n // d # maybe called by oblivious code in simplify
 
@@ -879,15 +888,11 @@ function Div(n,d, simplified=false; kw...)
     Div{promote_symtype((/), symtype(n), symtype(d))}(n,d, simplified; kw...)
 end
 
-function numerators(d::Div)
-    x = d.num
-    istree(x) && operation(x) == (*) ? arguments(x) : [x]
-end
+numerators(x) = istree(x) && operation(x) == (*) ? arguments(x) : [x]
+numerators(d::Div) = numerators(d.num)
 
-function denominators(d::Div)
-    x = d.den
-    istree(x) && operation(x) == (*) ? arguments(x) : [x]
-end
+denominators(x) = [1]
+denominators(d::Div) = numerators(d.den)
 
 TermInterface.istree(d::Type{Div}) = true
 
@@ -899,17 +904,7 @@ end
 
 Base.show(io::IO, d::Div) = show_term(io, d)
 
-function /(a::Union{SN,Number}, b::SN)
-    if a isa Div && b isa Div
-        Div(a.num * b.den, a.den * b.num)
-    elseif a isa Div
-        Div(a.num, a.den * b)
-    elseif b isa Div
-        Div(a * b.den, b.num)
-    else
-        Div(a,b)
-    end
-end
+/(a::Union{SN,Number}, b::SN) = Div(a,b)
 
 """
     Pow(base, exp)

test/polyform.jl

@@ -53,3 +53,16 @@ end
     @eqtest simplify_fractions(3*(x^2)*(y^3)/(3*(x^3)*(y^2))) == y/x
     @eqtest simplify_fractions(3*(x^x)/x*y) == 3*(x^x)/x*y
 end
+
+@testset "isone iszero" begin
+    @syms a b c d e f g h i
+    x = (f + ((((g*(c^2)*(e^2)) / d - e*h*(c^2)) / b + (-c*e*f*g) / d + c*e*i) /
+              (i + ((c*e*g) / d - c*h) / b + (-f*g) / d) - c*e) / b +
+         ((g*(f^2)) / d + ((-c*e*f*g) / d + c*f*h) / b - f*i) /
+         (i + ((c*e*g) / d - c*h) / b + (-f*g) / d)) / d
+
+    o = (d + (e*((c*(g + (-d*g) / d)) / (i + (-c*(h + (-e*g) / d)) / b + (-f*g) / d))) / b + (-f*(g + (-d*g) / d)) / (i + (-c*(h + (-e*g) / d)) / b + (-f*g) / d)) / d
+    @test SymbolicUtils.fraction_iszero(x)
+    @test !SymbolicUtils.fraction_isone(x)
+    @test SymbolicUtils.fraction_isone(o)
+end