Cubic Lattice
Contest: Codeforces Global Round 9
Difficulty rating: 3500
Tags: geometry, math, matrices, number theory

A cubic lattice $$$L$$$ in $$$3$$$-dimensional euclidean space is a set of points defined in the following way: $$$$$$L=\{u \cdot \vec r_1 + v \cdot \vec r_2 + w \cdot \vec r_3\}_{u, v, w \in \mathbb Z}$$$$$$ Where $$$\vec r_1, \vec r_2, \vec r_3 \in \mathbb{Z}^3$$$ are some integer vectors such that:

- $$$\vec r_1$$$, $$$\vec r_2$$$ and $$$\vec r_3$$$ are pairwise orthogonal: $$$$$$\vec r_1 \cdot \vec r_2 = \vec r_1 \cdot \vec r_3 = \vec r_2 \cdot \vec r_3 = 0$$$$$$ Where $$$\vec a \cdot \vec b$$$ is a dot product of vectors $$$\vec a$$$ and $$$\vec b$$$.
- $$$\vec r_1$$$, $$$\vec r_2$$$ and $$$\vec r_3$$$ all have the same length: $$$$$$|\vec r_1| = |\vec r_2| = |\vec r_3| = r$$$$$$

You have to find a cubic lattice $$$L$$$ such that $$$A \subset L$$$ and $$$r$$$ is the maximum possible.

Input format:
First line contains single integer $$$n$$$ ($$$1 \leq n \leq 10^4$$$) — the number of points in $$$A$$$.

The $$$i$$$-th of the following $$$n$$$ lines contains integers $$$x_i$$$, $$$y_i$$$, $$$z_i$$$ ($$$0 < x_i^2 + y_i^2 + z_i^2 \leq 10^{16}$$$) — coordinates of the $$$i$$$-th point.

It is guaranteed that $$$\gcd(g_1,g_2,\dots,g_n)=1$$$ where $$$g_i=\gcd(x_i,y_i,z_i)$$$.

Output format:
In first line output a single integer $$$r^2$$$, the square of maximum possible $$$r$$$.

In following $$$3$$$ lines output coordinates of vectors $$$\vec r_1$$$, $$$\vec r_2$$$ and $$$\vec r_3$$$ respectively.

If there are multiple possible answers, output any.

Constraints:
Time limit: 4.0 s
Memory limit: 256.0 MB

Example input:
2
1 2 3
1 2 1

Example output:
1
1 0 0
0 1 0
0 0 1

Example input:
1
1 2 2

Example output:
9
2 -2 1
1 2 2
-2 -1 2

Editorial:
$$$2$$$-dimensional version of this problem appeared as problem K in Makoto Soejima Contest 4 from Petrozavodsk Winter Training Camp 2016 (Also known as GP of Asia in OpenCup XVI). Solution for this case is based on the fact that the whole lattice equation may be represented as the product of two Gaussian integers $$$a_k = \alpha_k \cdot g$$$, where $$$g$$$ is the complex number representing $$$(\vec r_1, \vec r_2)$$$ basis. Thus, the solution was to find largest $$$g$$$ which divides every point $$$a_k$$$. This can be done with Euclidean algorithm.
Similar idea applies to this problem except for rotations in $$$3$$$-dimensional case are represented with quaternions rather than complex numbers. You may get some basic notion of quaternions from this article. Note that in $$$3$$$-dimensional space if integer vectors $$$\vec r_1$$$, $$$\vec r_2$$$ and $$$\vec r_3$$$ are pairwise orthogonal and have the same length, then their length is an integer number as well. This is due to the fact that $$$\vec r_3$$$ may be represented explicitly in the following form:
$$$$$$\vec r_3 = \pm\frac{\vec r_1 \times \vec r_2}{|r_1|}$$$$$$
We should note here that both $$$\vec r_3$$$ and $$$\vec r_1 \times \vec r_2$$$ have integer coordinates, thus $$$|r_1|$$$ is rational. But since $$$|r_1|^2$$$ is integer, we may conclude that $$$|r_1|$$$ is integer as well. This brings us to the conclusion that vectors $$$\frac{\vec r_1}{|\vec r_1|}$$$, $$$\frac{\vec r_2}{|\vec r_2|}$$$ and $$$\frac{\vec r_3}{|\vec r_3|}$$$ have rational coordinates, thus the transformation which maps basis vectors $$$\vec e_x$$$, $$$\vec e_y$$$ and $$$\vec e_z$$$ into $$$\vec r_1$$$, $$$\vec r_2$$$ and $$$\vec r_3$$$ may be represented as the combination of rational rotation and scaling, that is:
$$$$$$\vec v \mapsto k \cdot \frac{q \cdot \vec v \cdot \bar q}{q \cdot \bar q}$$$$$$
Here $$$q = s + ix + jy + kz$$$ is an integer quaternion representing the rotation, $$$\bar q = s - ix - jy - kz$$$ is its conjugate quaternion, and $$$k=r$$$ is some integer scaling factor.
If we write this transform explicitly in matrix form, its matrix will be as follows:
$$$$$$ \frac{k}{s^2+x^2+y^2+z^2} \cdot \begin{pmatrix} s^2+x^2-y^2-z^2 & 2xy-2sz & 2xz+2sy \\ 2xy+2sz & s^2-x^2+y^2-z^2 & 2yz-2sx \\ 2xz-2sy & 2yz+2sx & s^2-x^2-y^2+z^2 \end{pmatrix} $$$$$$
Since this transform maps $$$\vec e_x$$$, $$$\vec e_y$$$ and $$$\vec e_z$$$ to integer vectors $$$\vec r_1$$$, $$$\vec r_2$$$ and $$$\vec r_3$$$, all matrix elements should be integer. That is, if we reduce the fraction $$$\frac{k}{s^2+x^2+y^2+z^2}$$$, its denominator should divide every matrix element. Without loss of generality we may assume that $$$\gcd(s,x,y,z)=1$$$, because otherwise we may simply reduce all quaternion coordinates by this common divisor and matrix elements will stay same. Now if $$$\frac{k}{s^2+x^2+y^2+z^2}=\frac{P}{Q}$$$ then $$$Q$$$ should divide every element of the matrix and it also should divide $$$s^2+x^2+y^2+z^2$$$. Thus, it divides:
$$$$$$\gcd\begin{pmatrix} s^2+x^2+y^2+z^2\\ s^2+x^2-y^2-z^2\\ s^2-x^2+y^2-z^2\\ s^2-x^2-y^2+z^2 \end{pmatrix}=\gcd\begin{pmatrix} s^2+x^2+y^2+z^2\\ 2(y^2+z^2)\\ 2(x^2+z^2)\\ 2(x^2+y^2) \end{pmatrix}=\gcd\begin{pmatrix} s^2+x^2+y^2+z^2\\ 2(y^2+z^2)\\ 2(x^2+z^2)\\ 4z^2 \end{pmatrix}$$$$$$
Here we utilized the fact that $$$\gcd(a\pm b, b)=\gcd(a, b)=\gcd(-a,b)$$$. Next thing to note is that $$$\gcd(a, b)$$$ always divides $$$\gcd(k\cdot a, b)$$$, thus $$$Q$$$ should divide:
$$$$$$\gcd\begin{pmatrix} 4(s^2+x^2+y^2+z^2)\\ 4(y^2+z^2)\\ 4(x^2+z^2)\\ 4z^2 \end{pmatrix}=4\gcd(s^2,x^2,y^2,z^2)=4 $$$$$$
Which means that every transform we're looking for may be represented in the following form:
$$$$$$ \vec v \mapsto \frac{k}{4} \cdot q \cdot \vec v \cdot \bar q $$$$$$
Where $$$q$$$ is some integer quaternion with $$$\gcd(s,x,y,z)=1$$$ and $$$k$$$ is some integer scaling factor. Quaternions having all integer components are called Lipschitz quaternions. But in this problem we will need to find greatest common divisor of some quaternions and using Lipschitz quaternions won't be enough for us because they don't form an Euclidean domain. Instead, we'll stick to Hurwitz quaternions in which it's also allowed for quaternion to have all of its coordinates semi-integer. With such quaternions we may further divide quaternion by $$$2$$$ if all its components are odd, which will reduce the set of possible denominators to $$$\{1,2\}$$$ only, thus in Hurwitz quaternions any transform may be represented as:
$$$$$$ \vec v \mapsto \frac{k}{2} \cdot q \cdot \vec v \cdot \bar q $$$$$$
Now, given this knowledge, we have to find a Hurwitz quaternion $$$q$$$ and scaling factor $$$k$$$ such that every point from the set $$$A$$$ is the image of some integer vector under this transform. Basically, we assume here that $$$\vec r_1$$$, $$$\vec r_2$$$ and $$$\vec r_3$$$ are images of unit basis vectors $$$\vec e_x$$$, $$$\vec e_y$$$ and $$$\vec e_z$$$. If this is true, then $$$a_k = u_k \cdot \vec r_1 + v_k \cdot \vec r_2 + w_k \cdot \vec r_3$$$ is the image of the vector $$$b_k$$$ defined as $$$b_k = u_k \cdot \vec e_x + v_k \cdot \vec e_y + w_k \cdot \vec e_z$$$.
For simplicity we will further work with vectors $$$a_i$$$ multiplied by $$$2$$$, as in this way we may safely assume that they were obtained with $$$v \mapsto k \cdot q \cdot \vec v \cdot \bar q$$$ mapping. Now we should further work with the quaternion $$$g=\gcd(a_1,\dots,a_n)$$$ which is obtained as GCD of input vectors as if they were Hurwitz quaternions and with the number $$$G=\gcd(\|a_1\|,\dots,\|a_n\|)$$$, where $$$\|q\|=q\cdot\bar q = s^2+x^2+y^2+z^2$$$ is the norm of $$$q$$$. Among major properties of quaternionic norm we should note that it's multiplicative, that is $$$\|a\cdot b\| =\|a\| \cdot \|b\|$$$. Due to this we may conclude that $$$\|k\|\cdot \|q\|^2=k^2 \|q\|^2$$$ should always divide $$$G$$$. Thus with the number $$$G$$$ fixed there may only be at most $$$O(G^{1/6})$$$ possible values of $$$\|q\|$$$ which may be found in $$$O(G^{1/3})$$$ because if $$$a^2$$$ divides $$$b$$$ then either $$$a$$$ or $$$b/a$$$ is at most $$$b^{1/3}$$$.
In our problem $$$G$$$ may only be up to $$$10^{16}$$$, which makes it up to nearly $$$500$$$ candidate numbers for being equal to $$$\|q\|$$$ which may be found in $$$\approx 2.5\cdot 10^5$$$ arithmetic operations. Now to check if it's possible to obtain $$$q$$$ with $$$\|q\|$$$ being fixed we should look on the quaternion $$$g$$$. Under given constraints $$$g$$$ is primitive (not divisible by any integer constant larger than $$$1$$$). It may be proven thus that if $$$g$$$ is primitive and $$$\|g\|=ab$$$ then $$$g$$$ may be uniquely (up to units) represented as $$$g=qp$$$ where $$$\|q\|=a$$$ and $$$\|p\|=b$$$. Due to this if we fix $$$\|q\|$$$ we may find actual $$$q$$$ as $$$\gcd(g,\|q\|)$$$ because $$$q\bar q =\|q\|$$$. After this we only have to check if this $$$q$$$ actually produces all points from the input set, which may be done in $$$O(n)$$$. Therefore, the total complexity of described solution is $$$O(G^{1/3}+G^{1/6}n)$$$.