一些数论概念与算法——从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}} 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$）

算法代码：

inline long long findRoot(long long p)
{    
    long long s = ceil(sqrt(p - 1));
    
    for(long long i = 1; i < p; ++i)
    {
        if(pow(i, p - 1, p) == 1long long)
        {
            long long sign = 1;
            for(long long j = 2; j <= s; ++j)
                if(((p - 1) % j == 0) && (pow(i, (p - 1) / j, p) == 1))
                {
                    sign = 0; break;
                }
            if(sign)return i;
        }
    }
}

快速幂：

上面寻找原根的程序中调用了一个函数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)$

算法代码：

long long pow(long long n, long long k, long long p)
{
    if(k == 0long long)return 1long long;
    long long w = 1;
    if(k & 1)w = n % p;
    long long ans = pow(n * n % p, k >> 1, p);
    return ans * w % p;
}

接下来的工作就是求离散对数啦：我们采用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

算法代码：

inline long long ind(long long x, long long p, long long g)
{
    static map<long long, long long> tmp;
    long long s = ceil(sqrt(p));
    long long w = pow(g, s, p);
    
    for(long long r = 0; r <= s; ++r)tmp.insert(make_pair(pow(g, r, p), r));
    for(long long t = 0; t <= s; ++t)
    {
        long long gt = pow(w, t, p);
        long long anti = findv(gt);
        if(tmp.count(x * anti % p))return (tmp[x * anti % p] + t * s);
    }
    return -1;
}

注意到在这个算法中，查找合条件的$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)$，按照欧几里德算法的方式递归下去就好了，这样我们解决了求解逆元的问题。而这道题的最后一步求解同余方程也要使用扩展欧几里德算法解决。

算法代码：

long long exgcd(long long a, long long b, long long &x1, long long &y1)
{
    if(b == 0){x1 = 1; y1 = 0; return a;}
    long long x0, y0;
    long long G = exgcd(b, a % b, x0, y0);
    x1 = y0; y1 = x0 - a / b * y0;
    return G;
}

求解逆元的代码如下：

inline long long findv(long long x)
{
    long long d, t;
    exgcd(x, p, d, t);
    while(d < 0)d += p;
    return d;
}

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

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

//date 20140323
#include <cstdio>
#include <cstring>
#include <cmath>
#include <map>
#include <algorithm>
#include <iostream>

using namespace std;

typedef long long LL;

LL x, k, a, p;

LL pow(LL n, LL k, LL p)
{
    if(k == 0LL)return 1LL;
    LL w = 1;
    if(k & 1)w = n % p;
    LL ans = pow(n * n % p, k >> 1, p);
    return ans * w % p;
}

inline LL findRoot(LL p)
{    
    LL s = ceil(sqrt(p - 1));
    
    for(LL i = 1; i < p; ++i)
    {
        if(pow(i, p - 1, p) == 1LL)
        {
            LL sign = 1;
            for(LL j = 2; j <= s; ++j)
                if(((p - 1) % j == 0) && (pow(i, (p - 1) / j, p) == 1))
                {
                    sign = 0; break;
                }
            if(sign)return i;
        }
    }
}

LL exgcd(LL a, LL b, LL &x1, LL &y1)
{
    if(b == 0){x1 = 1; y1 = 0; return a;}
    LL x0, y0;
    LL G = exgcd(b, a % b, x0, y0);
    x1 = y0; y1 = x0 - a / b * y0;
    return G;
}

inline LL findv(LL x)
{
    LL d, t;
    exgcd(x, p, d, t);
    while(d < 0)d += p;
    return d;
}

inline LL ind(LL x, LL p, LL g)
{
    static map<LL, LL> tmp;
    LL s = ceil(sqrt(p));
    LL w = pow(g, s, p);
    
    for(LL r = 0; r <= s; ++r)tmp.insert(make_pair(pow(g, r, p), r));
    for(LL t = 0; t <= s; ++t)
    {
        LL gt = pow(w, t, p);
        LL anti = findv(gt);
        if(tmp.count(x * anti % p))return (tmp[x * anti % p] + t * s);
    }
    return -1;
}


LL ans[1000000], rans[1000000];
int nans;

int main()
{
    freopen("sgu.in", "r", stdin);
    freopen("sgu.out", "w", stdout);
    
    cin >> p >> k >> a;

    if(a >= p)
    {
        printf("0
");
        return 0;
    }
    
    if(a == 0)
    {
        printf("1
0
");return 0;
    }
    LL g = findRoot(p);
    LL q1, q2 = ind(a, p, g), w;
    
    if(w == -1)
    {
        printf("0
");
        return 0;
    }
    
    LL G = exgcd(k, p - 1, q1, w);

    if(q2 % G != 0)
    {
        printf("0
");
        return 0;
    }
    
    
    q1 = q1 * (q2 / G) % (p - 1);
    LL delta = (p - 1) / G;
    
    nans = 0;
    for(int i = 0 ; i < G; ++i)
    {
        q1 = ((q1 + delta) % (p - 1) + p - 1) % (p - 1);
        ans[++nans] = pow(g, q1, p);        
    }
    
    sort(ans + 1, ans + nans + 1);
    
    int tot = 0;
    for(int i = 1; i <= nans; ++i)
    {
        if(ans[i] != ans[i - 1])rans[++tot] = ans[i];
    }
    cout << tot << endl;
    for(int i = 1; i < tot; ++i)cout << rans[i] << ' ';
    cout << rans[tot] << endl;
    return 0;
}

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