ACM数论模版


在 ACM 程序设计大赛中, 数论占了一个很大的比例, 刚好博主当时在队伍里也是负责数论这个方便, 所以整理了一个数论模版分享给大家, 这个数论模版在博主已经翻烂了, 大家也可以把它当做一个数论总结大纲来看.
Waring : 比赛中模版固然很重要, 但是运用模版的同时, 大家也别忘记了将数论的知识理解透彻, 毕竟学到的都是自己的, 而且只有理解了, 才能灵活的运用好模版.

素数

素数的筛选:

一般用于:

  1. 为了求得某范围内的素数.
  2. 对求素质的因子加速.

用到的定理:

  1. 唯一分解定理: 任一自然数 N 一定可以唯一标示为素数之积. $$N={b_1}^{a_1} + {b_2}^{a_2} + … + {b_k}^{a_k}$$
  2. 如果 $$p^a || a!$$ 那么 $$a=[{n\over p^1}] + [{n\over p^2}] + [{n\over p^3}] + ……$$
  3. 整数 n 的正因数个数称为除法函数, 若 n 的标准分解式为: $$N={b_1}^{a_1} + {b_2}^{a_2} + … + {b_k}^{a_k}$$ 则有: $$d(n) = (a_1 + 1) * (a_2 + 1) …… (a_k + 1)$$

筛选素数的目的:

  1. 真的为了求得某范围内的素数个数:
  2. 对求数字的因子加速, 枚举素数列表.

素数的判断法:

  1. 暴力求解 O(sqrt(n)) .
  2. 预先打表判断.

素数筛选法模版:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
/* 求1-L内的所有素数, 结果保存在prime[L]中 */

const int L=1000005,inf=1<<30,maxn=1005;
int prime[L]; // 保存已经筛选过的素数.
bool is[L]; // 记录是否为素数.
void getPrime()
{
fill(is,is+L,1);
is[1]=0;
int np=0;
for(int i=2;i<L;i++)
if(is[i])
{
prime[++np]=i;
for(int j=2*i;j<L;j+=i) is[j]=0;
}
}

除法函数模板:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
/* 求N得因子个数, 以及最小因子. */

#include<stdio.h>
#include<string.h>
#define CC(m ,what) memset(m , what, sizeof(m))
int N ;
const int MAXN = 1010 ;
bool is_p[MAXN] ;
int p[MAXN] , cnt ;

void calc(){
cnt = 0 ;
for(int i=1;i<MAXN;i++) is_p[i] = 1 ;
is_p[1] = 0 ;
for(long long i=2;i<MAXN;i++){
if( is_p[i] == 0 ) continue ;
p[cnt++] = i ;
for(long long j=i;j*i<MAXN;j++){
is_p[i*j] = 0 ;
}
}
}

int main(){
calc() ;
while(scanf("%d",&N) == 1){
int num = N , res = 1 ;
int n = N ;
for(int i=0;i<cnt && p[i]*p[i]<=n;i++){
if( n % p[i] == 0 ){
int c = 0 ;
if( p[i] < num) num = p[i] ;
while( n % p[i] == 0 ){
c++ ;
n /= p[i] ;
}
res *= (c+1) ;
}
}
if(n > 1){
res *= 2 ;
}
res --;
printf("%d %d\n",res, N/num) ;
}
return 0 ;
}

模方程问题:

一般利用扩展欧几里得算法展开得到.

扩展的欧几里德算法模版:

1
2
3
4
void E_gcd(LL a, LL b, LL &d, LL &x, LL &y){
if(!b){ d = a; x = 1, y = 0 }
else { E_gcd(b, a % b, d, y, x); y -= x * (a / b); }
}

1 ) 求 ax + by = c 的整数解:

ax + by = gcd(a,b) = g 的一组解为 (x0, y0) , 则 ax + by = c 的一组解为 (x0 * c / g, y0 * c / g) 。当c不是g的倍数时无整数解
(x1, y1)ax + by = c 的一组解,则其任意整数解为 (x1 + k * bb, y1 - k * aa) ,其中 aa = a / gcd(a, b), bb = b / gcd(bb) ,k为任意整数.

模版:

1
2
3
4
5
6
7
8
9
bool linear_equation(int a,int b,int c,int &x,int &y)
{
int d=exgcd(a,b,x,y);
if(c%d) // 判断是否有解
return false;
int k=c/d;
x*=k; y*=k; // 求得的只是其中一组解
return true;
}

2 ) 用扩展欧几里德算法求解模线性方程:

同余方程 ax≡b (mod n) 对于未知数 x 有解,当且仅当 gcd(a,n) | b 。且方程有解时,方程有 gcd(a,n) 个解。

求解方程 ax≡b (mod n) 相当于求解方程 ax + ny = b , (x, y为整数)

ans = x * (b / d), s = n / d;
方程 ax≡b (mod n) 的最小整数解为:(ans % s + s) % s;

模版:

1
2
3
4
5
6
7
8
9
10
11
bool modular_linear_equation(int a, int b, int n)
{
int x, y, x0, i;
int d = exgcd(a, n, x, y);
if(b % d)
return false;
x0 = x * (b / d) % n; //特解
for(i = 1; i < d; i++)
printf("%d\n",(x0 + i * (n / d)) % n);
return true;
}

3 ) 求模乘法的逆:

模板:

1
2
3
4
5
6
7
//求得a在模n条件下的逆  
int inv(int a,int n)
{
int d,x,y;
gcd(a,n,d,x,y);
return d==1?(x+n)%n:-1;
}

求解模方程组(中国剩余定理):

模板(互质情况):

1
2
3
4
5
6
7
8
9
10
11
12
13
14
//中国剩余定理,求得M%A=a,M%B=b,...中的M,其中A,B,C...互质    
int china(int a[])
{
int i,j,k,d,ans=0,x,y,M=1;
for(i=0;i<n;i++)
M=M*x[i];
for(i=0;i<n;i++)
{
m=M/x[i];
gcd(x[i],m,d,x,y);
ans=(ans+y*m*a[i])%M;
}
return (ans+M)%M;
}

模板(不互质情况):

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
int China(int n)    
{
int m1,r1,m2,r2,flag=0,i,d,x,y,c,t;
scanf("%d%d",&m1,&r1);
flag=0;
for(i=1;i<n;i++)
{
scanf("%d%d",&m2,&r2);
if(flag)continue;
gcd(m1,m2,d,x,y);//d=gcd(m1,m2);x*m1+y*m2=d;
c=r2-r1;
if(c%d)//对于方程m1*x+m2*y=c,如果c不是d的倍数就无整数解
{
flag=1;
continue;
}
t=m2/d;//对于方程m1x+m2y=c=r2-r1,若(x0,y0)是一组整数解,那么(x0+k*m2/d,y0-k*m1/d)也是一组整数解(k为任意整数)
//其中x0=x*c/d,y0=x*c/d;
x=(c/d*x%t+t)%t;//保证x0是正数,因为x+k*t是解,(x%t+t)%t也必定是正数解(必定存在某个k使得(x%t+t)%t=x+k*t)
r1=m1*x+r1;//新求的r1就是前i组的解,Mi=m1*x+M(i-1)=r2-m2*y(m1为前i个m的最小公倍数);对m2取余时,余数为r2;
//对以前的m取余时,Mi%m=m1*x%m+M(i-1)%m=M(i-1)%m=r
m1=m1*m2/d;
}
if(flag)return -1;
else return r1;
}

2.3. 求解模方程a^x=b(mod n):

模版(n为素数。无解返回-1):

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
//求解模方程a^x=b(mod n)。n为素数,无解返回-1  
int log_mod (int a,int b,int n)
{
int m,v,e=1,i;
m=ceil(sqrt(n+0.5));
v=inv(pow_mod(a,m),n);
map<int,int>x;
x[1]=0;
for(i=1;i<m;i++)
{
e=(e*a)%n;
if(!x.count(e))x[e]=i;
}
for(i=0;i<m;i++)
{
if(x.count(b))return i*m+x[b];
b=(b*v)%n;
}
return -1;
}

模版(n不为素数):

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
//hdu 2815 Mod Tree  
#include <cstdio>
#include <cstring>
#include <cmath>
#include <map>
#include <iostream>
#include <algorithm>
using namespace std;
#define LL __int64
LL gcd(LL a,LL b)
{
return b==0?a:gcd(b,a%b);
}
//拓展欧几里得定理,求ax+by=gcd(a,b)的一组解(x,y),d=gcd(a,b)
void gcd_mod(LL a,LL b,LL &d,LL &x,LL &y)
{
if(!b){d=a;x=1;y=0;}
else{gcd_mod(b,a%b,d,y,x);y-=x*(a/b);}
}
//求解模方程d*a^(x-c)=b(mod n)。d,a和n互质,无解返回-1
LL log_mod (LL a,LL b,LL n,LL c,LL d)
{
LL m,v,e=1,i,x,y,dd;
m=ceil(sqrt(n+0.5)); //x=i*m+j
map<LL,LL>f;
f[1]=m;
for(i=1;i<m;i++) //建哈希表,保存a^0,a^1,...,a^m-1
{
e=(e*a)%n;
if(!f[e])f[e]=i;
}
e=(e*a)%n;//e=a^m
for(i=0;i<m;i++)//每次增加m次方,遍历所有1<=f<=n
{
gcd_mod(d,n,dd,x,y);//d*x+n*y=1-->(d*x)%n=1-->d*(x*b)%n==b
x=(x*b%n+n)%n;
if(f[x])
{
LL num=f[x];
f.clear();//需要清空,不然会爆内存
return c+i*m+(num==m?0:num);
}
d=(d*e)%n;
}

return -1;
}
int main()
{
LL a,b,n;
while(scanf("%I64d%I64d%I64d",&a,&n,&b)!=EOF)
{
if(b>=n)
{
printf("Orz,I can’t find D!\n");
continue;
}
if(b==0)
{
printf("0\n");
continue;
}
LL ans=0,c=0,d=1,t;
while((t=gcd(a,n))!=1)
{
if(b%t){ans=-1;break;}
c++;
n=n/t;
b=b/t;
d=d*a/t%n;
if(d==b){ans=c;break;}//特判下是否成立。
}
if(ans!=0)
{
if(ans==-1){printf("Orz,I can’t find D!\n");}
else printf("%I64d\n",ans);
}
else
{
ans=log_mod(a,b,n,c,d);
if(ans==-1)printf("Orz,I can’t find D!\n");
else printf("%I64d\n",ans);
}
}
return 0;
}
/*
求解模方程a^x=b(mod n),n不为素数。拓展Baby Step Giant Step
模板题。

方法:
初始d=1,c=0,i=0;
1.令g=gcd(a,n),若g==1则执行下一步。否则由于a^x=k*n+b;(k为某一整数),则(a/g)*a^k=k*(n/g)+b/g,(b/g为整除,若不成立则无解)
令n=n/g,d=d*a/g,b=b/g,c++则d*a^(x-c)=b(mod n),接着重复1步骤。
2.通过1步骤后,保证了a和d都与n互质,方程转换为d*a^(x-c)=b(mod n)。由于a和n互质,所以由欧拉定理a^phi(n)==1(mod n),(a,n互质)
可知,phi(n)<=n,a^0==1(mod n),所以构成循环,且循环节不大于n。从而推出如果存在解,则必定1<=x<n。(在此基础上我们就可以用
Baby Step Giant Step方法了)
3.令m=ceil(sqrt(n)),则m*m>=n。用哈希表存储a^0,a^1,..,a^(m-1),接着判断1~m*m-1中是否存在解。
4.为了减少时间,所以用哈希表缩减复杂度。分成m次遍历,每次遍历a^m长度。由于a和d都与n互质,所以gcd(d,n)=1,
所以用拓展的欧几里德定理求得d*x+n*y=gcd(d,n)=1,的一组整数解(x,y)。则d*x+n*y=1-->d*x%n=(d*x+n*y)%n=1-->d*(x*b)%n=((d*x)%n*b%n)%n=b。
所以若x*b在哈希表中存在,值为k,则a^k*d=b(mod n),答案就是ans=k+c+i*m。如果不存在,则令d=d*a^m,i++后遍历下一个a^m,直到遍历a^0到a^(m-1)还未找到,则说明不解并退出。
*/

欧拉phi函数:

求不超过n且与n互质的正整数个数.

模板:

1
2
3
4
5
6
7
8
9
10
11
12
13
int euler_phi(int n)//求单个欧拉函数  
{
int m=(int)sqrt(n+0.5);
int i,ans=n;
for(i=2;i<=m;i++)
if(n%i==0)
{
ans=ans/i*(i-1);
while(n%i==0)n/=i;
}
if(n>1)ans=ans/n*(n-1);
return ans;
}
int phi[maxn];    
void euler_phi()    
{    
    int i,j,k;    
    //欧拉函数,phi[i]表示不超过i的与i互质的整数个数    
    for(i=2;i<maxn;i++)phi[i]=0;    
    phi[1]=1;    
    for(i=2;i<maxn;i++)    
        if(!phi[i])    
        for(j=i;j<=maxn;j+=i){    
            if(!phi[j])phi[j]=j;    
            phi[j]=phi[j]/i*(i-1);    
        }    
}