这篇博客探讨了数论在ACM竞赛中的重要应用，包括辗转相除法、同余系运算、欧拉筛、素数因子分解、快速傅里叶变换等。文章详细阐述了裴蜀定理、欧拉函数及其性质，并介绍了如何利用这些概念解决乘法逆元、模方程等问题。此外，还提到了莫比乌斯反演、快速数论变换和Pell方程等高级数论技巧。
数论
辗转相除法与同余系
裴蜀定理：对任何
a
,
b
∈
Z
a,b\in Z
a,b∈Z和它们的最大公约数
d
d
d，关于未知数
x
,
y
x,y
x,y的线性不定方程（称为裴蜀等式）：
a
x
+
b
y
=
c
ax+by=c
ax+by=c当仅当
d
∣
c
d\mid c
d∣c，可知有无穷多解。特别地，
a
x
+
b
y
=
d
ax+by=d
ax+by=d一定有解。
推论：
a
,
b
a,b
a,b互质的充要条件是
a
x
+
b
y
=
1
ax+by=1
ax+by=1有整数解。
ll gcd(ll a,ll b)//a、b的最大公约数
{
return b?gcd(b,a%b):a;
}
ll lcm(ll a,ll b)//a、b的最小公倍数
{
return a/gcd(a,b)*b;
}
ll gcd(ll a,ll b,ll &x,ll &y)//扩展欧几里得，引用返回a*x+b*y=gcd(a,b)绝对值之和最小的解
{
if(!a)return x=0,y=1,b;
ll d=gcd(b%a,a,y,x);
return x-=y*(b/a),d;
}
同余系运算
求乘法逆元的另外一种方法是用欧拉定理
x
ϕ
(
m
)
≡
1
(
m
o
d
m
)
x^{\phi(m)}\equiv1\pmod m
xϕ(m)≡1(modm)，x的逆是
x
ϕ
(
m
)
−
1
x^{\phi(m)-1}
xϕ(m)−1。特别地，m为素数时
ϕ
(
m
)
=
m
−
1
\phi(m)=m-1
ϕ(m)=m−1，此时x的逆就是pow(x,m-2,m)。
log函数：m为素数时求解模方程
a
x
≡
b
(
m
o
d
m
)
a^x\equiv b\pmod m
ax≡b(modm)。设P为质数，G为P的原根，则
x
y
≡
b
(
m
o
d
P
)
x^y\equiv b\pmod P
xy≡b(modP)等价于
y
i
n
d
x
≡
b
(
m
o
d
P
−
1
)
y\ ind\ x\equiv b\pmod{P-1}
y ind x≡b(modP−1)，其中
G
i
n
d
x
≡
x
(
m
o
d
P
)
G\ ind\ x\equiv x\pmod P
G ind x≡x(modP)。
ll add(ll a,ll b,ll m)
{
if(a+=b,a<0)a=m-(-a%m);//负数取模和编译器有关，这样才是需要的同余系运算
return a%m;
}
ll mul(ll a,ll b,ll m)//根据a*b是否爆ll替换a*b%m
{
ll r=0;
for(a%=m; b; b>>=1,a=add(a,a,m))
if(b&1)r=add(r,a,m);
return r;
}
ll pow(ll a,ll b,ll m)
{
ll r=1;
for(a%=m; b; b>>=1,a=mul(a,a,m))
if(b&1)r=mul(r,a,m);
return r;
}
ll sub(ll a,ll b,ll m)
{
return add(a,-b,m);
}
ll inv(ll a,ll m)//模m下a的乘法逆元，不存在返回-1（m为素数时a不为0必有逆元）
{//return pow(a,phi(m)-1,m);
ll x,y,d=gcd(a,m,x,y);
return d==1?(x+m)%m:-1;
}
ll div(ll a,ll b,ll m)
{
return mul(a,inv(b,m),m);
}
ll log(ll a,ll b,ll m)
{
ll n=ceil(sqrt(m+0.5));
map<ll,ll> x;
for(ll i=0,e=1; i<n; e=mul(e,a,m),++i)
if(!x.count(e))x[e]=i;
for(ll i=0,v=inv(pow(a,n,m),m); i<n; ++i,b=mul(b,v,m))
if(x.count(b))return i*n+x[b];
return -1;
}
同余方程
void sol(ll a,ll b,ll n,vector<ll> &ans)//返回ax=b(mod n)循环节内所有解
{
ll x,y,d=gcd(a,n,x,y);
if(ans.clear(),b%d)return;
ans.push_back((b/d)%(n/d)*(x=(x%n+n)%n));
for(int i=1; i<d; ++i)
ans.push_back((ans[0]+i*n/d)%n);
}
同余方程组
ll sol(const vector<pair<ll,ll> > &v)//x%v[i].first==v[i].second，不存在返回-1
{
ll m=v[0].first,r=v[0].second,c,d,x,y,z;
for(int i=1; i<v.size(); ++i)
{
if(c=v[i].second-r,d=gcd(m,v[i].first,x,y),c%d)
return -1;
gcd(m/d,z=v[i].first/d,x,y),r+=c/d*x%z*m,r%=m*=z;
}
return r<0?r+m:r;
}
欧拉筛
欧拉函数
ϕ
(
n
)
\phi(n)
ϕ(n)是小于n的正整数中与n互素的数的数目。特别地，规定
ϕ
(
1
)
=
1
\phi(1)=1
ϕ(1)=1，易知n>2时都为偶数。
欧拉函数是积性函数，即对任意素数
p
,
q
p,q
p,q满足下列关系：
ϕ
(
p
q
)
=
ϕ
(
p
)
ϕ
(
q
)
=
(
p
−
1
)
(
q
−
1
)
\phi(pq)=\phi(p)\phi(q)=(p-1)(q-1)
ϕ(pq)=ϕ(p)ϕ(q)=(p−1)(q−1)对任何两个互质的正整数
x
,
m
(
m
≥
2
)
x, m(m\geq2)
x,m(m≥2)有欧拉定理：
x
ϕ
(
m
)
≡
1
(
m
o
d
m
)
x^{\phi(m)}\equiv1\pmod m
xϕ(m)≡1(modm)当m为素数p时，此式变为费马小定理：
x
p
−
1
≡
1
(
m
o
d
p
)
x^{p-1}\equiv1\pmod p
xp−1≡1(modp)利用欧拉函数和它本身不同质因数的关系，用筛法
O
(
N
)
O(N)
O(N)预处理某个范围内所有数的欧拉函数值，并求出素数表。同时，利用计算欧拉函数过程中求出的最小素因子m，可以实现
O
(
l
o
g
N
)
O(log N)
O(logN)的素因数分解。
更新：增加同时求莫比乌斯函数
μ
(
n
)
\mu(n)
μ(n)的代码，存在mu中
struct EulerSieve
{
vector<int> p,m,phi,mu;//素数序列，最小素因子，欧拉函数，莫比乌斯函数
EulerSieve(int N):m(N,0),phi(N,0),mu(N,0)
{
phi[1]=mu[1]=1;//m[1]=0
for(long long i=2,k; i<N; ++i)//防i*p[j]爆int
{
if(!m[i])p.push_back(m[i]=i),phi[i]=i-1,mu[i]=-1;//i是素数
for(int j=0; j<p.size()&&(k=i*p[j])<N; ++j)
{
phi[k]=phi[i]*p[j];
if((m[k]=p[j])==m[i])
{
mu[k]=0;
break;
}
phi[k]-=phi[i];
mu[k]=-mu[i];
}
}
}
};
常见数论函数变换
∑
d
∣
n
μ
(
d
)
=
[
n
=
1
]
\sum_{d|n}\mu(d)=[n=1]
∑d∣nμ(d)=[n=1]
ϕ
(
n
)
=
∑
i
=
1
n
[
gcd
(
i
,
n
)
=
1
]
=
∑
i
=
1
n
∑
k
∣
i
,
k
∣
n
μ
(
k
)
=
∑
k
∣
n
n
k
μ
(
k
)
\phi(n)=\sum_{i=1}^n[\gcd(i,n)=1]=\sum_{i=1}^n\sum_{k\mid i,k\mid n}\mu(k)=\sum_{k\mid n}\frac nk\mu(k)
ϕ(n)=∑i=1n[gcd(i,n)=1]=∑i=1n∑k∣i,k∣nμ(k)=∑k∣nknμ(k)
莫比乌斯反演
若
f
(
n
)
=
∑
d
∣
n
g
(
d
)
f(n)=\sum_{d|n}g(d)
f(n)=∑d∣ng(d)，则
g
(
n
)
=
∑
d
∣
n
μ
(
d
)
f
(
n
d
)
g(n)=\sum_{d|n}\mu(d)f(\frac{n}{d})
g(n)=∑d∣nμ(d)f(dn)
若
f
(
n
)
=
∑
i
=
1
n
t
(
i
)
g
(
n
i
)
f(n)=\sum_{i=1}^nt(i)g(\frac{n}{i})
f(n)=∑i=1nt(i)g(in)，则
g
(
n
)
=
∑
i
=
1
n
μ
(
i
)
t
(
i
)
f
(
n
i
)
g(n)=\sum_{i=1}^n\mu(i)t(i)f(\frac{n}{i})
g(n)=∑i=1nμ(i)t(i)f(in)（此时代
t
(
i
)
=
[
g
c
d
(
n
,
i
)
>
1
]
t(i)=[gcd(n,i)>1]
t(i)=[gcd(n,i)>1]可得上式）
举例（其中
T
=
k
d
T=kd
T=kd）：
∑
i
=
1
n
∑
j
=
1
m
gcd
(
i
,
j
)
=
∑
d
d
∑
i
=
1
n
∑
j
=
1
m
[
gcd
(
i
,
j
)
=
d
]
=
∑
d
d
∑
i
=
1
⌊
n
d
⌋
∑
j
=
1
⌊
m
d
⌋
[
gcd
(
i
,
j
)
=
1
]
=
∑
d
d
∑
i
=
1
⌊
n
d
⌋
∑
j
=
1
⌊
m
d
⌋
∑
k
∣
i
,
k
∣
j
μ
(
k
)
=
∑
d
d
∑
k
μ
(
k
)
∑
k
∣
i
⌊
n
d
⌋
∑
k
∣
j
⌊
m
d
⌋
=
∑
d
d
∑
k
μ
(
k
)
⌊
n
k
d
⌋
⌊
m
k
d
⌋
=
∑
T
⌊
n
T
⌋
⌊
m
T
⌋
∑
k
∣
T
T
k
μ
(
k
)
=
∑
T
⌊
n
T
⌋
⌊
m
T
⌋
φ
(
T
)
\begin{aligned} \sum_{i=1}^n\sum_{j=1}^m\gcd(i,j)&=\sum_d d\sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)=d] \\ &=\sum_{d}d\sum_{i=1}^{\lfloor\frac nd\rfloor}\sum_{j=1}^{\lfloor\frac md\rfloor}[\gcd(i,j)=1] \\ &=\sum_{d}d\sum_{i=1}^{\lfloor\frac nd\rfloor}\sum_{j=1}^{\lfloor\frac md\rfloor}\sum_{k\mid i,k\mid j}\mu(k) \\ &=\sum_d d\sum_k\mu(k)\sum_{k\mid i}^{\lfloor\frac nd\rfloor}\sum_{k\mid j}^{\lfloor\frac md\rfloor} \\ &=\sum_{d}d\sum_k\mu(k)\lfloor\frac n{kd}\rfloor\lfloor\frac m{kd}\rfloor\\&=\sum_{T}\lfloor\frac nT\rfloor\lfloor\frac mT\rfloor\sum_{k\mid T}\frac Tk\mu(k) \\ &=\sum_{T}\lfloor\frac nT\rfloor\lfloor\frac mT\rfloor\varphi(T) \end{aligned}
i=1∑nj=1∑mgcd(i,j)=d∑di=1∑nj=1∑m[gcd(i,j)=d]=d∑di=1∑⌊dn⌋j=1∑⌊dm⌋[gcd(i,j)=1]=d∑di=1∑⌊dn⌋j=1∑⌊dm⌋k∣i,k∣j∑μ(k)=d∑dk∑μ(k)k∣i∑⌊dn⌋k∣j∑⌊dm⌋=d∑dk∑μ(k)⌊kdn⌋⌊kdm⌋=T∑⌊Tn⌋⌊Tm⌋k∣T∑kTμ(k)=T∑⌊Tn⌋⌊Tm⌋φ(T)
φ
(
T
)
\varphi(T)
φ(T)可以使用线性筛预处理处理，我们就可以枚举
T
T
T求上式了，时间复杂度
O
(
n
)
O(n)
O(n)。多组数据
n
,
m
n,m
n,m询问上式，时间复杂度就变成了
O
(
T
n
)
O(Tn)
O(Tn)。事实上，
⌊
n
T
⌋
\lfloor\frac{n}{T}\rfloor
⌊Tn⌋是不会轻易变化的，是过了连续的一段后才发生变化的，那么我们就可以计算出这一段的结束位置，对
φ
\varphi
φ函数作前缀和，就可以直接分块了，这样的时间复杂度是
O
(
T
n
)
O(T\sqrt{n})
O(Tn)的。
前缀和
欧拉函数前缀和
S
ϕ
(
n
)
=
(
n
+
1
)
n
2
−
∑
d
=
1
n
S
ϕ
(
n
d
)
S_\phi(n)=\frac{(n+1)n}2-\sum_{d=1}^nS_\phi(\frac{n}{d})
Sϕ(n)=2(n+1)n−∑d=1nSϕ(dn)
莫比乌斯函数前缀和
S
μ
(
n
)
=
1
−
∑
d
=
1
n
S
μ
(
n
d
)
S_\mu(n)=1-\sum_{d=1}^nS_\mu(\frac{n}{d})
Sμ(n)=1−∑d=1nSμ(dn)
PollardRho大数素因子分解
时间复杂度 O ( N 1 / 4 ) O(N^{1/4}) O(N1/4)，数据多的时候可考虑欧拉筛优化。
struct PollardRho
{
bool isPrime(ll n,int S=12)//MillerRabin素数测试，S为测试次数，用前S个素数测一遍，S=12可保证unsigned long long范围内无错；n<2请特判
{
ll d,u,t,p[]= {2,3,5,7,11,13,17,19,23,29,31,37};
for(d=n-1; !(d&1);)d>>=1;//未对0，1做特判！
for(int i=0; i<S; ++i)
{
if(n%p[i]==0)return n==p[i];
for(u=d,t=pow(p[i],u,n); u!=n-1&&t!=n-1&&t!=1; u<<=1)
t=mul(t,t,n);
if(t!=n-1&&!(u&1))return 0;
}
return 1;
}
void fac(ll n,vector<ll> &factor)
{
if(isPrime(n))return factor.push_back(n);
for(ll c=1;; ++c)
for(ll i=0,k=1,x=rand()%(n-1)+1,y,p;;)
{
if(++i==k)y=x,k<<=1;
if(x=(mul(x,x,n)+c)%n,p=gcd(abs(x-y),n),p==n)break;
if(p>1)return fac(p,factor),fac(n/p,factor);
}
}
};
快速傅里叶变换（FFT）
struct Rader:vector<int>
{
Rader(int n):vector<int>(1<<int(ceil(log2(n))))
{
for(int i=at(0)=0; i<size(); ++i)
if(at(i)=at(i>>1)>>1,i&1)
at(i)+=size()>>1;
}
};
struct FFT:Rader
{
vector<complex<lf> > w;
FFT(int n):Rader(n),w(size(),polar(1.0,2*PI/size()))
{
w[0]=1;
for(int i=2; i<size(); ++i)w[i]*=w[i-1];
}
vector<complex<lf> > fft(const vector<complex<lf> > &a)const
{
vector<complex<lf> > x(size());
for(int i=0; i<a.size(); ++i)x[at(i)]=a[i];
for(int i=1; i<size(); i<<=1)
for(int j=0; j<i; ++j)
for(int k=j; k<size(); k+=i<<1)
{
complex<lf> &l=x[k],&r=x[k+i],t=w[size()/(i<<1)*j]*r;
r=l-t,l+=t;
}
return x;
}
vector<ll> ask(const vector<ll> &a,const vector<ll> &b)const
{
vector<complex<lf> > xa(a.begin(),a.end()),xb(b.begin(),b.end());
xa=fft(xa),xb=fft(xb);
for(int i=0; i<size(); ++i)xa[i]*=xb[i];
vector<ll> ans(size());
xa=fft(xa),ans[0]=xa[0].real()/size()+0.5;
for(int i=1; i<size(); ++i)ans[i]=xa[size()-i].real()/size()+0.5;
return ans;
}
};
优化高精度乘法
常规算法。
Wint& operator*=(Wint &a,const Wint &b)//fft优化乘法，注意double仅15位有效数字，调小Wint::width不超过2，计算自2*log2(base)+2*log2(len)<53
{
vector<ll> ax(a.begin(),a.end()),bx(b.begin(),b.end());
ax=FFT(a.size()+b.size()).ask(ax,bx);
for(int i=1; i<ax.size(); ++i)
ax[i]+=ax[i-1]/a.base,ax[i-1]%=a.base;
return a.assign(ax.begin(),ax.end()),a.trim(0);
}
快速数论变换（FNTT）
原理和FFT相同，解决特殊情况下FFT的浮点误差，并且可以在同余系进行变换。
对于形如
m
=
2
n
c
+
1
m=2^nc+1
m=2nc+1的费马素数，记其原根为g，则旋转因子为
g
(
m
−
1
)
/
n
g^{(m-1)/n}
g(m−1)/n，满足
g
m
−
1
=
1
g^{m-1}=1
gm−1=1且
2
n
∣
m
−
1
2^n|m-1
2n∣m−1。
常见素数的原根。
struct FNTT:Rader
{
ll M,G;
vector<ll> w;
FNTT(int N,ll M,ll G):Rader(N),M(M),G(G),w(size(),pow(G,(M-1)/size(),M))//M费马素数，G原根
{
for(int i=w[0]=1; i<size(); ++i)w[i]=mul(w[i],w[i-1],M);
}
vector<ll> fntt(const vector<ll> &a)const
{
vector<ll> x(size());
for(int i=0; i<a.size(); ++i)x[at(i)]=a[i];
for(int i=1,j; i<size(); i<<=1)
for(int j=0; j<i; ++j)
for(int k=j; k<size(); k+=i<<1)
{
ll t=mul(w[size()/(i<<1)*j],x[k+i],M);
if(x[k+i]=x[k]-t,x[k+i]<0)x[k+i]+=M;
if(x[k]+=t,x[k]>=M)x[k]-=M;
}
return x;
}
vector<ll> ask(vector<ll> a,vector<ll> b)const
{
a=fntt(a),b=fntt(b);
for(int i=0; i<size(); ++i)a[i]=mul(a[i],b[i],M);
a=fntt(a),reverse(a.begin()+1,a.end());
ll u=inv(size(),M);
for(int i=0; i<size(); ++i)a[i]=mul(a[i],u,M);
return a;
}
};
优化高精度乘法
常规算法。
Wint& operator*=(Wint &a,const Wint &b)//ntt优化，Wint::width不超过2
{
vector<ll> ax(a.begin(),a.end()),bx(b.begin(),b.end());
ax=FNTT(a.size()+b.size(),(7<<26)+1,3).ask(ax,bx);
for(int i=1; i<ax.size(); ++i)
ax[i]+=ax[i-1]/a.base,ax[i-1]%=a.base;
return a.assign(ax.begin(),ax.end()),a.trim(0);
}
快速沃尔什变换（FWT）
再给一个二分的代码便于理解。
void fwt(vector<ll> &x,void f(ll &l,ll &r))
{
for(int i=1; i<x.size(); i<<=1)
for(int j=0; j<i; ++j)
for(int k=j; k<x.size(); k+=i<<1)
f(x[k],x[k+i]);
}
void fwt(ll *b,ll *e,void f(ll &l,ll &r))
{
if(e-b<2)return;
ll *m=b+(e-b)/2;
fwt(b,m,f),fwt(m,e,f);
while(m<e)f(*(b++),*(m++));
}
XOR
可能要在同余系中进行运算，下面代码需要修改。
void tf(ll &l,ll &r)
{
ll tl=l+r,tr=l-r;
l=tl,r=tr;
}
void utf(ll &l,ll &r)
{
tf(l,r),l>>=1,r>>=1;
}
AND
void tf(ll &l,ll &r)
{
l+=r;
}
void utf(ll &l,ll &r)
{
l-=r;
}
OR
void tf(ll &l,ll &r)
{
r+=l;
}
void utf(ll &l,ll &r)
{
r-=l;
}
XNOR,NAND,NOR
直接用异或运算、与运算、或运算的方法求出来，然后将互反的两位交换即可。
Pell方程
形如
x
2
−
D
y
2
=
1
x^2-Dy^2=1
x2−Dy2=1（D为任意正整数）的方程称为佩尔方程，必有最小正整数解
(
x
0
,
y
0
)
(x_0,y_0)
(x0,y0)，用
x
n
=
x
0
x
n
−
1
+
D
y
0
y
n
−
1
,
y
n
=
y
0
x
n
−
1
+
x
0
y
n
−
1
x_n=x_0x_{n-1}+Dy_0y_{n-1},y_n=y_0x_{n-1}+x_0y_{n-1}
xn=x0xn−1+Dy0yn−1,yn=y0xn−1+x0yn−1可递推方程的第n小整数解（可用矩阵快速幂求），同时还有
2
x
0
x
n
=
x
n
−
1
+
x
n
+
1
,
2
x
0
y
n
=
y
n
−
1
+
y
n
+
1
2x_0x_n=x_{n-1}+x_{n+1},2x_0y_n=y_{n-1}+y_{n+1}
2x0xn=xn−1+xn+1,2x0yn=yn−1+yn+1
Jacobi’s Four Square Theorem
设 a 2 + b 2 + c 2 + d 2 = n a^2 + b^2 + c^2 + d^2 = n a2+b2+c2+d2=n 的自然数解个数为 r4(n)，d(n) 为 n 的约数和，由 Jacobi’s Four Square Theorem 可知，若 n 是奇数，则 r4(n) = 8d(n)，否则 r4(n) = 24d(k)，k 是 n 去除所有 2 后的结果。
扫一扫
466
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
