我爱数论
jycjycjyc
·
·
算法·理论
这篇文章的推导对有着一定的莫反印象，知道莫反的几个基本公式但不知道莫反的解题套路的人很友好，因为相邻两步之间的跨度并不远，认真阅读可以帮助加强对莫反解题套路的熟练度。
有一些题目的最后一步用到了杜教筛，没学过的读者可以自行学习。
[清华集训 2012] 模积和
题意
给定正整数 n,m，求
\sum_{i=1}^n \sum_{j=1}^m (n \bmod i)(m \bmod j),i \neq j
对 19940417 取模的值．
解析
不妨设 n \leq m，并将原式做如下变换：
\begin{aligned}
\sum_{i=1}^n \sum_{j=1}^m (n \bmod i)(m \bmod j),i \neq j
&= \sum_{i=1}^n \sum_{j=1}^m (n \bmod i)(m \bmod j) - \sum_{i=1}^n (n \bmod i)(m \bmod i) \\
&= \sum_{i=1}^n \sum_{j=1}^m (n- \lfloor \frac{n}{i} \rfloor i)(m - \lfloor \frac{m}{j} \rfloor j) -\sum_{i=1}^n (n- \lfloor \frac{n}{i} \rfloor i)(m- \lfloor \frac{m}{i} \rfloor i) \\
&=\sum_{i=1}^n \sum_{j=1}^m (nm-nj \lfloor \frac{m}{j} \rfloor-mi\lfloor \frac{n}{i} \rfloor+i\lfloor \frac{n}{i} \rfloor j \lfloor \frac{m}{j} \rfloor) - \sum_{i=1}^n(nm-ni \lfloor \frac{m}{i} \rfloor - mi \lfloor \frac{n}{i} \rfloor+i^2 \lfloor \frac{n}{i} \rfloor \lfloor \frac{m}{i} \rfloor) \\
&= n^2m^2-n^2\sum_{j=1}^m j \lfloor \frac{m}{j} \rfloor -m^2 \sum_{i=1}^n i \lfloor \frac{n}{i} \rfloor + \sum_{i=1}^n i \lfloor \frac{n}{i} \rfloor \sum_{j=1}^m j \lfloor \frac{m}{j} \rfloor -n^2m+n \sum_{i=1}^n i \lfloor \frac{m}{i} \rfloor +m\sum_{i=1}^n i \lfloor \frac{n}{i} \rfloor -\sum_{i=1}^n i^2 \lfloor \frac{n}{i} \rfloor \lfloor \frac{m}{i} \rfloor
\end{aligned}
设 val_1 = \sum_{i=1}^n i \lfloor \frac{n}{i} \rfloor,val_2 = \sum_{j=1}^m j \lfloor \frac{m}{j} \rfloor,val_3 = \sum_{i=1}^n i \lfloor \frac{m}{i} \rfloor,val_4 = \sum_{i=1}^n i^2 \lfloor \frac{n}{i} \rfloor \lfloor \frac{m}{i} \rfloor，则原式变为
n^2m^2-n^2 \times val_2-m^2 \times val_1+val_1 \times val_2-n^2m+n \times val_3+m \times val_1-val_4
:::success[AC code]
```cpp line-numbers
#include<bits/stdc++.h>
#define int long long
using namespace std;
constexpr int mod = 19940417,inv6 = 3323403;
int n,m;
inline int sum2(const int &l,const int &r){
return (r*(r+1)%mod*(2*r+1)%mod*inv6%mod-(l-1)*l%mod*(2*l-1)%mod*inv6%mod+mod)%mod;
}
signed main(){
ios::sync_with_stdio(false);
cin.tie(0),cout.tie(0);
cin >> n >> m;
if(n>m){
swap(n,m);
}
int val1 = 0,val2 = 0,val3 = 0,val4 = 0;
for(int l = 1,r = 1;l<=n;l = r+1){
r = min(n,n/(n/l));
val1 = (val1+(r-l+1)*(l+r)/2%mod*(n/l)%mod)%mod;
}
for(int l = 1,r = 1;l<=m;l = r+1){
r = min(m,m/(m/l));
val2 = (val2+(r-l+1)*(l+r)/2%mod*(m/l)%mod)%mod;
}
for(int l = 1,r = 1;l<=n;l = r+1){
r = min(n,m/(m/l));
val3 = (val3+(r-l+1)*(l+r)/2%mod*(m/l)%mod)%mod;
}
for(int l = 1,r = 1;l<=n;l = r+1){
r = min({n,n/(n/l),m/(m/l)});
val4 = (val4+sum2(l,r)*(n/l)%mod*(m/l)%mod)%mod;
}
cout << (((((((n*n%mod*m%mod*m%mod-n*n%mod*val2%mod+mod)%mod-m*m%mod*val1%mod+mod)%mod+val1*val2%mod)%mod-n*n%mod*m%mod+mod)%mod+n*val3%mod)%mod+m*val1%mod)%mod-val4+mod)%mod;
return 0;
}
```
:::
## [[HAOI2011] Problem b](https://www.luogu.com.cn/problem/P2522)
### 题意
$n$ 次询问，每次给定 $a,b,c,d,k$，求
$$
\sum_{x=a}^b \sum_{y=c}^d [\gcd(x,y) = k]
$$
### 解析
设 $f(a,b) = \sum_{x=1}^a \sum_{y=1}^b [\gcd(x,y) = k]$，则原问题答案为 $f(b,d)-f(a-1,d)-f(b,c-1)+f(a-1,c-1)$．
不妨设 $a \leq b$，考虑计算 $f(a,b)$．
$$
\begin{aligned}
f(a,b)
&= \sum_{x=1}^a \sum_{y=1}^b [\gcd(x,y) = k] \\
&= \sum_{x=1}^a [k \mid x] \sum_{y=1}^b [k \mid y] [\gcd(\frac{x}{k},\frac{y}{k}) = 1] \\
&= \sum_{x=1}^{\lfloor \frac{a}{k} \rfloor} \sum_{y=1}^{\lfloor \frac{b}{k} \rfloor} [\gcd(x,y) = 1] \\
&= \sum_{x=1}^{\lfloor \frac{a}{k} \rfloor} \sum_{y=1}^{\lfloor \frac{b}{k} \rfloor} \sum_{d \mid \gcd(x,y)} \mu(d) \\
&= \sum_{d=1}^{\lfloor \frac{a}{k} \rfloor} \mu(d)\sum_{x=1}^{\lfloor \frac{a}{k} \rfloor} [d \mid x] \sum_{y=1}^{\lfloor \frac{b}{k} \rfloor} [d \mid y] \\
&= \sum_{d=1}^{\lfloor \frac{a}{k} \rfloor} \mu(d) \lfloor \frac{\lfloor \frac{a}{k} \rfloor}{d} \rfloor \lfloor \frac{\lfloor \frac{b}{k} \rfloor}{d} \rfloor
\end{aligned}
$$
整除分块即可．
:::success[AC code]
```cpp line-numbers
#include<bits/stdc++.h>
#define int long long
using namespace std;
constexpr int N = 5e4+1;
int n,a,b,c,d,k,mu[N];
int primes[N],len = 0;
bitset<N>st;
void init(){
mu[1] = 1;
for(int i = 2;i<N;++i){
if(!st[i]){
primes[++len] = i;
mu[i] = -1;
}
for(int j = 1;j<=len&&i*primes[j]<N;++j){
st[i*primes[j]] = true;
if(i%primes[j] == 0){
mu[i*primes[j]] = 0;
break;
}
mu[i*primes[j]] = -mu[i];
}
}
for(int i = 1;i<N;++i){
mu[i]+=mu[i-1];
}
}
int f(int a,int b){
int res = 0;
for(int l = 1,r = 1;l<=min(a/k,b/k);l = r+1){
r = min({min(a/k,b/k),(a/k)/((a/k)/l),(b/k)/((b/k)/l)});
res+=((a/k)/l)*((b/k)/l)*(mu[r]-mu[l-1]);
}
return res;
}
signed main(){
ios::sync_with_stdio(false);
cin.tie(0),cout.tie(0);
init();
cin >> n;
while(n--){
cin >> a >> b >> c >> d >> k;
cout << f(b,d)-f(a-1,d)-f(b,c-1)+f(a-1,c-1) << '\n';
}
return 0;
}
```
:::
## [[SDOI2015] 约数个数和](https://www.luogu.com.cn/problem/P3327)
### 题意
设 $\tau(x)$ 为 $x$ 的约数个数．$T$ 次询问，每次给定 $n,m$，求
$$
\sum_{i=1}^n \sum_{j=1}^m \tau(ij)
$$
### 解析
对于 $\tau(ij)$，有如下式子：
$$
\tau(ij) = \sum_{x \mid i} \sum_{y \mid j} [\gcd(x,y) = 1]
$$
证明略．感性理解吧．
所以根据这个式子，可对原式做莫反之类的东西．
不妨设 $n \leq m$，则
$$
\begin{aligned}
\sum_{i=1}^n \sum_{j=1}^m \tau(ij)
&= \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 \mid \gcd(x,y)} \mu(d) \\
&= \sum_{i=1}^n \sum_{j=1}^m \sum_{d=1}^n \mu(d) (\sum_{d \mid x \mid i} 1)(\sum_{d \mid y \mid j} 1)
\end{aligned}
$$
最右面两个式子是等价的，所以我们只考虑 $\sum_{d \mid x \mid i} 1$．
思考这个式子的意义：固定 $d,i$，计数有多少个 $x$ 使得 $d \mid x$ 且 $x \mid i$．
设 $x = k_1d,i = k_2d$．实际上，$k_2 = \frac{i}{d}$．问题等价于计数有多少个 $k_1$ 使得 $k_1 \mid k_2$．我们发现这就是 $\tau(k_2)$．
所以原式可以做如下变换：
$$
\begin{aligned}
\sum_{i=1}^n \sum_{j=1}^m \sum_{d=1}^n \mu(d) (\sum_{d \mid x \mid i} 1)(\sum_{d \mid y \mid j} 1)
&= \sum_{d=1}^n \mu(d) \sum_{i=1}^n [d \mid i]\tau(\frac{i}{d}) \sum_{j=1}^m [d \mid j]\tau(\frac{j}{d})\\
&= \sum_{d=1}^n \mu(d) \sum_{i=1}^{\lfloor \frac{n}{d} \rfloor} \tau(i) \sum_{j=1}^{\lfloor \frac{m}{d} \rfloor} \tau(j)
\end{aligned}
$$
我们预处理 $\mu,\tau$ 的前缀和，然后直接整除分块即可。
:::success[AC code]
```cpp line-numbers
#include<bits/stdc++.h>
#define int long long
using namespace std;
constexpr int N = 5e4+1;
int T,n,m,tau[N];
int primes[N],len = 0;
int mu[N];
bitset<N>st;
void init(){
mu[1] = tau[1] = 1;
for(int i = 2;i<N;++i){
if(!st[i]){
primes[++len] = i;
mu[i] = -1;
tau[i] = 2;
}
for(int j = 1;i*primes[j]<N;++j){
st[i*primes[j]] = true;
if(i%primes[j] == 0){
int t = i,cnt = 0;
while(t%primes[j] == 0){
++cnt;
t/=primes[j];
}
tau[i*primes[j]] = tau[t]*(cnt+2);
break;
}
mu[i*primes[j]] = -mu[i];
tau[i*primes[j]] = tau[i]*tau[primes[j]];
}
}
for(int i = 1;i<N;++i){
mu[i]+=mu[i-1];
tau[i]+=tau[i-1];
}
}
signed main(){
ios::sync_with_stdio(false);
cin.tie(0),cout.tie(0);
init();
cin >> T;
while(T--){
cin >> 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/(n/l),m/(m/l)});
res+=(mu[r]-mu[l-1])*tau[n/l]*tau[m/l];
}
cout << res << '\n';
}
return 0;
}
```
:::
## [[SDOI2017] 数字表格](https://www.luogu.com.cn/problem/P3704)
### 题意
定义斐波那契数列
$$
f_i=
\begin{cases}
0 & i=0 \\
1 & i=1 \\
f_{i-1}+f_{i-2} & i \geq 2
\end{cases}
$$
求
$$
\prod_{i=1}^n \prod_{j=1}^m f_{\gcd(i,j)}
$$
对 $(10^9+7)$ 取模的值．
### 解析
不妨设 $n \leq m$，考虑计数每一个 $\gcd(i,j)$ 出现的次数．
有
$$
\begin{aligned}
\prod_{i=1}^n \prod_{j=1}^m f_{\gcd(i,j)}
&= \prod_{g=1}^n f_{g}^{\sum_{i=1}^n \sum_{j=1}^m [\gcd(i,j) = g]}
\end{aligned}
$$
指数上的式子就很好处理了。
$$
\begin{aligned}
\sum_{i=1}^n \sum_{j=1}^m [\gcd(i,j) = g]
&= \sum_{i=1}^n [g \mid i] \sum_{j=1}^m [g \mid j] [\gcd(\frac{i}{g},\frac{j}{g}) = 1] \\
&= \sum_{i=1}^{\lfloor \frac{n}{g} \rfloor} \sum_{j=1}^{\lfloor \frac{m}{g} \rfloor} [\gcd(i,j) = 1] \\
&= \sum_{i=1}^{\lfloor \frac{n}{g} \rfloor} \sum_{j=1}^{\lfloor \frac{m}{g} \rfloor} \sum_{d \mid \gcd(i,j)} \mu(d) \\
&= \sum_{d=1}^{\lfloor \frac{n}{g} \rfloor}\mu(d)\sum_{i=1}^{\lfloor \frac{n}{g} \rfloor} [d \mid i] \sum_{j=1}^{\lfloor \frac{m}{g} \rfloor} [d \mid j] \\
&= \sum_{d=1}^{\lfloor \frac{n}{g} \rfloor}\mu(d) \lfloor \frac{n}{dg} \rfloor \lfloor \frac{m}{dg} \rfloor
\end{aligned}
$$
代回原式，有
$$
\begin{aligned}
\prod_{g=1}^n f_{g}^{\sum_{i=1}^n \sum_{j=1}^m [\gcd(i,j) = g]}
&= \prod_{g=1}^n f_{g}^{\sum_{d=1}^{\lfloor \frac{n}{g} \rfloor}\mu(d) \lfloor \frac{n}{dg} \rfloor \lfloor \frac{m}{dg} \rfloor} \\
&= \prod_{g=1}^n \prod_{d=1}^{\lfloor \frac{n}{g} \rfloor} f_g^{\mu(d) \lfloor \frac{n}{dg} \rfloor \lfloor \frac{m}{dg} \rfloor}
\end{aligned}
$$
转而枚举 $T = dg$，有
$$
\prod_{g=1}^n \prod_{d=1}^{\lfloor \frac{n}{g} \rfloor} f_g^{\mu(d) \lfloor \frac{n}{dg} \rfloor \lfloor \frac{m}{dg} \rfloor}=\prod_{T=1}^n \prod_{g \mid T}f_g^{\mu(\frac{T}{g}) \lfloor \frac{n}{T} \rfloor \lfloor \frac{m}{T} \rfloor}
$$
设 $F(T) = \prod_{g \mid T} f_g^{\mu(\frac{T}{g})}$，则答案为
$$
\prod_{T=1}^n F^{\lfloor \frac{n}{T} \rfloor \lfloor \frac{m}{T} \rfloor}(T)
$$
我们发现 $F(T)$ 可以 $O(n \ln n)$ 预处理，所以这题直接整除分块即可．时间复杂度为 $O(n \ln n + q \sqrt n \log n)$。
代码中给出的实现在预处理 $F$ 时重复调用了 ``qpow``，所以代码中的实现是 $O(n \log n \log mod+q \sqrt n \log mod)$ 的，不过依然能过就是了。
:::success[AC code]
```cpp line-numbers
#include<bits/stdc++.h>
#define int long long
using namespace std;
constexpr int N = 1e6+1,mod = 1e9+7;
int T,n,m,f[N],mu[N],F[N],pre[N],inv[N];
int primes[N],len = 0;
bitset<N>st;
inline int qpow(int base,int e){
int res = 1;
while(e){
if(e&1){
res = res*base%mod;
}
base = base*base%mod;
e>>=1;
}
return res;
}
void init(){
f[0] = 0,f[1] = 1;
for(int i = 2;i<N;++i){
f[i] = (f[i-1]+f[i-2])%mod;
}
mu[1] = 1;
for(int i = 2;i<N;++i){
if(!st[i]){
primes[++len] = i;
mu[i] = -1;
}
for(int j = 1;i*primes[j]<N;++j){
st[i*primes[j]] = true;
if(i%primes[j] == 0){
break;
}
mu[i*primes[j]] = -mu[i];
}
}
for(int i = 1;i<N;++i){
F[i] = 1;
}
for(int g = 1;g<N;++g){
for(int T = g;T<N;T+=g){
int val = mu[T/g];
if(val == 1){
F[T] = F[T]*f[g]%mod;
}else if(val == -1){
F[T] = F[T]*qpow(f[g],mod-2)%mod;
}
}
}
pre[0] = 1;
for(int i = 1;i<N;++i){
pre[i] = pre[i-1]*F[i]%mod;
}
inv[N-1] = qpow(pre[N-1],mod-2);
for(int i = N-2;i>=0;--i){
inv[i] = inv[i+1]*F[i+1]%mod;
}
}
signed main(){
ios::sync_with_stdio(false);
cin.tie(0),cout.tie(0);
init();
cin >> T;
while(T--){
cin >> n >> m;
if(n>m){
swap(n,m);
}
int res = 1;
for(int l = 1,r;l<=n;l = r+1){
r = min({n,n/(n/l),m/(m/l)});
res = res*qpow(pre[r]%mod*inv[l-1]%mod,(n/l)*(m/l))%mod;
}
cout << res << '\n';
}
return 0;
}
```
:::
## [[SDOI2014] 数表](https://www.luogu.com.cn/problem/P3312)
### 题意
$Q$ 次询问，每次给定 $n,m,a$，计算
$$
\sum_{i=1}^n \sum_{j=1}^m [\sigma(\gcd(i,j)) \leq a] \sigma(\gcd(i,j))
$$
对 $2^{31}$ 取模的值。
### 解析
不妨设 $n \leq m$，并将原式套路地转化：
$$
\begin{aligned}
\sum_{i=1}^n \sum_{j=1}^m [\sigma(\gcd(i,j)) \leq a] \sigma(\gcd(i,j))
&= \sum_{g=1}^n \sum_{i=1}^n \sum_{j=1}^m [\gcd(i,j) = g] [\sigma(g) \leq a] \sigma(g) \\
&= \sum_{g=1}^n [\sigma(g) \leq a] \sigma(g) \sum_{i=1}^{\lfloor \frac{n}{g} \rfloor} \sum_{j=1}^{\lfloor \frac{m}{g} \rfloor} [\gcd(i,j) = 1] \\
&= \sum_{g=1}^n [\sigma(g) \leq a] \sigma(g) \sum_{i=1}^{\lfloor \frac{n}{g} \rfloor} \sum_{j=1}^{\lfloor \frac{m}{g} \rfloor} \sum_{d \mid \gcd(i,j)} \mu(d) \\
&= \sum_{g=1}^n [\sigma(g) \leq a] \sigma(g) \sum_{d=1}^{\lfloor \frac{n}{g} \rfloor} \mu(d) \lfloor \frac{n}{dg} \rfloor \lfloor \frac{m}{dg} \rfloor \\
&= \sum_{T=1}^n \sum_{g \mid T} [\sigma(g) \leq a] \sigma(g) \mu(\frac{T}{g}) \lfloor \frac{n}{T} \rfloor \lfloor \frac{m}{T} \rfloor \\
\end{aligned}
$$
设 $F(T) = \sum_{g \mid T} [\sigma(g) \leq a] \sigma(g) \mu(\frac{T}{g})$，可将原式化为：
$$
\sum_{T=1}^n F(T) \lfloor \frac{n}{T} \rfloor \lfloor \frac{m}{T} \rfloor
$$
考虑如何计算 $F(T)$。
一个很显然的思路是把所有的询问都离线下来，按 $a$ 升序排序。所以我们可以对每一个 $\sigma(g)$ 在 $(a_{i-1},a_i]$ 之间的 $g$ 枚举倍数 $T$ 然后直接加，用树状数组维护即可（单点加，区间求和）。但是“枚举合法的 $g$”这一步是 $O(n)$ 的，这会导致崩溃，我们需要更优的方法。
我们发现我们可以将所有二元组 $(\sigma(g),g)$ 按第一维升序排序，然后我们维护一个指针就好了。最后是一个树状数组单点加。时间复杂度 $O(Q \sqrt n \log n + n \log^2 n)$。
:::success[AC code]
```cpp line-numbers
#include<bits/stdc++.h>
#define int long long
using namespace std;
constexpr int N = 1e5+1,mod = 1ll<<31;
int T,mu[N],sigma[N],pe[N],ans[N];
struct node{
int n,m,a,id;
}que[N];
pair<int,int>pii[N];
inline bool cmp(const node &x,const node &y){
return x.a<y.a;
}
int primes[N],len = 0;
bitset<N>st;
void init(){
mu[1] = sigma[1] = pe[1] = 1;
for(int i = 2;i<N;++i){
if(!st[i]){
primes[++len] = i;
mu[i] = -1;
sigma[i] = pe[i] = i+1;
}
for(int j = 1;i*primes[j]<N;++j){
st[i*primes[j]] = true;
if(i%primes[j] == 0){
pe[i*primes[j]] = pe[i]*primes[j]+1;
sigma[i*primes[j]] = sigma[i]/pe[i]*pe[i*primes[j]];
break;
}
pe[i*primes[j]] = primes[j]+1;
sigma[i*primes[j]] = sigma[i]*sigma[primes[j]];
mu[i*primes[j]] = -mu[i];
}
}
}
namespace BIT{
int tree[N];
inline int lowbit(const int &x){
return x&-x;
}
inline void add(const int &x,const int &y){
for(int i = x;i<N;i+=lowbit(i)){
tree[i] = (tree[i]+y)%mod;
}
}
inline int query(const int &x){
int sum = 0;
for(int i = x;i;i^=lowbit(i)){
sum = (sum+tree[i])%mod;
}
return sum;
}
}using namespace BIT;
signed main(){
ios::sync_with_stdio(false);
cin.tie(0),cout.tie(0);
init();
for(int i = 1;i<N;++i){
pii[i].first = sigma[i];
pii[i].second = i;
}
sort(pii+1,pii+N);
cin >> T;
for(int i = 1;i<=T;++i){
cin >> que[i].n >> que[i].m >> que[i].a;
que[i].id = i;
}
sort(que+1,que+T+1,cmp);
int pos = 0;
for(int i = 1;i<=T;++i){
while(pos<N-1&&pii[pos+1].first<=que[i].a){
++pos;
for(int j = pii[pos].second;j<N;j+=pii[pos].second){
add(j,(pii[pos].first%mod*mu[j/pii[pos].second]%mod+mod)%mod);
}
}
if(que[i].n>que[i].m){
swap(que[i].n,que[i].m);
}
int res = 0;
for(int l = 1,r;l<=que[i].n;l = r+1){
r = min(que[i].n/(que[i].n/l),que[i].m/(que[i].m/l));
res = (res+(que[i].n/l)*(que[i].m/l)%mod*((query(r)-query(l-1)+mod)%mod)%mod)%mod;
}
ans[que[i].id] = res;
}
for(int i = 1;i<=T;++i){
cout << ans[i] << '\n';
}
return 0;
}
```
:::
## [P3768 简单的数学题](https://www.luogu.com.cn/problem/P3768)
### 题意
给定正整数 $n,p$，求
$$
\sum_{i=1}^n \sum_{j=1}^n ij \gcd(i,j)
$$
对 $p$ 取模的值。
$1 \leq n \leq 10^{10},1 \leq p \leq 1.1 \times 10^9$ 且 $p$ 为质数。
### 解析
我们套路地对原式转化：
$$
\begin{aligned}
\sum_{i=1}^n \sum_{j=1}^n ij \gcd(i,j)
&= \sum_{g=1}^n g\sum_{i=1}^n i\sum_{j=1}^n j [\gcd(i,j) = g] \\
&= \sum_{g=1}^n g\sum_{i=1}^n i [g \mid i]\sum_{j=1}^n j [g \mid j][\gcd(\frac{i}{g},\frac{j}{g}) = 1] \\
&= \sum_{g=1}^n g\sum_{i=1}^{\lfloor \frac{n}{g} \rfloor} ig \sum_{j=1}^{\lfloor \frac{n}{g} \rfloor} jg [\gcd(i,j) = 1] \\
&= \sum_{g=1}^n g^3 \sum_{i=1}^{\lfloor \frac{n}{g} \rfloor} i \sum_{j=1}^{\lfloor \frac{n}{g} \rfloor} j \sum_{d \mid \gcd(i,j)} \mu(d) \\
&= \sum_{g=1}^n g^3 \sum_{d=1}^{\lfloor \frac{n}{g} \rfloor} \mu(d) \sum_{i=1}^{\lfloor \frac{n}{g} \rfloor} i [d \mid i] \sum_{j=1}^{\lfloor \frac{n}{g} \rfloor} j [d \mid j] \\
&= \sum_{g=1}^n g^3 \sum_{d=1}^{\lfloor \frac{n}{g} \rfloor} \mu(d) d^2(\sum_{i=1}^{\lfloor \frac{n}{dg} \rfloor} i)^2 \\
&= \sum_{g=1}^n g^3 \sum_{d=1}^{\lfloor \frac{n}{g} \rfloor} \mu(d) d^2(\frac{(1+\lfloor \frac{n}{dg} \rfloor) \lfloor \frac{n}{dg} \rfloor}{2})^2 \\
&= \sum_{T=1}^n \sum_{g \mid T} g^3 \mu(\frac{T}{g}) (\frac{T}{g})^2 (\frac{(1+\lfloor \frac{n}{T} \rfloor) \lfloor \frac{n}{T} \rfloor}{2})^2 \\
&= \sum_{T=1}^n (\frac{(1+\lfloor \frac{n}{T} \rfloor) \lfloor \frac{n}{T} \rfloor}{2})^2 T^2 \sum_{g \mid T} g\mu(\frac{T}{g}) \\
&= \sum_{T=1}^n (\frac{(1+\lfloor \frac{n}{T} \rfloor) \lfloor \frac{n}{T} \rfloor}{2})^2 T^2 \varphi(T)
\end{aligned}
$$
然后杜教筛优化即可。
令 $f(n) = n^2 \varphi(n)$，取 $g = id_2$，可得 $h = f*g = id_3$，即
$$
(f*g)(n) = \sum_{d \mid n} d^2 \varphi(d) (\frac{n}{d})^2 = n^2 \sum_{d \mid n} \varphi(d) = n^3 = id_3(n) = h(n)
$$
$g,h$ 都是可以 $O(1)$ 计算前缀和的。所以这题就做完了。
:::success[AC code]
```cpp line-numbers
#include<bits/stdc++.h>
#define ll long long
using namespace std;
constexpr int N = 4641589;
ll mod,n,inv2,inv6;
int phi[N];
ll F[N];
int primes[N],len = 0;
bitset<N>st;
inline ll qpow(int base,ll e){
ll res = 1;
while(e){
if(e&1){
res = res*base%mod;
}
base = 1ll*base*base%mod;
e>>=1;
}
return res;
}
void init(){
phi[1] = 1;
for(int i = 2;i<N;++i){
if(!st[i]){
primes[++len] = i;
phi[i] = (i-1)%mod;
}
for(int j = 1;i*primes[j]<N;++j){
st[i*primes[j]] = true;
if(i%primes[j] == 0){
phi[i*primes[j]] = phi[i]*primes[j]%mod;
break;
}
phi[i*primes[j]] = phi[i]*phi[primes[j]]%mod;
}
}
for(int i = 1;i<N;++i){
F[i] = i%mod*i%mod*phi[i]%mod;
F[i] = (F[i]+F[i-1])%mod;
}
}
inline ll calc1(ll n){
return n%mod*(n+1)%mod*inv2%mod;
}
inline ll calc2(ll n){
return n%mod*(n+1)%mod*(2*n%mod+1)%mod*inv6%mod;
}
map<ll,ll>mp;
ll SumF(ll n){
if(n<N){
return F[n];
}else if(mp.count(n)){
return mp[n];
}
ll ans = calc1(n)*calc1(n)%mod;
for(ll l = 2,r;l<=n;l = r+1){
r = n/(n/l);
ans = (ans-(calc2(r)-calc2(l-1)+mod)%mod*SumF(n/l)%mod+mod)%mod;
}
return mp[n] = ans;
}
signed main(){
ios::sync_with_stdio(false);
cin.tie(0),cout.tie(0);
cin >> mod >> n;
inv2 = qpow(2,mod-2);
inv6 = qpow(6,mod-2);
init();
ll ans = 0;
for(ll l = 1,r;l<=n;l = r+1){
r = n/(n/l);
ans = (ans+calc1(n/l)*calc1(n/l)%mod*((SumF(r)-SumF(l-1)+mod)%mod)%mod)%mod;
}
cout << ans;
return 0;
}
```
:::
## [P2424 约数和](https://www.luogu.com.cn/problem/P2424)
### 题意
给定正整数 $x,y$，计算
$$
\sum_{i=x}^y \sigma(i)
$$
### 解析
我们首先把这个东西拆成差分，这样下标就从 $1$ 开始了。
只需计算 $f(n) = \sum_{i=1}^n \sigma(i)$。
我们对 $f(n)$ 做如下变换：
$$
\begin{aligned}
\sum_{i=1}^n \sigma(i)
&= \sum_{i=1}^n \sum_{d \mid i} d \\
&= \sum_{d=1}^n d \sum_{d \mid i,i \leq n} 1\\
&= \sum_{d=1}^n d \lfloor \frac{n}{d} \rfloor
\end{aligned}
$$
:::success[AC code]
```cpp line-numbers
#include<bits/stdc++.h>
#define ll long long
using namespace std;
int x,y;
ll f(int n){
ll ans = 0;
for(ll l = 1,r;l<=n;l = r+1){
r = n/(n/l);
ans+=(l+r)*(r-l+1)/2*(n/l);
}
return ans;
}
signed main(){
ios::sync_with_stdio(false);
cin.tie(0),cout.tie(0);
cin >> x >> y;
cout << f(y)-f(x-1);
return 0;
}
```
:::
## [[POI 2007] ZAP-Queries](https://www.luogu.com.cn/problem/P3455)
### 题意
给定正整数 $a,b,g,g \leq a,b$，计算
$$
\sum_{i=1}^a \sum_{j=1}^b [\gcd(i,j) = g]
$$
### 解析
不妨设 $a \leq b$，套路地对原式转化：
$$
\begin{aligned}
\sum_{i=1}^a \sum_{j=1}^b [\gcd(i,j) = g]
&= \sum_{i=1}^a [g \mid i] \sum_{j=1}^b [g \mid j] [\gcd(\frac{i}{g},\frac{j}{g}) = 1] \\
&= \sum_{i=1}^{\lfloor \frac{a}{g} \rfloor} \sum_{j=1}^{\lfloor \frac{b}{g} \rfloor} [\gcd(i,j) = 1] \\
&= \sum_{i=1}^{\lfloor \frac{a}{g} \rfloor} \sum_{j=1}^{\lfloor \frac{b}{g} \rfloor} \sum_{d \mid \gcd(i,j)} \mu(d) \\
&= \sum_{d=1}^{\lfloor \frac{a}{g} \rfloor} \mu(d) \sum_{i=1}^{\lfloor \frac{a}{g} \rfloor} [d \mid i] \sum_{j=1}^{\lfloor \frac{b}{g} \rfloor} [d \mid j] \\
&= \sum_{d=1}^{\lfloor \frac{a}{g} \rfloor} \mu(d) \lfloor \frac{\lfloor \frac{a}{g} \rfloor}{d} \rfloor \lfloor \frac{\lfloor \frac{b}{g} \rfloor}{d} \rfloor
\end{aligned}
$$
然后数论分块即可。
:::success[AC code]
```cpp line-numbers
#include<bits/stdc++.h>
#define ll long long
using namespace std;
constexpr int N = 5e4+1;
int T,a,b,g;
int mu[N];
ll premu[N];
int primes[N],len = 0;
bitset<N>st;
void init(){
mu[1] = 1;
for(int i = 2;i<N;++i){
if(!st[i]){
primes[++len] = i;
mu[i] = -1;
}
for(int j = 1;i*primes[j]<N;++j){
st[i*primes[j]] = true;
if(i%primes[j] == 0){
break;
}
mu[i*primes[j]] = -mu[i];
}
}
for(int i = 1;i<N;++i){
premu[i] = premu[i-1]+mu[i];
}
}
signed main(){
ios::sync_with_stdio(false);
cin.tie(0),cout.tie(0);
init();
cin >> T;
while(T--){
cin >> a >> b >> g;
if(a>b){
swap(a,b);
}
ll res = 0,n = a/g,m = b/g;
for(ll l = 1,r;l<=n;l = r+1){
r = min(n/(n/l),m/(m/l));
res+=(premu[r]-premu[l-1])*(n/l)*(m/l);
}
cout << res << '\n';
}
return 0;
}
```
:::
## [P5221 Product](https://www.luogu.com.cn/problem/P5221)
### 题意
给定正整数 $n$，计算
$$
\prod_{i=1}^n \prod_{j=1}^n \frac{\operatorname{lcm}(i,j)}{\gcd(i,j)}
$$
对 $104857601$ 取模，一个质数。
### 题解
我们对原式做如下转化：
$$
\begin{aligned}
\prod_{i=1}^n \prod_{j=1}^n \frac{\operatorname{lcm}(i,j)}{\gcd(i,j)}
&= \prod_{i=1}^n \prod_{j=1}^n \frac{ij}{\gcd^2(i,j)} \\
&= \frac{\prod_{i=1}^n \prod_{j=1}^n ij}{\prod_{g=1}^n (g^2)^{\sum_{i=1}^n \sum_{j=1}^n [\gcd(i,j) = g]}} \\
&= \frac{(n!)^{2n}}{\prod_{g=1}^n (g^2)^{\sum_{d=1}^{\lfloor \frac{n}{g} \rfloor} \mu(d) \lfloor \frac{n}{dg} \rfloor^2}} \\
&= \frac{(n!)^{2n}}{\prod_{T=1}^n \prod_{g \mid T} (g^2)^{\mu(\frac{T}{g}) \lfloor \frac{n}{T} \rfloor^2}}
\end{aligned}
$$
令 $F(T) = \prod_{g \mid T} (g^2)^{\mu(\frac{T}{g})}$，则原式变为：
$$
\frac{(n!)^{2n}}{\prod_{T=1}^n \prod_{g \mid T} (g^2)^{\mu(\frac{T}{g}) \lfloor \frac{n}{T} \rfloor^2}} = \frac{(n!)^{2n}}{\prod_{T=1}^n F^{\lfloor \frac{n}{T} \rfloor^2}(T)}
$$
我们发现 $F(T)$ 可以 $O(n \ln n)$ 预处理，所以这题就做完了。
注意卡常。
:::success[AC code]
```cpp line-numbers
#include<bits/stdc++.h>
#define ll long long
using namespace std;
constexpr int N = 1e6+1,mod = 104857601;
int n;
inline int qpow(int base,int e){
int res = 1;
while(e){
if(e&1){
res = 1ll*res*base%mod;
}
base = 1ll*base*base%mod;
e>>=1;
}
return res;
}
int mu[N],F[N],preF[N],invF[N],inv[N];
int primes[N],len = 0;
bitset<N>st;
void init(){
mu[1] = inv[1] = 1;
for(int i = 2;i<=n;++i){
inv[i] = 1ll*(mod-mod/i)*inv[mod%i]%mod;
}
for(int i = 2;i<=n;++i){
if(!st[i]){
primes[++len] = i;
mu[i] = -1;
}
for(int j = 1;i*primes[j]<=n;++j){
st[i*primes[j]] = true;
if(i%primes[j] == 0){
break;
}
mu[i*primes[j]] = -mu[i];
}
}
for(int i = 1;i<=n;++i){
F[i] = 1;
}
for(ll g = 1;g<=n;++g){
for(int T = g;T<=n;T+=g){
if(mu[T/g] == 1){
F[T] = 1ll*F[T]*g%mod*g%mod;
}else if(mu[T/g] == -1){
F[T] = 1ll*F[T]*inv[g]%mod*inv[g]%mod;
}
}
}
preF[0] = 1;
for(int i = 1;i<=n;++i){
preF[i] = 1ll*preF[i-1]*F[i]%mod;
}
invF[n] = qpow(preF[n],mod-2);
for(int i = n-1;i>=0;--i){
invF[i] = 1ll*invF[i+1]*F[i+1]%mod;
}
}
signed main(){
ios::sync_with_stdio(false);
cin.tie(0),cout.tie(0);
cin >> n;
init();
int ans = 1;
for(int i = 1;i<=n;++i){
ans = 1ll*ans*i%mod;
}
ans = qpow(ans,2*n);
for(int l = 1,r;l<=n;l = r+1){
r = n/(n/l);
ans = 1ll*ans*qpow(qpow(1ll*preF[r]*invF[l-1]%mod,1ll*(n/l)*(n/l)%(mod-1)),mod-2)%mod;
}
cout << ans;
return 0;
}
```
:::
## [P4450 双亲数](https://www.luogu.com.cn/problem/P4450)
同 [[POI 2007] ZAP-Queries](https://www.luogu.com.cn/problem/P3455)。
## [[集训队互测 2010] Crash的数字表格](https://www.luogu.com.cn/problem/P1829)
### 题意
给定正整数 $n,m$，计算
$$
\sum_{i=1}^n \sum_{j=1}^m \operatorname{lcm}(i,j)
$$
对 $20101009$ 取模的值。
### 解析
不妨设 $n \leq m$，并套路地对原式转化：
$$
\begin{aligned}
\sum_{i=1}^n \sum_{j=1}^m \operatorname{lcm}(i,j)
&= \sum_{i=1}^n \sum_{j=1}^m \frac{ij}{\gcd(i,j)} \\
&= \sum_{g=1}^n \frac{1}{g} \sum_{i=1}^n i \sum_{j=1}^m j [\gcd(i,j) = g] \\
&= \sum_{g=1}^n \frac{1}{g} \sum_{i=1}^n i [g \mid i] \sum_{j=1}^m j [g \mid j] [\gcd(\frac{i}{g},\frac{j}{g}) = 1] \\
&= \sum_{g=1}^n g \sum_{i=1}^{\lfloor \frac{n}{g} \rfloor} i \sum_{j=1}^{\lfloor \frac{m}{g} \rfloor} j [\gcd(i,j) = 1] \\
&= \sum_{g=1}^n g \sum_{i=1}^{\lfloor \frac{n}{g} \rfloor} i \sum_{j=1}^{\lfloor \frac{m}{g} \rfloor} j \sum_{d \mid \gcd(i,j)} \mu(d) \\
&= \sum_{g=1}^n g \sum_{d=1}^{\lfloor \frac{n}{g} \rfloor} \mu(d) \sum_{i=1}^{\lfloor \frac{n}{g} \rfloor} i [d \mid i] \sum_{j=1}^{\lfloor \frac{m}{g} \rfloor} j [d \mid j] \\
&= \sum_{g=1}^n g \sum_{d=1}^{\lfloor \frac{n}{g} \rfloor} \mu(d) d^2 \frac{(1+\lfloor \frac{n}{dg} \rfloor) \lfloor \frac{n}{dg} \rfloor}{2} \cdot \frac{(1+\lfloor \frac{m}{dg} \rfloor) \lfloor \frac{m}{dg} \rfloor}{2} \\
&= \sum_{T=1}^n \sum_{g \mid T} g \mu(\frac{T}{g}) (\frac{T}{g})^2 \frac{(1+\lfloor \frac{n}{T} \rfloor) \lfloor \frac{n}{T} \rfloor}{2} \cdot \frac{(1+\lfloor \frac{m}{T} \rfloor) \lfloor \frac{m}{T} \rfloor}{2} \\
&= \sum_{T=1}^n \frac{(1+\lfloor \frac{n}{T} \rfloor) \lfloor \frac{n}{T} \rfloor}{2} \cdot \frac{(1+\lfloor \frac{m}{T} \rfloor) \lfloor \frac{m}{T} \rfloor}{2} \cdot T^2 \sum_{g \mid T} \mu(\frac{T}{g}) \frac{1}{g} \\
&= \sum_{T=1}^n \frac{(1+\lfloor \frac{n}{T} \rfloor) \lfloor \frac{n}{T} \rfloor}{2} \cdot \frac{(1+\lfloor \frac{m}{T} \rfloor) \lfloor \frac{m}{T} \rfloor}{2} \cdot T^2 \sum_{g \mid T} \mu(g) \frac{1}{\frac{T}{g}} \\
&= \sum_{T=1}^n \frac{(1+\lfloor \frac{n}{T} \rfloor) \lfloor \frac{n}{T} \rfloor}{2} \cdot \frac{(1+\lfloor \frac{m}{T} \rfloor) \lfloor \frac{m}{T} \rfloor}{2} \cdot T \sum_{g \mid T} \mu(g) g
\end{aligned}
$$
函数 $F(T) = \sum_{g \mid T} \mu(g) g$ 可以线性预处理，其是一个积性函数。
具体地：
- $F(1) = 1
-
F(p) = 1-p
-
F(p^k) = F(p)(k \geq 2)
:::success[AC code]
#include<bits/stdc++.h>
#define ll long long
using namespace std;
constexpr int N = 1e7+1,mod = 20101009;
int n,m;
int inv[N],F[N];
int primes[N],len = 0;
bitset<N>st;
void init(){
inv[1] = F[1] = 1;
for(int i = 2;i<=max(m,6);++i){
inv[i] = 1ll*(mod-mod/i)*inv[mod%i]%mod;
}
for(int i = 2;i<=m;++i){
if(!st[i]){
primes[++len] = i;
F[i] = (1-i+mod)%mod;
}
for(int j = 1;i*primes[j]<=m;++j){
st[i*primes[j]] = true;
if(i%primes[j] == 0){
F[i*primes[j]] = F[i];
break;
}
F[i*primes[j]] = 1ll*F[i]*F[primes[j]]%mod;
}
}
for(int i = 1;i<=m;++i){
F[i] = 1ll*F[i]*i%mod;
F[i] = (F[i]+F[i-1])%mod;
}
}
inline int calc1(int n){
return 1ll*n*(n+1)%mod*inv[2]%mod;
}
inline int calc2(int n){
return 1ll*n*(n+1)%mod*(2*n+1)%mod*inv[6]%mod;
}
signed main(){
ios::sync_with_stdio(false);
cin.tie(0),cout.tie(0);
cin >> n >> m;
if(n>m){
swap(n,m);
}
init();
int res = 0;
for(int l = 1,r;l<=n;l = r+1){
r = min(n/(n/l),m/(m/l));
res = (res+1ll*calc1(n/l)*calc1(m/l)%mod*((F[r]-F[l-1]+mod)%mod)%mod)%mod;
}
cout << res;
return 0;
}
:::
P4449 于神之怒加强版
题意
给定 k，有 T 次询问，每次询问给定 n,m，计算
\sum_{i=1}^n \sum_{j=1}^m \gcd(i,j)^k
对 10^9+7 取模的值。
解析
不妨设 n \leq m，并套路地对原式转化：
\begin{aligned}
\sum_{i=1}^n \sum_{j=1}^m \gcd(i,j)^k
&= \sum_{g=1}^n g^k \sum_{i=1}^n \sum_{j=1}^m [\gcd(i,j) = g] \\
&= \sum_{g=1}^n g^k \sum_{i=1}^{\lfloor \frac{n}{g} \rfloor} \sum_{j=1}^{\lfloor \frac{m}{g} \rfloor} [\gcd(i,j) = 1] \\
&= \sum_{g=1}^n g^k \sum_{i=1}^{\lfloor \frac{n}{g} \rfloor} \sum_{j=1}^{\lfloor \frac{m}{g} \rfloor} \sum_{d \mid \gcd(i,j)} \mu(d) \\
&= \sum_{g=1}^n g^k \sum_{d=1}^{\lfloor \frac{n}{g} \rfloor} \mu(d) \lfloor \frac{n}{dg} \rfloor \lfloor \frac{m}{dg} \rfloor \\
&= \sum_{T=1}^n \lfloor \frac{n}{T} \rfloor \lfloor \frac{m}{T} \rfloor \sum_{g \mid T} g^k \mu(\frac{T}{g})
\end{aligned}
这题能 O(n \ln n)，就不继续优化了。
:::success[AC code]
#include<bits/stdc++.h>
#define ll long long
using namespace std;
constexpr int N = 5e6+1,mod = 1e9+7;
int t,k,n,m;
int primes[N],len = 0;
int mu[N],idk[N],F[N];
bitset<N>st;
inline int qpow(int base,int e){
int res = 1;
while(e){
if(e&1){
res = 1ll*res*base%mod;
}
base = 1ll*base*base%mod;
e>>=1;
}
return res;
}
void init(){
mu[1] = idk[1] = 1;
for(int i = 2;i<N;++i){
if(!st[i]){
primes[++len] = i;
mu[i] = mod-1;
idk[i] = qpow(i,k);
}
for(int j = 1;i*primes[j]<N;++j){
st[i*primes[j]] = true;
idk[i*primes[j]] = 1ll*idk[i]*idk[primes[j]]%mod;
if(i%primes[j] == 0){
break;
}
mu[i*primes[j]] = (-mu[i]+mod)%mod;
}
}
for(int g = 1;g<N;++g){
for(int T = g;T<N;T+=g){
F[T] = (F[T]+1ll*idk[g]*mu[T/g]%mod+mod)%mod;
}
}
for(int i = 1;i<N;++i){
F[i] = (F[i]+F[i-1])%mod;
}
}
signed main(){
ios::sync_with_stdio(false);
cin.tie(0),cout.tie(0);
cin >> t >> k;
init();
while(t--){
cin >> 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+1ll*(n/l)*(m/l)%mod*((F[r]-F[l-1]+mod)%mod)%mod)%mod;
}
cout << res << '\n';
}
return 0;
}
:::
P4917 天守阁的地板
题意
$$
\prod_{i=1}^n \prod_{j=1}^n \frac{\operatorname{lcm}(i,j)^2}{ij}
$$
对 $19260817$ 取模的值。
### 解析
$$
\prod_{i=1}^n \prod_{j=1}^n \frac{\operatorname{lcm}(i,j)^2}{ij}= \prod_{i=1}^n \prod_{j=1}^n \frac{ij}{\gcd(i,j)^2} \\
$$
然后同 [P5221 Product](https://www.luogu.com.cn/problem/P5221)。