2021 ccpc 新疆省赛

Problem A. balloon

将两个序列从小到大排序，从跳跃高度较低的人考虑到跳跃高度较高的人，使用一个变量 ( ext{now})​​​ 维护当前取到的气球，若气球能取就取，不能取就换成下一个人。

时间复杂度的瓶颈在于排序。

#include <iostream>
#include <algorithm>
#include <cstring>
#include <cstdio>
using namespace std;
const int N = 2e5 + 10;
int n, m;
struct Node {
	int h;
	int id;
} e[N];
int ans[N];
int h[N];
bool cmp(Node a, Node b) {
	return a.h < b.h;
}
void solve() {
	cin >> n >> m;
	for(int i = 1; i <= n ; i ++) {
		cin >> e[i].h;
		e[i].id = i;
	}
	sort(e + 1, e + 1 + n, cmp);

	for(int i = 1; i <= m; i ++)
		cin >> h[i];
	int now  = 1;
	sort(h + 1, h + 1 + m);
	for(int i = 1; i <= n; i ++) {
		int num = 0;
		while(now <= m && h[now] <= e[i].h)
			num ++, now ++;
		ans[e[i].id] = num ;
	}
	for(int i = 1; i <= n; i ++)
		cout << ans[i] << endl;
	return ;
}
signed main() {
	int t = 1;
	while(t--)
		solve();
	return 0;
}

Problem B. sophistry

可以从后往前倒序 dp。

设 (f_i) 表示：考虑到了第 (i) 天到第 (n) 天时，造成的最大伤害。

• 当 (a_i leq m) 时，有转移：

[f_i = f_{i + 1} + a_i ]

• 当 (a_i > m) 时，有转移：

[f_i = max(f_{i + 1}, f_{i + d + 1} + a_i) ]

最后答案即为 (f_1)​，时间复杂度 (mathcal{O}(n))​。

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int maxn=2e5+10;
ll dp[maxn];
ll a[maxn];
int main()
{
    int n,d,m;
    scanf("%d%d%d",&n,&d,&m);
    for(int i=1;i<=n;i++) scanf("%lld",&a[i]);
    
    for(int i=n;i>=1;i--)
    {
        if(a[i]>m)
        {
            dp[i]=max(dp[i+d+1]+a[i],dp[i+1]);
        }
        else
        {
            dp[i]=dp[i+1]+a[i];
        } 
    }
    cout<<dp[1]<<endl;
}

Problem C. bomb

先考虑一下没有环的情况。

结论：答案为有向图的最长链长度。

证明：

• 由于最长链上的点两两不能在同一轮染色。所以最优解 (ge) 最长链长度。

• 令 (dist_x)​​ 表示从 (x)​​ 开始的最长链长度，可以在第 (i)​​ 轮中对所有 (dist_x = i)​​ 的点 (x)​​ 进行染色，恰好可以得到一个轮数为最长链长度的方案。所以最优解 (leq)​​ 最长链长度。

故答案为有向图的最长链长度。

再考虑一下有环的情况。

考虑使用 tarjan 将这张有向图中的所有强连通分量求出，那么对于任意一个强连通分量，如果在一轮染色中对这个强连通分量中的其中一个点进行染色，那么这个强连通分量中的其他点、以及它所能到达的强连通分量里的所有点都不能染色。

我们发现这与题目的定义类似，对于一个大小为 (s) 的强连通分量，可以将其缩为一个需要被染色 (s) 次的点。这样就得到了一张新图，答案即为这张新图中的最长链，使用拓扑排序即可解决。

时间复杂度 (mathcal{O}(n + m))​。

#include <iostream>
#include <algorithm>
#include <cstring>
#include <cstdio>
#include <map>
#include <queue>
using namespace std;
typedef long long ll;
const int N = 2e6 + 10;
int e[N], h[N], ne[N], idx;
int cnt;
int stk[N];
int top;
int in_stk[N];
int dfn[N];
int low[N];
int scc_cnt;
int in_scc_num[N];
int scc_size[N];
int n, m;
void add(int a, int b) {
	e[idx] = b;
	ne[idx] = h[a];
	h[a] = idx ++;
}
void tarjan(int u) {
	dfn[u] = low[u] = ++ cnt;
	stk[++ top] = u;
	in_stk[u] = 1;
	
	for(int i = h[u]; i != -1; i = ne[i]) {
		int v = e[i];
		
		if(!dfn[v]) {
			tarjan(v);
			low[u] = min(low[u], low[v]);
			
		} else if(in_stk[v]) {
			low[u] = min(low[u], dfn[v]);
		}
	}

	if(dfn[u] == low[u]) {
		scc_cnt ++;
		int x;
		do {
			x = stk[top --];
			in_stk[x] = 0;
			in_scc_num[x] = scc_cnt;
			scc_size[scc_cnt] ++;
		} while(x != u);
	}

}



int h1[N], e1[N], ne1[N], idx1;
int du[N];
void add1(int a, int b) {
	
	e1[idx1] = b;
	ne1[idx1] = h1[a];
	h1[a] = idx1 ++;
	du[b] ++;
	
}


ll dist[N];


void topsort() {
	queue<int>q;
	
	
	for(int i = 1; i <= scc_cnt; i ++)
		if(du[i] == 0)
			q.push(i), dist[i] = scc_size[i];

	while(q.size()) {
		int u = q.front();
		q.pop();
		for(int i = h1[u]; i != -1; i = ne1[i]) {
			int v = e1[i];
			dist[v] = max(dist[v], dist[u] + scc_size[v]);
			if(--du[v] == 0)
				q.push(v);
		}
	}
}
void solve() {
	cin >> n >> m;
	idx = idx1 = 0;
	memset(h, -1, sizeof h);
	memset(h1, -1, sizeof h1);
	for(int i = 1; i <= m; i ++) {
		int x, y;
		cin >> x >> y;
		add(x, y);
	}
	for(int i = 1; i <= n ; i++) {
		if(!dfn[i])
			tarjan(i);
	}

	for(int u = 1; u <= n; u ++) {
		for(int i = h[u]; i != -1; i = ne[i]) {
			int v = e[i];
			if(in_scc_num[u] == in_scc_num[v])
				continue;

			add1(in_scc_num[u], in_scc_num[v]);

		}
	}

	topsort();
	ll ans = 0;
	for(int i = 1; i <= scc_cnt; i ++)
		ans = max(ans, dist[i]);

	cout << ans << "
";
	
	return ;
}
signed main() {
	int t = 1;
	while(t--)
		solve();
	return 0;
}

Problem D. maxsum

题目要求的就是前 (w) 大子段和。

因为 (a_i) 为非负整数，所以如果当前最大的是 ([l,r]) 子段，那么易得 ([l,r+1]) 子段和 ([l-1,r]) 子段一定之前就已经取出。

最大的子段一定为 ([1,n])，所以首先将 ([1,n]) 加入堆。设堆顶为 ([l,r])，则 ([l+1,r]) 和 ([l,r-1]) 可以成为备选答案，加入堆并用 map 判重即可。

时间复杂度 (O(w log ⁡n))。

#include <iostream>
#include <algorithm>
#include <cstring>
#include <cstdio>
#include <map>
#include <queue>
using namespace std;
typedef long long ll;
const int N = 1e5 + 10;

typedef pair<int, int> pii;
int n, w;
ll a[N];
ll sum;
map<pii, bool>mp;
struct node {
	int l, r;
	ll sum;
	friend bool operator < (node x, node y) {
		return x.sum < y.sum;
	}

};
priority_queue<node>q;
void solve() {
	cin >> n >> w;
	for(int i = 1; i <= n; i++)
		cin >> a[i], sum += a[i];
	q.push({1, n, sum});
	mp[ {1, n}] = 1;
	while(w) {
		auto now = q.top();
		q.pop();
		w--;
		cout << now.sum << " ";
		if(now.l != now.r) {
			if(!mp[ {now.l + 1, now.r}])
				q.push({now.l + 1, now.r, now.sum - a[now.l]});
			if(!mp[ {now.l, now.r - 1}])
				q.push({now.l, now.r - 1, now.sum - a[now.r]});
			mp[ {now.l + 1, now.r}] = mp[ {now.l, now.r - 1}] = 1;
		}
	}
	return ;
}
signed main() {
	int t = 1;
	while(t--)
		solve();
	return 0;
}

Problem E. array

令 (A, B) 分别为两个数组的最大值，则 (c_i ≥ max(A, B))。
另一方面，如果 (c_i > max(A, B))，则对应的 (a_i) 和 (b_j) 都非零。
根据题意，两个数组最多只有 (5000) 个数非零，暴力计算即可。

时间复杂度 (mathcal O(n + 5000^2))。

#include <iostream>
#include <algorithm>
#include <cstring>
#include <cstdio>
#include <map>
#include <queue>
using namespace std;
typedef long long ll;
const int N = 2e5 + 10;
const int mod = 1e9 + 7;
int a[N];
int b[N];
int c[N];
int n;
int nexta[N];
int nextb[N];
int maxv = 0;
void solve() {
	cin >> n;
	for(int i = 0; i <= n - 1; i ++)
		cin >> a[i], maxv = max(maxv, a[i]);

	for(int i = 0; i <= n - 1; i ++)
		cin >> b[i], maxv = max(maxv, b[i]);

	for(int i = 0; i < n; i ++)
		c[i] = maxv;
	nexta[n - 1] = nextb[n - 1] = n;
	for(int i = n - 1; i >= 1; i --) {
		if(a[i])
			nexta[i - 1] = i;
		else
			nexta[i - 1] = nexta[i];
	}

	for(int i = n - 1; i >= 1; i --) {
		if(b[i])
			nextb[i - 1] = i;
		else
			nextb[i - 1] = nextb[i];
	}


	for(int i = 0 ; i < n; i = nexta[i])
		for(int j = 0; j < n; j = nextb[j])
			c[(i + j) % n] = max(c[(i + j) % n], a[i] + b[j]);

	for(int i = 0; i < n ; i++)
		cout << c[i] << " ";
	cout << "
";
	return ;

}
signed main() {
	int t = 1;
	while(t--)
		solve();
	return 0;
}

Problem F. fare

设 (u)​​ 是树的一个结点，设 (num1(u),num2(u),dp1(u))​​ 分别表示以 (u)​​ 为根的子树内所有参加活动的选手到 (u)​​ 的距离的零次方和，一次方和，二次方和，再设 (tmp1(u),tmp2(u),dp2(u))​​ 分别表示除去以 (u)​​ 为根的子树内所有参加活动的选手到 (u)​​ 的距离的零次方和，一次方和，二次方和。

遍历整棵树就可以求出所有结点的 (num1,num2,dp1)​​​，再次遍历整棵树，我们使用换根 (dp) 的方法，就可以求出所有结点的 (tmp1,tmp2,dp2)​​​，最后的答案便为

(minlimits_{1 leq u leq n} ⁡{ dp1(u)+dp2(u) })

时间复杂度 (O(n))​。

#include <iostream>
#include <algorithm>
#include <cstring>
#include <cstdio>
#include <map>
#include <queue>
using namespace std;
typedef long long ll;
const int N = 4e5 + 10;
const int mod = 1e9 + 7;
int h[N], e[N], ne[N], idx, w[N];
ll c[N];
ll num1[N];// num
ll num2[N];// edge*num
ll dp1[N];// 1 as root

ll tmp1[N];
ll tmp2[N];
ll dp2[N];//change root

ll n;

void add(int a, int b, int c) {
	e[idx] = b;
	w[idx] = c;
	ne[idx] = h[a];
	h[a] = idx ++;
}
void dfs1(int u, int fa) {
	num1[u] = c[u];
	
	for(int i = h[u]; i != -1; i = ne[i]) {
		int v = e[i];
		if(v == fa)
			continue;

		dfs1(v, u);

		num1[u] += num1[v];
		num2[u] += num2[v] + w[i] * num1[v];
		dp1[u] += dp1[v] + 2 * w[i] * num2[v] + w[i] * w[i] * num1[v];
	}
	
	


}

void dfs2(int u, int fa) {
	
	for(int i = h[u]; i != -1; i = ne[i]) {
		int v = e[i];
		if(v == fa)
			continue;

		tmp1[v] = num1[u] - num1[v] + tmp1[u];
                   
		tmp2[v] = num2[u] - num1[v] * w[i] - num2[v] + tmp2[u] + tmp1[v] * w[i];

		dp2[v] = dp1[u] - dp1[v] - 2 * w[i] * num2[v] - w[i] * w[i] * num1[v] + dp2[u] +
				 2 * w[i] * (num2[u] - num1[v] * w[i] - num2[v] + tmp2[u] ) 
				 + w[i] * w[i] * tmp1[v];
		dfs2(v, u);


	}
}
void solve() {
	cin >> n;
	for(int i = 1; i <= n; i ++)
		cin >> c[i];
	idx = 0;
	memset(h, -1 , sizeof h);

	for(int i = 1; i < n; i ++) {
		int a, b, c;
		cin >> a >> b >> c;
		add(a, b, c);
		add(b, a, c);
	}

	dfs1(1, 0);
	
//	cout << dp1[1] << " ";
	dfs2(1, 0);

	ll ans = 1e18;
	for(int i = 1; i <= n; i++)
		ans = min(ans , dp1[i] + dp2[i]);

	cout << ans << "
";
	return ;
}
signed main() {
	int t = 1;
	while(t--)
		solve();
	return 0;
}

Problem G. road

注意到异或不进位，每位之间是相互独立的，于是我们可以考虑计算一下每一位的贡献。

比如说我们现在要计算一下第 ( ext{index}) 位的贡献。

设原图中的邻接矩阵为 (A)，将第 ( ext{index}) 位提取出来可以得到一个 (0 / 1) 矩阵 (A')。我们让点在新的邻接矩阵上走，注意到按位异或的性质，一条非简单路径只有经过了奇数个边，才会对答案有 (2^{ ext{index}}) 的贡献。

考虑 dp。

设 (f(i, x, 0 / 1)) 表示：有多少条 (S) 到 (x) 的边数为 (i) 的路径，且经过了 偶数 / 奇数 个 (1)。

则有状态转移方程：
( f(i, v, 0) = sumlimits_{(u, v) in E, a'_{u, v} = 0} f(i - 1, u, 0) + sumlimits_{(u, v) in E, a'_{u, v} = 1} f(i - 1, u, 1) \ f(i, v, 1) = sumlimits_{(u, v) in E, a'_{u, v} = 0} f(i - 1, u, 1) + sumlimits_{(u, v) in E, a'_{u, v} = 1} f(i - 1, u, 0) \ )
则第 ( ext{index}) 位对答案的贡献即为 (f_{k, T, 1} ast 2^{ ext{index}})。

直接做的话显然会 TLE，但是这个转移方程对于所有的阶段 (i) 转移都一样，于是矩阵快速幂加速即可。

单次回答询问的时间复杂度是 (mathcal{O}(n^3 log^2 ext{size})) 的，其中 ( ext{size}) 表示值域大小。
还不够优秀。

因为回答的询问只关注单源信息，注意到 (1 imes n) 的矩阵和 (n imes n) 的矩阵相乘的复杂度是 (mathcal{O}(n^2)) 的，(n imes n) 的矩阵和 (n imes n) 的矩阵相乘的复杂度是 (mathcal{O}(n^3)) 的。

于是我们可以先把 (A^1, A^2, A^4, A^8, ...) 都预处理出来。
然后将 (k) 二进制分解，若 (k) 的第 (i) 位上为 (1)，则令 (f = f ast A^{2^i})。
这样做的话，单次询问的时间复杂度就成功降到了 (mathcal{O}(n^2 log^2 n))。

总时间复杂度 (mathcal{O}(n^3 log^2 ext{size} + Q ast n^2 log^2 ext{size}))​。

#include <cstdio>
#include <cstring>
#include <algorithm>

using namespace std;

inline int read() {
	int x = 0, f = 1; char s = getchar();
	while (s < '0' || s > '9') { if (s == '-') f = -f; s = getchar(); }
	while (s >= '0' && s <= '9') { x = x * 10 + s - '0'; s = getchar(); }
	return x * f;
}

const int N = 45;
const int mod = 1e9 + 7;

int n, m, Q;
int a[N][N];
int f[31][31][N][N][2];

void mul(int f[N][N][2], int a[N][N][2]) {
	static int c[N][N][2]; memset(c, 0, sizeof(c));

	for (int i = 1; i <= n; i ++)
		for (int j = 1; j <= n; j ++)
			for (int k = 1; k <= n; k ++) {
				c[i][j][0] = (c[i][j][0] + 1ll * a[i][k][0] * a[k][j][0]) % mod;
				c[i][j][0] = (c[i][j][0] + 1ll * a[i][k][1] * a[k][j][1]) % mod;
				c[i][j][1] = (c[i][j][1] + 1ll * a[i][k][0] * a[k][j][1]) % mod;
				c[i][j][1] = (c[i][j][1] + 1ll * a[i][k][1] * a[k][j][0]) % mod;
			}

	memcpy(f, c, sizeof(c));
}

void mulstar(int a[N][2], int b[N][N][2]) {
	static int c[N][2]; memset(c, 0, sizeof(c));

	for (int j = 1; j <= n; j ++)
		for (int k = 1; k <= n; k ++) {
			c[j][0] = (c[j][0] + 1ll * a[k][0] * b[k][j][0]) % mod;
			c[j][0] = (c[j][0] + 1ll * a[k][1] * b[k][j][1]) % mod;
			c[j][1] = (c[j][1] + 1ll * a[k][0] * b[k][j][1]) % mod;
			c[j][1] = (c[j][1] + 1ll * a[k][1] * b[k][j][0]) % mod;
		}

	memcpy(a, c, sizeof(c));
}

void calc(int index) {
	for (int u = 1; u <= n; u ++)
		for (int v = 1; v <= n; v ++) {
			if (a[u][v] == 0x3f3f3f3f) continue;
			int w = a[u][v] >> index & 1;
			f[index][0][u][v][w] = 1;
		}

	for (int i = 1; i <= 30; i ++)
		mul(f[index][i], f[index][i - 1]);
}

void work() {
	int S = read(), T = read(), k = read() - 1;
	int ans = 0;

	for (int index = 0; index <= 30; index ++) {
		int d[N][2]; memset(d, 0, sizeof(d));
		for (int i = 1; i <= n; i ++) {
			if (a[S][i] == 0x3f3f3f3f) continue;
			int w = a[S][i] >> index & 1;
			d[i][w] = 1; 
		}

		for (int i = 0; i <= 30; i ++)
			if (k >> i & 1) mulstar(d, f[index][i]);

		ans = (ans + 1ll * (1 << index) * d[T][1]) % mod;
	}

	printf("%d
", ans);
}

int main() {
	n = read(), m = read();

	memset(a, 0x3f, sizeof(a));
	for (int i = 1; i <= m; i ++) {
		int u = read(), v = read(), w = read();
		a[u][v] = a[v][u] = w;
	}

	for (int index = 0; index <= 30; index ++)
		calc(index);

	Q = read();

	while (Q --)    work();

	return 0;
}

Problem H. xor

将 (a_i) 和 (a_j) 按二进制每一位展开：

(a_i =sumlimits_{x=0}^{29}2^xa_{i, x})，(a_j =sumlimits_{y=0}^{29}2^ya_{j, y})

(a_i ⊕ a_j=sumlimits_{k=0}^{29} 2^k[a_{i, k} ot=a_{j, k}])

将原式展开，

(sumlimits_{i=1}^{n}sumlimits_{j=1}^{n} (a_i ⊕a_j)^2)

(=sumlimits_{i=1}^{n}sumlimits_{j=1}^{n}(sumlimits_{x=0}^{29} 2^x[a_{i, x} ot=a_{j, x}])(sumlimits_{y=0}^{29} 2^y[a_{i, y} ot=a_{j, y}]))

(=sumlimits_{i=1}^{n}sumlimits_{j=1}^{n}sumlimits_{x=0}^{29}sumlimits_{y=0}^{29} 2^{x+y}[a_{i, x} ot=a_{j, x}][a_{i, y} ot=a_{j, y}])

(=sumlimits_{x=0}^{29}sumlimits_{y=0}^{29} 2^{x+y}sumlimits_{i=1}^{n}sumlimits_{j=1}^{n}[a_{i, x} ot=a_{j, x}][a_{i, y} ot=a_{j, y}])

枚举 (x, y, a_{i,x}, a_{j,x}, a_{i,y}, a_{j,y}) 的取值，乘以对应的方案数（即有多少个 (i) 满足第 (x) 位为 (a_{i,x}) 且第 (y) 位为 (a_{i,y})）即可。
预处理出 (f_{i,j},x,y) 表示有多少个数满足第 (i) 位为 (x) 且第 (j) 位为 (y)，则方案数可以查表 (O(1)) 得到。
时间复杂度 (O(n log^2 a))。

#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;
const int N = 50010;
const int mod = 1e9 + 7;
typedef long long ll;
int n;
ll a[N];
ll inv[N];
ll f[40][40][2][2];
signed main() {
	
	/*
5
12133
544123
47778
69885
11234
	*/
	
	cin >> n;

	for(int i = 1; i <= n; i ++)
		cin >> a[i];


//	ll res = 0 ;
//	for(int i = 1; i <= n; i ++)
//		for(int j = 1; j <= n; j ++)
//			res = (res + (a[i] ^ a[j]) * (a[i] ^ a[j]) % mod) % mod;
//	cout << res << endl;
	
	inv[0] = 1;
	for(int i = 1; i <= 63; i ++)
		inv[i] = 1ll * 2 * inv[i - 1] % mod;


	for(int i = 1; i <= n; i ++)
		for(int j = 0; j <= 31; j ++)
			for(int k = 0; k <= 31; k ++)
				f[j][k][a[i] >> j & 1][a[i] >> k & 1]  ++;


	ll ans = 0;
	for(int i = 0; i <= 31; i ++)
		for(int j = 0; j <= 31; j ++)
			for(int k = 0; k <= 1; k ++)
				for(int z = 0; z <= 1; z ++)
					ans = (ans + 1ll * f[i][j][k][z] * f[i][j][k ^ 1][z ^ 1] % mod * inv[i + j] % mod) % mod;

	cout << ans << '
';
	
	return 0;
}

Problem I. Fibonacci sequence

存在递推式

(sqrt{3 + f_n f_{n - 1}} = sqrt{3 + f_{n - 1} f_{n - 2}} + f_{n - 1})

证明如下

((sqrt{3 + f_{n - 1} f_{n - 2}} + f_{n - 1})^2)

(=3 + f_{n - 1} f_{n - 2} + 2 f_{n - 1} sqrt{3 + f_{n - 1} f_{n - 2}} + f_{n - 1}^2)

(=3 + (f_{n - 1} + f_{n - 2} + 2 sqrt{3 + f_{n - 1}f_{n - 2}}) f_{n - 1})

(=3 + f_n f_{n - 1})

于是可以构建矩阵来优化递推，但是考虑到有 (T(T leq 10^6)) 组，我们可以考虑光速幂。

具体来说，我们想要知道 (G^t)（(G) 为转移矩阵），我们预处理出矩阵 (G) 的 (G^0, G^1, G^2 cdots, G^s, G^{2s},G^{3s}, G^{⌈frac{MAX}{s}⌉s})（其中，(s=sqrt{MAX})）这样，我们就可以 (O(3^3))（即矩阵乘法的复杂度） 求出一个 (G^t (t = ks + b))，进而优化复杂度。

#include <bits/stdc++.h>
#define int long long
const int N = 3 + 10;
const int RN = 1e6 + 10;

using namespace std;

struct matrix
{
	int l, r;
	long long val[N][N];
} a, G, t[100000], t2[100000];

int x[RN];

int G0, G1, M, n, ans = 1;

inline matrix mul(matrix A, matrix B)
{
	matrix C; C.l = A.l, C.r = B.r;
	memset(C.val, 0, sizeof(C.val));

	for (int i = 1; i <= C.l; ++ i)
		for (int j = 1; j <= C.r; ++ j)
			for (int k = 1; k <= A.r; ++ k)
				C.val[i][j] = (C.val[i][j] + A.val[i][k] * B.val[k][j]) % M;

	return C;
}

inline matrix pow(matrix A, int b)
{
	matrix res = A; b --;
	for (; b; b >>= 1)
	{
		if (b & 1) res = mul(res, A);
		A = mul(A, A);
	}
	return res;
}

inline matrix solve(int b)
{
	a.l = 3, a.r = 3;
	a.val[1][1] = 1, a.val[1][2] = 1, a.val[1][3] = 1;
	a.val[2][1] = 1, a.val[2][2] = 0, a.val[3][3] = 0;
	a.val[3][1] = 2, a.val[3][2] = 0, a.val[3][3] = 1;
	
	return pow(a, b);
}

int Max;

signed main()
{
	
	cin >> G0 >> G1 >> M >> n;


	G.l = 1, G.r = 3;
	G.val[1][1] = G1, G.val[1][2] = G0, G.val[1][3] = (int)sqrt(3 + G1 * G0) % M;
		
	a.l = 3, a.r = 3;
	a.val[1][1] = 1, a.val[1][2] = 1, a.val[1][3] = 1;
	a.val[2][1] = 1, a.val[2][2] = 0, a.val[3][3] = 0;
	a.val[3][1] = 2, a.val[3][2] = 0, a.val[3][3] = 1;

	for (int i = 1; i <= n; ++ i)
	{
		cin >> x[i];
		Max = max(Max, x[i]);	
	}
	Max = sqrt(Max) + 1;
	t[1] = solve(1);
	for (int i = 2; i <= Max; ++ i)
		t[i] = mul(t[i - 1], t[1]);
	t2[1] = solve(Max);
	for (int i = 2; i <= Max; ++ i)
		t2[i] = mul(t2[i - 1], t2[1]);
	for (int i = 1; i <= n; ++ i)
	{
		if (x[i] == 0) 
			ans = ans * G0 % M;
		else if (x[i] == 1) 
			ans = ans * G1 % M;
		else
		{
			x[i] -= 1;
			int s = x[i] / Max;
			matrix c1 = t2[s], c2 = t[x[i] - Max * s];
			matrix res = mul(c1, c2);
			if (s == 0) res = c2;
			if (x[i] - Max * s == 0) res = c1;
			ans = ans * mul(G, res).val[1][1] % M;
		}
	}
	cout << ans << endl;
	
	return 0;
} 

Problem J. mesh

直接上莫比乌斯反演：
( egin{aligned} sumlimits_{i = 1}^n sumlimits_{j = 1}^m varphi(gcd(i, j)) & = sumlimits_{d = 1} varphi(d) cdot sumlimits_{i = 1}^n sumlimits_{j = 1}^m [gcd(i, j) = d]\ & = sumlimits_{d = 1} varphi(d) cdot sumlimits_{i = 1}^{leftlfloorfrac{n}{d} ight floor} sumlimits_{j = 1}^{leftlfloorfrac{m}{d} ight floor} [gcd(i, j) = 1]\ & = sumlimits_{d = 1} varphi(d) cdot sumlimits_{i = 1}^{leftlfloorfrac{n}{d} ight floor} sumlimits_{j = 1}^{leftlfloorfrac{m}{d} ight floor} sumlimits_{p mid gcd(i, j)} mu(p) \ & = sumlimits_{d = 1} varphi(d) cdot sumlimits_{p = 1} mu(p) cdot leftlfloorfrac{n}{pd} ight floor cdot leftlfloorfrac{m}{pd} ight floor end{aligned} )

记 (T = pd)，则式子化简为：
( sumlimits_{T = 1}^{min(n, m)} leftlfloorfrac{n}{T} ight floor cdot leftlfloorfrac{m}{T} ight floor cdot left(sumlimits_{d mid T} varphi(d) cdot muleft(frac{T}{d} ight) ight) )
设 (g = varphi ast mu)（这里的 (ast) 是狄利克雷卷积，下文同），则式子化简为：
( sumlimits_{T = 1}^{min(n, m)} leftlfloorfrac{n}{T} ight floor cdot leftlfloorfrac{m}{T} ight floor cdot g(T) )
注意到函数 (g) 显然为积性函数，可以使用线性筛将 (g) 的前缀和筛出。

配合数论分块即可做到 (mathcal{O}(n)) 预处理、(mathcal{O}(sqrt{n})) 回答询问。但还不够优秀。

---

可以考虑杜教筛，设：
( S(n) = sumlimits_{i = 1}^n g(i) )
记 (d = 1 ast 1)，其中函数 (1(n) = 1)，显然函数 (d) 为约数个数函数。

那么可以得到 (S(n)) 有关 (Sleft(leftlfloorfrac{n}{d} ight floor ight)) 的递推式：
( d(1)S(n) = sumlimits_{i = 1}^n (g ast d)(i) - sumlimits_{i = 2}^n d(i) cdot Sleft(leftlfloorfrac{n}{i} ight floor ight) \ )
化简：
( d(1)S(n) = sumlimits_{i = 1}^n (varphi ast 1 ast mu ast 1)(i) - sumlimits_{i = 2}^n d(i) cdot Sleft(leftlfloorfrac{n}{i} ight floor ight) \ d(1)S(n) = frac{n cdot (n + 1)}{2} - sumlimits_{i = 2}^n d(i) cdot Sleft(leftlfloorfrac{n}{i} ight floor ight) )
该式显然需要数论分块来做。取一个阈值 (K)，分段处理：

• 对于 (n leq K) 的部分，使用线性筛求出每个 (S(n), sum_{i = 1}^n d(i))。
• 对于 (n > K) 的部分，使用数论分块求解 (sum_{i = 1}^n d(i))，使用上述递推式求解 (S(i))。

可以使用积分近似证明当 (K = mathcal{O}(n^{frac{2}{3}})) 时，取到杜教筛的理论最优复杂度 (mathcal{O}(n^{frac{2}{3}}))。

#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <map>
#include <ctime>

using namespace std;

const int N = 1000100, SIZE = 1e6;
const int mod = 1e9 + 7;

int n, m;

int k, prime[N], low[N];
int h[N], d[N];
int Sd[N], Sh[N];

void GetPrimes(int N) {
	h[1] = 1, d[1] = 1;

	for (int i = 2; i <= N; i ++) {
		if (!low[i]) {
			prime[++ k] = i;
			low[i] = i;
			d[i] = 2;
			h[i] = i - 2;

			if (1ll * i * i <= N) {
				low[i * i] = i * i;
				d[i * i] = 3;
				h[i * i] = (i * i - i) - (i - 1);

				long long v = 1ll * i * i * i;
				int c = 3;

				while (v <= N) {
					low[v] = v;
					d[v] = c + 1;
					h[v] = h[v / i] * i;

					v *= i, c ++;
				}
			}
		}

		for (int j = 1; j <= k; j ++) {
			if (prime[j] > N / i) break;

			if (i % prime[j] == 0)
				low[i * prime[j]] = low[i] * prime[j];
			else
				low[i * prime[j]] = prime[j];

			d[i * prime[j]] = d[(i * prime[j]) / low[i * prime[j]]] * d[low[i * prime[j]]];
			h[i * prime[j]] = h[(i * prime[j]) / low[i * prime[j]]] * h[low[i * prime[j]]];

			if (i % prime[j] == 0) break;
		}
	}

	for (int i = 1; i <= N; i ++)
		Sd[i] = (Sd[i - 1] + d[i]) % mod,
		Sh[i] = (Sh[i - 1] + h[i]) % mod;
}

map<int, int> Gd, Gh;

int sig(int n) {
	if (n <= SIZE) return Sd[n];
	if (Gd[n]) return Gd[n];
	int ans = 0;
	for (int x = 1, Nx; x <= n; x = Nx + 1) {
		Nx = n / (n / x);
		ans = (ans + 1ll * (Nx - x + 1) * (n / x)) % mod; 
	}
	return Gd[n] = ans;
}

int S(int n) {
	if (n <= SIZE) return Sh[n];
	if (Gh[n]) return Gh[n];
	int ans = (1ll * n * (n + 1) / 2) % mod;
	int sub = 0; 
	for (int x = 2, Nx, cur; x <= n; x = Nx + 1) {
		Nx = n / (n / x);

		cur = ((sig(Nx) - sig(x - 1)) % mod + mod) % mod;
		cur = 1ll * cur * S(n / x) % mod;

		sub = (sub + cur) % mod;
	}
	ans = ((ans - sub) % mod + mod) % mod;
	return Gh[n] = ans;
}

int main() {
	cin >> n >> m;

	GetPrimes(SIZE);

	int ans = 0;
	for (int x = 1, Nx, cur; x <= min(n, m); x = Nx + 1) {
		Nx = min(n / (n / x), m / (m / x));

		cur = ((S(Nx) - S(x - 1)) % mod + mod) % mod;
		cur = 1ll * cur * (n / x) % mod;
		cur = 1ll * cur * (m / x) % mod;
 
		ans = (ans + cur) % mod;
	}

	cout << ans << endl;

	return 0;
}