There are four lines A, B, C, D in 3D space, they all go through (0,0,0). There is a point Ma on a line A.
I need to find points Mb, Mc, Md on lines B, C, D, such that Ma, Mb, Mc, Md forms a square in space (equal edge lengths, right angles ...). We expect that such square exists.
Is there an explicit solution? Can it be computed with an iterative algorithm?
Let $\vec a=M_a$, $\vec b$, $\vec c$, $\vec d$ be four vectors giving the directions of lines $A$, $B$, $C$, $D$. We must find three numbers $r$, $s$, $t$ such that $M_b=r\vec b$, $M_c=s\vec c$, $M_d=t\vec d$ and $M_aM_bM_cM_d$ is a square.
A square is (1) a parallelogram with (2) equal sides and (3) right angles. Condition (1) is satisfied if and only if diagonals $M_aM_c$ and $M_bM_d$ have the same midpoint, that is if and only if: $$ \tag{1} r\vec b + t\vec d - s\vec c = \vec a. $$ We can get rid of two terms on the LHS by taking the scalar product of both sides with $\vec b\times\vec d$ or $\vec b\times\vec c$ or $\vec c\times\vec d$. We can then solve for the remaining unknown: $$ s=-{(\vec b\times\vec d)\cdot\vec a\over(\vec b\times\vec d)\cdot\vec c},\quad r={(\vec c\times\vec d)\cdot\vec a\over(\vec c\times\vec d)\cdot\vec b},\quad t={(\vec b\times\vec c)\cdot\vec a\over(\vec b\times\vec c)\cdot\vec d}. $$ Note that all denominators give the same result (sign apart), for the circular shift property of scalar triple product, and do not vanish if lines $B$, $C$, $D$ are not coplanar (the case when all lines lie in the same plane requires a different solution, see edit below). Hence you can always find a parallelogram but of course there is no guarantee it is also a square, if you take four generic lines. On the other hand the solution is unique, hence if you know it was a square you'll obtain that.
EDIT.
In the case all lines lie in the same plane, then equation $(1)$ is equivalent to 2 scalar equations and its solution gives the values for only two unknowns, for instance $r$ and $t$. We can solve it by the same technique used above, taking the vector product of both sides with $\vec b$ or $\vec d$. We obtain: $$ r={(\vec d\times\vec a)+s(\vec d\times\vec c) \over(\vec d\times\vec b)},\quad t={(\vec b\times\vec a)+s(\vec b\times\vec c) \over(\vec b\times\vec d)}. $$ Parameter $s$ is then available to enforce one of the other conditions, for instance $(r\vec b-\vec a)\cdot(t\vec d-\vec a)=0$ to obtain a rectangle.