数论
Tags:Noip前的大抱佛脚
知识点
Exgcd
(O(logn))求解(Ax+By=C)的问题
1、若(C\%gcd(A,B)!=0)则无解
2、(Gcd=gcd(A,B);A/=Gcd,B/=Gcd,C/=Gcd)
3、代入下面代码求(Ax+By=1)
4、(x*C),得到一组特解
5、通解为(egin{cases}x=x_0+k*B \y=y_0+k*Aend{cases})
void Exgcd(ll a,ll b,ll &x,ll &y)
{
if(!b){x=1;y=0;return;}
Exgcd(b,a%b,y,x);y-=a/b*x;
}
逆元
在(gcd(A,P)==1)时,(A)在模(P)意义下存在逆元(证明可以符合Exgcd有解的证明),其余情况不存在逆元
通常(P)为质数时就是(A^{p-2})作为逆元(费马小定理)
-
欧拉定理(费马小定理) (x^{-1}equiv x^{p-2}(mod p)),要求(p)为质数
-
解方程(Exgcd)(Px+Ay=1,A=frac{1}{y}(mod P))
-
线性递推逆元:
令(a=lfloorfrac{p}{i} floor,b=p\%i),则有(ai+b=p,iequiv-frac{b}{a}(mod p),frac{1}{i}equiv-frac{a}{b}(mod p))
由于(b<i),所以(b)的逆元已经求出,直接可以得到:
inv[i]=p-(p/i*inv[p%i])%p;
gcd
具有一些奇妙的性质,如可合并
- 往数集中加入一个数,要么gcd不变,要么至少变为原来的(frac{1}{2})(所以可以用来分块了(同样的xor也每次只会增加log次))
欧拉函数(varphi(x))
表示小于x且与x互质的数的个数
计算公式
降幂公式
- (x,p)互质:(x^kequiv x^{k\%varphi(p)} (mod p))
- (x,p)不互质:(x^kequiv egin{cases}x^{k} (mod p),klevarphi(p) \x^{k\%varphi(p)+varphi(p)} (mod p),k>varphi(p)end{cases})
CRT&EXCRT
(O(nlogn))求解一系列同余方程,如
代码如下(未判无解,无解即exgcd无解)
ll gcd(ll a,ll b) {return !b?a:gcd(b,a%b);}
void Exgcd(ll a,ll b,ll &x,ll &y)
{
if(!b) {x=1;y=0;return;}
Exgcd(b,a%b,y,x);y-=a/b*x;
}
int EXCRT()
{
for(int i=2;i<=n;i++)
{
ll g=gcd(P[i-1],P[i]),C=(A[i]-A[i-1])/g,x,y;
Exgcd(P[i-1]/g,P[i]/g,x,y);
P[i]=P[i-1]*P[i]/g;
A[i]=A[i-1]+P[i-1]*x*C;
A[i]=(A[i]%P[i]+P[i])%P[i];
}
return A[n];
}
方法是每次合并两个方程,手推式子(P[i]x+A[i]=P[i+1]y+A[i+1] -> P[i]x+P[i+1]y=A[i+1]-A[i])
可以发现模数为原来的(lcm)(先除以(gcd)再乘,很容易爆(long long)),余数为原来余数加上模数的(x)倍
注意经常x算出来是负数所以时刻(+mod)\%mod !)
有的时候题目并没有那么裸
- (x)前带系数(屠龙勇士)
(Axequiv B(mod p) -> Ax+pk=B(当且仅当B\%gcd(A,p)==0时有解) -> x=x_0+tp -> xequiv x_0(mod p))
- CRT合并答案
对于一些公式,只适用于模数是质数的情况,而题目要求的模数要是任意正数,所以可以求得
然后用CRT合并答案,求得(Ansequiv OUTANS (mod Mod))
但是就所见到的题目来说(任意模数NTT),只有在答案不太大的时候适用,例如该题答案不会超过(10^9*10^9*10^5=10^{23}),方法是选取三个乘积大于(10^{23})的NTT模数,得到
令(M=P_1*P_2)于是(Ans=kM+A_1=k_3P_3+A_2),即(kM+A_1equiv A_2(mod P_3))
由于(Ansle10^{23}<P_1P_2P_3=MP_3),所以(k<P_3),根据此同余方程可以求得(k)在(mod P_3)的解也就是(k)的实际值
然后就可以带入求得答案(Ans=kM+A_1)了,不过会爆(long long),可以采用这个
ll mul(ll x,ll y,ll p) {return (x*y-(ll)((long double)x/p*y+0.5)*p+p)%p;}
BSGS&EXBSGS
快速((O(sqrt n)))求解(A^x=B(mod P))的(x)的解,板子题:[SDOI2013]随机数生成器
普通的BSGS要求P是质数,拓展的可以不用是质数
通过降幂公式,能够知道(xle P-1),令(M=sqrt P),则(x=iM+j),可以枚举(i)的值,得到(t=A^{(i+1)M}),查表看是否存在(T=A^{M-j}B)满足条件,得到解就返回,这样就能保证解是最小的了
int BSGS(int A,int B,int P)
{
if(!A%P) return -1;
int M=sqrt(P)+1;Hash.reset();
for(int i=0,t=B;i<M;i++,t=1ll*t*A%P) Hash.Add(t%Mo,t,i);//存入B*A^i的哈希值
for(int i=1,bs=ksm(A,M,P),t=bs;i<=M;i++,t=1ll*t*bs%P)//t=A^{(i+1)M}
if(Hash.Query(t)!=-1) return i*M-Hash.Query(t);
return -1;
}
FFT/NTT/MTT/FWT
板子见最下方
原根相关
定义:
原根类似于FFT的单位复根,使得能够进行取模操作
NTT模数是指有原根且是(r×2^k+1)的形式的模数,如(998244353(3))、(1004535809(3))
存在性及判定:
一个数有原根当且仅当它为(2,4,p,2p,p^r)(p为奇质数)
(g^1,g^2...g^{phi(p)})在(mod p)下各不相同,且(g^{phi (p)}=1),当p为质数时,也就是说(g^1...g^{p-1})分别对应(1...p-1)
在(p)为质数时,令(p-1=p_1^{a1}p_2^{a_2}...p_n^{a_n}),若存在(g^{frac{p-1}{p_i}}=1)则(g)不是原根,否则是原根(证明)
注意点
- 两个多项式最高次项分别为(l1,l2),那么需要开的长度(包括0)是(>l1+l2)的第一个形如(2^k)的数
组合公式
-
从((0,0))走到((n,m)),不碰到直线(y=x+b)的方案数:(C(n+m,n)-C(n+m,n-b))
意义为沿着直线翻折,越过直线的路径对称域从翻着顶点到终点的路径
斯特林数
第一类斯特林数
第一类斯特林数:(n)个人坐(m)张圆桌的方案数(人不同,圆桌相同,n个元素构成m个圆排列)
第二类斯特林数
第二类斯特林数:(n)个球放入(m)个盒子的方案数(球不同,盒子相同,n个元素分成m个非空集合)
卡塔兰数
常用数学公式
- 平方和的前缀和 (sum_{i=1}^{n}i^2=frac{n(n+1)(2n+1)}{6})
- 立方和的前缀和 (sum_{i=1}^{n}i^3=(frac{n(n+1)}{2})^2)
技巧经验
容斥
-
从((0,0))到((n,m))且不经过一些禁点的方案数(BZOJ两双手)
(f[i]=sum_{j=0}^{j<i}-f[j]*way(j,i),f[0]=-1,way(j,i))表示从j点到i点的方案数,按照从((0,0))走到该点的所需步数排序
组合计数
- 可以画成二维平面上的点去移动,辅助理解、推式子
区间筛
-
区间筛质数个数
luogu1835素数密度:筛出长度为(10^6)的一段区间内的质数个数,左右端点在int以内
方法:枚举到(sqrt r)然后用(j=(l-1)/i*i+i)暴力在这个区间内筛掉合数
-
区间筛(varphi())之和
luogu3601签到题:筛出长度为(10^6)的一段区间的(varphi(i))之和,左右端点在(10^{12})以内
方法:枚举到(sqrt r)然后(ans[i]=ans[i]/p*(p-1))
-
区间筛(mu())之和
luogu3653小清新数学题:筛出长度为(10^6)的一段区间的(mu(i))之和,左右端点在(10^{18})以内
方法:枚举到(10^6)然后用质数暴力搞,除剩下的数如果不是(1)那么只有三种情况:
-
质数的平方:用(sqrt)判掉
-
一个质数
-
两个质数之积
上述两种可以用Miller—Rabin素性判断判掉,不过更方便的是(check):(2^{p-1}\%p==1)、(3^{p-1}\%p==1)。这题判两次可以过,当然是有可能不准确的。。
博弈
- Nim游戏:异或和为0则先手必败,否则必胜
- 威佐夫博弈:约定(a<b),(a=frac{sqrt 5+1}{2}(a-b))时是奇异局势,先手必败。
有趣的式子
gcd有关
-
(sum_{i=1}^{n}gcd(i,n))
(=sum_{d=1}^{n}(dsum_{i=1}^{n}[gcd(i,n)==d]))
(=sum_{d=1}^{n}dvarphi(frac{n}{d}))(当然到这步就可以(O(nsqrt n))做了)
(=sum_{d=1}^{n}frac{n}{d}varphi(d))
(=sum_{d=1}^{n}frac{n}{d}dprod_{p_i|d}frac{p_i-1}{p_i})
(=nsum_{d=1}^{n}prod_{p_i|d}frac{p_i-1}{p_i})
考虑每个质数选或者不选的生成函数:((b_ifrac{p_i-1}{p_i}+1)),其中(b_i)表示(p_i)的指数,含义为选了这个质数就会有(b_i)种指数
所以最后答案就是(nprod_{p_i|n}(b_ifrac{p_i-1}{p_i}+1))
数论模板库
黑科技
(long long)相乘取模
ll mul(ll x,ll y,ll P) {return (x*y-(ll)((long double)x/P*y+0.5)*P+P)%P;}
子集枚举
(0)~(2^n)的子集个数之和是(3^n)(一个位置可以有三种情况:没被选到、选到但是没有被子集枚举到、选到并被枚举到)
以下可以快速枚举子集
for(int i=s;;i=(i-1)&s)
{
cout<<i<<endl;
if(!i) break;
}
高维前缀和
统计子集
for(int p=1;p<=1<<20;p<<=1)
for(int i=0;i<1<<21;i++)
if(i&p) f[i]+=f[i^p];
统计超集
for(int p=1;p<=1<<20;p<<=1)
for(int i=0;i<1<<21;i++)
if(!(i&p)) f[i]+=f[i|p];
各种线性筛
用接近线性的方法筛取一些函数,首先要先找到其值与质数的关联才能做
积性函数都能线性筛
线性筛质数
【模板】线性筛素数
void Get_pri()
{
ispri[1]=1;
for(int i=2;i<=N;i++)
{
if(!ispri[i]) pri[++tot]=i;
for(int j=1;j<=tot&&i*pri[j]<=N;j++)
{
ispri[i*pri[j]]=1;
if(i%pri[j]==0) break;
}
}
}
线性筛(mu)
void Get_mu()
{
ispri[1]=mu[1]=1;
for(int i=2;i<=N;i++)
{
if(!ispri[i]) pri[++tot]=i,mu[i]=-1;
for(int j=1;j<=tot&&i*pri[j]<=N;j++)
{
ispri[i*pri[j]]=1;
if(i%pri[j]) mu[i*pri[j]]=-mu[i];
else break;
}
}
}
线性筛(varphi)
void Get_phi()
{
ispri[1]=phi[1]=1;
for(int i=2;i<=N;i++)
{
if(!ispri[i]) pri[++tot]=i,phi[i]=i-1;
for(int j=1;j<=tot&&i*pri[j]<=N;j++)
{
ispri[i*pri[j]]=1;
if(i%pri[j]) phi[i*pri[j]]=phi[i]*(pri[j]-1);
else {phi[i*pri[j]]=phi[i]*pri[j];break;}
}
}
}
高级算法
Exgcd
上文有讲用法,注意传入的时候要求都除以了(gcd)
void Exgcd(int a,int b,int &x,int &y)
{
if(!b) {x=1;y=0;return;}
Exgcd(b,a%b,y,x);y-=a/b*x;
}
Lucas
【模板】卢卡斯定理
用于计算模数比较小的时候的组合数,公式看就是(C(n,k)=C(n/p,k/p)*C(n\%p,k\%p))
然后继续递归,边界条件是(k)为(0)时答案是(1)
int C(int n,int k) {return 1ll*jc[n]*inv[k]%P*inv[n-k]%P;}
int Lucas(int n,int k) {return !k?1:1ll*Lucas(n/P,k/P)*C(n%P,k%P)%P;}
int pre()
{
inv[0]=inv[1]=jc[0]=1;
for(int i=2;i<P;i++) inv[i]=(P-1ll*P/i*inv[P%i]%P)%P;
for(int i=1;i<P;i++) jc[i]=1ll*jc[i-1]*i%P,inv[i]=1ll*inv[i]*inv[i-1]%P;
}
EXCRT
【模板】扩展中国剩余定理(EXCRT)
合并同余方程,复杂度瓶颈在于exgcd
ll mul(ll x,ll y,ll p) {return (x*y-(ll)((long double)x/p*y+0.5)*p+p)%p;}
ll gcd(ll a,ll b) {return !b?a:gcd(b,a%b);}
void Exgcd(ll a,ll b,ll &x,ll &y) {!b?(x=1,y=0):(Exgcd(b,a%b,y,x),y-=a/b*x);}
ll n,B1,P1,B2,P2;
int main()
{
cin>>n>>P1>>B1;n--;
while(n--)
{
cin>>P2>>B2;
ll g=gcd(P1,P2),C=(B2-B1)/g,x,y,P;
Exgcd(P1/g,P2/g,x,y); P=P1/g*P2;
B1=(mul(mul(x,C,P),P1,P)+B1+P)%P; P1=P;
}
cout<<B1<<endl;
}
BSGS
[CQOI2018]破解D-H协议
(g^a equiv A(mod P)),给定(A),求(a)。
方法:
- (a=im-j)所以(g^{im}equiv A*g^j),其中(m=sqrt P)能够保证最优复杂度
- 然后从(0)到(m)枚举(i),把(g^{im})存进哈希表
- 暴力判断对于每一个(j),哈希表里是否有(A*g^j)的值
复杂度(O(sqrt PlogP))
const int Mo=100003;
int g,A,a,b,P,n;
int ksm(int x,int k)
struct HashTable
{
struct Line{int next,v1,v2;}a[Mo];
int head[Mo],cnt;
void reset() {memset(head,0,sizeof(head));cnt=0;}
void Add(int p,int v1,int v2) {a[++cnt]=(Line){head[p],v1,v2};head[p]=cnt;}
int Query(int x)
{
for(int i=head[x%Mo];i;i=a[i].next)
if(a[i].v1==x) return a[i].v2;return -1;
}
}Hash;
int BSGS(int g,int A,int P)
{
int M=sqrt(P)+1;Hash.reset();
for(int i=0,t=A;i<M;i++,t=1ll*t*g%P) Hash.Add(t%Mo,t,i);
for(int i=1,bs=ksm(g,M),t=bs;i<=M;i++,t=1ll*t*bs%P)
if(Hash.Query(t)!=-1) return i*M-Hash.Query(t);return -1;
}
高斯消元
【模板】高斯消元法
(O(n^3))求解线性方程组。用线性基的方式实现
const int N=110;
const double eps=1e-7;
int n;
double F[N][N],A[N];
int main()
{
cin>>n;
for(int w=1;w<=n;w++)
{
for(int i=1;i<=n+1;i++) cin>>A[i];
for(int i=1;i<=n;i++)
{
if(fabs(A[i])<eps) continue;
if(fabs(F[i][i])<eps)
{
for(int j=i;j<=n+1;j++) F[i][j]=A[j]/A[i];
F[i][i]=1;break;
}
else for(int j=i+1;j<=n+1;j++) A[j]-=F[i][j]*A[i];
}
}
for(int i=1;i<=n;i++)
if(fabs(F[i][i])<eps)
return puts("No Solution"),0;
for(int i=n;i>=1;i--)
for(int j=i-1;j>=1;j--)
F[j][n+1]-=F[j][i]*F[i][n+1],F[j][i]=0;
for(int i=1;i<=n;i++) printf("%.2f
",F[i][n+1]);
}
线性基
【模板】线性基
很方便地求解一个数集中能够异或出来的数
ll n,x,a[60],ans;
int main()
{
for(cin>>n;n;n--)
{
cin>>x;
for(int i=50;i>=0;i--)
if(x&(1ll<<i))
{
if(!a[i]) {a[i]=x;break;}
else x^=a[i];
}
}
for(int i=50;i>=0;i--)
if((ans^a[i])>ans) ans^=a[i];
return cout<<ans<<endl,0;
}
裴蜀定理
【模板】裴蜀定理
问题:给定(A_1,A_2...A_k),求(A_1x_1+A_2x_2+...A_kx_k)所能取到的最小值
结论:所有数的绝对值的(gcd)
int g,x,n;
int main()
{
cin>>n;while(n--) cin>>x,g=__gcd(g,abs(x));
cout<<g<<endl;
}
FFT
【模板】多项式乘法(FFT)
(O(nlogn))求解两多项式卷积
//FFT
for(int i=1;i<l;i++) r[i]=(r[i>>1]>>1)|((i&1)<<tt);
for(int i=0;i<l;i++) w[i].rl=cos(Pi/l*i),w[i].im=sin(Pi/l*i);
void FFT(Complex *P,int op)
{
for(int i=1;i<l;i++) if(r[i]>i) swap(P[i],P[r[i]]);
for(int i=1;i<l;i<<=1)
for(int j=0,p=i<<1;j<l;j+=p)
for(int k=0;k<i;k++)
{
Complex W=w[l/i*k];W.im*=op;
Complex X=P[j+k],Y=P[j+k+i]*W;
P[j+k]=X+Y;P[j+k+i]=X-Y;
}
}
for(int i=0;i<=m+n;i++) printf("%d ",(int)(A[i].rl/l+0.5));
拉格朗日插值
(n-1)次多项式带入(x=k)求值。
int Lagrange()
{
cin>>n>>k;
for(int i=1;i<=n;i++) cin>>x[i]>>y[i];
for(int i=1;i<=n;i++)
{
int fz=1,fm=1;
for(int j=1;j<=n;j++)
if(i!=j) fz=1ll*fz*(k-x[j])%P,fm=1ll*fm*(x[i]-x[j])%P;
(ans+=1ll*y[i]*fz%P*ksm(fm,P-2)%P)%=P;
}
cout<<(ans+P)%P<<endl;
}
NTT
FWT
这两个一个用NTT代替,一个用高维前缀和。考到了话我直播CS!