In some sense, the reason you are running into these 'problems' is that you are working with the ring of real polynomials rather than the ring of real-valued analytic functions, which is also a UFD.  For example, it is not hard to show that there is a (unique) real-valued, real-analytic function $f$ defined in a neighborhood of $0$ and satisfying $f(0)=\frac12$ such that
$$
y^3 + 2x^2y-x^4 = \bigl(y - x^2f(x^2)\bigr)\bigl(y^2 + x^2f(x^2) y + x^2/f(x^2)\bigr),
$$
so $y^3 + 2x^2y-x^4$ is reducible in the ring of real-valued, real analytic functions defined on a neighborhood of the origin.  The curve $y = x^2f(x^2)$ is smooth (in fact, real-analytic, of course), but the $y$-discriminant of the quadratic factor is $x^4f(x^2)^2-4x^2/f(x^2) = -8x^2 + \cdots$, so the only real point of $y^2 + x^2f(x^2) y + x^2/f(x^2)=0$ near the origin is the origin itself.  (The quadratic factor is irreducible in the ring of real-valued analytic functions defined on a neighborhood of the origin.)

Thus, a more easily approached question might be:  **Suppose that the origin is a singular zero of an *irreducible* element $f$ in the ring of real-valued analytic functions defined on a neighborhood of the origin.  Can the zero locus of $f$ be a nonsingular real-analytic hypersurface near the origin?**

The answer to this question is 'no' because near a nonsingular point of such a hypersurface, there will always be a real-analytic function $g$ that vanishes on it whose differential at that point is nonvanishing.  However, this will imply that $g$ is a factor of $f$, which is assumed to be irreducible.

**Added Remark:**  The case of a curve is not hard, using some standard facts about resolution of curve singularities:  If $f(x,y)$ is a nonzero, real-valued analytic function defined on a neighborhood of the origin in $\mathbb{R}^2$ that is irreducible in this ring and satisfies $f_x(0,0)=f_y(0,0)=0$, then the locus $f(x,y)=0$ cannot be a smooth curve in a neighborhood of the origin.  A sketch of a proof is as follows:  If the origin is not isolated, then $f(z,w)$ is a $\mathbb{C}$-valued analytic function defined on a neighborhood of the origin in $\mathbb{C}^2$ that is also irreducible in this larger ring, and hence there is a neighborhood of the origin in $\mathbb{C}^2$ such that, in this neighborhood, the locus $f(z,w)=0$ can be parametrized by an embedded disc in $\mathbb{C}$ in the form $(z,w) = (a(\tau),b(\tau))$ where $a$ and $b$ are analytic functions of $\tau$ for $|\tau| < 1$.  We can assume that neither $a$ nor $b$ is constant, and so we can reparametrize so that $a(\tau) = \tau^k$ for some $k>1$.  In particular, the real locus will be parametrized by some curves of the form $\tau = \omega t$ where $t$ is real and $\omega^k = \pm 1$.  By replacing $t$ by $t/\omega$, we can assume that $(a(t),b(t))$ is real for all small real $t$, and that this parametrizes a 'branch' of the real locus.  In particular, the coefficients of $b$ are real, so our curve is parametrized in the form
$$
(x,y) = \bigl(\ t^k,\ b_l t^l + b_{l+1} t^{l+1} + \cdots\ \bigr)
$$
where, because of the embeddedness property of the disk, the greatest common divisor of $k$ and those $l$ for which $b_l\not=0$ must be $1$.  By a (real) rotation, we can assume that $l>k$, so the curve is expressed in the form
$$
 y = b_l\ x^{l/k} + b_{l+1}\ x^{(l+1)/k} + \cdots
$$
where at least one of the exponents in this series is not an integer.  It follows that the function on the right hand side of this equation cannot be smooth at $x=0$, even though, since $l>k$, it is $C^1$.

Probably, there is a similar argument in more variables that somehow uses either the Mather Division Theorem or something like it that will prove what you want.