模板：数论

最新推荐文章于 2025-01-25 18:45:25 发布

wu-kan

最新推荐文章于 2025-01-25 18:45:25 发布

阅读量754

点赞数

分类专栏： acm 模板文章标签： acm 模板数论

本文及后续更新已转到https://wu-kan.cn 本文为博主原创文章，转载需要注明来源。

本文链接：https://blog.csdn.net/w_weilan/article/details/78867002

版权

acm 同时被 2 个专栏收录

54 篇文章

订阅专栏

模板

16 篇文章

订阅专栏

这篇博客探讨了数论在ACM竞赛中的重要应用，包括辗转相除法、同余系运算、欧拉筛、素数因子分解、快速傅里叶变换等。文章详细阐述了裴蜀定理、欧拉函数及其性质，并介绍了如何利用这些概念解决乘法逆元、模方程等问题。此外，还提到了莫比乌斯反演、快速数论变换和Pell方程等高级数论技巧。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

数论

辗转相除法与同余系

裴蜀定理：对任何 $a,b\in Z$ 和它们的最大公约数 $d$ ，关于未知数 $x, y$ 的线性不定方程（称为裴蜀等式）： $a x + b y = c$ 当仅当 $d\mid c$ ，可知有无穷多解。特别地， $a x + b y = d$ 一定有解。
推论： $a, b$ 互质的充要条件是 $a x + b y = 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^{\phi(m)}\equiv1\pmod m$ ，x的逆是 $x^{\phi(m)-1}$ 。特别地，m为素数时 $\phi(m)=m-1$ ，此时x的逆就是pow(x,m-2,m)。
log函数：m为素数时求解模方程 $a^x\equiv b\pmod m$ 。设P为质数，G为P的原根，则 $x^y\equiv b\pmod P$ 等价于 $y\ ind\ x\equiv b\pmod{P-1}$ ，其中 $G\ ind\ x\equiv x\pmod P$ 。

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

欧拉筛

欧拉函数 $\phi(n)$ 是小于n的正整数中与n互素的数的数目。特别地，规定 $\phi(1)=1$ ，易知n>2时都为偶数。
欧拉函数是积性函数，即对任意素数 $p, q$ 满足下列关系： $\phi(pq)=\phi(p)\phi(q)=(p-1)(q-1)$ 对任何两个互质的正整数 $m(m\geq2)$ 有欧拉定理： $x^{\phi(m)}\equiv1\pmod m$ 当m为素数p时，此式变为费马小定理： $x^{p-1}\equiv1\pmod p$ 利用欧拉函数和它本身不同质因数的关系，用筛法 $O (N)$ 预处理某个范围内所有数的欧拉函数值，并求出素数表。同时，利用计算欧拉函数过程中求出的最小素因子m，可以实现 $O (l o g N)$ 的素因数分解。
更新：增加同时求莫比乌斯函数 $\mu(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];
			}
		}
	}
};

常见数论函数变换

$\sum_{d|n}\mu(d)=[n=1]$
$\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)$

莫比乌斯反演

若 $f(n)=\sum_{d|n}g(d)$ ，则 $g(n)=\sum_{d|n}\mu(d)f(\frac{n}{d})$
若 $f(n)=\sum_{i=1}^nt(i)g(\frac{n}{i})$ ，则 $g(n)=\sum_{i=1}^n\mu(i)t(i)f(\frac{n}{i})$ （此时代 $t (i) = [g c d (n, i) > 1]$ 可得上式）
举例（其中 $T = k d$ ）： $\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}$
$\varphi(T)$ 可以使用线性筛预处理处理，我们就可以枚举 $T$ 求上式了，时间复杂度 $O (n)$ 。多组数据 $n, m$ 询问上式，时间复杂度就变成了 $O (T n)$ 。事实上， $\lfloor\frac{n}{T}\rfloor$ 是不会轻易变化的，是过了连续的一段后才发生变化的，那么我们就可以计算出这一段的结束位置，对 $\varphi$ 函数作前缀和，就可以直接分块了，这样的时间复杂度是 $O(T\sqrt{n})$ 的。

前缀和

欧拉函数前缀和 $S_\phi(n)=\frac{(n+1)n}2-\sum_{d=1}^nS_\phi(\frac{n}{d})$
莫比乌斯函数前缀和 $S_\mu(n)=1-\sum_{d=1}^nS_\mu(\frac{n}{d})$

PollardRho大数素因子分解

时间复杂度 $O(N^{1/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^nc+1$ 的费马素数，记其原根为g，则旋转因子为 $g^{(m-1)/n}$ ，满足 $g^{m-1}=1$ 且 $2^n|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-Dy^2=1$ （D为任意正整数）的方程称为佩尔方程，必有最小正整数解 $x_0,y_0)$ ，用 $x_n=x_0x_{n-1}+Dy_0y_{n-1},y_n=y_0x_{n-1}+x_0y_{n-1}$ 可递推方程的第n小整数解（可用矩阵快速幂求），同时还有
$2x_0x_n=x_{n-1}+x_{n+1},2x_0y_n=y_{n-1}+y_{n+1}$

Jacobi’s Four Square Theorem

设 $a^2 + b^2 + c^2 + d^2 = n$ 的自然数解个数为 r4(n)，d(n) 为 n 的约数和，由 Jacobi’s Four Square Theorem 可知，若 n 是奇数，则 r4(n) = 8d(n)，否则 r4(n) = 24d(k)，k 是 n 去除所有 2 后的结果。