随笔- 26  评论- 12  文章- 0 

一些数论概念与算法——从SGU261谈起

话说好久没来博客上面写过东西了,之前集训过于辛苦了,但有很大的收获,我觉得有必要把它们拿出来总结分享。之前一直是个数论渣(小学初中没好好念过竞赛的缘故吧),经过一道题目对一些基础算法有了比较深刻的理解,在这里我打算系统地讲出这道题目涉及的大部分内容,希望可以帮到大家。

原题地址:http://acm.sgu.ru/problem.php?contest=0&problem=261

题目大意:给出质数$p$、$k$和一个自然数$a$,求关于$x$的同余方程$x^k \equiv a \pmod p$在区间$[0,\,p-1]$内所有解

数据范围:$2 \le p \le 10^9, 2 \le k \le 100000$

这种同余方程叫做“k次剩余”,这题也是一个模板题,但它涉及了比较多的数论算法,我会一一在这里尽可能讲清楚,这些概念和算法包括:

原根、缩系、离散对数、乘法逆元、线性同余方程、快速幂算法、扩展欧几里德算法、Baby Step - Gaint Step算法

虽然我希望近可能地系统化讲解,但是以下基本概念不会做任何讲解,不懂的同学可以查到大量的相关资料,它们包括:

素数、带余除法、模运算、幂运算、欧拉函数$\varphi(n)$、渐进复杂度分析

 

 首先我们通过引入一些概念对所求进行一些变形,最后再讲述求解所需的算法,前面的理论定义部分可能较为枯燥:

阶、原根、缩系与离散对数:

定义 1(阶):若$gcd(a,\,m)=1$,满足$a^r \equiv 1 \pmod m$的最小正整数$r$称为$a$模$m$的

性质 1(阶的判定定理):$gcd(a,\,m)=1$时,$a$模$m$的阶是$r$等价于以下两点:

  1:$a^r \equiv 1 \pmod m$

  2:对于$r$的每一个质因子$p$,均有$a^{\frac{r}{p}} \not\equiv 1 \pmod m$

  证明:第一条是阶的定义,第二条应用反证法,若存在质因子$p$使得其不成立,则与$r$的最小性矛盾

定义 2(原根):若整数$g$模$m$的阶为$\varphi(m)$,则称$g$为$m$的原根

性质 2(原根的存在情况):正整数$m$存在原根,当且仅当$m=2,\,4,\,p^a,\,2p^a$,其中$p$是奇素数且$a \ge 1$

由性质2我们可以得知所有的质数都有原根。

性质 3(原根的分布):一个数可能对应多个原根,而且最小的原根一般都很小(基本没例外)

定义 3(缩系):若模$m$有原根$g$,则$\{g^0,\,g^1,\,\dots,\,g^{\varphi(m)-1}\}$称为模$m$的缩系。显然对于每个与$m$互质的正整数$a$,在$m$的缩系中仅存在唯一的一个整数$k$,使得$g^r \equiv a \pmod m$

定义 4(离散对数):我们将形如$a^k \equiv n \pmod m$(m是质数)的式子中的$k$称为$n$以$a$为底在模$m$下的离散对数,特别地,在$a$为$m$的原根$g$时,我们称$k$为$n$模$m$的离散对数,由于$k$的值与原根$g$有关,我们将它记作$ind_{g}n$

由于我们讨论的模是质数,所以我们可以将不超过$m$的运算全部放到它的缩系里面来做了。而且不超过$m$的整数$n$模$m$的离散对数正巧与$m$的缩系中$n$对应的值相等。

性质 4(离散对数运算性质):

  类似于对数的运算性质,离散对数的性质大致有以下几条:

  1:若$gcd(a,\,m)=gcd(b,\,m)$,$a \equiv b \pmod m$等价于$ind_{g}a=ind_{g}b$

  2:$ind_{g}ab=ind_{g}a+ind_{g}b \pmod{\varphi(m)}$

  3:$ind_{g}a^n=n\,ind_{g}a \pmod{\varphi(m)}$

以上涉及的内容全在整数范围内

 

好了,进行了那么多定义和基础知识的灌输,题目涉及的理论基础也差不多说完了,让我们来回归题目本身吧:

$x^k \equiv a \pmod p$等价于$k\,ind_{g}x \equiv ind_{g}a \pmod{\varphi(p)}$,其中$g$是模$p$的原根(这个过程强烈建议读者手推)。这样,我们直接求解这个线性同余方程就好了,问题被分解为了求原根、离散对数以及求解线性同余方程,接下来我们开始讲解每步的求解过程,可以选择不会的地方跳着阅读

求解原根:

由上面的性质3可以知道最小的原根一般都很小,由于目前不存在寻找原根的比暴力枚举更有效的方法,所以我们直接枚举原根然后进行判定就好了

算法流程:

1.输入一个质数$p$

2.枚举1到$p-1$的所有整数$k$,判断$p-1$是否为$k$的阶(即$k^{p-1}$是否同余于1,并枚举$p-1$的所有质因子$r$判断$k^{\frac{p-1}{r}}$是否不同余于1,若有任何一条不满足,则$k$不是原根,如果两条都满足,则返回$k$)

算法代码:

 1 inline long long findRoot(long long p)
 2 {    
 3     long long s = ceil(sqrt(p - 1));
 4     
 5     for(long long i = 1; i < p; ++i)
 6     {
 7         if(pow(i, p - 1, p) == 1long long)
 8         {
 9             long long sign = 1;
10             for(long long j = 2; j <= s; ++j)
11                 if(((p - 1) % j == 0) && (pow(i, (p - 1) / j, p) == 1))
12                 {
13                     sign = 0; break;
14                 }
15             if(sign)return i;
16         }
17     }
18 }

快速幂:

上面寻找原根的程序中调用了一个函数pow(n, k, p),用来计算$n^k\mod p$的值。这个算法用到了指数运算的一些简单技巧

算法流程:

1.若k = 0, 返回 1

2.设w = pow(n*n mod p, k / 2, p)

3.若k为奇数,返回w * n mod p, 否则返回w

这个算法过于简单了,不多说

时间复杂度:$O(\lg k)$

算法代码:

1 long long pow(long long n, long long k, long long p)
2 {
3     if(k == 0long long)return 1long long;
4     long long w = 1;
5     if(k & 1)w = n % p;
6     long long ans = pow(n * n % p, k >> 1, p);
7     return ans * w % p;
8 }

接下来的工作就是求离散对数啦:我们采用Baby Step-Gaint Step算法(小步大步)

Baby Step - Gaint Step算法:

这一算法在$O(\sqrt{n})$的时间内解决求解离散对数的问题(关于$k$的方程$g^k \equiv n \pmod p$的解)

首先我们设$s=\sqrt{\lceil p \rceil}$,由于$k$的范围始终是$[0,\,p-1]$,所以可以设$k=st+r$

预先计算出所有小于等于$s$的$i$的$g^i$的值,然后从0到$s$枚举$t$,之后检查是否有满足条件的$r < s$存在,如果有,则$ts+r$则是符合条件的解。

注意寻找是否有满足条件的$r$的时候,我们有两种处理方式:第一种是使用哈希表或树形结构维护,第二种是将它们预处理之后排序再二分。这里为了方便我使用了标准库中的map

算法代码:

 1 inline long long ind(long long x, long long p, long long g)
 2 {
 3     static map<long long, long long> tmp;
 4     long long s = ceil(sqrt(p));
 5     long long w = pow(g, s, p);
 6     
 7     for(long long r = 0; r <= s; ++r)tmp.insert(make_pair(pow(g, r, p), r));
 8     for(long long t = 0; t <= s; ++t)
 9     {
10         long long gt = pow(w, t, p);
11         long long anti = findv(gt);
12         if(tmp.count(x * anti % p))return (tmp[x * anti % p] + t * s);
13     }
14     return -1;
15 }

注意到在这个算法中,查找合条件的$r$时我们使用了除法,但是在模意义下除法会出问题,我们该为乘以除数模意义下的乘法逆元。

乘法逆元:

定义 5(乘法逆元):若$xk \equiv 1 \pmod p$,则称$k$和$x$互为模$p$意义下的乘法逆元

解法:将同余方程$xk \equiv 1 \pmod p$转化为等式方程$xk-np=1\quad(n \in N)$,使用下面的扩展欧几里德算法求解

扩展欧几里德算法:

这个是比较基础但是讲起来比较麻烦的数论算法了,这里只讲推导过程和做法,严谨的证明请查阅相关书籍。

此算法用于求解形如$ax+by=gcd(a,\,b)$的不定方程的解$(x,\,y)$。根据欧几里德算法我们知道有$bx_0+(a \mod b)y_0=gcd(a,\,b)$

经过一些整理我们可以发现由于$a \mod b = a - b \times \lfloor \frac{a}{b} \rfloor$,$y_0a+\left(x_0+y_0 \times \lfloor \frac{a}{b} \rfloor \right)b=gcd(a,\,b)$

所以有$x=y_0, y = \left(x_0+y_0 \times \lfloor \frac{a}{b} \rfloor \right)$,按照欧几里德算法的方式递归下去就好了,这样我们解决了求解逆元的问题。而这道题的最后一步求解同余方程也要使用扩展欧几里德算法解决。

算法代码:

1 long long exgcd(long long a, long long b, long long &x1, long long &y1)
2 {
3     if(b == 0){x1 = 1; y1 = 0; return a;}
4     long long x0, y0;
5     long long G = exgcd(b, a % b, x0, y0);
6     x1 = y0; y1 = x0 - a / b * y0;
7     return G;
8 }

 

求解逆元的代码如下:

1 inline long long findv(long long x)
2 {
3     long long d, t;
4     exgcd(x, p, d, t);
5     while(d < 0)d += p;
6     return d;
7 }

 

(今晚太累了所以后面讲的可能不是很清楚,同学们有不明白的地方可以查阅相关资料或者给我留言,我会尽力完善)

到这里我们总结了一些常用的模运算下的算法,并完整地讲完了这道题,最后附上这道题的完整代码以供查错等需要

  1 //date 20140323
  2 #include <cstdio>
  3 #include <cstring>
  4 #include <cmath>
  5 #include <map>
  6 #include <algorithm>
  7 #include <iostream>
  8 
  9 using namespace std;
 10 
 11 typedef long long LL;
 12 
 13 LL x, k, a, p;
 14 
 15 LL pow(LL n, LL k, LL p)
 16 {
 17     if(k == 0LL)return 1LL;
 18     LL w = 1;
 19     if(k & 1)w = n % p;
 20     LL ans = pow(n * n % p, k >> 1, p);
 21     return ans * w % p;
 22 }
 23 
 24 inline LL findRoot(LL p)
 25 {    
 26     LL s = ceil(sqrt(p - 1));
 27     
 28     for(LL i = 1; i < p; ++i)
 29     {
 30         if(pow(i, p - 1, p) == 1LL)
 31         {
 32             LL sign = 1;
 33             for(LL j = 2; j <= s; ++j)
 34                 if(((p - 1) % j == 0) && (pow(i, (p - 1) / j, p) == 1))
 35                 {
 36                     sign = 0; break;
 37                 }
 38             if(sign)return i;
 39         }
 40     }
 41 }
 42 
 43 LL exgcd(LL a, LL b, LL &x1, LL &y1)
 44 {
 45     if(b == 0){x1 = 1; y1 = 0; return a;}
 46     LL x0, y0;
 47     LL G = exgcd(b, a % b, x0, y0);
 48     x1 = y0; y1 = x0 - a / b * y0;
 49     return G;
 50 }
 51 
 52 inline LL findv(LL x)
 53 {
 54     LL d, t;
 55     exgcd(x, p, d, t);
 56     while(d < 0)d += p;
 57     return d;
 58 }
 59 
 60 inline LL ind(LL x, LL p, LL g)
 61 {
 62     static map<LL, LL> tmp;
 63     LL s = ceil(sqrt(p));
 64     LL w = pow(g, s, p);
 65     
 66     for(LL r = 0; r <= s; ++r)tmp.insert(make_pair(pow(g, r, p), r));
 67     for(LL t = 0; t <= s; ++t)
 68     {
 69         LL gt = pow(w, t, p);
 70         LL anti = findv(gt);
 71         if(tmp.count(x * anti % p))return (tmp[x * anti % p] + t * s);
 72     }
 73     return -1;
 74 }
 75 
 76 
 77 LL ans[1000000], rans[1000000];
 78 int nans;
 79 
 80 int main()
 81 {
 82     freopen("sgu.in", "r", stdin);
 83     freopen("sgu.out", "w", stdout);
 84     
 85     cin >> p >> k >> a;
 86 
 87     if(a >= p)
 88     {
 89         printf("0\n");
 90         return 0;
 91     }
 92     
 93     if(a == 0)
 94     {
 95         printf("1\n0\n");return 0;
 96     }
 97     LL g = findRoot(p);
 98     LL q1, q2 = ind(a, p, g), w;
 99     
100     if(w == -1)
101     {
102         printf("0\n");
103         return 0;
104     }
105     
106     LL G = exgcd(k, p - 1, q1, w);
107 
108     if(q2 % G != 0)
109     {
110         printf("0\n");
111         return 0;
112     }
113     
114     
115     q1 = q1 * (q2 / G) % (p - 1);
116     LL delta = (p - 1) / G;
117     
118     nans = 0;
119     for(int i = 0 ; i < G; ++i)
120     {
121         q1 = ((q1 + delta) % (p - 1) + p - 1) % (p - 1);
122         ans[++nans] = pow(g, q1, p);        
123     }
124     
125     sort(ans + 1, ans + nans + 1);
126     
127     int tot = 0;
128     for(int i = 1; i <= nans; ++i)
129     {
130         if(ans[i] != ans[i - 1])rans[++tot] = ans[i];
131     }
132     cout << tot << endl;
133     for(int i = 1; i < tot; ++i)cout << rans[i] << ' ';
134     cout << rans[tot] << endl;
135     return 0;
136 }
SGU 261

结束语:数论在计算机科学中占据着重要的地位,这道简单的题也只能算是管中窥豹,后面我会学习新的知识并于大家分享

posted on 2014-03-24 21:57  SnowyJone  阅读(...)  评论(...编辑  收藏