P2019 四平方和定理题解

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

雅可比四平方和定理

聪明的雅可比证明了：自然数 n 表示为四个整数平方和的有序表示方式的数量为：
$$8\sum _{d\mid n,4\nmid n}d$$

证明过程在这里：https://en.wikipedia.org/wiki/Jacobi%27s_four-square_theorem。

芝士铺垫

Miller Rabin

好解析：OIwiki
& dalao 好文章

我们小学二年级都学过费马小定理：

如果 $p$ 是一个质数，而 $a$ 不是 $p$ 的质数，那么 $a^{p-1}\equiv 1 \pmod{p} $

原理有点复杂，可以看看上面的文章，因为篇幅原因不多展开。

代码如下：

int qmul(int a, int b, int p) { // 快速乘
	ull c = (ull)a * (ull)b - (ull)((long double)a / p * b + 0.5L) * (ull)p;
	if (c < (ull)p) return c;
	return c + p;
}
int qpow(int a, int b, int p) { // 快速幂
	int res = 1;
	while (b > 0) {
		if (b & 1) {
			res = qmul(res, a, p);
		}
		a = qmul(a, a, p);
		b /= 2;
	}
	return res;
}
bool MillerRabin(int n) {
	if (n < 2) {
		return false;
	}
	if (n == 2 || n == 3) {
		return true;
	}
	int d = n - 1, r = 0;
	while (!(d & 1)) {
		r++;
		d /= 2;
	}
	for (int k = 0; k < 10; ++k) {
		int a = rand() % (n - 2) + 2;
		int x = qpow(a, d, n);
		if (x == 1 || x == n - 1) {
			continue;
		}
		for (int i = 0; i < r - 1; i++) {
			x = qmul(x, x, n);
			if (x == n - 1) {
				break;
			}
		}
		if (x != n - 1) {
			return false;
		}
	}
	return true;
}

（虽然这个代码放到高精度里会炸，但既然都判断高精度是不是质数了……出题人非人哉）

Pollard-Rho

这个太过于奇妙了，导致我一点都解释不了，可以自行 BD。

代码如下：

int qmul(int a, int b, int p) {
	// 省略实现
}
int qpow(int a, int b, int p) { 
	// 省略实现
}
bool MillerRabin(int n) {
	// 省略实现
}

int gcd(int a, int b) {
	if (b == 0) {
		return a;
	}
	return gcd(b, a % b);
}
int PollardRho(int x) {
	int s = 0, t = 0, c = rand() % (x - 1) + 1;
	for (int g = 1; ; g *= 2, s = t) {
		int val = 1;
		for (int k = 1; k <= g; k++) {
			t = (qmul(t, t, x) + c) % x;
			val = qmul(val, abs(t - s), x);
			if ((k % 127) == 0) {
				int d = gcd(val, x);
				if (d > 1) {
					return d;
				}
			}
		}
		int d = gcd(val, x);
		if (d > 1) {
			return d;
		}
	}
}
int maxfac;
void fac(int x) {  
	if (x <= maxfac || x < 2) return;
	if (MillerRabin(x)) {             
		maxfac = max(maxfac, x); 
		return;
	}
	int p = x;
	while (p >= x) {
		p = PollardRho(x); 
	}
	while ((x % p) == 0) {
		x /= p;
	}
	fac(x);
	fac(p);
}

所以有了上面的芝士，我们可以先 AC 这道题

AC code：

#include <bits/stdc++.h>
#define int long long
#define ull unsigned long long
using namespace std;
const int MOD = 1e9 + 7;
int qmul(int a, int b, int p) { 
	// 省略实现
}
int qpow(int a, int b, int p) {
	// 省略实现
}
bool MillerRabin(int n) {
	// 省略实现
}

int gcd(int a, int b) {
	// 省略实现
}
int PollardRho(int x) {
	// 省略实现
}
int maxfac;
void fac(int x) {  
	// 省略实现
}

signed main() { 
	int t;
	cin >> t;
	while (t--) {
		maxfac = 0;
		int n;
		cin >> n;
		fac(n);
		if (maxfac == n) {
			cout << "Prime\n";
		}
		else {
			cout << maxfac << endl;
		}
	}
	return 0;
}

此题写法

我们可以先根据雅可比四平方和定理敲出暴力代码：

#include <bits/stdc++.h>
#define int long long
#define ull unsigned long long
using namespace std;
const int MOD = 1e9 + 7;
int solve(int n) {
	int ans = 0;
	for (int i = 1; i * i <= n; i++) {
		if(n % i == 0) {
			if (i % 4 != 0) {
				ans += i;
			}
			if (i * i != n && (n / i) % 4 != 0) {
				ans += (n / i);
			}
			ans %= MOD;
		}
	}
	return 8 * ans % MOD;
}
signed main() {
	int t; 
	cin >> t;
	while (t--) {
		int n;
		cin >> n;
		cout << solve(n) << endl; 
	}
	return 0;
}

直接 60 分拿到手！

那么我们根据前面的 MillerRabin 和 PollardRho，可以写出一下代码：

#include <bits/stdc++.h>
#define int long long
#define ull unsigned long long
using namespace std;
const int MOD = 1e9 + 7;
int qmul(int a, int b, int p) { // 快速乘
	ull c = (ull)a * (ull)b - (ull)((long double)a / p * b + 0.5L) * (ull)p;
	if (c < (ull)p) return c;
	return c + p;
}
int qpow(int a, int b, int p) { // 快速幂
	int res = 1;
	while (b > 0) {
		if (b & 1) {
			res = qmul(res, a, p);
		}
		a = qmul(a, a, p);
		b /= 2;
	}
	return res;
}
bool MillerRabin(int n) {
	if (n < 2) {
		return false;
	}
	if (n == 2 || n == 3) {
		return true;
	}
	int d = n - 1, r = 0;
	while (!(d & 1)) {
		r++;
		d /= 2;
	}
	for (int k = 0; k < 10; ++k) {
		int a = rand() % (n - 2) + 2;
		int x = qpow(a, d, n);
		if (x == 1 || x == n - 1) {
			continue;
		}
		for (int i = 0; i < r - 1; i++) {
			x = qmul(x, x, n);
			if (x == n - 1) {
				break;
			}
		}
		if (x != n - 1) {
			return false;
		}
	}
	return true;
}

int gcd(int a, int b) {
	if (b == 0) {
		return a;
	}
	return gcd(b, a % b);
}
int PollardRho(int x) {
	int s = 0, t = 0, c = rand() % (x - 1) + 1;
	for (int g = 1; ; g *= 2, s = t) {
		int val = 1;
		for (int k = 1; k <= g; k++) {
			t = (qmul(t, t, x) + c) % x;
			val = qmul(val, abs(t - s), x);
			if ((k % 127) == 0) {
				int d = gcd(val, x);
				if (d > 1) {
					return d;
				}
			}
		}
		int d = gcd(val, x);
		if (d > 1) {
			return d;
		}
	}
}
int maxfac;
void fac(int x) {  
	if (x <= maxfac || x < 2) return;
	if (MillerRabin(x)) {             
		maxfac = max(maxfac, x); 
		return;
	}
	int p = x;
	while (p >= x) {
		p = PollardRho(x); 
	}
	while ((x % p) == 0) {
		x /= p;
	}
	fac(x);
	fac(p);
}

int solve_60(int n) {
	int ans = 0;
	for (int i = 1; i * i <= n; i++) {
		if(n % i == 0) {
			if (i % 4 != 0) {
				ans += i;
			}
			if (i * i != n && (n / i) % 4 != 0) {
				ans += (n / i);
			}
			ans %= MOD;
		}
	}
	return 8 * ans % MOD;
}

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;
}

AC！

完结撒花！

纪念我的首黑！