Simplify div of polynomial denominator (#118)

Co-authored-by: Benoît Legat <[email protected]>

Author: odow

src/nl_to_polynomial.jl

@@ -94,21 +94,38 @@ function _to_polynomial!(
         return a^b
     elseif _is_operator(expr, :/) && length(operands) == 2
         a, b = _to_polynomial!.(Ref(d), T, operands)
-        return a / b
-        # TODO(odow): see PR#113
-        # divisor, remainder = Base.divrem(a, b)
-        # if iszero(remainder)
-        #     return divisor
-        # else
-        #     return a / b
-        # end
+        return _checked_div(a, b)
     elseif _is_variable(expr)
         return _to_polynomial!(d, T, operands[1])
     else
         throw(InvalidNLExpression("Cannot convert `$(expr)` into a polynomial"))
     end
 end
 
+_checked_div(a, b) = a / b
+
+function _checked_div(α, p::MP.AbstractPolynomialLike)
+    if !MP.isconstant(p)
+        throw(InvalidNLExpression("Cannot convert `α / ($p)` to a polynomial."))
+    end
+    return MP.constant_term(α / MP.convert_to_constant(p), p)
+end
+
+function _checked_div(
+    a::MP.AbstractPolynomialLike,
+    b::MP.AbstractPolynomialLike,
+)
+    divisor, remainder = Base.divrem(a, b)
+    if !iszero(remainder)
+        throw(
+            InvalidNLExpression(
+                "Cannot convert `($(a)) / ($b)` into a polynomial as the Euclidean division has a nonzero remainder `$remainder` (the divisor is `$divisor`).",
+            ),
+        )
+    end
+    return divisor
+end
+
 function _to_polynomial(expr, ::Type{T}) where {T}
     d = Dict{MOI.VariableIndex,VarType}()
     poly = _to_polynomial!(d, T, expr)

test/qcqp.jl

@@ -215,6 +215,52 @@ function test_scalar_nonlinear_function(x, y, T)
     return
 end
 
+function test_scalar_nonlinear_function_div_rem_zero(x, y, T)
+    inner = Model{T}()
+    model = PolyJuMP.JuMP.GenericModel{T}() do
+        return PolyJuMP.QCQP.Optimizer{T}(MOI.Utilities.MockOptimizer(inner))
+    end
+    PolyJuMP.@variable(model, x)
+    PolyJuMP.@objective(model, Min, x^3 / x)
+    PolyJuMP.optimize!(model)
+    @test MOI.get(inner, MOI.NumberOfVariables()) == 1
+    @test isempty(MOI.get(inner, MOI.ListOfConstraintTypesPresent()))
+    @test PolyJuMP.objective_function(model) isa PolyJuMP.GenericNonlinearExpr
+    F = MOI.get(inner, MOI.ObjectiveFunctionType())
+    @test F <: MOI.ScalarQuadraticFunction{T}
+    return
+end
+
+function test_scalar_nonlinear_function_div_rem_err(x, y, T)
+    inner = Model{T}()
+    model = PolyJuMP.JuMP.GenericModel{T}() do
+        return PolyJuMP.QCQP.Optimizer{T}(MOI.Utilities.MockOptimizer(inner))
+    end
+    PolyJuMP.@variable(model, x)
+    PolyJuMP.@variable(model, y)
+    PolyJuMP.@objective(model, Min, x^3 / y)
+    @test_throws PolyJuMP.InvalidNLExpression PolyJuMP.optimize!(model)
+    PolyJuMP.@objective(model, Min, 2 / y)
+    @test_throws PolyJuMP.InvalidNLExpression PolyJuMP.optimize!(model)
+    PolyJuMP.@objective(model, Min, 2 / (x * y^2 - y^2 * x - 1))
+    PolyJuMP.optimize!(model)
+    return
+end
+
+function test_scalar_nonlinear_function_div_rem_number(x, y, T)
+    inner = Model{T}()
+    model = PolyJuMP.JuMP.GenericModel{T}() do
+        return PolyJuMP.QCQP.Optimizer{T}(MOI.Utilities.MockOptimizer(inner))
+    end
+    PolyJuMP.@variable(model, x)
+    PolyJuMP.@objective(model, Min, 4x^3 / 2)
+    PolyJuMP.optimize!(model)
+    @test MOI.get(inner, MOI.NumberOfVariables()) == 2
+    F, S = MOI.ScalarQuadraticFunction{T}, MOI.EqualTo{T}
+    @test (F, S) in MOI.get(inner, MOI.ListOfConstraintTypesPresent())
+    return
+end
+
 function test_variable_primal(x, y, T)
     inner = Model{T}()
     optimizer = MOI.Utilities.MockOptimizer(inner)