The "**short answer**" is as follows.

The point $P$ is in general **not constructible**, for instance, for the triangle with sides $6$, $9$, $13$
constructing $P$ implies we can solve a polynomial equation $\Pi(K)=0$, where
$\Pi=a_7K^7+a_6K^6+a_5K^5+a_4K^4+a_3K^3+a_2K^2 + a_1K+a_0$
is an irreducible polynomial of degree seven with rational coefficients,
and a root $K$ of it is the square of the common distance of $P$ to the vertices.
I will use $K$ for both the indeterminate when writing $\Pi$,
and also the particular root value(s) $K\in\overline{\Bbb Q}$ with corresponding to real positive numbers.
A simpler (particular) case is when $\Delta ABC$ is isosceles.
Then $\Pi$ factorizes as $\Pi=\Pi_2^2\cdot\Pi_3$ with $\deg \Pi_2=2$, $\deg \Pi_3=3$,
$\Pi_2$ without roots in $\Bbb R$.
So the possible $K$ values come from the real positive roots of $\Pi_3$. In explicit cases we can easier check $\Pi_3$
has no rational roots. So the constructibility **fails**.

The point $P$ fails also to be unique **in the plane** - as stated by the OP - in infinitely many cases.
To see which cases lead to a failure,
explicit computations of the coefficients of $\Pi$ in terms of the parameters $a,b,c$ are needed.
It turns out that in "most cases" there is only one sign change in the list of coefficients of $\Pi$,
the first coefficients $a_7$, $\color{blue}{a_6}$, $a_5$, $a_4$, $a_3$ being positive, and $a_2$, $a_1$, $a_0$ being negative,
so Decartes' Rule of Signs insures an unique $K$.
However, in "few cases" the coefficient $\color{blue}{a_6}=\color{blue}{-p^4} + 12 p^2 S^2 E + 12 S^4 E^2 + 192 S^6$
is negative, Decartes' Rule of Signs predicts at most three positive roots, and indeed, we have three such roots $K>0$,
leading to three solutions.
Here $p=abc$, $E=a^2+b^2+c^2$, and $S$ is the area of $\Delta ABC$.
Counterexamples should thus have a "small area" $S$, but "big product" $p$.
It is natural to search then for a counterexample with $c$ being "almost" $a+b$, and this produces them quickly, as checked by computer.

We have however a uniqueness of $P$ **in the interior** of the given triangle.
This, together with the existence, may be seen by geometric arguments as in the other answers.
For the existence of a solution i have a deformation argument.
For existence and uniqueness **in the interior** this answer gives some hints / details for the deformation.
(Making things concrete goes beyond the question of the OP.)

Given this verdict, a / my proof cannot work by geometric means,
so analytic tools are needed, my choice was to use barycentric coordinates.
I will try to follow in presentation a minimal path, however providing the needed details
makes minimal not really compact.

The question raised a big echo, it is indeed an interesting question.
(I'm afraid this answer shows things are not so simple, but please give it a chance.)
To control the computations, i need to use computer algebra support (CAS), my choice of weapons is sage,
it also has a great merit in making the code pretty readable for a mathematician, since the methods used
are named to be easily digested by a mathematician.

As reference for barycentric coordinates i will use

Barycentric Coordinates for the Impatient, Max Schindler, Evan Chen, July 13, 2012

The detailed answer starts now.

Let $a,b,c$ be the lengths of the sides of a general non-degenerated triangle $\Delta ABC$.

Let $x,y,z$ with $1=x+y+z$ be the coordinates of the unknown point $P$. So $P=(x,y,z)$ in notation.

Then $A,B,C;P$ and the points $D=AP\cap BC$, $E=BP\cap CA$, $F=CP\cap AB$.
have explicit barycentric descriptions:
$$
\begin{aligned}
A &= (1,0,0)\ ,\\
B &= (0,1,0)\ ,\\
C &= (0,0,1)\ ,\\[2mm]
P &= (x,y,z)\ ,\\
D &= [0:y:z]=\left(0,\frac {y}{y+z},\frac {z}{y+z}\right)\ ,\\
E &= [x:0:z]=\left(\frac {x}{x+z},0,\frac {z}{x+z}\right)\ ,\\
F &= [x:y:0]=\left(\frac {x}{x+y},\frac {y}{x+y},0\right)\ ,\\[3mm]
&\qquad\text{and corresponding displacement vectors are}\\[3mm]
\overrightarrow{DP} &= D - P = \left(-x,\frac {xy}{y+z},\frac {xz}{y+z}\right) = \frac x{y+z}(-(y+z),y,z) \ ,\\
\overrightarrow{EP} &= E - P = \left(\frac {xy}{x+z},-y,\frac {yz}{x+z}\right) = \frac y{z+x}(x,-(z+x),z) \ ,\\
\overrightarrow{FP} &= F - P = \left(\frac {xz}{x+y},\frac {yz}{x+y},-z\right) = \frac z{x+y}(x,y,-(x+y)) \ ,\\[3mm]
&\qquad\text{and corresponding squared lengths are}\\[3mm]
|DP|^2 &= \frac{x^2}{(y+z)^2}\Big(\ -a^2yz + (b^2z+c^2y)(y+z)\ \Big) = \frac{x^2}{(1-x)^2}\Big(\ Q + b^2z + c^2y\ \Big)\ ,\\
|EP|^2 &= \frac{y^2}{(z+x)^2}\Big(\ -b^2zx + (c^2x+a^2z)(z+x)\ \Big) = \frac{y^2}{(1-y)^2}\Big(\ Q + c^2x + a^2z\ \Big)\ ,\\
|FP|^2 &= \frac{z^2}{(x+y)^2}\Big(\ -c^2xy + (a^2y+b^2x)(x+y)\ \Big) = \frac{z^2}{(1-z)^2}\Big(\ Q + a^2y + b^2x\ \Big)\ ,\\
&\qquad\text{where}\\
Q &= -a^2 yz -b^2 zx -c^2xy\ .
\end{aligned}
$$
Let $K>0$ be the common value of the squared distances from $P$
to each of the points $D,E,F$, i.e. $K=|DP|^2=|EP|^2=|FP|^2$.
Then we have to solve the following system of equations in the unknowns $x,y,z;K,Q\in \Bbb R$, $K>0$:
$$
\tag{$\dagger$}
$$
$$
\left\{
\begin{aligned}
1 &= x+y+z\ ,\\
Q &= -a^2yz -b^2zx -c^2xy\ ,\\[2mm]
K(1-x)^2 &= x^2(Q + b^2 z + c^2 y)\ ,\\
K(1-y)^2 &= y^2(Q + c^2 x + a^2 z)\ ,\\
K(1-z)^2 &= z^2(Q + a^2 y + b^2 x)\ ,
\end{aligned}
\right.
$$
(Solutions were introduced by eliminating denomiators, for instance $A=(1,0,0)$ is now a solution.)

Let $J$ be the ideal generated by the above equations (rewritten in terms vanishing of expressions).
From now on, we will work algebraically with the system $(\dagger)$ above.
(Since the parameters $a,b,c$ appear only through their squares, it may be convenient to
use short hand notations $A=a^2$, $B=b^2$, $C=c^2$ for them.)
(There should be no confusion with the vertices
$A,B,C$ of the given triangle.)

Now i am trying to address the points from the question in terms of this system.

**Constructibility of $P$, special case of the $6$, $9$, $13$ triangle**:

This triangle is used for the purpose of (numerically) searching for knwon (general) centers in a triangle, see also

ETC, Search_6_9_13 .

If $P$ is somehow constructible (starting from some fixed values $a,b,c$), then $K$ is also constructible (by rule and compass constructions),
so starting with $6,9,13\in\Bbb Q$
one should obtain that $K$ is an algebraic number in an extension of degree among $1,2,4,8,16,\dots$ -
however, it turns out that $K$ is the root of an irreducible polynomial $\Pi$ of degree seven. The sage code computing $\Pi$
is postponed.
In this special case:
$$
\tag{$1$}
\Pi =
K^7
+ \frac{2420849677}{40884480}\; K^6
+ \frac{2635885}{2106 }\; K^5
+ \frac{920646335}{82134 }\; K^4
+ \frac{349438600}{9477 }\; K^3
- \frac{10186414000}{123201 }\; K^2
- \frac{8451520000}{9477 }\; K
- \frac{215129600000}{123201}
\ .
$$
So the constructibility fails in this case, so it fails.

Let us observe that there is exactly one sign change in the coefficients of $\Pi$.
This happens "often", as seen in the next section.
Below there will be an other example ($61,61,120$ isosceles triangle) where the
above polynomial splits as a product of a squared quadratic and a cubic polynomial,
the cubic part only has real roots, we have an explicit formula for them, and it is easy to check
there is no rational root.

**The minimal polynomial of $K$** in the general case.

In the code part, a polynomial $\Pi$ of degree seven in $R[K]$ over the ring $R=\Bbb Q[a,b,c]$ is computed, so that
the special value $K=|DP|^2=|EP|^2=|FP|^2$ (for each solution $P$ of our problem) satisfies $\Pi(K)=0$.
Explicitly,
$$
\tag{$2$}
$$
$$
\begin{aligned}
\Pi &=
a_7 K^7 +
a_6 K^6 +
a_5 K^5 +
a_4 K^4 +
a_3 K^3 +
a_2 K^2 +
a_1 K +
a_0\ ,\\[3mm]
&\qquad\text{where}\\[3mm]
a_7 &= 64 S^2 p^2\ ,\\
a_6 &= \color{blue}{-p^4} + 12 p^2 S^2 E + 12 S^4 E^2 + 192 S^6\ ,\\
a_5 &= 2 S^4 \;(E^3 + 13p^2 + 34 S^2 E)\ ,\\
a_4 &= S^4 \;(p^2E + 10 S^2 E^2 + 103 S^4)\ ,\\
a_3 &= S^6 \;(2p^2 + 17 S^2 E)\ ,\\
a_2 &= -\frac 14 S^8 \;(E^2 - 40 S^2)\ ,\\
a_1 &= -E \;S^{10} \ ,\\
a_0 &= -S^{12}\ .
\end{aligned}
$$
The coefficients are homogeneous of degrees $\deg a_7=10$, $\deg a_6=12$, ... , $\deg a_0=24$,
when considering the weights $\deg S=\deg E=\deg K=2$, $\deg p=3$.
If $K$ is considered with $\deg K=2$, then $\Pi$ is homogeneous of degree $24$.

These explicit formulas allow to control the signs of the coefficients of $\Pi=\Pi(K)$.
We have (at least) one sign change obtained when passing from $a_3$ to $a_2$.
To see this note that
$$
E^2 - 40 S^2 > E^2 - 48 S^2 = 4(a^4+b^4+c^4-a^2b^2-b^2c^2-c^2a^2)>0\ .
$$
An other change of signs may occur only around $a_6$.
In case $a_6>0$ there is exactly one change of signs at all,
so by Decartes' rule of signs we have **exactly one real positive root**.
So uniqueness of $K$ is insured in case of $a_6>0$.
Let us analyze this in detail.

**Further relations:**

Let $x,y,z;K$ be the solutions of the system $(\dagger)$. Then
there is a relation joining $K$ with each of the variables $x,y,z$:
$$
\tag{$3$}
$$
$$
\begin{aligned}
0 &= 4S^2 x^3 - 2S^2x^2 +K(b^2+c^2-a^2)x^2 + 8 xK^2 - 2K^2 \ ,\\
0 &= 4S^2 y^3 - 2S^2y^2 +K(c^2+a^2-b^2)y^2 + 8 yK^2 - 2K^2 \ ,\\
0 &= 4S^2 z^3 - 2S^2z^2 +K(a^2+b^2-c^2)z^2 + 8 zK^2 - 2K^2 \ ,
\end{aligned}
$$
and relations for $(y,K)$ and $(z,K)$ can be written analogously.
We also have:
$$
\tag{$4$}
$$
$$
\begin{aligned}
0 &= a^2\;y^2z^2 + 4K\;yz(y+z) - K(y+z)^2\ ,\\
0 &= b^2\;z^2x^2 + 4K\;zx(z+x) - K(z+x)^2\ ,\\
0 &= c^2\;x^2y^2 + 4K\;xy(x+y) - K(x+y)^2\ ,
\end{aligned}
$$
as found in the code section #4.

(These relations can be used maybe in some way to show algebraically the existence
and uniqueness in the interior of $\Delta ABC$
for the point matching a fixed $K$, solution of the above polynomial $\Pi$.
I tried to make it, but it made me tired.)

**Uniqueness** of $P$ **fails in the plane**:

We have uniqueness,
in case there is a unique positive root $K$
of $\Pi$.
For instance, in the "often" case of a coefficient $a_6>0$ of $\Pi$,
we know from above that there exists exactly one $K>0$ so that for a $P$ with $|PD|^2=|PE|^2=|PF|^2$ (in case of its existence)
the common value of the above squared distances is $K$.

However cases can be found where such a $K$ is not unique, and we can check that multiple solutions arise.
In the code section #5 there is the following solution found for the triangle with sides $a=b=61$, $c=120$.
The height $h_c$ is $11$. The solutions are roughly:
$$
\begin{aligned}
x = y & = -1.749464267458191 \dots & z &= 4.498928534916382 \dots & K &= 2449.08331\dots\\
x = y & = -0.5905364374782210\dots & z &= 2.181072874956442 \dots & K &= 575.606545\dots\\
x = y & = 0.1209924404736021\dots & z &= 0.7580151190527959\dots & K &= 69.5250174\dots\\
\end{aligned}
$$
Only the last point lives inside $\Delta ABC$.
$K$ is a root of the polynomial $(121K^3 - 374400K^2 + 196020000 K - 11859210000)$.

(For the point $P$ inside $\Delta ABC$ the corresponding $F$ is the mid point of $AB$,
the height $CF$ has length $11$ and $PF$ is approximatively $\sqrt{69.5250174\dots}$. The other two points
have the same $F$. Using this point as origin,
$$\overrightarrow{FP}=x(\overrightarrow{FA}+\overrightarrow{FB})+z\overrightarrow{FC}
=z\overrightarrow{FC}
$$
in all three listed cases. And indeed, $121\;z^2=|CF|^2 \; z^2=|PF|^2=K$ matches the values for $K$
for each $z$ in the list.

For these $x,y,z$ we can also write down exact polynomials having them as roots.
The three $x$-values are for instance roots of the polynomial
$\displaystyle x^3 + \frac{537}{242} x^2 + \frac 34 x - \frac 18$.
none of them is rational.

Of course, "small deformations" lead also to situations with three solutions.

For the existence, geometric ideas may work better.

**Existence and uniqueness** of $P$ in the **interior** of $\Delta ABC$:

So we restrict the values for the barycentric coordinates $x,y,z$ to positive values,
$x, y, z>0$, $x+y+z=1$.
Consider the following picture:

(We do not know right now that an interior $P^*$ exists so that for the cevians $AD^*$, $BE^*$, $CF^*$ through $P^*$ we have $|P^*D^*|=|P^*E^*|=|P^*F^*|$.
Also its uniqueness is still an issue.)
(Also the other answers suggest a possibility to "get closer to the / a solution $P^*$ - if it exists.)

I propose to use two geometrical schemes to "get" (a) $P^*$ in the limit.

*First scheme:*

Start with a point $P=P_0$. Compare the distances from $P$ to $D,E,F$ and take the maximal one, say it is $PE$.
Let a point $Q=P_1$ slide from $P$ to $E$.
Then we stop when
$f_P(Q):=2|QE_1|-|QD_1|-|QF_1|$
vanishes first $f_P(Q)=0$.
Then iterate.
For the convergence we need some argument. (Is this procedure / function $P\to Q$ a contraction?)

*Second scheme:*

Start with an interior point $P=(x,y,z)$, and consider the squared distances $PD^2$, $PE^2$, $PF^2$.
Then consider the point $Q=[X:Y:Z]$ where $X,Y,Z$ are weighted versions of $(x,y,z)$ as follows:
$$
\begin{aligned}
X &=x\cdot \frac{PE^2+PF^2}{2(PD^2+PE^2+PF^2)}\ ,\\
Y &=y\cdot \frac{PF^2+PD^2}{2(PD^2+PE^2+PF^2)}\ ,\\
Z &=z\cdot \frac{PD^2+PE^2}{2(PD^2+PE^2+PF^2)}\ .
\end{aligned}
$$
(Intuition: If for instance $PE^2$ is maximal among $PD^2$, $PE^2$, $PF^2$, then we would like to move $P$ along $PE$ towards $E$.
This corresponds to making $y$ smaller. We can try to multiply $y$ with some "controlled" subunitary factor - and
$\frac{PF^2+PD^2}{2(PD^2+PE^2+PF^2)}$ is such a factor.
At the end, the point $[X:Y:Z]$ has to be normed again, so the result is $Q=\frac 1{X+Y+Z}(X,Y,Z)$.
Of course, this is only an Ansatz, one has to show it works.
I only have numerical support so far - see the code section #6.)

*Third scheme:*
This is less explicit, but we do not have issues as the ones described above.

Let $K\ge0$ be a parameter, we let it variate from $0$ to $\infty$. For each such $K$ we draw three curves.
The $D$-curve is the locus of all points $P$ obtained as follows. Let $D$ slide along the line $BC$.
Draw the cevian $AD$. Let $P$ be on the ray $[DA$ such that $PD^2=K$. Then the $D$-curve is the locus of the points $P$ as $D$ runs on $BC$.
The $D$-curve has the line $BC$ as double asymptote, it is symmetric w.r.t. the height from $A$ and for $K\to 0_+$ it tends to $BC$.
Use as orientation for it the direction that in the limit corresponds to the direction of $BC$ from $B$ to $C$.

Similarly consider the $E$-curve, and the $F$-curve.
They are oriented corresponding to the directions that in the limit give the direction of $CA$ from $C$ to $A$, and of $AB$ from $A$ to $B$.

For $K\to 0_+$, in the limit $K=0$, the three curves become $BC$, $CA$, $AB$, and the oriented area between these curves
is $S>0$. Now let $K$ grow. At some point, e.g. when $K$ is greater than the square of the biggest height, the oriented area becomes negative.
So at some point $K^*$ this area is zero. This corresponds to the the case when the three curves are passing through one point.
(I do not have a simple argument to show this is an interior point.)
This elucidates the existence.

$\square$

A rough picture for the three curves at their intersection in the zoom of our objective is as follows:

I hope it is clear how the three "lens" of the zoom are moving. To fix idea, assume $\Delta ABC$ has all angles $<90^\circ$,
so the orthocenter $H$ is in its interior. Let $AA'$, $BB'$, $CC'$ be the heights, intersecting in $H$.

- For $K=0$ the area between the three curves is $\Delta ABC$.
- Arrange that $HA'\le HB'\le HC'$ by permuting $A,B,C$.
- when $K=A'H^2$, we still have a positive area between the three curves.
- when $K=C'H^2$, we have a negative area between the three curves.
- so the value $K^*$ is squared between these values.

In case of an obtuse angle - say in $A$ - use an alternative zooming objective
with mobile lens for the $E$- and $F$-curves and keep the third lens constant to be the line $BC$ (instead of the $D$-curve)
to see that the intersection of the three curves still has to be an interior point.

Computer algebra support.

**Code #1:** The value of $K=|PD|^2 =|PD|^2 =|PD|^2$ in the case of the $6,9,13$ triangle.

```
a, b, c = 6, 9, 13
R.<x,y,z,K> = PolynomialRing(QQ)
def eq(a, b, c, x, y, z, K):
Q = - a^2*y*z - b^2*z*x - c^2*x*y
return x^2 * ( Q + b^2*z + c^2*y ) - K*(1 - x)^2
J = R.ideal([ x + y + z - 1,
eq(a, b, c, x, y, z, K),
eq(b, c, a, y, z, x, K),
eq(c, a, b, z, x, y, K), ])
JK = J.elimination_ideal([x, y, z])
print("Generator(s) of the elimination ideal JK after eliminating x, y, z:")
for g in JK.groebner_basis():
print(f'{g}\n')
```

And we obtain:

```
Generator(s) of the elimination ideal JK after eliminating x, y, z:
K^9 + 2420849677/40884480*K^8 + 2635885/2106*K^7 + 920646335/82134*K^6 + 349438600/9477*K^5
- 10186414000/123201*K^4 - 8451520000/9477*K^3 - 215129600000/123201*K^2
```

It turns out that the above polynomial is $K^2$ times some irreducible polynomial of degree seven.
To have the entry that should be compared with ETC, we ask for the points in the ring of real algebraic numbers:

```
sage: J.variety(ring=AA)
[{K: 0, z: 0, y: 0, x: 1},
{K: 0, z: 0, y: 1, x: 0},
{K: 0, z: 1, y: 0, x: 0},
{K: 3.973344192056688?,
z: 0.5361103522937736?, y: 0.2667643469973961?, x: 0.1971253007088303?}]
```

(New roots were introduced after multiplication with denominators.)
Only the last real point is significant. In ETC we have the value for
$2\cdot\operatorname{Area}(\Delta ABC)\cdot \frac xa$. The area of the $6,9,13$ triangle is $\sqrt{14(14-6)(14-9)(14-13)}=\sqrt{560}$.
So we try to match the following value with the 6-9-13-search:

```
sage: x0 = J.variety(ring=AA)[-1][x]
sage: etc_match = 2 * sqrt(560) * x0/6
sage: etc_match.n(200)
1.5549453416812577768807502448833414277833445395269911652987
```

This was already done in the comment of
Peter Taylor. To have again a convincing statement against constructibility,
the value `x0`

used above is an algebraic number, root of an irreducible polynomial in $\Bbb Q[x]$
of degree seven:

```
sage: x0.minpoly()
x^7 - 28472791/10221120*x^6 + 46797133/15331680*x^5 - 52400213/30663360*x^4
+ 28733/54756*x^3 - 5173/54756*x^2 + 280/13689*x - 35/13689
```

**Code #2:** The value of $K=|PD|^2 =|PD|^2 =|PD|^2$ in the case of an isosceles, general triangle with sides $a,a,b$.
To avoid the factor $K^2$, i will assume $K$ invertible below.

```
R.<a,b, x,y,z, K,K_inv> = PolynomialRing(QQ)
def eq(a, b, c, x, y, z, K):
Q = - a^2*y*z - b^2*z*x - c^2*x*y
return x^2 * ( Q + b^2*z + c^2*y ) - K*(1 - x)^2
J = R.ideal([ x + y + z - 1, K*K_inv - 1,
eq(a, a, b, x, y, z, K),
eq(a, b, a, y, z, x, K),
eq(b, a, a, z, x, y, K), ])
JK = J.elimination_ideal([x, y, z, K_inv])
print("Generator(s) of the elimination ideal JK after eliminating x, y, z:")
for g in JK.groebner_basis():
print(f'{g.factor()}\n')
```

And we obtain:

```
Generator(s) of the elimination ideal JK after eliminating x, y, z:
(1/256) * (16*a^4*b^2 - 8*a^2*b^4 + b^6 + 64*a^4*K + 32*a^2*b^2*K - 12*b^4*K + 256*a^2*K^2)
* (16*a^4*b^4 - 8*a^2*b^6 + b^8 - 32*a^2*b^4*K + 8*b^6*K
- 1024*a^2*b^2*K^2 + 272*b^4*K^2 - 4096*a^2*K^3 + 1024*b^2*K^3)
```

(Result was manually rearranged.)
Even in this case, the constructibility fails. For instance for the sides $2,2,3$ we have:

```
a, b, c = 2, 2, 3
R.<x,y,z, K,K_inv> = PolynomialRing(QQ)
J = R.ideal([ x + y + z - 1, K*K_inv - 1,
eq(a, b, c, x, y, z, K),
eq(b, c, a, y, z, x, K),
eq(c, a, b, z, x, y, K), ])
JK = J.elimination_ideal([x, y, z, K_inv])
print("Generator(s) of the elimination ideal JK after eliminating x, y, z:")
for g in JK.groebner_basis():
print(f'{g.factor()}\n')
```

This gives:

```
Generator(s) of the elimination ideal JK after eliminating x, y, z:
(1/7340032) * (1024*K^2 + 1204*K + 441) * (7168*K^3 + 14832*K^2 + 4536*K - 3969)
```

The quadratic factor has no real roots.
The cubic factor has exactly one positive root,
$$
\frac3{448}
\left( -103 + \sqrt[3]{5(111205 + 6272\sqrt{105}} + \sqrt[3]{5(111205 - 6272\sqrt{105}} \right)
\\ \approx 0.3643855138832992\dots
$$
Constructibility fails.

**Code #3:**

We consider $K=|PD|^2 =|PE|^2 =|PF|^2$ (if the needed $P$ exists) as a function of the sides $a,b,c$
of a triangle. Then there is a polynomial $\Pi\in\Bbb Q[a,b,c]\;[\kappa]$ of degree seven w.r.t. $\kappa$, which vanishes in $K$, $\Pi(a,b,c;K)=0$.
Its coefficients, expressed in terms of the quantities $E = a^2+b^2+c^2$, $S^2=s(s-a)(s-b)(s-c)$ (squared area), and $p=abc$ (product) are obtained as follows:

```
R.<a,b,c, x,y,z, K,K_inverse, p,SS,E> = PolynomialRing(QQ)
J = R.ideal([ x + y + z - 1,
K * K_inverse - 1,
SS - 1/16 * (a+b+c) * (a+b-c) * (b+c-a) * (c+a-b),
E - (a^2 + b^2 + c^2),
p - a*b*c,
eq(a, b, c, x, y, z, K),
eq(b, c, a, y, z, x, K),
eq(c, a, b, z, x, y, K), ])
for g in J.elimination_ideal([x, y, z, K_inverse, a, b, c]).groebner_basis():
for k in range(g.degree(K), 0, -1):
print(f'Coefficient of K^{k} is {g.coefficient(K**k).factor()}')
print(f'Coefficient of K^0 is {g.subs({K : 0}).factor()}')
```

And the coefficient are explicitly:

```
Coefficient of K^7 is (-64) * SS * p^2
Coefficient of K^6 is (-1) * (-p^4 + 12*p^2*SS*E + 12*SS^2*E^2 + 192*SS^3)
Coefficient of K^5 is (-2) * SS^2 * (E^3 + 13*p^2 + 34*SS*E)
Coefficient of K^4 is (-1) * SS^2 * (p^2*E + 10*SS*E^2 + 103*SS^2)
Coefficient of K^3 is (-1) * SS^3 * (2*p^2 + 17*SS*E)
Coefficient of K^2 is (1/4) * SS^4 * (E^2 - 40*SS)
Coefficient of K^1 is E * SS^5
Coefficient of K^0 is SS^6
```

so after changing the sign of each coefficient we obtain the values from $(2)$.

**Code #4:** If we know $K$, can we "easily" obtain $x,y,z$?
In other words, are there any "simple relations" among the variables, that would lead to
a quick determination of $x$, when $a,b,c,K$ are known?

Here is the elimination of variables, all but $a,b,c;x,K$.

```
R.<x,y,z, a,b,c, K,K_inv> = PolynomialRing(QQ)
J = R.ideal([ x + y + z - 1, K*K_inv - 1,
eq(a, b, c, x, y, z, K),
eq(b, c, a, y, z, x, K),
eq(c, a, b, z, x, y, K), ])
JxK = J.elimination_ideal([y, z, K_inv])
print("Generator(s) of the elimination ideal JxK after eliminating y, z:")
for g in JxK.groebner_basis():
print(f'{g}\n')
```

And among the many relations shown there is also the following one:

```
x^3*a^4 - 2*x^3*a^2*b^2 + x^3*b^4 - 2*x^3*a^2*c^2 - 2*x^3*b^2*c^2 + x^3*c^4
- 1/2*x^2*a^4 + x^2*a^2*b^2 - 1/2*x^2*b^4 + x^2*a^2*c^2 + x^2*b^2*c^2 - 1/2*x^2*c^4
+ 4*x^2*a^2*K - 4*x^2*b^2*K - 4*x^2*c^2*K
- 32*x*K^2
+ 8*K^2
```

This leads to the relation $(3)$.

We can also try to eliminate $K,z$ and obtain relations among $y,z$.
Here are such relations:

```
sage: for g in J.elimination_ideal([z, K_inv, a, b]).groebner_basis(): print(g, '\n')
x^2*y^2*c^2 + 4*x^2*y*K + 4*x*y^2*K - x^2*K - 2*x*y*K - y^2*K
sage: for g in J.elimination_ideal([z, K_inv, c]).groebner_basis(): print(g, '\n')
x^2*y^2*a^2 + 2*x*y^3*a^2 - 2*x^3*y*b^2 - x^2*y^2*b^2 - 2*x*y^2*a^2 + 2*x^2*y*b^2
+ 4*x^2*y*K - 4*x*y^2*K - x^2*K + y^2*K + 2*x*K - 2*y*K
y^4*a^2 + 2*x^3*y*b^2 + x^2*y^2*b^2 - 2*y^3*a^2 - 2*x^2*y*b^2 + y^2*a^2
+ 8*x*y^2*K - 8*x*y*K - 5*y^2*K + 6*y*K - K
x^4*b^2 + 2*x^3*y*b^2 + x^2*y^2*b^2 - 2*x^3*b^2 - 2*x^2*y*b^2 + x^2*b^2
+ 4*x^2*y*K + 4*x*y^2*K - 4*x^2*K - 8*x*y*K - y^2*K + 4*x*K + 2*y*K - K
```

**Code #5:**

Let us see what happens in the case of an isosceles triangle with sides $a=61$, $b=61$, $c=120$.
Its height is $11$, since $61^2-60^2=9^2$.

```
a, b, c = 61, 61, 120
R.<x,y,z, K,K_inv> = PolynomialRing(QQ)
def eq(a, b, c, x, y, z, K):
Q = - a^2*y*z - b^2*z*x - c^2*x*y
return x^2 * ( Q + b^2*z + c^2*y ) - K*(1 - x)^2
J = R.ideal([ x + y + z - 1, K*K_inv - 1,
eq(a, b, c, x, y, z, K),
eq(b, c, a, y, z, x, K),
eq(c, a, b, z, x, y, K), ])
points = J.variety(ring=AA)
print(f'{a} {b} {c}')
for dic in points:
print(f'x = {dic[x]} y = {dic[y]} z = {dic[z]} K = {dic[K]}')
g = J.elimination_ideal([x, y, z, K_inv]).groebner_basis()[0]
print(f'K is a root of the following polynomial:\n{g.factor()}')
```

And we obtain:

```
61 61 120
x = -1.749464267458191? y = -1.749464267458191? z = 4.498928534916382? K = 2449.083313436469?
x = -0.5905364374782210? y = -0.5905364374782210? z = 2.181072874956442? K = 575.6065451903618?
x = 0.1209924404736021? y = 0.1209924404736021? z = 0.7580151190527959? K = 69.52501740622754?
K is a root of the following polynomial:
(1/1800964) * (14884*K^2 + 1757041*K + 52707600) * (121*K^3 - 374400*K^2 + 196020000*K - 11859210000)
```

A note on the degree of the above polynomial. It is five, not seven.
However, using the formulas for the coefficients in this special case,

```
a, b, c = 61, 61, 120
s, p, E = (a + b + c)/2, a*b*c, a^2 + b^2 + c^2
SS = s*(s - a)*(s - b)*(s - c)
a7 = 64 * SS * p^2
a6 = -p^4 + 12*p^2*SS*E + 12*SS^2*E^2 + 192*SS^3
a5 = 2 * SS^2 * (E^3 + 13*p^2 + 34*SS*E)
a4 = SS^2 * (p^2*E + 10*SS*E^2 + 103*SS^2)
a3 = SS^3 * (2*p^2 + 17*SS*E)
a2 = (-1/4) * SS^4 * (E^2 - 40*SS)
a1 = - E * SS^5
a0 = - SS^6
var('K')
PI = a7*K^7 + a6*K^6 + a5*K^5 + a4*K^4 + a3*K^3 + a2*K^2 + a1*K + a0
print(PI.factor())
print(PI.roots(ring=AA, multiplicities=False))
```

we obtain

```
207360000 * (121*K^3 - 374400*K^2 + 196020000*K - 11859210000) * (14884*K^2 + 1757041*K + 52707600)^2
[69.52501740622754?, 575.6065451903618?, 2449.083313436469?]
```

So the quadratic factor appears squared.

**Code #6:**

Let us implement the recursion described in the theoretical section.

```
def recursion(a, b, c, bit_precision=40):
IR = RealField(bit_precision) # IR is "real field" - precision is given as argument
a, b, c = IR(a), IR(b), IR(c) # pass to the numerical world with the given sides a, b, c
def d2(P, Q):
x1, y1, z1 = P
x2, y2, z2 = Q
x, y, z = x1 - x2, y1 - y2, z1 - z2
return - a^2*y*z - b^2*z*x - c^2*x*y
P = vector(IR, 3, [ IR(1/3), IR(1/3), IR(1/3) ]) # start recursion in centroid
for k in range(10):
print(P)
x, y, z = P
D = vector(IR, 3, [ 0, y/(y+z), z/(y+z)])
E = vector(IR, 3, [ x/(x+z), 0, z/(x+z)])
F = vector(IR, 3, [ x/(x+y), y/(x+y), 0])
PD2, PE2, PF2 = d2(P, D), d2(P, E), d2(P, F)
d2sum = PD2 + PE2 + PF2
X = x * (PE2 + PF2) / 2. / d2sum
Y = y * (PF2 + PD2) / 2. / d2sum
Z = z * (PD2 + PE2) / 2. / d2sum
X, Y, Z = X/(X+Y+Z), Y/(X+Y+Z), Z/(X+Y+Z)
P = vector(IR, 3, [X, Y, Z])
```

Then calling `recursion(6, 9, 13)`

we return the following result:

```
sage: recursion(6, 9, 13)
(0.33333333333, 0.33333333333, 0.33333333333)
(0.22960372960, 0.30827505828, 0.46212121212)
(0.20266594188, 0.26478637847, 0.53254767965)
(0.19806794747, 0.26656776311, 0.53536428942)
(0.19727171883, 0.26670350678, 0.53602477439)
(0.19715195703, 0.26675849714, 0.53608954584)
(0.19712950067, 0.26676259568, 0.53610790365)
(0.19712606659, 0.26676417869, 0.53610975472)
(0.19712542143, 0.26676429665, 0.53611028192)
(0.19712532272, 0.26676434216, 0.53611033512)
```

This "seems to converge" to the algebraic exact value found by code #1.

2more comments