做题计划

唐。

Luogu P3455 [POI2007] ZAP-Queries

莫比乌斯反演。

令：\(a\le b\)

求：

\[\sum\limits_{i=1}^a\sum\limits_{j=1}^b[\gcd(i,j)=x] \]

消掉 \(x\)：

\[=\sum\limits_{i=1}^{\left\lfloor\frac{a}{x}\right\rfloor}\sum\limits_{j=1}^{\left\lfloor\frac{b}{x}\right\rfloor}[\gcd(i,j)=1] \]

根据莫比乌斯的定义：

\[=\sum\limits_{i=1}^{\left\lfloor\frac{a}{x}\right\rfloor}\sum\limits_{j=1}^{\left\lfloor\frac{b}{x}\right\rfloor}\sum_{d|\gcd(i,j)} \mu(d) \]

改成和的形式：

\[=\sum\limits_{i=1}^{\left\lfloor\frac{a}{x}\right\rfloor}\sum\limits_{j=1}^{\left\lfloor\frac{b}{x}\right\rfloor}\sum_{d=1}^{\left\lfloor\frac{a}{x}\right\rfloor} \mu(d) \times [d \mid \gcd(i,j)] \]

移项：

\[=\sum\limits_{d=1}^{\left\lfloor\frac{a}{x}\right\rfloor} \mu(d)\sum\limits_{i=1}^{\left\lfloor\frac{b}{x}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{a}{x}\right\rfloor} \times [d \mid \gcd(i,j)] \]

显然 \(d|i,d|j\)，又因为满足条件的 \(i\) 在 \(1 \sim \left\lfloor\frac{a}{x}\right\rfloor\) 中出现 \(\left\lfloor\frac{a}{x \times d}\right\rfloor\) 次，\(j\) 在 \(1 \sim \left\lfloor\frac{b}{x}\right\rfloor\) 中出现 \(\left\lfloor\frac{b}{x \times d}\right\rfloor\) 次，所以：

\[=\sum\limits_{d=1}^{\left\lfloor\frac{a}{x}\right\rfloor} \mu(d)\left\lfloor\frac{a}{x \times d}\right\rfloor\left\lfloor\frac{b}{x \times d}\right\rfloor \]

然后可以预处理出 \(\mu\) 的前缀和，另外的整除分块即可。时间复杂度 \(\mathcal{O(m \log \sqrt{n})}\)。

点击查看代码

#include <bits/stdc++.h>
#define int long long
	
using namespace std;

const int N = 5e4;

int pri[N + 100], cnt, mu[N + 100], a, b, x, n, sumu[N + 100];
bool vis[N + 100];

void solve() {
	cin >> a >> b >> x;
	if (a > b) {
		swap(a, b);
	}
	a /= x;
	b /= x;
	int ans = 0;
	for (int l = 1, r; l <= a; l = r + 1) {
		r = min(a / (a / l), b / (b / l));
		ans += (a / l) * (b / l) * (sumu[r] - sumu[l - 1]);
	}
	cout << ans << endl;
}

void init() {
	mu[1] = 1;
	vis[1] = 1;
	for (int i = 2; i <= N; i++) {
		if (!vis[i]) {
			pri[++cnt] = i;
			mu[i] = -1;
		}
		for (int j = 1; j <= cnt && i * pri[j] <= N; j++) {
			vis[i * pri[j]] = 1;
			if (i % pri[j] == 0) {
				mu[i * pri[j]] = 0;
				break;
			} else {
				mu[i * pri[j]] = mu[i] * mu[pri[j]];
			}
		}
	}
	for (int i = 1; i <= N; i++) {
		sumu[i] = sumu[i - 1] + mu[i];
	}
}

signed main() {
	init();
	cin >> n;
	while (n--) {
		solve();
	}
	return 0;
}

Luogu P2261 [CQOI2007] 余数求和

求：

\[\sum\limits_{i=1}^nk\bmod i \]

根据 \(\bmod\) 的性质：

\[a \bmod b = a - b \times \left\lfloor\frac{a}{b}\right\rfloor \]

得：

\[\sum\limits_{i=1}^n k - i \times \left\lfloor\frac{k}{i}\right\rfloor \]

整除分块即可，注意一下边界。

如果 \(k \le l\)，那么直接跳出循环，否则 \(r = \min(n, \left\lfloor\frac{\left\lfloor\frac{k}{l}\right\rfloor}{l}\right\rfloor)\)。

点击查看代码

#include <bits/stdc++.h>
#define int long long

using namespace std;

int n, k, ans;

signed main() {
	cin >> n >> k;
	for (int l = 1, r; l <= n; l = r + 1) {
		if (k <= l) {
			break;
		} else {
			r = min(n, k / (k / l));
		}
		ans += 1ll * (r - l + 1) * (k / l) * (l + r) / 2ll;
	}
	cout << n * k - ans;
	return 0;
}

Luogu P2257 YY 的 GCD

求：

\[\sum\limits_{i=1}^n\sum\limits_{j=1}^m[\gcd(i,j)\in \text{prime}] \]

枚举 \(k\)

\[=\sum\limits_{k = 1, k\in \text{prime}}^n\sum\limits_{i=1}^n\sum\limits_{j=1}^m[\gcd(i,j)=k] \]

除以 \(k\)：

\[=\sum\limits_{k = 1,k \in \text{prime}}^n\sum\limits_{i=1}^{\left\lfloor\frac{n}{k}\right\rfloor}\sum\limits_{j=1}^{\left\lfloor\frac{m}{k}\right\rfloor}[\gcd(i,j)=1] \]

莫比乌斯定义转换：

\[=\sum\limits_{k = 1,k \in \text{prime}}^n\sum\limits_{i=1}^{\left\lfloor\frac{n}{k}\right\rfloor}\sum\limits_{j=1}^{\left\lfloor\frac{m}{k}\right\rfloor}\sum\limits_{p|\gcd(i,j)} \mu(p) \]

套路：

\[=\sum\limits_{k =1,k \in \text{prime}}^n\sum\limits_{p=1}^{\left\lfloor\frac{n}{k}\right\rfloor} \mu(p) \times \left\lfloor\frac{n}{kp}\right\rfloor \times \left\lfloor\frac{m}{kp}\right\rfloor \]

设 \(T=kp\)

原式：

\[=\sum\limits_{k =1,k \in \text{prime}}^n\sum\limits_{p=1}^{\left\lfloor\frac{n}{k}\right\rfloor} \mu(p) \times \left\lfloor\frac{n}{T}\right\rfloor \times \left\lfloor\frac{m}{T}\right\rfloor \]

枚举 \(T\)：

\[=\sum\limits_{T = 1}^{n} \left\lfloor\frac{n}{T}\right\rfloor \times \left\lfloor\frac{m}{T}\right\rfloor \sum\limits_{k \mid T,k \in \text{prime}} \mu(\frac{T}{k}) \]

线性筛预处理 + 整除分块即可。

点击查看代码

#include <bits/stdc++.h>

using namespace std;

const int N = 1e7 + 100;

int n, m, mu[N], cnt, pri[N], f[N], t;
long long sumu[N];
bool vis[N];

void init() {
	mu[1] = 1;
	vis[1] = 1;
	for (int i = 2; i <= N - 100; ++i) {
		if (!vis[i]) {
			pri[++cnt] = i;
			mu[i] = -1;
		}
		for (int j = 1; j <= cnt && i * pri[j] <= N - 100; ++j) {
			vis[i * pri[j]] = 1;
			if (i % pri[j] == 0) {
				mu[i * pri[j]] = 0;
				break;
			} else {
				mu[i * pri[j]] = mu[i] * mu[pri[j]];
			}
		}
	}
	for (int i = 1; i <= cnt; ++i) {
		for (int j = 1; j * pri[i] <= N - 100; ++j) {
			f[j * pri[i]] += mu[j];
		}
	}
	for (int i = 1; i <= N - 100; ++i) {
		sumu[i] = sumu[i - 1] + f[i];
//		cout << sumu[i] << ' '; 
	}
}

void solve() {
	cin >> n >> m;
	long long sum = 0;
	for (int l = 1, r; l <= min(n, m); l = r + 1) {
		r = min(n / (n / l), m / (m / l));
		sum += (sumu[r] - sumu[l - 1]) * (n / l) * (m / l);
	}
	cout << sum << '\n';
}

int main() {
	ios::sync_with_stdio(0);
	cin.tie(0), cout.tie(0);
	init();
	cin >> t;
	while (t--) {
		solve();
	}
	return 0;
}

SPOJ SP26017

求：

\[\sum\limits_{i=1}^n\sum\limits_{j=1}^m\gcd(i,j) \]

令：\(n \le m\)

\(\texttt{F1}：\)

枚举 \(d\)，除以 \(d\)：

\[=\sum\limits_{d=1}^{n}d\sum\limits_{j=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum\limits_{j=1}^{\left\lfloor\frac{m}{d}\right\rfloor}[\gcd(i,j)=1] \]

莫比乌斯反演：

\[=\sum\limits_{d=1}^{n}d\sum\limits_{j=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum\limits_{j=1}^{\left\lfloor\frac{m}{d}\right\rfloor}\sum\limits_{p \mid \gcd(i,j)} \mu(p) \]

枚举 \(p\)：

\[=\sum\limits_{d=1}^{n}d\sum\limits_{p=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum\limits_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}[p \mid i]\sum\limits_{j = 1}^{\left\lfloor\frac{m}{d}\right\rfloor}[p \mid j] \mu(p) \]

套路：

\[=\sum\limits_{d=1}^{n}d\sum\limits_{p=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\left\lfloor\frac{n}{pd}\right\rfloor\left\lfloor\frac{m}{pd}\right\rfloor \mu(p) \]

设 \(T=pd\)：

原式：

\[=\sum\limits_{d=1}^{n}d\sum\limits_{p=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\left\lfloor\frac{n}{T}\right\rfloor\left\lfloor\frac{m}{T}\right\rfloor \mu(p) \]

枚举 \(T\)：

\[=\sum\limits_{T=1}^{n}\sum\limits_{p|T}\left\lfloor\frac{n}{T}\right\rfloor\left\lfloor\frac{m}{T}\right\rfloor \mu(p) \times \frac{T}{p} \]

调转顺序：

\[=\sum\limits_{T=1}^{n}\left\lfloor\frac{n}{T}\right\rfloor\left\lfloor\frac{m}{T}\right\rfloor\sum\limits_{p|T} \mu(p) \times \frac{T}{p} \]

\(\texttt{F2}\)：

有：

\[\sum\limits_{d|n}\varphi(d)=n \]

带入原式：

\[=\sum\limits_{i=1}^n\sum\limits_{j=1}^m\sum\limits_{d|\gcd(i,j)}\varphi(d) \]

套路：

\[=\sum\limits_{d=1}^n\sum\limits_{i=1}^n\sum\limits_{j=1}^m\varphi(d) \times [d \mid i] \times [d|j] \]

接着套路：

\[=\sum\limits_{d=1}^n\varphi(d) \times \left\lfloor\frac{n}{d}\right\rfloor \times \left\lfloor\frac{m}{d}\right\rfloor \]

给出 \(\texttt{F2}\) 代码:

点击查看代码

#include <bits/stdc++.h>
#define int long long

using namespace std;

const int N = 5e4;
const int mod = 1e9 + 7;

int n, m, n1, m1, T, pri[N + 100], phi[N + 100], cnt, sum[N + 100], p1, p2;
bool vis[N + 100];

int num(int n, int m) {
	if (n > m) {
		swap(n, m);
	}
	int ans = 0;
	for (int l = 1, r; l <= n; l = r + 1) {
		r = min(n / (n / l), m / (m / l));
//		cout << "l = " << l << ' ' << "r = " << r << '\n';
//		cout << n << ' ' << m << ' ' << min(n / (n / l), m / (m / l)) << endl;
		ans = (ans + (phi[r] - phi[l - 1]) * (n / l) % mod * (m / l) % mod) % mod;
	}
	return ans;
}

void solve() {
	cin >> n >> m >> n1 >> m1;
//	cout << num(n1, m1) << ' ' << num(n - 1, m1) << ' ' << num(n1, m - 1) << ' ' << num(n - 1, m - 1) << endl;
	cout << ((num(n1, m1) - num(n - 1, m1) - num(n1, m - 1) + num(n - 1, m - 1)) % mod + mod) % mod << endl;
}

void init() {
	phi[1] = 1;
	vis[1] = 1;
	for (int i = 2; i <= N; i++) {
		if (!vis[i]) {
			pri[++cnt] = i;
			phi[i] = i - 1;
		}
		for (int j = 1; j <= cnt && i * pri[j] <= N; j++) {
			vis[i * pri[j]] = 1;
			if (i % pri[j] == 0) {
				phi[i * pri[j]] = phi[i] * pri[j];
				break;
			} else {
				phi[i * pri[j]] = phi[i] * (pri[j] - 1);
			}
		}
	}
	for (int i = 1; i <= N; i++) {
		phi[i] += phi[i - 1];
		phi[i] %= mod;
//		sum[i] = sum[i - 1] + phi[i];
	}
//	cout << sum[100] << endl;
}

signed main() {
	init();
	cin >> T;
	cin >> p1 >> p2;
	while (T--) {
		solve();
	}
	return 0;  
}

Luogu P6055 [RC-02] GCD

求：

\[\sum\limits_{i=1}^n\sum\limits_{j=1}^n\sum\limits_{p=1}^{\left\lfloor\frac{n}{j}\right\rfloor}\sum\limits_{q=1}^{\left\lfloor\frac{n}{j}\right\rfloor}[\gcd(i,j)=1][\gcd(p,q)=1] \]

左右同乘 \(j\)：

\[\sum\limits_{i=1}^n\sum\limits_{j=1}^n\sum\limits_{p=1}^n\sum\limits_{q=1}^n[\gcd(i,j)=1][\gcd(p,q)=j] \]

所以 \(j=\gcd(p,q)\)，带入 \(\gcd(i,j)=1\)，得：

\[=\sum\limits_{i=1}^n\sum\limits_{p=1}^n\sum\limits_{q=1}^n[\gcd(i,\gcd(p,q))=1] \]

根据 \(\gcd\) 的性质（\(\gcd(a,\gcd(b,c))=\gcd(a,b,c)\)）：

\[=\sum\limits_{i=1}^n\sum\limits_{p=1}^n\sum\limits_{q=1}^n[\gcd(i,p,q)=1] \]

莫比乌斯反演：

\[=\sum\limits_{i=1}^n\sum\limits_{p=1}^n\sum\limits_{q=1}^n\sum\limits_{d|\gcd(i,p,q)}\mu(d) \]

套路：

\[=\sum\limits_{d=1}^n{\left\lfloor\frac{n}{d}\right\rfloor}^3\mu(d) \]

杜教筛维护 \(\mu\) 函数，整除分块即可。

点击查看代码

#include <bits/stdc++.h>
#define int long long

using namespace std;

const int N = 5e5;
const int mod = 998244353;

int cnt, n, num, sum[N + 100], pri[N + 100];
bool vis[N + 100];
unordered_map<int, int> mu;

void init() {
	mu[1] = 1;
	for (int i = 2; i <= N; i++) {
		if (!vis[i]) {
			pri[++cnt] = i;
			mu[i] = -1;
		}
		for (int j = 1; j <= cnt && i * pri[j] <= N; j++) {
			vis[i * pri[j]] = 1;
			if (i % pri[j]) {
				mu[i * pri[j]] = -mu[i];  
			} else {
				mu[i * pri[j]] = 0;
				break;
			}
		}
	}
	for (int i = 1; i <= N; i++) {
		sum[i] = sum[i - 1] + mu[i];
		sum[i] %= mod;
	}
}

int suma(int n) {
	if (n <= N) {
		return sum[n];
	}
	if (mu[n]) {
		return mu[n];
	}
	int ans = 1ll;
	for (int l = 2, r; l <= n; l = r + 1) {
		r = n / (n / l);
		ans -= (r - l + 1) * suma(n / l);
	}
	return mu[n] = ans;
}

signed main() {
	ios::sync_with_stdio(0);
	cin.tie(0);
	cout.tie(0);
	cin >> n;
	init();
	for (int l = 1, r; l <= n; l = r + 1) {
		r = n / (n / l);
		int w = (n / l) * (n / l) % mod * (n / l) % mod;
	//		cout << l << ' ' << r << endl;
		num = (num + (suma(r) - suma(l - 1) + mod) % mod * w % mod) % mod;
	}	
	cout << num;
	return 0;
}

Luogu P4449 于神之怒加强版

求：

\[\sum\limits_{i=1}^n\sum\limits_{j=1}^m\gcd(i,j)^k \]

枚举 \(d\)：

\[=\sum\limits_{d=1}^nd^k\sum\limits_{i=1}^n\sum\limits_{j=1}^m[\gcd(i,j)=k] \]

除以 \(d\)：

\[=\sum\limits_{d=1}^nd^k\sum\limits_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum\limits_{j=1}^{\left\lfloor\frac{m}{d}\right\rfloor}[\gcd(i,j)=1] \]

莫比乌斯反演：

\[=\sum\limits_{d=1}^nd^k\sum\limits_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum\limits_{j=1}^{\left\lfloor\frac{m}{d}\right\rfloor}\sum\limits_{p=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\mu(p)[d|i][d|j] \]

套路：

\[=\sum\limits_{d=1}^nd^k\sum\limits_{p=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\mu(p)\left\lfloor\frac{n}{pd}\right\rfloor\left\lfloor\frac{m}{pd}\right\rfloor \]

设 \(T=pd\)：

\[=\sum\limits_{d=1}^nd^k\sum\limits_{p=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\mu(p)\left\lfloor\frac{n}{T}\right\rfloor\left\lfloor\frac{m}{T}\right\rfloor \]

枚举 \(T\)：

\[=\sum\limits_{T=1}^n\left\lfloor\frac{n}{T}\right\rfloor\left\lfloor\frac{m}{T}\right\rfloor\sum\limits_{d|T}\mu(\frac{T}{d})d^k \]

线性筛即可。

在线性筛中：

令 \(g(x)=\sum\limits_{d|x}\mu(\frac{x}{d})d^k\)：

显然 \(g(x)=id_{k} \ * \ \mu\)（\(*\) 为范德蒙德卷积），为积性函数（若 \(f(x),g(x)\) 均为积性函数，那么 \(h(x)=\sum\limits_{d|x}f(d)g(\frac{x}{d})\) 也为积性函数）。

有：

\(g(x)=\begin{cases}-(T^{x-1})^k&d=T^{x-1}\\(T^x)^k&d=T^x\\0& \texttt{other}\end{cases}\)

那么 \(g(p^{k+1})=p^{(k+1)x}-p^{kx}=p^{x}(p^k-p^{(k-1)x})=p^xg(p^x)\)

有两种情况：

\(f_{i\times pri_j}=\begin{cases}f_i \times f_{pri_j} &\gcd(i,pri_j)=1&线性函数的性质\\f_i \times pri_j^k&pri_j |i&pri_j 的指数增加 1，即 g(pri_j^{k+1})\end{cases}\)

点击查看代码

#include <bits/stdc++.h>
#define int long long

using namespace std;

const int N = 5e6;
const int mod = 1e9 + 7;

int n, m, k, pri[N + 100], phi[N + 100], cnt, T;
bool vis[N + 100];

int qpow(int a, int b) {
	if (!b) {
		return 1;
	}
	int tmp = qpow(a, b / 2);
	if (b & 1) {
		return tmp % mod * tmp % mod * a % mod;
	}
	return tmp % mod * tmp % mod;
}

void solve() {
cin >> n >> m;
	if (n > m) {
		swap(n, m);
	}
	int ans = 0;
	for (int l = 1, r; l <= n; l = r + 1) {
		r = min(n / (n / l), m / (m / l));
		ans = (ans + (phi[r] - phi[l - 1] + mod) % mod * (n / l) % mod * (m / l) % mod) % mod;
	}
	cout << ans << endl;
	return ;
}

void init() {
	phi[1] = 1;
	for (int i = 2; i <= N; i++) {
		if (!vis[i]) {
			pri[++cnt] = i;
//			cout << i << ' ' << k << ' ' << qpow(i, k) << endl;
			phi[i] = (qpow(i, k) - 1 + mod) % mod;
		}
		for (int j = 1; j <= cnt && i * pri[j] <= N; j++) {
			vis[i * pri[j]] = 1;
			if (i % pri[j] == 0) {
				phi[i * pri[j]] = phi[i] * qpow(pri[j], k) % mod;
				break;
			} else {
				phi[i * pri[j]] = phi[i] * phi[pri[j]] % mod;
			}	
		}
	}
	for (int i = 1; i <= N; i++) {
		phi[i] += phi[i - 1];
		phi[i] %= mod;
		//		sum[i] = sum[i - 1] + phi[i];
	}
	//	cout << phi[10] << endl;
}

signed main() {
	ios::sync_with_stdio(0);
	cin.tie(0), cout.tie(0);
	cin >> T >> k;
	init();
	while (T--) {
		solve();
	}
	return 0;  
}

Luogu P3172 [CQOI2015] 选数

求：

\[\sum\limits_{a_1=L}^R\sum\limits_{a_2=L}^R \cdots \sum\limits_{a_n=L}^R[\gcd(a_1,a_2 \cdots a_n)=k] \]

除以 \(k\)，令 \(l=\left\lfloor\frac{R}{k}\right\rfloor,r=\left\lfloor\frac{L}{k}\right\rfloor\)：

\[\sum\limits_{a_1=l}^r\sum\limits_{a_2=l}^r \cdots \sum\limits_{a_n=l}^r[\gcd(a_1,a_2 \cdots a_n)=1] \]

枚举 \(p\)：

\[\sum\limits_{a_1=l}^r\sum\limits_{a_2=l}^r \cdots \sum\limits_{a_n=l}^r\sum\limits_{p|\gcd(a_1,a_2 \cdots a_n)}\mu(p) \]

套路：

\[\sum\limits_{p=1}^r\mu(p){(\left\lfloor\frac{r}{p}\right\rfloor-\left\lfloor\frac{l-1}{p}\right\rfloor)}^n \]

杜教筛预处理，整除分块即可。

点击查看代码

#include <bits/stdc++.h>
#define int long long

using namespace std;

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

int cnt, n, num, sum[N + 100], pri[N + 100], n1, k1, l1, h1;
bool vis[N + 100];
unordered_map<int, int> mu;

void init() {
mu[1] = 1;
for (int i = 2; i <= N; i++) {
	if (!vis[i]) {
		pri[++cnt] = i;
		mu[i] = -1;
	}
	for (int j = 1; j <= cnt && i * pri[j] <= N; j++) {
		vis[i * pri[j]] = 1;
		if (i % pri[j]) {
			mu[i * pri[j]] = -mu[i];  
		} else {
			mu[i * pri[j]] = 0;
			break;
		}
	}
}
for (int i = 1; i <= N; i++) {
//		cout << i << ' ';
	sum[i] = sum[i - 1] + mu[i];
	sum[i] %= mod;
}
}

int qpow(int a, int b) {
	if (!b) {
		return 1;
	}
	int tmp = qpow(a, b / 2);
	if (b & 1) {
		return tmp % mod * tmp % mod * a % mod;
	}
	return tmp % mod * tmp % mod;
}

int suma(int n) {
	if (n <= N) {
		return sum[n];
	}
	if (mu[n]) {
		return mu[n];
	}
	int ans = 1ll;
	for (int l = 2, r; l <= n; l = r + 1) {
		r = n / (n / l);
		ans -= (r - l + 1) * suma(n / l);
	}
	return mu[n] = ans;
}

signed main() {
	ios::sync_with_stdio(0);
	cin.tie(0);
	cout.tie(0);
	cin >> n1 >> k1 >> l1 >> h1;
	l1--;
	l1 /= k1;
	h1 /= k1;
	//	cout << l1 << ' ' << h1 << endl; 
	// 2 2 2 4
	init();
	for (int l = 1, r; l <= h1; l = r + 1) {
		if (l1 / l) {
			r = min(h1 / (h1 / l), l1 / (l1 / l));
		} else {
			r = h1 / (h1 / l); 
		}
	//		cout << l << ' ' << r << endl;
		num = (num + (suma(r) - suma(l - 1) + mod) % mod * qpow((h1 / l - l1 / l + mod) % mod, n1) % mod) % mod;
	}
	cout << num;
	return 0;
}

Luogu P6825 「EZEC-4」求和

毒瘤卡常题。

求：

\[\sum\limits_{i=1}^n\sum\limits_{j=1}^n\gcd(i,j)^{i+j} \]

枚举 \(d\)：

\[\sum\limits_{d=1}^{n}\sum\limits_{i=1}^n\sum\limits_{j=1}^nd^{i+j}[\gcd(i,j)=d] \]

除以 \(d\)：

\[\sum\limits_{d=1}^{n}\sum\limits_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum\limits_{j=1}^{\left\lfloor\frac{n}{d}\right\rfloor}d^{d(i+j)}[\gcd(i,j)=1] \]

反演：

\[\sum\limits_{d=1}^{n}\sum\limits_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum\limits_{j=1}^{\left\lfloor\frac{n}{d}\right\rfloor}d^{d(i+j)}\sum\limits_{p|\gcd(i,j)}\mu(p) \]

套路：

\[\sum\limits_{d=1}^n\sum\limits_{p=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\mu(p)\sum\limits_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum\limits_{j=1}^{\left\lfloor\frac{n}{d}\right\rfloor}d^{d(i+j)}[p|i][p|j] \]

不考虑 \([p|i]\) 和 \([p|j]\) 了，直接就是 \(\left\lfloor\frac{n}{d}\right\rfloor\) 里面有 \(\left\lfloor\frac{n}{pd}\right\rfloor\) 个 \(p\) 的倍数，但是 \(i,j\) 同除多除了一个 \(p\)，在 \(d^{d(i+j)}\) 中补上一个 \(p\) 得 \(d^{pd(i+j)}\)：

\[=\sum\limits_{d=1}^n\sum\limits_{p=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\mu(p)\sum\limits_{i=1}^{\left\lfloor\frac{n}{pd}\right\rfloor}\sum\limits_{j=1}^{\left\lfloor\frac{n}{pd}\right\rfloor}d^{pd(i+j)} \]

设 \(T=pd\)：

\[=\sum\limits_{d=1}^n\sum\limits_{p=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\mu(p)\sum\limits_{i=1}^{\left\lfloor\frac{n}{T}\right\rfloor}\sum\limits_{j=1}^{\left\lfloor\frac{n}{T}\right\rfloor}d^{T(i+j)} \]

拆开：

\[=\sum\limits_{d=1}^n\sum\limits_{p=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\mu(p)\sum\limits_{i=1}^{\left\lfloor\frac{n}{T}\right\rfloor}(d^{T})^i\sum\limits_{j=1}^{\left\lfloor\frac{n}{T}\right\rfloor}(d^{T})^j \]

发现 \(i,j\) 两个循环一样：

\[=\sum\limits_{d=1}^n\sum\limits_{p=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\mu(p)\Large(\small\sum\limits_{i=1}^{\left\lfloor\frac{n}{T}\right\rfloor}(d^{T})^i\Large)^{\small2} \]

发现括号里面为等比数列，有：

\(s(x)=\begin{cases}s(\frac{n}{2})(1+q^{\frac{n}{2}})&2|x\\s(\frac{n-1}{2})(1+q^{\frac{n+1}{2}})&2\nmid x\end{cases}\)

\(\mathcal{O(\log n)}\) 求和即可。

点击查看代码

#include <bits/stdc++.h>

using namespace std;

const int N = 1.5e6;

int T, n, p, pri[N + 100], mu[N + 100], cnt;
bool vis[N + 100];

namespace FastIO
{
	
	class FastIOBase
	{
	protected:
#ifdef OPENIOBUF
		static const int BUFSIZE=1<<16;
		char buf[BUFSIZE+1];
		int buf_p=0;
#endif
		FILE *target;
	public:
#ifdef OPENIOBUF
		virtual void flush()=0;
#endif
		FastIOBase(FILE *f): target(f){}
		~FastIOBase()=default;
	};
	
	class FastOutput final: public FastIOBase
	{
#ifdef OPENIOBUF
	public:
		inline void flush()
		{ fwrite(buf,1,buf_p,target),buf_p=0; }
#endif
	private:
		inline void __putc(char x)
		{
#ifdef OPENIOBUF
			if(buf[buf_p++]=x,buf_p==BUFSIZE)flush();
#else
			putc(x,target);
#endif
		}
		template<typename T>
		inline void __write(T x)
		{
			static char stk[64],*top;top=stk;
			if(x<0) return __putc('-'),__write(-x);
			do *(top++)=x%10,x/=10; while(x);
			for(;top!=stk;__putc(*(--top)+'0'));
		}
	public:
		FastOutput(FILE *f=stdout): FastIOBase(f){}
#ifdef OPENIOBUF
		~FastOutput(){ flush(); }
#endif
		inline FastOutput &operator <<(char x)
		{ return __putc(x),*this; }
		inline FastOutput &operator <<(const char *s)
		{ for(;*s;__putc(*(s++)));return *this; }
		inline FastOutput &operator <<(const string &s)
		{ return (*this)<<s.c_str(); }
		template<typename T>
		inline enable_if_t<is_integral<T>::value,FastOutput&> operator <<(const T &x)
		{ return __write(x),*this; }
		template<typename ...T>
		inline void writesp(const T &...x)
		{ initializer_list<int>{(this->operator<<(x),__putc(' '),0)...}; }
		template<typename ...T>
		inline void writeln(const T &...x)
		{ initializer_list<int>{(this->operator<<(x),__putc('\n'),0)...}; }
		template<typename Iter>
		inline void writesp(Iter begin,Iter end)
		{ while(begin!=end) (*this)<<*(begin++)<<' '; }
		template<typename Iter>
		inline void writeln(Iter begin,Iter end)
		{ while(begin!=end) (*this)<<*(begin++)<<'\n'; }
		
	}qout;
	
	class FastInput final: public FastIOBase
	{
#ifdef OPENIOBUF
	public:
		inline void flush()
		{ buf[fread(buf,1,BUFSIZE,target)]=EOF,buf_p=0; }
#endif
	private:
		bool __eof,__ok;
		inline char __getc()
		{
			if(__eof) return EOF;
#ifdef OPENIOBUF
			if(buf_p==BUFSIZE) flush();
			return ~buf[buf_p]?buf[buf_p++]:(__eof=true,EOF);
#else
			char ch=getc(target);
			return ~ch?ch:(__eof=true,EOF);
#endif
		}
	public:
		FastInput(FILE *f=stdin): FastIOBase(f),__eof(false),__ok(true)
#ifdef OPENIOBUF
		{ buf_p=BUFSIZE; }
#else
		{}
#endif
		inline bool eof()const { return __eof; }
		inline char peek() { return __getc(); }
		explicit inline operator bool()const { return __ok; }
		inline FastInput &operator >>(char &x)
		{ while(isspace(x=__getc()));return *this; }
		template<typename T>
		inline enable_if_t<is_integral<T>::value,FastInput&> operator >>(T &x)
		{
			char ch,sym=0;
			while(isspace(ch=__getc()));
			if(ch=='-') sym=1,ch=__getc();
			if(__eof) return __ok=false,*this;
			for(x=0;isdigit(ch);x=(x<<1)+(x<<3)+(ch^48),ch=__getc());
			return sym?x=-x:x,*this;
		}
		inline FastInput &operator >>(char *s)
		{
			while(isspace(*s=__getc()));
			if(__eof) return __ok=false,*this;
			for(;!isspace(*s) && !__eof;*(++s)=__getc());
			return *s='\0',*this;
		}
		inline FastInput &operator >>(string &s)
		{
			char str_buf[(1<<8)+1],*p=str_buf;
			char *const buf_end=str_buf+(1<<8);
			while(isspace(*p=__getc()));
			if(__eof) return __ok=false,*this;
			for(s.clear(),p++;;p=str_buf)
			{
				for(;p!=buf_end && !isspace(*p=__getc()) && !__eof;p++);
				*p='\0',s.append(str_buf);
				if(p!=buf_end) break;
			}
			return *this;
		}
		template<typename ...T>
		inline void read(T &...x)
		{ initializer_list<int>{(this->operator>>(x),0)...}; }
		template<typename Iter>
		inline void read(Iter begin,Iter end)
		{ while(begin!=end) (*this)>>*(begin++); }
	}qin;
	
} // namespace FastIO
using FastIO::qin,FastIO::qout;

inline long long qpow(int a, int b, int p) {
	if (!b) {
		return 1;
	}
	long long tmp = qpow(a, b / 2, p);
	if (b & 1) {
		return tmp % p * tmp % p * a % p;
	}
	return tmp % p * tmp % p;
}

inline void init() {
	mu[1] = 1;
	for (int i = 2; i <= N; ++i) {
		if (!vis[i]) {
			pri[++cnt] = i;
			mu[i] = -1;
		}
		for (int j = 1; j <= cnt && i * pri[j] <= N; ++j) {
			vis[i * pri[j]] = 1;
			if (i % pri[j] == 0) {
				mu[i * pri[j]] = 0;
				break;
			}
			mu[i * pri[j]] = -mu[i];
		}
	}
}

inline long sum(int a, int n, int mod) {
	if (n == 1) {
		return a;
	}
	if (!(n & 1)) {
		return 1ll * (1ll + qpow(a, n / 2, mod)) * sum(a, n / 2, mod) % mod;
	}
	return (1ll * (1ll + qpow(a, n / 2, mod)) * sum(a, n / 2, mod) % mod + qpow(a, n, mod)) % mod;
}

int main() {
//	freopen("1.txt", "w", stdout);
	qin >> T;
	init();
	while (T--) {		
		qin >> n >> p;
		long long suma = 0;
		for (int i = 1; i <= n; ++i) {
			long long ans = qpow(i, i, p);
			long long t = ans;
			for (int j = 1; j <= n / i; ++j) {
//				suma = (suma + mu[j] * sum(n / (i * j), t, p) % p * sum(n / (i * j), t, p) % p) % p;
//				t = (t * ans) % p; 
//				qout << sum(t, n / (i * j), p) << endl;
				suma = 1ll * (suma + (1ll * sum(t, n / (i * j), p) % p * sum(t, n / (i * j), p) % p) * (mu[j] + p) % p) % p;
				t = t * ans % p;
			}
		}
		qout << suma % p << '\n';
	}
	return 0;
}

P1829 [国家集训队] Crash的数字表格 / JZPTAB

求：

\[\sum\limits_{i=1}^{n}\sum\limits_{j=1}^m\operatorname{lcm}(i,j) \]

令：\(n\le m\)：

有：

\[\sum\limits_{i=1}^m\sum\limits_{j=1}^m\frac{ij}{\gcd(i,j)} \]

枚举 \(d\)：

\[\sum\limits_{d=1}^n\sum\limits_{i=1}^m\sum\limits_{j=1}^m[\gcd(i,j)=d]\frac{ij}{d} \]

除以 \(d\)：

\[\sum\limits_{d=1}^nd\sum\limits_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum\limits_{j=1}^{\left\lfloor\frac{n}{d}\right\rfloor}[\gcd(i,j)=1]ij \]

莫比乌斯反演：

\[\sum\limits_{d=1}^nd\sum\limits_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum\limits_{j=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum\limits_{p|\gcd(i,j)}\mu(p)ij \]

移相：

\[\sum\limits_{d=1}^nd\sum\limits_{p=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\mu(p)\sum\limits_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}i[p\mid i]\sum\limits_{j=1}^{\left\lfloor\frac{m}{d}\right\rfloor}j[p\mid j] \]

枚举 \(p\) 的整数数倍：

\[\sum\limits_{d=1}^nd\sum\limits_{p=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\mu(p)\sum\limits_{up=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum\limits_{vp=1}^{\left\lfloor\frac{m}{d}\right\rfloor}uvp^2 \]

提出 \(p^2\)：

\[\sum\limits_{d=1}^nd\sum\limits_{p=1}^{\left\lfloor\frac{n}{d}\right\rfloor}p^2\mu(p)(\sum\limits_{u=1}^{\left\lfloor\frac{n}{pd}\right\rfloor}u)(\sum\limits_{v=1}^{\left\lfloor\frac{m}{pd}\right\rfloor}v) \]

化一下：

\[\sum\limits_{d=1}^nd\sum\limits_{p=1}^{\left\lfloor\frac{n}{d}\right\rfloor}p^2\mu(p)S(\left\lfloor\frac{n}{pd}\right\rfloor)S(\left\lfloor\frac{m}{pd}\right\rfloor) \]

其中：

\[S(n)=\frac{n(n+1)}{2} \]

线性筛预处理，两次整除分块即可。

点击查看代码

#include <bits/stdc++.h>
#define int long long

using namespace std;

const int N = 1e7;
const int mod = 20101009;

int pri[N + 100], mu[N + 100], cnt, n, m, sum[N + 100], ans;
bool vis[N + 100];

void init() {
	mu[1] = 1;
	for (int i = 2; i <= N; i++) {
		if (!vis[i]) {
			pri[++cnt] = i;
			mu[i] = -1;
		}
		for (int j = 1; j <= cnt && i * pri[j] <= N; j++) {
			vis[i * pri[j]] = 1;
			if (i % pri[j] == 0) {
				mu[i * pri[j]] = 0;
				break;
			}
			mu[i * pri[j]] = -mu[i];
		}
	}
	for (int i = 1; i <= N; i++) {
		sum[i] = (sum[i - 1] + i * i % mod * mu[i] % mod) % mod;
	}
}

int sol(int n) {
	return (1ll * n * (n + 1) / 2ll) % mod;
} 

int suma(int n, int m) {
	if (n > m) {
		swap(n, m);
	}
	int ans = 0;
	for (int l = 1, r; l <= n; l = r + 1) {
		r = min(n / (n / l), m / (m / l));
		ans = (ans + (sum[r] - sum[l - 1] + mod) % mod * sol(n / l) % mod * sol(m / l) % mod) % mod; 	
	}
	return ans % mod;
}

void solve(int n, int m) {
	if (n > m) {
		swap(n, m);
	}
	int ans = 0;
	for (int l = 1, r; l <= n; l = r + 1) {
		r = min(n / (n / l), m / (m / l));
		ans = (ans + 1ll * (r - l + 1) * (l + r) / 2ll % mod * suma(n / l, m / l) % mod) % mod;
	}
	cout << ans << endl;
}

signed main() {
	init();
	cin >> n >> m;
	solve(n, m);
	return 0;
}

Luogu P3327 [SDOI2015] 约数个数和

求：

\[\sum\limits_{i=1}^n\sum\limits_{j=1}^md(ij) \]

其中：

\[d(x)=\sum\limits_{p|x}1 \]

有：

\[d(i\times j)=\sum\limits_{x|i}\sum\limits_{y|j}[\gcd(x,y)=1] \]

带入：

\[\sum\limits_{i=1}^n\sum\limits_{j=1}^m\sum\limits_{x|i}\sum\limits_{y|j}[\gcd(x,y)=1] \]

套路，消掉 \(i,j\)：

\[\sum\limits_{x=1}^n\sum\limits_{y=1}^m\left\lfloor\frac{n}{x}\right\rfloor\left\lfloor\frac{m}{y}\right\rfloor[\gcd(x,y)=1] \]

反演：

\[\sum\limits_{x=1}^n\sum\limits_{y=1}^m\left\lfloor\frac{n}{x}\right\rfloor\left\lfloor\frac{m}{y}\right\rfloor\sum\limits_{p|\gcd(x,y)}\mu(p) \]

套路：

\[\sum\limits_{x=1}^n\left\lfloor\frac{n}{x}\right\rfloor\sum\limits_{y=1}^m\left\lfloor\frac{m}{y}\right\rfloor\sum\limits_{p=1}^n[p\mid x][p\mid y]\mu(p) \]

改 \(x\) 为 \(ip\)，\(y\) 为 \(jp\)，满足 \(x\mid p\)，\(y\mid p\)：

\[\sum\limits_{p=1}^n\mu(p)\sum\limits_{i=1}^{\left\lfloor\frac{n}{p}\right\rfloor}\sum\limits_{j=1}^{\left\lfloor\frac{m}{p}\right\rfloor}\left\lfloor\frac{n}{ip}\right\rfloor\left\lfloor\frac{m}{jp}\right\rfloor \]

移项：

\[\sum\limits_{p=1}^n\mu(p)\sum\limits_{i=1}^{\left\lfloor\frac{n}{p}\right\rfloor}\left\lfloor\frac{n}{ip}\right\rfloor\sum\limits_{j=1}^{\left\lfloor\frac{m}{p}\right\rfloor}\left\lfloor\frac{m}{jp}\right\rfloor \]

设：

\[f(x)=\sum\limits_{i=1}^x\left\lfloor\frac{x}{i}\right\rfloor \]

有：\(\dfrac{n}{ip}=\dfrac{\left\lfloor\dfrac{n}{p}\right\rfloor}{i}=f(\left\lfloor\dfrac{n}{ip}\right\rfloor)\)，\(\dfrac{m}{jp}\) 同理：

\[\sum\limits_{p=1}^n\mu(p)f(\left\lfloor\frac{n}{p}\right\rfloor)f(\left\lfloor\frac{m}{p}\right\rfloor) \]

预处理 \(\mu\) 函数和 \(f\) 数组，查询的时候整除分块即可。

时间复杂度 \(\mathcal{O(N \log N)}\)，\(\mathcal{N=5\times 10^4}\)。

点击查看代码

#include <bits/stdc++.h>
#define int long long

using namespace std;

const int N = 5e4 + 100;

int pri[N + 100], mu[N + 100], cnt, t, n, m, sum[N + 100], s[N + 100];
bool vis[N + 100];

void init() {
	mu[1] = 1;
	for (int i = 2; i <= N; i++) {
		if (!vis[i]) {
			pri[++cnt] = i;
			mu[i] = -1;
		}
		for (int j = 1; j <= cnt && i * pri[j] <= N; j++) {
			vis[i * pri[j]] = 1;
			if (i % pri[j] == 0) {
				mu[i * pri[j]] = 0;
				break;
			} 
			mu[i * pri[j]] = -mu[i];
		}
 	}
	for (int i = 1; i <= N; i++) {
		sum[i] = sum[i - 1] + mu[i];
	}
	for (int i = 1; i <= N; i++) {
		for (int l = 1, r; l <= i; l = r + 1) {
			r = i / (i / l);
			s[i] += (r - l + 1) * (i / l);
		}
	}
}

void solve() {
	cin >> n >> m;
	if (n > m) {
		swap(n, m);
	}
	int ans = 0;
	for (int l = 1, r; l <= n; l = r + 1) { // sqrt(n)
		r = min(n / (n / l), m / (m / l));
		ans += (sum[r] - sum[l - 1]) * s[n / l] * s[m / l];
	}
	cout << ans << endl;
	return ;
}

signed main() {
	init();
	cin >> t;
	while (t--) {
		solve();
	}
	return 0;
}

Luogu P1891 疯狂 LCM & SPOJ LCMSUM - LCM Sum

求：

\[\sum\limits_{i=1}^n \operatorname{lcm}(i,n) \]

根据 \(\operatorname{lcm}\) 的定义：

\[\sum\limits_{i=1}^n\frac{in}{\gcd(i,n)} \]

把 \(n\) 放在前面：

\[n\sum\limits_{i=1}^n\frac{i}{\gcd(i,n)} \]

枚举 \(d\)：

\[n\sum\limits_{i=1}^n\sum\limits_{d\mid n}[\gcd(i,n)=d]\frac{i}{d} \]

改 \(i\) 为 \(\left\lfloor\frac{i}{d}\right\rfloor\)：

\[n\sum\limits_{d\mid n}\sum\limits_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}[\gcd(id, n)=d]\frac{i}{d} \]

除以 \(d\)：

\[n\sum\limits_{d\mid n}\sum\limits_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}[\gcd(i,\frac{n}{d})=1]i \]

因为 \(d\mid n\)，所以 \(d\) 和 \(\frac{n}{d}\) 对称：

\[n\sum\limits_{d\mid n}\sum\limits_{i=1}^d[\gcd(i,d)=1]i \]

有：

\[n\sum\limits_{d\mid n}\frac{\varphi(d)d}{2} \]

线性筛预处理即可

点击查看代码

#include <bits/stdc++.h>
#define int long long

using namespace std;

const int N = 1e6;

int t, n, pri[N + 100], phi[N + 100], cnt, f[N + 100];
bool vis[N + 100];

void init() {
	phi[1] = 1;
	for (int i = 2; i <= N; i++) {
		if (!vis[i]) {
			pri[++cnt] = i;
			phi[i] = i - 1;
		}
		for (int j = 1; j <= cnt && i * pri[j] <= N; j++) {
			vis[i * pri[j]] = 1;
			if (i % pri[j] == 0) {
				phi[i * pri[j]] = pri[j] * phi[i];
				break;
			}
			phi[i * pri[j]] = (pri[j] - 1) * phi[i];
		}
	}
	for (int i = 1; i <= N; i++) {
		for (int j = 1; i * j <= N; j++) {
			if (i == 1) {
				f[i * j] = 1;
			} else {
				f[i * j] += phi[i] * i / 2ll;
			}
		}
	} 
}

signed main() {
	cin >> t;
	init();
	while (t--) {
		cin >> n;
		cout << f[n] * n << '\n';
	}
	return 0;
}

Luogu P3911 最小公倍数之和

求：

\[\sum\limits_{i=1}^n\sum\limits_{j=1}^n\operatorname{lcm}(a_i,a_j) \]

定义：

\[c_i=\sum\limits_{j=1}^n[a_j=i] \]

转化：

\[\sum\limits_{i=1}^n\sum\limits_{j=1}^nc_i\times c_j\times \operatorname{lcm}(i,j) \]

\(\operatorname{lcm}\) 的定义：

\[\sum\limits_{i=1}^n\sum\limits_{j=1}^nc_i\times c_j\times \frac{ij}{\gcd(i,j)} \]

枚举 \(d\)：

\[\sum\limits_{i=1}^n\sum\limits_{j=1}^n\sum\limits_{d=1}^n[\gcd(i,j)=d]c_i\times c_j\times \frac{ij}{d} \]

除以 \(d\)，反演：

\[\sum\limits_{d=1}^n\sum\limits_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum\limits_{j=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum\limits_{p|\gcd(i,j)}\mu(p)\times c_{id}\times c_{jd}\times i\times j\times d \]

枚举 \(p\)：

\[\sum\limits_{d=1}^n\sum\limits_{p=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\mu(p)p^2\sum\limits_{i=1}^{\left\lfloor\frac{n}{pd}\right\rfloor}\sum\limits_{j=1}^{\left\lfloor\frac{n}{pd}\right\rfloor}c_{ipd}\times c_{jpd}\times i\times j\times d \]

令：

\[T=pd \]

枚举 \(T\)：

\[\sum\limits_{T=1}^nT(\sum\limits_{i=1}^{\left\lfloor\frac{n}{T}\right\rfloor}i\times c_{iT})^2\sum\limits_{p|T}\mu(p)p \]

前面的暴力算，后面预处理即可。

时间复杂度 \(\mathcal{O(N \log N)}\)。

点击查看代码

#include <bits/stdc++.h>
#define int long long

using namespace std;

const int N = 5e4 + 100;

int n, a[N], c[N], pri[N], mu[N], cnt, f[N];
bool vis[N];

void init() {
	mu[1] = 1;
	for (int i = 2; i <= N - 100; i++) {
		if (!vis[i]) {
			pri[++cnt] = i;
			mu[i] = -1;
		}
		for (int j = 1; j <= cnt && i * pri[j] <= N; j++) {
			vis[i * pri[j]] = 1;
			if (i % pri[j] == 0) {
				mu[i * pri[j]] = 0;
				break;
			}
			mu[i * pri[j]] = -mu[i];
		}
	}
	for (int i = 1; i <= N - 100; i++) {
		for (int j = 1; i * j <= N - 100; j++) {
			f[i * j] += i * mu[i];
		}
	}
}

int qpow(int a, int b) {
	if (!b) return 1;
	int tmp = qpow(a, b / 2);
	if (b & 1) return tmp * tmp * a;
	return tmp * tmp;
}

signed main() {
	cin >> n;
	init();
	for (int i = 1; i <= n; i++) {
		cin >> a[i];
		c[a[i]]++;
	}
	int sum = 0;
	for (int T = 1; T <= N; T++) {
		int ans = 0;
		for (int i = 1; i <= N / T; i++) {
			ans += i * c[i * T];
		}	
		sum += T * qpow(ans, 2) * f[T];
	}
	cout << sum;
	return 0;
}

SPOJ PROD1GCD - Product it again

求：

\[\prod\limits_{i=1}^n\prod\limits_{j=1}^m\gcd(i,j) \]

有：

\[\prod\limits_{d=1}^nd^{\sum\limits_{i=1}^n\sum\limits_{j=1}^m[\gcd(i,j)=d]} \]

看指数部分：

\[\sum\limits_{i=1}^n\sum\limits_{j=1}^m[\gcd(i,j)=d] \]

反演：

\[\sum\limits_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum\limits_{j=1}^{\left\lfloor\frac{m}{d}\right\rfloor}\sum\limits_{p|\gcd(i,j)}\mu(p) \]

枚举 \(p\)：

\[\sum\limits_{p=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\left\lfloor\frac{n}{pd}\right\rfloor\left\lfloor\frac{m}{pd}\right\rfloor\mu(p) \]

代回：

\[\prod\limits_{d=1}^nd^{\sum\limits_{p=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\left\lfloor\frac{n}{pd}\right\rfloor\left\lfloor\frac{m}{pd}\right\rfloor\mu(p)} \]

两次整除分块即可。

点击查看代码

#include <bits/stdc++.h>
#define int long long

using namespace std;

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

int t, inv[N + 100], pri[N + 100], mu[N + 100], n, m, cnt;
bool vis[N + 100];

void init() {
  mu[1] = 1;
  inv[0] = inv[1] = 1;
  for (int i = 2; i <= N; i++) {
    if (!vis[i]) {
      pri[++cnt] = i;
      mu[i] = -1;
    }
    for (int j = 1; j <= cnt && i * pri[j] <= N; j++) {
      vis[i * pri[j]] = 1;
      if (i % pri[j] == 0) {
        mu[i * pri[j]] = 0;
        break;
      }
      mu[i * pri[j]] = -mu[i];
    }
    mu[i] = (mu[i] + mu[i - 1]) % mod; 
    inv[i] = (inv[i - 1] * i) % mod;
  }
}

int qpow(int a, int b) {
  if (!b) {
    return 1;
  }
  int tmp = qpow(a, b / 2);
  if (b & 1) {
    return tmp % mod * tmp % mod * a % mod;
  }
  return tmp % mod * tmp % mod; 
}

int solve(int n, int m) {
  if (n > m) {
    swap(n, m);
  } 
  int ans = 0;
  for (int l = 1, r; l <= n; l = r + 1) {
    r = min(n / (n / l), m / (m / l));
    ans = (ans + (mu[r] - mu[l - 1] + (mod - 1)) % (mod - 1) * (n / l) % (mod - 1) * (m / l) % (mod - 1)) % (mod - 1);
  }
  return ans % (mod - 1);
}

signed main() {
  cin >> t;
  init();
  while (t--) {
    cin >> n >> m;
    if (n > m) {
      swap(n, m);
    }
    int ans = 1;
    for (int l = 1, r; l <= n; l = r + 1) {
      r = min(n / (n / l), m / (m / l));
      ans = (ans * qpow(inv[r] * qpow(inv[l - 1], mod - 2) % mod, solve(n / l, m / l) % mod) % mod) % mod;
    }
    cout << ans << endl;
  }
  return 0;
}

Luogu P2260 [清华集训2012] 模积和

求：

\[\sum\limits_{i=1}^n\sum\limits_{j=1}^m(n\bmod i)\times(m\bmod j)(i\neq j) \]

容斥：

\[\sum\limits_{i=1}^n\sum\limits_{j=1}^m(n\bmod i)\times(m\bmod j)-\sum\limits_{i=1}^n(n\bmod i)\times(m\bmod i) \]

\(\bmod\) 的定义：

\[\sum\limits_{i=1}^n\sum\limits_{j=1}^m(n-\left\lfloor\frac{n}{i}\right\rfloor\times i)(m - \left\lfloor\frac{m}{j}\right\rfloor\times j)-\sum\limits_{i=1}^n(n-\left\lfloor\frac{n}{i}\right\rfloor\times i)(m - \left\lfloor\frac{m}{i}\right\rfloor\times i) \]

移相：

\[\sum\limits_{i=1}^n(n-\left\lfloor\frac{n}{i}\right\rfloor\times i)\sum\limits_{j=1}^m(m - \left\lfloor\frac{m}{j}\right\rfloor\times j)-\sum\limits_{i=1}^n(n-\left\lfloor\frac{n}{i}\right\rfloor\times i)(m - \left\lfloor\frac{m}{i}\right\rfloor\times i) \]

展开：

\[(n^2-\sum\limits_{i=1}^ni\times\left\lfloor\frac{n}{i}\right\rfloor)(m^2-\sum\limits_{i=1}^mi\times\left\lfloor\frac{m}{i}\right\rfloor)-\sum\limits_{i=1}^n(nm-mi\times \left\lfloor\frac{m}{i}\right\rfloor-ni\times \left\lfloor\frac{n}{i}\right\rfloor+i^2\times \left\lfloor\frac{m}{i}\right\rfloor\times \left\lfloor\frac{n}{i}\right\rfloor) \]

公式：

\[\sum\limits_{i=1}^ni^2=\frac{n(n+1)(2n+1)}{6} \]

点击查看代码

#include <bits/stdc++.h>
#define int __int128

using namespace std;

const int mod = 19940417;

int n, m, ans1, ans2, ans3, p, ans4, ans5;

int len(int x) {
  return x * (x + 1) / 2;
}

int ipow(int x) {
  return x * (x + 1) * (2 * x + 1) / 6;
}

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

void write(int x) {
  if (x < 0) {
    putchar('-');
    x = -x;
  }
  if (x > 9) {
    write(x / 10);
  }
  putchar(x % 10 + '0');
  return ;
}

signed main() {
	n = read();
	m = read();
	if (n > m) {
	  swap(n, m);
  }
  for (int l = 1, r; l <= n; l = r + 1) {
    r = n / (n / l);
    ans1 = (ans1 + (len(r) - len(l - 1) + mod) % mod * (n / l) % mod) % mod;
  }
  ans1 = n * n - ans1;
  for (int l = 1, r; l <= m; l = r + 1) {
    r = m / (m / l);
    ans2 = (ans2 + (len(r) - len(l - 1) + mod) % mod * (m / l) % mod) % mod;
  }
  ans2 = m * m - ans2;
  p = n * n * m; 
  for (int l = 1, r; l <= n; l = r + 1) {
    r = min(n / (n / l), m / (m / l));
    ans3 = (ans3 + (len(r) - len(l - 1) + mod) % mod * (n / l) % mod) % mod;
  }
  p -= ans3 * m;
  for (int l = 1, r; l <= n; l = r + 1) {
    r = min(n / (n / l), m / (m / l));
    ans4 = (ans4 + (len(r) - len(l - 1) + mod) % mod * (m / l) % mod) % mod;
  }
  p -= ans4 * n;
  for (int l = 1, r; l <= n; l = r + 1) {
    r = min(n / (n / l), m / (m / l));
    ans5 = (ans5 + (ipow(r) - ipow(l - 1) + mod) % mod * (n / l) % mod * (m / l)) % mod;
  }
  p += ans5;
  write((ans1 * ans2 % mod - p % mod + mod) % mod);
  return 0;
}

UVA Count LCM

求：

\[\sum\limits_{i=1}^n\sum\limits_{j=1}^m[\operatorname{lcm}(i,j)=ij] \]

有：

\[\sum\limits_{i=1}^n\sum\limits_{j=1}^m[\frac{ij}{\gcd(i,j)}=ij] \]

即求：

\[\sum\limits_{i=1}^n\sum\limits_{j=1}^m[\gcd(i,j)=1] \]

反演：

\[\sum\limits_{i=1}^n\sum\limits_{j=1}^m\sum\limits_{p|\gcd(i,j)}\mu(p) \]

套路：

\[\sum\limits_{p=1}^n\sum\limits_{i=1}^n[p\mid i]\sum\limits_{j=1}^m[p\mid j]\mu(p) \]

套路：

\[\sum\limits_{p=1}^n\left\lfloor\frac{n}{p}\right\rfloor\left\lfloor\frac{m}{p}\right\rfloor\mu(p) \]

点击查看代码

#include <bits/stdc++.h>
#define int long long

using namespace std;

const int N = 1e6;

int pri[N + 100], mu[N + 100], cnt, t, n, m;
bool vis[N + 100];

void init() {
    mu[1] = 1;
    for (int i = 2; i <= N; i++) {
        if (!vis[i]) {
            pri[++cnt] = i;
            mu[i] = -1;
        }
        for (int j = 1; j <= cnt && i * pri[j] <= N; j++) {
            vis[i * pri[j]] = 1;
            if (i % pri[j] == 0) {
                mu[i * pri[j]] = 0;
                break;
            }
            mu[i * pri[j]] = -mu[i];
        }
        mu[i] += mu[i - 1];
    }
}

void solve() {
    cin >> n >> m;
    if (n > m) {
        swap(n, m);
    }
    int ans = 0;
    for (int l = 1, r; l <= n; l = r + 1) {
        r = min(n / (n / l), m / (m / l));
        ans += (mu[r] - mu[l - 1]) * (n / l) * (m / l);
    }
    cout << ans << '\n';
    return ;
}

signed main() {
    cin >> t;
    init();
    while (t--) {
        solve();
    }
    return 0;
}

SPOJ VLATTICE - Visible Lattice Points

转化题意：

\[\sum\limits_{i=1}^n\sum\limits_{j=1}^n\sum\limits_{k=1}^n[\gcd(i,j,k)=1]+3\sum\limits_{i=1}^n\sum\limits_{j=1}^n[\gcd(i,j)=1]+3 \]

反演：

\[\sum\limits_{i=1}^n\sum\limits_{j=1}^n\sum\limits_{k=1}^n\sum\limits_{d|\gcd(i,j,k)}\mu(d)+3\sum\limits_{i=1}^n\sum\limits_{j=1}^n\sum\limits_{d|\gcd(i,j)}\mu(d)+3 \]

套路：

\[\sum\limits_{d=1}^n\sum\limits_{i=1}^n[d\mid i]\sum\limits_{j=1}^n[d\mid j]\sum\limits_{k=1}^n[d\mid k]\mu(d)+3\sum\limits_{d=1}^n\sum\limits_{i=1}^n[d\mid i]\sum\limits_{j=1}^n[d\mid j]\mu(d)+3 \]

套路：

\[\sum\limits_{d=1}^{n}(\left\lfloor\frac{n}{d}\right\rfloor)^3\mu(d)+3\sum\limits_{d=1}^n(\left\lfloor\frac{n}{d}\right\rfloor)^2\mu(d)+3 \]

合并同类项：

\[\sum\limits_{d=1}^n\mu(d)(\left\lfloor\frac{n}{d}\right\rfloor)^2(\left\lfloor\frac{n}{d}\right\rfloor+3)+3 \]

点击查看代码

#include <bits/stdc++.h>
#define int long long

using namespace std;

const int N = 1e6;

int t, n, mu[N + 100], pri[N + 100], cnt;
bool vis[N + 100];

void init() {
	mu[1] = 1;
	for (int i = 2; i <= N; i++) {
		if (!vis[i]) {
			pri[++cnt] = i;
			mu[i] = -1;
		}
		for (int j = 1; j <= cnt && i * pri[j] <= N; j++) {
			vis[i * pri[j]] = 1;
			if (i % pri[j] == 0) {
				mu[i * pri[j]] = 0;
				break;
			}
			mu[i * pri[j]] = -mu[i];
		}
		mu[i] += mu[i - 1];
	}
}

int qpow(int a, int b) {
	if (!b) {
		return 1;
	}
	int tmp = qpow(a, b / 2);
	if (b & 1) {
		return tmp * tmp * a;
	}
	return tmp * tmp;
}

signed main() {
//	cout << qpow(1, 1)
	cin >> t;
	init();
	while (t--) {
		cin >> n;
		int ans = 0;
		for (int l = 1, r; l <= n; l = r + 1) {
			r = n / (n / l);
			ans += (mu[r] - mu[l - 1]) * qpow(n / l, 2) * (n / l + 3);
		}
		cout << (ans + 3) << endl;
	}
	return 0;
}

SPOJ GCDMAT2 - GCD OF MATRIX (hard)

此题 的升级版。

发现时间复杂度正确，但是常数过大，考虑优化：

• 根号分治

• 把除法变成乘这个数的倒数

稍微卡一卡就可以了。

代码太丑，不放了

Luogu P5221 Product

求：

\[\prod\limits_{i=1}^n\prod\limits_{j=1}^n\frac{\operatorname{lcm}(i,j)}{\gcd(i,j)} \]

有：

\[\prod\limits_{i=1}^n\prod\limits_{j=1}^n\frac{ij}{\gcd(i,j)^2} \]

拆开：

\[\prod\limits_{i=1}^n\prod\limits_{j=1}^nij\times\prod\limits_{i=1}^n\prod\limits_{j=1}^n\gcd(i,j)^{-2} \]

反演：

\[(n!)^{2n}\times(\prod\limits_{d=1}^nd^{\sum\limits_{p=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\left\lfloor\frac{n}{pd}\right\rfloor\left\lfloor\frac{n}{pd}\right\rfloor\mu(p)})^{-2} \]

点击查看代码

#include <bits/stdc++.h>

using namespace std;

const int N = 1e6;
const long long mod = 104857601;

int t, pri[N + 100], mu[N + 100], n, cnt;
long long fac;
bool vis[N + 100];

void init() {
	mu[1] = 1;
	for (int i = 2; i <= N; i++) {
		if (!vis[i]) {
			pri[++cnt] = i;
			mu[i] = -1;
		}
		for (int j = 1; j <= cnt && i * pri[j] <= N; j++) {
			vis[i * pri[j]] = 1;
			if (i % pri[j] == 0) {
				mu[i * pri[j]] = 0;
				break;
			}
			mu[i * pri[j]] = -mu[i];
		}
		mu[i] = (mu[i] + mu[i - 1]) % mod;
	}
	fac = 1ll;
	for (int i = 1; i <= n; i++) {
		fac = (long long)fac * i % mod;
	}
}

long long qpow(int a, int b) {
	if (!b) {
		return 1;
	}
	int tmp = qpow(a, b / 2);
	if (b & 1) {
		return tmp % mod * tmp % mod * a % mod;
	}
	return tmp % mod * tmp % mod; 
}

int solve(int n) {
	long long ans = 0;
	for (int l = 1, r; l <= n; l = r + 1) {
		r = n / (n / l);
		ans = (ans + (mu[r] - mu[l - 1] + (mod - 1)) % (mod - 1) * (n / l) % (mod - 1) * (n / l) % (mod - 1)) % (mod - 1);
	}
	return ans % (mod - 1);
}

int main() {
	cin >> n;
	init();
	long long sum = qpow((long long)fac * fac % mod, n);
	long long ans = 1ll;
	for (int l = 1, r; l <= n; l = r + 1) {
		r = n / (n / l);
		fac = 1ll;
		for (int i = l; i <= r; i++) {
			fac = (long long)fac * i % mod;
		}
		long long t = qpow((long long)fac * fac % mod, solve(n / l));
		ans = (long long)ans * t % mod;
	}
	cout << (long long)sum * qpow(ans, mod - 2) % mod << endl;
	return 0;
}

Atcoder [ABC162E] Sum of gcd of Tuples (Hard)

求：

\[\sum\limits_{i=1}^ki\sum\limits_{a_1=1}^k\sum\limits_{a_2=1}^k\cdots\sum\limits_{a_n=1}^k[\gcd(\{a_n\})=i] \]

令：

\[p=\left\lfloor\frac{k}{i}\right\rfloor \]

有：

\[\sum\limits_{i=1}^ki\sum\limits_{a_1=1}^p\sum\limits_{a_2=1}^p\cdots\sum\limits_{a_n=1}^p[\gcd(\{a_n\})=1] \]

变成 Luogu P3172 [CQOI2015] 选数（见上文）：

\[\sum\limits_{i=1}^ki\sum\limits_{t=1}^p\mu(t)(\left\lfloor\frac{p}{t}\right\rfloor)^n \]

代回：

\[\sum\limits_{i=1}^ki\sum\limits_{t=1}^{\left\lfloor\frac{k}{i}\right\rfloor}\mu(t)(\left\lfloor\frac{\left\lfloor\frac{k}{i}\right\rfloor}{t}\right\rfloor)^n \]

时间复杂度 \(\mathcal{O(k \log k)}\)。

点击查看代码

#include <bits/stdc++.h>
#define int long long

using namespace std;

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

int cnt, n, k, pri[N + 100], mu[N + 100], sum;
bool vis[N + 100];

void init() {
	mu[1] = 1;
	for (int i = 2; i <= N; i++) {
		if (!vis[i]) {
			pri[++cnt] = i;
			mu[i] = -1;
		}
		for (int j = 1; j <= cnt && i * pri[j] <= N; j++) {
			vis[i * pri[j]] = 1;
			if (i % pri[j]) {
				mu[i * pri[j]] = -mu[i];  
			} else {
				mu[i * pri[j]] = 0;
				break;
			}
		}
	}
}

int qpow(int a, int b) {
	if (!b) {
		return 1;
	}
	int tmp = qpow(a, b / 2);
	if (b & 1) {
		return tmp % mod * tmp % mod * a % mod;
	}
	return tmp % mod * tmp % mod;
}

signed main() {
	cin >> n >> k;
	init();
	for (int i = 1; i <= k; i++) {
		int ans = 0;
		for (int j = 1; j <= k / i; j++) {
			ans = (ans + mu[j] * qpow(k / i / j, n) % mod) % mod;
		}
		sum = (sum + i * ans % mod) % mod; 
	}
	cout << sum << endl;
	return 0;
}