定义三维空间中的一个格 (Lattice) 是满足下列条件的三维整点的集合：$$ L = \left\{ u \cdot \mathbf r_1 + v \cdot \mathbf r_2 + w \cdot \mathbf r_3 \mid u, v, w \in \mathbb Z \right\} $$ 其中 $\mathbf r_1, \mathbf r_2, \mathbf r_3$ 为三维空间中的整系数坐标向量，满足：
现在有 $n$ 个三维整点 (向量) 构成的集合 $A = \left\{ \mathbf a_1, \mathbf a_2, \cdots, \mathbf a_n \right\}$，其中 $\mathbf a_i = \left( x_i, y_i, z_i \right)$。
你需要找到一个格 $L$ 满足 $A \subset L$ 且组成 $L$ 的三个基向量的长度 (即 $r$) 尽可能大。
第一行包含一个正整数 $n$ ($n \leq 10000$)，表示三维向量的个数。
接下来 $n$ 行，每行三个整数 $x_i, y_i, z_i$ ($0 < x_i^2 + y_i^2 + z_i^2 \leq 10^{16}$)，表示第 $i$ 个向量的坐标。保证 $\gcd \left( x_1, y_1, z_1, x_2, y_2, z_2, \cdots, x_n, y_n, z_n \right) = 1$。
第一行输出一个整数，表示基向量的长度的最大可能值的平方 $r^2$。
接下来三行，每行三个整数 $x, y, z$，描述其中一个基向量的坐标。
如果有多组可能的解，输出任意一组均可。
首先，由于 $\mathbf r_1, \mathbf r_2, \mathbf r_3$ 两两正交，因此 $\mathbf r_3$ 和 $\mathbf r_1 \times \mathbf r_2$ 是共线向量。
而 $\left| \mathbf r_1 \times \mathbf r_2 \right| = r^2, \left| \mathbf r_3 \right| = r$，说明 $\mathbf r_1 \times \mathbf r_2 = \pm r \cdot \mathbf r_3 \Rightarrow r \in \mathbb N^*$。
其次，如果 $\mathbf a = u \cdot \mathbf r_1 + v \cdot \mathbf r_2 + w \cdot \mathbf r_3$，那么有 $\left| \mathbf a \right|^2 = \left( u^2 + v^2 + w^2 \right) r^2 \Rightarrow r^2 \mid \left| \mathbf a \right|^2$。
也就是说，最终的 $r$ 需要满足一个必要条件：$r^2 \mid \gcd \left( x_1^2 + y_1^2 + z_1^2, x_2^2 + y_2^2 + z_2^2, \cdots, x_n^2 + y_n^2 + z_n^2 \right)$。
事实上，这样的 $r$ 的个数是不超过 $\max\limits_{n \leq 10^9} d \left( n \right) = 1344$ 个的。这为我们对 $r$ 逐一检验带来了可能性。
下面考虑对于一组满足必要条件的 $r$，检验它是否可能成为真正的答案 (即是否存在这样一个满足条件的格基)，如果这个检验过程可以在一个可接受的时间复杂度内完成，那么我们只需要按照所有候选 $r$ 从大到小依次检验即可。
众所周知，二维平面上的变换我们一般用复数来刻画，其中相当于整数在实数中的地位的，被称为 Gauss 整数 (即 $a + b \mathrm i$ 满足 $a, b \in \mathbb Z$)。所有的 Gauss 整数关于 $+$ 和 $\times$ 运算构成一个 Euclid 整环。
那么，在三维空间中的变换，我们一般使用四元数 (Quaternion) $a + b \mathrm i + c \mathrm j + d \mathrm k$ 来刻画，而在四元数中，相当于整数的地位的，被称为 Hurwitz 四元数。
(ps: 四元数的运算满足结合律，不满足交换律，虚数单位之间的运算满足 $\mathrm i^2 = \mathrm j^2 = \mathrm k^2 = \mathrm i \mathrm j \mathrm k = -1$)
一个最直观的定义是，$a + b \mathrm i + c \mathrm j + d \mathrm k$，其中 $a, b, c, d \in \mathbb Z$)。可惜的是，这样定义出来的四元数集合不构成 Euclid 整环，而这道题中的操作就是类似于 $\gcd$ 的一个过程，不能 Euclid 看起来就很难受 QAQ (ps: 其实这样定义出来的四元数被称为 Lipschitz 四元数)。
因此我们需要扩充一下它的定义，考虑增加部分的有理四元数。在复数中，$\left| a + b \mathrm i \right| = \sqrt{a^2 + b^2}$，考虑 $z = a + b \mathrm i$，且满足 $a^2 + b^2 \in \mathbb Z$，这样可能会有一些 $\dfrac 45 + \dfrac 35 \mathrm i$ 这类玩意，虽然它的模长为整数，但是它可以和 $1, \mathrm i$ 整系数表出 $\dfrac {1 + \mathrm i} 5$，模长就不为整数。
于是在 $\mathbb C$ 中，如果能扩充的话，我们肯定要考虑 $a + b \mathrm i$，其中 $\dfrac 1a, \dfrac 1b \in \mathbb Z$。那么，在这种情况下，除了 $a = b = 1$，其余所有的有理复数的模长均不等于 $1$，因此就没办法扩充了。
但是，这件事情一旦搬到四元数 $\mathbb H$ 来就不一样了：因为容易验证当 $a = b = c = d = 2$ 时，$\left| \dfrac {1 + \mathrm i + \mathrm j + \mathrm k} 2 \right| = 1$！(事实上它也是唯一的可扩充候选)
因此，如果我们将 $\omega = \dfrac {1 + \mathrm i + \mathrm j + \mathrm k} 2$ 加入刚才定义进的 Lipschitz 四元数集合，并取它们的正系数表示，就得到了 Hurwitz 四元数：$$ H = \left\{ a + b \mathrm i + c \mathrm j + d \mathrm k \mid a, b, c, d \in \mathbb Z \vee a, b, c, d \in \mathbb Z + \frac 12 \right\} $$ (事实上有 $\omega^6 = 1$，即 $\omega$ 还是四元数的一个单位呢！因此怎么能抛弃 $\omega$ 小可爱呢~)
可以证明，虽然没有交换律 (因此 $H$ 肯定不是整环了 qwq)，但是可以证明，$H$ 中依然有类似地唯一分解定理 (但是表述可能会有点区别)，甚至依然可以做 Euclid 算法！(ps: 注意 $H$ [Hurwitz 四元数] 和 $\mathbb H$ [一般四元数] 的区别)
具体地，在 $H$ 中，对于两个数 $A, B$，可以定义它的左商、左余数，右商、右余数，以及左 $\gcd$ 和右 $\gcd$。如，左商和左余数对应表达式 $A = q B + r$，而左 $\gcd$ 就表示模最大的 $d$ 使得 $A = d u, B = d v$。
有了 Hurwitz 四元数后，整个问题就瞬间清晰了很多：
首先，对于三维空间中的向量 $\mathbf a = \left( x, y, z \right)$，我们将其对应到四元数 $A = x \mathrm i + y \mathrm j + z \mathrm k$。那么三维空间中的其中一组标准正交基就可以表示为 $\mathrm i, \mathrm j, \mathrm k$。
由立体几何的基础知识可得，三维空间中的旋转变换可以这样表示：$$ f \left( w \right) = \frac {q \cdot w \cdot \bar q} {q \cdot \bar q} = \frac {q \cdot w \cdot \bar q} {\left| q \right|^2} $$ 其中 $w$ 表示实部为 $0$ 的四元数 (通常表示向量)，$q = a + b \mathrm i + c \mathrm j + d \mathrm k$ 为一个特定的四元数，而 $\bar q = a - b \mathrm i - c \mathrm j - d \mathrm k$ 表示 $q$ 的共轭。
于是，$\mathbf r_1, \mathbf r_2, \mathbf r_3$ 可以看成标准正交基经过某个旋转并缩放后的结果，因此有 \begin{cases} \mathbf r_1 = r \cdot \dfrac {q \cdot \mathrm i \cdot \bar q} {\left| q \right|^2} \\ \mathbf r_2 = r \cdot \dfrac {q \cdot \mathrm j \cdot \bar q} {\left| q \right|^2} \\ \mathbf r_3 = r \cdot \dfrac {q \cdot \mathrm k \cdot \bar q} {\left| q \right|^2} \end{cases}
记比例系数 $c = \dfrac r {\left| q \right|^2} = \dfrac r {a^2 + b^2 + c^2 + d^2}$，则可以证明 $c$ 的分母是 $4$ 的约数 (因为可以不妨假设 $\gcd \left( a, b, c, d \right) = 1$)。
当然，如果我们放宽要求令 $c \in H$，那么进一步可以证明 $c$ 的分母是 $2$ 的约数 (即 $2 c \in \mathbb Z$)。通过适当调整使得 $q$ 为 Lipschitz 四元数，那么由题目条件可知 $c = 1$。
定义 $f \left( w \right) = q \cdot w \cdot \bar q$，那么存在一组向量 $\mathbf b_1, \mathbf b_2, \cdots, \mathbf b_n$ 满足 $$ \mathbf a_i = f \left( \mathbf b_i \right) = q \cdot \mathbf b_i \cdot \bar q $$ (ps: 这里我们直接将一个向量表示为对应的四元数)
于是，此时有 $q \mid \mathbf a_i$ (这里指的是作为左因子)，因此 $q \mid \gcd \left( \mathbf a_1, \mathbf a_2 \cdots, \mathbf a_n \right)$。
而且，我们还枚举了 $r$，因此必须有 $\left| q \right|^2 = r$，由唯一分解定理可知 $q \mid r$ (作为左因子)。
于是，$q \mid \gcd \left( \mathbf a_1, \mathbf a_2 \cdots, \mathbf a_n, r \right)$。事实上，设 $q_0 = \gcd \left( \mathbf a_1, \mathbf a_2 \cdots, \mathbf a_n, r \right)$，则 $\left| q_0 \right|^2 \mid r$，因此只有当 $\left| q_0 \right|^2 = r$ 时才有必要继续验证下去。
此时我们只需要反解出 $\mathbf b_i = \dfrac {\bar q \cdot \mathbf a_i \cdot q} {r^2}$ 并检验它们是否是整数即可。故这部分的总时间复杂度为 $O \left( n + \log W \right)$ 的 ($\log W$ 指的是 $\gcd$ 复杂度)。
于是时间复杂度的一个上界为 $O \left( W^{1/3} + \max d \cdot \left( n + \log W \right) \right)$，可以通过。
#include <bits/stdc++.h>
#define gcd std::__gcd
using std::cin;
using std::cout;
using std::vector;
typedef long long ll;
typedef long double ld;
const int N = 10054;
struct quaternion {
ld s, x, y, z;
quaternion (ld s0 = 0, ld x0 = 0, ld y0 = 0, ld z0 = 0) : s(s0), x(x0), y(y0), z(z0) {}
inline quaternion operator - () const {return quaternion(-s, -x, -y, -z);}
inline quaternion conj() const {return quaternion(s, -x, -y, -z);}
inline quaternion operator + (const quaternion &B) const {return quaternion(s + B.s, x + B.x, y + B.y, z + B.z);}
inline quaternion operator - (const quaternion &B) const {return quaternion(s - B.s, x - B.x, y - B.y, z - B.z);}
inline quaternion operator * (const ld k) const {return quaternion(s * k, x * k, y * k, z * k);}
inline quaternion operator * (const quaternion &B) const {
return quaternion(
s * B.s - x * B.x - y * B.y - z * B.z,
s * B.x + x * B.s + y * B.z - z * B.y,
s * B.y - x * B.z + y * B.s + z * B.x,
s * B.z + x * B.y - y * B.x + z * B.s
);
}
inline operator void * () const {return s || x || y || z ? (void *)this : 0;}
inline ld norm2() const {return s * s + x * x + y * y + z * z;}
inline ll intnorm2() const {return llroundl(norm2());}
inline quaternion inv() const {return conj() * (1.l / norm2());}
inline quaternion round() const {
quaternion A(roundl(s), roundl(x), roundl(y), roundl(z)), B(floorl(s) + .5l, floorl(x) + .5l, floorl(y) + .5l, floorl(z) + .5l );
return (*this - A).norm2() <= (*this - B).norm2() ? A : B;
}
inline friend quaternion Ldiv(const quaternion &A, const quaternion &B) {return (B.inv() * A).round();}
inline friend quaternion Rdiv(const quaternion &A, const quaternion &B) {return (A * B.inv()).round();}
inline friend quaternion Lmod(const quaternion &A, const quaternion &B) {return A - B * Ldiv(A, B);}
inline friend quaternion Rmod(const quaternion &A, const quaternion &B) {return A - Rdiv(A, B) * B;}
// find (Hurwitz) quaternion d where A = d u, B = d v.
inline friend quaternion Lgcd(quaternion A, quaternion B) {
static quaternion tmp;
for (; B; tmp = Lmod(A, B), A = B, B = tmp);
return A;
}
inline bool is_int() const {return fabsl(s - roundl(s)) < 1e-10l && fabsl(x - roundl(x)) < 1e-10l && fabsl(y - roundl(y)) < 1e-10l && fabsl(z - roundl(z)) < 1e-10l;}
} p[N], q[N], v;
int n, cnt = 0;
int potential[2003731];
bool final(int r) {
// unique factorization for Hurwitz quaternion.
quaternion u = Lgcd(v, quaternion(r)), x, y, z;
ll r1 = u.norm2(); ld cf = 1. / ((ll)r * r);
if (assert(r % r1 == 0), r1 != r) return false;
for (int i = 0; i < n; ++i) {
q[i] = u.conj() * p[i] * u * cf;
if (!q[i].is_int()) return false;
}
x = u * quaternion(0, 1, 0, 0) * u.conj(),
y = u * quaternion(0, 0, 1, 0) * u.conj(),
z = u * quaternion(0, 0, 0, 1) * u.conj();
return cout << (ll)r * r << '\n'
<< llroundl(x.x) << ' ' << llroundl(x.y) << ' ' << llroundl(x.z) << '\n'
<< llroundl(y.x) << ' ' << llroundl(y.y) << ' ' << llroundl(y.z) << '\n'
<< llroundl(z.x) << ' ' << llroundl(z.y) << ' ' << llroundl(z.z) << '\n', true;
}
int main() {
int i, x, y, z, g = 0; ll d, D = 0;
std::ios::sync_with_stdio(false), cin.tie(NULL);
cin >> n;
for (i = 0; i < n; ++i)
cin >> x >> y >> z, g = gcd(gcd(gcd(g, x), y), z),
v = Lgcd(v, p[i] = quaternion(0, x, y, z)),
D = gcd(D, (ll)x * x + (ll)y * y + (ll)z * z);
assert(abs(g) == 1), d = v.intnorm2(), z = cbrtl(D);
for (i = 1; i <= z; ++i) if (!(D % i)) {
if (!(D % ((ll)i * i) || d % i)) potential[cnt++] = i;
if (y = sqrtl(D / i), (ll)y * y * i == D && !(d % y)) potential[cnt++] = y;
}
std::sort(potential, potential + cnt, std::greater <int> ());
for (i = 0; i < cnt; ++i) if (final(potential[i])) return 0;
return -1;
}
坑1：寻找所有候选 $r$ 的过程类似分解质因数，可以在 $O \left( W^{1/3} \right)$ 时间内完成，注意不要写错或超时。
坑2：注意 std::__gcd 可能会返回负数，判断时需要注意。