浅谈 OI 中的数论
xwh_hh
·
·
算法·理论
一些记号
$\gcd(a,b)$ 表示 $a,b$ 的最大公约数。
$[x=k]$ 表示函数 $f(x)=\begin{cases}1&x=k\\0&x\ne k\end{cases}
素数的基本知识
素数定义
因数个数为 2 的数为素数（质数）。
素数判定
朴素判定
也可以先打素数表，单次查询可降至 $O(\frac{\sqrt n}{\ln n})$。
#### Miller-Rabin 素数测试
对于素数 $p>2$ 而言，若 $x^2\equiv 1\pmod p$，则有 $x\equiv \pm1\pmod p$。
证明：等价于 $p|(x-1)(x+1)$，而 $p>(x+1)-(x-1)\ge \gcd(x+1,x-1)$，故只能 $p|x-1$ 或 $p|x+1$。
我们知道 $a^{p-1}\equiv1\pmod p$，可将 $p-1$ 分解为 $2^kq$ 的形式，其中 $q$ 为奇数，先计算 $a^q$，然后不断平方，看它是否由 $-1$ 变为 $1$。
[模板题](https://loj.ac/p/143)代码实现（`long long` 范围内 $100\%$ 正确）：
::::info[Code]
```cpp
#include<bits/stdc++.h>
using namespace std;
long long d[]={2,3,5,7,11,13,17,19,23,29,31,37};
long long pw(long long a,long long b,long long p){
long long res=1;
while(b){
if(b&1) res=(__int128)(res)*a%p;
a=(__int128)(a)*a%p;
b>>=1;
}
return res;
}
bool is_prime(long long x){
for(long long i=0;i<12;i++){
long long u=d[i];
if(x==u) return true;
if(x<u) return false;
if(pw(u,x-1,x)!=1) return false;
long long pos=x-1;
while(!(pos&1)) pos>>=1;
long long now=pw(u,pos,x);
while(pos!=x-1){
if(now==x-1) break;
if(now==1) break;
if((__int128)(now)*now%x==1) return false;
now=(__int128)(now)*now%x;
pos<<=1;
}
}
return true;
}
signed main(){
long long p;
while(cin>>p){
cout<<(is_prime(p)?'Y':'N')<<'\n';
}
return 0;
}
```
::::
## Lucas 定理、exLucas 定理
### Lucas
$$\Large C_{n}^{m}\equiv C_{[n/p]}^{[m/p]}C_{n\bmod p}^{m\bmod p}\pmod p$$
这里认为 $m>n\Rightarrow C_n^m=0$。
### exLucas
假如模数不是质数怎么办？
考虑质因数分解，只需处理模数为素数幂次的问题。
预处理每个数阶乘，将其表示为 $p^kq$ 的形式，其中 $\gcd(p,q)=1$，后略，可参考 [oi-wiki](https://oi-wiki.org/math/number-theory/lucas/#exlucas-%E7%AE%97%E6%B3%95)。
[exLucas 模板题](https://www.luogu.com.cn/problem/P4720)代码：
::::info[Code]
```cpp
#include<bits/stdc++.h>
using namespace std;
typedef __int128 LL;
namespace io{
inline __int128 read(){
__int128 ans=0,f=1;
char c=getchar();
while(c<'0'||c>'9'){
if(c=='-') f=-f;
c=getchar();
}
while(c>='0'&&c<='9'){
ans=(ans<<3)+(ans<<1)+c-48;
c=getchar();
}
return ans*f;
}
inline char readchar(){
bk:
char c=getchar();
if(c>32&&c<127) return c;
if(c==EOF) return 0;
goto bk;
}
void write(__int128 x){
if(x<0){
putchar('-');
x=-x;
}
if(x<10) putchar(x+48);
else{
write(x/10);
putchar(x%10+48);
}
}
}
namespace exLucas{
LL gcd(LL x,LL y){
if(y==0) return x;
return gcd(y,x%y);
}
LL pw(LL a,LL b,LL p=0x3f3f3f3f3f3f3f3f){
LL res=1,now=a;
while(b){
if(b&1) res=res*now%p;
now=now*now%p;
b>>=1;
}
return res;
}
LL mod(LL a,LL b,LL m){
if(m==0) exit(3);
a%=m,b%=m;
if(a==1) return b;
if(a==0&&b==0) return 0;
if(a==0) exit(4);
LL k=mod((a-1)*m,b,a);
return ((b+m*k)/a)%m;
}
LL excrt(vector<LL>a,vector<LL>b){
LL xa=0,xb=1;
LL pos=0;
while(pos!=a.size()){
LL na,nb,ee;
na=a[pos++]; nb=b[pos-1];
LL xx=xb,yy=((na-xa)%nb+nb)%nb,zz=nb;
ee=gcd(xx,zz);
if(ee==0) exit(2);
if(yy%ee!=0){
exit(1);
}
xx/=ee;
yy/=ee;
zz/=ee;
LL k=mod(xx,yy,zz);
xa=xa+xb*k;
xb*=nb/gcd(xb,nb);
xa%=xb;
}
return xa;
}
typedef pair<LL,LL>P;
vector<P>t;
LL f(LL n,LL p){
if(n<p) return 0;
return n/p+f(n/p,p);
}
LL g(LL n,LL p,LL k){
if(n==0) return 1;
LL t=pw(p,k),now=1;
for(LL i=1;i<=t;i++){
LL j=i;
if(j%p!=0)
now=now*j%t;
}
now=pw(now,n/t,t);
now=now*g(n/p,p,k)%t;
for(LL i=1;i<=n%t;i++){
LL j=i;
if(j%p!=0)
now=now*j%t;
}
assert(now);
return now;
}
LL C(LL x,LL y,LL p,LL k){
LL t=f(x,p)-f(y,p)-f(x-y,p);
if(t>=k) return 0;
k-=t;
LL tmp=pw(p,k);
assert(tmp);
LL ans=(g(x,p,k)*mod(g(y,p,k),1,tmp)%tmp*mod(g(x-y,p,k),1,tmp)%tmp)*pw(p,t);
return ans;
}
LL exlucas(LL n,LL m,LL p){
t.clear();
LL tmp=p;
for(LL i=2;i*i<=tmp;i++){
LL cnt=0;
while(tmp%i==0){
tmp/=i;
cnt++;
}
if(cnt) t.push_back(P(i,cnt));
}
if(tmp!=1) t.push_back(P(tmp,1));
vector<LL>a,b;
for(LL i=0;i<t.size();i++){
LL u=t[i].first,v=t[i].second;
a.push_back(C(n,m,u,v)),b.push_back(pw(u,v));
}
return excrt(a,b);
}
}
signed main(){
LL n=io::read(),m=io::read(),p=io::read();
io::write(exLucas::exlucas(n,m,p)); putchar(10);
return 0;
}
```
::::
## 一些数论函数
### 定义
$1(x)$ 表示函数值恒为 $1$ 的函数。
$d(x)$ 表示 $x$ 的正因数个数。
$\varepsilon (x)=[x=1]$。
$\mu(x)=\begin{cases}1&x\ \text{的质因子个数为偶数且互不相同}\\0&x\ \text{存在平方因子}\\-1&x\ \text{的质因子个数为奇数且互不相同}\end{cases}
$\varphi(x)$ 表示 $1$ 到 $x$ 中与 $x$ 互质的数的个数。
### 举例
以防你看不懂，下面列举部分函数值：
::::info[表格]
| $x$ | $1(x)$ | $\varepsilon(x)$ | $\mu(x)$ | $id(x)$ | $\varphi(x)$ |$d(x)$|
|:-:|:-:|:-:|:-:|:-:|:-:|:-:|
| $1$ | $1$ | $1$ | $1$ | $1$ | $1$ | $1$ |
| $2$ | $1$ | $0$ | $-1$ | $2$ | $1$ | $2$ |
| $3$ | $1$ | $0$ | $-1$ | $3$ | $2$ | $2$ |
| $4$ | $1$ | $0$ | $0$ | $4$ | $2$ | $3$ |
| $5$ | $1$ | $0$ | $-1$ | $5$ | $4$ | $2$ |
| $6$ | $1$ | $0$ | $1$ | $6$ | $2$ | $4$ |
| $7$ | $1$ | $0$ | $-1$ | $7$ | $6$ | $2$ |
| $8$ | $1$ | $0$ | $0$ | $8$ | $4$ | $4$ |
| $9$ | $1$ | $0$ | $0$ | $9$ | $6$ | $3$ |
| $10$ | $1$ | $0$ | $1$ | $10$ | $4$ | $4$ |
::::
### 性质
如果对于任意**互质**的正整数 $x,y$，有 $f(x)f(y)=f(xy)$，称 $f$ 为积性函数。
如果对于任意正整数 $x,y$，有 $f(x)f(y)=f(xy)$，称 $f$ 为完全积性函数。
上述函数均为积性函数，其中 $1(x),\varepsilon(x),id(x)$ 为完全积性函数。
## 莫比乌斯反演
### 狄利克雷卷积定义
若有
$$f(x)=\sum_{d|x}g(d)h(\frac{x}{d})$$
恒成立，则认为 $f$ 为 $g$ 与 $h$ 的卷积，记作 $f=g\ast h$。
### 重要性质
$f,g,h$ 中有任意两个为积性函数，可推出第三个亦为积性函数。
$$\mu\ast1=\varepsilon$$
$$1\ast1=d$$
$$1\ast\varphi=id$$
$$d(xy)=\sum_{i|x}\sum_{j|y}[\gcd(i,j)=1]$$
## 筛法
### 素数筛
$O(n\log\log n)$ 的埃氏筛和 $O(n)$ 的欧拉筛，想必大佬们都会。
### 杜教筛
$O(n^{\frac 23})$ 求积性函数前缀和。
$$f\ast g=h$$
$$F(n)=\sum_{i=1}^nf(i)$$
$$G(n)=\sum_{i=1}^ng(i)$$
$$H(n)=\sum_{i=1}^nh(i)$$
$$g(1)F(n)=g(1)\sum_{i=1}^nf(i)\\=\sum_{i=1}^nh(i)-\sum_{i=2}^nF([\frac ni])g(i)$$
整除分块。
预处理可平衡复杂度至 $O(n^\frac 23)$。
[模板题](https://www.luogu.com.cn/problem/P4213)代码：
::::info[Code]
```cpp
#include<bits/stdc++.h>
using namespace std;
long long sumI(long long l,long long r){
return r-l+1;
}
long long sumid(long long l,long long r){
return (l+r)*(r-l+1)/2;
}
long long sumx(long long l,long long r){
return l<=1&&r>=1;
}
long long hashh(long long t){
return t%47224349;
}
long long mp[50000005];
long long a[15],b[15],c[15];
long long s[1500005],t[1500005],pf[1500005],u[1500005];
long long pr[100005],cnt;
long long solve(long long n,long long (*g)(long long l,long long r),long long(*h)(long long l,long long r)){//f*g=h
long long w=hashh(n);
if(n<=1500000) return h(159,159)==159?t[n]:u[n];
if(mp[w]) return mp[w];
mp[w]=h(1,n);
long long u=1;
for(long long i=2;i*i<=n;i++){
u=i;
mp[w]-=g(i,i)*solve(n/i,g,h);
}
for(long long i=1;i<=u;i++){
long long l=n/(i+1)+1,r=n/i;
if(l<=u) l=u+1;
mp[w]-=g(l,r)*solve(i,g,h);
}
return mp[w];
}
void del(long long n){//f*g=h
long long w=hashh(n);
if(!mp[w]) return;
mp[w]=0;
long long u=1;
for(long long i=2;i*i<=n;i++){
u=i;
del(n/i);
del(i);
}
return ;
}
signed main(){
ios::sync_with_stdio(false); cin.tie(0); cout.tie(0);
for(long long i=2;i<=1500000;i++){
if(!s[i]){
pr[++cnt]=i;
s[i]=i;
long long now=i;
while(now<=1500000){
t[now]=now/i*(i-1);
now*=i;
}
}
if(i/s[i]%s[i]==0) pf[i]=1;
for(long long j=1;j<=cnt;j++){
if(i*pr[j]>1500000) break;
s[i*pr[j]]=pr[j];
if(s[i]==pr[j]) break;
}
}
for(long long i=2;i<=1500000;i++){
if(t[i]) continue;
long long u=s[i],v=i/s[i];
while(v%s[i]==0){
u*=s[i];
v/=s[i];
}
t[i]=t[u]*t[v];
}
u[1]=1;
t[1]=1;
for(long long i=2;i<=1500000;i++){
if(!pf[i]) u[i]=-u[i/s[i]];
else u[i]=0;
}
// for(long long i=1;i<=100;i++){
// cout<<i<<' '<<s[i]<<' '<<t[i]<<' '<<u[i]<<'\n';
// }
for(long long i=2;i<=1500000;i++){
u[i]+=u[i-1];
t[i]+=t[i-1];
}
// exit(0);
long long t;
cin>>t;
for(long long i=1;i<=t;i++) cin>>a[i];
for(long long i=1;i<=t;i++){
b[i]=solve(a[i],sumI,sumid);
}
// memset(mp,0,sizeof(mp));
for(long long i=1;i<=t;i++) del(a[i]);
// for(long long i=1;i<=50000;i++) solve(i,sumI,sumx);
for(long long i=1;i<=t;i++){
c[i]=solve(a[i],sumI,sumx);
}
for(long long i=1;i<=t;i++){
cout<<b[i]<<' '<<c[i]<<'\n';
}
/*
while(t--){
cin>>n;
cout<<solve(n,sumI,sumid)<<' ';
del(n);
cout<<solve(n,sumI,sumx)<<'\n';
del(n);
}
*/
return 0;
}
/*
n/l<i+1 n/(l-1)>=i+1 n/(i+1)<l n/(i+1)>=l-1
*/
```
::::
### PN 筛
定义形如 $a^2b^3$ 的数为 Powerful Number（即不含一次素因子的数）。
结论：PN 个数为 $O(\sqrt n)$。
想求 $\sum_{i=1}^n f(i)$，构造 $g$ 使得 $\forall p\in prime$ 有 $f(p)=g(p)$，构造 $h\ast g=f$，则有 $h$ 为积性函数，且 $\forall i\notin PN,h(i)=0$，构造 $G(u)=\sum_{i=1}^ug(i)$。
$$\large\sum_{i=1}^nf(i)=\sum_{i=1}^n\sum_{d|i}g(d)h(\frac id)\\=\sum_{i=1}^n\sum_{j=1}^{\frac ni}g(j)h(i)\\=\sum_{i=1,i\in PN}^nh(i)G(\frac ni)$$
## 例题选讲
### [P2257](https://www.luogu.com.cn/problem/P2257)
#### 题面
计算：
$$\large \sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)=prime]$$
#### 解法
莫比乌斯反演，下面式子中 $p$ 均为质数。
$$\large \sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)=prime]\\=\sum_{p=2}^{\min(n,m)}\sum_{i=1}^{\frac np}\sum_{j=1}^{\frac mp}[\gcd(i,j)=1]\\=\sum_{p=2}^{\min(n,m)}\sum_{i=1}^{\frac np}\sum_{j=1}^{\frac mp}\sum_{k|i,k|j}\mu(k)\\=\sum_{p=2}^{\min(n,m)}\sum_{k=1}^{\frac{\min(n,m)}p}[\frac n{pk}][\frac m{pk}]\mu(k)\\=\sum_{T=2}^{\min(n,m)}\sum_{p|T}[\frac nT][\frac mT]\mu(\frac Tp)\\=\sum_{T=2}^{\min(n,m)}a_T[\frac nT][\frac mT]$$
整除分块即可，$O(T\sqrt{\max(n,m)})$。
预处理 $a_T=\sum_{p|T}\mu(\frac Tp)$， $O(n\log\log n)$。
本题轻微卡常。
::::info[Code]
```cpp
#include<bits/stdc++.h>
using namespace std;
long long a[10000005];
long long p[10000005];
long long sum[10000005];
long long mu[10000005];
signed main(){
ios::sync_with_stdio(false); cin.tie(0); cout.tie(0);
mu[1]=1;
for(long long i=2;i<=10000000;i++){
if(a[i]==0){
mu[i]=-1;
p[i]=i;
for(long long j=i*i;j<=10000000;j+=i){
p[j]=i;
a[j]++;
}
}
}
for(long long i=2;i<=10000000;i++){
a[i]=0;
if(!mu[i]){
if(i/p[i]%p[i]==0) continue;
mu[i]=-mu[i/p[i]];
}
}
for(long long i=1;i<=10000000;i++){
if(p[i]==i){
for(long long j=i;j<=10000000;j+=i){
a[j]+=mu[j/i];
}
}
sum[i]=sum[i-1]+a[i];
}
long long t;
cin>>t;
while(t--){
long long n,m,ans=0;
cin>>n>>m;
for(long long i=2;i<=n&&i<=m&&i<=2000;i++) ans+=(n/i)*(m/i)*a[i];
long long l=2001,r;
while(l<=n&&l<=m){
r=min(n/(n/l),m/(m/l));
ans+=(n/l)*(m/l)*(sum[r]-sum[l-1]);
l=r+1;
}
cout<<ans<<'\n';
}
return 0;
}
```
::::
### [P4240](https://www.luogu.com.cn/problem/P4240)
#### 题面
$T$ 次询问，每次给定 $n, m$，Salamander 需要回答出 $\left( \sum_{i=1}^n \sum_{j=1}^m \varphi(ij) \right)\! \bmod 998244353$。
对于 $100\%$ 的数据，$1 \le T
\le {10}^4$，$1 \le n, m \le {10}^5$。
#### 解法
$$\varphi(xy)=\varphi(x)\varphi(y)\frac{\gcd(x,y)}{\varphi(\gcd(x,y))}$$
$$\large\sum_{i=1}^n \sum_{j=1}^m \varphi(ij)=\sum_{x=1}^n\sum_{y=1}^m\varphi(x)\varphi(y)\frac{\gcd(x,y)}{\varphi(\gcd(x,y))}\\=\sum_{T=1}^{\min(n,m)}\frac T{\varphi(T)}\sum_{i=1}^\frac nT\sum_{j=1}^\frac mT\varphi(iT)\varphi(jT)[\gcd(i,j)=1]\\=\sum_{T=1}^{\min(n,m)}\frac T{\varphi(T)}\sum_{t=1}^\frac nT\mu(t)\sum_{i=1}^\frac n{Tt}\varphi(iTt)\sum_{j=1}^\frac m{Tt}\varphi(jTt)\\=\sum_{k=1}^{\min(n,m)}\sum_{i=1}^{\frac{n}{k}}\varphi(ik)\sum_{j=1}^{\frac{m}{k}}\varphi(jk)\sum_{d|k}\frac{d\mu(\frac kd)}{\varphi(d)}$$
记 $\large a_k=\sum_{d|k}\frac{d\mu(\frac kd)}{\varphi(d)}$（可 $O(n\log n)$ 预处理），$\large b_{t,w}=\sum_{i=1}^{t}\varphi(wi)$，则即求：
$$\large\sum_{k=1}^{\min(n,m)}b_{\frac nk,k}b_{\frac mk,k}a_k$$
这个式子 $b$ 里面带 $k$，不能整除分块，我们令
$$\large t_{x,y,z}=\sum_{i=1}^xa_ib_{y,i}b_{z,i}$$
根号分治，后略。
::::info[Code]
```cpp
#include<bits/stdc++.h>
using namespace std;
const long long p=998244353;
const long long B=18;
long long phi[100005],mu[100005];
long long a[100005],sum[100005];
vector<long long>b[100005];
long long t[100005][20][20];
long long pw(long long a,long long b=998244351,long long p=998244353){
long long res=1;
while(b){
if(b&1) res=res*a%p;
a=a*a%p;
b>>=1;
}
return res;
}
signed main(){
//O(nlogn) 筛出phi,mu
for(long long i=1;i<=100000;i++) phi[i]=i;
mu[1]=1;
for(long long i=1;i<=100000;i++){
for(long long j=i+i;j<=100000;j+=i){
mu[j]-=mu[i];
phi[j]-=phi[i];
}
}
for(long long i=1;i<=100000;i++){
for(long long j=i;j<=100000;j+=i){//i向j贡献
a[j]+=i*mu[j/i]*pw(phi[i]);
a[j]%=p;
}
a[i]+=p;
a[i]%=p;
sum[i]=(sum[i-1]+a[i])%p;
}
for(long long i=1;i<=100000;i++) b[i].push_back(0);
for(long long w=1;w<=100000;w++){
for(long long t=1;w*t<=100000;t++){
b[w].push_back((b[w][b[w].size()-1]+phi[w*t])%p);
}
}
for(long long x=1;x<=100000;x++){
for(long long y=1;y<=B&&y*x<=100000;y++){
for(long long z=1;z<=B&&z*x<=100000;z++){
assert(b[y].size()>x);
assert(b[z].size()>x);
t[x][y][z]=t[x-1][y][z]+a[x]*b[x][y]%p*b[x][z]%p;
t[x][y][z]%=p;
}
}
}
ios::sync_with_stdio(false); cin.tie(0); cout.tie(0);
long long q;
cin>>q;
while(q--){
long long n,m;
cin>>n>>m;
if(min(n,m)<5600){
long long ans=0;
for(long long i=1;i<=min(n,m);i++){
ans+=a[i]*b[i][n/i]%p*b[i][m/i]%p;
ans%=p;
}
cout<<ans<<'\n';
}else{
long long ans=0;
for(long long i=1;i<=5600;i++){
ans+=a[i]*b[i][n/i]%p*b[i][m/i]%p;
// cout<<a[i]*b[i][n/i]%p*b[i][m/i]<<'\n';
ans%=p;
}
long long l=5601,r;
while(l<=n&&l<=m){
r=min(n/(n/l),m/(m/l));
ans+=t[r][n/l][m/l]-t[l-1][n/l][m/l];
// cout<<t[r][n/l][m/l]-t[l-1][n/l][m/l]<<'\n';
ans%=p;
ans+=p;
ans%=p;
l=r+1;
}
cout<<ans<<'\n';
}
}
return 0;
}
```
::::
### [P3726](https://www.luogu.com.cn/problem/P3726)
记小 A 向上 $s$ 次，向下 $(a-s)$ 次，小 B 向上 $(b-t)$ 次，向下 $t$ 次，只需 $s+t>b$ 即可。
$$Ans=\sum_{i=b+1}^{a+b} C_{a+b}^i=\frac 12\sum_{i=b+1}^{a+b}C_{a+b}^{i}+\frac 12\sum_{i=0}^{a-1}C_{a+b}^i\\=\frac{2^{a+b}+\sum_{i=b+1}^{a-1}C_{a+b}^i}2$$
exLucas 求解即可。
喜提洛谷第四劣解，并喜提 $2$ 发 $95$ 分 $1.01s\ TLE$。
::::info[Code]
```cpp
#include<bits/stdc++.h>
using namespace std;
typedef __int128 LL;
namespace exLucas{
LL id[]={0,4,5,8,16,25,32,64,125,128,256,512,625,1024,3125,15625,78125,390625,1953125};
LL len=18;
vector<LL>a[20];
LL gcd(LL x,LL y){
if(y==0) return x;
return gcd(y,x%y);
}
void init(){
for(LL i=1;i<=len;i++){
LL now=id[i];
a[i].resize(now+1);
a[i][0]=1;
for(LL j=1;j<=now;j++){
a[i][j]=a[i][j-1]*(now%5==0&&j%5==0||now%2==0&&j%2==0?1:j)%now;
}
}
}
LL pw(LL a,LL b,LL p=0x3f3f3f3f3f3f3f3f){
LL res=1,now=a;
while(b){
if(b&1) res=res*now%p;
now=now*now%p;
b>>=1;
}
return res;
}
LL mod(LL a,LL b,LL m){
if(m==0) exit(3);
a%=m,b%=m;
if(a==1) return b;
if(a==0&&b==0) return 0;
if(a==0) exit(4);
LL k=mod((a-1)*m,b,a);
return ((b+m*k)/a)%m;
}
LL excrt(vector<LL>a,vector<LL>b){
LL xa=0,xb=1;
LL pos=0;
while(pos!=a.size()){
LL na,nb,ee;
na=a[pos++]; nb=b[pos-1];
LL xx=xb,yy=((na-xa)%nb+nb)%nb,zz=nb;
ee=gcd(xx,zz);
if(ee==0) exit(2);
if(yy%ee!=0){
exit(1);
}
xx/=ee;
yy/=ee;
zz/=ee;
LL k=mod(xx,yy,zz);
xa=xa+xb*k;
xb*=nb/gcd(xb,nb);
xa%=xb;
}
return xa;
}
typedef pair<LL,LL>P;
vector<P>t;
LL f(LL n,LL p){
if(n<p) return 0;
return n/p+f(n/p,p);
}
LL g(LL n,LL p,LL k){
if(n==0) return 1;
LL t=pw(p,k),now=1;
// for(LL i=1;i<=t;i++){
// if(i%p!=0) now=now*i%t;
// }
LL d=lower_bound(id+1,id+1+len,t)-id;
now=now*a[d][t]%t;
now=pw(now,n/t,t);
now=now*g(n/p,p,k)%t;
// for(LL i=1;i<=n%t;i++){
// LL j=i;
// if(j%p!=0)
// now=now*j%t;
// }
now=now*a[d][n%t]%t;
assert(now);
return now;
}
LL C(LL x,LL y,LL p,LL k){
LL t=f(x,p)-f(y,p)-f(x-y,p);
if(t>=k) return 0;
k-=t;
LL tmp=pw(p,k);
assert(tmp);
LL ans=(g(x,p,k)*mod(g(y,p,k),1,tmp)%tmp*mod(g(x-y,p,k),1,tmp)%tmp)*pw(p,t);
return ans;
}
LL exlucas(LL n,LL m,LL p){
t.clear();
LL tmp=p;
for(LL i=2;i*i<=tmp;i++){
LL cnt=0;
while(tmp%i==0){
tmp/=i;
cnt++;
}
if(cnt) t.push_back(P(i,cnt));
}
if(tmp!=1) t.push_back(P(tmp,1));
vector<LL>a,b;
for(LL i=0;i<t.size();i++){
LL u=t[i].first,v=t[i].second;
a.push_back(C(n,m,u,v)),b.push_back(pw(u,v));
}
return excrt(a,b);
}
}
LL p[10]={1,10,100,1000,10000,100000,1000000,10000000,100000000,1000000000};
signed main(){
ios::sync_with_stdio(false); cin.tie(0); cout.tie(0);
exLucas::init();
long long a,b,k;
while(cin>>a>>b>>k){
LL mod=p[k]*2;
LL ans=exLucas::pw(2,a+b,mod);
for(LL i=b+1;i<=a-1;i++){
ans+=exLucas::exlucas(a+b,i,mod);
ans%=mod;
}
if(a==b) ans+=mod-exLucas::exlucas(a+b,b,mod);
ans%=mod;
ans/=2;
for(LL i=k-1;~i;i--) cout<<(int)(ans/p[i]%10); cout<<'\n';
}
return 0;
}
```
::::
### [P5325](https://www.luogu.com.cn/problem/P5325)
Min_25 筛的板子题，复杂度 $O(\frac {n^{\frac 34}}{\log n})$，但可以用杜教筛和 PN 筛在 $O(n^{\frac 23})$ 时间内通过。
令 $g(x)=x\varphi(x),t(x)=x^2$，有 $id\ast g=t$。
易证 $p$ 为质数时 $f(x)=g(x)$，可以 PN 筛，求 $G(x)$ 部分用杜教筛。
此时推导可得 $h(p^k)=(k-1)(p-1)p^k$，后略。
::::info[Code]
```cpp
#include<bits/stdc++.h>
using namespace std;
const long long B=15000000;
const long long mod=1000000007;
const long long d2=500000004;
const long long d6=166666668;
//g(x)=x\varphi(x),g*id=id2
long long G[B+5];//预处理g函数前缀和
long long p[(B>>3)+5],cnt;
long long s[B+5];
map<long long,long long>mp;
long long getG(long long h){//杜教筛
if(h<=B) return G[h];
if(mp[h]) return mp[h];
long long ans=h%mod*((h+1)%mod)%mod*((h+h+1)%mod)%mod*d6%mod;
long long l=2,r;
while(l*l<=h){
ans-=getG(h/l)*l%mod;
l++;
}
while(l<=h){
r=h/(h/l);
ans-=getG(h/l)*((r-l+1)%mod)%mod*((l+r)%mod)%mod*d2%mod;
l=r+1;
}
ans=(ans%mod+mod)%mod;
return mp[h]=ans;
}
struct edge{
long long num,p;
bool operator<(const edge &b)const{
return num<b.num;
}
bool operator==(const edge &b)const{
return num==b.num;
}
};
edge pn[1000005];
long long m;//Powerful Numbers
long long h[1000005];//i=PN处的h值
long long solve(long long n){//PN筛
//先筛Powerful Number
for(long long a=1;a*a<=n;a++){
for(long long b=1;a*a*b*b*b<=n;b++){
pn[++m]={a*a*b*b*b,s[max(a,b)]};
}
}
sort(pn+1,pn+m+1);
m=unique(pn+1,pn+m+1)-pn-1;
//再筛h
for(long long i=1;i<=m;i++){
long long k=pn[i].num;
if(k==1){
h[1]=1;
}
else{
long long p=pn[i].p;
assert(~p);
long long tmp=k;
long long w=0;
while(tmp%p==0){
tmp/=p;
w++;
}
assert(w>=2);
h[i]=h[lower_bound(pn+1,pn+m+1,(edge){tmp,0})-pn]*(w-1)*(p-1)*(k/tmp)%mod;
}
}
long long ans=0;
for(long long i=1;i<=m;i++){
ans+=h[i]*getG(n/pn[i].num);
ans%=mod;
}
return ans;
}
signed main(){
//线性预处理g函数前缀和
G[1]=1;
for(long long i=2;i<=B;i++){
if(!s[i]){
s[i]=i;
p[++cnt]=i;
}
for(long long j=1;j<=cnt&&i*p[j]<=B;j++){
s[i*p[j]]=p[j];
if(i%p[j]==0) break;
}
}
for(long long i=1;i<=cnt;i++){
long long d=p[i];
for(long long j=d;j<=B;j*=d){
G[j]=j/d*(d-1);
}
}
for(long long i=2;i<=B;i++) if(!G[i]){
if(i/s[i]%s[i]==0) G[i]=G[i/s[i]]*s[i];
else G[i]=G[i/s[i]]*G[s[i]];
}
for(long long i=1;i<=B;i++) G[i]*=i;
for(long long i=1;i<=B;i++) G[i]%=mod;
for(long long i=1;i<=B;i++) G[i]+=G[i-1];
for(long long i=1;i<=B;i++) G[i]%=mod;
long long n;
cin>>n;
cout<<solve(n);
return 0;
}
```
::::