初等数论入门

简介

数论是纯粹数学的分支之一，主要研究整数的性质。整数可以是方程式的解。有些解析函数中包括了一些整数、质数的性质，透过这些函数也可以了解一些数论的问题。

数论被高斯成为“数学皇后”，是一大必不可少的分支。

初等数论主要研究的是整除理论和同余理论。下面将会给出本博客介绍的东西：

• 逆元。
• 欧几里得和扩展欧几里得。
• 中国剩余定理和扩展中国剩余定理。
• BSGS（大步小步算法）算法。
• Lucas 定理。

模意义与逆元

考虑这样一个整数除法 $27\div 4 = 6\cdots 3$。这样的式子可以用模运算记为 $27\bmod 4 = 3$ 或 $27\equiv 3\pmod 4$。它们的意义有略微的区别：前者表示 $27$ 除以 $4$ 的余数是 $3$，后者是 $27$ 在模 $4$ 意义下和 $3$ 是相同的。

对于三个数 $a,b,p$，若 $a\div p,b\div p$ 余数相同，就称 $a\equiv b\pmod p$。

模运算有如下性质：

$(a+b)\bmod p = (a\bmod p + b\bmod p)\bmod p$

$(a-b)\bmod p = (a\bmod p - b\bmod p)\bmod p$

$ab\bmod p = [(a\bmod p) (b\bmod p)]\bmod p$

注意 $\frac{a}{b}\mod p \neq \frac{a\bmod p}{b\mod p}\bmod p$。例如 $24\bmod 13 = 9, 12\bmod 13 =12 , 2\mod 13 =2$。

如果我们得到的式子需要用除法该怎么办？不妨把除法转化成乘法。这个时候就相当于，我们需要找到一个数在模意义下的“倒数”，也就是找到一个数 $k$，满足 $a\div b\equiv ak\pmod p$。这个 $k$ 就叫做 $b$ 在模 $p$ 意义下的逆元。

如何求逆元？当 $p$ 是质数是，可以用费马小定理：

$a^{p-1} \equiv 1\pmod p$

于是有 $a^{p-2}$ 和 $a$ 在模 $p$ 意义下互为倒数。因此 $a^{p-2}$ 是 $a$ 的一个逆元。

若 $p$ 不是质数，后面还会提到用欧几里得求逆元。

扩展欧几里得

数论题中经常会要求你求出一个这样的不定方程的整数解（其中参数也是整数）：

$ax+by=c$

这时如果直接暴力枚举显然不现实。因此我们需要一个优秀的算法解决这类问题。

裴蜀定理

对于不定方程 $ax+by=c$，$x,y$ 有整数解当且仅当 $\gcd(a,b)\mid c$，即 $\gcd(a,b)$ 是 $c$ 的约数。

欧几里得算法即辗转相除法，用于计算 $\gcd(a,b)$。需要用到一个公式：

$\gcd(a,b) = \gcd(b,a\bmod b)$

当 $b=0$ 时，$\gcd(a,b) = a$。因此直接递归就可以求解：

inline int gcd(int a,int b){
    if(!b) return a;
    return gcd(b,a%b);
}

当然也有非递归版本：

inline int gcd(int x,int y){
    while(y^=x^=y^=x%=y);
    return x;
}

这个看起来比较阴间，其实就是每次 $x\bmod y$ 后交换 $x,y$，直到 $y$ 为 $0$。

这时我们已经可以对这个方程有一个无解与否的判断了：如果原方程没有整数解，则 $\gcd(a,b)\nmid c$。然而我们还没有解这个方程。这就需要用到扩展欧几里得了。

扩展欧几里得是一个非常强大的算法，本质是边算 $\gcd(a,b)$ 边求 $ax+by=c$ 的解。

扩展欧几里得算法的基本思想，就是在通过辗转相除法求得 $\gcd(a,b)$ 的同时，通过辗转相处过程中的一些式子，倒回去推导出 $x,y$。

假设我们已经找到了一组特殊解 $(x_0,y_0)$ 以及 $k=\gcd(a,b)$。则可以用这一组特解表示所有的解：

$\begin{cases} x = x_0 + \frac{b}{k}t \\ y = y_0 - \frac{a}{k}t \end{cases}$

考虑如何求出一组特解 $(x_0,y_0)$，考虑辗转相除法的结束状态 $a'=\gcd,b'=0$。此时令 $x=1$ 有 $a'+yb'=c'$，因为 $b'=0$，所以 $y$ 是多少都没关系，当 $x=1$ 时，$c'= a' = \gcd(a,b)$，此时一定有解。

考虑推回去，对于一个状态 $(a,b)$，我们已经求出了 $(b,a\bmod b)$ 的解，也就是一组 $x_0,y_0$ 满足 $bx_0 + (a\bmod b)y_0 = \gcd$。

而 $a\bmod b = a-\left\lfloor\frac{a}{b}\right\rfloor b$：

$\begin{aligned} \gcd &= bx_0 + (a-\left\lfloor\frac{a}{b}\right\rfloor b)y_0 \\ &= bx_0 + ay_0 - \left\lfloor\frac{a}{b}\right\rfloor b y_0 \end{aligned}$

对于当前状态 $a,b$，要让 $ax+by=\gcd$，其实就是：

$x\gets y_0$

$y\gets x_0 - \left\lfloor\frac{a}{b}\right\rfloor y_0$

于是就得到一个状态怎么用另一个状态得到了，有代码：

int x1,x2,res;
inline int exgcd(int a,int b){ // 扩展欧几里得
    if(!b){
        x1=1,x2=0;
        return a;
    }
    res=exgcd(b,a%b);
    int x3=x1;
    x1=x2,x2=x3-x2*(a/b);
    return res;
}

还有一种传址实现：

inline int exgcd(int a,int b,int &x,int &y){
    if(!b) return x=1,y=0,a;
    int res=exgcd(b,a%b,y,x);
    y-=a/b*x;
    return res;
}

给定 $a,p$，求 $a$ 模 $p$ 意义下的逆元。

注意到若 $x$ 为 $a$ 的逆元，有 $ax\equiv 1\pmod p$，于是可以转化成求 $ax+py=1$ 的整数解。就可以用扩欧了。

扩展中国剩余定理

来趁热打铁地了解一下扩展中国剩余定理。不过在此之前还要看一下什么叫做中国剩余定理（CRT）：

孙子定理是中国古代求解一次同余式组（见同余）的方法。是数论中一个重要定理。又称中国余数定理。一元线性同余方程组问题最早可见于中国南北朝时期（公元5世纪）的数学著作《孙子算经》卷下第二十六题，叫做“物不知数”问题。

中国剩余定理需要解决这样一类问题：求解下面这一个一元线性同余方程组的整数解：

$\begin{cases} x \equiv a_1 \pmod{p_1} \\ x \equiv a_2 \pmod{p_2} \\ \quad\vdots \\ x \equiv a_n \pmod{p_n} \end{cases}$

其中 $p_1,p_2\cdots,p_n$ 两两互质。

首先需要知道的是，若 $a\equiv b \pmod p$，则 $a\div k \equiv b\div k \pmod p,a+pk\equiv b\pmod p$。

下面举一个例子，这是一个小学奥数题：计算一个数，它除以 $3$ 余 $2$，除以 $5$ 余 $7$，除以 $7$ 余 $2$。

先给每个方程依次写出一个解 $x_1,x_2,x_3$。考虑从 $x_1$ 除法，找一个 $x_2$ 使得 $(x_1+x_2)\bmod 3 =2$。根据前面的内容得到满足条件需要 $3\mid x_2$。同理，再找一个数 $x_3$ 满足 $(x_1+x_2+x_3)\bmod 3 =2$ 需要 $3\mid x_3$。同理总共能写出如下条件：

$5\mid x_1, 7\mid x_1,3\mid x_2,7 \mid x_2, 3\mid x_3, 5\mid x_3$

如果找到了符合条件的 $x_1,x_2,x_3$，把它们加起来就是原方程的一个解。

考虑 $x_1$ 有关的条件，相当于在 $5\times7=35$ 的倍数里找 $x_1$ 满足 $x_1\bmod 3 = 2$。为 $2$ 不好求，考虑求一个余数为 $1$，然后让这个解乘上 $2$。而求一个 $x_1$ 满足 $x_1\equiv 1\pmod 3$ 的形式和上面的扩展欧几里得求逆元是一样的，简单转化即可。

对于其他情况也用类似方法分析即可。

然而如果 $p_1,p_2,\cdots,p_n$ 不保证两两互质，上面的做法就有问题了，不信可以手玩一下。这时候就要用扩展中国剩余定理（exCRT）。

这时假设单独提出前两个方程：

$\begin{cases} x \equiv a_1 \pmod{p_1} \\ x \equiv a_2 \pmod{p_2} \end{cases}$

可以简单转化为：

$\begin{cases} x = a_1 + c_1 p_1 \\ x = a_2 + c_2 p_2 \end{cases}$

因此 $p_1 c_1 -p_2 c_2 = a_2 - a_1$。这是一个不定方程的形式，可以用扩欧求出一组 $(c_1,c_2)$。将 $c_1$ 代入求 $x'=p_1 c_1 + a_1$，这个 $x'$ 对于前两个方程成立。为了让它对后面的方程成立，探求其和后面解的关系，把第三个方程联立进来，求解新的方程的解 $x$。合并前两个方成后联立：

$\begin{cases} x \equiv x' \pmod{\rm{lcm}(p_1,p_2)} \\ x \equiv a_3 \pmod{p_3} \end{cases}$

然后就又求出了新的 $x$。再接着联立第四个方程……依次类推即可。

int x1,x2,res;
inline int exgcd(int a,int b){ // 扩展欧几里得
    if(!b){
        x1=1,x2=0;
        return a;
    }
    res=exgcd(b,a%b);
    int x3=x1;
    x1=x2,x2=x3-x2*(a/b);
    return res;
}
inline int china(){ // 扩展中国剩余定理
    int a1=a[1],b1=b[1];
    for(register int i=2;i<=n;++i){
        int a2=a[i],b2=b[i],c=b2-b1;
        exgcd(a1,a2);

        int mod=a2/res;
        x1=((x1*c/res)%mod+mod)%mod;

        int LCM=lcm(a1,a2);
        b1=((a1*x1+b1)%LCM+LCM)%LCM;
        a1=LCM;
    }
    return b1%a1;
}

BSGS

考虑这样一个问题，给定 $a,b,p$，$p$ 是质数，求下面方程的最小非负整数解：

$a^x\equiv b\pmod p$

求法并不复杂，将 $x$ 拆分为 $im-j$ 的形式，其中 $m=\left\lceil\sqrt p\right\rceil$，也就是求出：

$a^{im-j}\equiv b\pmod p$

也就是：

$a^{im}\equiv ba^j \pmod p$

开一个 Hash，从 $0\sim m$ 枚举 $j$，把 $ba^j \bmod p$ 丢进去。然后 $1\sim m$ 枚举 $i$，检查 Hash 表里面 $a^{im}\bmod p$ 是否存在，如果存在则取出其对应的 $j$。显然取出的第一个解就是原方程的最小非负整数解。

如果每个数的值域很小，快速幂能预处理，一次询问的复杂度就是 $\sqrt p$ 的。否则是 $\sqrt p\log n$ 的（要算快速幂）。如果用 STL 实现就是 $\sqrt p \log (n\sqrt p)$ 的，

// mp 是哈希数组
// mi 是快速幂，都会打吧，没必要放出来
inline int bsgs(int x,int z,int mod){ 
// 求 x 的 k 次方模 mod 等于 z 的最小的 k
    mp.clear();
    x%=mod,z%=mod;
    int p=sqrt(1.0*mod)+1;

    if(!z && !x) return 1;
    if(!x && z) return -1;
    // 注意特殊情况的判断
    for(register int i=0;i<=p;++i)
        mp[z*mi(x,i,mod)%mod]=i;
    x=mi(x,p,mod);
    for(register int i=0;i<=p;++i){
        int a=mi(x,i,mod);
        if(mp[a]){
            int b=mp[a];
            if(i*p>=b) return i*p-b;
        }
    }
    return -1;
}

Lucas 定理

Lucas 定理用中文表述就是卢卡斯定理。一般用于快速求一个在模意义下的组合数。

虽然把它放在最后一个（写这篇文章时的最后一个，如果后面还有别的知识都是本篇文章的更新版本加的了），但是它可能比上面讲的很多知识都简单许多。确切地说，卢卡斯定理就是：

设 $p$ 为质数，有（这里为了简化式子就将斜杠当做整除了）：

$\binom{n}{m} \equiv \binom{n\div p}{m\div p}\binom{n\bmod p}{m\bmod p} \pmod p$

其中有除法的那一个，可以递归实现，另一个直接套组合数公式求即可（很简单吧）：

/// fac[i] 表示 i 的阶乘
// inv[i] 表示 fac[i] 的逆元
inline int C(int n,int m,int p){
    if(m>n) return 0;
    return (fac[n]*inv[m]%p*inv[n-m]%p)%p;
}
inline int lucas(int n,int m,int p){
    if(!m) return 1;
    return lucas(n/p,m/p,p)*C(n%p,m%p,p)%p;
}