Existence of solution for a system of multi-variate polynomial equations and in-equations
The question seems to be based on false premises, so let me try to deconfuse it.


Solvability of systems of polynomial equations with integer coefficients is NP-hard over any fixed field (or integral domain). For instance, it is straightforward to reduce the satisfiability of a 3-CNF $\phi(x_1,\dots,x_k)$ to the solvability of a polynomial system including $\{x_i^2-x_i:i=1,\dots,k\}$ and a bunch of polynomials translating each 3-clause.
Solvability of general polynomial systems with integer (or rational) coefficients is polynomially reducible to solvability of systems of multilinear degree-$2$ polynomials with coefficients in $\{-1,0,1\}$ and at most $3$ nonzero monomials each. To see this, express each polynomial by an arithmetical circuit (aka straightline program) using operations $+$, $\cdot$, variables, and constants $0,1,-1$. Introduce a new variable for each node in the circuit, and take the set of polynomials describing the value of each node in terms of its inputs.
The multivariate polynomial ring $K[x_1,\dots,x_k]$ (for a field $K$) is indeed a GCD domain, but for $k\ge2$ it is not a Bézout domain, hence the gcd’s do not behave the way you think they do. For example, the polynomials $x_1$ and $x_2$ have a common root $(0,0)$, even though $\gcd(x_1,x_2)=1$.


Thus, your problem is NP-hard, and as such it has no known polynomial-time algorithm (deterministic or randomized). Under plausible assumptions, it even has no subexponential-time randomized algorithm.

Furthermore:


Over a field, inequations can be replaced by equations: instead of $q(x_1,\dots,x_k)\ne0$, use $q(x_1,\dots,x_k)y=1$, where $y$ is a fresh variable.
Solvability of systems of polynomial equations (hence also inequations, by above) in $\mathbb C$ can be ckecked in PSPACE, due to Canny [1]. This is the best known unconditional result, but:
Assuming the Riemann hypothesis for Dedekind $\zeta$ functions, solvability of polynomial systems in $\mathbb C$ can be checked in the class AM, due to Koiran [2]. Moreover, we can derandomize AM to NP under plausible Nisan–Widgerson-style complexity assumptions, thus assuming both, the problem becomes NP-complete.


[1] John Canny, Some algebraic and geometric computations in PSPACE, Proc. 20th Annual ACM Symposium on Theory of Computing, 1988, pp. 460–467.

[2] Pascal Koiran, Hilbert’s Nullstellensatz is in the Polynomial Hierarchy, Journal of Complexity 12 (1996), no. 4, pp. 273–286.