Japanese theorem for cyclic quadrilaterals Proof Inversion

Two weeks ago our professor taught us the Japanese theorem for cyclic quadrilaterals. It states that the inscribed centres of the four triangles formed by two sides and a diagonal of a cyclic quadrilateral always form a rectangle.

After the proof he mentioned that he was convinced that the reverse also holds, that is, that whenever the incircle centres form a rectangle the outer quadrilateral is a cyclic quadrilateral, but he never managed to prove it. After finding absolutely nothing on the internet about the inversion of the Japanese theorem, I set out to prove it myself.

The big problem is that the construction of the rectangle is not unique, many different chordal quadrilaterals have the same rectangle of incircle centres. Since there is thus no construction rule in the opposite direction, all my elementary geometric attempts have failed.

So I came up with the idea of a purely arithmetic proof for the Japanese theorem. In the hope that, unlike the geometric proof, it works in both directions. Without limiting generality, I have given the four vertices of the chord quadrilateral the following coordinates: $A (0, 0), B (1, 0), C(x_c, y_c), D(x_d, y_d)$. This is possible because by moving and scaling in the coordinate system alone, any quadrilateral can be brought into this shape. The simple values for A and B simplify the calculation enormously. After much effort I have found formulas to represent all vectors between the incircle centres only dependent from $x_c, y_c, x_d, y_d$:

$$\overrightarrow{M_1M_2} = \left(\begin{array}{c} \frac{\sqrt{x_c^2 + y_c^2} + x_c}{\sqrt{1 + x_c^2 - 2x_c + y_c^2} + \sqrt{x_c^2 + y_c^2} + 1} - \frac{\sqrt{x_d^2 + y_d^2}+ x_d}{1 + \sqrt{1 + x_d^2 - 2x_d + y_d^2} + \sqrt{x_d^2 + y_d^2}}\\ \frac{y_c}{\sqrt{1 + x_c^2 - 2x_c + y_c^2} + \sqrt{x_c^2 + y_c^2} + 1} - \frac{y_d}{1 + \sqrt{1 + x_d^2 - 2x_d + y_d^2} + \sqrt{x_d^2 + y_d^2}} \end{array}\right)$$

$$\overrightarrow{M_2M_3} = \left(\begin{array}{c} \frac{\sqrt{x_c^2 + y_c^2 + x_d^2 + y_d^2 - 2x_cx_d - 2y_cy_d} + \sqrt{1 + x_d^2 - 2x_d + y_d^2} \cdot x_c + \sqrt{1 + x_c^2 - 2x_c + y_c^2} \cdot x_d}{\sqrt{x_c^2 + y_c^2 + x_d^2 + y_d^2 - 2x_cx_d - 2y_cy_d} + \sqrt{1 + x_d^2 - 2x_d + y_d^2} + \sqrt{1 + x_c^2 - 2x_c + y_c^2}} - \frac{\sqrt{x_c^2 + y_c^2} + x_c}{\sqrt{1 + x_c^2 - 2x_c + y_c^2} + \sqrt{x_c^2 + y_c^2} + 1}\\ \frac{\sqrt{1 + x_d^2 - 2x_d + y_d^2} \cdot y_c + \sqrt{1 + x_c^2 - 2x_c + y_c^2} \cdot y_d}{\sqrt{x_c^2 + y_c^2 + x_d^2 + y_d^2 - 2x_cx_d - 2y_cy_d} + \sqrt{1 + x_d^2 - 2x_d + y_d^2} + \sqrt{1 + x_c^2 - 2x_c + y_c^2}} - \frac{y_c}{\sqrt{1 + x_c^2 - 2x_c + y_c^2} + \sqrt{x_c^2 + y_c^2} + 1} \end{array}\right)$$

$$\overrightarrow{M_3M_4} = \left(\begin{array}{c} \frac{\sqrt{x_d^2 + y_d^2} \cdot x_c + \sqrt{x_c^2 + y_c^2} \cdot x_d}{\sqrt{x_c^2 + y_c^2 + x_d^2 + y_d^2 - 2x_cx_d - 2y_cy_d} + \sqrt{x_d^2 + y_d^2} + \sqrt{x_c^2 + y_c^2}} - \frac{\sqrt{x_c^2 + y_c^2 + x_d^2 + y_d^2 - 2x_cx_d - 2y_cy_d} + \sqrt{1 + x_d^2 - 2x_d + y_d^2} \cdot x_c + \sqrt{1 + x_c^2 - 2x_c + y_c^2} \cdot x_d}{\sqrt{x_c^2 + y_c^2 + x_d^2 + y_d^2 - 2x_cx_d - 2y_cy_d} + \sqrt{1 + x_d^2 - 2x_d + y_d^2} + \sqrt{1 + x_c^2 - 2x_c + y_c^2}}\\ \frac{\sqrt{x_d^2 + y_d^2} \cdot y_c + \sqrt{x_c^2 + y_c^2} \cdot y_d}{\sqrt{x_c^2 + y_c^2 + x_d^2 + y_d^2 - 2x_cx_d - 2y_cy_d} + \sqrt{x_d^2 + y_d^2} + \sqrt{x_c^2 + y_c^2}} - \frac{\sqrt{1 + x_d^2 - 2x_d + y_d^2} \cdot y_c + \sqrt{1 + x_c^2 - 2x_c + y_c^2} \cdot y_d}{\sqrt{x_c^2 + y_c^2 + x_d^2 + y_d^2 - 2x_cx_d - 2y_cy_d} + \sqrt{1 + x_d^2 - 2x_d + y_d^2} + \sqrt{1 + x_c^2 - 2x_c + y_c^2}} \end{array}\right)$$

$$\overrightarrow{M_4M_1} = \left(\begin{array}{c} \frac{\sqrt{x_d^2 + y_d^2}+ x_d}{1 + \sqrt{1 + x_d^2 - 2x_d + y_d^2} + \sqrt{x_d^2 + y_d^2}} - \frac{\sqrt{x_d^2 + y_d^2} \cdot x_c + \sqrt{x_c^2 + y_c^2} \cdot x_d}{\sqrt{x_c^2 + y_c^2 + x_d^2 + y_d^2 - 2x_cx_d - 2y_cy_d} + \sqrt{x_d^2 + y_d^2} + \sqrt{x_c^2 + y_c^2}}\\ \frac{y_d}{1 + \sqrt{1 + x_d^2 - 2x_d + y_d^2} + \sqrt{x_d^2 + y_d^2}} - \frac{\sqrt{x_d^2 + y_d^2} \cdot y_c + \sqrt{x_c^2 + y_c^2} \cdot y_d}{\sqrt{x_c^2 + y_c^2 + x_d^2 + y_d^2 - 2x_cx_d - 2y_cy_d} + \sqrt{x_d^2 + y_d^2} + \sqrt{x_c^2 + y_c^2}} \end{array}\right)$$

As I feared, the formulae that result from the geometric construction are so complicated that I could not work with them even with the help of Mathematica and Matlab. So I discarded the arithmetic approach. After many more failed attempts, I wanted to rule out the possibility that a counterexample exists and that the reverse of the Japanese theorem is not true in the first place. I have written a programme in C++ which tries out millions and millions of constellations and indeed all points where the centres of the incircles form a rectangle always lie on a circle.

But I've run out of ideas how I could still attempt the proof. And apparently no one has ever cared about the validity of the reversion, I've only found one old StackExchange post. But it is very sparse and the counterexample presented in the only answer is invalid because the vertices of the rectangle do not coincide with the inscribed centres.

Easy angle chasing shows that lines $CM_3, CM_4$ are isogonal in angle $DCM_2$ and lines $DM_3, DM_4$ are isogonal in angle $M_1DC$. Hence, letting $X$ be the intersection of $CM_2$ and $DM_1$, we get that $M_3$ and $M_4$ are isogonal conjugates in triangle $CDX$. Hence $\angle M_2XM_3 = \angle M_4XM_1$.

Since $M_1M_2M_3M_4$ is a rectangle and $\angle M_2XM_3 = \angle M_4XM_1$, an easy symmetry argument shows that $X$ lies on the perpendicular bisector $\ell$ of $M_1M_2$ (which, of course, coincides with the perpendicular bisector of $M_3M_4$).

Similarly, letting $Y$ be the intersection of $AM_4$ and $BM_3$, we get that $Y$ lies on $\ell$.

Let $P$ be the intersection of $AC$ and $BD$. Let $Z$ be the intersection of lines $BY, CX$. Let $T$ be the intersection of lines $DX, AY$. Then $Z$ and $T$ are incenters of $BCP$ and $DAP$, respectively. Again, a symmetry argument shows that $Z$ and $T$ are symmetric wrt to $\ell$. In particular $\angle PZX = \angle XTP$. Therefore $\angle CZP = \angle PTD$, so $\frac\pi 2 + \frac 12\angle CBD = \frac \pi 2 + \frac 12 \angle CAD$, hence $\angle CBD = \angle CAD$, hence $ABCD$ is cyclic.