ACM 数论基本模板

1.欧几里得

求最大公约数,最小公倍数

(1)递归的写法:int

gcd(int a,int b) {return b?gcd(b,a%b):a;}

(2)辗转相除法:

int gcd(int

a,int b)

{

if(a

int

r;

while(b)

{r=a%b;a=b;b=r;}

return

a;

}

(3)stein+欧几里得

快速求解大数的最大公约数

i64 stein(i64

a,i64 b)

{

if(a

if(b==0)

return a;

if((a&1)==0&&(b&1)==0)

return

2*stein(a>>1,b>>1);//a

and b are even

if((a&1)==0)

return stein(a>>1,b); // only a is even

if((b&1)==0)

return stein(a,b>>1); // only b is even

return stein((a+b)>>1,(a-b)>>1); // a and b are odd

}

最小公倍数: int

lcm(int a,int b) {return a/gcd(a,b)*b;}

2.扩展欧几里得 求ax=b

(mod m) ax+my=b

如果r=gcd(a,m)且b%r==0,则同余方程有解,其最小解为x*(b/r);

ax+by=c

如r=gcd(a,b),则存在x,y,使xa+yb=r;当x+=b,y-=a后仍然成立

因为xa+yb+ab-ab=r;==>(x+b)a+(y-a)b=r

int exgcd(int

a,int b,int &x,int

&y)

{

if(b==0)

{x=1;y=0;return a;}

int

r=exgcd(b,a%b,y,x);

y-=x*(a/b);

return

r;

}

3.素数判定

(1)试除法:

bool isprime(int n)

{

int

i;

for(i=2;i<=(int)sqrt(n*1.0);i++)

if(n%i==0)  return false;

return true;

}

bool isprime(int

n)

{

if(n==2)

return true;

if(n==1||(n&1)==0)

return false;

for(int

i=3;i*i<=n;i+=2)if(n%i==0) return fals;

return

true;

}(2)miller-rabin

算法

bool

witness(i64 a,i64 n)

{

i64

x,d=1,i=ceil(log(n-1.0)/log(2.0))-1;

for(;i>=0;i--)

{

x=d; d=(d*d)%n;

if(d==1&&x!=1&&x!=n-1)

return 1;

if(((n-1)&(1<0)

d=(d*a)%n;

}

return d==1?0:1;

}

bool

miller_rabin(i64 n)

{

if(n==2) return

1;

if(n==1||(n&1)==0) return 0;

i64

j,a;

for(j=0;j<50;j++)

{

a=rand()*(n-2)/RAND_MAX+1;

if(witness(a,n)) return 0;

}

return

1;

}

另一种写法，更好理解

bool witness(i64 a,i64

n)

{

int

i,j=0;

i64

m=n-1,x,y;

while(m%2==0)

{

m>>=1;

j++;

}

x=pow(a,m,n);///快速幂取模

for(i=1;i<=j;i++)

{

y=pow(x,2,n);

if(y==1&&x!=1&&x!=n-1)

return true;

x=y;

}

return

y==1?false:true;

}

bool miller_rabin(i64

n)

{

if(n==2)

return true;

if(n==1||n%2==0)

return false;

for(int

i=1;i<=10;i++)

{

i64

a=rand()%(n-1)+1;

if(witness(a,n)) return false;

}

return

true;

}

4.素数筛法 //前17个素数

prime[18]={17,2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59}

bool

f[100002];//保存判断是否是素数的结果，p[i]=1 是素数，p[i]=0 则不是素数

int

prime[78499];//保存素数prime[0]为素数的个数

void PRIME(int

M)

{

int

i,i2,k;

for(i=0;i<=M;i+=2)

f[i]=0;

for(i=1;i<=M;i+=2)

f[i]=1;

f[1]=0; f[2]=1;

for(i=3;i<=(int)sqrt(1.0*M);i+=2)

if(p[i])

{

i2=i+i; k=i*i;

while(k<=M) {f[k]=0; k+=i2;}

}

prime[1]=2;

k=1;

for(i=3;i<=M;i+2)

if(f[i]) prime[++k]=i;

prime[0]=k;

}

(2)

void PRIME(int

M)

{

int

i,j,k;

prime[1]=2;

prime[2]=3;

for(i=5;i<=M;i+=2)

{

for(j=1;prime[j]*prime[j]<=i;j++)

if(i%prime[j]==0) goto loop;

prime[++k]=i;

loop:;

}

prime[0]=k;

}

5.整数分解

(1)

void split(int

n,int *p,int *t)

{

int

i,s,top=0;

for(i=1;i<=prime[0];i++)

{

s=0;

while(n%prime[i]==0) {s++;n/=prime[i];}

if(s) {p[++top]=prime[i];t[top]=s;}

if(n==1)

break;

if(n

{p[++top]=n;t[top]=1;n=1;break;}

}

p[0]=t[0]=top;

}

(2)分解1-100000的因子,且由prime[n][]保存n的素因子(prime[n][0]为质因子的个数):

void split(int

n)//p[]为素数表

{

int

i,x=n;

prime[n][0]=0;

for(i=1;i<=p[0];i++)if(x%p[i]==0)严重坑爹的bug

{

prime[n][++prime[n][0]]=p[i];

while(x%p[i]==0) x/=p[i];

if(x==1)

break;

}

if(x>1)

prime[n][++prime[n][0]]=x;

}

(3)Pollard-pho大数分解

i64 Pollard(i64 n,int

c)

{

i64

i=1,k=2,x=rand()%n,y=x,d;

srand(time(NULL));

while(true)

{

i++;

x=(mod_mult(x,x,n)+c)%n;

d=gcd(y-x,n);

if(d>1&&d

return d;

if(y==x)

return n;

if(i==k){y=x;k<<=1;}

}

}

6.求因子和与因子个数(包含1和本身)

因子和s是积性函数,即:gcd(a,b)=1==>s(a*b)=s(a)*s(b);

如果p是素数==>s(p^X)=1+p+p^2+...+p^X=(p^(X+1)-1)/(p-1);s(p^2x)=1+p+p^2+...+p^2x=(p^(2x+1)-1)/(p-1);

求因子和:

(1)ans=1+n;

for(i=2;i<=n/2;i++)

if(n%i==0)

{

if(n/i>i)

ans+=i+n/i;

else if(n/i==i) ans+=i;

else break;

}

(2)另一种递推的写法:

for(i=1;i<=lmax;i++) num[i]=1+i;

for(i=2;i<=lmax/2;i++)

for(j=i<<1;j<=lmax;j+=i) num[j]+=i;

求因子个数:

n=p1^t1*p2^t2*p3^t3***pk^tk;

因子个数为:(t1+1)(t2+1)***(tk+1)

for(ret=i=1;i<=prime[0]&&prime[i]<=(int)sqrt(1.0*n);i++)

if(n%prime[i]==0)

{

k=0;

while(n%prime[i]==0)

{k++,n/=prime[i];}

ret*=k+1;

}

if(n>1)

ret*=2;

当求n^2的因子个数的时候:n^2=p1^(2*t1)*p2^(2*t2)***pk^(2*tk); 因子个数为:(2*t1+1)(2*t2+1)***(2*tk+1)

for(ret=i=1;i<=prime[0]&&prime[i]<=(int)sqrt(1.0*n);i++)

if(n%prime[i]==0)

{

k=0;

while(n%prime[i]==0)

{k++,n/=prime[i];}

ret*=2*k+1;

}

if(n>1)

ret*=3;

快速求出一个比较大的区间内的所有因子和:

const int

lmax=50005;//求出[1,50005]区间内每一个数的因子和(不包括本身),并用facsum[]数组保存

i64 facsum[lmax];

for(i=0;i<=lmax;i++) facsum[i]=1;

for(i=2;i*i<=lmax;i++)

{

for(j=i+1;j*i<=lmax;j++) facsum[i*j]+=i+j;

facsum[i*i]+=i;

}

7.欧拉函数

(1)单独求欧拉函数

int eular(int

n)

{

int

ret=1,i;

for(i=2;i*i<=n;i++)

if(n%i==0)

{

n/=i; ret*=i-1;

while(n%i==0) {n/=i;ret*=i;}

}

if(n>1)

ret*=n-1;

return ret;

}

int euler(int x)

{

int i,

res=x,tmp=(int)sqrt(x*1.0)+1;

for(i=2;i

if(x%i==0)

{

res=res/i*(i-1);

while(x%i==

0) x/=i;

}

if(x>1)

res=res/x*(x-1);

return

res;

}

int eular(int

n)

{

int

ret=n,i;

for(i=2;i*i<=n;i++)

if(n%i==0)

{

ret=ret/i*(i-1);

while(n%i==0) n/=i;

}

if(n>1)

ret=ret/n*(n-1);

return

ret;

}

先素数筛法在用欧拉函数(在此仅写其中的一个)

void eular(int

n)

{

int

ret=n,i;

for(i=1;i<=prime[0]&&prime[i]<=(int)sqrt(1.0*n);i++)

if(n%prime[i]==0)

{

ret=ret/prime[i]*(prime[i]-1);

while(n%prime[i]==0) n/=prime[i];

}

if(n>1)

ret=ret/n*(n-1);

return

ret;

}

(2)递推求欧拉函数

const

int lmax=300000;

int PHI(int

lmax)

{

int

i,j;

for(i=1;i<=lmax;i++)phi[i]=i&1?i:i/2;

for(i=3;i<=lmax;i+=2)

if(phi[i])

for(j=i;j<=lmax;j+=i)phi[j]=phi[j]/i*(i-1);

}

(3)同时求出欧拉值和素数

int

prime[lmax][25],num[lmax],eular[lmax];//prime[n][i]表示n的第i+1个素数因子,num[n]表示n的因子个数,eular[n]表示n的欧拉值

///下面这个貌似不对

void

eular_prime()

{

int

i,j;

for(i=1;i<=lmax;i++)

eular[i]=i,prime[i][0]=0;

for(i=2;i<=lmax;i++)

if(eular[i]==i)

{

eular[i]=i-1;

for(j=i<<1;j<=lmax;j+=i)

{

eular[j]=eular[j]/i*(i-1);

prime[j][++prime[j][0]]=i;//prime[j][0]=num[j]表示j的素数因子个数

}

}

}

/上面这个貌似不对

void

eular_prime()//每个数的欧拉函数值及筛选法得到数的素因子num[i]为i的因子个数

{

eular[1]=1;

for(int

i=2;i

{

if(eular[i]==0)

for(int

j=i;j

{

if(eular[j]==0) eular[j]=j;

eular[j]=eular[j]*(i-1)/i;

prime[j][num[j]++]=i;

}

//eular[i]+=eular[i-1];//进行累加(法里数列长度)

}

}

void

eular_prime()

{

int

i,j;

eular[1]=1;

for(i=2;i<=lmax;i++)

if(eular[i]==0)

for(j=i;j<=lmax;j+=i)

{

if(eular[j]==0) eular[j]=j;

eular[j]=eular[j]/i*(i-1);

prime[j][++prime[j][0]]=i;

}

}

欧拉定理的一个重要应用:A^x

mod m=A^(x%phi(m)+phi(m)) mod m (当x>=phi(m)时)

8.求逆元ax=1 (mod

m) x是a的逆元

(1)用扩展欧几里得求

int Inv(int

a,int m)

{

int

r,x,y;

r=exgcd(a,m,x,y);

if(r==1)

return (x%m+m)%m;

return

-1;

}

(2)用快速幂取模求a*a^(p-2)=a^(p-1)=1

(mod p) p必须为素数,a的逆元是a^(p-2);

int pow(int

a,int n)//a^n%p (n=p-2)

{//这里的做法会让a的值变化,可令t=a;用t代替a计算

int

r=1;

while(n)

{

if(n&1) r=r*a%p;

a=a*a%p;

n>>=1;

}

return

r;

}

9.快速模乘

a*b%p

int mul(int

a,int b)

{

int

r=0;

while(b)

{

if(b&1) r=(r+a)%p;

a=(a<<1)%p;

b>>=1;

}

return

r;

}

10.求解模线性方程组(中国剩余定理)x=a1 mod

m1

x=a2 mod

m2

......

x=an mod mn其中,a[],m[]已知,m[i]>0且m[i]与m[j]互质,求x.

设m1,m2,...,mn是两两互素的正数,则对任意的整数a1,a2,...,an,同余方程组其解为:X=((M_1*M1*a1)+(M_2*M2*a2)+...+(M_n*Mn*an))

mod m;

其中m=m1*m2*...*mn; Mi=m/mi; M_i是Mi的逆元

int china(int *a,int *m,int n)

{

int M=1,ans=0,mi,i,x,y;

for(i=0;i

M*=m[i];

for(i=0;i

{

mi=M/m[i];

exgcd(m[i],mi,x,y);

ans=(ans+mi*y*a[i])%M;

}

return (ans%M+M)%M;

}