【模板】数学
数论
筛法
线性筛
//P3383 【模板】线性筛素数
#include <bits/stdc++.h>
using namespace std;
const int N=1e8+1,M=6e6+10;
int prime[M],cnt=0;
bool vis[N];
void get_prime( int n )
{
memset( vis,0,sizeof(vis) ); cnt=0;
for ( int i=2; i<=n; i++ )
{
if ( !vis[i] ) prime[++cnt]=i;
for ( int j=1; j<=cnt && prime[j]*i<=N-1; j++ )
{
vis[prime[j]*i]=1;
if ( i%prime[j]==0 ) break;
}
}
}
int main()
{
int n,q; scanf( "%d%d",&n,&q );
get_prime( n );
while ( q-- )
{
int x; scanf( "%d",&x );
printf( "%d\n",prime[x] );
}
return 0;
}
杜教筛
//Author: RingweEH
//P4213 【模板】杜教筛
const int N=5e5+10;
int pri[N],tot=0,n;
ll mu[N],phi[N];
bool is[N];
void Sieve()
{
mu[1]=phi[1]=1; is[1]=1;
for ( int i=2; i<=N-10; i++ )
{
if ( !is[i] ) { pri[++tot]=i,mu[i]=-1,phi[i]=i-1; }
for ( int j=1; j<=tot && (i*pri[j]<=(N-10)); j++ )
{
is[i*pri[j]]=1;
if ( i%pri[j]==0 ) { mu[i*pri[j]]=0,phi[i*pri[j]]=phi[i]*pri[j]; break; }
mu[i*pri[j]]=-mu[i]; phi[i*pri[j]]=phi[i]*phi[pri[j]];
}
}
for ( int i=2; i<=(N-10); i++ )
mu[i]+=mu[i-1],phi[i]+=phi[i-1];
}
map<int,ll> mp_mu,mp_phi;
ll Sieve_Phi( int n )
{
if ( n<=N-10 ) return phi[n];
if ( mp_phi[n] ) return mp_phi[n];
ll _res=0;
for ( int l=2,r; l<=n; l=r+1 )
{
r=n/(n/l); _res+=(r-l+1)*Sieve_Phi(n/l);
if ( r>=2147483647 ) break;
}
ll res=(ull)n*(n+1ll)/2ll-_res;
mp_phi[n]=res; return res;
}
ll Sieve_Mu( int n )
{
if ( n<=N-10 ) return mu[n];
if ( mp_mu[n] ) return mp_mu[n];
ll res=1;
for ( int l=2,r; l<=n; l=r+1 )
{
r=n/(n/l); res-=(r-l+1)*Sieve_Mu(n/l);
if ( r>=2147483647 ) break;
}
return mp_mu[n]=res;
}
int main()
{
Sieve();
int T=read();
while ( T-- )
{
int n=read();
printf("%lld %lld\n",Sieve_Phi(n),Sieve_Mu(n) );
}
return 0;
}
数论
BSGS
#include <bits/stdc++.h>
using namespace std;
int p,x,y;
map<int,int> M;
int power( int x,int y )
{
int ans=1;
for ( ; y; y>>=1,x=(1ll*x*x)%p )
if ( y&1 ) ans=1ll*ans*x%p;
return ans;
}
int main()
{
scanf( "%d%d%d",&p,&x,&y );
int m=sqrt(p)+1,s=y;
for ( int i=0; i<m; i++ )
M[s]=i,s=1ll*s*x%p;
int pow=power( x,m ); s=1;
for ( int i=1; i<=m; i++ )
{
s=1ll*s*pow%p;
if ( M.count(s) )
{
printf( "%d\n",i*m-M[s] ); return 0;
}
}
printf( "no solution" );
}
exBSGS
//Author: RingweEH
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <map>
#include <cmath>
#define ll long long
#define db double
using namespace std;
int Mod,x,y;
map<int,int> mp;
int power(int a,int b){int res=1;for(;b;b>>=1,a=1ll*a*a%Mod)if(b&1)res=1ll*res*a%Mod;return res;}
int gcd(int a,int b){ return (!b) ? a : gcd(b,a%b); }
void exBSGS( int x,int y )
{
if ( y==1 ) { puts("0"); return ; }
int d=gcd(x,Mod),k=1,t=0;
while ( d^1 )
{
if ( y%d ) { puts("No Solution"); return; }
t++; y/=d; Mod/=d; k=(x/d)*1ll*k%Mod;
if ( y==k ) { printf("%d\n",t ); return; }
d=gcd(x,Mod);
}
int s=y; mp.clear(); int m=sqrt(Mod)+1;
for ( int i=0; i<m; i++ ) mp[s]=i,s=1ll*s*x%Mod;
s=k; k=power(x,m);
for ( int i=1; i<=m; i++ )
{
s=1ll*s*k%Mod;
if ( mp[s] ) { printf("%d\n",i*m-mp[s]+t ); return; }
}
puts("No Solution");
}
int main()
{
//freopen( "exam.in","r",stdin );
while ( scanf( "%d%d%d",&x,&Mod,&y ) && (x||Mod||y) )
x%=Mod,y%=Mod,exBSGS(x,y);
return 0;
}
exCRT
#include <cstdio>
#include <cmath>
#define ll long long
using namespace std;
int n;
ll exgcd( ll a, ll b, ll &x, ll &y )
{
if ( b==0 ) { x=1,y=0; return a; }
ll d=exgcd( b,a%b,y,x ); y-=a/b*x;
return d;
}
ll inline mod( ll a,ll b )
{
return ((a % b)+b)%b;
}
int main()
{
scanf( "%d",&n );
ll a1,m1; scanf( "%lld%lld",&a1,&m1 );
for ( int i=1; i<n; i++ )
{
ll a2,m2,k1,k2; scanf( "%lld%lld",&a2,&m2 );
ll d=exgcd( a1,-a2,k1,k2 );
if ( (m2-m1)%d )
{ printf( "-1\n" ); return 0; }
k1=mod( k1*(m2-m1)/d,abs( a2/d ) );
m1=k1*a1+m1; a1=abs( a1/d*a2 );
}
printf( "%lld\n",m1 );
return 0;
}
exLucas
#include <bits/stdc++.h>
#define ll long long
using namespace std;
const int N=11;
ll y,z,p,w[N];
ll power( ll a,ll b,ll p )
{
ll res=1;
for ( ; b; b>>=1,a=a*a%p )
if ( b&1 ) res=res*a%p;
return res;
}
ll exgcd( ll a,ll b,ll &x,ll &y )
{
if ( !b ) { x=1; y=0; return a; }
ll res=exgcd( b,a%b,x,y );
ll t=x; x=y; y=t-a/b*y; return res;
}
ll fac( ll n,ll pi,ll pk )
{
if ( !n ) return 1;
ll res=1;
for ( ll i=2; i<=pk; i++ )
if ( i%pi ) (res*=i)%=pk;
res=power( res,n/pk,pk );
for ( ll i=2; i<=n%pk; i++ )
if ( i%pi ) (res*=i)%=pk;
return res*fac(n/pi,pi,pk)%pk;
}
ll inv( ll n,ll mod )
{
ll x,y; exgcd( n,mod,x,y );
return (x+=mod)>mod ? x-mod : x;
}
ll CRT( ll b,ll mod )
{
return b*inv(p/mod,mod)%p*(p/mod)%p;
}
ll C( ll n,ll m,ll pi,ll pk )
{
ll up=fac( n,pi,pk ),t1=fac( m,pi,pk ),t2=fac( n-m,pi,pk ),k=0;
for ( ll i=n; i; i/=pi )
k+=i/pi;
for ( ll i=m; i; i/=pi )
k-=i/pi;
for ( ll i=n-m; i; i/=pi )
k-=i/pi;
return up*inv(t1,pk)%pk*inv(t2,pk)%pk*power(pi,k,pk)%pk;
}
ll exlucas( ll n,ll m )
{
ll res=0,tmp=p,pk;
int lim=sqrt(p)+5;
for ( int i=2; i<=lim; i++ )
if ( tmp%i==0 )
{
pk=1;
while ( tmp%i==0 ) pk*=i,tmp/=i;
res=(res+CRT(C(n,m,i,pk),pk))%p;
}
if ( tmp>1 ) res=(res+CRT(C(n,m,tmp,tmp),tmp))%p;
return res;
}
int main()
{
scanf( "%lld%lld%lld",&y,&z,&p );
printf( "%lld",exlucas(y,z) );
}
二次剩余
//Author: RingweEH
//P5491 【模板】二次剩余
#include <cstdio>
#include <cstring>
#include <algorithm>
#define ll long long
#define db double
using namespace std;
int min( int a,int b ) { return a<b ? a : b; }
int max( int a,int b ) { return a>b ? a : b; }
void bmax( int &a,int b ) { a=(a>b) ? a : b; }
void bmin( int &a,int b ) { a=(a<b) ? a : b; }
int read()
{
int x=0,w=1; char ch=getchar();
while ( ch>'9' || ch<'0' ) { if ( ch=='-' ) w=-1; ch=getchar(); }
while ( ch<='9' && ch>='0' ) x=x*10+ch-'0',ch=getchar();
return x*w;
}
ll Mod,w,n;
struct Complex
{
ll x,y;
Complex( ll _x=0,ll _y=0 ) : x(_x),y(_y) {}
Complex operator * ( Complex t ) { return Complex(x*t.x+y*t.y%Mod*w,x*t.y+y*t.x); }
Complex operator % ( ll Mod ) { return Complex( (x%Mod+Mod)%Mod,(y%Mod+Mod)%Mod ); }
};
ll power( ll a,ll b,ll Mod ) { ll res=1; for ( ; b; b>>=1,a=a*a%Mod ) if ( b&1 ) res=res*a%Mod; return res; }
Complex power( Complex a,ll b,ll Mod ) { Complex res(1,0); for ( ; b; b>>=1,a=a*a%Mod ) if ( b&1 ) res=res*a%Mod; return res; }
ll Solve( ll n,ll P )
{
n%=P;
if ( power(n,(P-1)>>1,P)==P-1 ) return -1; //判断无解
ll a=0;
while ( 1 )
{
a=rand()%P; w=(a*a+P-n)%P;
if ( power(w,(P-1)>>1,P)==P-1 ) break;
}
return power( Complex(a,1),(P+1)>>1,P ).x%P;
}
int main()
{
int T=read();
while ( T-- )
{
n=read(); Mod=read();
if ( n==0 ) { puts("0"); continue; }
ll ans1=Solve(n,Mod),ans2=Mod-ans1;
if ( ans1==-1 ) { puts("Hola!"); continue; }
if ( ans1>ans2 ) swap( ans1,ans2 );
printf( "%lld %lld\n",ans1,ans2 );
}
return 0;
}
线性代数
矩阵和行列式
高斯消元
//P3389 【模板】高斯消元法
#include <bits/stdc++.h>
#define ll long long
using namespace std;
const int N=2010;
int n,cnt;
double a[N][N];
int main()
{
scanf( "%d",&n );
for ( int i=1; i<=n; i++ )
for ( int j=1; j<=n+1; j++ )
scanf( "%lf",&a[i][j] );
for ( int i=1; i<=n; i++ )
{
cnt=i;
while ( a[cnt][i]==0 && cnt<=n ) cnt++;
if ( cnt==n+1 ) { printf( "No Solution" ); return 0; }
for ( int j=1; j<=n+1; j++ )
swap( a[i][j],a[cnt][j] );
double k=a[i][i];
for ( int j=1; j<=n+1; j++ )
a[i][j]/=k;
for ( int j=1; j<=n; j++ )
if ( i!=j )
{
double kt=a[j][i];
for ( int m=1; m<=n+1; m++ )
a[j][m]-=kt*a[i][m];
}
}
for ( int i=1; i<=n; i++ )
printf( "%.2lf\n",a[i][n+1] );
}
矩乘加速
//P1939 【模板】矩阵加速（数列）
#include <bits/stdc++.h>
#define ll long long
using namespace std;
const ll mod=1e9+7;
struct matrix
{
ll mt[5][5];
matrix() { memset(mt,0,sizeof(mt)); }
matrix operator * ( matrix b )
{
matrix c;
for ( int i=1; i<=3; i++ )
for ( int j=1; j<=3; j++ )
for ( int k=1; k<=3; k++ )
c.mt[i][j]=(c.mt[i][j]+(mt[i][k]*b.mt[k][j])%mod)%mod;
return c;
}
}ans,base;
matrix power( ll b )
{
matrix res=ans,a=base;
for ( ; b; b>>=1,a=a*a )
if ( b&1 ) res=res*a;
return res;
}
int main()
{
int T; scanf( "%d",&T );
ans.mt[1][1]=1; ans.mt[1][2]=1; ans.mt[1][3]=1;
base.mt[1][1]=1; base.mt[1][2]=1; base.mt[2][3]=1; base.mt[3][1]=1;
while ( T-- )
{
ll n; scanf( "%lld",&n );
if ( n<=3 ) { printf( "1\n"); continue; }
matrix res=power( n-3 );
printf( "%lld\n",res.mt[1][1] );
}
}
线性基
异或最大值
//Author:RingweEH
//P3812
const int N=55,MX=50;
int n;
ll a[N];
void Insert( ll x )
{
for ( int i=MX; ~i; i-- )
if ( x>>i&1 )
{
if ( !a[i] ) { a[i]=x; return; }
x^=a[i];
}
}
ll Query_mx( ll res=0 )
{
for ( int i=MX; ~i; i-- )
res=max( res,res^a[i] );
return res;
}
int main()
{
n=read();
for ( int i=1; i<=n; i++ )
{
ll x=read(); Insert( x );
}
printf( "%lld\n",Query_mx() );
return 0;
}
异或第 k 小
//Author:RingweEH
const int N=63,M=65;
int n,cnt;
ll a[M],b[M];
bool fl=0;
void Insert( ll x )
{
for ( int i=N; ~i; i-- )
if ( x>>i&1 )
{
if ( a[i] ) x^=a[i];
else
{
a[i]=x;
for ( int j=i-1; ~j; j-- )
if ( a[j] && (a[i]>>j&1) ) a[i]^=a[j];
for ( int j=N; j>i; j-- )
if ( a[j]>>i&1 ) a[j]^=a[i];
return;
}
}
if ( x==0 ) fl=1;
}
int main()
{
int T=read();
for ( int cas=1; cas<=T; cas++ )
{
memset( a,0,sizeof(a) ); fl=0; cnt=0;
n=read();
for ( int i=1; i<=n; i++ )
{
ll x=read(); Insert( x );
}
for ( int i=0; i<=N; i++ )
if ( a[i] ) b[cnt++]=a[i];
int q=read(); printf( "Case #%d:\n",cas );
while ( q-- )
{
ll k=read(),ans=0; k-=fl;
if ( k>=(1ll<<cnt) ) { printf( "-1\n" ); continue; }
for ( int i=0; i<cnt; i++ )
if ( k>>i&1 ) ans^=b[i];
printf( "%lld\n",ans );
}
}
return 0;
}
实数线性基
//Author:RingweEH
//P3265 [JLOI2015]装备购买
const int N=510;
const db eps=1e-5;
struct Vector
{
db a[N]; int val;
db &operator [] ( const int &x ) { return a[x]; }
bool operator < ( const Vector&tmp ) const { return val<tmp.val; }
}c[N];
int n,m,a[N];
int main()
{
n=read(); m=read();
for ( int i=1; i<=n; i++ )
for ( int j=1; j<=m; j++ )
scanf( "%lf\n",&c[i][j] );
for ( int i=1; i<=n; i++ )
c[i].val=read();
sort( c+1,c+1+n ); int cnt=0; ll ans=0;
for ( int i=1; i<=n; i++ )
for ( int j=1; j<=m; j++ )
{
if ( fabs(c[i][j])<eps ) continue;
if ( !a[j] ) { a[j]=i; cnt++; ans+=c[i].val; break; }
db tmp=1.0*c[i][j]/c[a[j]][j];
for ( int k=j; k<=m; k++ )
c[i][k]-=tmp*c[a[j]][k];
}
printf( "%d %lld\n",cnt,ans );
return 0;
}
组合数学
Prufer序列
Prufer序列模板
//Author: RingweEH
const int N=5e6+10;
int n,opt,deg[N],fa[N],prufer[N];
ll ans=0;
void TreeToPrufer()
{
memset( deg,0,sizeof(deg) );
for ( int i=1; i<n; i++ )
fa[i]=read(),deg[fa[i]]++;
for ( int i=1,pos=1; i<=n-2; i++ )
{
while ( deg[pos] ) pos++; //找一个编号最小的叶节点
prufer[i]=fa[pos]; //将其父亲加入序列
for ( ; i<=n-2; )
{
deg[prufer[i]]--;
if ( deg[prufer[i]] ) break; //减完之后不是叶
if ( prufer[i]>=pos ) break; //是叶节点但是编号在后面
prufer[i+1]=fa[prufer[i]]; i++; //新产生的点加入序列
}
pos++;
}
for ( int i=1; i<=n-2; i++ )
ans=ans^(1ll*i*prufer[i]);
printf( "%lld\n",ans );
}
void PruferToTree()
{
memset( deg,0,sizeof(deg) );
for ( int i=1; i<=n-2; i++ )
prufer[i]=read(),deg[prufer[i]]++;
prufer[n-1]=n;
for ( int i=1,pos=1; i<n; i++ )
{
while ( deg[pos] ) pos++; //自增找一个叶节点
fa[pos]=prufer[i]; //与当前序列首相连
for ( ; i<n-1; )
{
deg[prufer[i]]--;
if ( deg[prufer[i]] ) break;
if ( prufer[i]>=pos ) break;
fa[prufer[i]]=prufer[i+1]; i++;
}
pos++;
}
for ( int i=1; i<n; i++ )
ans=ans^(1ll*i*fa[i]);
printf( "%lld\n",ans );
}
int main()
{
n=read(); opt=read();
if ( opt==1 ) TreeToPrufer();
else PruferToTree();
return 0;
}
大地也该是从一片类似的光明中冒出来的。
