莫比乌斯反演学习笔记

莫比乌斯反演

一、前置技能

积性函数
定义$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;
}

posted @ 2020-07-04 23:56 _kiko 阅读(317) 评论(1) 编辑收藏举报

会员力量，点亮园子希望

刷新页面返回顶部

_kiko

莫比乌斯反演 学习笔记

一、前置技能

二、莫比乌斯函数

三、莫比乌斯反演

四、例题：

公告

莫比乌斯反演学习笔记