2018冬令营模拟测试赛（十六）

2018冬令营模拟测试赛（十六）

[Problem A]欧拉子图

试题描述

给定一个 (n) 个点，(m) 条边无向图（可能有重边），对于这个图的边集的子集（一共有 (2^m) 个），如果其导出的子图的每个联通块内都存在欧拉回路，我们就把答案加上这个子图的边数的平方，答案对 (10^9+7) 取模。

输入

第一行两个数 (n) 和 (m)。

接下来 (m) 行，每行两个数 (a) 和 (b)，表示 (a) 和 (b) 之间有一条路，(a) 不等于 (b)。

输出

一行输出答案。

输入示例

5 8
1 2
1 3
1 4
3 5
5 4
4 2
3 4
3 5

输出示例

296

数据规模及约定

对于 (10\%) 的数据 (n le 15, m le 20)

对于 (20\%) 的数据 (n le 15, m le 100)

对于 (40\%) 的数据 (n le 60, m le 100)

对于 (60\%) 的数据 (n le 300, m le 300)

对于 (100\%) 的数据 (n, m le 200000)

题解

首先有这么个结论：一张 (n) 个点 (m) 条边的无向连通图的欧拉子图的个数等于 (2^{m-n+1})，即找出任意一个生成树，然后非树边任意组合，每选择一条非树边，就将它所覆盖的链上的所有边的选择情况改变一次；这样就能对于每种非树边的组合给出一个且唯一的对应的欧拉子图的方案了。同理我们可以将这个推广到不连通的无向图去：一张 (n) 个点 (m) 条边的无向图的欧拉子图的个数等于 (2^{m-n+c})，其中 (c) 是连通块个数，因为树边的条数为 (n-c)。

然后对于这道题，我们还需要变个形，考虑它要求什么（令原无向图的边集为 (E)）

[sum_{S in E, S 的生成子图是欧拉子图} |S|^2 \\ = sum_{X in E, E - X 的生成子图是欧拉子图} (|E| - |X|)^2 ]

就是说我们考虑删除哪些边，这样会好统计。

展开一下得到

[sum_{X in E, E - X 的生成子图是欧拉子图} |E|^2 - 2|E| cdot |X| + |X|^2 \\ = |E|^2 cdot 欧拉子图方案数 + 2|E| sum_{X in E, E - X 的生成子图是欧拉子图} |X| + sum_{X in E, E - X 的生成子图是欧拉子图} |X|^2 \\ = |E|^2 cdot 2^{n-m+c} + 2|E| sum_{X in E, E - X 的生成子图是欧拉子图} |X| + sum_{X in E, E - X 的生成子图是欧拉子图} |X|^2 ]

(sum_{X in E, E - X 的生成子图是欧拉子图} |X|) 这个东西其实可以换个顺序求（以下的 (e)、(f) 都代表一条无向边）：

[sum_{X in E, E - X 的生成子图是欧拉子图} |X| \\ = sum_{X in E, E - X 的生成子图是欧拉子图} sum_{e in X} 1 \\ = sum_{e in E} sum_{e in X, X in E, E - x 的生成子图是欧拉子图} 1 ]

可以看到我只需要求出对于每条边有多少删边方案包含它就可以了。现在讨论边 (e) 在多少删除方案中，如果 (e) 是桥边，那么删掉它会导致连通块数量必定 (+1)，同时总边数 (-1)，所以方案数还是 (2^{m-n+c})（事实上欧拉子图中不可能包含桥边，所以所有的方案都会删掉桥边）；若 (e) 不是桥边，那么删掉它会导致边数 (-1) 连通块数不变，即 (2^{m-1-n+c})。

对于平方项我们也可以如此求

[sum_{X in E, E - X 的生成子图是欧拉子图} |X|^2 \\ = sum_{X in E, E - X 的生成子图是欧拉子图} sum_{e in X} sum_{f in X} 1 \\ = sum_{e, f in E} sum_{X in E, E - X 的生成子图是欧拉子图} 1 ]

就是要求对于任意两条边 (e, f)，同时包含这两条边的方案数。讨论若两条边均为桥，则连通块数量 (+2)，边数 (-2)，方案数为 (2^{m-n+c})；若其中一条边为桥，则连通块数量 (+1)，边数 (-2)，方案数为 (2^{m-n+c-1})；若都不是桥，那么连通块数量不一定不变，如果两条边被一模一样的非树边覆盖（特别地，非树边只被自己覆盖），删掉它们会导致连通块数目 (+1)，否则连通块数目不变，当联通卡数目不变时，方案数为 (2^{m-n+c-2})。

如何判断是否被一模一样的非树边覆盖呢？给每条非树边随机一个权值，一条树边的权值就是将所有覆盖它的非树边的权值异或起来的值，若两条边权值相等则认为被一模一样的非树边覆盖。用 unsigned long long 存权值，错误率小到可以 AC。（但是如果用 int 错误率就大到会 WA 两三个点）

#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cctype>
#include <algorithm>
using namespace std;
#define rep(i, s, t) for(int i = (s), mi = (t); i <= mi; i++)
#define dwn(i, s, t) for(int i = (s), mi = (t); i >= mi; i--)

int read() {
	int x = 0, f = 1; char c = getchar();
	while(!isdigit(c)){ if(c == '-') f = -1; c = getchar(); }
	while(isdigit(c)){ x = x * 10 + c - '0'; c = getchar(); }
	return x * f;
}

#define maxn 200010
#define maxm 400010
#define MOD 1000000007
#define LL long long
#define UL unsigned long long

UL RNDUL() {
	UL res = 0;
	rep(i, 0, 7) res = res << 8 | (rand() & 255);
	return res;
}

int n, m, head[maxn], nxt[maxm], to[maxm], id[maxm], pow2[maxn];
UL val[maxn];

void AddEdge(int a, int b, int i) {
	to[++m] = b; id[m] = i; nxt[m] = head[a]; head[a] = m;
	swap(a, b);
	to[++m] = b; id[m] = i; nxt[m] = head[a]; head[a] = m;
	return ;
}

int clo, dfn[maxn], low[maxn], bcno[maxn], cntb, S[maxn], top, con, bridge;
UL Xor[maxn];
bool isbridge[maxn], isback[maxn];
void dfs(int u, int fai) {
	dfn[u] = low[u] = ++clo;
	S[++top] = u;
	// printf("in: %d
", u);
	for(int e = head[u]; e; e = nxt[e])if(id[e] != fai) {
		if(dfn[to[e]]) {
			low[u] = min(low[u], dfn[to[e]]);
			if(dfn[to[e]] < dfn[u]) val[id[e]] = RNDUL(), Xor[u] ^= val[id[e]], Xor[to[e]] ^= val[id[e]], isback[id[e]] = 1;
		}
		else dfs(to[e], id[e]), low[u] = min(low[u], low[to[e]]);
	}
	// printf("out: %d %d
", u, Xor[u]);
	if(low[u] == dfn[u]) {
		cntb++;
		while(S[top] != u) bcno[S[top--]] = cntb;
		bcno[S[top--]] = cntb;
		isbridge[u] = 1;
	}
	return ;
}
void dfs2(int u, int fa) {
	for(int e = head[u]; e; e = nxt[e]) if(to[e] != fa && !isback[id[e]]) dfs2(to[e], u), Xor[u] ^= Xor[to[e]];
	return ;
}

int Pow2(int a) {
	return a < 0 ? 0 : pow2[a];
}

UL getXor[maxn];
int cntg;

int main() {
	srand(20162109);
	
	n = read(); int m = read();
	pow2[0] = 1;
	rep(i, 1, m) {
		pow2[i] = pow2[i-1] << 1; if(pow2[i] >= MOD) pow2[i] -= MOD;
		int a = read(), b = read();
		AddEdge(a, b, i);
	}
	
	rep(i, 1, n) if(!bcno[i]) con++, dfs(i, 0), dfs2(i, 0);
	bridge = cntb - con;
	
	int ans = (LL)m * m % MOD * Pow2(m - n + con) % MOD, tmp = 0;
	ans -= (LL)(m << 1) * ((LL)bridge * Pow2(m - n + con) % MOD + (LL)(m - bridge) * Pow2(m - 1 - n + con) % MOD) % MOD;
	if(ans < 0) ans += MOD;
	tmp += ((LL)bridge * (bridge - 1) >> 1) % MOD * Pow2(m - n + con) % MOD; if(tmp >= MOD) tmp -= MOD;
	tmp += (LL)bridge * (m - bridge) % MOD * Pow2(m - 1 - n + con) % MOD; if(tmp >= MOD) tmp -= MOD;
	rep(i, 1, n) if(!isbridge[i]) getXor[++cntg] = Xor[i];
	rep(i, 1, m) if(isback[i]) getXor[++cntg] = val[i];
	sort(getXor + 1, getXor + cntg + 1);
	int cnow = 0;
	rep(i, 1, cntg) {
		cnow++;
		if(i == cntg || getXor[i] != getXor[i+1]) {
			tmp += (LL)(cntg - i) * cnow % MOD * Pow2(m - 2 - n + con) % MOD; if(tmp >= MOD) tmp -= MOD;
			tmp += ((LL)cnow * (cnow - 1) >> 1) % MOD * Pow2(m - 1 - n + con) % MOD; if(tmp >= MOD) tmp -= MOD;
			cnow = 0;
		}
	}
	tmp <<= 1; if(tmp >= MOD) tmp -= MOD;
	tmp += ((LL)bridge * Pow2(m - n + con) + (LL)(m - bridge) * Pow2(m - 1 - n + con)) % MOD; if(tmp >= MOD) tmp -= MOD;
	ans += tmp; if(ans >= MOD) ans -= MOD;
	
	printf("%d
", ans);
	
	return 0;
}

[Problem B]Ball

试题描述

一行有 (n) 个球，现在将这些球分成 (K) 组，每组可以有一个球或相邻两个球。一个球只能在至多一个组中（可以不在任何组中）。求对于 (1 le K le m) 的所有 (K) 分别有多少种分组方法。答案对 (998244353) 取模。

输入

一行两个正整数 (n)，(m)。

输出

一行 (m) 个整数，第 (i) 个整数表示组数为 (i) 时的答案。

输入示例1

3 3

输出示例1

5 5 1

输入示例2

5 10

输出示例2

9 25 25 9 1 0 0 0 0 0

数据规模及约定

测试点编号	$n$	$m$
1	$leq 10$	$< 2^3$
2	$leq 1000$	$< 2^{10}$
3	$leq 10^9$
4		$< 2^{12}$
5		$< 2^{14}$
6		$< 2^{15}$
7		$< 2^{16}$
8,9,10		$< 2^{17}$

题解

暴力做的话就是一个 dp，令 (f(i, j)) 表示 (i) 个球分 (j) 组的方案数。转移如下

[f(i, j) = f(i-1, j) + f(i-1, j-1) + f(i-2, j-1) ]

如果我们搞出 (f(i, j)) 生成函数 (F_i(x) = sum_{j=0}^m f(i, j)x^j)，那么上面的转移可以写成如下形式

[egin{matrix} F_i(x) equiv F_{i-1}(x) + xF_{i-1}(x) + xF_{i-2}(x) & (mathrm{mod} x^{m+1}) \\ F_i(x) equiv (x+1) F_{i-1}(x) + xF_{i-2}(x) & (mathrm{mod} x^{m+1}) end{matrix} ]

最后要求的就是 (F_n(x))。下面略去所有多项式的 ((x))

由于 (n) 很大，转化成多项式之后又是一个递推式，自然想到矩阵快速幂，把转移矩阵构造出来是这样的

[left( egin{matrix} 0 & 1 \ x & x + 1 end{matrix} ight) ^n left( egin{matrix} F_0 \ F_1 end{matrix} ight) = left( egin{matrix} F_n \ F_{n+1} end{matrix} ight) ]

然后我们可以求出这个矩阵的特征值（就是特征方程的解，特征方程就是特征多项式 (p_A(lambda) = |lambda I - A|) 等于 (0)，(A) 为该矩阵），它的特征矩阵如下

[lambda I - A = left( egin{matrix} lambda & 0 \ 0 & lambda end{matrix} ight) - left( egin{matrix} 0 & 1 \ x & x + 1 end{matrix} ight) = left( egin{matrix} lambda & -1 \ -x & lambda - x - 1 end{matrix} ight) ]

于是特征方程就是

[left| egin{matrix} lambda & -1 \ -x & lambda - x - 1 end{matrix} ight| = lambda(lambda-x-1) - (-x)(-1) = lambda^2 - (x+1)lambda - x = 0 ]

解得 (lambda = frac{x+1 pm sqrt{x^2+6x+1}}{2})。

然后可以将这个转移矩阵 (A) 对角化，即找到可逆的矩阵 (P) 和对角矩阵 (B)，使得 (P^{-1}BP = A)，其中那个对角矩阵是

[left( egin{matrix} lambda_1 & 0 \ 0 & lambda_2 end{matrix} ight) ]

可以证明能够这样操作，因为矩阵 (A) 有两个不同的特征值，那么就必然有两个线性无关的特征向量（详见《线性代数》）。

同时还可以证明

[A^n = P^{-1}B^nP = P^{-1} left( egin{matrix} lambda_1^n & 0 \ 0 & lambda_2^n end{matrix} ight) P ]

这是由于 (A^n) 的特征值是 (lambda^n)。

然后原矩阵快速幂就可以转化成下面这种形式

[left( egin{matrix} F_n \ F_{n+1} end{matrix} ight) = left( egin{matrix} 0 & 1 \ x & x + 1 end{matrix} ight) ^n left( egin{matrix} F_0 \ F_1 end{matrix} ight) = P^{-1} left( egin{matrix} lambda_1^n & 0 \ 0 & lambda_2^n end{matrix} ight) P left( egin{matrix} F_0 \ F_1 end{matrix} ight) ]

我们不用关心 (P) 具体是什么，因为我们知道矩阵乘法永远是线性变换，所以我们知道 (F_n = c_1 lambda_1^n + c_2 lambda_2^n)，现在的任务是解出 (c_1) 和 (c_2)。

我们将 (F_0) 和 (F_1) 分别代入上式就好啦。

[F_0 = 1 = c_1 + c_2 \\ F_1 = 1 + x = c_1 lambda_1 + c_2 lambda_2 ]

解得

[c_1 = frac{x+1-lambda_2}{lambda_1-lambda_2} \\ c_2 = frac{lambda_1-x-1}{lambda_1-lambda_2} ]

不妨令 (lambda_1 = frac{x+1 + sqrt{x^2+6x+1}}{2})，(lambda_2 = frac{x+1 - sqrt{x^2+6x+1}}{2})，你发现 (lambda_2) 的常数项为 (0)，所以 (lambda_2^n equiv 0(mathrm{mod} x^{m+1}))（注意有意义的 (m) 总是 (le n)，对于 (> n) 的部分直接输出 (0) 就好了）。这告诉我们 (F_n equiv c_1 lambda_1^n(mathrm{mod} x^{m+1}))。

带进去，把 (c_1) 彻底求出来

[c_1 = frac{x+1-frac{x+1-sqrt{x^2+6x+1}}{2}}{sqrt{x^2+6x+1}} ]

于是只需要写一个多项式逆元、开方、exp、ln、乘方就好啦（exp 和 ln 在算乘方时要用）。

然后我们发现乘方的对象 (lambda_1) 的常数项是 (1)，于是问题又简化了许多呢！

#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cctype>
#include <algorithm>
using namespace std;
#define rep(i, s, t) for(int i = (s), mi = (t); i <= mi; i++)
#define dwn(i, s, t) for(int i = (s), mi = (t); i >= mi; i--)
#define LL long long

LL read() {
	LL x = 0, f = 1; char c = getchar();
	while(!isdigit(c)){ if(c == '-') f = -1; c = getchar(); }
	while(isdigit(c)){ x = x * 10 + c - '0'; c = getchar(); }
	return x * f;
}

#define maxn 262144
#define Groot 3
#define MOD 998244353

int Pow(int a, int b) {
	int ans = 1, t = a;
	while(b) {
		if(b & 1) ans = (LL)ans * t % MOD;
		t = (LL)t * t % MOD; b >>= 1;
	}
	return ans;
}

int brev[maxn];
void FFT(int *a, int len, int tp) {
	int n = 1 << len;
	rep(i, 0, n - 1) if(i < brev[i]) swap(a[i], a[brev[i]]);
	rep(i, 1, len) {
		int wn = Pow(Groot, MOD - 1 >> i);
		if(tp < 0) wn = Pow(wn, MOD - 2);
		for(int j = 0; j < n; j += 1 << i) {
			int w = 1;
			rep(k, 0, (1 << i >> 1) - 1) {
				int la = a[j+k], ra = (LL)w * a[j+k+(1<<i>>1)] % MOD;
				a[j+k] = (la + ra) % MOD;
				a[j+k+(1<<i>>1)] = (la - ra + MOD) % MOD;
				w = (LL)w * wn % MOD;
			}
		}
	}
	if(tp < 0) {
		int invn = Pow(n, MOD - 2);
		rep(i, 0, n - 1) a[i] = (LL)a[i] * invn % MOD;
	}
	return ;
}
void Mul(int *A, int *B, int n, int m) {
	int N = 1, len = 0;
	while(N <= n + m) N <<= 1, len++;
	rep(i, 0, N - 1) brev[i] = (brev[i>>1] >> 1) | ((i & 1) << len >> 1);
	rep(i, n + 1, N - 1) A[i] = 0; rep(i, m + 1, N - 1) B[i] = 0;
	FFT(A, len, 1); FFT(B, len, 1);
	rep(i, 0, N - 1) A[i] = (LL)A[i] * B[i] % MOD;
	FFT(A, len, -1);
	return ;
}

int tmp[maxn];
void inverse(int *f, int *g, int n) {
	if(n == 1) return (void)(f[0] = Pow(g[0], MOD - 2));
	inverse(f, g, n + 1 >> 1);
	int N = 1, len = 0;
	while(N <= (n << 1)) N <<= 1, len++;
	rep(i, 0, N - 1) brev[i] = (brev[i>>1] >> 1) | ((i & 1) << len >> 1);
	rep(i, 0, n - 1) tmp[i] = g[i]; rep(i, n, N - 1) tmp[i] = f[i] = 0;
	FFT(f, len, 1); FFT(tmp, len, 1);
	rep(i, 0, N - 1) f[i] = (2ll - (LL)tmp[i] * f[i] % MOD + MOD) * f[i] % MOD;
	FFT(f, len, -1);
	rep(i, n, N - 1) f[i] = 0;
	return ;
}

int t1[maxn], t2[maxn], t3[maxn];
void p_sqrt(int *f, int *g, int n) { // g[0] = 1;
	if(n == 1) return (void)(f[0] = 1);
	p_sqrt(f, g, n + 1 >> 1);
	rep(i, 0, n - 1) t1[i] = g[i];
	rep(i, 0, n - 1) t2[i] = (f[i] << 1) % MOD;
	inverse(t3, t2, n);
	Mul(t1, t3, n - 1, n - 1);
	int inv2 = Pow(2, MOD - 2); rep(i, 0, n - 1) t2[i] = (LL)f[i] * inv2 % MOD;
	rep(i, 0, n - 1) f[i] = (t1[i] + t2[i]) % MOD;
	return ;
}

int inv[maxn];
void getinv(int n) {
	inv[1] = 1;
	rep(i, 2, n) inv[i] = (LL)(MOD - MOD / i) * inv[MOD%i] % MOD;
	return ;
}

void logarithm(int *f, int *g, int n) { // g[0] = 1;
	rep(i, 1, n - 1) t1[i-1] = (LL)g[i] * i % MOD;
	inverse(t2, g, n);
	Mul(t1, t2, n - 2, n - 1);
	dwn(i, n - 1, 0) f[i+1] = (LL)t1[i] * inv[i+1] % MOD; f[0] = 0;
	return ;
}

void exponential(int *f, int *g, int n) { // g[0] = 0;
	if(n == 1) return (void)(f[0] = 1);
	exponential(f, g, n + 1 >> 1);
	logarithm(t1, f, n);
	rep(i, 0, n - 1) t1[i] = (g[i] - t1[i] + MOD) % MOD; t1[0]++; if(t1[0] >= MOD) t1[0] -= MOD;
	int N = 1; while(N <= (n - 1 << 1)) N <<= 1;
	Mul(f, t1, n - 1, n - 1);
	rep(i, n, N - 1) f[i] = 0;
	return ;
}

void p_pow(int *f, int *g, int k, int n) { // g[0] = 1;
	logarithm(t3, g, n);
	rep(i, 0, n - 1) t3[i] = (LL)t3[i] * k % MOD;
	exponential(f, t3, n);
	return ;
}

int A[maxn], B[maxn], tB[maxn], lmd[maxn], lmdn[maxn];
int main() {
	LL n = read(); int om = read(), m = min((LL)om, n);
	getinv(m + 2);
	
	tB[0] = tB[2] = 1; tB[1] = 6;
	p_sqrt(B, tB, m + 1);
	rep(i, 0, m) A[i] = MOD - B[i];
	A[0]++; A[1]++; if(A[0] >= MOD) A[0] -= MOD; if(A[1] >= MOD) A[1] -= MOD;
	rep(i, 0, m) A[i] = MOD - ((LL)A[i] * inv[2] % MOD);
	A[0]++; A[1]++; if(A[0] >= MOD) A[0] -= MOD; if(A[1] >= MOD) A[1] -= MOD;
	
	inverse(tB, B, m + 1);
	Mul(A, tB, m, m); rep(i, m + 1, m << 1) A[i] = 0;
	
	rep(i, 0, m) lmd[i] = B[i];
	lmd[0]++; lmd[1]++; if(lmd[0] >= MOD) lmd[0] -= MOD; if(lmd[1] >= MOD) lmd[1] -= MOD;
	rep(i, 0, m) lmd[i] = (LL)lmd[i] * inv[2] % MOD;
	p_pow(lmdn, lmd, n % MOD, m + 1);
	
	Mul(A, lmdn, m, m);
	
	rep(i, 1, m) printf("%d%c", A[i], i < om ? ' ' : '
');
	rep(i, m + 1, om) printf("%d%c", 0, i < om ? ' ' : '
');
	
	return 0;
}

[Problem C]Div

试题描述

定义复数 (a+bi) 为整数 (k) 的约数，当且仅当 (a) 和 (b) 为整数且存在整数 (c) 和 (d) 满足 ((a+bi)(c+di)=k)

定义复数 (a+bi) 的实部为 (a)，虚部为 (b)

定义 (f(n)) 为整数 (n) 的所有实部大于 (0) 的约数的实部之和

给定正整数 (n)，求出 (sum_{i=1}^n f(i))

对 (1004535809) 取模后得到的值

即求

[sum_{i=1}^n sum_{z in C, z|i, Re[z] > 0} Re[z] ]

（(C) 表示复数集合，(Re[z]) 表示 (z) 的实部，(Im[z]) 表示 (z) 的虚部）

输入

一行一个正整数 (n)

输出

一行一个整数表示答案

输入示例1

5

输出示例1

35

输入示例2

1000

输出示例2

1752541

输入示例3

1000000

输出示例3

636408476

数据规模及约定

数据编号	$n$
1	$le 10$
2	$le 100$
3	$le 200$
4	$le 300$
5	$le 400$
6	$le 500$
7	$le 3000$
8	$le 5000$
9	$le 3 imes 10^6$
10	$le 5 imes 10^6$
11	$le 7 imes 10^6$
12	$le 10^7$
13	$le 3 imes 10^8$
14	$le 5 imes 10^8$
15	$le 7 imes 10^8$
16	$le 10^9$
17	$le 3 imes 10^9$
18	$le 5 imes 10^9$
19	$le 7 imes 10^9$
20	$le 10^{10}$

题解

又是一道大大的反演 + 杜教筛题……

我们就是要求满足 ((a+bi)(c+di) le n(a, c > 0, (a+bi)(c+di) 是实数)) 的所有的 (a) 的和。

就是 (ac-bd le n, ad+bc = 0)，移项发现 (b e 0) 时，(frac{a}{b} = -frac{c}{d})，我们设

[frac{p}{q} = frac{a}{b} = -frac{c}{d} ]

其中 (p, q > 0, mathrm{gcd}(p,q) = 1)，由于我们忽略了 (b = 0) 的情况，且虚部可正可负，所以最后答案需要 ( imes 2)，再加上虚部等于 (0) 的情况（就是 (1 sim n) 的所有数的约数和）。

那么我们现在要求这个东西

[sum_{k=1}^n sum_{p^2+q^2=k} { [mathrm{gcd}(p,q)=1] sum_{x,y in [1, n]} [kxy le n]px } \\ = sum_{k=1}^n sum_{p^2+q^2=k} { [mathrm{gcd}(p,q)=1]p sum_{x=1}^{lfloor frac{n}{k} floor} x lfloor frac{n}{kx} floor } ]

令 (g(n) = sum_{x=1}^n x lfloor frac{n}{x} floor)，其实 (g(n)) 也是 (1 sim n) 的所有约数的和；再令 (t(n) = sum_{p^2+q^2=n} [mathrm{gcd}(p,q)=1]p)，式子变成如下形式

[sum_{k=1}^n t(k)g(lfloor frac{n}{k} floor) ]

分块求它。现在考虑求 (T(n) = sum_{i=1}^n t(i))，可以用杜教筛（注意，杜教筛是一个广泛使用的词，只要满足可以预处理前 (n^{frac{2}{3}})，后面较大的部分能够 (sqrt{n}) 求出即可满足最终复杂度是 (O(n^{frac{2}{3}}))，就能叫做杜教筛；有了这个性质，杜教筛可以任意套杜教筛，复杂度不变，待会会看到，求 (T(n)) 的超过 (n^{frac{2}{3}}) 部分的过程中可以再套一个杜教筛）。

现在考虑求 (T(n))

[T(n) \\ = sum_{k=1}^n t(k) \\ = sum_{k=1}^n sum_{p^2+q^2=k} [mathrm{gcd}(p,q)=1]p \\ = sum_{d=1}^{lfloor sqrt{n} floor} { mu(d) sum_{k=1}^n sum_{d|p} sum_{d|q} [p^2+q^2=k]p } \\ = sum_{d=1}^{lfloor sqrt{n} floor} { mu(d) sum_{d|p} sum_{d|q} [p^2+q^2 le n]p } \\ = sum_{d=1}^{lfloor sqrt{n} floor} { mu(d) sum_{p=1}^{lfloor frac{sqrt{n}}{d} floor} sum_{q=1}^{lfloor frac{sqrt{n}}{d} floor} [p^2+q^2 le lfloor frac{n}{d^2} floor]pd } \\ = sum_{d=1}^{lfloor sqrt{n} floor} { mu(d) cdot d sum_{p=1}^{lfloor sqrt{frac{n}{d^2}} floor} p lfloor sqrt{lfloor frac{n}{d^2} floor - p^2} floor } ]

然后这个 (d) 可以直接暴力枚举，做到 (O(sqrt{n}))，接下来令 (h(n) = sum_{p=1}^{lfloor sqrt{n} floor} p lfloor sqrt{n - p^2} floor)，那么

[T(n) = sum_{d=1}^{lfloor sqrt{n} floor} { mu(d) cdot d cdot h(lfloor frac{n}{d^2} floor) } ]

(h(n)) 如果暴力求的话，最后复杂度会带一个 (ln n)，因为它的复杂度是 (int_0^{sqrt{n}} sqrt{frac{n}{x^2}} = sqrt{n} ln sqrt{n} = O(sqrt{n} ln n))（注意这里分析的是求 (T(n)) 超出 (n^{frac{2}{3}}) 的部分的复杂度），要想不带这个 (ln n)，就得用 (h(n) = sum_{p=1}^{lfloor frac{sqrt{n}}{d} floor} sum_{q=1}^{lfloor frac{sqrt{n}}{d} floor} [p^2+q^2 le lfloor frac{n}{d^2} floor]p) 这个形式预处理 (h(n)) 的前 (n^{frac{2}{3}}) 项（这就是我说的杜教筛套杜教筛了）；当然不预处理也能过。

【好像做完了？】

#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cctype>
#include <algorithm>
#include <cmath>
using namespace std;
#define rep(i, s, t) for(int i = (s), mi = (t); i <= mi; i++)
#define dwn(i, s, t) for(int i = (s), mi = (t); i >= mi; i--)
#define LL long long

LL read() {
	LL x = 0, f = 1; char c = getchar();
	while(!isdigit(c)){ if(c == '-') f = -1; c = getchar(); }
	while(isdigit(c)){ x = x * 10 + c - '0'; c = getchar(); }
	return x * f;
}

#define maxn 4700000
#define maxsq 2200
#define MOD 1004535809

bool vis[maxn];
int prime[maxn], cp, mu[maxn], _sigma[maxn], sigma[maxn], sum[maxn], funt[maxn], sumt[maxn], gcd[maxsq][maxsq];
void init() {
	int n = maxn - 1;
	mu[1] = sigma[1] = sum[1] = 1;
	rep(i, 2, n) {
		if(!vis[i]) prime[++cp] = i, mu[i] = -1, sigma[i] = 1 + i, _sigma[i] = 1;
		for(int j = 1; j <= cp && i * prime[j] <= n; j++) {
			vis[i*prime[j]] = 1;
			if(i % prime[j] == 0) {
				sigma[i*prime[j]] = ((LL)sigma[i] * prime[j] + _sigma[i]) % MOD;
				_sigma[i*prime[j]] = _sigma[i];
				mu[i*prime[j]] = 0;
				break;
			}
			sigma[i*prime[j]] = (LL)sigma[i] * (prime[j] + 1) % MOD;
			_sigma[i*prime[j]] = sigma[i];
			mu[i*prime[j]] = -mu[i];
		}
		sum[i] = sum[i-1] + sigma[i]; if(sum[i] >= MOD) sum[i] -= MOD;
	}
	rep(i, 1, maxsq - 1) gcd[i][0] = gcd[0][i] = i;
	rep(i, 1, maxsq - 1) rep(j, 1, maxsq - 1) if(i * i + j * j < maxn) {
		gcd[i][j] = gcd[j][i] = gcd[i][j%i];
		if(gcd[i][j] == 1) (funt[i*i+j*j] += i) %= MOD;
	}
	rep(i, 1, n) {
		sumt[i] = sumt[i-1] + funt[i];
		if(sumt[i] >= MOD) sumt[i] -= MOD;
	}
	return ;
}

const int HMOD = 1000037, maxtot = 60000;
struct Hash {
	int ToT, head[HMOD], nxt[maxtot], val[maxtot];
	LL key[maxtot];
	int Find(LL x) {
		int u = x % HMOD;
		for(int e = head[u]; e; e = nxt[e]) if(key[e] == x) return val[e];
		return -1;
	}
	void Insert(LL x, int v) {
		int u = x % HMOD;
		nxt[++ToT] = head[u]; key[ToT] = x; val[ToT] = v; head[u] = ToT;
		return ;
	}
} gh, Th, hh;

int g(LL n) {
	if(n < maxn) return sum[n];
	int ans = gh.Find(n);
	if(ans >= 0) return ans;
	ans = 0;
	for(LL i = 1; i <= n; ) {
		LL r = min(n / (n / i), n);
		ans += ((__int128)(i + r) * (r - i + 1) >> 1) % MOD * (n / i) % MOD;
		if(ans >= MOD) ans -= MOD;
		i = r + 1;
	}
	gh.Insert(n, ans);
	return ans;
}

int h(LL n) {
	int ans = hh.Find(n);
	if(ans >= 0) return ans;
	ans = 0;
	int m = sqrt(n + .5);
	rep(i, 1, m) {
		ans += (LL)i * (LL)sqrt(n - (LL)i * i) % MOD;
		if(ans >= MOD) ans -= MOD;
	}
	hh.Insert(n, ans);
	return ans;
}

int T(LL n) {
	if(n < maxn) return sumt[n];
	int ans = Th.Find(n);
	if(ans >= 0) return ans;
	ans = 0;
	int m = sqrt(n + .5);
	rep(i, 1, m) {
		ans += ((LL)mu[i] * i + MOD) * h(n / ((LL)i * i)) % MOD;
		if(ans >= MOD) ans -= MOD;
	}
	Th.Insert(n, ans);
	return ans;
}

int main() {
	init();
	LL n = read();
	
	int ans = 0;
	for(LL i = 1; i <= n; ) {
		LL r = min(n / (n / i), n);
		ans += (LL)(T(r) - T(i - 1) + MOD) * g(n / i) % MOD;
		if(ans >= MOD) ans -= MOD;
		i = r + 1;
	}
	// printf("ans: %d
", ans);
	ans <<= 1; if(ans >= MOD) ans -= MOD;
	// printf("g(n): %d
", g(n));
	ans += g(n); if(ans >= MOD) ans -= MOD;
	
	printf("%d
", ans);
	
	return 0;
}