扩展GCD
给定
a,b
a
,
b
,求使得
ax+by=gcd(a,b)
a
x
+
b
y
=
g
c
d
(
a
,
b
)
的一组解(x,y)。
我们知道,
gcd(a,b)=gcd(b,a%b)
g
c
d
(
a
,
b
)
=
g
c
d
(
b
,
a
%
b
)
。
所以通过辗转相除法,可以得到这些方程:
a=b∗m+n,m=a div b,n=a%b
a
=
b
∗
m
+
n
,
m
=
a
d
i
v
b
,
n
=
a
%
b
m=n∗m′+n′
m
=
n
∗
m
′
+
n
′
通过这样的轮换,最后得到
gcd(a,b)∗x+0∗y=gcd(a,b)
g
c
d
(
a
,
b
)
∗
x
+
0
∗
y
=
g
c
d
(
a
,
b
)
,返回
x=1,y=0
x
=
1
,
y
=
0
。
ax+by=gcd(a,b)
a
x
+
b
y
=
g
c
d
(
a
,
b
)
bx′+a%b∗y′=gcd(a,b)
b
x
′
+
a
%
b
∗
y
′
=
g
c
d
(
a
,
b
)
合并同类项,最后得出
x=y′,y=x′−y′∗⌊ab⌋
x
=
y
′
,
y
=
x
′
−
y
′
∗
⌊
a
b
⌋
自己推一下式子。
int extgcd(int a,int b,int &x,int &y){
if (!b){
x=1;y=0;
return a;
}
int c=extgcd(b,a%b,x,y),t;
t=x;
x=y;
y=t-y*(a/b);
return c;
}
矩阵快速幂
差点忘了矩乘怎么打…….
struct matrix{
int mat[N][N];
void clear(){
memset(mat,0,sizeof(mat));
}
void unit(){
clear();
for(int i=1;i<=n;i++)mat[i][i]=1;
}
};matrix a,b,c;
matrix mul(matrix a,matrix b){
matrix c;c.clear();
for(int k=1;k<=n;k++)
for(int i=1;i<=n;i++)if(a.mat[i][k])//这个优化很强大
for(int j=1;j<=n;i++)if(b.mat[k][j])
c.mat[i][j]+=a.mat[i][k]*b.mat[k][j];
return c;
}
matrix mksm(matrix a,int p){
matrix b;b.unit();
for(;p;p/=2,a=mul(a,a))
if(p&1) b=mul(b,a);//注意顺序,矩乘没有交换律
return b;
}
线筛求φ(N)和μ(N)
通过线筛质数来求。
void getprime(){
int i,j;
for(i=2;i<=N;i++){
if(!bz[i]){
prime[++prime[0]]=i;
phi[i]=i-1;
mu[i]=-1;
}
fo(j,1,prime[0]){
if(i*prime[j]>N)break;
bz[i*prime[j]]=0;
if(i%prime[j]==0){
mu[i*prime[j]]=0;
phi[i*prime[j]]=phi[i]*prime[j];
break;
} else{
mu[i*prime[j]]=-mu[i];
phi[i*prime[j]]=phi[i]*(prime[j]-1);
}
}
}
}
中国剩余定理
就是给定这么几个方程:
A≡a1(mod p1)
A
≡
a
1
(
m
o
d
p
1
)
A≡a2(mod p2)
A
≡
a
2
(
m
o
d
p
2
)
……
A≡an(mod pn)
A
≡
a
n
(
m
o
d
p
n
)
求最小的A。(一般任意两个p之间互质,如果不互质自己对号入座)
考虑
A=p1∗k1+a1=p2∗k2+a2
A
=
p
1
∗
k
1
+
a
1
=
p
2
∗
k
2
+
a
2
移项,得
p1∗k1−p2∗k2=a2−a1
p
1
∗
k
1
−
p
2
∗
k
2
=
a
2
−
a
1
,即
p1∗k1a2−a1+p2∗−k2a2−a1=1
p
1
∗
k
1
a
2
−
a
1
+
p
2
∗
−
k
2
a
2
−
a
1
=
1
,那么用
extgcd
e
x
t
g
c
d
求出
k1a2−a1
k
1
a
2
−
a
1
和
−k2a2−a1
−
k
2
a
2
−
a
1
就大功告成,合并了两个同余方程。以此类推,将n个同余方程合并到一个的时候,所得到的
a′n
a
n
′
就是最小的A。
void get(){
int i,ans,last,l1,p1,p2,a1,a2;
last=a[1];l1=p[1];
fo(i,2,p[0]){
p1=l1,p2=p[i];
a1=last,a2=a[i];
extgcd(p1,p2,k1,k2);
k1=k1*(a2-a1);
l1=l1*p2;
last=k1*p1+a1;
}
ans=(last%P+P)%P;//P为lcm(p1,p2,...,pn)
return ans;
}
Lucas定理
求
C(n,m)%p
C
(
n
,
m
)
%
p
。
Lucas(n,m,p)=Lucas(⌊np⌋,⌊mp⌋,p)∗C(n%p,m%p,p)
L
u
c
a
s
(
n
,
m
,
p
)
=
L
u
c
a
s
(
⌊
n
p
⌋
,
⌊
m
p
⌋
,
p
)
∗
C
(
n
%
p
,
m
%
p
,
p
)
。
证明:设
n=ap+b,m=cp+d
n
=
a
p
+
b
,
m
=
c
p
+
d
,现证明
Ccp+dap+b=Cca∗Cdb%p)
C
a
p
+
b
c
p
+
d
=
C
a
c
∗
C
b
d
%
p
)
∀x∈Z
∀
x
∈
Z
,
(1+x)ap+b=Σap+bi=0xi∗Cin
(
1
+
x
)
a
p
+
b
=
Σ
i
=
0
a
p
+
b
x
i
∗
C
n
i
(1+x)p≡1+xp (mod p)
(
1
+
x
)
p
≡
1
+
x
p
(
m
o
d
p
)
(1+x)ap+b≡[(1+x)p]a∗(1+x)b≡(1+xp)a∗(1+x)b≡[Σai=0Cia∗xip]∗[Σbi=0Cib∗xi]
(
1
+
x
)
a
p
+
b
≡
[
(
1
+
x
)
p
]
a
∗
(
1
+
x
)
b
≡
(
1
+
x
p
)
a
∗
(
1
+
x
)
b
≡
[
Σ
i
=
0
a
C
a
i
∗
x
i
p
]
∗
[
Σ
i
=
0
b
C
b
i
∗
x
i
]
那么二项式中
xcp+d项的系数≡Cca∗Cdb≡Ccp+dap+b(modp)
x
c
p
+
d
项
的
系
数
≡
C
a
c
∗
C
b
d
≡
C
a
p
+
b
c
p
+
d
(
m
o
d
p
)
∴得证。
(LL)即为long long
LL lucas(LL n,LL m,LL p){
if(!m)return 1;
return lucas(n/p,m/p,p)*C(n%p,m%p);
}
高斯消元
运用加减消元法进行消元,一次消一个,总复杂度 O(n3) O ( n 3 )
void gauss(){
int i,j,k;
double temp;
fo(i,1,n){
j=i;
while(!mat[j][i]&&j<=n)j++;
if(j==n+1)continue;
if(j!=i)fo(k,1,n+1)swap(mat[i][k],mat[j][k]);
fo(j,1,n)
if(j!=i && mat[j][i]){
temp=mat[i][i]/mat[j][i];
fo(k,1,n+1)mat[j][k]=mat[j][k]*temp-mat[i][k];
}
}
fd(i,n,1)
if(mat[i][i]!=0){
fo(j,i+1,n)mat[i][n]-=mat[i][j]*x[j];
x[i]=mat[i][n]/mat[i][i];
}
}
FFT
这里就不详细讲了。
虚数的定义。
struct Cpx{
DB r,i;
Cpx(){}
Cpx(DB R,DB I){r=R,i=I;}
Cpx operator +(const Cpx&x)const{return Cpx(r+x.r,i+x.i);}
Cpx operator -(const Cpx &x)const{return Cpx(r-x.r,i-x.i);}
Cpx operator *(const Cpx &x)const{return Cpx(r*x.r-i*x.i,r*x.i+i*x.r);}
};
void FFT(Cpx *y,int opz){
int i,h,k;
fo(i,0,m-1)if(i<Rev[i])swap(y[i],y[Rev[i]]);
for(h=2;h<=m;h<<=1){
Cpx wn(cos(2*Pi/h),opz*sin(2*Pi/h));
for(i=0;i<m;i+=h){
Cpx w(1,0);
fo(k,i,i+(h>>1)-1){
Cpx u=y[k],t=w*y[k+(h>>1)];
y[k]=u+t;
y[k+(h>>1)]=u-t;
w=w*wn;
}
}
}
if(opz==-1){
fo(i,0,m-1)y[i].r/=m;
}
}
fo(i,0,m-1)Rev[i]=Rev[i>>1]>>1|(i&1)<<len-1;//求i的反过来的二进制。
//意思是,假设目前要加一个位,已经知道不加位的反过来的。如果最后一位是0,反过来后第一位是0.
//最后一位是1同理。
NTT(有原根版)
#include<iostream>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<cstdio>
#define N 200010
#define M 525000
#define LL long long
#define mo 924844033
#define P(a) putchar(a)
#define fo(i,a,b) for(i=a;i<=b;i++)
#define fd(i,a,b) for(i=a;i>=b;i--)
using namespace std;
LL i,j,k,l,n,m,ans;
LL u,v,L,len;
LL jc[N],ny[N],fa[N];
LL siz[N],f[N],h[N];
LL f1[M],f2[M];
LL w[M],_w[M],Rev[M];
LL Ans[N];
LL read(){
LL fh=1,rs=0;char ch;
while((ch<'0'||ch>'9')&&(ch^'-'))ch=getchar();
if(ch=='-')fh=-1,ch=getchar();
while(ch>='0'&&ch<='9')rs=(rs<<3)+(rs<<1)+(ch^'0'),ch=getchar();
return fh*rs;
}
LL ksm(LL x,LL y){
LL rs=1;
for(;y;y>>=1,x=(x*x)%mo)
if(y&1)rs=(rs*x)%mo;
return rs;
}
void pre_(){
int i;
jc[0]=jc[1]=ny[0]=ny[1]=1;
fo(i,2,N-10)jc[i]=(jc[i-1]*i)%mo;
ny[N-10]=ksm(jc[N-10],mo-2);
fd(i,N-11,2)ny[i]=(ny[i+1]*(i+1))%mo;
}
void NTT(LL *y,LL *w){
LL i,h,k;
fo(i,0,L-1)if(i<Rev[i])swap(y[i],y[Rev[i]]);
for(h=2;h<=L;h<<=1){
for(i=0;i<L;i+=h){
fo(k,i,i+(h>>1)-1){
LL u=y[k],t=(w[L/h*(k-i)]*y[k+(h>>1)])%mo;
y[k]=(u+t)%mo;
y[k+(h>>1)]=((u-t)%mo+mo)%mo;
}
}
}
}
int main(){
n=read();
//...
len=0;
for(L=1;L<(n+1)*2;L<<=1)len++;
w[0]=w[L]=1;
w[1]=ksm(5,(mo-1)/L);
fo(i,1,L-1) w[i]=w[i-1]*1ll*w[1]%mo;
fo(i,0,L) _w[i]=w[L-i];
fo(i,0,L-1)Rev[i]=Rev[i>>1]>>1|(i&1)<<len-1;
fo(i,0,n)f1[i]=f[n-i];
fo(i,n+1,L-1)f1[i]=0;
fo(i,0,n)f2[i]=h[i];
fo(i,n+1,L-1)f2[i]=0;
NTT(f1,w);
NTT(f2,w);
fo(i,0,L-1)f1[i]=f1[i]*f2[i]%mo;
NTT(f1,_w);
LL inv=ksm(L,mo-2);
fo(i,0,L-1)f1[i]=(f1[i]*inv)%mo;
}