数论专题（二）数论基础知识

转载自：https://blog.csdn.net/whereisherofrom/article/details/78922798

二、数论基础知识

       1、欧几里德算法(辗转相除法)

       2、扩展欧几里德定理

           a.线性同余

           b.同余方程求解

           c.逆元

       3、中国剩余定理（孙子定理）

       4、欧拉函数

           a.互素

           b.筛选法求解欧拉函数

           c.欧拉定理和费马小定理

       5、容斥原理

二、数论基础知识

      1、欧几里德定理（辗转相除法）

      定理：gcd(a, b) = gcd(b, a % b)。

      证明：a = kb + r = kb + a%b，则a % b = a - kb。令d为a和b的公约数，则d|a且d|b 根据整除的组合性原则，有d|(a-kb)，即d|(a%b)。

      这就说明如果d是a和b的公约数，那么d也一定是b和a%b的公约数，即两者的公约数是一样的，所以最大公约数也必定相等。

      这个定理可以直接用递归实现，代码如下：

    int gcd(int a, int b) {
        return b ? gcd(b, a%b) : a;
    }

      这个函数揭示了一个约定俗成的概念，即任何非零整数和零的最大公约数为它本身。

     【例题8】f[0] = 0, 当n>1时，f[n] = (f[n-1]+a) % b，给定a和b，问是否存在一个自然数k (0 <= k< b)，是f[n]永远都取不到的。

      永远有多远？并不是本题的范畴。

      但是可以发现的是这里的f[...]一定是有循环节的，如果在某个循环节内都无法找到那个自然数k，那么必定是永远都找不到了。

      求出f[n]的通项公式，为f[n] = an % b，令an = kb + r，那么这里的r = f[n]，如果t = gcd(a, b)，r = an-kb = t ( (a/t)n - (b/t)k )，则有t|r，要满足所有的r使得t|r,只有当t = 1的时候，于是这个问题的解也就出来了，只要求a和b的gcd，如果gcd(a, b) > 1，则存在一个k使得f[n]永远都取不到，直观的理解是当gcd(a, b) > 1，那么f[n]不可能是素数。

 

      2、扩展欧几里德定理

      a、线性同余

      线性同余方程（也可以叫模线性方程）是最基本的同余方程，即ax≡b (mod n)，其中a、b、n都为常量，x是未知数，这个方程可以进行一定的转化，得到：ax = kn + b,这里的k为任意整数，于是我们可以得到更加一般的形式即：ax + by + c = 0，这个方程就是二维空间中的直线方程，但是x和y的取值为整数，所以这个方程的解是一些排列成直线的点集。

      b、同余方程求解

      求解同余方程第一步是转化成一般式：ax + by = c，这个方程的求解步骤如下：

      i) 首先求出a和b的最大公约数d = gcd(a, b)，那么原方程可以转化成d(ax/d + by/d) = c，容易知道(ax/d + by/d)为整数，如若d不能整除b，方程必然无解，算法结束；否则进入ii)。

      ii) 由i)可以得知，方程有解则一定可以表示成 ax + by = c = gcd(a, b)*c'，那么我们先来看如何求解d = gcd(a, b) = ax + by，根据欧几里德定理，有：

      d = gcd(a, b) = gcd(b, a%b) = bx' + (a%b)y' = bx' + [a-b*(a/b)]y' = ay' + b[x' - (a/b)y']

      于是有x = y',  y = x' - (a/b)y'。

      由于gcd(a, b)是一个递归的计算，所以在求解(x, y)时，(x', y')其实已经利用递归计算出来了，递归出口为b == 0的时候（对比辗转相除，也是b == 0的时候递归结束），那么这时方程的解x0 = 1, y0 = 0。代码如下：

    #define LL __int64
    LL Extend_Euclid(LL a, LL b, LL &X, LL &Y) {
        LL q, temp;
                   if( !b ) {
            X = 1; Y = 0;
                            return a;
        }else {
            q = Extend_Euclid(b, a % b, X, Y);
            temp = X;
            X = Y;
            Y = temp - (a / b) * Y;
                             return q;
        }
    }

      扩展欧几里德算法和欧几里德算法的返回值一致，都是gcd(a, b)，传参多了两个未知数X, Y，采用引用的形式进行传递，对应上文提到的x, y，递归出口为b == 0，这时返回值为当前的a，因为gcd(a, 0) = a，(X, Y)初值为(1, 0)，然后经过回溯不断计算新的(X, Y)，这个计算是利用了之前的(X, Y)进行迭代计算的，直到回溯到最上层算法终止。最后得到的(X, Y)就是方程gcd(a, b) = ax + by的解。

      通过扩展欧几里德求的是ax + by = gcd(a, b)的解，令解为(x0, y0)，代入原方程，得：ax0 + by0 = gcd(a, b)，如果要求ax + by = c = gcd(a, b)*c'，可以将上式代入，得：ax + by = c = (ax0 + by0)c'，则x = x0c', y = y0c'，这里的(x, y)只是这个方程的其中一组解，x的通解为 { x0c' + kb/gcd(a, b) | k为任意整数 }，y的通解可以通过x通解的代入得出。

    【例题9】有两只青蛙，青蛙A和青蛙B，它们在一个首尾相接的数轴上。设青蛙A的出发点坐标是x，青蛙B的出发点坐标是y。青蛙A一次能跳m米，青蛙B一次能跳n米，两只青蛙跳一次所花费的时间相同。数轴总长L米。要求它们至少跳了几次以后才会碰面。 

      假设跳了t次后相遇，则可以列出方程：(x + mt) % L = (y + nt) % L

      将未知数t移到等式左边，常数移到等式右边，得到模线性方程：(m-n)t%L = (y-x)%L   （即 ax≡b (mod n) 的形式）

      利用扩展欧几里德定理可以求得t的通解{ t0 + kd | k为任意整数 }，由于这里需要求t的最小正整数，而t0不一定是最小的正整数，甚至有可能是负数，我们发现t的通解是关于d同余的，所以最后的解可以做如下处理：ans = (t0 % d + d) % d。

      c、逆元

      模逆元的最通俗含义可以效仿乘法，a*x = 1，则称x为a在乘法域上的逆（倒数）；同样，如果ax≡1 (mod n)，则称b为a模n的逆，简称逆元。求a模n的逆元，就是模线性方程ax≡b (mod n)中b等于1的特殊形式，可以用扩展欧几里德求解。并且在gcd(a, n) > 1时逆不存在。

      3、中国剩余定理

      上文提到了模线性方程的求解，再来介绍一种模线性方程组的求解，模线性方程组如图二-3-1所示，其中(ai, mi)都是已知量，求最小的x满足以下n个等式：

图二-3-1

      将模数保存在mod数组中，余数保存在rem数组中，则上面的问题可以表示成以下几个式子，我们的目的是要求出一个最小的正整数K满足所有等式：

      K = mod[0] * x[0] + rem[0]     (0)

      K = mod[1] * x[1] + rem[1]     (1)

      K = mod[2] * x[2] + rem[2]     (2)

      K = mod[3] * x[3] + rem[3]     (3)

      ... ...

     这里给出我的算法，大体的思想就是每次合并两个方程，经过n-1次合并后剩下一个方程，方程的自变量取0时得到最小正整数解。算法描述如下：

      i)   迭代器i = 0

      ii)  x[i] = (newMod[i]*k + newRem[i])            (k为任意整数)

      iii) 合并(i)和(i+1)，得 mod[i] * x[i] - mod[i+1] * x[i+1] = rem[i+1] - rem[i]

           将x[i]代入上式，有 newMod[i]*mod[i]*k - mod[i+1] * x[i+1] = rem[i+1] - rem[i] - newRem[i]*mod[i]

      iv)  那么产生了一个形如 a*k + b*x[i+1] = c的同余方程，

           其中a = newMod[i]*mod[i], b = - mod[i+1], c = rem[i+1] - rem[i] - newRem[i]*mod[i]

                       求解同余方程，如果a和b的gcd不能整除c，则整个同余方程组无解，算法结束；

                       否则，利用扩展欧几里德求解x[i+1]的通解，通解可以表示成 x[i+1] = (newMod[i+1]*k + newRem[i+1])            (k为任意整数)

      v)  迭代器i++,如果i == n算法结束，最后答案为 newRem[n-1] * mod[n-1] + rem[n-1]；否则跳转到ii)继续迭代计算。

      4、欧拉函数

      a、互素

      两个数a和b互素的定义为：gcd(a, b) = 1，那么如何求不大于n且与n互素的数的个数呢？

      朴素算法，枚举i从1到n，当gcd(i, n)=1时计数器++，算法时间复杂度O(n)。

      这里引入一个新的概念：用φ(n)表示不大于n且与n互素的数的个数，该函数以欧拉的名字命名，称为欧拉函数。

      如果n是一个素数，即n = p，那么φ(n) = p-1（所有小于n的都互素）；

      如果n是素数的k次幂，即n = p^k，那么φ(n) = p^k - p^(k-1) （除了p的倍数其它都互素）；

      如果m和n互素，那么φ(mn) = φ(m)φ(n)（可以利用上面两个性质进行推导）。

      将n分解成如图二-4-1的素因子形式，那么利用上面的定理可得φ(n)如图二-4-2所示：

图二-4-1

图二-4-2

      前面已经讲到n的因子分解复杂度为O(k)，所以欧拉函数的求解就是O(k)。

      b、筛选法求解欧拉函数

       由于欧拉函数的表示法和整数的素数拆分表示法很类似，都可以表示成一些素数的函数的乘积，所以同样可以利用筛选法进行求解。伪代码如下：

    #define MAXP 2000010
    #define LL __int64
    void Eratosthenes_Phi() {
        notprime[1] = true;
        for(int i = 1; i < MAXP; i++) phi[i] = 1;
        for(int i = 2; i < MAXP; i++) {
            if( !notprime[i] ) {
                phi[i] *= i - 1;
                // 和传统素数筛法的区别在于这个i+i
                for(int j = i+i; j < MAXP; j += i) {
                    notprime[j] = true;
                    int n = j / i;
                    phi[j] *= (i - 1);
                    while(n % i == 0) n /= i, phi[j] *= i;
                } 
            }
        }
    }

       这里的phi[i]保存了i这个数的欧拉函数，还是利用素数筛选将所有素数筛选出来，然后针对每个素因子计算它的倍数含有该素因子的个数，利用欧拉公式计算该素因子带来的欧拉函数分量，整个筛选过程可以参考素数筛选。

      c、欧拉定理和费马小定理

      欧拉定理：若n,a为正整数，且n,a互素，则: 。

      费马小定理：若p为素数，a为正整数且和p互素，则： 。

      由于当n为素数时φ(n) = p-1，可见费马小定理是欧拉定理的特殊形式。

      证明随处可见，这里讲一下应用。

     【例题10】整数a和n互素，求a的k次幂模n，其中k = X^Y, 正整数a,n,X,Y（X,Y<=10^9）为给定值。

      问题要求的是a^(X^Y) % n，指数上还是存在指数，需要将指数化简，注意到a和n互素，所以可以利用欧拉定理，令X^Y = kφ(n) + r，那么kφ(n)部分并不需要考虑，问题转化成求r = X^Y % φ(n)，可以采用快速幂取模，二分求解，得到r后再采用快速幂取模求解a^r % n。

      5、容斥原理

      容斥原理是应用在集合上的，来看图二-5-1，要求图中两个圆的并面积，我们的做法是先将两个圆的面积相加，然后发现相交的部分多加了一次，予以减去；对于图二-5-2的三个圆的并面积，则是先将三个圆的面积相加，然后减去两两相交的部分，而三个圆相交的部分被多减了一次，予以加回。

图二-5-1

图二-5-2

     这里的“加”就是“容”，“减”就是“斥”，并且“容”和“斥”总是交替进行的（一个的加上，两个的减去，三个的加上，四个的减去），而且可以推广到n个元素的情况。

    【例题11】求小于等于m(m < 2^31)并且与n(n < 2^31)互素的数的个数。

     当m等于n，就是一个简单的欧拉函数求解。

     但是一般情况m都是不等于n的，所以可以直接摈弃欧拉函数的思路了。

     考虑将n分解成素数幂的乘积，来看一种最简单的情况，当n为素数的幂即n = p^k时，显然答案等于m - m/p（m/p表示的是p的倍数，去掉p的倍数，则都是和n互素的数了）；然后再来讨论n是两个素数的幂的乘积的情况，即n = p1^k1 * p2^k2，那么我们需要做的就是找到p1的倍数和p2的倍数，并且要减去p1和p2的公公倍数，这个思想其实已经是容斥了，所以这种情况下答案为：m - ( m/p1 + m/p2 - m/(p1*p2) )。

      类比两个素因子，如果n分解成s个素因子，也同样可以用容斥原理求解。

      容斥原理其实是枚举子集的过程，常见的枚举方法为dfs，也可以采用二进制法（0表示取，1表示不取）。这里给出一版dfs版本的容斥原理的伪代码，用于求解小于等于m且与n互素的数的个数。

    #define LL __int64
    void IncludeExclude(int depth, LL m, LL mul, int op, int* p, LL &ans) {
        if(m < mul) return ; 
        if(depth == p[0]) { 
            ans += (op ? -1 : 1) * (m / mul);
            return ;
        }
        for(int i = 0; i < 2; i++) { 
            // 0 表示不取, 1表示取
            IncludeExclude( depth+1, m, mul * (i?p[depth+1]:1), op^i, p, ans );
        }
    }

      p[ 1 : p[0] ]存储的是n的所有素因子，p[0]表示数组长度，mul表示该次的素因子子集的乘积，op表示子集的奇偶性，ans存储最后的答案。

      例如求[1, 9]中和6互素的数的个数，这时p = [2, 2, 3] （注意p[0]是存素数的个数的，6分解的素因子为2和3）。

 

      ans = 9/1 - （9/2 + 9/3） + 9/6 = 3，ans分为三部分，0个数的组合，1个数的组合，2个数的组合。