**NOTE :** Incomplete answer : solution framework, suitable for numerical computation (any precision), but no formal solution.

If two quantities are equal, their squares are also equal. You can therefore get *candidate* solutions to `eqn`

by solving `eqn.lhs()^2==eqn.rhs()^2`

, which van be abbreviated as `eqn^2`

.

But by doing this, you may introduce *spurious* solutions : `x==y`

implies that `x^2==y^2`

, but so does `x==-y`

. So you must filter your set of candidate solutions by checking against `eqn`

.

In your case, your candidate solutions are :

```
sage: (eqn^2).solve(p)
[p == 1/2*(2*h*n - sqrt(-4*(h^2 - h)*n + 1) + 1)/(n + 1), p == 1/2*(2*h*n + sqrt(-4*(h^2 - h)*n + 1) + 1)/(n + 1)]
```

$$ \left[p = \frac{2 \, h n - \sqrt{-4 \, {\left(h^{2} - h\right)} n + 1} + 1}{2 \, {\left(n + 1\right)}}, \quad p = \frac{2 \, h n + \sqrt{-4 \, {\left(h^{2} - h\right)} n + 1} + 1}{2 \, {\left(n + 1\right)}}\right] $$

but *formally* checking these two roots isn't specially easy (numerically, it's a piece of cake...) : it amounts to check that :

```
sage: [eqn.subs(u) for u in (eqn^2).solve(p)]
[-h + 1/2*(2*h*n - sqrt(-4*(h^2 - h)*n + 1) + 1)/(n + 1) == 1/2*sqrt(-(2*h*n - sqrt(-4*(h^2 - h)*n + 1) + 1)*((2*h*n - sqrt(-4*(h^2 - h)*n + 1) + 1)/(n + 1) - 2)/((n + 1)*n)),
-h + 1/2*(2*h*n + sqrt(-4*(h^2 - h)*n + 1) + 1)/(n + 1) == 1/2*sqrt(-(2*h*n + sqrt(-4*(h^2 - h)*n + 1) + 1)*((2*h*n + sqrt(-4*(h^2 - h)*n + 1) + 1)/(n + 1) - 2)/((n + 1)*n))]
```

$$ -h + \frac{2 \, h n - \sqrt{-4 \, {\left(h^{2} - h\right)} n + 1} + 1}{2 \, {\left(n + 1\right)}} = \frac{1}{2} \, \sqrt{-\frac{{\left(2 \, h n - \sqrt{-4 \, {\left(h^{2} - h\right)} n + 1} + 1\right)} {\left(\frac{2 \, h n - \sqrt{-4 \, {\left(h^{2} - h\right)} n + 1} + 1}{n + 1} - 2\right)}}{{\left(n + 1\right)} n}} $$

$$ -h + \frac{2 \, h n + \sqrt{-4 \, {\left(h^{2} - h\right)} n + 1} + 1}{2 \, {\left(n + 1\right)}} = \frac{1}{2} \, \sqrt{-\frac{{\left(2 \, h n + \sqrt{-4 \, {\left(h^{2} - h\right)} n + 1} + 1\right)} {\left(\frac{2 \, h n + \sqrt{-4 \, {\left(h^{2} - h\right)} n + 1} + 1}{n + 1} - 2\right)}}{{\left(n + 1\right)} n}} $$

which is *more complex* than the original problem.

Nevertheless, HTH,

**EDIT :** A bit of numerical exploration *suggests* that for real `h`

and `p`

, the second solution is the only one valid :

```
sage: implicit_plot((h-p+sqrt(p*(1-p)/n)).subs(n==10),(p,-2,2),(h,-2,2))
Launched png viewer for Graphics object consisting of 1 graphics primitive
```

**EDIT 2:** `sympy`

ofers a tool (`unrad`

) which can eliminate radicals from expressions, returning an expression whose set of roots contains those of the original expression, as well as a possiblke set of change of variables.

Let's try it. Put our equation in an "expression" form (preferred by sympy) :

```
sage: E0 = eqn.lhs()-eqn.rhs()
```

Eliminate the radical from our equation

```
sage: from sympy.solvers.solvers import unrad
sage: E1 = unrad(E0._sympy_()) ; E1
(h**2*n - 2*h*n*p + n*p**2 + p**2 - p, [])
```

Solve it

```
sage: Cand = E1[0]._sage_().solve(p) ; Cand
[p == 1/2*(2*h*n - sqrt(-4*(h^2 - h)*n + 1) + 1)/(n + 1), p == 1/2*(2*h*n + sqrt(-4*(h^2 - h)*n + 1) + 1)/(n + 1)]
```

This set of roots contains the roots of E0, if any. But substituting `p`

in E0 gives expressions containing radicals again :

```
sage: [E0.subs(u) for u in Cand]
[-h - 1/2*sqrt(-(2*h*n - sqrt(-4*(h^2 - h)*n + 1) + 1)*((2*h*n - sqrt(-4*(h^2 - h)*n + 1) + 1)/(n + 1) - 2)/((n + 1)*n)) + 1/2*(2*h*n - sqrt(-4*(h^2 - h)*n + 1) + 1)/(n + 1),
-h - 1/2*sqrt(-(2*h*n + sqrt(-4*(h^2 - h)*n + 1) + 1)*((2*h*n + sqrt(-4*(h^2 - h)*n + 1) + 1)/(n + 1) - 2)/((n + 1)*n)) + 1/2*(2*h*n + sqrt(-4*(h^2 - h)*n + 1) + 1)/(n + 1)]
```

**EDIT :** what follows between lines is wrong, *wrong*, **wrong**. See after the second line for an explanation :

================================================================================

So eliminate them the same way :

```
sage: [unrad(E0.subs(u)._sympy_()) for u in Cand]
[(0, []), (0, [])]
```

It turns out that *both* solutions are also solutions of our original equation.

An interesting question is "why `implicit_plot`

suggests only one real solution ?", for which I have not (yet) an answer.

================================================================================

**EDIT :** My error above is to re-use `unrad`

, which may re-introduce spurious answers (and *does*, in the present case).

A numerical check is (again) in order :

```
sage: [E0.subs(u).subs({h:1/2, n:10}).numerical_approx() for u in Sol]
[-0.301511344577764, 0.000000000000000]
```

Only one candidate solution can possibly a root of the original system.

HTH,