题解：P2019 四平方和定理

温馨提示：本题的实现基于 Miller Rabin 和 Pollard-Rho，不了解的读者可以自行学习，本篇题解不作赘述，此处是我的实现。

---

此题可以直接套用“雅可比四平方和定理”解决。

雅可比四平方和定理：自然数 \(n\) 表示为四个整数平方和的有序表示方式的数量为 \(8\sum_{d \mid n,4\nmid n}^{} d\)，证明过程在这里。

怎么得出这个结论（不严谨的打表大法）？

我们可以光速写出一个打表程序，可以发现所有的方案数都是 \(8\) 的倍数。我们可以把所有的数都除以 \(8\) 来观察规律。

来看看此时的数据：

number	有序整数对方案数	有序整数对方案数除以 \(8\)
\(10\)	\(144\)	\(18\)
\(5\)	\(48\)	\(6\)
\(13\)	\(112\)	\(14\)
\(6\)	\(96\)	\(12\)
\(7\)	\(64\)	\(8\)

因为这是一道数论题目，我们把他们的因数和列出来，容易发现有序整数对方案数就是 \(8\sum_{d \mid n}^{} d\)。可是我们看看当 \(n\) 为 \(4\) 的倍数时……

number	有序整数对方案数除以 \(8\)	\(\sum_{d \mid n}^{} d\)
\(8\)	\(3\)	\(13\)
\(4\)	\(3\)	\(7\)
\(24\)	\(12\)	\(60\)

出现了 bug。

但实际上我们可以把他们的因数都列出来：

number	因数
\(4\)	\(1,2,4\)
\(8\)	\(1,2,4,8\)
\(24\)	\(1,2,3,4,6,8,12,24\)

似乎我们凭借着我们优秀的注意力，发现 \(4\) 与 \(8\) 之的因数间多了一个 \(8\)，可有序整数对方案数除以 \(8\) 的值却还是一样的。我们自己再通过打表找规律（如 \(2\) 和 \(4\)），就不难发现方案数除以 \(8\) 其实是 \(n\) 的所有非 \(4\) 的倍数的因数的和，也就是说，自然数 \(n\) 表示为四个整数平方和的有序表示方式的数量为 \(8\sum_{d \mid n,4\nmid n}^{} d\)。

---

我们可以先根据雅可比四平方和定理敲出暴力代码，会有四个点超时。

我们可以从另外一个角度入手：

首先还补充一个小奥的知识：

对于任意大于一的整数 \(n\) 可以唯一的表示为 \(p_1^{u_1}p_2^{u_2}p_3^{u_3}\cdots p_m^{u_m}\)（其中 \(p_i\) 都是质数，\(u_i\) 不为 \(0\)），此时 \(\sum_{d \mid n}^{}d=\prod_{i=1}^{m}\sum_{j=0}^{u_i}p_i^j\)。

我们可以先将 \(n\) 疯狂除以 \(2\) 直至其不为 \(4\) 的倍数，因为这样子我们找到它的所有因数时，这些因数都不是 \(4\) 的倍数了。

然后方案数式子 \(8\sum_{d \mid n,4\nmid n}^{} d\) 就简化为了 \(8\sum_{d \mid n}^{} d\)，套用小奥的知识得到方案数为 \(8\prod_{i=1}^{m}\sum_{j=0}^{u_i}p_i^j\)（其中 \(p_i\) 都是质数）。

那么我们找质因数的过程可以怎么实现呢？没错！可以直接套用 Pollard-Rho。此时到这里你应该已经知道怎么写了。

AC 代码如下：

#include <bits/stdc++.h>
#define int long long
#define ull unsigned long long
using namespace std;
const int MOD = 1e9 + 7;

	// 省略 qmul（快速乘）、qpow（快速幂）、MillerRabin、gcd、PollardRho、fac（找出 x 的最大质因数） 的实现

int solve(int n) {
	// 让 n 不为 4 的倍数
	while (n % 4 == 0) {
		n /= 2;
	} 
	
	int ans = 1;
	
	while (n != 1) {
		maxfac = 0;
		fac(n);
		int there = 1;
		int temp = maxfac;
		while(n % maxfac == 0) {
			there += temp;
			temp *= maxfac;
			there %= MOD;
			n /= maxfac;
		}
		ans *= there;
		ans %= MOD;
	}
	
	return 8 * ans % MOD;
}
signed main() {
	int t; 
	cin >> t;
	while (t--) {
		int n;
		cin >> n;
		cout << solve(n) << endl; 
	}
	return 0;
}