浅谈 OI 中的数论

一些记号

a|b

表示

b

为

a

的倍数。

\gcd(a,b)

表示

a,b

的最大公约数。

[x=k]

表示函数

f(x)=\begin{cases}1&x=k\\0&x\ne k\end{cases}

素数的基本知识

素数定义

因数个数为

2

的数为素数（质数）。

素数判定

朴素判定

O(\sqrt n)

，应该都会。
也可以先打素数表，单次查询可降至

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

。
模板题代码实现（long long 范围内

100\%

正确）：

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。
exLucas 模板题代码：

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}

id(x)=x

表示函数值恒为自变量本身的函数。

\varphi(x)

表示

1

到

x

中与

x

互质的数的个数。

举例

以防你看不懂，下面列举部分函数值：

表格

$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)

。模板题代码：

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

题面

计算：

\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)

。
本题轻微卡常。

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

题面

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}

根号分治，后略。

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

记小 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

。

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

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

，后略。

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;
}

文章操作

一些记号

素数的基本知识

素数定义

素数判定

朴素判定

Miller-Rabin 素数测试

Lucas 定理、exLucas 定理

Lucas

exLucas

一些数论函数

定义

举例

性质

莫比乌斯反演

狄利克雷卷积定义

重要性质

筛法

素数筛

杜教筛

PN 筛

例题选讲

P2257

题面

解法

P4240

题面

解法

P3726

P5325

相关推荐

评论

浅谈 OI 中的数论