[数论第二节]欧拉函数/快速幂/扩展欧几里得算法

• 欧拉函数

  • 欧拉函数$\phi (N)$ : 1-N中与N互质的数的个数
  • 若$N=p_1^{a_1}·p_2^{a_2}·p_3^{a_3}····p_n^{a_n}$ 其中p为N的所有质因子
  • 则$\phi (N)=N(1-\frac{1}{p_1})(1-\frac{1}{p_2})···(1-\frac{1}{p_n})$

  • 证明：

    • 互质：两数的公共因子只有1
    • 去掉所有与N有(大于1的)公共因子的数，剩下的数就是与N互质的数
    • 对N的所有质因子$p_k$，去掉所有$\underset{―}{\mathrm{质}\mathrm{数}{p}_k\mathrm{的}\mathrm{倍}\mathrm{数}}$(与N有公共因子的数)，$\underset{―}{\mathrm{每}\mathrm{个}\mathrm{质}\mathrm{数}\mathrm{的}\mathrm{倍}\mathrm{数}}$个数为 $\frac{N}{p_k}$，即$$N-\frac{N}{p_1}-\frac{N}{p_2}-···-\frac{N}{p_n}$$
    • 对于数 $p_i·p_j$，在去掉 $p_i$ 的倍数和去掉 $p_j$ 的倍数的过程中去除了两次，所以要加上一次，即$$N-(\frac{N}{p_1}+\frac{N}{p_2}+···+\frac{N}{p_n})+(\frac{N}{p_1{p}_2}+\frac{N}{p_1{p}_3}+···+\frac{N}{p_{n-1}{p}_n})$$
    • 对于数 $p_i·p_j·p_k$，在去掉$p_k$的倍数的过程中被去掉了三次，在加上合数 $p_i·p_j$ 的过程中被加上了三次，所以合数 $p_i·p_j·p_k$ 没有被去掉，因此要去掉它，即$$N-(\frac{N}{p_1}-\frac{N}{p_2}-···-\frac{N}{p_n})+(\frac{N}{p_1{p}_2}+\frac{N}{p_1{p}_3}+···+\frac{N}{p_{n-1}{p}_n})-(\frac{N}{p_1{p}_2{p}_3}+\frac{N}{p_1{p}_2{p}_4}+···+\frac{N}{p_{n-2}{p}_{n-1}{p}_n})$$
    • 对于合数 $p_i·p_j·p_k·p_m$ 同理，归纳递推可知，所有质数个数为：$$N-{\displaystyle \sum _{i=1}^n\frac{N}{p_i}+{\displaystyle \sum _{1<=i<j<=n}\frac{N}{p_i{p}_j}-{\displaystyle \sum _{1<=i<j<k<=n}\frac{N}{p_i{p}_j{p}_k}+···+(-1{)}^{2n-1}{\displaystyle \sum _{1<=i<j<k<···<=n}\frac{N}{p_i{p}_j{p}_k···p}}}}}$$
    • 同样基于容斥原理：

      • $$|A\cup B\cup C|=|A|+|B|+|C|-|A\cap B|-|A\cap C|-|B\cap C|+|A\cap B\cap C|$$

      • $$|{\displaystyle {\cup }_{i=1}^n{A}_i|=\sum _i|A_i|-\sum _{i,j}|A_i\cap A_j|+\dots +(-1{)}^{n+1}|{\cap }_{i=1}^n{A}_i|}$$

      • 与N有(大于1的)公共因子的个数 =| $p_i$ 的倍数的个数| - |$p_i{p}_j$ 的倍数的个数| + |$p_i{p}_j{p}_k$ 的倍数的个数 ··· $(-1{)}^{n+1}$ |$p_i{p}_j···p_n$ 的倍数的个数|
    • 将上式因式分解后：$$N(1-\frac{1}{p_1})(1-\frac{1}{p_2})···(1-\frac{1}{p_n})$$

  • 容斥原理：

    int res = a;
    for(int i = 2; i <= a / i; ++ i){
        if(a % i == 0){//i为质因子
            res = res / i * (i - 1);//套公式
            while(a % i == 0) a /= i;//把因子除干净
        }
    }
    if(a > 1) res = res / a * (a - 1);//最后一个因子可能大于sqrt(a)
    

  • 筛选法：

    • 利用线性筛选质数的过程求出每个数的欧拉函数
    • 欧拉函数为积性函数，当a与b互质时有$$\phi (a·b)=\phi (a)·\phi (b)$$
    • 当i为质数时，$$\phi (i)=i-1$$
    • 当$\frac{i}{p_j}=0$时，$p_j$为i的质因子，此时$p_j·i$的质因子与i的质因子完全相同，所以 $$\phi (pj·i)=pj·i·(1-\frac{1}{p_1})(1-\frac{1}{p_2})···(1-\frac{1}{p_n})=p_j·\phi (i)$$
    • 当$\frac{i}{p_j}\ne 0$时，$p_j$不为i的质因子，此时$p_j·i$的质因子比i的质因子多一个$p_j$，所以 $$\phi (pj·i)=pj·i·(1-\frac{1}{p_1})(1-\frac{1}{p_2})···(1-\frac{1}{p_n})(1-\frac{1}{p_j})=p_j·\phi (i)·(1-\frac{1}{p_j})=(p_j-1)·\phi (i)$$
    • 代码：

      const int N = 1e6 + 10;
      typedef long long LL;
      int primes[N], cnt;//质数数组，下标
      int st[N];//标记为合数
      int phi[N];//欧拉函数
      int n;
      
      LL get_eulers(int n){
          phi[1] = 1;
          for(int i = 2; i <= n; ++ i){
              if(!st[i]){
                  primes[ ++ cnt] = i;
                  phi[i] = i - 1;//质数i的欧拉函数为i - 1
              }
              for(int j = 1; primes[j] <= n / i; ++ j){
                  int pj = primes[j];
                  st[pj * i] = 1;
                  if(i % pj == 0){//i是合数
                      phi[pj * i] = pj * phi[i];//pj为i的质因子
                      break;
                  }
                  else phi[pj * i] = (pj - 1) * phi[i];//pj不是i的因子
              }
          }
          
          LL ans = 0;
          for(int i = 1; i <= n; ++ i) ans += phi[i];
          return ans;
      }
      

  • 欧拉函数的应用

    • 欧拉定理

      • 若a与n互质，则$$a^{\phi (n)}(\mathrm{mod}n)\equiv 1$$
      • 证明：
        • 若a与b互质，且$a(\mathrm{mod}b)\ne 0$则$a(\mathrm{mod}b)$ 也与b互质
        • 设$a_1,a_2,a_3···a_{\phi (n)}$ 为1~n中的所有与n互质的数
        • 每个数同时乘上a得$a·a_1,a·a_2,a·a_3···a·a_{\phi (n)}$ 由于a也与n互质，所以乘a后所得的数也全部与n互质
        • 将每个数mod n，取模后的每个数都在1~n范围内，并且仍然都与n互质，显然取模后的这几个数就是原来的几个数，只不过数的顺序可能发生了变化
        • 将所有数相乘，则有$$a·a_1·a·a_2·a·a_3···a·a_{\phi (n)}(\mathrm{mod}n)=a_1·a_2·a_3···a_{\phi (n)}$$
        • 所以有：$$a^{\phi (n)}·(a_1{a}_2{a}_3···a_{\phi (n)})(\mathrm{mod}n)\equiv (a_1·a_2·a_3···a_{\phi (n)})$$即$$a^{\phi (n)}(\mathrm{mod}n)\equiv 1$$
        • 特别地，当n为质数时，$a^{n-1}(\mathrm{mod}n)\equiv 1$ 为费马定理

• 快速幂

  • 快速地求出$a^k(\mathrm{mod}b)$，时间复杂度为$O(log(k))$
  • 将k拆分为二进制相加，即

  • $$k=c_1·2^1+c_2·2^2+c_3·2^3+···+c_n·2^n$$

  • 其中 $c_i$ 是k的二进制的每一位(0/1)，n最大为k的二进制的位数
  • 所以

  • $$a^k(\mathrm{mod}b)=a^{c_1·2^1+c_2·2^2+c_3·2^3+···+c_n·2^n}(\mathrm{mod}b)=a^{c_1·2^1}·a^{c_2·2^2}·a^{c_3·2^3}···a^{c_n·2^n}(\mathrm{mod}b)$$

  • 所以每次遍历k的所有二进制位，同时预处理出$a^i$，根据k的二进制的取值将答案累乘
  • 代码：

    typedef long long LL;
    //求a^k mod b
    int qmi(int a, int k, int b){
    	int res = 1;
    	while(k){
    		if(k & 1) res = (LL)res * a % b;//k的当前位非0，则将a累乘到答案
    		k >>= 1;//k右移一位
    		a = (LL)a * a % b;//a的幂倍增
    	}
    	return res;
    }
    

  • 快速幂求逆元

    • 若b与p互质，且 $b|a$，$\frac{a}{b}=a·x(\mathrm{mod}p)$ ，则x为b的逆元
    • 两边乘b有$$a=a·b·x(\mathrm{mod}p)$$
    • 所以有$$b·x=1(\mathrm{mod}p)$$
    • 若p是质数，有费马定理$$b^{p-1}=1(\mathrm{mod}p)$$则$$x=b^{p-2}$$
    • 若p不是质数，有欧拉定理$$b^{\phi (p)}=1(\mathrm{mod}p)$$则$$x=b^{\phi (p)-1}$$
    • 代码：

      int qmi(int a, int b, int p){
      	略...
      }
      if(b与q互质){ //互质是有解的前提
      	if(p为质数) cout << qmi(b, p - 2, p);
      	else cout << qmi(b, phi[p] - 1, p);
      }else 无解
      

• 扩展欧几里得定理

  • 裴蜀定理：对于任意正整数a,b, 都一定存在整数x,y, 使得$ax+by=gcd(a,b)$, 并且 $gcd(a,b)$ 一定是a与b能构造出来的最小公约数

  • 扩展欧几里得算法就是求这样的x，y

  • 用于求解方程 $ax+by=gcd(a,b)$ 的解

  • 当 $b=0$ 时 $ax+by=a$ 故而 $x=1,y=0$

  • 当 $b\ne 0$ 时,因为

  • $$gcd(a,b)=gcd(b,a\mathrm{\%}b)$$

  • 而

  • $$b·x_1+(a\mathrm{\%}b)·y_1=gcd(b,a\mathrm{\%}b)$$

  • $$b·x_1+(a-\lfloor a/b\rfloor ·b)·y_1=gcd(b,a\mathrm{\%}b)$$

  • $$a{y}_1+b·(x_1-\lfloor a/b\rfloor )=gcd(b,a\mathrm{\%}b)=gcd(a,b)$$

  • 故而 $$x=y_1,y=x_1-\lfloor a/b\rfloor ·y_1$$

  • 代码：

    int exgcd(int a, int b, int &x, int &y){
        if(!b){ //b = 0时，ax + by = gcd(a, 0) = a,所以x = 1, y = 0;
            x = 1, y = 0;
            return a;
        }else {
            int d = exgcd(b, a % b, x, y);//交换a与b，x与y
            int t = x;
            x = y;//x1 = y2
            y = t - a / b * y;//y1 = x2 - a / b * y2
            return d;
        }
    }
    
    //简化版
    void exgcd(int a, int b, int &x, int &y){
        if(!b){
            x = 1, y = 0;
            return a;
        }else{
    	    int d = exgcd(b, a % b, y, x);
    	    y -= a / b * x;
    	    return d;
        }    
    }
    

  • 扩展欧几里得求解线性同余方程

    • 线性同余方程：$$ax\equiv b(\mathrm{mod}m)$$ 其中a，x，b，m均为整数，x为未知数，方程等价于 $ax=km+b$，也等价于 $$ax+mk=b$$其中x,k均为未知数
    • 方程形如：$$ax+by=c$$ 的直线方程，解(x,y)为直线上散列的点
    • 求出a，b的最大公约数$d=gcd(a,b)$，方程化为 $$d(x\frac{a}{d}+y\frac{b}{d})=c$$易知其中 $x\frac{a}{d}+y\frac{b}{d}$ 为整数，所以要求 $d|c$，故 $$c=k·d=k·gcd(a,b)$$
    • 所以线性同余方程$ax\equiv b(\mathrm{mod}m)$有解的充要条件为b为$gcd(a,m)$的倍数
    • 当 $b\ne k^{\prime }·gcd(a,m)$ 时($k^{\prime }$ 为整数)，方程无解
    • 当 $b=k^{\prime }·gcd(a,m)$ 时，先用扩展欧几里得求出方程 $ax+mk=gcd(a,m)$ 的特解$(x,k)$，可知 $ax+mk=b$ 方程的解为 $\frac{b}{gcd(a,m)}·(x,k)$ ，所以原式的解为 $\frac{b}{gcd(a,m)}·x$
    • 代码：

      const int N = 1e5 + 10;
      typedef long long LL;
      int n; 
      int a, b, m;
      int x, y;
      //扩展欧几里得算法
      int exgcd(int a, int b , int &x, int &y){
          if(!b){
              x = 1, y = 0;
              return a;
          }else{
              int t = exgcd(b, a % b, y, x);
              y -= a / b * x;
              return t;
          }
      }
      
      int main(){
          cin >> n;
          while(n -- ){
              cin >> a >> b >> m;
              int t = exgcd(a, m, x, y);
              //b是gcd(a,m)的倍数才有解
              if(b % t != 0) cout << "impossible" << endl;
              else cout << (LL)x * (b / t) % m << endl;//特解乘上倍数
          }
          return 0;
      }
      

• 中国剩余定理

  • 中国剩余定理 (Chinese Remainder Theorem, CRT) 可求解如下形式的一元线性同余方程组（其中 $n_1,n_2,n_3,···n_k$ 两两互质）：

  • $$\{\begin{array}{ll}x& \equiv a_1(\mathrm{mod}{n}_1)\\ x& \equiv a_2(\mathrm{mod}{n}_2)\\ & ⋮\\ x& \equiv a_k(\mathrm{mod}{n}_k)\end{array}$$

  • 过程：
    • 计算所有模数的积
    • 对于第 $i$ 个方程：
      • 计算 $m_i=\frac{n}{n_i}$
      • 计算 $m_i$ 在模 $n_i$ 意义下的逆元 $m_i^{-1}$
      • 计算 $c_i=m_i{m}_i^{-1}$ （不要对 $n_i$ 取模）
    • 方程组在模 $n$ 意义下的唯一解为：$x=\sum _{i=1}^k{a}_i{c}_i(\mathrm{mod}n)$

  • 扩展版中国剩余定理：

    • 当不满足 $n_1,n_2,n_3,···n_k$ 两两互质时，求解线性同余方程组

    • $$\mathrm{①}x\equiv a_1(\mathrm{mod}{n}_1)\iff x=k_1{n}_1+a_1$$

    • $$\mathrm{②}x\equiv a_2(\mathrm{mod}{n}_2)\iff x=k_2{n}_2+a_2$$

    • 由①②知$$k_1{n}_1+a_1=k_2{n}_2+a_2$$
    • 移项$$k_1{n}_1-k_2{n}_2=a_2-a_1$$
    • 其中 $n_1,n_2,a_2,a_1$ 已知，$k_1,k_2$ 未知，可以用扩展欧几里得算法求出
    • 当 $d=\frac{a_2-a_1}{gcd(n_1,n_2)}\ne 0$ 时无解
    • 当 $d=\frac{a_2-a_1}{gcd(n_1,n_2)}=0$ 时，先求出方程

    • $$k_1^{‘}{n}_1-k_2^{‘}{n}_2=gcd(n_1,n_2)$$

    • 的特解，求出 $k_1^{‘},k_2^{‘}$ 后

    • $$k_1^{\ast }=k_1^{‘}·d,k_2^{\ast }=k_2^{‘}·d$$

    • 此时 $k_1^{\ast },k_2^{\ast }$ 就是原方程的特解
    • 观察原方程可知通解为

    • $$k_1=k_1^{\ast }+k·\frac{n_2}{d},k_2=k_2^{\ast }+k·\frac{n_1}{d}(\mathrm{其}\mathrm{中}k\mathrm{为}\mathrm{整}\mathrm{数})$$

    • 将通解代入①式得

    • $$x=(k_1^{\ast }+k·\frac{n_2}{d})n_1+a_1=(k_1^{\ast }{n}_1)+a_1+k·\frac{n_1{n}_2}{d}=(k_1^{\ast }{n}_1)+a_1+k·lcm(n_1,n_2)$$

    • 所以

    • $$\mathrm{③}\mathrm{③}x=k·lcm(n_1,n_2)+(k_1^{\ast }{n}_1)+a_1$$

    • 其中 $k_1^{\ast },n_1,a_1,lcm(n_1,n_2)$ 已知，$k$ 未知
    • 令 $n_2^{‘}=lcm(n_1,n_2),a_2^{‘}=(k_1^{\ast }{n}_1)+a_1$ 所以③式化为

    • $$\mathrm{③}x=k·n_2^{‘}+a_2^{‘}$$

    • 所以将①②式合并后的③式仍与方程组方程相似，于是可以循环将方程组合并为一个式子，通过最后一个式子即可用扩展欧几里得算法求出x
    • 代码：

      const int N = 30;
      typedef long long LL;
      int n;
      
      LL exgcd(LL a, LL b, LL &x, LL &y){
          if(!b){
              x = 1, y = 0;
              return a;
          }
          LL d = exgcd(b, a % b, y, x);
          y -= a / b * x;
          return d;
      }
      int main(){
          cin >> n;
          LL a1, m1;
          bool flag = true;
          cin >> a1 >> m1;//用于存储更新的a与m
          
          for(int i = 1; i < n; ++ i){ //合并n-1次
              LL an, mn;//接受新的a与m
              LL k1, k2; //存储用欧几里得求的系数
              cin >> an >> mn;
              LL d = exgcd(a1, an, k1, k2); //求出gcd(a1, an)以及系数k1, k2
              if((mn - m1) % d){ //无解判断
                  flag = false;
                  break;
              }
              LL t = an / d; //通解k1 = k1 + k * (an / d), k2 = k2 + k * (a1 / d),题目为了求最小的x，故要让k1为最小解故要k1不断膜上t
              k1 = k1 * (mn - m1) / d; //更新k1
              k1 = (k1 % t + t) % t; //取k1的最小解
              m1 = k1 * a1 + m1; //m更新为k1 * a1 * m1
              a1 = abs(a1 / d * an);  //a1更新为gcd(a1, an)
          }
          
          if(flag){
              cout << (m1 % a1 + a1) % a1; //最后一个式子的就是余数m
          }else cout << -1;
          return 0;
      }