证明：

• 互质：两数的公共因子只有1
• 去掉所有与N有(大于1的)公共因子的数，剩下的数就是与N互质的数
• 对N的所有质因子\(p_k\)，去掉所有\(\underline{质数p_k的倍数}\)(与N有公共因子的数)，\(\underline{每个质数的倍数}\)个数为 \(\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_1p_2} + \frac{N}{p_1p_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_1p_2} + \frac{N}{p_1p_3} + ··· +\frac{N}{p_{n-1}p_n}) - (\frac{N}{p_1p_2p_3} + \frac{N}{p_1p_2p_4} + ··· +\frac{N}{p_{n-2}p_{n-1}p_n})$$
• 对于合数 \(p_i·p_j·p_k·p_m\) 同理，归纳递推可知，所有质数个数为：$$N - \displaystyle\sum^n_{i=1}{\frac{N}{p_i}}+\displaystyle \sum_{1<=i<j<=n}{\frac{N}{p_ip_j}}-\displaystyle \sum_{1<=i<j<k<=n}{\frac{N}{p_ip_jp_k}}+···+(-1)^{2n-1}\displaystyle \sum_{1<=i<j<k<···<=n}{\frac{N}{p_ip_jp_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|+\ldots +(-1)^{n+1}|\cap_{i=1}^n A_i | \]

  • 与N有(大于1的)公共因子的个数 =| \(p_i\) 的倍数的个数| - |\(p_ip_j\) 的倍数的个数| + |\(p_ip_jp_k\) 的倍数的个数 ··· \((-1)^{n+1}\) |\(p_ip_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互质时有$$\varphi(a·b)=\varphi(a)·\varphi(b)$$
• 当i为质数时，$$\varphi(i)=i-1$$
• 当\(\frac{i}{p_j} = 0\)时，\(p_j\)为i的质因子，此时\(p_j·i\)的质因子与i的质因子完全相同，所以 $$\varphi(pj·i)=pj·i·(1-\frac{1}{p_1})(1-\frac{1}{p_2})···(1-\frac{1}{p_n})=p_j·\varphi(i)$$
• 当\(\frac{i}{p_j} \neq 0\)时，\(p_j\)不为i的质因子，此时\(p_j·i\)的质因子比i的质因子多一个\(p_j\)，所以 $$\varphi(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·\varphi(i)·(1-\frac{1}{p_j})=(p_j - 1)·\varphi(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^{\varphi(n)}\pmod n\equiv1$$
  • 证明：
    • 若a与b互质，且\(a\pmod b\neq 0\)则\(a\pmod b\) 也与b互质
    • 设\(a_1,a_2,a_3···a_{\varphi(n)}\) 为1~n中的所有与n互质的数
    • 每个数同时乘上a得\(a·a_1,a·a_2,a·a_3···a·a_{\varphi(n)}\) 由于a也与n互质，所以乘a后所得的数也全部与n互质
    • 将每个数mod n，取模后的每个数都在1~n范围内，并且仍然都与n互质，显然取模后的这几个数就是原来的几个数，只不过数的顺序可能发生了变化
    • 将所有数相乘，则有$$a·a_1·a·a_2·a·a_3···a·a_{\varphi(n)} \pmod n=a_1·a_2·a_3···a_{\varphi(n)}$$
    • 所以有：$$a^{\varphi(n)}·(a_1a_2a_3···a_{\varphi(n)})\pmod n \equiv (a_1·a_2·a_3···a_{\varphi(n)})$$即$$a^{\varphi(n)}\pmod n\equiv1$$
    • 特别地，当n为质数时，\(a^{n-1}\pmod n\equiv1\) 为费马定理

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 \pmod p\) ，则x为b的逆元
• 两边乘b有$$a=a·b·x \pmod p$$
• 所以有$$b·x=1 \pmod p$$
• 若p是质数，有费马定理$$b^{p-1} = 1 \pmod p$$则$$x=b^{p-2}$$
• 若p不是质数，有欧拉定理$$b^{\varphi(p)}=1 \pmod p$$则$$x=b^{\varphi(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 \neq 0\) 时,因为

而

故而 $$x=y_1,\quad 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 \pmod 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 \pmod m\)有解的充要条件为b为\(gcd(a,m)\)的倍数
• 当 \(b \neq k'·gcd(a,m)\) 时(\(k'\) 为整数)，方程无解
• 当 \(b = k'·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;
  }
  

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

扩展版中国剩余定理：

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

• \[①\quad x\equiv a_1\pmod {n_1} \Leftrightarrow x=k_1n_1+a_1 \]

• \[②\quad x \equiv a_2\pmod {n_2} \Leftrightarrow x=k_2n_2+a_2 \]

• 由①②知$$k_1n_1+a_1=k_2n_2+a_2$$
• 移项$$k_1n_1-k_2n_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)} \neq 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^*=k_1^`·d,\quad k_2^*=k_2^`·d \]

• 此时 \(k_1^*,k_2^*\) 就是原方程的特解
• 观察原方程可知通解为

• \[k_1=k_1^*+k·\frac{n_2}{d},\quad k_2=k_2^*+k·\frac{n_1}{d}(其中k为整数) \]

• 将通解代入①式得

• \[x=(k_1^*+k·\frac{n_2}{d})n_1+a_1=(k_1^*n_1)+a_1+k·\frac{n_1n_2}{d}=(k_1^*n_1)+a_1+k·lcm(n_1,n_2) \]

• 所以

• \[③③\quad x=k·lcm(n_1,n_2)+(k_1^*n_1)+a_1 \]

• 其中 \(k_1^*,n_1,a_1,lcm(n_1,n_2)\) 已知，\(k\) 未知
• 令 \(n_2^`=lcm(n_1,n_2),\quad a_2^`=(k_1^*n_1)+a_1\) 所以③式化为

• \[③\quad 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;
  }