数学相关模板

行列式和矩阵树定理一直没学通，学了又忘，只得重学一下。

行列式：

$det(A)=det(A^T)$  => 行列式的行和列是等价的

两行互换行列式变号

上(下)三角矩阵行列式为对角线乘积

斜上(下)三角矩阵为斜对角线乘积$ imes (-1)^{frac{n(n-1)}{2}}$

一行乘k，行列式乘k

每行每列和为0的矩阵行列式为0

矩阵树定理：

Kirchhoff矩阵=度数矩阵-邻接矩阵 $C=D-A$

G的生成树个数=C[G]的任意n-1阶矩阵主子式行列式绝对值。

模意义下的高斯消元：

和辗转相除法有点像。

因为$D[i][i]$和$D[j][i]$一定需要消掉一个，就和$gcd(x,y)$最后算到$y=0$的时候才停止一样。

$gcd(x,y)=gcd(y,x \% y)$ 的x、y变成这里的$D[i]$和$D[j]$

$ans(D[i],D[j])=ans(D[j],D[i] \% D[j])=ans(D[j],D[i]-D[j] imes lfloor frac{D[i][i]}{D[j][i]} floor)$

矩阵树定理(bzoj2467)

//Serene
#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdlib>
#include<cstdio>
#include<cmath>
using namespace std;
#define ll long long
const int maxn=400+7,mod=2007;
int T,n,D[maxn][maxn];
 
char cc;ll ff;
template<typename T>void read(T& aa) {
    aa=0;ff=1; cc=getchar();
    while(cc!='-'&&(cc<'0'||cc>'9')) cc=getchar();
    if(cc=='-') ff=-1,cc=getchar();
    while(cc>='0'&&cc<='9') aa=aa*10+cc-'0',cc=getchar();
    aa*=ff;
}
 
int Guess(int n) {
    n--;
    for(int i=1;i<=n;++i) for(int j=1;j<=n;++j) {
        D[i][j]%=mod;
        if(D[i][j]<0) D[i][j]+=mod;
    }
    int a,b,r,fl=0,ans=1;
    for(int i=1;i<=n;++i) {
        for(int j=i+1;j<=n;++j) {
            a=D[i][i];b=D[j][i];
            while(b) {
                r=a/b; a=a%b; swap(a,b);
                for(int k=i;k<=n;++k) D[i][k]=(D[i][k]-r*D[j][k]%mod+mod)%mod;
                for(int k=i;k<=n;++k) swap(D[i][k],D[j][k]);
                fl^=1;
            }
        }
        if(!D[i][i]) return 0;
        (ans*=D[i][i])%=mod;
    }
    if(fl) ans=(mod-ans)%mod;
    return ans;
}
 
int main() {
    read(T); 
    while(T--) {
        read(n);
        memset(D,0,sizeof(D));
        for(int i=1;i<=(n<<2);++i) {
            if(i%4==0) D[i][i]=4;
            else D[i][i]=2;
        }
        for(int i=1;i<n;++i) --D[i<<2][(i+1)<<2],--D[(i+1)<<2][i<<2];
        for(int i=1;i<(n<<2);++i) --D[i][i+1],--D[i+1][i];
        --D[n<<2][1]; --D[1][n<<2];
        --D[n<<2][4]; --D[4][n<<2];
        printf("%d
",Guess(n<<2));
    }
    return 0;
}

Miller-Rabin + Pollard Rho(poj1811)

//Serene
#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdlib>
#include<cstdio>
#include<cmath>
using namespace std;
#define ll long long
#define db double
#define For(i,a,b) for(int i=(a);i<=(b);++i)
#define Rep(i,a,b) for(int i=(a);i>=(b);--i)
const int M=120,W=20;
ll Td,n,ans;

char cc; ll ff;
template<typename T>void read(T& aa) {
	aa=0;cc=getchar();ff=1;
	while((cc<'0'||cc>'9')&&cc!='-') cc=getchar();
	if(cc=='-') ff=-1,cc=getchar();
	while(cc>='0'&&cc<='9') aa=aa*10+cc-'0',cc=getchar();
	aa*=ff;
}

ll gcd(ll x,ll y){return y==0? x:gcd(y,x%y);}

ll qc(ll x,ll k,ll mod) {
	ll rs=0;
	while(k) {
		if(k&1) rs=(rs+x)%mod;
		k>>=1; x=(x+x)%mod;
	}
	return rs;
}

ll qp(ll x,ll k,ll mod) {
	ll rs=1;
	while(k) {
		if(k&1) rs=qc(rs,x,mod);
		k>>=1; x=qc(x,x,mod);
	}
	return rs;
}

bool isprime(ll n) {
	if(n==2||n==3) return 1;
	if(n<2||(n%6!=1&&n%6!=5)) return 0;
	ll m=n-1,k=0,x,y;
	while((m&1)==0) m>>=1,k++;
	For(i,1,W) {
		x=qp(rand()%(n-1)+1,m,n);
		For(j,1,k) {
			y=qc(x,x,n);
			if(y==1&&x!=1&&x!=n-1) return 0;
			x=y;
		}
		if(y!=1) return 0;
	}
	return 1;
}

ll prho(ll n,ll m) {
	ll x=rand()%(n-1)+1,y=0;
	for(ll i=1,k=1,d;y!=x;++i) {
		if(i==k) {y=x;k<<=1;}
		x=(qc(x,x,n)+m)%n;
		d=gcd((y-x+n)%n,n);
		if(d>1&&d<n) return d;
	}
	return n;
}

void find(ll n,ll m) {
	if(isprime(n)) {
		ans=min(ans,n);
		return;
	}
	ll p; for(p=n;p>=n&&m;m--) p=prho(n,m--);
	find(p,M); find(n/p,M);
}

int main() {
	read(Td);
	while(Td--) {
		read(n); ans=n;
		find(n,M);
		if(ans==n) printf("Prime
");
		else printf("%lld
",ans);
	}
	return 0;
}

FFT模板(洛谷）

#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdlib>
#include<cstdio>
#include<cmath>
using namespace std;
#define ll long long
#define db double
const db PI=acos(-1);
const int maxn=2097152+10;
int n,m; 

char cc;ll ff;
template<typename T>void read(T& aa) {
    aa=0;ff=1; cc=getchar();
    while(cc!='-'&&(cc<'0'||cc>'9')) cc=getchar();
    if(cc=='-') ff=-1,cc=getchar();
    while(cc>='0'&&cc<='9') aa=aa*10+cc-'0',cc=getchar();
    aa*=ff;
}

struct Virt{
    db r,i;
    Virt(db r=0.0,db i=0.0) {
        this->r=r;
        this->i=i;
    }
    Virt operator + (const Virt& x) {
        return Virt(r+x.r,i+x.i);
    }
    Virt operator - (const Virt& x) {
        return Virt(r-x.r,i-x.i);
    }
    Virt operator * (const Virt& x) {
        return Virt(r*x.r-i*x.i,r*x.i+i*x.r);
    }
}a[maxn],b[maxn];

void Rader(Virt F[],int len) {
    int i,j,k;
    for(i=1,j=len/2;i<len-1;++i) {
        if(i<j) swap(F[i],F[j]);
        k=len/2;
        while(j>=k) {
            j-=k;
            k>>=1;
        }
        if(j<k) j+=k;
    }
}

void FFT(Virt F[],int len,int on) {
    Rader(F,len);
    for(int h=2;h<=len;h<<=1) {
        Virt wn(cos(-on*2*PI/h),sin(-on*2*PI/h));
        for(int j=0;j<len;j+=h) {
            Virt w(1,0);
            for(int k=j;k<j+h/2;++k) {
                Virt u=F[k];
                Virt t=w*F[k+h/2];
                F[k]=u+t;
                F[k+h/2]=u-t;
                w=w*wn;
            }
        }
    }
    if(on==-1)
        for(int i=0;i<len;++i) F[i].r/=len;
}

int main() {
    read(n); read(m);
    for(int i=0;i<=n;++i) read(a[i].r);
    for(int i=0;i<=m;++i) read(b[i].r);
    m+=n; for(n=1;n<=m+1;n<<=1);
    FFT(a,n,1); FFT(b,n,1);
    for(int i=0;i<n;++i) a[i]=a[i]*b[i];
    FFT(a,n,-1);
    for(int i=0;i<=m;++i) printf("%d ",(int)(a[i].r+0.5));
    return 0;
}

拉格朗日插值法(51nod1258）

//Serene
#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdlib>
#include<cstdio>
#include<cmath>
using namespace std;
#define ll long long
#define db double
#define For(i,a,b) for(int i=(a);i<=(b);++i)
#define Rep(i,a,b) for(int i=(a);i>=(b);--i)
const int maxn=5e4+7,W=5e4+7;
const ll mod=1e9+7;
ll Td,n,k,mi[maxn],y[maxn];

char cc; ll ff;
template<typename T>void read(T& aa) {
	aa=0;cc=getchar();ff=1;
	while((cc<'0'||cc>'9')&&cc!='-') cc=getchar();
	if(cc=='-') ff=-1,cc=getchar();
	while(cc>='0'&&cc<='9') aa=aa*10+cc-'0',cc=getchar();
	aa*=ff;
}

ll qp(ll x,ll k) {
	ll rs=1;
	while(k) {
		if(k&1) rs=rs*x%mod;
		k>>=1; x=x*x%mod;
	}
	return rs;
}

ll inv(ll x) {return qp(x,mod-2);}

void mo(ll &x) {if(x>=mod) x-=mod;}

struct T{
	ll x,y;
	T(){}
	T(ll x,ll y):x(x),y(y){}
	T operator * (const ll& b) const{
		if(b==0) return T(x,y+1);
		return T(x*b%mod,y);
	}
	T operator / (const ll& b) const{
		if(b==0) return T(x,y-1);
		return T(x*inv(b)%mod,y);
	}
}t;

ll solve() {
	ll rs=0; y[0]=0; 
	For(i,1,k+2) y[i]=y[i-1]+qp(i,k),mo(y[i]);
	if(n<=k+2) return y[n];
	t=T(1,0); T s; ll now;
	For(i,0,k+1) t=t*((n%mod-i+mod)%mod);
	For(i,1,k+1) {
		s=t/((n%mod-i+mod)%mod)*y[i];
		if(s.y) continue;
		now=s.x*mi[i]%mod;
		now=now*mi[k+1-i]%mod;
		if((k+1-i)&1) now=now*(mod-1)%mod;
		rs+=now; mo(rs);
	}
	return rs;
}

int main() {
	read(Td);
	mi[0]=1; For(i,1,W+1) mi[i]=mi[i-1]*i%mod;
	mi[W+1]=inv(mi[W+1]); Rep(i,W+1,1) mi[i-1]=mi[i]*i%mod;
	while(Td--) {
		read(n); read(k);
		printf("%lld
",solve());
	}
	return 0;
} 

ud20180620: bzoj3453

//Serene
#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdlib>
#include<cstdio>
#include<cmath>
using namespace std;
#define ll long long
#define db double
#define For(i,a,b) for(ll i=(a);i<=(b);++i)
#define Rep(i,a,b) for(ll i=(a);i>=(b);--i)
const int maxn=2000+7;
const ll mod=1234567891;
ll Td,n,K,A,D,g[maxn],f[maxn];
 
char cc;ll ff;
template<typename T>void read(T& aa) {
    aa=0;ff=1; cc=getchar();
    while(cc!='-'&&(cc<'0'||cc>'9')) cc=getchar();
    if(cc=='-') ff=-1,cc=getchar();
    while(cc>='0'&&cc<='9') aa=aa*10+cc-'0',cc=getchar();
    aa*=ff;
}
 
ll qp(ll x,ll k) {
    ll rs=1;
    while(k) {
        if(k&1) rs=rs*x%mod;
        k>>=1; x=x*x%mod;
    }
    return rs;
}

ll finv(ll x) {return qp(x,mod-2);}
 
ll get_g(ll n) {
    if(n<=K+3) return g[n];
    ll x,y,rs=0;
    For(i,1,K+3) {
        x=g[i]; y=1;
        For(j,1,K+3) if(j!=i) {
            x=x*(n-j+mod)%mod;
            y=y*(i-j+mod)%mod;
        }
        rs+=x*finv(y)%mod;
    }
    return rs%mod;
}
 
ll get_f(ll n) {
    if(n<=K+5) return f[n];
    ll x,y,rs=0;
    For(i,1,K+5) {
        x=f[i]; y=1;
        For(j,1,K+5) if(j!=i) {
            x=x*(n-j+mod)%mod;
            y=y*(i-j+mod)%mod;
        }
        rs+=x*finv(y)%mod;
    }
    return rs%mod;
}
 
int main() {
    read(Td);
    while(Td--) {
        read(K); read(A); read(n); read(D);
        For(i,1,K+3) g[i]=(g[i-1]+qp(i,K))%mod;
        For(i,1,K+3) g[i]=(g[i-1]+g[i])%mod;
        f[0]=get_g(A);
        For(i,1,K+5) f[i]=(f[i-1]+get_g((A+D*i)%mod))%mod;
        printf("%lld
",get_f(n));
    }
    return 0;
}
//bzoj3453

　　

一直没怎么打过杜教筛的题，去bzoj找了个板子打，bzoj3944，卡常，不会卡，弃疗。

//Serene
#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdlib>
#include<cstdio>
#include<cmath>
#include<map>
using namespace std;
#define ll long long
#define For(i,a,b) for(int i=(a);i<=(b);++i)
#define Rep(i,a,b) for(int i=(a);i>=(b);--i)
const int maxn=2e6+7,W=2e6;
int T,n;
 
char cc; ll ff;
template<typename T>void read(T& aa) {
    aa=0;cc=getchar();ff=1;
    while((cc<'0'||cc>'9')&&cc!='-') cc=getchar();
    if(cc=='-') ff=-1,cc=getchar();
    while(cc>='0'&&cc<='9') aa=aa*10+cc-'0',cc=getchar();
    aa*=ff;
}
 
ll phi[maxn],mu[maxn]; int prime[maxn],totp;
bool ok[maxn];
void getp() {
    phi[1]=mu[1]=1;
    For(i,2,W) {
        if(!ok[i]) {
            prime[++totp]=i;
            phi[i]=i-1; mu[i]=-1;
        }
        for(int j=1;j<=totp&&prime[j]<=W/i;++j) {
            ok[i*prime[j]]=1;
            if(i%prime[j]) {
                phi[i*prime[j]]=phi[i]*(prime[j]-1);
                mu[i*prime[j]]=-mu[i];
            }
            else {
                phi[i*prime[j]]=phi[i]*prime[j];
                mu[i*prime[j]]=0;
                break;
            }
        }
    }
    For(i,2,W) phi[i]+=phi[i-1],mu[i]+=mu[i-1];
}
 
map<ll,ll> P,M;
 
inline ll get_phi(ll n) {
    if(n<=W) return phi[n];
    if(P[n]) return P[n];
    ll rs=n*(n+1)/2;
    for(ll i=2,j;i<=n;i=j+1) {
        j=n/(n/i);
        rs-=get_phi(n/i)*(j-i+1);
    }
    return P[n]=rs;
}
 
inline ll get_mu(ll n) {
    if(n<=W) return mu[n];
    if(M[n]) return M[n];
    ll rs=1;
    for(ll i=2,j;i<=n;i=j+1) {
        j=n/(n/i);
        rs-=get_mu(n/i)*(j-i+1);
    }
    return M[n]=rs;
}
 
int main() {
    getp();
    read(T);
    while(T--) {
        read(n);
        printf("%lld %lld
",get_phi(n),get_mu(n));
    }
    return 0;
}
/*
334
345287345
*/

　　

扩展欧拉定理(bzoj4869)

>=phi的情况我直接在快速幂里面判断了

//Serene
#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdlib>
#include<cstdio>
#include<cmath>
using namespace std;
#define ll long long
#define db double
#define For(i,a,b) for(int i=(a);i<=(b);++i)
#define Rep(i,a,b) for(int i=(a);i>=(b);--i)
const int maxn=2e5+7;
ll n,m,num,phi[maxn],a[maxn],W=2e5;

char cc; ll ff;
template<typename T>void read(T& aa) {
	aa=0;cc=getchar();ff=1;
	while((cc<'0'||cc>'9')&&cc!='-') cc=getchar();
	if(cc=='-') ff=-1,cc=getchar();
	while(cc>='0'&&cc<='9') aa=aa*10+cc-'0',cc=getchar();
	aa*=ff;
}

ll get_phi(ll x) {
	ll rs=1,r=x;
	for(int i=2;i<=r/i;++i) if(x%i==0){
		rs*=(i-1); x/=i;
		while(x%i==0) x/=i,rs*=i;
	}
	if(x>1) rs*=(x-1);
	return rs;
}

ll sum[4*maxn],minnum[4*maxn],ql,qr;

ll qp(ll x,ll k,ll mod) {
	ll rs=1,t=0;
	while(k) {
		if(k&1) {
			rs=rs*x;
			if(rs>=mod) t=mod,rs%=mod;
		}
		x=x*x; if(x>=mod) t=mod,x%=mod;
		k>>=1;
	}
	return rs+t;
}

ll get_num(ll x,ll t,ll tot) {
	if(t==tot) return x>=phi[t]? x%phi[t]+phi[t]:x;
	ll r=get_num(x,t+1,tot);
	return qp(num,r,phi[t]);
}

void bld(int pos,int l,int r) {
	if(l==r) {
		sum[pos]=a[l];
		return;
	}
	int mid=(l+r)>>1;
	bld(pos<<1,l,mid);
	bld(pos<<1|1,mid+1,r);
	sum[pos]=(sum[pos<<1]+sum[pos<<1|1])%phi[0];
}

void chge(int pos,int l,int r) {
	if(minnum[pos]>W) return;
	if(l==r) {
		++minnum[pos];
		sum[pos]=qp(num,get_num(a[l],1,minnum[pos]),phi[0])%phi[0];
		return;
	} 
	int mid=(l+r)>>1;
	if(ql<=mid) chge(pos<<1,l,mid); 
	if(qr>mid) chge(pos<<1|1,mid+1,r);
	minnum[pos]=min(minnum[pos<<1],minnum[pos<<1|1]);
	sum[pos]=(sum[pos<<1]+sum[pos<<1|1])%phi[0];
}

ll q(int pos,int l,int r) {
	if(l>=ql&&r<=qr) return sum[pos];
	int mid=(l+r)>>1;ll rs=0;
	if(ql<=mid) rs+=q(pos<<1,l,mid);
	if(qr>mid) rs+=q(pos<<1|1,mid+1,r);
	return rs%phi[0]; 
}

int main() {
	ll op;
	read(n); read(m); read(phi[0]); read(num);
	For(i,1,W) {
		phi[i]=get_phi(phi[i-1]);
		if(phi[i]==1&&phi[i-1]==1) { W=i; break; }
	}
	For(i,W,n) phi[i]=1;
	For(i,1,n) read(a[i]);
	bld(1,1,n);
	For(i,1,m) {
		read(op); read(ql); read(qr);
		if(op==0) chge(1,1,n);
		else printf("%lld
",q(1,1,n));
	}
	return 0;
}

　　

勾股数定理 (bzoj1041)

根据勾股数定理，如果有PPT(a,b,c)(a,b,c互质)，必对应一组s,t使得

a=2*s*t ，b=s*s-t*t ，c=s*s+t*t (s与t互质并且s、t 一奇一偶)

对于这道题，我们先枚举导出组(ak,bk,ck)的k，然后就知道了c，就枚举s就可以了

//Serene
#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdlib>
#include<cstdio>
#include<cmath>
using namespace std;
#define ll long long
#define db double
#define For(i,a,b) for(ll i=(a);i<=(b);++i)
#define Rep(i,a,b) for(ll i=(a);i>=(b);--i)
ll n;

char cc;ll ff;
template<typename T>void read(T& aa) {
	aa=0;ff=1; cc=getchar();
	while(cc!='-'&&(cc<'0'||cc>'9')) cc=getchar();
	if(cc=='-') ff=-1,cc=getchar();
	while(cc>='0'&&cc<='9') aa=aa*10+cc-'0',cc=getchar();
	aa*=ff;
}

ll gcd(ll x,ll y) {
	return y==0? x:gcd(y,x%y);
}

ll cal(ll x) {//x=s*s+t*t,a=2*s*t,b=s*s-t*t
	ll r=sqrt(x),rs=0,t;
	For(s,1,r) {
		t=sqrt(x-s*s);
		if(s*s+t*t!=x||t==0||gcd(s,t)!=1||(s&t&1)) continue;
		rs++;
	}
	return rs;
}

ll get_ans() {
	ll r=sqrt(n),rs=0;
	For(i,1,r) if(n%i==0){//枚举导出组的k
		rs+=cal(n/i);
		if(i!=n/i) rs+=cal(i);
	}
	return rs;
}

int main() {
	read(n);
	ll ans=get_ans();
	printf("%lld",ans*4+4);//四个象限，加上坐标轴上的4个点.
	return 0;
}

中国剩余定理(51Nod1079)

注意exgcd之后x可能是负数，模b加b就可以了

//Serene
#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdlib>
#include<cstdio>
#include<cmath>
using namespace std;
#define ll long long
#define db double
#define For(i,a,b) for(int i=(a);i<=(b);++i)
#define Rep(i,a,b) for(int i=(a);i>=(b);--i)
const ll maxn=17;
ll n,P[maxn],M[maxn],N;

char cc; ll ff;
template<typename T>void read(T& aa) {
	aa=0;cc=getchar();ff=1;
	while((cc<'0'||cc>'9')&&cc!='-') cc=getchar();
	if(cc=='-') ff=-1,cc=getchar();
	while(cc>='0'&&cc<='9') aa=aa*10+cc-'0',cc=getchar();
	aa*=ff;
}

void exgcd(ll a,ll b,ll& x,ll& y) {
	if(b==0) {x=1;y=0;return;}
	exgcd(b,a%b,y,x); y-=a/b*x;
}

int main() {
	read(n); N=1;
	For(i,1,n) {
		read(P[i]); read(M[i]);
		N=N*P[i];
	}
	ll num=0,a,b,x,y;
	For(i,1,n) {
		a=N/P[i]; b=P[i];
		exgcd(a,b,x,y);
		x=(x%b+b)%N*a%N*M[i]%N;
		num=(num+x)%N;
	}
	printf("%lld
",num);
	return 0;
}

原根(51Nod1135)

求$n$的原根$g$:令$r=phi (n)$ （当$n$为素数的时候就是$n-1$啦）

把$r$分解质因数，分别是$p_1$到$p_k$。

那么$g ^ {frac{r}{p_i}} != 1 (mod n)$

//Serene
#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdlib>
#include<cstdio>
#include<cmath>
using namespace std;
#define ll long long
#define db double
#define For(i,a,b) for(int i=(a);i<=(b);++i)
#define Rep(i,a,b) for(int i=(a);i>=(b);--i)
const int maxn=1e6+7;
ll n,x,mod,r,prime[maxn],totp;

char cc; ll ff;
template<typename T>void read(T& aa) {
	aa=0;cc=getchar();ff=1;
	while((cc<'0'||cc>'9')&&cc!='-') cc=getchar();
	if(cc=='-') ff=-1,cc=getchar();
	while(cc>='0'&&cc<='9') aa=aa*10+cc-'0',cc=getchar();
	aa*=ff;
}

ll qp(ll x,ll k) {
	ll rs=1;
	while(k) {
		if(k&1) rs=rs*x%mod;
		k>>=1; x=x*x%mod;
	}
	return rs;
}

int main() {
	read(n); mod=n; n--;
	r=sqrt(n); x=n;
	For(i,2,r) if(x%i==0){
		prime[++totp]=i;
		while(x%i==0) x/=i;
	}
	int fl;
	For(i,2,n) {
		fl=0;
		For(j,1,totp) if(qp(i,n/prime[j])==1) {
			fl=1; break;
		}
		if(!fl) {
			printf("%d
",i);
			break;
		}
	}
	return 0;
}

多项式求逆、多项式除法(bzoj3456)

有两种方法证明：http://blog.miskcoo.com/2015/05/polynomial-inverse、https://www.cnblogs.com/joyouth/p/5587648.html

注意每次FFT前后多出去的部分要清零，否则会wa，还要注意79行是N<<1而不是N

//Serene
#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdlib>
#include<cstdio>
#include<cmath>
using namespace std;
#define ll long long
#define db double
#define For(i,a,b) for(int i=(a);i<=(b);++i)
#define Rep(i,a,b) for(int i=(a);i>=(b);--i)
const int maxn=5e6+7;
const ll mod=1004535809,B=3;
ll n,mi[maxn],inv[maxn],C[maxn],G[maxn],Ginv[maxn],X[maxn];

char cc; ll ff;
template<typename T>void read(T& aa) {
	aa=0;cc=getchar();ff=1;
	while((cc<'0'||cc>'9')&&cc!='-') cc=getchar();
	if(cc=='-') ff=-1,cc=getchar();
	while(cc>='0'&&cc<='9') aa=aa*10+cc-'0',cc=getchar();
	aa*=ff;
}

void debug() {
	printf("Ginv:
");
	For(i,0,n) printf("%lld ",Ginv[i]);
	printf("
"); 
}

ll qp(ll x,ll k) {
	ll rs=1;
	while(k) {
		if(k&1) rs=rs*x%mod;
		k>>=1; x=x*x%mod;
	}
	return rs;
}

ll qp1(ll x,ll k) {
	if(k<0) return qp(qp(x,mod-2),-k);
	else return qp(x,k);
}

void Rader(ll F[],ll len) {
	for(int i=1,j=len/2,k;i<len-1;++i) {
		if(i<j) swap(F[i],F[j]);
		k=len>>1;
		while(j>=k) {j-=k,k>>=1;}
		if(j<k) j+=k;
	}
}

void FFT(ll F[],ll len,ll on) {
	Rader(F,len);
	for(int h=2;h<=len;h<<=1) {
		ll wn=qp1(B,(mod-1)*on/h);
		for(int j=0;j<len;j+=h) {
			ll w=1;
			for(int i=j;i<j+h/2;++i) {
				ll u=F[i],v=w*F[i+h/2]%mod;
				F[i]=(u+v)%mod;
				F[i+h/2]=(u-v+mod)%mod;
				w=w*wn%mod;
			}
		}
	}
	ll x=qp1(len,mod-2);
	if(on==-1) For(i,0,len) F[i]=F[i]*x%mod;
}

void get_Ginv(ll N) {
	if(N==1) {
		Ginv[0]=qp(G[0],mod-2);
		return;
	}
	get_Ginv((N+1)>>1);
	ll n=1; for(;n<=(N<<1);n<<=1);
	For(i,0,N-1) X[i]=G[i]; For(i,N,n) X[i]=0;
	FFT(X,n,1); FFT(Ginv,n,1);
	For(i,0,n-1) Ginv[i]=Ginv[i]*(2-Ginv[i]*X[i]%mod+mod)%mod;
	FFT(Ginv,n,-1);
	For(i,N,n) Ginv[i]=0;
//	printf("N=%lld:
",N); debug();
}

int main() {
	read(n);
	mi[0]=1; For(i,1,n) mi[i]=mi[i-1]*i%mod;
	inv[n]=qp(mi[n],mod-2); Rep(i,n,1) inv[i-1]=inv[i]*i%mod;
	For(i,1,n) C[i]=qp(2,(ll)i*(i-1)/2)*inv[i-1]%mod;
	For(i,0,n) G[i]=qp(2,(ll)i*(i-1)/2)*inv[i]%mod;
	get_Ginv(n+1);
//	debug();
	ll ans=0,len=1; for(;len<=n;len<<=1);
	For(i,n+1,len) C[i]=Ginv[i]=0;
	FFT(Ginv,len,1); FFT(C,len,1);
	For(i,0,len) Ginv[i]=Ginv[i]*C[i]%mod;
	FFT(Ginv,len,-1);
	ans=Ginv[n]*mi[n-1]%mod;
	printf("%lld
",ans);
//	cerr<<clock()<<"
";
	return 0;
}

这道题也可以用多项式求ln做：

//Serene
#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdlib>
#include<cstdio>
#include<cmath>
using namespace std;
#define ll long long
#define db double
#define For(i,a,b) for(int i=(a);i<=(b);++i)
#define Rep(i,a,b) for(int i=(a);i>=(b);--i)
const int maxn=1e6+7;
ll mod=1004535809,B=3;
ll n,mi[maxn],inv[maxn],G[maxn],Ginv[maxn],Gln[maxn],X[maxn];

char cc; ll ff;
template<typename T>void read(T& aa) {
	aa=0;cc=getchar();ff=1;
	while((cc<'0'||cc>'9')&&cc!='-') cc=getchar();
	if(cc=='-') ff=-1,cc=getchar();
	while(cc>='0'&&cc<='9') aa=aa*10+cc-'0',cc=getchar();
	aa*=ff;
}

ll qp(ll x,ll k) {
	ll rs=1;
	while(k) {
		if(k&1) rs=rs*x%mod;
		k>>=1; x=x*x%mod;
	}
	return rs;
}

ll finv(ll x) {return qp(x,mod-2);}

ll qp1(ll x,ll k) {
	if(k<0) return qp(finv(x),-k);
	return qp(x,k);
}

void Rader(ll F[],ll len) {
	for(int i=1,j=len>>1,k;i<len-1;++i) {
		if(i<j) swap(F[i],F[j]);
		k=len>>1;
		while(j>=k) {j-=k;k>>=1;}
		if(j<k) j+=k;
	}
}

void FFT(ll F[],ll len,ll on) {
	Rader(F,len);
	for(int h=2;h<=len;h<<=1) {
		ll wn=qp1(B,(mod-1)*on/h);
		for(int j=0;j<len;j+=h) {
			ll w=1;
			for(int i=j;i<j+h/2;++i) {
				ll u=F[i],v=w*F[i+h/2]%mod;
				F[i]=(u+v)%mod;
				F[i+h/2]=(u-v+mod)%mod;
				w=w*wn%mod;
			}
		}
	}
	ll x=finv(len);
	if(on==-1) For(i,0,len) F[i]=F[i]*x%mod;
}

void get_inv(ll N,ll* A,ll *B) {
	if(N==1) {
		B[0]=finv(A[0]);
		return;
	}
	get_inv((N+1)>>1,A,B);
	ll n=1;for(;n<=(N<<1);n<<=1);
	For(i,0,N-1) X[i]=A[i]; For(i,N,n) X[i]=B[i]=0;
	FFT(X,n,1); FFT(B,n,1);
	For(i,0,n) B[i]=B[i]*(2-B[i]*X[i]%mod+mod)%mod;
	FFT(B,n,-1); For(i,N,n) B[i]=0;
}

void get_ln(ll n,ll* A,ll *B) {
	get_inv(n+1,A,Ginv);
	For(i,0,n) B[i]=A[i];
	For(i,1,n) B[i-1]=B[i]*i%mod; B[n]=0;
	ll len=1;for(;len<=(n<<1);len<<=1);
	FFT(B,len,1); FFT(Ginv,len,1);
	For(i,0,len) B[i]=B[i]*Ginv[i]%mod;
	FFT(B,len,-1);
	Rep(i,n,1) B[i+1]=B[i]*finv(i+1)%mod; B[0]=0;
}

int main() {
	read(n);
	mi[0]=1; For(i,1,n) mi[i]=mi[i-1]*i%mod;
	inv[n]=qp(mi[n],mod-2);
	Rep(i,n,1) inv[i-1]=inv[i]*i%mod;
	For(i,0,n) G[i]=qp(2,(ll)i*(i-1)/2)*inv[i]%mod;
	get_ln(n,G,Gln);
	printf("%lld
",Gln[n]*mi[n]%mod);
	return 0;
}
/*
Method 2: get ln(G(x))
*/

多项式开方(cf round 200 E - The Child and Binary Tree)

题解：https://blog.csdn.net/q582116859/article/details/77804903

bz上很卡常，没卡过，只能去cf上交了。

注意每次FFT之前清零。

//Serene
#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdlib>
#include<cstdio>
#include<cmath>
using namespace std;
#define ll long long
#define db double
#define For(i,a,b) for(int i=(a);i<=(b);++i)
#define Rep(i,a,b) for(int i=(a);i>=(b);--i)
const int maxn=1e6+7;
ll mod=998244353,B=3,I2=499122177;
ll n,m,G[maxn],Gsq[maxn],Ginv[maxn],X[maxn],C[maxn];

char cc; ll ff;
template<typename T>void read(T& aa) {
	aa=0;cc=getchar();ff=1;
	while((cc<'0'||cc>'9')&&cc!='-') cc=getchar();
	if(cc=='-') ff=-1,cc=getchar();
	while(cc>='0'&&cc<='9') aa=aa*10+cc-'0',cc=getchar();
	aa*=ff;
}

void debug(ll *F) {
    For(i,0,n) printf("%lld ",F[i]);
    printf("
");
}

ll qp(ll x,ll k) {
	ll rs=1;
	while(k) {
		if(k&1) rs=rs*x%mod;
		k>>=1; x=x*x%mod;
	}
	return rs;
}

ll inv(ll x) {return qp(x,mod-2);}

ll qp1(ll x,ll k) {
	if(k<0) return qp(inv(x),-k);
	return qp(x,k);
}

bool cansq(ll x) {
	return qp((x+mod)%mod,(mod-1)>>1)==1;
}
/*
ll sq(ll x) {
	if(!cansq(x)) cerr<<"!!!
";
	ll a=rand()%(2*mod),wpf;
	while(!a||cansq(a*a-x)) a=rand()%(2*mod);
	wpf=(a*a-x+mod)%mod;
	ll rs1=1,rs2=0,x1=a,x2=1,k=(mod+1)>>1,y1,y2;
	while(k) {
		if(k&1) {
			y1=rs1*x1%mod+rs2*x2%mod*wpf%mod;
			y2=rs1*x2%mod+rs2*x1%mod;
			rs1=y1%mod; rs2=y2%mod;
		}
		k>>=1;
		y1=x1*x1%mod+x2*x2%mod*wpf%mod;
		y2=2*x1*x2%mod;
		x1=y1%mod; x2=y2%mod;
	}
	return rs1;
}
*/
void Rader(ll F[],ll len) {
	for(int i=1,j=len>>1,k;i<len-1;++i) {
		if(i<j) swap(F[i],F[j]);
		k=len>>1;
		while(j>=k) {j-=k;k>>=1;}
		if(j<k) j+=k;
	}
}

void FFT(ll F[],ll len,ll on) {
	Rader(F,len);
	for(int h=2;h<=len;h<<=1) {
		ll wn=qp1(B,(mod-1)*on/h);
		for(int j=0;j<len;j+=h) {
			ll w=1;
			for(int i=j;i<j+h/2;++i) {
				ll u=F[i],v=w*F[i+h/2]%mod;
				F[i]=(u+v)%mod;
				F[i+h/2]=(u-v+mod)%mod;
				w=w*wn%mod;
			}
		}
	}
	ll x=qp(len,mod-2);
	if(on==-1) For(i,0,len) F[i]=F[i]*x%mod;
}

void get_inv(ll N,ll* A,ll* B) {
	if(N==1) {
		B[0]=inv(A[0]);
		return;
	}
	get_inv((N+1)>>1,A,B);
	ll n=1;for(;n<=(N<<1);n<<=1);
	For(i,0,N-1) X[i]=A[i]; For(i,N,n) X[i]=B[i]=0;
	FFT(X,n,1); FFT(B,n,1);
	For(i,0,n) B[i]=B[i]*(2-B[i]*X[i]%mod+mod)%mod;
	FFT(B,n,-1); For(i,N,n) B[i]=0;
//	printf("get_inv %lld:
",N); debug(B);
}

void get_sq(ll N,ll* A,ll* B) {
	if(N==1) {
		B[0]=1;//B[0]=sq(A[0]);
		return;
	}
	get_sq((N+1)>>1,A,B);//B=(B^2+A)/(2B)
//	printf("move into get_inv
");
	get_inv(N,B,Ginv);
//	printf("return get_sq
");
	ll n=1;for(;n<=(N<<2);n<<=1);
	For(i,0,N-1) X[i]=A[i]; For(i,N,n) X[i]=B[i]=Ginv[i]=0;
	FFT(B,n,1); FFT(Ginv,n,1); FFT(X,n,1);
	For(i,0,n) B[i]=(B[i]*B[i]%mod+X[i]+mod)*Ginv[i]%mod*I2%mod;
	FFT(B,n,-1); For(i,N,n) B[i]=0;
//	printf("get_sq %lld:
",N); debug(B);
}

int main() {//F=2/(1+sqrt(1-4G))
	read(m); read(n); int x;
	For(i,1,m) { read(x); if(x<=n) G[x]=1; }
	For(i,0,n) if(G[i]) G[i]=mod-4;G[0]=1;
	get_sq(n+1,G,Gsq);
	memcpy(G,Gsq,sizeof(Gsq));
	G[0]=(G[0]+1)%mod;
	memset(Ginv,0,sizeof(Ginv));
	get_inv(n+1,G,Ginv);
	For(i,0,n) Ginv[i]=Ginv[i]*2%mod;
	For(i,1,n) printf("%lld
",Ginv[i]);
	return 0;
}

多项式快速幂、拉格朗日反演(bzoj3684)

拉格朗日反演：

llj大佬都不会证我怎么会证

题解链接：https://blog.csdn.net/PoPoQQQ/article/details/46366549

一些细节：

FFT前后清零！FFT前后清零！FFT前后清零！

求$frac{x}{G(x)}$这个东西，你发现$G(x)$的0次项系数为0，所以直接把系数平移就好了，就不用算了inv又去*x了。

get_pow套get_exp套get_ln套get_inv的时候，注意不要为了节dai省ma空hao间xie用相同的数组互相影响，宁愿多开几个数组，多copy几次。

//Serene
#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdlib>
#include<cstdio>
#include<cmath>
using namespace std;
#define ll long long
#define db double
#define For(i,a,b) for(int i=(a);i<=(b);++i)
#define Rep(i,a,b) for(int i=(a);i>=(b);--i)
const int maxn=1e6+7;
const ll mod=950009857,B=7;
ll n,m,G[maxn],Ginv[maxn],Gln[maxn],Gexp[maxn],X[maxn];

char cc; ll ff;
template<typename T>void read(T& aa) {
	aa=0;cc=getchar();ff=1;
	while((cc<'0'||cc>'9')&&cc!='-') cc=getchar();
	if(cc=='-') ff=-1,cc=getchar();
	while(cc>='0'&&cc<='9') aa=aa*10+cc-'0',cc=getchar();
	aa*=ff;
}

ll qp(ll x,ll k) {
	ll rs=1;
	while(k) {
		if(k&1) rs=rs*x%mod;
		k>>=1; x=x*x%mod;
	}
	return rs;
}

ll finv(ll x) {return qp(x,mod-2);}

ll qp1(ll x,ll k) {
	if(k<0) return qp(finv(x),-k);
	return qp(x,k);
}

void Rader(ll F[],ll len) {
	for(int i=1,j=len>>1,k;i<len-1;++i) {
		if(i<j) swap(F[i],F[j]);
		k=len>>1;
		while(j>=k) {j-=k;k>>=1;}
		if(j<k) j+=k;
	}
}

void FFT(ll F[],ll len,ll on) {
	Rader(F,len);
	for(int h=2;h<=len;h<<=1) {
		ll wn=qp1(B,(mod-1)*on/h);
		for(int j=0;j<len;j+=h) {
			ll w=1;
			for(int i=j;i<j+h/2;++i) {
				ll u=F[i],v=w*F[i+h/2]%mod;
				F[i]=(u+v)%mod;
				F[i+h/2]=(u-v+mod)%mod;
				w=w*wn%mod;
			}
		}
	}
	ll x=finv(len);
	if(on==-1) For(i,0,len) F[i]=F[i]*x%mod;
}

void get_inv(ll N,ll* A,ll* B) {
	if(N==1) {
		B[0]=finv(A[0]);
		return;
	}
	get_inv((N+1)>>1,A,B);
	ll n=1;for(;n<=(N<<1);n<<=1);
	For(i,0,N-1) X[i]=A[i]; For(i,N,n) X[i]=B[i]=0;
	FFT(B,n,1); FFT(X,n,1);
	For(i,0,n) B[i]=B[i]*(2-B[i]*X[i]%mod+mod)%mod;
	FFT(B,n,-1); For(i,N,n) B[i]=0;
}

void get_ln(ll N,ll* A,ll *B) {
	For(i,0,N-1) B[i]=A[i];
	get_inv(N,A,Ginv);
	For(i,1,N) B[i-1]=B[i]*i%mod; B[N]=0;
	ll n=1;for(;n<=(N<<1);n<<=1);
	For(i,N,n) B[i]=Ginv[i]=0;
	FFT(B,n,1); FFT(Ginv,n,1);
	For(i,0,n) B[i]=B[i]*Ginv[i]%mod;
	FFT(B,n,-1); For(i,N,n) B[i]=0;
	Rep(i,N,1) B[i]=B[i-1]*finv(i)%mod; B[0]=0;
//	printf("get_ln %lld
",N);
//	For(i,0,N-1) printf("%lld ",B[i]);
//	printf("
");
}

void get_exp(ll N,ll *A,ll* B) {
	if(N==1) {
		B[0]=1;
		return;
	}
	get_exp((N+1)>>1,A,B);
	get_ln(N,B,Gln);
	ll n=1;for(;n<=(N<<1);n<<=1);
	For(i,0,N-1) X[i]=(A[i]-Gln[i]+mod)%mod; X[0]=(X[0]+1)%mod;
	For(i,N,n) X[i]=B[i]=0;
	FFT(B,n,1); FFT(X,n,1);
	For(i,0,n) B[i]=B[i]*X[i]%mod;
	FFT(B,n,-1); For(i,N,n) B[i]=0;
//	printf("get_exp %lld
",N);
//	For(i,0,N-1) printf("%lld ",B[i]);
//	printf("
");
}

ll get_qpow(ll N,ll* A,ll* B,ll k) {
	get_ln(N,A,Gln);
	For(i,0,N-1) A[i]=Gln[i]*k%mod;
	get_exp(N,A,B);
}

int main() {
	read(n); read(m);
	ll x;
	For(i,1,m) {
		read(x);
		G[x-1]=mod-1;
	}
	G[0]=(G[0]+1)%mod;
	get_inv(n+1,G,Ginv);
	For(i,0,n) G[i]=Ginv[i];
	get_qpow(n+1,G,Gexp,n);
	printf("%lld
",Gexp[n-1]*finv(n)%mod);
	return 0;
}

Min-Max容斥

loj2542

//Serene
#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdlib>
#include<cstdio>
#include<cmath>
using namespace std;
#define ll long long
#define db double
#define For(i,a,b) for(int i=(a);i<=(b);++i)
#define Rep(i,a,b) for(int i=(a);i>=(b);--i)
const int maxn=23,maxt=(1<<18)+7;
const ll mod=998244353;
ll n,m,root,inv[maxn],U,cnt[maxt];
ll f[maxt],g[maxt],A[maxn],B[maxn];

char cc;ll ff;
template<typename T>void read(T& aa) {
	aa=0;cc=getchar();ff=1;
	while((cc<'0'||cc>'9')&&cc!='-') cc=getchar();
	if(cc=='-') ff=-1,cc=getchar();
	while(cc>='0'&&cc<='9') aa=aa*10+cc-'0',cc=getchar();
	aa*=ff;
}

ll qp(ll x,ll k) {
	ll rs=1;
	while(k) {
		if(k&1) rs=rs*x%mod;
		k>>=1; x=x*x%mod;
	}
	return rs;
}

int fir[maxn],nxt[2*maxn],to[2*maxn],e=0;
void add(int x,int y) {
	to[++e]=y;nxt[e]=fir[x];fir[x]=e;
	to[++e]=x;nxt[e]=fir[y];fir[y]=e;
}

int fa[maxn];ll d[maxn];//d:sum of son
void dfs(int pos) {
	int y,z;
	for(y=fir[pos];y;y=nxt[y]) {
		if((z=to[y])==fa[pos]) continue;
		fa[z]=pos; dfs(z); ++d[pos];
	}
}

void get_f(int pos,int t) {
	if((t>>(pos-1))&1) {
		A[pos]=B[pos]=f[t]=0;
		return;
	}
	if(!d[pos]) {
		A[pos]=1; B[pos]=1;
		return;
	}
	int y,z;ll a=0,b=0;
	for(y=fir[pos];y;y=nxt[y]) {
		if((z=to[y])==fa[pos]) continue;
		get_f(z,t);
		a+=A[z]; b+=B[z];
	}
	if(pos!=root) {
		a=(d[pos]+1-a%mod+mod)%mod; 
		a=qp(a,mod-2);
		b=(b+d[pos]+1)%mod;
		A[pos]=a; B[pos]=b*a%mod;
	}
	else {
		a=(d[pos]-a%mod+mod)%mod;
		a=qp(a,mod-2);
		b=(b+d[pos])%mod;
		f[t]=b*a%mod;
	}
}

int main() {
	read(n); read(m); read(root);
	inv[0]=1; For(i,1,n+1) inv[i]=qp(i,mod-2);
	int x,y,t;
	For(i,1,n-1) {
		read(x); read(y);
		add(x,y);
	}
	U=1<<n; dfs(root);
	For(i,1,U-1) {
		get_f(root,i);
		cnt[i]=cnt[i&(i-1)]^1;
		if(!cnt[i]) f[i]=(mod-f[i])%mod;
	}
	For(i,1,U-1) {
		for(int j=i;j;j=(j-1)&i) g[i]+=f[j];
		g[i]%=mod;
	}
	For(i,1,m) {
		read(x); t=0;
		For(j,1,x) {
			read(y);
			t|=(1<<y-1);
		}
		printf("%lld
",g[t]);
	}
	return 0;
}

bzoj4555

计算

$S(i, j)$表示第二类斯特林数

$n leq 1e5$

第二类斯特林数？我不会呀。

第二类Stirling数实际上是集合的一个拆分，表示将n个不同的元素拆分成m个集合的方案数

推式子？麻烦，我们来看定义。$S(n,m) imes (m!)$就是把n个不同的元素拆分成$m$个不同的有序集合的方案数。

令$F(n) = sumlimits_{m=0}^{n} S(n,m) imes (m!) $，这个的含义就是将$i$个不同的数分到任意多个有序集合的方案数。

那么$F$的递推公式就是$F(n)= sumlimits_{i=1}^{n} F(n-i) C(n,i)$。

如果我们令$F(n) = sumlimits_{m=0}^{n} S(n,m) imes (m!) imes 2^m$

那么$F$的递推公式就是$F(n) = sumlimits_{i=1}^{n} 2 F(n-i) C(n,i) $

令$f(n)=frac{F(n)}{n!}$，那么式子变成$f(n) = sumlimits_{i=1}^{n} f(n-i) imes frac{2}{i!}$

这是一个卷积的形式，如果令$g(n)=frac{2}{n!}$，那么$f=f*g+1$，因为$f(0)=1$

所以$f=frac{1}{1-g}$，就是多项式求逆的裸题啦。

//Serene
#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdlib>
#include<cstdio>
#include<cmath>
using namespace std;
#define ll long long
#define db double
#define For(i,a,b) for(int i=(a);i<=(b);++i)
#define Rep(i,a,b) for(int i=(a);i>=(b);--i)
const int maxn=1e6+7;
const ll mod=998244353,B=3;
ll n,mi[maxn],inv[maxn],G[maxn],Ginv[maxn],X[maxn];

char cc; ll ff;
template<typename T>void read(T& aa) {
	aa=0;cc=getchar();ff=1;
	while((cc<'0'||cc>'9')&&cc!='-') cc=getchar();
	if(cc=='-') ff=-1,cc=getchar();
	while(cc>='0'&&cc<='9') aa=aa*10+cc-'0',cc=getchar();
	aa*=ff;
}

void debug() {
    printf("Ginv:
");
    For(i,0,n) printf("%lld ",Ginv[i]);
    printf("
");
}

ll qp(ll x,ll k) {
	ll rs=1;
	while(k) {
		if(k&1) rs=rs*x%mod;
		k>>=1; x=x*x%mod;
	}
	return rs;
}

ll qp1(ll x,ll k) {
	if(k<0) return qp(qp(x,mod-2),-k);
	else return qp(x,k);
}

void Rader(ll F[],ll len) {
	for(int i=1,j=len/2,k;i<len-2;++i) {
		if(i<=j) swap(F[i],F[j]);
		k=len>>1;
		while(j>=k) {j-=k;k>>=1;}
		if(j<k) j+=k;
	}
}

void FFT(ll F[],ll len,ll on) {
	Rader(F,len);
	for(int h=1;h<=len;h<<=1) {
		ll wn=qp1(B,(mod-1)*on/h);
		for(int j=0;j<len;j+=h) {
			ll w=1;
			for(int i=j;i<j+h/2;++i) {
				ll u=F[i],v=w*F[i+h/2]%mod;
				F[i]=(u+v)%mod;
				F[i+h/2]=(u-v+mod)%mod;
				w=w*wn%mod;
			}
		}
	}
	ll x=qp1(len,mod-2);
	if(on==-1) For(i,0,len) F[i]=F[i]*x%mod;
}

void get_Ginv(ll N) {
	if(N==1) {
		Ginv[0]=qp(G[0],mod-2);
		return;
	}
	get_Ginv((N+1)>>1);
	ll n=1;for(;n<=(N<<1);n<<=1);
	For(i,0,N-1) X[i]=G[i]; For(i,N,n) X[i]=0;
	FFT(X,n,1); FFT(Ginv,n,1);
	For(i,0,n) Ginv[i]=Ginv[i]*(2-Ginv[i]*X[i]%mod+mod)%mod;
	FFT(Ginv,n,-1);
	For(i,N,n) Ginv[i]=0;
//	printf("N=%lld:
",N);debug();
}

int main() {
	read(n); 
	mi[0]=1; For(i,1,n) mi[i]=mi[i-1]*i%mod;
	inv[n]=qp(mi[n],mod-2);
	Rep(i,n,1) inv[i-1]=inv[i]*i%mod;
	G[0]=1; For(i,1,n) G[i]=(mod-2*inv[i]%mod)%mod;
//	For(i,0,n) printf("%lld ",G[i]);
//	printf("
");
	get_Ginv(n+1);
	ll ans=0;
//	debug();
	For(i,0,n) ans+=Ginv[i]*mi[i]%mod;
	printf("%lld
",ans%mod);
	return 0;
}

bzoj1406密码箱

枚举一下>=sqrt(n)的因子，然后暴力往set里面插就可以了

下午改了一下午上午看了一上午都没有de出bug来的t2，心情特别不爽。

写这个没有什么脑子的题都三番五次的WA。

//Serene
#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdlib>
#include<cstdio>
#include<cmath>
#include<set>
using namespace std;
#define ll long long
#define db double
#define For(i,a,b) for(int i=(a);i<=(b);++i)
#define Rep(i,a,b) for(int i=(a);i>=(b);--i)
ll n;
set<ll> G;

char cc; ll ff;
template<typename T>void read(T& aa) {
	aa=0;cc=getchar();ff=1;
	while((cc<'0'||cc>'9')&&cc!='-') cc=getchar();
	if(cc=='-') ff=-1,cc=getchar();
	while(cc>='0'&&cc<='9') aa=aa*10+cc-'0',cc=getchar();
	aa*=ff;
}

int main() {
	read(n);ll x,y,t=sqrt(n);
	for(ll i=1;i<=t;++i) if(n%i==0) {
		x=n/i;
		for(y=x;y<=n;y+=x) {
			if(y*(y-2)%n==0) G.insert(y-1);
			if(y*(y+2)%n==0) G.insert((y+1)%n);
		}
	}
	if(G.empty()) printf("None
");
	while(!G.empty()) {
		printf("%lld
",*G.begin());
		G.erase(G.begin());
	}
	return 0;
}

　　

bzoj4173 数学

 N,M<=10^15

一眼看去不可做

如何能神奇地发现$k$当且仅当满足$frac{n+m}{k} - frac{n}{k} - frac{m}{k}$的时候$m mod k + n mod k geq k$？

$frac{n+m}{k} = frac{n}{k} + frac{m}{k} + frac{n mod k + m mod k}{k}$

于是，我们需要求的东西转化成

$varphi(n) imes varphi(m) imes ( sum varphi(k) imes lfloor frac{n+m}{k} floor -sum varphi(k) imes lfloor frac{n}{k} floor -sum varphi(k) imes lfloor frac{m}{k} floor )$

如何才能神奇地发现$sum varphi(k) imes lfloor frac{n}{k} floor = sumlimits_{i=1}^{n} i$？

$sumlimits_{i=1}^{n} i = sumlimits_{i=1}^{n} sumlimits_{d|i} varphi(d)$

ARC076F - Exhausted?

hall定理

反过来，对于一个L,R，有一些人是必须在[1,L]或[R,m]里面的，用人数-区间并的长度的max来更新答案。

如果枚举L、R复杂度是n^2的，所以把l放到线段树里面。

对于L>=R的情况，就找到最多的人-m，那么就是n-m

//Serene
#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdlib>
#include<cstdio>
#include<cmath>
using namespace std;
#define ll long long
#define db double
#define For(i,a,b) for(int i=(a);i<=(b);++i)
#define Rep(i,a,b) for(int i=(a);i>=(b);--i)
const int maxn=2e5+7;
int n,m,p,ans;

char cc;ll ff;
template<typename T>void read(T& aa) {
	aa=0;ff=1; cc=getchar();
	while(cc!='-'&&(cc<'0'||cc>'9')) cc=getchar();
	if(cc=='-') ff=-1,cc=getchar();
	while(cc>='0'&&cc<='9') aa=aa*10+cc-'0',cc=getchar();
	aa*=ff;
}

struct Node{
	int l,r;
	bool operator < (const Node& b) const{return r>b.r;}
}node[maxn];

int num[4*maxn],laz[4*maxn],ql,qr;

int pd(int pos,int l,int r) {
	int mid=(l+r)>>1;
	if(!laz[pos]) return mid;
	laz[pos<<1]+=laz[pos];
	laz[pos<<1|1]+=laz[pos];
	num[pos<<1]+=laz[pos];
	num[pos<<1|1]+=laz[pos];
	laz[pos]=0;
	return mid;
}

void bld(int pos,int l,int r) {
	if(l==r) {
		num[pos]=-l;
		return;
	}
	int mid=(l+r)>>1;
	bld(pos<<1,l,mid);
	bld(pos<<1|1,mid+1,r);
	num[pos]=max(num[pos<<1],num[pos<<1|1]);
}

void chge(int pos,int l,int r) {
	if(l>=ql&&r<=qr) {
		num[pos]++;
		laz[pos]++;
		return;
	}
	int mid=pd(pos,l,r);
	if(ql<=mid) chge(pos<<1,l,mid);
	if(qr>mid) chge(pos<<1|1,mid+1,r);
	num[pos]=max(num[pos<<1],num[pos<<1|1]);
}

int main() {
	read(n); read(m);
	For(i,1,n) {
		read(node[i].l);
		read(node[i].r);
		if(node[i].l==0&&node[i].r>m) p++,i--,n--;
	}
	sort(node+1,node+n+1);
	bld(1,1,m+1);
	For(i,1,n) {
		ql=node[i].l+1; qr=m+1;
		chge(1,1,m+1);
		ans=max(ans,num[1]+1-(m-node[i].r+1));
	}
	printf("%d",max(ans,n-m)+p);
	return 0;
}

bzoj4259残缺的字符串

很久很久以前，在你刚刚学习字符串匹配的时候，有两个仅包含小写字母的字符串A和B，其中A串长度为m，B串长度为n。

可当你现在再次碰到这两个串时，这两个串已经老化了，每个串都有不同程度的残缺。

你想对这两个串重新进行匹配，其中A为模板串，那么现在问题来了，请回答，对于B的每一个位置i，从这个位置开始连续m个字符形成的子串是否可能与A串完全匹配?

1<=m<=n<=300000

FFT

我们考虑构造一个函数，可以根据函数值区分匹配和不匹配。

把'*'的位置先设为0

$dis(A,B)=sumlimits_{i=0}^{n-1}(A[i]-B[i])^2A[i]B[i]$

当dis为0的时候两个匹配

明显，把A翻转之后就是卷积形式，可以把$(A[i]-B[i])^2A[i]B[i]$拆成$A[i]^3B[i]-2A[i]^2B[i]^2+A[i]B[i]^3$

//Serene
#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdlib>
#include<cstdio>
#include<cmath>
using namespace std;
#define ll long long
#define db double
#define For(i,a,b) for(int i=(a);i<=(b);++i)
#define Rep(i,a,b) for(int i=(a);i>=(b);--i)
const int maxn=11e5+7;
const db PI=acos(-1);
int n,m,len,A[maxn],B[maxn],ans;
char s[maxn];
db f[maxn];

char cc; ll ff;
template<typename T>void read(T& aa) {
	aa=0;cc=getchar();ff=1;
	while((cc<'0'||cc>'9')&&cc!='-') cc=getchar();
	if(cc=='-') ff=-1,cc=getchar();
	while(cc>='0'&&cc<='9') aa=aa*10+cc-'0',cc=getchar();
	aa*=ff;
}

struct Virt{
	db r,i;
	Virt(db r=0.,db i=0.){this->r=r; this->i=i;}
	Virt operator + (const Virt& b) const{return Virt(r+b.r,i+b.i);}
	Virt operator - (const Virt& b) const{return Virt(r-b.r,i-b.i);}
	Virt operator * (const Virt& b) const{return Virt(r*b.r-i*b.i,r*b.i+i*b.r);}
}a[maxn],b[maxn];

void Rader(Virt F[],int len) {
	for(int i=1,j=len/2,k;i<len-1;++i) {
		if(i<j) swap(F[i],F[j]);
		k=len>>1;
		while(j>=k) { j-=k; k>>=1; }
		if(j<k) j+=k;
	}
}

void FFT(Virt F[],int len,int on) {
	Rader(F,len);
	for(int h=2;h<=len;h<<=1) {
		Virt wn(cos(-on*2*PI/h),sin(-on*2*PI/h));
		for(int j=0;j<len;j+=h) {
			Virt w(1,0);
			for(int k=j;k<j+h/2;++k) {
				Virt u=F[k],v=w*F[k+h/2];
				F[k]=u+v;
				F[k+h/2]=u-v;
				w=w*wn;
			}
		}
	}
	if(on==-1) For(i,0,len-1) F[i].r/=len;
}

void get_FFT(db x) {
	FFT(a,len,1);
	FFT(b,len,1);
	For(i,0,len) a[i]=a[i]*b[i];
	FFT(a,len,-1);
	For(i,0,len) f[i]+=x*a[i].r;
}

int main() {
	read(m); read(n);
	scanf("%s",s+1);
	For(i,1,m) if(s[i]!='*')
		A[m-i]=s[i]-'a'+1;
	scanf("%s",s+1);
	For(i,1,n) if(s[i]!='*')
		B[i-1]=s[i]-'a'+1;
	for(len=1;len<=n+m;len<<=1);
	For(i,0,len) {
		a[i]=Virt(A[i]*A[i]*A[i],0);
		b[i]=Virt(B[i],0);
	}
	get_FFT(1);
	For(i,0,len) {
		a[i]=Virt(A[i]*A[i],0);
		b[i]=Virt(B[i]*B[i],0);
	}
	get_FFT(-2);
	For(i,0,len) {
		a[i]=Virt(A[i],0);
		b[i]=Virt(B[i]*B[i]*B[i],0);
	}
	get_FFT(1);
	For(i,m-1,n-1) if((int)(f[i]+0.49)==0) ans++;
	printf("%d
",ans);
	For(i,m-1,n-1) if((int)(f[i]+0.49)==0) printf("%d ",i-m+2);
	printf("
");
	return 0;
}

bzoj2005能量采集

栋栋有一块长方形的地，他在地上种了一种能量植物，这种植物可以采集太阳光的能量。

在这些植物采集能量后，栋栋再使用一个能量汇集机器把这些植物采集到的能量汇集到一起。

栋栋的植物种得非常整齐，一共有n列，每列有m棵，植物的横竖间距都一样，

因此对于每一棵植物，栋栋可以用一个坐标(x, y)来表示，其中x的范围是1至n，表示是在第x列，y的范围是1至m，表示是在第x列的第y棵。

由于能量汇集机器较大，不便移动，栋栋将它放在了一个角上，坐标正好是(0, 0)。 能量汇集机器在汇集的过程中有一定的能量损失。

如果一棵植物与能量汇集机器连接而成的线段上有k棵植物，则能量的损失为2k + 1。

例如，当能量汇集机器收集坐标为(2, 4)的植物时，由于连接线段上存在一棵植物(1, 2)，会产生3的能量损失。

注意，如果一棵植物与能量汇集机器连接的线段上没有植物，则能量损失为1。

现在要计算总的能量损失。

1 ≤ n, m ≤ 100,000。

所有我能独立做出来的数学题都是大水题

我们知道，令一个位置(a,b)到(0,0)之间的植物数量(如果包括(a,b)这棵植物的话)是f(a,b)

>令$r=gcd(a,b)$，$x=frac{a}{r}$，$y=frac{b}{r}$

那么$f(a,b)=frac{a}{x}=frac{b}{x}$

也就是$f(a,b)=frac{a}{frac{a}{r}}=r$

问题转化成求$sumlimits_{i=1}^{n}sumlimits_{j=1}^{m}gcd(i,j)$

设$n leq m$

套路推式子……

那么$ans=sumlimits_{T=1}^{n} lfloor frac{n}{T} floor lfloor frac{m}{T} floor sumlimits_{d|T} frac{T}{d} mu (d)$

令$g(n)=sumlimits_{T=1}^{n} sumlimits_{d|T} frac{T}{d} mu (d) = sumlimits_{d=1}^{n} mu (d) sumlimits_{i=1}^{frac{n}{d}}i$

//Serene
#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdlib>
#include<cstdio>
#include<cmath>
using namespace std;
#define ll long long
#define db double
#define For(i,a,b) for(int i=(a);i<=(b);++i)
#define Rep(i,a,b) for(int i=(a);i>=(b);--i)
const int maxn=1e5+7;
ll n,m,g[maxn];

char cc; ll ff;
template<typename T>void read(T& aa) {
	aa=0;cc=getchar();ff=1;
	while((cc<'0'||cc>'9')&&cc!='-') cc=getchar();
	if(cc=='-') ff=-1,cc=getchar();
	while(cc>='0'&&cc<='9') aa=aa*10+cc-'0',cc=getchar();
	aa*=ff;
}

int prime[maxn],mu[maxn],totp;
bool ok[maxn];
void get_p() {
	mu[1]=1;
	For(i,2,n) {
		if(!ok[i]) {
			mu[i]=-1;
			prime[++totp]=i;
		}
		For(j,1,totp) {	
			if(i*prime[j]>n) break;
			ok[i*prime[j]]=1;
			if(i%prime[j]) mu[i*prime[j]]=-mu[i];
			else {
				mu[i*prime[j]]=0;
				break;
			}
		}
	}
	For(i,2,n) mu[i]+=mu[i-1];
}

ll get_g(ll n) {
	if(~g[n]) return g[n];
	ll rs=0;
	for(ll i=1,j;i<=n;i=j+1) {
		j=n/(n/i);
		rs+=(mu[j]-mu[i-1])*(n/i)*(n/i+1)/2;
	}
	return g[n]=rs;
}

ll solve() {
	ll rs=0;
	for(ll i=1,j;i<=n;i=j+1) {
		j=min(n/(n/i),m/(m/i));
		rs+=(n/i)*(m/i)*(get_g(j)-get_g(i-1));
	}
	return rs;
}

int main() {
	read(n); read(m);
	if(n>m) swap(n,m);
	get_p();
	memset(g,-1,sizeof(g)); g[0]=0;
	printf("%lld
",2*solve()-n*m);
	return 0;
}

bzoj4332 分零食

校长从幼儿园旁边的小吃店买了大量的零食决定分给同学们。

同学们依次排成了一列，其中有A位小朋友，有三个共同的欢乐系数O，S和U。如果有一位小朋友得到了x个糖果，那么她的欢乐程度就是$F(x)=O*x^2+S*x+U$。

现在校长开始分糖果了，一共有M个糖果。有些小朋友可能得不到糖果，对于那些得不到糖果的小朋友来说，欢乐程度就是1。

如果一位小朋友得不到糖果，那么在她身后的小朋友们也都得不到糖果。(即这一列得不到糖果的小朋友一定是最后的连续若干位)

求所有方案下小朋友欢乐程度乘积的总和对P取膜的结果

$M leq 10000 , P leq 255 , A leq 10^8 , O leq 4 , S leq 300 , U leq 100$

终于知道传说中的分治FFT是什么了，感觉什么多项式求逆，多项式……都是分治FFT呐

先给一个xxy大佬的题解的链接：https://www.cnblogs.com/TheRoadToTheGold/p/8960712.html

因为如果一位小朋友得不到糖果，那么在她身后的小朋友们也都得不到糖果。

所以设g[i][j] 表示前i位小朋友，分到j个糖果，且前i位小朋友都分到糖果的方案数

令F(x) 表示分到x个糖果的欢乐程度 ∴g[i][j] = ∑ g[i-1][j-k]*F(k)

记g[i]=g[i-1]*F，则 g[i]=F ^ i

但是要求的是 Σ g[i][m]

记f[n]=Σ g[i] i∈[1,n] ，那么ans=f[n][m]

f[n]=Σ g[i] i∈[1,n]

=Σ f(n/2)+Σ g[i] i∈[n/2+1,n]

=Σ f(n/2)+Σ F^i i∈[n/2+1,n]

=Σ f(n/2)+Σ F^(n/2+i) i∈[1,n/2]

=Σ f(n/2)+F^(n/2) * Σ F^i i∈[1,n/2]

=Σ f(n/2)+g(n/2)*f(n/2)

然后可以分治解决

如果n是奇数，f(n)=f(n-1)+g[n]=f(n-1)+g(n-1)*F

边界条件：g[][0]=1

//Serene
#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdlib>
#include<cstdio>
#include<cmath>
using namespace std;
#define ll long long
#define db double
#define For(i,a,b) for(int i=(a);i<=(b);++i)
#define Rep(i,a,b) for(int i=(a);i>=(b);--i)
const db PI=acos(-1),eps=1e-2;
const int maxn=(1<<17)+7;
ll n,m,len,mod;

char cc;ll ff;
template<typename T>void read(T& aa) {
	aa=0;ff=1; cc=getchar();
	while(cc!='-'&&(cc<'0'||cc>'9')) cc=getchar();
	if(cc=='-') ff=-1,cc=getchar();
	while(cc>='0'&&cc<='9') aa=aa*10+cc-'0',cc=getchar();
	aa*=ff;
}

struct Virt{
	db r,i;
	Virt(db r=0.0,db i=0.0) {this->r=r;this->i=i;}
	Virt operator + (const Virt& b) const{return Virt(r+b.r,i+b.i);}
	Virt operator - (const Virt& b) const{return Virt(r-b.r,i-b.i);}
	Virt operator * (const Virt& b) const{return Virt(r*b.r-i*b.i,r*b.i+i*b.r);}
}F[maxn],f[maxn],g[maxn],X[maxn];

void Rader(Virt F[],ll len) {
	for(int i=1,j=len>>1,k;i<len-1;++i) {
		if(i<j) swap(F[i],F[j]);
		k=len>>1;
		while(j>=k) {j-=k;k>>=1;}
		if(j<k) j+=k;
	}
}

void FFT(Virt F[],ll len,ll on) {
	Rader(F,len);
	for(int h=2;h<=len;h<<=1) {
		Virt wn=Virt(cos(-on*2*PI/h),sin(-on*2*PI/h));
		for(int j=0;j<len;j+=h) {
			Virt w=Virt(1,0);
			for(int i=j;i<j+h/2;++i) {
				Virt u=F[i],v=w*F[i+h/2];
				F[i]=u+v;
				F[i+h/2]=u-v;
				w=w*wn;
			}
		}
	}
	if(on==-1) For(i,0,len) F[i].r/=len;
}

void debug() {
	printf("g: "); For(i,0,m) printf("%lld ",(ll)g[i].r);printf("
");
	printf("f: "); For(i,0,m) printf("%lld ",(ll)f[i].r);printf("
");
}

void get_mod() {
	For(i,0,m) f[i].r=((ll)(f[i].r+eps))%mod,f[i].i=0;
	For(i,0,m) g[i].r=((ll)(g[i].r+eps))%mod,g[i].i=0;
}

void solve(ll N) {
	if(!N) {g[0].r=1;return;}
	if(N&1) {
		solve(N-1);
		FFT(g,len,1);
		For(i,0,len) g[i]=g[i]*F[i];
		FFT(g,len,-1); For(i,m+1,len) f[i].r=f[i].i=g[i].r=g[i].i=0;
		For(i,0,m) f[i]=f[i]+g[i];
	}
	else {
		solve(N>>1);
		For(i,0,m) X[i]=f[i];
		FFT(X,len,1); FFT(g,len,1);
		For(i,0,len) X[i]=X[i]*g[i];
		For(i,0,len) g[i]=g[i]*g[i];
		FFT(X,len,-1);FFT(g,len,-1);
		For(i,m+1,len) X[i].r=X[i].i=g[i].r=g[i].i=0;
		For(i,0,m) f[i]=f[i]+X[i];
	}
	get_mod();
//	printf("solve %lld:
",N); debug();
}

int main() {
	ll x,y,z;
	read(m); read(mod); read(n);
	read(x); read(y); read(z);
	For(i,1,m) F[i].r=((ll)x*i*i+y*i+z)%mod;
	for(len=1;len<=(m<<1);len<<=1);
	FFT(F,len,1);
	solve(n);
	printf("%lld
",(ll)f[m].r);
	return 0;
}