官方题解只能做到 $O(\sqrt{n \log n})$？玩的太差了。这题可以把关于 $n$ 的幂次完全消掉，做到确定性的多对数复杂度。官方题解最后给出的复杂度是 $O(\sqrt{n\log n})$；下面不再沿着阈值分治继续优化，而是将问题化为固定维数的加权整数点计数。([QOJ][1])
记 $L=1+\lceil\log_2(n+1)\rceil$。下面实现的复杂度为 $O(L^{1+\log_{3/2}3})$ 次 $O(L)$ 位整数运算，即约 $O(\log^{3.71}n)$ 次算术操作。按朴素大整数运算计费，位复杂度为 $O(L^{3+\log_{3/2}3})$，约 $O(\log^{5.71}n)$；空间为 $O(L\log L)$ 比特。这里的指数来自三维锥分解，是算法的上界，不是已经证明的最优下界。
将问题化为四个加权求和
交换三次染色的顺序，使 $a_1\le a_2\le a_3$，并记 $a=a_1$。
先只考虑第一条等差数列。相邻两个黑点 $x,x+a$ 之间有 $a-1$ 个白点，贡献为 $\binom a2$。最关键的观察是：另外两条等差数列在这个开区间内都至多出现一个点。
若加入一个黑点 $x+u$，其中 $1\le u< a$，贡献减少 $u(a-u)$。
若另外两条等差数列分别出现于 $x+u,x+v$，且 $u\le v$，把两个减少量直接相加会多减 $u(a-v)$，需要将它加回来。因此，这个间隔的贡献是 $\binom a2-u(a-u)-v(a-v)+u(a-v)$。这个式子在 $u=v$ 时也成立，恰好处理了两次染色落在同一个位置的情况。
设第一条数列在 $[0,n)$ 内有 $K+1$ 个黑点，则它们之间有 $K$ 个完整间隔，其左端点为 $x_k=b_1+ak$，其中 $0\le k< K$。对于第 $j$ 条数列，定义相对位置 $d_j=a_jt_j-ak+b_j-b_1$。
记 $T_j$ 为所有满足 $0\le k< K$、$1\le d_j\le a-1$ 的整数对 $(k,t_j)$ 上，权值 $d_j(a-d_j)$ 的总和。
再定义两个修正量：$C_{23}$ 在所有满足 $0\le k< K$、$1\le d_2\le d_3\le a-1$ 的整数三元组上求和，权值为 $d_2(a-d_3)$；$C_{32}$ 在所有满足 $0\le k< K$、$1\le d_3< d_2\le a-1$ 的整数三元组上求和，权值为 $d_3(a-d_2)$。
于是答案为
$$ \mathrm{Ans}=\mathrm{Boundary}+K\binom a2-T_2-T_3+C_{23}+C_{32}. $$
这里 $\mathrm{Boundary}$ 是第一个黑点之前、最后一个黑点之后的贡献。加上哨兵 $-1,n$ 后，这两个间隔的长度都不超过 $a$，另外每条数列仍然至多出现一个点，直接计算即可。
至此，原问题只剩下四个求和：$T_2,T_3$ 是二维平行四边形内的二次加权整数点求和；$C_{23},C_{32}$ 是三维三棱柱内的二次加权整数点求和。它们的面数、顶点数、维数和权值次数都是常数。
整个推导没有使用互质条件，因此也适用于题目中公差不互质、黑点重合的情况。题目允许 $n\le10^{13}$、$a_i\le2n$，下面也不会枚举公差或最小公倍数。([QOJ][2])
多对数计算这些加权和
我们实现一个统一的计数器，计算 $\sum_{\mathbf z\in P\cap\mathbb Z^d}(\mathbf u\cdot\mathbf z+c_u)(\mathbf v\cdot\mathbf z+c_v)$，其中 $d$ 只可能是 $2$ 或 $3$。
首先处理退化情况。所有约束的左侧在整数点上都是整数，因此可以将整数下界减少 $1/2$、整数上界增加 $1/2$，而不改变整数点集合。例如，将 $0\le k
这样，即使 $K=1$，得到的多面体仍然是满维的。$a=1$ 时答案直接为零；其余情况下，代码中的两个三棱柱也都是满维的。所有顶点均用分子、分母表示，不需要浮点数。
接下来使用顶点锥分解。 Brion 定理说明，一个有界有理多面体的整数点生成函数，等于各个顶点处切锥的整数点生成函数之和，因此只需要处理常数个锥。([arXiv][3])
为了避免处理锥分解时的公共边界，我们在对偶空间进行分解：对偶空间中的低维锥，变回原空间后包含直线，其有理生成函数为零，可以直接丢弃。这是 Barvinok 算法中使用的对偶化技巧。([arXiv][3])
设一个对偶锥的生成向量是整数矩阵 $B$ 的列，记 $\Delta=|\det B|$。当 $\Delta=1$ 时，它是幺模锥，可以直接求生成函数；否则需要找到一个整数向量 $\mathbf w=B\boldsymbol\alpha$，使每个 $|\alpha_i|$ 都很小。
为此，在格 $\operatorname{adj}(B)\mathbb Z^d$ 上运行精确 LLL。这个格的行列式是 $\Delta^{d-1}$，LLL 可以找到非零向量 $\mathbf h$，满足 $|\mathbf h|_\infty\le 2^{(d-1)/4}\Delta^{(d-1)/d}$。代码同时维护其对应的整数向量 $\mathbf w$。这个短向量界来自 LLL 的标准行列式界。([LSA Technology Services][4])
然后从 $\boldsymbol\alpha$ 的每个坐标减去最近的整数，相当于从 $\mathbf w$ 中减去若干列向量的整数倍。这样同时保证 $|\alpha_i|\le1/2$。必要时将 $\mathbf w$ 取反，使至少一个 $\alpha_i>0$。
用 $\mathbf w$ 替换 $B$ 的第 $i$ 列得到子锥，其符号是 $\operatorname{sgn}(\alpha_i)$，行列式绝对值为 $|\alpha_i|\Delta$。因此，在三维情况下，每次至多产生三个子锥，并满足 $\Delta'\le\min(\Delta/2,\sqrt2,\Delta^{2/3})$。
所以递归深度为 $O(\log\log\Delta)$，递归树大小为 $O((\log\Delta)^{\log_{3/2}3})$。本题初始行列式至多为 $a_2a_3=O(n^2)$，因而总节点数为 $O(L^{\log_{3/2}3})$。
每个节点的精确 LLL 和最大公约数计算需要 $O(L)$ 次整数运算。又因为分解向量的坐标每层至多放大常数倍、深度只有 $O(\log L)$，所有中间整数始终只有 $O(L)$ 位。这就得到了开头给出的复杂度。代码深度优先处理各个锥，不保存整棵递归树。
还剩下最后一步：如何在幺模锥上计算二次权值，而不只是计算点数。
设顶点为 $\boldsymbol\xi$，对偶锥矩阵 $B$ 满足 $\det B=\pm1$，令 $U=B^{-T}$。原空间锥的整数点恰好可以写成 $\mathbf p+\sum_i m_i\mathbf r_i$，其中 $\mathbf r_i$ 是 $U$ 的列、$m_i\ge0$，并且 $\mathbf p=U\lceil B^T\boldsymbol\xi\rceil$。因此生成函数是
$$ G(\mathbf z)=\frac{\mathbf z^{\mathbf p}}{\prod_{i=1}^d(1-\mathbf z^{\mathbf r_i})}. $$
引入两个满足 $\varepsilon^2=\eta^2=0$ 的形式变量。对整数点 $\mathbf x$，令 $E(\mathbf x)=\boldsymbol\lambda\cdot\mathbf x+\varepsilon(\mathbf u\cdot\mathbf x+c_u)+\eta(\mathbf v\cdot\mathbf x+c_v)$。那么 $\exp(tE(\mathbf x))$ 中 $t^2\varepsilon\eta$ 的系数，恰好就是需要的权值。
对锥的每条射线，记 $q_i=\boldsymbol\lambda\cdot\mathbf r_i+\varepsilon(\mathbf u\cdot\mathbf r_i)+\eta(\mathbf v\cdot\mathbf r_i)$。利用 $B_0(z)=z/(e^z-1)=1-z/2+z^2/12-z^4/720+O(z^6)$，只需要求
$$ [t^{d+2}\varepsilon\eta]\; \frac{(-1)^d}{\prod_iq_i}\, e^{tE(\mathbf p)}\prod_i B_0(tq_i). $$
由于 $d\le3$，多项式最高只保留到五次，每个叶子只需常数次代数运算。
这里还有一个必须处理的问题：不能随便选择 $\boldsymbol\lambda$，否则某条射线可能导致分母为零。下面的实现不使用随机试值，而是在三次扩域 $\mathbb F_p[\theta]/(\theta^3-\theta-1)$ 中选取 $\boldsymbol\lambda=(1,\theta,\theta^2)$，其中 $p=998244353$。
多项式 $X^3-X-1$ 在这个模数下不可约，验证信息写在代码注释中。幺模矩阵的每一列都是本原向量，不可能所有坐标都被 $p$ 整除，因此 $\boldsymbol\lambda\cdot\mathbf r_i$ 一定非零。这样，所有求逆都有保证。最终结果属于原来的 $\mathbb F_p$，$\theta,\theta^2$ 的系数会严格消失。
代码包含精确 LLL、带符号对偶锥分解、三次扩域以及二次权值提取。对黑点出现次数不超过固定常数 $512$ 的情况直接枚举，不影响多对数渐近复杂度。
提交记录：https://qoj.ac/submission/2939184 。数据范围开太小了，谴责一下出题人。