CF449E Jzzhu and Squares

题目大意：有一个$n imes m$的方格图，求其中所有的格点正方形完整包含的小方格个数，多组询问。$n,mleqslant 10^6$

题解：令$nleqslant m$。有一个显然的式子：
$$
ans=sumlimits_{i=1}^n(n-i+1)(m-i+1)f(i)
$$
$f(i)$表示可以完整包含在$i imes i$的正方形中且顶点在这个正方形边上的正方形所包含的小方格总数。可以分选的正方形和$i imes i$正方形重合和边转动$j$格来计算
$$
f(n)=n^2+sumlimits_{i=1}^{n-1}[(n-2max{i,n-1})^2+g(i,n-i)]
$$
其中$(n-2max{i,n-1})^2$是转动$i$格后正中间的完整小方格，$g(i,n-i)$表示四周的$4$个小直角三角形中包含的完整小方格个数。

比如计算左上角的直角三角形，发现完整的小方格的左上角满足在直角三角形的斜边上或三角形内。三角形内的点可以用皮克定理来解决，令三角形两直角边为$n,m$，则三角形内的点个数为$dfrac{nm-n-m+2-gcd(n,m)}2$，再加上在斜边上的点，则一个直角三角形内的方格数为$dfrac{nm-m-n+gcd(n,m)}2$。

把$n=i,m=n-i$带入式子
$$
egin{align*}
f(n)=&n^2+sumlimits_{i=1}^{n-1}[(n-2max{i,n-1})^2+g(i,n-i)]\
    =&n^2+sumlimits_{i=1}^{n-1}(n-2max{i,n-1})^2+sumlimits_{i=1}^{n-1}g(i,n-i)\
    =&n^2+sumlimits_{i=1}^{leftlfloorfrac{n-1}2 ight floor}(n-2i)^2+sumlimits_{i=leftlfloorfrac{n-1}2 ight floor+1}^{n-1}(2i-n)^2\
    &+sumlimits_{i=1}^{n-1}dfrac{i(n-i)-(n-i)-i+gcd(i,n-i)}2\
    =&n^2+2sumlimits_{i=1}^{leftlfloorfrac{n-1}2 ight floor}(n-2i)^2+sumlimits_{i=1}^{n-1}dfrac{in-i^2-n+gcd(i,n)}2
end{align*}
$$
其中$n^2$可以快速计算，$sumlimits_{i=1}^{leftlfloorfrac{n-1}2 ight floor}(n-2i)^2$可以前缀和解决，问题在$sumlimits_{i=1}^{n-1}dfrac{in-i^2-n+gcd(i,n)}2$部分。令$h(n)=sumlimits_{i=1}^{n-1}dfrac{in-i^2-n+gcd(i,n)}2$，$sgcd(n)=sumlimits_{i=1}^{n-1}gcd(i,n)$
$$
h(n)=dfrac12[sumlimits_{i=1}^{n-1}(in-i^2-n)+sgcd(n)]\
egin{align*}
h(n-1)&=dfrac12[sumlimits_{i=1}^{n-2}(i(n-1)-i^2-(n-1))+sgcd(n-1)]\
    &=dfrac12[sumlimits_{i=1}^{n-2}(in-i^2-n-i+1)+sgcd(n-1)]\
end{align*}\
$$
$$
egin{align*}
h(n)=&h(n-1)+dfrac12[sumlimits_{i=1}^{n-2}(i-1)+(n-1)n-(n-1)^2-n\
    &-sgcd(n-1)+sgcd(n)]\
    =&h(n-1)+dfrac12[left(sumlimits_{i=0}^{n-3}i ight)-1-sgcd(n-1)+sgcd(n)]\
    =&h(n-1)+dfrac12[left(dfrac{(n-2)(n-3)}2 ight)-1-sgcd(n-1)+sgcd(n)]
end{align*}
$$

其他部分都可以快速求出，问题在求$sgcd(n)$
$$
egin{align*}
sgcd(n)&=sumlimits_{i=1}^{n-1}gcd(i,n)\
    &=left(sumlimits_{d|n}dfrac ndvarphi(d) ight)-n
end{align*}
$$
全部预处理出来即可。

卡点：无

C++ Code：

#include <cstdio>
#include <algorithm>
#include <iostream>
#define mul(a, b) (static_cast<long long> (a) * (b) % mod)
const int mod = 1e9 + 7, maxn = 1e6 + 10, half = (mod + 1) / 2;
inline void reduce(int &x) { x += x >> 31 & mod; }

int gcd(int a, int b) {
	if (!b) return a;
	return gcd(b, a % b);
}
inline int sqr(int x) { return mul(x, x); }

int g[maxn], sumgcd[maxn], phi[maxn], plist[maxn / 2], ptot;
int pre[maxn], R[maxn], T[maxn], preF[maxn];
bool notp[maxn];
int Q;

int F(int n) {
	int ans = g[n];
	reduce(ans += sqr(n) - mod);
	reduce(ans += pre[n - 2] - mod);
	reduce(ans += pre[n - 2] - mod);
	return ans;
}
void sieve(int N) {
	phi[1] = 1;
	for (int i = 2; i <= N; i++) {
		if (!notp[i]) phi[plist[ptot++] = i] = i - 1;
		for (int j = 0, t; j < ptot && (t = i * plist[j]) <= N; j++) {
			notp[t] = true;
			if (i % plist[j] == 0) {
				phi[t] = phi[i] * plist[j];
				break;
			}
			phi[t] = phi[i] * phi[plist[j]];
		}
	}
	for (int i = 1; i <= N; ++i) {
		for (int j = i + i; j <= N; j += i) reduce(sumgcd[j] += mul(phi[j / i], i));
	}
	for (int i = 4; i <= N; ++i) {
		g[i] = (1ll * (i - 3) * (i - 2) / 2 - 1 - sumgcd[i - 1] + sumgcd[i]) % mod;
		reduce(g[i]);
		g[i] = mul(g[i], half);
		reduce(g[i] += g[i - 1] - mod);
	}
	for (int i = 1; i <= N; ++i) g[i] = mul(g[i], 4);
	pre[1] = 1;
	for (int i = 2; i <= N; ++i) reduce(pre[i] = pre[i - 2] + sqr(i) - mod);
	for (int i = 1; i <= N; ++i) reduce(preF[i] = preF[i - 1] + F(i) - mod);
	for (int i = 1; i <= N; ++i) {
		reduce(R[i] = R[i - 1] + preF[i - 1] - mod);
		reduce(R[i] += F(i) - mod);
	}
	for (int i = 1; i <= N; ++i) {
		reduce(T[i] = T[i - 1] + preF[i - 1] - mod);
		reduce(T[i] += R[i - 1] - mod);
		reduce(T[i] += R[i - 1] - mod);
		reduce(T[i] += F(i) - mod);
	}
}
int solve(int n, int m) {
	int ans = T[n];
	reduce(ans += mul(R[n], m - n) - mod);
	return ans;
}

int main() {
	std::ios::sync_with_stdio(false), std::cin.tie(0), std::cout.tie(0);
	sieve(1000005);
	std::cin >> Q;
	while (Q --> 0) {
		static int n, m;
		std::cin >> n >> m;
		if (n > m) std::swap(n, m);
		std::cout << solve(n, m) << '
';
	}
	return 0;
}