话说好久没来博客上面写过东西了，之前集训过于辛苦了，但有很大的收获，我觉得有必要把它们拿出来总结分享。之前一直是个数论渣（小学初中没好好念过竞赛的缘故吧），经过一道题目对一些基础算法有了比较深刻的理解，在这里我打算系统地讲出这道题目涉及的大部分内容，希望可以帮到大家。
原题地址：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}} otequiv 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^kmod 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 ceil}$，由于$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=1quad(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 imes lfloor frac{a}{b} floor$，$y_0a+left(x_0+y_0 imes lfloor frac{a}{b} floor ight)b=gcd(a,\,b)$
所以有$x=y_0, y = left(x_0+y_0 imes lfloor frac{a}{b} floor ight)$，按照欧几里德算法的方式递归下去就好了，这样我们解决了求解逆元的问题。而这道题的最后一步求解同余方程也要使用扩展欧几里德算法解决。
算法代码：
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 "); 90 return 0; 91 } 92 93 if(a == 0) 94 { 95 printf("1 0 ");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 "); 103 return 0; 104 } 105 106 LL G = exgcd(k, p - 1, q1, w); 107 108 if(q2 % G != 0) 109 { 110 printf("0 "); 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 }
结束语：数论在计算机科学中占据着重要的地位，这道简单的题也只能算是管中窥豹，后面我会学习新的知识并于大家分享