数学相关模板

Posted 2020-10-21 shixinyi

矩阵树定理(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\n",Guess(n<<2));
    }
    return 0;
}

 

Miller-Rabin + Pollard Rho(poj1811)

#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdlib>
#include<cstdio>
#include<cmath>
using namespace std;
#define ll long long
const int TM=25,maxn=5500;
const ll INF=1e17;
ll T,n,cnt,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;
}

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

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

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

bool isprime(ll n) {
	if(n==2||n==3||n==5||n==7) return 1;
	if(n<2||(n%6!=1&&n%6!=5)) return 0;
	ll t=0,a,m=n-1,x,y;
	while((m&1)==0) t++,m>>=1;
	for(int i=0;i<TM;++i) {
		a=rand()%(n-2)+2;
		x=qp(a,m,n); 
		for(int j=0;j<t;++j) {
			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 c) {
	ll i=1,k=2,d,x=rand()%(n-1)+1;
	ll y=x;
	while(1) {
		i++;
		x=(qc(x,x,n)+c)%n;
		d=gcd((y-x+n)%n,n);
		if(1<d&&d<n) return d;
		if(y==x) return n;
		if(i==k) {
			y=x; k<<=1;
		}
	}
}

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

int main() {
	read(T);
	while(T--) {
		read(n); ans=INF;
		find(n,120);
		if(ans==n) printf("Prime\n");
		else printf("%lld\n",ans);
	}
	return 0;
}
/*
100
138
5
10
*/

 

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）

#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdlib>
#include<cstdio>
#include<cmath>
using namespace std;
#define ll long long
const int maxn=50000+7;
const ll mod=1e9+7;
ll T,n,k;

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

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

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

ll solve() {
	y[0]=0;
	for(int i=1;i<k+2;++i) {
		y[i]=y[i-1]+qp(i,k);
		mo(y[i]); 
	}
	if(n<k+2) return y[n];
	ll rs=0,t=1,nowf,tot=0;
	for(int i=0;i<k+2;++i) 
	if(n%mod-i) t=t*(n%mod-i+mod)%mod;
	else tot++;
	if(tot>1) return 0;
	for(int i=1;i<k+2;++i) {
		if(!tot) nowf=y[i]*t%mod*qp(n-i,mod-2)%mod;
		else if(n%mod-i==0) nowf=y[i]*t%mod;
		else continue;
		nowf=nowf*mi[i]%mod;
		if(i<k+1) nowf=nowf*mi[k+1-i]%mod;
		if((k+1-i)&1) nowf=nowf*(mod-1)%mod;
		rs+=nowf; mo(rs);
	}
	return rs;
}

int main() {
	read(T);
	mi[1]=1;
	for(int i=2;i<=50000+3;++i) mi[i]=mi[i-1]*qp(i,mod-2)%mod;
	while(T--) {
		read(n); read(k);
		printf("%lld\n",solve());
	}
	return 0;
}

 

杜教筛留坑