Using the Positivstellensatz, you know that the set `{(x, y) | f1(x, y) <= 0, f2(x, y) <= 0, x = a}`

is empty if and only if there exists `f`

in the preorder generated by `-f1, -f2`

and `h`

in the ideal generated by `x - a`

such that `f + h = -1`

.

You can formulate it as a Sum-of-Squares program with

`s_1 * f1 + s_2 * f2 + s_3 * f1 * f2 + p * (x - a) = -1`

where `s_i`

are SOS polynomials and `p`

is any polynomial.

To solve it numerically, you will need to fix the maximal degree of the expression but you can build a hierarchy of SOS programs with higher and higher maximal degree.

An issue with this is that if the coefficients of `p`

are JuMP decision variables and `a`

is also a JuMP decision variable, then `p * (x - a)`

will be bilinear in the JuMP decision variables. Since you have SDP constraints (to constraint the `s_i`

to be SOS) and bilinear constraint, it becomes tricky to solve. You will need to mix an SDP solver with a meta-solver like https://github.com/lanl-ansi/Alpine.jl to handle bilinearity.

However, if you fix `a`

, e.g. you are commuting `g(a)`

and use a separate method that find the root of `g`

given an oracle for `g(a)`

then you don’t have any bilinearity issues.

Here is how it would look like. First define your polynomials:

```
using DynamicPolynomials
@polyvar x y
f1 = -(1 - x)* y + 1/2* (1 - x - x^2 - x^3 - 2* (1 - x)* y)
f2 = (1 - 2*x)*(1 + x + x^2) - 2*(-1 + x)^2*(1 + x)*y
```

You want `s * p`

to have degree `d`

, hence `s`

to have degree `d - deg(p)`

. As `s`

is SOS, it has even degree. Taking all this into account, you get something like this to create `s`

:

```
function s(model, p, d)
deg = div(d - maxdegree(p, y) + 1, 2)
monos = monomials([y], 0:deg)
return @variable(model, variable_type=SOSPoly(monos))
end
```

Now you can combine all this to create your oracle for `g`

.

```
function g(a, d, solver = CSDP.Optimizer)
f_1 = subs(f1, x => a)
f_2 = subs(f2, x => a)
model = Model(solver)
set_silent(model)
p = sum(s(model, f, d) * f for f in [f_1, f_2, f_1 * f_2])
@constraint(model, p == -1)
optimize!(model)
status = termination_status(model)
if status == MOI.OPTIMAL
return false
elseif status == MOI.INFEASIBLE
return true
else
error("Termination status: $status.")
end
end
```

To test it, suppose you want to evaluate `g`

at the given values of `x`

: `x_values = range(0, stop=1, length=10)`

.

With degree `d`

, you get `g.(x_values, 2)`

which gives `true`

for all values except for `x = 1`

.

Climbing up the hierarchy, with `d = 6`

: `g.(x_values, 6)`

, you get the same except that you also get `false`

for an additional values before `x = 1`

.

With `d = 10`

, you get `false`

for the last three values of `x`

.

Accepting the status `MOI.ALMOST_OPTIMAL`

as `false`

, with `d = 20`

, you get `false`

for the last 4 values of `g`

.

**EDIT**: I just realized that in the implementation, I forgot we need to take `-f1`

and `-f2`

insted of `f1`

and `f2`

as mentioned in the beginning.