莫比乌斯反演练习归纳
这里主要是讲题目，然后提炼出一些有用的经验。统一钦定 \(n\leq m\)。
HDU - 1695 GCD
题意概述
多次询问（大约 \(3000\) 次），每次给出 \(a,b,c,d,k\)，求：
其中，\(1\leq a\leq b\leq 10^5,1\leq c\leq d\leq 10^5,0\leq k\leq 10^5\)。
分析
这是一道经典题目，容斥求一下就行了。
转化成求：
这个是好求的，容易写成（钦定 \(n<m\)）：
不难得到：
后面那个是好求的，可以 \(\mathcal{O}(1)\)。
然后数论分块即可，或者这道题目直接暴力就可以了，不过需要注意这个 \(k=0\)，特判一下就可以了。
代码
#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <stdlib.h>
#include <vector>
#define int long long
#define N 100005
using namespace std;
int cnt,prime[N],mu[N];
bool vis[N];
void init(int n) {
mu[1] = 1;
for (int i = 2;i <= n;i ++) {
if (!vis[i]) prime[++cnt] = i,mu[i] = -1;
for (int j = 1;j <= cnt && prime[j] * i <= n;j ++) {
vis[i * prime[j]] = 1;
if (i % prime[j] == 0) {
mu[i * prime[j]] = 0;
break;
}
mu[i * prime[j]] = -mu[i];
}
}
// for (int i = 1;i <= 10;i ++) cout << mu[i] << ' ';
}
int gcd(int a,int b) {
return b ? gcd(b,a % b) : a;
}
int T,k;
int calc(int n,int m) {
if (k == 0) return 0;
if (n > m) swap(n,m);//n<=m
int res = 0;
// for (int i = 1;i <= n / k;i ++)
// for (int j = i;j <= m / k;j ++)
// res += (gcd(i,j) == 1);
// if (k)
for (int i = 1;i <= n / k;i ++) {
int tn = n / k,tm = m / k;
tn /= i,tm /= i;
res += mu[i] * ((2 * tm - tn + 1) * tn / 2);
}
return res;
}
signed main(){
init(1e5);
cin >> T;
for (int cntt = 1;cntt <= T;cntt ++) {
int a,b,c,d;
scanf("%lld%lld%lld%lld%lld",&a,&b,&c,&d,&k);
// int res = 0;
// for (int i = a;i <= b;i ++)
// for (int j = max(c,i);j <= d;j ++) res += (gcd(i,j) == k);
// printf("%lld\n",res);
printf("Case %lld: %lld\n",cntt,calc(b,d) - calc(a - 1,d) - calc(b,c - 1) + calc(a - 1,c - 1));
}
return 0;
}
HDU - 4746 Mophues
题目概述
设一个数为 \(x=\prod_i p_i\)，如果 \(x\) 分解出来的最多质因数的个数为 \(k\) 且不大于 \(P\) 则称此 \(x\) 为 \(P\) 幸运数。
有 \(Q\) 次询问，每次询问给出 \(n,m,P\)，并定义 \(f(x)=k\)（含义见上），求出：
其中 \(1\leq n,m,P\leq 5\times 10^5,1\leq Q\leq 5000\)。
分析
显然地：
想要将套着的东西拿出来就需要枚举嘛，得到：
不难得到：
套路地令 \(T=dp\) 则：
将 \(T\) 提到前面：
设 \(g_P(t)=\sum_{d\mid T}[f(d)\leq P]\mu(\frac t d)\)。
考虑可以按 \(P\) 升序处理询问，考虑如何 \(g_P(t)\rightarrow g_{P+1}(t)\)。
可以有：\(g_{P+1}(t)=g_{P}(t)+\sum_{d\mid t}\mu(\frac t d)[f(d)=P+1]\)。
考虑到用的是尺取法，而且每个数字对应唯一一个 \(f(x)\)，且每个数字最多只会被用一次，所以说，这里复杂度总共为 \(\mathcal{O}(n\log n)\)。
如何加这个 \(\mu(x)\) 呢？直接存几个桶 \(rtj_x\) 表示满足 \(f(y)=x\) 的 \(y\) 组成的集合，然后直接枚举再枚举倍数就可以了。
但是我们求答案的时候是前缀和，所以说我们这种散点求和不好维护，因此我们使用单点修改，区间查询的树状数组就好了。
于是这题就做完了。
代码
时间复杂度 \(\mathcal{O}(n\log^2 n)\)。
#include <iostream>
#include <cstdio>
#include <cstring>
#include <stdlib.h>
#include <algorithm>
#include <vector>
#define int long long
#define N 500005
using namespace std;
int gcd(int a,int b) {
return b ? gcd(b,a % b) : a;
}
int cnt,prime[N],pri[N],mu[N],tj[N];
vector<int> rtj[25];
int ans[N];
bool vis[N];
void init(int n){
mu[1] = 1;
for (int i = 2;i <= n;i ++) {
if (!vis[i]) prime[++cnt] = i,mu[i] = -1;
for (int j = 1;j <= cnt && prime[j] * i <= n;j ++) {
vis[i * prime[j]] = 1;
pri[i * prime[j]] = prime[j];
if (i % prime[j] == 0) break;
mu[i * prime[j]] = -mu[i];
}
}
for (int i = 2;i <= n;i ++) {
int p = i;
tj[i] = 1;
while(pri[p]) {
tj[i] ++;
p /= pri[p];
}
}
for (int i = 1;i <= n;i ++) rtj[tj[i]].push_back(i);
}
struct node{
int n,m,p,id;
}q[N];
int tr[N];
void update(int x,int val,int n) {
for (;x <= n;x += x & -x) tr[x] += val;
}
int query(int x) {
int res = 0;
for (;x;x -= x & -x) res += tr[x];
return res;
}
int querysum(int l,int r) {
return query(r) - query(l - 1);
}
signed main(){
init(5e5);
int T;
cin >> T;
for (int i = 1;i <= T;i ++) scanf("%lld%lld%lld",&q[i].n,&q[i].m,&q[i].p),q[i].id = i;
stable_sort(q + 1,q + 1 + T,[](node x,node y) {
return x.p < y.p;
});
int nowp = -1;
for (int i = 1;i <= T;i ++) {
int n = q[i].n,m = q[i].m,qid = q[i].id,qp = q[i].p;
while(nowp < qp) {
nowp ++;
for (auto j : rtj[nowp])
for (int t = j;t < N;t += j) update(t,mu[t / j],N - 1);
}
if (n > m) swap(n,m);
int res = 0;
for (int l = 1,r;l <= n;l = r + 1) {
r = min(n / (n / l),m / (m / l));
int ret = 0;
// for (int j = l;j <= r;j ++) ret += g[j];
res += querysum(l,r) * (n / l) * (m / l);
}
ans[qid] = res;
}
for (int i = 1;i <= T;i ++) printf("%lld\n",ans[i]);
return 0;
}
HDU - 4675 GCD of Sequence
题目概述
给出 \(n\) 个元素 \(a_1,a_2,\dots,a_n\) 和 \(m\)，并满足 \(a_i\in[1,m]\)。
询问有多少个序列 \(b\) 满足：
- \(1\leq b_i\leq m\)。
- \(\gcd\{b_n\}=d\)。
- 恰好有 \(k\) 个位置满足 \(a_i\ne b_i\)。
多测。数据范围：\(1\leq n,m\leq 3\times 10^5\)。
分析
感觉经典问题，先设 \(f_{i}\) 表示当序列 \(b\) 的 \(\gcd\) 为 \(i\) 的倍数的方案。
现在考虑第 \(3\) 个限制，记 \(cnt_i\) 表示 \(a\) 中有多少数是 \(i\) 的倍数。
那么我们就可以得到 \(f_i\) 的表达式了，首先我们肯定要选择 \(n-k\) 个位置与 \(a\) 相同，那么方案为 \(\binom{cnt_i}{n-k}\)，这些只有一个选择，其他位置的选择考虑分成当前位置的数值是否是 \(i\) 的倍数。
如果是 \(i\) 的倍数，那么我们要使得 \(b\) 在这个位置上不同，选择为 \(\left\lfloor\frac{m}i\right\rfloor\)，而这样的位置有 \((cnt_i-n+k)\) 个，因此方案为 \(\left(\left\lfloor\frac{m}i-1\right\rfloor\right)^{cnt_i-n+k}\)。
如果不是 \(i\) 的倍数，那么选择就不需要减去 \(1\) 了，而这样的位置有 \(n-cnt_i\) 个。
因此：
考虑如何求答案 \(ans_i\) 表示 \(\gcd\) 恰好为 \(i\) 的方案。
那么显然满足：
然后从大到小求 \(ans\) 就可以了。
代码
时间复杂度 \(\mathcal{O}(\sum m\log m)\)。
#include <iostream>
#include <cstring>
#include <cstdio>
#include <algorithm>
#include <stdlib.h>
#include <vector>
#define int long long
#define N 300005
using namespace std;
const int mod = 1e9 + 7;
int qpow(int a,int b) {
int res = 1;
while(b) {
if (b & 1) res = res * a % mod;
a = a * a % mod;
b >>= 1;
}
return res;
}
int n,m,k,mp[N],cnt[N],inv[N],jc[N];
int C(int a,int b) {
if (a < 0 || b < 0 || a < b) return 0;
return jc[a] * inv[b] % mod * inv[a - b] % mod;
}
int a[N],ans[N];
signed main(){
jc[0] = jc[1] = inv[0] = inv[1] = 1;
for (int i = 2;i < N;i ++) jc[i] = jc[i - 1] * i % mod,inv[i] = (mod - mod / i) * inv[mod % i] % mod;
for (int i = 2;i < N;i ++) inv[i] = inv[i - 1] * inv[i] % mod;
while(~scanf("%lld%lld%lld",&n,&m,&k)) {
for (int i = 1;i <= n;i ++) scanf("%lld",&a[i]),mp[a[i]] ++;
for (int d = 1;d <= m;d ++)
for (int t = d;t <= m;t = t + d) cnt[d] += mp[t];
for (int i = m;i;i --) {
if (cnt[i] >= n - k) ans[i] = C(cnt[i],n - k) * qpow(m / i - 1,cnt[i] - n + k) % mod * qpow(m / i,n - cnt[i]) % mod;//注意条件
for (int j = i + i;j <= m;j += i) ans[i] -= ans[j],ans[i] = (ans[i] + mod) % mod;
}
for (int i = 1;i <= m;i ++) printf("%lld ",ans[i]);
putchar('\n');
for (int i = 1;i <= m;i ++) ans[i] = cnt[i] = mp[i] = 0;
}
return 0;
}
拓展
事实上，若 \(f_i\) 表示的是倍数，而 \(F_i\) 表示的是恰好，那么满足：
这由莫比乌斯反演定理直接得到。
这个相较于上述方法的优势在于可能优化到 \(\mathcal{O}(n)\) 以下。
HDU - 4676 Sum of Gcd
题目概述
给出 \(a\) 以及 \(Q\) 组 \((L,R)\) 求：
多测大约 \(10\) 组，数据范围：\(1\leq n,Q \leq2\times 10^4\)。
分析
这道题目确实难想。
首先看到完整的 \(\gcd\) 考虑用欧拉函数代替，于是乎：
考虑提到前面去，将最大公约数转化成约数：
观察后面的式子，可以发现其本质就是要知道含有 \(d\) 这个约数的 \(a_i(i\in[L,R])\) 有多少个，然后任选两个。
于是我们处理出 \(cnt_i\) 表示在当前区间 \([L,R]\) 中，有多少个数能被 \(i\) 整除，那么答案就为：
现在将目光放在询问区间上，发现可以离线，考虑莫队。
当我们增加一个数 \(x\) 时，如果增加了一个约数 \(d\)，那么我们想要将 \(\varphi(d)\binom{cnt_d}2\) 变成 \(\varphi(d)\binom{cnt_d+1}2\)，而这里只需加一个 \(\varphi(d)\binom{cnt_d}1\) 即可，这是杨辉三角。
少一个数差不多。
于是这道题目就做完了。
代码
分析枚举其因数，考虑均摊，那么 \(\frac 1 n\sum_{i=1}^nd(i)=\frac 1 n\sum_{i=1}^n\frac n i\approx\log n\)，两个都可以得到。
于是时间复杂度为 \(\mathcal{O}(n\sqrt n\log n)\)。
#include <iostream>
#include <cstdio>
#include <cstring>
#include <stdlib.h>
#include <algorithm>
#include <vector>
#define int long long
#define N 20004
using namespace std;
int cntp,prime[N],phi[N];
vector<int> d[N];
bool vis[N];
void init(int n) {
phi[1] = 1;
for (int i = 2;i <= n;i ++) {
if (!vis[i]) prime[++cntp] = i,phi[i] = i - 1;
for (int j = 1;j <= cntp && prime[j] * i <= n;j ++) {
vis[i * prime[j]] = 1;
if (i % prime[j] == 0) {
phi[i * prime[j]] = phi[i] * prime[j];
break;
}
phi[i * prime[j]] = (prime[j] - 1) * phi[i];
}
}
// for (int i = 1;i <= 10;i ++) cout << phi[i] << ' ';
for (int i = 1;i <= n;i ++)
for (int j = i;j <= n;j += i) d[j].push_back(i);
}
int gcd(int a,int b) {
return b ? gcd(b,a % b) : a;
}
int T,n,a[N];
struct node{
int l,r,p,id;
}q[N];
int cnt[N],sum,ans[N];
int m;
void add(int x) {
x = a[x];
for (auto i : d[x]) {
sum += cnt[i] * phi[i];
cnt[i] ++;
}
}
void del(int x) {
x = a[x];
for (auto i : d[x]) {
sum -= (cnt[i] - 1) * phi[i];
cnt[i] --;
}
}
void solve(){
scanf("%lld",&n);
for (int i = 1;i <= n;i ++) scanf("%lld",&a[i]);
scanf("%lld",&m);
int len = ceil(sqrt(n));
for (int i = 1;i <= m;i ++) scanf("%lld%lld",&q[i].l,&q[i].r),q[i].p = q[i].l / len,q[i].id = i;
stable_sort(q + 1,q + 1 + m,[](node x,node y) {
if (x.p == y.p) return (x.p & 1) ? x.r > y.r : x.r < y.r;
return x.p < y.p;
});
sum = 0;
for (int i = 1;i <= n;i ++) cnt[i] = 0;
int l = 1,r = 0;
for (int i = 1;i <= m;i ++) {
int ql = q[i].l,qr = q[i].r,qid = q[i].id;
while(l < ql) del(l ++);
while(l > ql) add(-- l);
while(r < qr) add(++ r);
while(r > qr) del(r --);
ans[qid] = sum;
}
for (int i = 1;i <= m;i ++) printf("%lld\n",ans[i]);
}
signed main(){
init(2e4);
cin >> T;
for (int cntt = 1;cntt <= T;cntt ++) {
printf("Case #%lld:\n",cntt);
solve();
}
return 0;
}
P2257 YY的GCD
题目概述
给定 \(N, M\)，求 \(1 \leq x \leq N\)，\(1 \leq y \leq M\) 且 \(\gcd(x, y)\) 为质数的 \((x, y)\) 有多少对。
数据范围：\(T=10^4,N,M\leq 10^7\)。
分析
经典题目，简单写写。
考虑枚举质数，转化成艾佛森括号形式从而莫比乌斯反演。
于是乎：
不难得到：
套路地令 \(T=dP\) 则：
后面那个可以预处理，然后整除分块就行了。
代码
时间复杂度 \(\mathcal{O}(n+T\sqrt n)\)。
#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <stdlib.h>
#include <vector>
#define int long long
#define N 10000007
#define M 5000006
using namespace std;
bool vis[N];
int prime[M],mu[N],cnt,g[N];
void init(int n){
mu[1] = 1;
for (int i = 2;i <= n;i ++) {
if (!vis[i]) prime[++cnt] = i,mu[i] = -1;
for (int j = 1;j <= cnt && prime[j] * i <= n;j ++) {
vis[i * prime[j]] = 1;
if (i % prime[j] == 0) break;
mu[i * prime[j]] = -mu[i];
}
}
for (int j = 1;j <= cnt;j ++)
for (int t = prime[j];t <= n;t += prime[j]) g[t] += mu[t / prime[j]];
for (int i = 1;i <= n;i ++) g[i] += g[i - 1];
}
int T;
int n,m;
signed main(){
init(1e7);
cin >> T;
for (;T--;) {
scanf("%lld%lld",&n,&m);
if (n > m) swap(n,m);
int res = 0;
for (int l = 1,r;l <= n;l = r + 1) {
r = min(n / (n / l),m / (m / l));
res = (res + (n / l) * (m / l) * (g[r] - g[l - 1]));
}
printf("%lld\n",res);
}
return 0;
}
P3312 [SDOI2014] 数表
题目概述
给定 \(n,m,a\) 求：
数据范围：多测，\(1\leq n,m\leq 10^5,1\leq Q\leq 2\times 10^4\)。
分析
先不考虑 \(a\) 的限制，因此我们只需要保留 \(\sigma(d)\) 即可，在最后考虑。
套路地：
于是有：
套路地令 \(T=dp\) 则：
现在加入 \(a\) 的限制，我们发现对 \(a\) 只需要尺取法即可，跟幸运数那题差不多，还是需要树状数组。
于是这题就做完了。
代码
时间复杂度 \(\mathcal{O}(Q\log^2 n + Q\sqrt n)\)。
#include <iostream>
#include <cstdio>
#include <algorithm>
#include <stdlib.h>
#include <cstring>
#include <vector>
#define int long long
#define N 100005
using namespace std;
const int mod = 2147483648;
int cnt,prime[N],mu[N],sp[N],sigma[N];
bool vis[N];
struct Node{
int id,sigma;
}box[N];
void init(int n) {
mu[1] = sp[1] = sigma[1] = 1;
for (int i = 2;i <= n;i ++) {
if (!vis[i]) prime[++cnt] = i,mu[i] = -1,sp[i] = sigma[i] = i + 1;
for (int j = 1;j <= cnt && prime[j] * i <= n;j ++) {
vis[i * prime[j]] = 1;
int t = i * prime[j];
if (i % prime[j] == 0) {
mu[t] = 0;
sp[t] = sp[i] * prime[j] + 1;
sigma[t] = sigma[i] / sp[i] * sp[t];
break;
}
mu[t] = -mu[i];
sp[t] = prime[j] + 1;
sigma[t] = sigma[i] * (prime[j] + 1);
}
}
for (int i = 1;i <= n;i ++) box[i] = (Node){i,sigma[i]};
stable_sort(box + 1,box + 1 + n,[](Node x,Node y) {
return x.sigma < y.sigma;
});
}
int tr[N];
void update(int x,int val,int n) {
for (;x <= n;x += x & -x) tr[x] += val;
}
int query(int x) {
int res = 0;
for (;x;x -= x & -x) res += tr[x];
return res;
}
struct node{
int n,m,a,id;
}q[N];
int ans[N];
signed main(){
int mxn = 1e5;
init(1e5);
int T;
cin >> T;
for (int i = 1;i <= T;i ++) scanf("%lld%lld%lld",&q[i].n,&q[i].m,&q[i].a),q[i].id = i;
stable_sort(q + 1,q + 1 + T,[](node x, node y) {
return x.a < y.a;
});
int pa = 1;
for (int i = 1;i <= T;i ++) {
int n = q[i].n,m = q[i].m,qa = q[i].a,qid = q[i].id;
if (n > m) swap(n,m);
while(pa <= mxn && box[pa].sigma <= qa) {
int id = box[pa].id;
for (int t = id;t <= mxn;t += id) update(t,(mu[t / id] * sigma[id] % mod + mod) % mod,mxn);
pa ++;
}
int res = 0;
for (int l = 1,r;l <= n;l = r + 1) {
r = min(n / (n / l),m / (m / l));
int s = (query(r) - query(l - 1)) % mod;
s = (s + mod) % mod;
res = (res + s * (n / l) % mod * (m / l) % mod) % mod;
}
ans[qid] = res;
}
for (int i = 1;i <= T;i ++) printf("%lld\n",ans[i]);
return 0;
}
BZOJ - 2693 jzptab
题目概述
求：
数据范围：\(T\leq10^4,1\leq n,m\leq 10^7\).
分析
展开：
想要把 \(\gcd\) 提出，显然枚举：
得到：
其中 \(S(x)\) 代表 \(\sum_{i=1}^x i\)。
套路地令 \(T=dp\) 则：
由对称性化简得到：
设 \(g(T)=T\sum_{d\mid T}\mu(d)d\)，不难发现是一个积性函数。
证明：
若 \(a,b\) 互质则：
而：
于是：
得证。
知道：\(g(p)=p-p^2,g(p^k)p=g(p^{k+1})\)，因此可以直接线性筛。
于是这道题目做完了。
代码
时间复杂度 \(\mathcal{O}(T\sqrt n)\)。
#include <iostream>
#include <cstdio>
#include <stdlib.h>
#include <cstring>
#include <algorithm>
#include <vector>
#define int long long
#define N 10000007
#define M 5000005
using namespace std;
const int mod = 1e8 + 9;
int prime[M],cnt,sum[N],f[N],sum2[N];
bool vis[N];
void init(int n) {
f[1] = 1;
for (int i = 2;i <= n;i ++) {
if (!vis[i]) prime[++cnt] = i,f[i] = i * (1 - i) % mod;
for (int j = 1;j <= cnt && prime[j] * i <= n;j ++) {
vis[prime[j] * i] = 1;
if (i % prime[j] == 0) {
f[i * prime[j]] = f[i] * prime[j] % mod;
break;
}
f[i * prime[j]] = f[prime[j]] * f[i] % mod;
}
}
for (int i = 1;i <= n;i ++) sum[i] = (sum[i - 1] + f[i]) % mod;
for (int i = 1;i <= n;i ++) sum2[i] = (sum2[i - 1] + i) % mod;
}
int n,m;
signed main(){
init(1e7);
int T;
cin >> T;
for (;T--;) {
scanf("%lld%lld",&n,&m);
if (m < n) 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 + sum2[n / l] * sum2[m / l] % mod * (sum[r] - sum[l - 1] + mod) % mod + mod) % mod;
}
printf("%lld\n",ans);
}
return 0;
}
拓展：线性筛求积性函数 \(g(x)\) 的充要条件
我们需要考虑三个部分：
- 为质数，直接代入（已知）。
- \(i,p\) 互质，直接用积性函数。
- \(p\) 为 \(i\) 的因子，那么此时设 \(i=p^km\)，则 \(g(i)=g(p^k)g(m),g(ip)=g(p^{k+1})g(m)\)，那么 \(g(ip)=g(i)\frac{g(p^{k+1})}{g(p^k)}\)。
所以只需知道 \(\frac{g(p^{k+1})}{g(p^k)}\) 是一个什么值就可以了（或者说你可以用其他方法得到他）。
P2231 [HNOI2002] 跳蚤
之前写过题解：https://www.cnblogs.com/high-sky/p/19162197。
P3172 [CQOI2015] 选数
题目概述
选 \(n\) 个数满足在 \([L,H]\) 中，请回答选出数的方案满足他们的 \(\gcd\) 为 \(k\) 的答案，并对 \(10^9+7\) 取模。
数据范围：\(1\leq L,H,n,k\leq 10^9,1\leq H-L\leq 10^5\)。
分析
我之前写的题解：https://www.cnblogs.com/high-sky/p/19154610
但是这里有一点需要补充，为什么要减去全部相同的情况呢？
这是因为我们把递推范围变小了，所以变换后有些 \(\gcd\) 不在这个区间里面，所以容斥不掉，因此减去即可，在最后如果有 \(1\) 就加上即可。
举个例子：假如变换后的区间为 \(L'=100,H'=200\)，而我们的递推是 \([1,100]\)，如果此时有 \(\gcd=150\)，所以容斥不掉，而这种方案对答案影响不大，故可以单独处理。
代码
时间复杂度 \(\mathcal{O}((H-L+1)\log (H-L+1))\)。
#include <iostream>
#include <cstring>
#include <stdlib.h>
#include <cstdio>
#include <algorithm>
#include <vector>
#define int long long
#define N 100005
using namespace std;
const int mod = 1e9 + 7;
int qpow(int a,int b) {
int res = 1;
while(b) {
if (b & 1) res = res * a % mod;
a = a * a % mod;
b >>= 1;
}
return res;
}
int n,k,l,h,f[N],ans[N];
signed main(){
cin >> n >> k >> l >> h;
h /= k;
l = (l + k - 1) / k;
if (l > h) return cout << 0,0;
k = 1;
for (int i = 1;i <= h - l;i ++) {
int L = (l + i - 1) / i,R = h / i;
if (L > R) continue;
f[i] = qpow(R - L + 1,n) - (R - L + 1);
f[i] %= mod;
f[i] = (f[i] + mod) % mod;
}
for (int i = h - l;i;i --) {
ans[i] = f[i];
for (int j = i + i;j <= h - l;j += i) ans[i] = (ans[i] - ans[j] + mod) % mod;
}
cout << (ans[1] + (l == 1)) % mod;
return 0;
}
P3327 [SDOI2015] 约数个数和
题意概述
设 \(d(x)\) 为 \(x\) 的约数个数，给定 \(n,m\)，求：
数据范围：\(1\le T,n,m \le 50000\)。
分析
考虑这个 \(d(ij)\) 怎么分析。
设 \(i=\prod_k p_k^{\alpha_k},j=\prod_k p_k^{\beta_k}\)，那么：\(d(ij)=\prod_k(\alpha_k+\beta_k+1)\)。
我们考虑什么时候还会变成 \(\prod_k(\alpha_k+\beta_k+1)\)（尽管这有点从答案推过程）。
我们往 \(\gcd\) 那方面想，对于 \(i,j\) 中的 \(p_k\)，我可以选 \(1,2,\dots,\alpha_k\) 个，此时另外一个不选，那么就是 \((\alpha_k+\beta_k)\)，那为什么 \(+1\)，两个都不选呗。
这种情况的产生可以在 \(ij\) 的约数分成两个部分，一个从 \(i\) 来一个从 \(j\) 来，那么他们两个必须互质，所以：
转化一下不难得到：
经典式子变换一下：
后面的那个预处理一下就行了，然后整除分块即可。
代码
时间复杂度 \(\mathcal{O}(T\sqrt n)\)。
#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <stdlib.h>
#include <vector>
#define int long long
#define N 50005
using namespace std;
int T;
int gcd(int a,int b) {
return b ? gcd(b,a % b) : a;
}
int prime[N],cnt,mu[N],sums[N],summu[N];
bool vis[N];
int S(int n) {
int res = 0;
for (int l = 1,r;l <= n;l = r + 1) {
r = n / (n / l);
res += (n / l) * (r - l + 1);
}
return res;
}
void init(int n) {
mu[1] = 1;
for (int i = 2;i <= n;i ++) {
if (!vis[i]) prime[++cnt] = i,mu[i] = -1;
for (int j = 1;j <= cnt && prime[j] * i <= n;j ++) {
vis[i * prime[j]] = 1;
if (i % prime[j] == 0) break;
mu[i * prime[j]] = -mu[i];
}
}
for (int i = 1;i <= n;i ++) summu[i] = summu[i - 1] + mu[i];
for (int i = 1;i <= n;i ++) sums[i] = S(i);
}
signed main(){
init(5e4);
cin >> T;
for (int n,m;T--;) {
scanf("%lld%lld",&n,&m);
int ans = 0;
if (n > m) swap(n,m);
for (int l = 1,r;l <= n;l = r + 1) {
r = min(n / (n / l),m / (m / l));
int t1 = n / l,t2 = m / l;
ans += (summu[r] - summu[l - 1]) * sums[t1] * sums[t2];
}
printf("%lld\n",ans);
}
return 0;
}
浙公网安备 33010602011771号