[ACM]数论模版

欢迎访问我的主页www.hh-yzm.com

最大公约数(GCD)/最小公倍数(LCM)

/*非递归版本求最大公约数*/
int gcd(int a,int b)
{
	if(0==b)	return a;
	while(b>0)
	{
		int temp=a%b;
		a=b;
		b=temp;
	}
	return a;
	//while(b^=a^=b^=a%=b)
}
/*递归版本求最大公约数*/
int gcd(int a,int b)
{
	return b?gcd(b,a%b):a;
}
/*求最小公倍数*/
int lcm(int a,int b)
{
	return a*b/gcd(a,b);
}

素数判断及打表

/*判断n为素数返回1,合数返回0*/
int IsPrime(int n)
{
	if(n==2)	return 1;
	if(n%2==0||n<2)	return 0;
	int l=sqrt(n+1);
	for(int i=3;i<=l;i+=2)
		if(n%i==0)	return 0;
	return 1;	
}

/*素数表打表 欧拉筛法*/
/*primes[]  [2,N] 之间的素数 primes[0]第0位素数2*/
/*isprime[] */
const int maxn=1000007;
int primes[maxn];
int isprime[maxn];

void euler_sieve()
{
	int tot=0;
	memset(isprime,1,sizeof(isprime));
	isprime[0]=isprime[1]=0;
	for(int i=2;i<=maxn;i++)
	{
		if(isprime[i]) primes[tot++]=i;
		for(int j=0;j<tot;j++)
		{
			if(i*primes[j]>maxn) break;
			isprime[i*primes[j]]=0;
			if(i%primes[j]==0) break;
		}
	}
}

快速幂/乘取模

/*快速幂取模 (a^b)%p*/
int quickpow(int a,int b,int p)
{
	int ret=1;
	a%=p;
	while(b)
	{
		if(b&1)	ret=(a*ret)%p;
		a=(a*a)%p;
		b>>=1;
	}
	return ret;
}
/*快速乘法取模 (a*b)%p*/
int quickmul(int a,int b,int p)
{
	int ans=0;
	while(b)
	{
		if(b&1)	{b--;ans=(ans+a)%p;}
		b>>=1;
		a=(a+a)%p;
	}
	return ans;
}

拓展欧几里得

/*拓展欧几里德 ax+by=1*/
int exgcd(int a,int b,int &x,int &y)
{
	if(b==0)
	{
		x=1;
		y=0;
		return a;
	}
	int ans=exgcd(b,a%b,y,x);
	y-=a/b*x;
	return ans;
}

求乘法逆元

/*求乘法逆元(数论倒数) ax=1(modb) -> ax-1=by -> ax+by=1*/ 
int modinverse(int a,int b)
{
	int x,y;
	int d=exgcd(a,b,x,y);
	return d==1?(x%b+b)%b:-1;
}

中国剩余定理(解一元线性同余方程组)

/*中国剩余定理(互质) 一元线性同余方程组 x=a(modm)*/
/*设ai,Mi为除该数外的摸数乘积,ti为Mi模mi数论倒数*/
/*ans(ai*ti*Mi)*/
/*a[],p[]模,n,长度*/
int crt(int a[],int p[],int n)
{
	int muls=1;
	int ret=0;
	for(int i=0;i<n;i++)	muls*=p[i];
	for(int i=0;i<n;i++)
	{
		int x,y;
		int mi=muls/p[i];
		exgcd(mi,p[i],x,y);
		ret=(ret+x*a[i]*mi)%muls;
	}
	return (ret+muls)%muls;
}
/*非互质,无解返回-1*/
int crt(int a[],int p[],int n)
{
	if(n==1)
	{
		return p[0]>a[0]?a[0]:-1;
	}
	int x,y,d;
	for(int i=1;i<n;i++)
	{
		if(p[i]<=a[i])	return -1;
		d=exgcd(p[0],p[i],x,y);
		if((a[i]-a[0])%d!=0)	return -1;
		int t=p[i]/d;
		x=((a[i]-a[0])/d*x%t+t)%t;		
		a[0]=x*p[0]+a[0];
		p[0]=p[0]*p[i]/d;
		a[0]=(a[0]%p[0]+p[0])%p[0];
	}
	return a[0];
}

求解方程ax=b(modn)

/*求解ax=b(modn), 解的个数为gcd(a,n)，返回x为vector<>*/
verctor<int> lmodeq(int a,int b,int n)
{
	int x,y;
	int d=exgcd(a,n,x,y);
	vector<int> ans;
	ans.clear();
	if(b%d==0)
	{
		x=(x%n+n)%n;
		x%=(n/d);
		ans.push_back(x*(b/d)%(n/d));
		for(int i=1;i<d;i++)
			ans.push_back((ans[0]+i*n/d)%n);
	}
	return ans;
}

莫比乌斯函数求解

/*求莫比乌斯函数mu
 mu={
 1	     ,u=1
 (-1)^k  ,能分解成k个不同的质因数因子
 0       ,u包含平方因子
 }
 */
int Mobius(int n)
{
	int cnt,k=0;
	for(int i=2;i*i<=n;i++)
	{
		if(n%i) continue;
		cnt=0;
		k++;
		while(n%i==0)
		{
			n/=i;
			cnt++;
		}
		if(cnt>=2)	return 0;
	}
	if(n)	k++;
	return k%2?-1:1;
}

组合数求解(逆元求组合数/Lucas定理)

const int MOD=1e9+7;
/*利用逆元求组合数*/
DataType exgcd(DataType a,DataType b,DataType &x,DataType &y){
    if(b==0){
        x=1;y=0;return a;
    }
    DataType ans=exgcd(b,a%b,y,x);
    y-=a/b*x;
    return ans;
}
/*计算阶乘*/
DataType fac(DataType n){
    DataType sum=1;
    for(int i=1;i<=n;i++){
        sum=(sum*(DataType)i)%MOD;
    }
    return sum;
}
/*逆元求组合数*/
DataType comb(DataType n,DataType m){
    DataType a=fac(m)*fac(n-m)%MOD;
    DataType x,y;
    exgcd(a,MOD,x,y);
    return ((fac(n)*x)%MOD+MOD)%MOD;
}
/*Lucas定理,运用于MOD数较小时*/
DataType lucas(DataType n,DataType m){
    if(m==0) return 1;
    return comb(n%MOD,m%MOD)*lucas(n/MOD,m/MOD)%MOD;
}