We have the Euler's Identity:
$$ \begin{align} n & = (a_1^2 + a_2^2 + a_3^2 + a_4^2) (b_1^2 + b_2^2 + b_3^2 + b_4^2) \newline & = (a_1 b_1-a_2 b_2-a_3 b_3-a_4 b_4)^2 \newline & + (a_1 b_2+a_2 b_1+a_3 b_4-a_4 b_3)^2 \newline & + (a_1 b_3-a_2 b_4+a_3 b_1+a_4 b_2)^2 \newline & + (a_1 b_4+a_2 b_3-a_3 b_2+a_4 b_1)^2 \end{align} \tag{1} $$
Let
$$ \begin{align} a_1 b_1-a_2 b_2-a_3 b_3-a_4 b_4 & = w \newline a_1 b_2+a_2 b_1+a_3 b_4-a_4 b_3 & = x \newline a_1 b_3-a_2 b_4+a_3 b_1+a_4 b_2 & = y \newline a_1 b_4+a_2 b_3-a_3 b_2+a_4 b_1 & = z \end{align} \tag{2} $$
Integers $w, x, y, z$ exist for all $n \in \mathbb{Z}$ by Lagrange’s Four-square theorem and can be found efficiently using methods described in [PT2018] or [RS1986].
Let $n$ be an odd integer that we want to factor for which $n = w^2 + x^2 + y^2 + z^2$ is a 4-square representation. Let $r_4(n)$ denote the number of representations of a natural number $n$ as the sum of four squares. There are $r_4(n) = 8 \sigma(n)$ representations for odd $n$ where $\sigma(n)$ is the sum of divisors of $n$.
Is there any hope of solving Eqn. (2) and thereby factoring $n$? The closest method I found was [LB1968] for solving a linear diophantine equation in $n$ variables. But Eqn. (2) is not a linear diophantine equation and we also have 4 simultaneous equations.
We could set $-a_2 b_2-a_3 b_3-a_4 b_4 = 1$ and get $a_1 b_1 = w - 1$ and then find $a_1, b_1$ by factoring $w - 1$. This is a smaller factorization problem. But how do we proceed from there to find $a_2, b_2, a_3, b_3, a_4, b_4$?
Update (Aug 12, 2023):
Since $n$ is odd, it has a representation of the form $w^2 + x^2 + 2y^2$. Since $n = ab$, $a, b$ also have a representation in the same form. We will use this.
Take $a_1 = a_2, b_1 = b_2$ in the Euler identity.
$$ \begin{align} (a_1 b_1-a_1 b_1)-a_3 b_3-a_4 b_4 & = -a_3 b_3-a_4 b_4 = w \newline (a_1 b_1+a_1 b_1)+a_3 b_4-a_4 b_3 & = 2a_1 b_1 + a_3 b_4-a_4 b_3 = x \newline a_1 b_3-a_1 b_4+a_3 b_1+a_4 b_1 & = y \newline a_1 b_4+a_1 b_3-a_3 b_1+a_4 b_1 & = z \end{align} \tag{3} $$
Let $y=z$.
Therefore,
$$ a_1 b_3-a_1 b_4+a_3 b_1+a_4 b_1 = a_1 b_4+a_1 b_3-a_3 b_1+a_4 b_1\\ \implies -a_1 b_4+a_3 b_1 = a_1 b_4-a_3 b_1\\ \implies 2a_3 b_1 = 2a_1 b_4\\ \implies a_3 b_1 = a_1 b_4\\ $$
Substituting $a_3 b_1 = a_1 b_4$, we get
$$ \begin{align} -a_3 b_3-a_4 b_4 & = w \newline 2a_1 b_1 + a_3 b_4-a_4 b_3 & = x \newline a_1 b_3+a_4 b_1 & = y \newline \end{align} \tag{4} $$
By change of signs on $w$, we can take
$$ \begin{align} a_3 b_3+a_4 b_4 & = -w \newline 2a_1 b_1 + a_3 b_4-a_4 b_3 & = x \newline a_1 b_3+a_4 b_1 & = y \newline \end{align} \tag{4} $$
We have $\max(a_1^2, a_3^2, a_4^2) \lt \sqrt{n}$. Therefore, $\max(a_1, a_3, a_4) \lt \sqrt[4]{n}$.
What else could we do to solve this cleverly?
References:
[PT2018]: Pollack, P.; Treviño, E. (2018). "Finding the four squares in Lagrange's theorem" (PDF). Integers. 18A: A15.
[RS1986]: Rabin, M. O.; Shallit, J. O. (1986). "Randomized Algorithms in Number Theory". Communications on Pure and Applied Mathematics. 39 (S1): S239–S256. doi:10.1002/cpa.3160390713.
[LB1968]: Bernstein, L.; (1968) "The Linear Diophantine Equation in n Variables and its Application to Generalized Fibonacci Numbers".