莫比乌斯反演 学习笔记
一、前置技能
- 积性函数
定义\(f(n)\)是积性函数，如果满足\(f(nm)=f(n)f(m), gcd(n,m) = 1\)
\(f(n)\)称为完全积性函数，当且仅当满足\(f(nm) = f(n)f(m), n,m\in Z\)
数论函数\(f(n)\)是积性函数的充分必要条件是：
\(f(1) = 1\) 且 \(f(n) = f(p_1^{a_1})f(p_2^{a_2})\cdots f(p_s^{a_s})\)
数论函数\(f(n)\)是完全积性函数的充分必要条件是：
\(f(1) = 1\) 且 \(f(n) = f^{a_1}(p_1)f^{a_2}(p_2) \cdots f^{a_s}(p_s)\)
若\(f(x)\)和\(g(x)\)是积性函数则：
- \(h(x) = f(x^p)\)
- \(h(x) = f^p(x)\)
- \(h(x) = f(x)g(x)\)
- \(h(x) = \sum_{d|x}f(d)g(\frac{x}{d})\)
也是是积性函数
- 常用数论函数
- 单位函数：\(\epsilon(n)=[n=1]\)（当且仅当\(n=1\)时值为\(1\)，其他情况值为\(0\)）
- 恒等函数：\(ID(n) = n\)（就是自身）
- 常数函数：\(1(n) = 1\)（都是\(1\)）
- 除数函数：\(d(n) = (a_1+1)\cdots(a_s+1)\)（其中\(n=\prod_{i=1}^{s}p_i^{a_i}\)，表示\(n\)的因子数量）
- 欧拉函数：\(\phi(n) = \sum_{i=1}^{n}[gcd(i,n)=1]\)（欧拉函数表示\(1\)到\(n\)之间和\(n\)的最大公约数为\(1\)的数的数量）
- 莫比乌斯函数：\(\mu(n) = \begin{cases} 1 & n=1 \\ 0 & \exists p > 1\ and\ p^2|n \\ (-1)^{\omega(n)} & otherwise\end{cases}\)
（\(\omega(n)\)表示\(n\)的不同质因子数量,\(p\)为素因子）
-
\(Dirichlet\)卷积
\((f*g)(n) = \sum_{d|n}f(d)g(\frac{n}{d})\)
\(Dirichlet\)卷积满足交换律和结合律
\(\epsilon 为Dirichlet卷积单位元，满足f*\epsilon = f\) -
常用\(Dirichlet\)卷积函数
- \(\epsilon = \mu * 1 \Leftrightarrow \epsilon(n) = \sum_{d|n}\mu(d)\)
- \(f = f * \epsilon \Leftrightarrow f(n) = \sum_{d|n}f(d)\cdot \epsilon(\frac{n}{d})\)（任意函数\(f\)和单位函数卷积都为自身，因为单位函数只有在等于\(1\)的时候才是\(1\)）
- \(\phi = \mu * id \Leftrightarrow \phi(n) = \sum_{d|n}d\cdot \mu(\frac{n}{d})\)
- \(id = \phi * 1 \Leftrightarrow n = \sum_{d|n}\phi(d)\)
- \(d = 1 * 1 \Leftrightarrow d(n) = \sum_{d|n} 1\)
- \(\frac{\phi(n)}{n} = \frac{\mu(n)}{n} * 1 \Leftrightarrow \sum_{d|n}\frac{\mu(d)}{d} \cdot \frac{n}{d}\)
- 关于整除
- \(\lfloor\frac{n}{a\cdot b} \rfloor = \lfloor\frac{\lfloor\frac{n}{a} \rfloor}{b} \rfloor\)
证明：
\[\lfloor \frac n{ab} \rfloor = \lfloor \frac na\cdot \frac1b \rfloor = \lfloor (\lfloor \frac na \rfloor + r) \cdot \frac1b \rfloor = \lfloor\frac{\lfloor\frac{n}{a} \rfloor}{b} + \frac rb \rfloor =\lfloor\frac{\lfloor\frac{n}{a} \rfloor}{b} \rfloor
\]
- 数论分块：对于任意\(i\le n\)，最大的满足\(\lfloor \frac ni \rfloor =\lfloor \frac nj\rfloor\)的\(j=\lfloor\frac{n}{\lfloor\frac ni \rfloor}\rfloor\)
证明：
\[\lfloor \frac ni \rfloor \le \frac ni\Rightarrow \frac n{\lfloor \frac ni \rfloor} \ge \frac n{\frac ni} = i \Rightarrow \lfloor \frac n{\lfloor \frac ni \rfloor}\rfloor \ge \lfloor i \rfloor = i
\]
\[\Rightarrow i \le \lfloor\frac{n}{\lfloor\frac ni \rfloor}\rfloor
\]
\[\Rightarrow j = \lfloor\frac{n}{\lfloor\frac ni \rfloor}\rfloor
\]
二、莫比乌斯函数
定义莫比乌斯函数为：\(\mu(n) = \begin{cases}1 & n = 1 \\ 0 & \exists\ p>1\ and\ p^2|n \\ (-1)^{\omega(n)} & otherwise \end{cases}\)
\(\omega(n)\)为\(n\)的不同的素因子数量,\(p\)为素因子
一个最基本、最重要的性质:
\[\sum_{d|n}\mu(d) = \lfloor \frac{1}{n}\rfloor = \begin{cases} 1 & n=1 \\ 0 & n\ne 1\end{cases}
\]
\[即\mu * 1 = \epsilon
\]
证明：
\[n = 1时显然成立，现设n=p_1^{a_1}\cdots p_s^{a_s},a_i\ge 1
\]
\[\sum_{d|n}\mu(d) = \sum_{i=0}^{s}C(s,i)(-1)^i = (1-1)^s = 0
\]
欧拉反演：
\[\phi(n) = \sum_{d|n}\mu(d)\frac{n}{d}
\]
\[即\mu * id = \phi
\]
证明：
\[\phi(n) = \sum_{i=1}^n[gcd(i.n)=1] = \sum_{i=1}^{n}\sum_{d|gcd(i,n)}\mu(d)
\]
\[=\sum_{d|n}\mu(d)\sum_{i=1}^{n}[d\ |\ i] = \sum_{d|n}\mu(d)\frac{n}{d}
\]
莫比乌斯函数是个积性函数，所以可以直接用线性筛筛出来：
view code
int mu[MAXN];
void sieve(){
vector<bool> pm(MAXN,true);
vector<int> prime;
mu[1] = 1;
for(int i = 2; i < MAXN; i++){
if(pm[i]) prime.push_back(i), mu[i] = -1;
for(int j = 0; i * prime[j] <= MAXN; j++){
mu[i*prime[j]] = -mu[i];
pm[i*prime[j]] = false;
if(i%prime[j]==0){
mu[i*prime[j]] = 0;
break;
}
}
}
}
三、莫比乌斯反演
莫比乌斯反演定理：设\(f(m),F(n)\)是数论函数，那么
\[F(n)=\sum_{d|n}f(d)
\]
成立的充分必要条件是
\[f(n)=\sum_{d|n}\mu(d)F(\frac{n}{d})
\]
利用Dirichlet卷积证明
证明：
\[已知\mu * 1 = \epsilon, f * \epsilon = f
\]
\[F = f * 1\Leftrightarrow F*\mu = f*1*\mu \Leftrightarrow F*\mu = f*(1*\mu) = f
\]
利用定义证明
充分性：
\[如果f(n)=\sum_{d|n}\mu(d)F(\frac{n}{d})成立
\]
\[则F(n)=\sum_{d|n}f(d) = \sum_{d|n}\sum_{k|d}\mu(k)F(\frac{d}{k}) = \sum_{k|n}\mu(k)\sum_{k|d,d|n}F(\frac{d}{k})
\]
\[令d = kl
\]
\[F(n)=\sum_{k|n}\mu(k)\sum_{l|\frac{n}{k}}F(l)=\sum_{l|n}F(l)\sum_{k|\frac{n}{l}}\mu(k)
\]
\[\because \sum_{k|\frac{n}{1}}\mu(k) = [\frac{n}{l}=1] = [l=n]
\]
\[\therefore \sum_{l|n}F(l)\sum_{k|\frac{n}{l}}\mu(k) = F(n)
\]
必要性：
\[若F(n)=\sum_{d|n}f(d)成立
\]
\[则f(n)=\sum_{d|n}\mu(d)F(\frac{n}{d})=\sum_{d|n}\mu(d)\sum_{k|\frac{n}{d}}f(k)=\sum_{k|n}f(k)\sum_{d|\frac{n}{k}}\mu(d)
\]
\[\because \sum_{d|\frac{n}{k}}\mu(d)=[\frac{n}{k}=1]=[k=n]
\]
\[\therefore \sum_{k|n}f(k)\sum_{d}{\frac{n}{k}}\mu(d)=f(n)
\]
另一形式：
\[F(n)=\sum_{n|d}f(d)\Leftrightarrow f(n) = \sum_{n|d}\mu(\frac{d}{n})F(d)
\]
证明：
充分性：
\[如果f(n)=\sum_{n|d}\mu(\frac{d}{n})F(d)成立
\]
\[则F(n) = \sum_{n|d}f(d) = \sum_{n|d}\sum_{d|t}\mu(\frac{t}{d})F(t)=\sum_{n|t}F(t)\sum_{n|d,d|t}\mu(\frac{t}{d})
\]
\[令k=\frac{t}{n}, d = ln
\]
\[\sum_{n|t}F(t)\sum_{n|d,d|t}\mu(\frac{t}{d}) = \sum_{n|t}F(t)\sum_{l|k}\mu(\frac{kn}{ln}) = \sum_{n|t}F(t)\sum_{l|k}\mu(\frac{k}{l})=\sum_{n|t}F(t)\sum_{l|k}\mu(l)
\]
\[\because \sum_{l|k}\mu(l) = [k=1] = [\frac{t}{n}=1]
\]
\[\therefore \sum_{n|t}F(t)\sum_{l|k}\mu(l) = F(n)
\]
必要性
\[如果F(n)=\sum_{n|d}f(d)成立
\]
\[则f(n)=\sum_{n|d}\mu(\frac{d}{n})F(d)=\sum_{n|d}\mu(\frac{d}{n})\sum_{d|t}f(t)
\]
\[=\sum_{n|t}f(t)\sum_{n|d,d|t}\mu(\frac{d}{n})
\]
\[令k = \frac{t}{n}, d = ln
\]
\[\Rightarrow \sum_{n|t}f(t)\sum_{n|d,d|t}\mu(\frac{d}{n}) = \sum_{n|t}f(t)\sum_{l|k}\mu(l)
\]
\[\because \sum_{l|k}\mu(l) = [k=1] = [\frac{t}{n}=1]
\]
\[\therefore \sum_{n|t}f(t)\sum_{l|k}\mu(l) = f(n)
\]
四、例题：
- LuoguP3935 Calculating🔗
题解：
\[f(n)就其实表示除数函数
\]
\[所以就是计算\sum_{i=1}^n\sum_{d|i}1
\]
\[转换枚举顺序\sum_{i=1}^n\sum_{d|i}1 = \sum_{d=1}^n\sum_{d|x,x\le n}1=\sum_{d=1}^n \lfloor\frac{n}{d} \rfloor
\]
数论分块来做就好了
复杂度\(O(\sqrt{r})\)
view code
//#pragma GCC optimize("O3")
//#pragma comment(linker, "/STACK:1024000000,1024000000")
#include<bits/stdc++.h>
using namespace std;
function<void(void)> ____ = [](){ios_base::sync_with_stdio(false); cin.tie(0); cout.tie(0);};
typedef long long int LL;
const LL MOD = 998244353;
LL l, r;
LL solve(LL n){
LL ret = 0;
for(LL l = 1, r; l <= n; l = r + 1){
r = n / (n/l);
ret = (ret + (r-l+1) * (n/l)) % MOD;
}
return ret;
}
int main(){
____();
cin >> l >> r;
cout << (solve(r) - solve(l-1) + MOD) % MOD << endl;
return 0;
}
- LuoguP3455 [POI2007]ZAP-Queries🔗
题解：
\[要求计算\sum_{x=1}^a\sum_{y=1}^b[gcd(x,y)=k]
\]
\[这里假定a\le b
\]
\[\sum_{x=1}^a\sum_{y=1}^b[gcd(x,y)=k] = \sum_{x=1}^{\lfloor \frac ak \rfloor}\sum_{y=1}^{\lfloor \frac bk\rfloor}[gcd(x,y)=1]= \sum_{x=1}^{\lfloor \frac ak \rfloor}\sum_{y=1}^{\lfloor \frac bk\rfloor}\epsilon(gcd(x,y))
\]
\[= \sum_{x=1}^{\lfloor \frac ak \rfloor}\sum_{y=1}^{\lfloor \frac bk\rfloor}\sum_{d|gcd(x,y)}\mu(d)
\]
\[转变为枚举d
\]
\[=\sum_{d=1}^{\lfloor \frac ak \rfloor}\mu(d)\sum_{d|x,x\le \lfloor \frac ak\rfloor}\sum_{d|y,y\le \lfloor \frac bk\rfloor} 1
\]
\[= \sum_{d=1}^{\lfloor \frac ak \rfloor}\mu(d) \lfloor \frac{\lfloor \frac ak\rfloor}{d} \rfloor \cdot \lfloor \frac{\lfloor \frac bk\rfloor}{d} \rfloor
\]
然后处理出来莫比乌斯函数的前缀和用整除分块来做就好了
复杂度\(O(Q\sqrt{max\{\bar{a},\bar{b}\}})\)
view code
//#pragma GCC optimize("O3")
//#pragma comment(linker, "/STACK:1024000000,1024000000")
#include<bits/stdc++.h>
using namespace std;
function<void(void)> ____ = [](){ios_base::sync_with_stdio(false); cin.tie(0); cout.tie(0);};
const int MAXN = 5e4+7;
typedef long long int LL;
LL mu[MAXN];
void sieve(){
vector<bool> pm(MAXN,true);
vector<int> prime;
mu[1] = 1;
for(int i = 2; i < MAXN; i++){
if(pm[i]) prime.push_back(i), mu[i] = -1;
for(int j = 0; i * prime[j] < MAXN; j++){
pm[i*prime[j]] = false;
if(i%prime[j]==0){
mu[i*prime[j]] = 0;
break;
}
mu[i*prime[j]] = -mu[i];
}
}
for(int i = 1; i < MAXN; i++) mu[i] += mu[i-1];
}
void solve(){
int a,b,k;
LL ret = 0;
scanf("%d %d %d",&a,&b,&k);
a /= k, b /= k;
if(a>b) swap(a,b);
for(int l = 1, r; l <= a; l = r + 1){
r = min(a / (a / l), b / (b / l));
ret += (mu[r] - mu[l-1]) * (a/l) * (b/l);
}
printf("%lld\n",ret);
}
int main(){
int T;
sieve();
for(scanf("%d",&T); T; T--) solve();
return 0;
}
- LuoguP4450 双亲数（第二题的简化版）🔗
题解：
和上一题一样
复杂度\(O(\sqrt{max(A,B)})\)
view code
//#pragma GCC optimize("O3")
//#pragma comment(linker, "/STACK:1024000000,1024000000")
#include<bits/stdc++.h>
using namespace std;
function<void(void)> ____ = [](){ios_base::sync_with_stdio(false); cin.tie(0); cout.tie(0);};
const int MAXN = 1e6+7;
typedef long long int LL;
LL mu[MAXN];
void sieve(){
vector<bool> pm(MAXN,true);
vector<int> prime;
mu[1] = 1;
for(int i = 2; i < MAXN; i++){
if(pm[i]) prime.push_back(i), mu[i] = -1;
for(int j = 0; i * prime[j] < MAXN; j++){
pm[i*prime[j]] = false;
if(i%prime[j]==0){
mu[i*prime[j]] = 0;
break;
}
mu[i*prime[j]] = -mu[i];
}
}
for(int i = 1; i < MAXN; i++) mu[i] += mu[i-1];
}
void solve(){
int a,b,k;
LL ret = 0;
scanf("%d %d %d",&a,&b,&k);
a /= k, b /= k;
if(a>b) swap(a,b);
for(int l = 1, r; l <= a; l = r + 1){
r = min(a / (a / l), b / (b / l));
ret += (mu[r] - mu[l-1]) * (a/l) * (b/l);
}
printf("%lld\n",ret);
}
int main(){
sieve();
solve();
return 0;
}
- LuoguP2257 YY的GCD🔗
题解：
我们枚举每个质数\(p\)然后用前两题的办法来计算有多少个\(x,y\)满足\(gcd(x,y)=p\)
这样复杂度是很高的，所以得要继续推柿子
下面假设\(n\le m\)
\[要计算的就是\sum_{p\in prime}\sum_{x = 1}^n\sum_{y=1}^m[gcd(x,y)=p]
\]
\[\sum_{p\in prime}\sum_{x = 1}^n\sum_{y=1}^m[gcd(x,y)=p] = \sum_{p\in prime}\sum_{x=1}^{\lfloor \frac np\rfloor}\sum_{y=1}^{\lfloor \frac mp\rfloor}[gcd(x,y)=1]
\]
\[=\sum_{p\in prime}\sum_{x=1}^{\lfloor \frac np\rfloor}\sum_{y=1}^{\lfloor \frac mp\rfloor}\epsilon(gcd(x,y))=\sum_{p\in prime}\sum_{x=1}^{\lfloor \frac np\rfloor}\sum_{y=1}^{\lfloor \frac mp\rfloor}\sum_{d|gcd(x,y)}\mu(d)
\]
\[=\sum_{p\in prime}\sum_{d=1}^{\lfloor \frac np\rfloor}\mu(d)\sum_{d|x,x\le \lfloor \frac np\rfloor}\sum_{d|y,y\le \lfloor \frac mp\rfloor}1 = \sum_{p\in prime}\sum_{d=1}^{\lfloor \frac np\rfloor}\mu(d)\lfloor \frac n{pd}\rfloor \lfloor \frac m{pd}\rfloor
\]
\[令T = pd，然后枚举T
\]
\[\sum_{p\in prime}\sum_{d=1}^{\lfloor \frac np\rfloor}\mu(d)\lfloor \frac n{pd}\rfloor \lfloor \frac m{pd}\rfloor = \sum_{T=1}^{n}\lfloor \frac nT\rfloor \lfloor \frac mT\rfloor\sum_{p\in prime,p|T}\mu(\frac Tp)
\]
\[发现\sum_{p\in prime,p|T}\mu(\frac Tp)是可以预处理出来的
\]
\[那么对于每次询问就可以数论分块做了
\]
复杂度\(O(MAXN + T\sqrt{max(n,m)})\)
view code
//#pragma GCC optimize("O3")
//#pragma comment(linker, "/STACK:1024000000,1024000000")
#include<bits/stdc++.h>
using namespace std;
function<void(void)> ____ = [](){ios_base::sync_with_stdio(false); cin.tie(0); cout.tie(0);};
const int MAXN = 1e7+7;
typedef long long int LL;
int mu[MAXN];
LL pre[MAXN];
bool pm[MAXN];
vector<int> prime;
void sieve(){
fill(begin(pm),end(pm),true);
mu[1] = 1;
for(int i = 2; i < MAXN; i++){
if(pm[i]) prime.push_back(i), mu[i] = -1;
for(int j = 0; i * prime[j] < MAXN; j++){
pm[i*prime[j]] = false;
if(i%prime[j]==0){
mu[i*prime[j]] = 0;
break;
}
mu[i*prime[j]] = -mu[i];
}
}
for(int p : prime) for(int i = 1; i * p < MAXN; i++) pre[i*p] += mu[i];
for(int i = 1; i < MAXN; i++) pre[i] += pre[i-1];
}
void solve(){
int n,m;
LL ret = 0;
scanf("%d %d",&n,&m);
if(n>m) swap(n,m);
for(int l = 1, r; l <= n; l = r + 1){
r = min(n/(n/l),m/(m/l));
ret += (LL)(n/l) * (m/l) * (pre[r] - pre[l-1]);
}
printf("%lld\n",ret);
}
int main(){
sieve();
int T; for(scanf("%d",&T); T; T--) solve();
return 0;
}
- LuoguP3768 简单的数学题🔗
题解：
不多BB直接开始推柿子
\[\sum_{i=1}^n\sum_{j=1}^ni\cdot j\cdot gcd(i,j)
\]
\[枚举g = gcd(i,j)
\]
\[\Rightarrow \sum_{g=1}^n g\sum_{g|i}\sum_{g|j}i\cdot j\cdot[gcd(\frac ig,\frac jg)=1]=\sum_{g=1}^ng\sum_{g|i}\sum_{g|j}i\cdot j\cdot \epsilon(gcd(\frac ig,\frac jg))
\]
\[=\sum_{g=1}^ng\sum_{g|i}\sum_{g|j}i\cdot j\sum_{d|gcd(\frac ig,\frac jg)}\mu(d)
\]
\[令x=\frac ig,y=\frac jg,枚举d
\]
\[\Rightarrow \sum_{g=1}^ng\sum_{d=1}^{\lfloor \frac ng\rfloor}\mu(d)\sum_{d|x}^{\lfloor \frac ng\rfloor}\sum_{d|y}^{\lfloor \frac ng\rfloor}g^2xy
\]
\[令a = \frac xd,b = \frac yd
\]
\[\Rightarrow \sum_{g=1}^ng^3\sum_{d=1}^{\lfloor\frac ng \rfloor}\mu(d)\sum_{a=1}^{\lfloor \frac n{gd} \rfloor}\sum_{b=1}^{\lfloor \frac n{gd} \rfloor}d^2ab
\]
\[\sum_{i=1}^n\sum_{j=1}^ni\cdot j = \sum_{i=1}^ni\sum_{j=1}^nj = (1+\cdots+n)^2
\]
\[令sum(n) = 1+\cdots+n
\]
\[\Rightarrow \sum_{g=1}^ng^3\sum_{d=1}^{\lfloor\frac ng \rfloor}\mu(d)\sum_{a=1}^{\lfloor \frac n{gd} \rfloor}\sum_{b=1}^{\lfloor \frac n{gd} \rfloor}d^2ab = \sum_{g=1}^ng^3\sum_{d=1}^{\lfloor\frac ng \rfloor}d^2\mu(d)sum^2(\lfloor \frac n{gd} \rfloor)
\]
\[令T = gd,枚举T
\]
\[\Rightarrow \sum_{T=1}^n sum^2(\lfloor \frac n{T} \rfloor)\sum_{d|T}d^2\mu(d)\frac{T^3}{d^3}=\sum_{T=1}^n T^2 sum^2(\lfloor \frac n{T} \rfloor)\sum_{d|T}\mu(d)\frac{T}{d}
\]
\[=\sum_{T=1}^n sum^2(\lfloor \frac n{T} \rfloor)T^2\phi(T)
\]
\[考虑数论分块，对于一块区间,sum^2(\lfloor \frac n{T} \rfloor)很好算，然后就是计算T^2\phi(T)的前缀和了，这是个积性函数，上杜教筛
\]
\[接下来是杜教筛的内容
\]
\[假设F = g*f
\]
\[\sum_{i=1}^nF(i) = \sum_{i=1}^n\sum_{d|i}g(d)f(\frac id)=\sum_{d=1}^ng(d)\sum_{d|i}f(i)=\sum_{d=1}^ng(d)\sum_{x=1}^{\lfloor\frac nd \rfloor}f(x)
\]
\[令S(n)=\sum_{i=1}^{n}f(i)
\]
\[\sum_{d=1}^ng(d)\sum_{x=1}^{\lfloor\frac nd \rfloor}f(x) = \sum_{d=1}^ng(d)S(\lfloor\frac nd \rfloor)
\]
\[\Rightarrow \sum_{i=1}^nF(i)=\sum_{d=1}^ng(d)S(\lfloor\frac nd \rfloor)
\]
\[\Rightarrow g(1)S(n) = \sum_{i=1}^nF(i) - \sum_{d=2}^n g(d)S(\lfloor\frac nd \rfloor)
\]
\[如果能快速得到\sum_{i=1}^nF(i)的值，那么后面的值可以数论分块算出来
\]
\[回到刚才要算的式子：\sum_{T=1}^n sum^2(\lfloor \frac n{T} \rfloor)T^2\phi(T)
\]
\[现在f(n) = n^2\phi(n)，我们需要构造一个g(n)使得\sum_{d|n}d^2\phi(d)g(\frac nd)是个比较好算的式子
\]
\[根据欧拉反演 \phi=\mu*ID\Rightarrow ID=1*\phi
\]
\[我们构造希望把\sum_{d|n}d^2\phi(d)g(\frac nd)中的d约掉，然后剩下的就可以用欧拉反演来做了
\]
\[那么构造出g(n) = n^2
\]
\[那么f*g = \sum_{d|n}d^2\phi(d)\frac{n^2}{d^2} = n^2\sum_{d|n}\phi(d)=n^3
\]
\[所以我们的杜教筛是这样的：S(n)=\sum_{i=1}^ni^3 - \sum_{d=2}^nd^2S(\lfloor\frac nd \rfloor)
\]
\[我们小学二年级的时候就学过：\sum_{i=1}^ni^3=\frac{n^2(n+1)^2}{4}
\]
\[所以最终S(n)=\frac{n^2(n+1)^2}{4} - \sum_{d=2}^nd^2S(\lfloor\frac nd \rfloor)
\]
\[注意后面需要整除分块来算的时候需要计算\sum_{i=l}^r i^2
\]
\[这个我们小学二年级也学过：\sum_{i=1}^ni^2 = \frac{n(n+1)(2n+1)}{6}
\]
\[然后就可以做了！
\]
\[其实有个简单的方法
\]
\[\sum_{i=1}^n\sum_{i=1}^ni\cdot j\cdot gcd(i,j)
\]
\[=\sum_{i=1}^n\sum_{j=1}^ni\cdot j\sum_{d|i,d|j}\phi(d)
\]
\[令x=\frac id,y = \frac id
\]
\[=\sum_{d=1}^nd^2\phi(d)\sum_{x=1}^{\lfloor \frac nd\rfloor}\sum_{y=1}^{\lfloor \frac nd\rfloor}xy
\]
\[\sum_{d=1}^nd^2\phi(d)sum^2(\lfloor \frac nd\rfloor)
\]
view code
//#pragma GCC optimize("O3")
//#pragma comment(linker, "/STACK:1024000000,1024000000")
#include<bits/stdc++.h>
using namespace std;
function<void(void)> ____ = [](){ios_base::sync_with_stdio(false); cin.tie(0); cout.tie(0);};
typedef long long int LL;
const int MAXN = 5e6+7;
int phi[MAXN], prime[MAXN], cnt;
LL p,S[MAXN],inv2,inv6;
unordered_map<LL,LL> MS;
void sieve(){
phi[1] = 1;
for(int i = 2; i < MAXN; i++){
if(!phi[i]) phi[i] = i-1, prime[cnt++] = i;
for(int j = 0; i * prime[j] < MAXN; j++){
if(i%prime[j]==0){
phi[i*prime[j]] = phi[i] * prime[j];
break;
}
phi[i*prime[j]] = phi[i] * phi[prime[j]];
}
}
for(int i = 1; i < MAXN; i++) S[i] = (S[i-1] + 1ll * i * i % p * phi[i] % p) % p;
}
LL ksm(LL a, LL b){
a %= p;
LL ret = 1;
while(b){
if(b&1) ret = ret * a % p;
b >>= 1;
a = a * a % p;
}
return ret;
}
LL inv(LL x){ return ksm(x,p-2); }
LL sum1(LL n){ n %= p; return n * (n + 1) % p * inv2 % p; }
LL sum2(LL n){ n %= p; return n * (n + 1) % p * (2 * n + 1) % p * inv6 % p; }
LL sum3(LL n){ return sum1(n) * sum1(n) % p; }
LL calS(LL n){
if(n<MAXN) return S[n];
if(MS.count(n)) return MS[n];
LL ret = sum3(n);
for(LL l = 2, r; l <= n; l = r + 1){
r = n / (n / l);
ret = (ret - calS(n / l) * (sum2(r) - sum2(l - 1) + p) % p + p ) % p;
}
return MS[n] = ret;
}
void solve(LL n){
LL ret = 0;
for(LL l = 1, r; l <= n; l = r + 1){
r = n / (n / l);
ret = (ret + sum1(n/l) * sum1(n/l) % p * (calS(r) - calS(l-1) + p) % p) % p;
}
printf("%lld\n",ret);
}
int main(){
LL n;
scanf("%lld %lld",&p,&n);
sieve();
inv2 = inv(2), inv6 = inv(6);
solve(n);
return 0;
}
- LuoguP3172 [CQOI2015]选数🔗
题解：
\[计算\sum_{i_1=L}^H\cdots\sum_{i_n=L}^H[gcd(i_1,\cdots,i_n)=K]
\]
\[\Rightarrow \sum_{i_1={\lfloor \frac{L-1}{k}+1\rfloor}}^{\lfloor\frac{H}{K}\rfloor}\cdots\sum_{i_n=\lfloor \frac{L-1}{k}+1\rfloor}^{\lfloor\frac{H}{K}\rfloor}[gcd(i_1,\cdots,i_n)=1]=\sum_{i_1={\lfloor \frac{L-1}{k}+1\rfloor}}^{\lfloor\frac{H}{K}\rfloor}\cdots\sum_{i_n=\lfloor \frac{L-1}{k}+1\rfloor}^{\lfloor\frac{H}{K}\rfloor}\sum_{d|gcd(i_1,\cdots,i_n)}\mu(d)
\]
\[令\lfloor\frac {L-1}K+1 \rfloor = l,\lfloor\frac HK \rfloor = r,枚举d
\]
\[\Rightarrow \sum_{d=1}^{\lfloor\frac HK \rfloor }\mu(d)\sum_{d|i_1}\cdots\sum_{d|i_n}1=\sum_{d=1}^{\lfloor\frac HK \rfloor }\mu(d)(\lfloor\frac H{Kd} \rfloor-\lfloor\frac {L-1}{Kd}+1 \rfloor+1)^n=\sum_{d=1}^{\lfloor\frac HK \rfloor }\mu(d)(\lfloor\frac H{Kd} \rfloor-\lfloor\frac {L-1}{Kd} \rfloor)^n
\]
\[用数论分块做，但是\mu的范围很大，所以用杜教筛来处理\mu的前缀和
\]
\[f(1)S(n) = \sum_{i=1}^nF(i)-\sum_{d=2}^nf(d)S(\lfloor\frac nd \rfloor)
\]
\[其中F=\epsilon,g = \mu,f = 1
\]
view code
//#pragma GCC optimize("O3")
//#pragma comment(linker, "/STACK:1024000000,1024000000")
#include<bits/stdc++.h>
using namespace std;
function<void(void)> ____ = [](){ios_base::sync_with_stdio(false); cin.tie(0); cout.tie(0);};
const int MAXN = 2e6+7;
typedef long long int LL;
const LL MOD = 1e9+7;
int N,K,L,H,mu[MAXN];
vector<int> prime;
bool pm[MAXN];
LL S[MAXN];
unordered_map<LL,LL> MS;
void sieve(){
mu[1] = 1;
fill(begin(pm),end(pm),true);
for(int i = 2; i < MAXN; i++){
if(pm[i]) prime.push_back(i), mu[i] = -1;
for(int j = 0; i * prime[j] < MAXN; j++){
pm[i*prime[j]] = false;
if(i%prime[j]==0){
mu[i*prime[j]] = 0;
break;
}
mu[i*prime[j]] = -mu[i];
}
}
for(int i = 1; i < MAXN; i++) S[i] = S[i-1] + mu[i];
}
LL ksm(LL a, LL b){
LL ret = 1;
while(b){
if(b&1) ret = ret * a % MOD;
b >>= 1;
a = a * a % MOD;
}
return ret;
}
LL calS(LL n){
if(n<MAXN) return S[n];
if(MS.count(n)) return MS[n];
LL ret = 1;
for(LL l = 2, r; l <= n; l = r + 1){
r = n / (n / l);
ret = (ret - (r-l+1) * calS(n/l) % MOD + MOD) % MOD;
}
return MS[n] = ret;
}
void solve(){
LL ret = 0;
LL ll = (L-1)/K, rr = H/K;
for(LL l = 1, r; l <= rr; l = r + 1){
if(ll/l==0) r = rr / (rr / l);
else r = min(rr / (rr / l),ll / (ll / l));
ret = (ret + (calS(r) - calS(l-1) + MOD) * ksm(rr/l-ll/l,N) % MOD) % MOD;
}
cout << ret << endl;
}
int main(){
scanf("%d %d %d %d",&N,&K,&L,&H);
sieve();
solve();
return 0;
}
- LuoguP2522 [HAOI2011]Problem b🔗
题解：
\[考虑容斥之后，套用第二题的做法即可
\]
\[设F(n,m)表示1\le x\le n,1\le y\le m的情况下gcd(n,m)=k的情况数
\]
\[Ans = F(b,d) - F(a-1,d) - F(b,c-1) + F(a-1,c-1)
\]
view code
//#pragma GCC optimize("O3")
//#pragma comment(linker, "/STACK:1024000000,1024000000")
#include<bits/stdc++.h>
using namespace std;
function<void(void)> ____ = [](){ios_base::sync_with_stdio(false); cin.tie(0); cout.tie(0);};
const int MAXN = 2e6+7;
int mu[MAXN],S[MAXN];
vector<int> prime;
bool pm[MAXN];
void sieve(){
mu[1] = 1;
fill(begin(pm),end(pm),true);
for(int i = 2; i < MAXN; i++){
if(pm[i]) prime.push_back(i), mu[i] = -1;
for(int j = 0; i * prime[j] < MAXN; j++){
pm[i*prime[j]] = false;
if(i%prime[j]==0){
mu[i*prime[j]] = 0;
break;
}
mu[i*prime[j]] = -mu[i];
}
}
for(int i = 1; i < MAXN; i++) S[i] = S[i-1] + mu[i];
}
int rua(int n, int m, int k){
if(n>m) swap(n,m);
n /= k, m /= k;
int ret = 0;
for(int l = 1, r; l <= n; l = r + 1){
r = min(n / (n / l),m / (m / l));
ret = ret + (S[r] - S[l-1]) * (n / l) * (m / l);
}
return ret;
}
void solve(){
int a, b, c, d, k;
cin >> a >> b >> c >> d >> k;
cout << rua(b,d,k) - rua(b,c-1,k) - rua(a-1,d,k) + rua(a-1,c-1,k) << endl;
}
int main(){
sieve();
____();
int T;
for(cin >> T; T; T--) solve();
return 0;
}
- LuoguP3327 [SDOI2015]约数个数和🔗
题解：
\[d(n)为除数函数,计算\sum_{i=1}^n\sum_{j=1}^md(i\cdot j)
\]
\[首先有个公式：d(ij)=\sum_{x\mid i}\sum_{y\mid j}[gcd(x,y)=1]
\]
证明：
\[d(ij)=\sum_{x\mid i}\sum_{y\mid j}[gcd(x,y)=1]
\]
考虑把\(i\cdot j\)的所有因子一一映射，考虑其中一个因子中存在\(p^x\)，而\(i\)中存在\(p^a\)，\(j\)中存在\(p^b\)
\(\begin{cases}x\le a时只取i中的p^x \\ x>a时取i中的p^a再取j中的p^{x-a} \end{cases}\)
比如\(a=2,b=3\)
\(\begin{cases}(p^0,p^0) & x=0 \\ (p^1,p^0)&x=1\\ (p^2,p^0) & x=2 \\ (p^0,p^1) & x=3\\ (p^0,p^2) & x=4 \\ (p^0,p^3) & x=5 \end{cases}\)
对于所有质因子都这样处理
\[下面假设n\le m
\]
\[\sum_{i=1}^n\sum_{j=1}^md(i\cdot j) = \sum_{i=1}^n\sum_{j=1}^m\sum_{x\mid i}\sum_{y\mid j}[gcd(x,y)=1]
\]
\[=\sum_{i=1}^n\sum_{j=1}^m\sum_{x\mid i}\sum_{y\mid j}\sum_{d|gcd(x,y)}\mu(d)
\]
\[枚举x和y
\]
\[\Rightarrow \sum_{x=1}^n\sum_{y=1}^m\sum_{x\mid i}^n\sum_{y\mid j}^m\sum_{d\mid gcd(x,y)}\mu(d)=\sum_{x=1}^n\sum_{y=1}^m\lfloor \frac nx\rfloor\lfloor\frac my \rfloor \sum_{d\mid gcd(x,y)}\mu(d)
\]
\[枚举d
\]
\[\Rightarrow \sum_{d=1}^n \mu(d)\sum_{d\mid x}^n\sum_{d\mid y}^m\lfloor \frac nx\rfloor\lfloor\frac my \rfloor
\]
\[令a = \frac xd,b=\frac yd
\]
\[\Rightarrow \sum_{d=1}^n\mu(d)\sum_{a}^{\lfloor\frac nd \rfloor}\sum_{b}^{\lfloor\frac md \rfloor}\lfloor \frac n{ad}\rfloor\lfloor\frac m{bd} \rfloor=\sum_{d=1}^n\mu(d)(\sum_{a}^{\lfloor\frac nd \rfloor}\lfloor \frac n{ad}\rfloor)(\sum_{b}^{\lfloor\frac md \rfloor}\lfloor\frac m{bd} \rfloor)
\]
\[令\sum_{i=1}^n \lfloor\frac{n}{i} \rfloor=P(n)
\]
\[\Rightarrow \sum_{d=1}^n\mu(d)(\sum_{a}^{\lfloor\frac nd \rfloor}\lfloor \frac n{ad}\rfloor)(\sum_{b}^{\lfloor\frac md \rfloor}\lfloor\frac m{bd} \rfloor)=\sum_{d=1}^n\mu(d)P(\lfloor\frac nd \rfloor)P(\lfloor\frac md \rfloor)
\]
\[预处理出来P就可以数论分块解决了
\]
\[P也可以数论分块来计算
\]
view code
//#pragma GCC optimize("O3")
//#pragma comment(linker, "/STACK:1024000000,1024000000")
#include<bits/stdc++.h>
using namespace std;
function<void(void)> ____ = [](){ios_base::sync_with_stdio(false); cin.tie(0); cout.tie(0);};
const int MAXN = 5e4+7;
typedef long long int LL;
int mu[MAXN],S[MAXN];
LL pre[MAXN];
vector<int> prime;
bool pm[MAXN];
void preprocess(){
mu[1] = 1;
fill(begin(pm),end(pm),true);
for(int i = 2; i < MAXN; i++){
if(pm[i]) prime.push_back(i), mu[i] = -1;
for(int j = 0; i * prime[j] < MAXN; j++){
pm[i*prime[j]] = false;
if(i%prime[j]==0){
mu[i*prime[j]] = 0;
break;
}
mu[i*prime[j]] = -mu[i];
}
}
for(int i = 1; i < MAXN; i++) S[i] = S[i-1] + mu[i];
for(int i = 1; i < MAXN; i++){
for(int l = 1, r; l <= i; l = r + 1){
r = i / (i / l);
pre[i] += (i / l) * (r - l + 1);
}
}
}
void solve(){
static LL n, m;
scanf("%d %d",&n,&m);
if(n>m) swap(n,m);
LL ret = 0;
for(int l = 1, r; l <= n; l = r + 1){
r = min(n / (n / l), m / (m / l));
ret += (S[r] - S[l-1]) * (pre[n/l]*pre[m/l]);
}
printf("%lld\n",ret);
}
int main(){
preprocess();
int T; for(scanf("%d",&T); T; T--) solve();
return 0;
}
- LuoguP3911 最小公倍数之和🔗
题解：
\[A_i\le 5\times 10^4 所以可以先把每个数出现的次数记录下来，记为c数组
\]
\[转化之后就是计算：\sum_{i=1}^n\sum_{j=1}^nlcm(i,j)\cdot c_i\cdot c_j
\]
\[\Rightarrow \sum_{i=1}^n\sum_{j=1}^n\frac{i\cdot j\cdot c_i\cdot c_j}{gcd(i,j)}
\]
\[枚举g=gcd(i,j)
\]
\[\Rightarrow \sum_{g=1}^n\sum_{g|i}\sum_{g|j}[gcd(\frac ig,\frac jg)=1]\frac{i\cdot j\cdot c_i\cdot c_j}{g}
\]
\[令x=\frac ig,y = \frac jg
\]
\[\Rightarrow \sum_{g=1}^n\sum_{x=1}^{\lfloor \frac ng \rfloor}\sum_{y=1}^{\lfloor \frac ng\rfloor}[gcd(x,y)=1]g\cdot x\cdot y\cdot c_{xg}\cdot c_{yg}
\]
\[=\sum_{g=1}^n\sum_{x=1}^{\lfloor \frac ng \rfloor}\sum_{y=1}^{\lfloor \frac ng\rfloor}\sum_{d\mid gcd(x,y)}\mu(d)\cdot g\cdot x\cdot y\cdot c_{xg}\cdot c_{yg}
\]
\[枚举d
\]
\[\Rightarrow \sum_{g=1}^n\sum_{d=1}^{\lfloor \frac ng\rfloor}\mu(d)\sum_{d\mid x}^{\lfloor \frac ng\rfloor }\sum_{d\mid y}^{\lfloor \frac ng\rfloor}g\cdot x\cdot y\cdot c_{xg}\cdot c_{yg}
\]
\[令a = \frac xd,b = \frac yd
\]
\[\Rightarrow \sum_{g=1}^n\sum_{d=1}^{\lfloor \frac ng\rfloor}\mu(d)\sum_{a=1}^{\lfloor \frac{n}{gd}\rfloor}\sum_{b=1}^{\lfloor \frac{n}{gd}\rfloor}g\cdot a\cdot b\cdot d^2\cdot c_{adg}\cdot c_{bdg}
\]
\[令T=gd，枚举T
\]
\[\Rightarrow \sum_{T=1}^nT\sum_{d\mid T}d\cdot\mu(d)\sum_{a=1}^{\lfloor \frac nT\rfloor}\sum_{b=1}^{\lfloor \frac nT\rfloor}a\cdot b\cdot c_{aT}\cdot c_{bT}
\]
\[=\sum_{T=1}^nT(\sum_{i=1}^{\lfloor \frac nT\rfloor}i\cdot c_{iT})^2\sum_{d|T}d\cdot \mu(d)
\]
\[令S(T)=\sum_{d\mid T}d\cdot \mu(d)，P(T)=\sum_{i=1}^{\lfloor \frac nT\rfloor}i\cdot c_{iT}
\]
\[\Rightarrow \sum_{T=1}^nT\cdot P^2(T)\cdot S(T)
\]
\[其中S和P都能n\log n预处理出来
\]
view code
//#pragma GCC optimize("O3")
//#pragma comment(linker, "/STACK:1024000000,1024000000")
#include<bits/stdc++.h>
using namespace std;
function<void(void)> ____ = [](){ios_base::sync_with_stdio(false); cin.tie(0); cout.tie(0);};
const int MAXN = 5e4+7;
typedef long long int LL;
int n,c[MAXN],mu[MAXN];
LL S[MAXN],P[MAXN];
vector<int> prime;
bool pm[MAXN];
void preprocess(){
fill(begin(pm),end(pm),true);
mu[1] = 1;
for(int i = 2; i < MAXN; i++){
if(pm[i]) prime.push_back(i), mu[i] = -1;
for(int j = 0; i * prime[j] < MAXN; j++){
pm[i*prime[j]] = false;
if(i%prime[j]==0){
mu[i*prime[j]] = 0;
break;
}
mu[i*prime[j]] = -mu[i];
}
}
for(int i = 1; i <= n; i++)
for(int j = 1; i * j <= n; j++)
S[i*j] += i * mu[i];
for(int T = 1; T <= n; T++)
for(int i = 1; i <= n / T; i++)
P[T] += i * c[i * T];
}
int main(){
scanf("%d",&n);
for(int i = 1; i <= n; i++){
int x; scanf("%d",&x);
c[x] += 1;
}
n = 50000;
preprocess();
LL ret = 0;
for(int i = 1; i <= n; i++) ret += i * S[i] * P[i] * P[i];
cout << ret << endl;
return 0;
}
- LuoguP1829 [国家集训队]Crash的数字表格 / JZPTAB🔗
题解：
\[计算\sum_{i=1}^n\sum_{j=1}^mlcm(i,j)
\]
\[开始推柿子吧，假设n\le m
\]
\[\sum_{i=1}^n\sum_{j=1}^mlcm(i,j) = \sum_{i=1}^n\sum_{j=1}^m\frac{i\cdot j}{gcd(i,j)}
\]
\[枚举g=gcd(i,j)
\]
\[\Rightarrow \sum_{g=1}^n\sum_{g\mid i}\sum_{g\mid j}[gcd(\lfloor \frac ig\rfloor,\lfloor \frac jg \rfloor)=1]\frac{i\cdot j}{g}
\]
\[令x = \frac ig, y = \frac jg
\]
\[\Rightarrow \sum_{g=1}^n\sum_{x=1}^{\lfloor \frac ng\rfloor }\sum_{y=1}^{\lfloor \frac mg\rfloor }[gcd(x,y)=1]x\cdot y\cdot g = \sum_{g=1}^n\sum_{x=1}^{\lfloor \frac ng\rfloor }\sum_{y=1}^{\lfloor \frac mg\rfloor }\sum_{d|gcd(x,y)}\mu(d)\cdot x\cdot y\cdot g
\]
\[枚举d
\]
\[\Rightarrow \sum_{g=1}^n\sum_{x=1}^{\lfloor \frac ng\rfloor }\sum_{y=1}^{\lfloor \frac mg\rfloor }\sum_{d|gcd(x,y)}\mu(d)\cdot x\cdot y\cdot g = \sum_{g=1}^ng\sum_{d=1}^{\lfloor \frac ng\rfloor}\mu(d)\sum_{d\mid x}^{\lfloor \frac ng\rfloor }\sum_{d\mid y}^{\lfloor \frac mg\rfloor }x\cdot y
\]
\[令a = \frac xd, b = \frac yd
\]
\[\Rightarrow \sum_{g=1}^ng\sum_{d=1}^{\lfloor \frac ng\rfloor}\mu(d)\sum_{d\mid x}^{\lfloor \frac ng\rfloor }\sum_{d\mid y}^{\lfloor \frac mg\rfloor }x\cdot y = \sum_{g=1}^ng\sum_{d=1}^{\lfloor \frac ng\rfloor}\mu(d)\sum_{a=1}^{\lfloor\frac n{gd}\rfloor}\sum_{b=1}^{\lfloor\frac m{gd}\rfloor}d^2\cdot a\cdot b
\]
\[=\sum_{g=1}^ng\sum_{d=1}^{\lfloor \frac ng\rfloor}d^2\cdot \mu(d)\sum_{a=1}^{\lfloor\frac n{gd}\rfloor}\sum_{b=1}^{\lfloor\frac m{gd}\rfloor}a\cdot b=\sum_{g=1}^ng\sum_{d=1}^{\lfloor \frac ng\rfloor}[d^2\cdot \mu(d)(\sum_{a=1}^{\lfloor\frac n{gd}\rfloor}a)(\sum_{b=1}^{\lfloor\frac m{gd}\rfloor}b)]
\]
\[到这里，首先对最外层进行数论分块，现在只考虑里面的式子，令N=\lfloor\frac ng\rfloor,M=\lfloor \frac mg\rfloor
\]
\[\sum_{d=1}^{\lfloor \frac ng\rfloor}[d^2\cdot \mu(d)(\sum_{a=1}^{\lfloor\frac n{gd}\rfloor}a)(\sum_{b=1}^{\lfloor\frac m{gd}\rfloor}b)]=\sum_{d=1}^{N}d^2\cdot \mu(d)(\sum_{a=1}^{\lfloor\frac Nd\rfloor}a)(\sum_{b=1}^{\lfloor\frac Md\rfloor}b)
\]
\[令sum(n) = \sum_{i=1}^ni
\]
\[\sum_{d=1}^{N}d^2\cdot \mu(d)(\sum_{a=1}^{\lfloor\frac Nd\rfloor}a)(\sum_{b=1}^{\lfloor\frac Md\rfloor}b)=\sum_{d=1}^{N}d^2\cdot \mu(d)sum(\lfloor\frac Nd\rfloor)sum(\lfloor\frac Md\rfloor)
\]
\[发现这个式子也是可以数论分块做的，所以总的就是再一个分块里再做一次分块，复杂度差不多就是O(n+m)
\]
view code
//#pragma GCC optimize("O3")
//#pragma comment(linker, "/STACK:1024000000,1024000000")
#include<bits/stdc++.h>
using namespace std;
function<void(void)> ____ = [](){ios_base::sync_with_stdio(false); cin.tie(0); cout.tie(0);};
const int MAXN = 1e7+7;
typedef long long int LL;
const LL MOD = 20101009;
const LL inv2 = 10050505;
vector<int> prime;
bool pm[MAXN];
int mu[MAXN];
LL S[MAXN],pre[MAXN];
void preprocess(){
fill(begin(pm),end(pm),true);
mu[1] = 1;
for(int i = 2; i < MAXN; i++){
if(pm[i]) mu[i] = -1, prime.push_back(i);
for(int j = 0; i * prime[j] < MAXN; j++){
pm[i*prime[j]] = false;
if(i%prime[j]==0){
mu[i*prime[j]] = 0;
break;
}
mu[i*prime[j]] = -mu[i];
}
}
for(LL i = 1; i < MAXN; i++){
S[i] = (S[i-1] + i * i * mu[i]) % MOD;
pre[i] = (pre[i-1] + i) % MOD;
}
}
LL sum(LL x){ return x * (x + 1) % MOD * inv2 % MOD; }
void solve(){
int n, m;
scanf("%d %d",&n,&m);
if(n>m) swap(n,m);
LL ret = 0;
for(int l = 1, r; l <= n; l = r + 1){
r = min(n / (n / l), m / (m / l));
int N = n / l, M = m / l;
LL tp = 0;
for(int L = 1, R; L <= N; L = R + 1){
R = min(N / (N / L), M / (M / L));
tp = (tp + (S[R] - S[L-1] + MOD) * sum(N / L) % MOD * sum(M / L) % MOD) % MOD;
}
ret = (ret + (pre[r] - pre[l-1] + MOD) * tp) % MOD;
}
cout << ret << endl;
}
int main(){
preprocess();
solve();
return 0;
}
- CF435F Cowslip Collections🔗
题解：
\[要求计算的是所有k元组的gcd的和
\]
\[那么假设f(n)表示gcd为n的k元组的数量，那么答案就是\sum_{i=1}^Ni\cdot f(i)
\]
\[假设F(n)表示gcd为n的倍数的k元组的数量，可以得到f(n) = \sum_{n|d}\mu(\frac dn)F(d)
\]
\[设cnt_x表示x的倍数出现的次数，那么F(n)=C(cnt_n,k)
\]
\[那么答案可以改写为：\sum_{i=1}^Ni\sum_{i|d}\mu(\frac di)C(cnt_d,k)
\]
\[由于是多组修改询问，考虑每次增加一个数对答案的贡献，假设这个数为x，容易每次修改会对\sqrt{x}个cnt进行修改
\]
\[考虑每次改变一个cnt_x对答案的贡献是(C(cnt_x,k)-C(cnt_x-1,k))\cdot \sum_{d|x}d\cdot\mu(\frac xd)
\]
\[发现后式可以n\log n预处理，所以每次修改可以做到O(\sqrt{x})
\]
\[总复杂度为O(n\sqrt{n})
\]
view code
//#pragma GCC optimize("O3")
//#pragma comment(linker, "/STACK:1024000000,1024000000")
#include<bits/stdc++.h>
using namespace std;
function<void(void)> ____ = [](){ios_base::sync_with_stdio(false); cin.tie(0); cout.tie(0);};
typedef long long int LL;
const int MAXN = 1e6+7;
const LL MOD = 1e9+7;
int n,k,q,cnt[MAXN],mu[MAXN];
LL S[MAXN],f[MAXN],rf[MAXN];
bool pm[MAXN];
vector<int> prime;
LL ksm(LL a, LL b){
LL ret = 1;
while(b){
if(b&1) ret = ret * a % MOD;
b >>= 1;
a = a * a % MOD;
}
return ret;
}
void sieve(){
mu[1] = 1;
fill(begin(pm),end(pm),true);
for(int i = 2; i < MAXN; i++){
if(pm[i]) prime.push_back(i), mu[i] = -1;
for(int j = 0; i * prime[j] < MAXN; j++){
pm[i*prime[j]] = false;
if(i%prime[j]==0){
mu[i*prime[j]] = 0;
break;
}
mu[i*prime[j]] = -mu[i];
}
}
for(int i = 1; i < MAXN; i++)
for(int j = i; j < MAXN; j += i)
S[j] += i * mu[j/i];
}
vector<int> getFact(int x){
vector<int> vec;
for(int i = 1; i * i <= x; i++){
if(x%i!=0) continue;
vec.push_back(i);
if(i!=x/i) vec.push_back(x/i);
}
return vec;
}
LL C(int n, int m){ return m > n ? 0 : f[n] * rf[m] % MOD * rf[n-m] % MOD; }
int main(){
sieve();
scanf("%d %d %d",&n,&k,&q);
f[0] = 1;
for(int i = 1; i < MAXN; i++) f[i] = f[i-1] * i % MOD;
rf[MAXN-1] = ksm(f[MAXN-1],MOD-2);
for(int i = MAXN - 2; i >= 0; i--) rf[i] = rf[i+1] * (i+1) % MOD;
LL ret = 0;
for(int i = 1; i <= n + q; i++){
int x; scanf("%d",&x);
auto fact = getFact(x);
for(int m : fact){
cnt[m]++;
ret = (ret + (C(cnt[m],k) - C(cnt[m]-1,k) + MOD) * S[m]) % MOD;
}
if(i>n) printf("%I64d\n",ret);
}
return 0;
}
- CF1043F Make It One🔗
题解：
\[给出n个数，问最少选几个数可以使得这些数的gcd等于1，不存在输出-1
\]
\[首先如果所有数的gcd不等于1，那就直接输出-1，否则必然存在解
\]
\[我们设最少选k个数可以满足条件，可以发现k是可以二分的
\]
\[现在问题转化为判断选定k个数能否得到gcd=1
\]
\[令F(n)表示gcd为n的倍数的k元组的方案数，f(n)表示gcd为n的k元组的方案数
\]
\[有f(n) = \sum_{n\mid d}\mu(\frac dn)F(d)
\]
\[我们要判断的就是f(1)=\sum_{i=1}^N\mu(i)F(i)是否不为0
\]
\[设cnt_x表示x的倍数的数的个数，F(x) = C(cnt_x,k)
\]
\[cnt数组可以O(n\log n)计算出来，二分之后每次可以O(n)计算f(1)
\]
\[总复杂度O(n\log n)
\]
view code
//#pragma GCC optimize("O3")
//#pragma comment(linker, "/STACK:1024000000,1024000000")
#include<bits/stdc++.h>
using namespace std;
function<void(void)> ____ = [](){ios_base::sync_with_stdio(false); cin.tie(0); cout.tie(0);};
typedef long long int LL;
const int MAXN = 3e5+7;
const LL MOD = 998244353;
LL f[MAXN],rf[MAXN];
int n,cnt[MAXN],mu[MAXN];
vector<int> prime;
bool pm[MAXN];
LL ksm(LL a, LL b){
LL ret = 1;
while(b){
if(b&1) ret = ret * a % MOD;
a = a * a % MOD;
b >>= 1;
}
return ret;
}
void sieve(){
mu[1] = 1;
fill(begin(pm),end(pm),true);
for(int i = 2; i < MAXN; i++){
if(pm[i]) prime.push_back(i), mu[i] = -1;
for(int j = 0; i * prime[j] < MAXN; j++){
pm[i*prime[j]] = false;
if(i%prime[j]==0){
mu[i*prime[j]] = 0;
break;
}
mu[i*prime[j]] = -mu[i];
}
}
}
LL C(LL n, LL m){ return m>n?0:f[n]*rf[m]%MOD*rf[n-m]%MOD; }
bool check(int m){
LL tot = 0;
for(int i = 1; i < MAXN; i++) tot = (tot + mu[i] * C(cnt[i],m)) % MOD;
return tot!=0;
}
int main(){
sieve();
f[0] = 1;
for(int i = 1; i < MAXN; i++) f[i] = f[i-1] * i % MOD;
rf[MAXN-1] = ksm(f[MAXN-1],MOD-2);
for(int i = MAXN - 2; i >= 0; i--) rf[i] = rf[i+1] * (i+1) % MOD;
scanf("%d",&n);
int g = 0;
for(int i = 1; i <= n; i++){
int x; scanf("%d",&x);
g = __gcd(x,g);
for(int j = 1; j * j <= x; j++){
if(x%j!=0) continue;
cnt[j]++;
if(j!=x/j) cnt[x/j]++;
}
}
if(g!=1){
cout << -1 << endl;
return 0;
}
int l = 1, r = n;
while(l<=r){
int mid = (l+r) >> 1;
if(check(mid)) r = mid - 1;
else l = mid + 1;
}
cout << l << endl;
return 0;
}
- CF235E. Number Challenge🔗
题解：
\[首先有个公式：d(ijk)=\sum_{x\mid i}\sum_{y\mid j}\sum_{z\mid k}[(x,y)=1][(x,z)=1][(y,z)=1]
\]
\[所以就是计算\sum_{i=1}^a\sum_{j=1}^b\sum_{k=1}^c\sum_{x\mid i}\sum_{y\mid j}\sum_{z\mid k}[(x,y)=1][(x,z)=1][(y,z)=1]
\]
\[枚举x,y,z
\]
\[\Rightarrow \sum_{x=1}^a\sum_{y=1}^b\sum_{z=1}^c\sum_{x\mid i}\sum_{y\mid j}\sum_{z\mid k}[(x,y)=1][(x,z)=1][(y,z)=1]
\]
\[=\sum_{x=1}^a\sum_{y=1}^b\sum_{z=1}^c\lfloor \frac ax \rfloor\lfloor \frac by\rfloor \lfloor \frac cz\rfloor [(x,y)=1][(x,z)=1][(y,z)=1]
\]
\[把[(x,y)=1]用\sum_{d\mid (x,y)}\mu(d)代替
\]
\[\Rightarrow \sum_{z=1}^c\sum_{y=1}^b\sum_{z=1}^c\lfloor \frac ax \rfloor\lfloor \frac by\rfloor \lfloor \frac cz\rfloor \sum_{d\mid (x,y)}\mu(d)[(x,z)=1][(y,z)=1]
\]
\[枚举d
\]
\[\Rightarrow \sum_{z=1}^c\lfloor \frac cz\rfloor \sum_{d=1}^{a}\mu(d) \sum_{d\mid x}\lfloor \frac ax \rfloor\sum_{d\mid y}\lfloor \frac by\rfloor [(x,z)=1][(y,z)=1]
\]
\[令i = \frac xd,j = \frac yd
\]
\[\Rightarrow \sum_{z=1}^c\lfloor \frac cz\rfloor \sum_{d=1}^{a}\mu(d) \sum_{i=1}^{\lfloor \frac ad \rfloor}\lfloor \frac a{id} \rfloor[(id,z)=1]\sum_{j=1}^{\lfloor \frac bd \rfloor}\lfloor \frac b{jd} \rfloor [(jd,z)=1]
\]
\[令f(n,m) = \sum_{i=1}^n\lfloor \frac ni \rfloor [(i,z)=1]
\]
\[\Rightarrow \sum_{z=1}^c\lfloor \frac cz\rfloor \sum_{d=1}^{a}\mu(d) [(d,z)=1]f(\lfloor \frac ad\rfloor,z)f(\lfloor \frac bd\rfloor,z)
\]
\[这个式子可以在n^2\log n的时间内算出来
\]
view code
//#pragma GCC optimize("O3")
//#pragma comment(linker, "/STACK:1024000000,1024000000")
#include<bits/stdc++.h>
using namespace std;
function<void(void)> ____ = [](){ios_base::sync_with_stdio(false); cin.tie(0); cout.tie(0);};
typedef long long int LL;
const int MAXN = 2222;
const LL MOD = 1<<30;
int a,b,c,mu[MAXN],g[MAXN][MAXN];
bool pm[MAXN];
vector<int> prime;
void sieve(){
fill(begin(pm),end(pm),true);
mu[1] = 1;
for(int i = 2; i < MAXN; i++){
if(pm[i]) prime.push_back(i), mu[i] = -1;
for(int j = 0; i * prime[j] < MAXN; j++){
pm[i*prime[j]] = false;
if(i%prime[j]==0){
mu[i*prime[j]] = 0;
break;
}
mu[i*prime[j]] = -mu[i];
}
}
}
int gcd(int a, int b){
if(g[a][b]) return g[a][b];
return g[a][b] = g[b][a] = __gcd(a,b);
}
LL f(int n, int m){
LL ret = 0;
for(int i = 1; i <= n; i++) if(gcd(i,m)==1) ret += n / i;
return ret;
}
void solve(){
scanf("%d %d %d",&a,&b,&c);
if(a>b) swap(a,b);
LL ret = 0;
for(int m = 1; m <= c; m++){
for(int d = 1; d <= a; d++){
if(gcd(d,m)!=1) continue;
ret = (ret + c / m * (mu[d] * f(a/d,m) % MOD * f(b/d,m) % MOD) % MOD) % MOD;
}
}
cout << (ret + MOD) % MOD << endl;
}
int main(){
sieve();
solve();
return 0;
}
- CF1139D. Steps to One🔗
题解：
\[f[x]表示当前gcd为x的情况下期望还要加多少个数能使得gcd变成1
\]
\[那么f[x] = 1+\frac{\sum_{i=1}^nf[gcd(i,x)]}{n}
\]
\[n\cdot f[x] = n + \sum_{i=1}^nf[gcd(i,x)]
\]
\[现在要计算\sum_{i=1}^nf[gcd(i,x)],考虑枚举gcd
\]
\[\sum_{d\mid x}f[d]\sum_{i=1}^n[gcd(i,x)=d]=\sum_{d\mid x}f[d]\sum_{i=1}^{\lfloor \frac nd\rfloor}[gcd(i,\frac xd)=1]
\]
\[=\sum_{d\mid x}f[d]\sum_{i=1}^{\lfloor \frac nd\rfloor}\sum_{p\mid i,p\mid \frac xd}\mu(p)=\sum_{d\mid x}f[d]\sum_{p\mid \frac xd}\mu(p)\lfloor \frac n{pd}\rfloor
\]
\[代回原式把f[x]消掉可以得到：
\]
\[(n-\lfloor \frac nx\rfloor )f[x] = n + \sum_{d\mid x,d\ne x}f(d)\sum_{p\mid \frac xd}\mu(p)\lfloor \frac n{pd}\rfloor
\]
\[解一下方程就好了再全部统计一遍就好了
\]
view code
//#pragma GCC optimize("O3")
//#pragma comment(linker, "/STACK:1024000000,1024000000")
#include<bits/stdc++.h>
using namespace std;
function<void(void)> ____ = [](){ios_base::sync_with_stdio(false); cin.tie(0); cout.tie(0);};
typedef long long int LL;
const LL MOD = 1e9+7;
const int MAXN = 1e5+7;
LL ksm(LL a, LL b){
LL ret = 1;
while(b){
if(b&1) ret = ret * a % MOD;
b >>= 1;
a = a * a % MOD;
}
return ret;
}
vector<int> prime,fact[MAXN];
bool pm[MAXN];
int mu[MAXN];
LL f[MAXN];
void sieve(){
fill(begin(pm),end(pm),true);
mu[1] = 1;
for(int i = 2; i < MAXN; i++){
if(pm[i]) prime.push_back(i), mu[i] = -1;
for(int j = 0; i * prime[j] < MAXN; j++){
pm[i*prime[j]] = false;
if(i%prime[j]==0){
mu[i*prime[j]] = 0;
break;
}
mu[i*prime[j]] = -mu[i];
}
}
for(int i = 1; i < MAXN; i++)
for(int j = i; j < MAXN; j += i)
fact[j].push_back(i);
}
void solve(){
int n;
scanf("%d",&n);
f[1] = 0;
for(int i = 2; i <= n; i++){
f[i] = n;
for(int d : fact[i]){
if(d==i) continue;
for(int p : fact[i / d]) f[i] = (f[i] + f[d] * mu[p] * (n/p/d)) % MOD;
}
f[i] = f[i] * ksm(n-n/i,MOD-2) % MOD;
}
LL ans = 0, invn = ksm(n,MOD-2);
for(int i = 1; i <= n; i++) ans = (ans + f[i]) % MOD;
ans = ans * invn % MOD;
ans = (ans + 1 + MOD) % MOD;
cout << ans << endl;
}
int main(){
sieve();
solve();
return 0;
}
- CF439E. Devu and Birthday Celebration🔗
题解：
$$f[x]表示当前gcd为x的情况下期望还要加多少个数能使得gcd变成1$$\[F(x)表示分组之后的gcd为x的倍数的方案数
\]
\[f(x)表示分组之后的gcd为x的方案数
\]
\[f(x)=\sum_{x|d}\mu(\frac xd)F(d)
\]
\[答案就是f(1) = \sum_{i=1}^n\mu(i)F[i]
\]
\[而F(x)只有在x\mid n且\frac nx >m的情况下才有值，值为C(\frac nx-1,m-1)
\]
\[对于每次询问计算所有因子的F值，然后计算f(1)即可
\]
view code
//#pragma GCC optimize("O3")
//#pragma comment(linker, "/STACK:1024000000,1024000000")
#include<bits/stdc++.h>
using namespace std;
function<void(void)> ____ = [](){ios_base::sync_with_stdio(false); cin.tie(0); cout.tie(0);};
typedef long long int LL;
const int MAXN = 1e5+7;
const LL MOD = 1e9+7;
int n,m,mu[MAXN];
LL f[MAXN],rf[MAXN],g[MAXN];
vector<int> prime;
bool pm[MAXN];
LL ksm(LL a, LL b){
LL ret = 1;
while(b){
if(b&1) ret = ret * a % MOD;
a = a * a % MOD;
b >>= 1;
}
return ret;
}
void sieve(){
fill(begin(pm),end(pm),true);
mu[1] = 1;
for(int i = 2; i < MAXN; i++){
if(pm[i]) prime.push_back(i), mu[i] = -1;
for(int j = 0; i * prime[j] < MAXN; j++){
pm[i*prime[j]] = false;
if(i%prime[j]==0){
mu[i*prime[j]] = 0;
break;
}
mu[i*prime[j]] = -mu[i];
}
}
f[0] = 1;
for(int i = 1; i < MAXN; i++) f[i] = f[i-1] * i % MOD;
rf[MAXN-1] = ksm(f[MAXN-1],MOD-2);
for(int i = MAXN - 2; i >= 0; i--) rf[i] = rf[i+1] * (i+1) % MOD;
}
LL C(int n, int m){ return m>n?0:f[n]*rf[m]%MOD*rf[n-m]%MOD; }
void solve(){
scanf("%d %d",&n,&m);
vector<int> fact;
for(int i = 1; i * i <= n; i++){
if(n%i!=0) continue;
if(n/i>=m) fact.push_back(i);
if(i==n/i) continue;
if(i>=m) fact.push_back(n/i);
}
LL ret = 0;
for(int x : fact) ret = (ret + mu[x] * C(n/x-1,m-1)) % MOD;
printf("%I64d\n",(ret+MOD) % MOD);
}
int main(){
sieve();
int T;
for(scanf("%d",&T); T; T--) solve();
return 0;
}
