[数论第三节]高斯消元法/求组合数/卡特兰数

• 高斯消元

  • 求解含有n个未知数，n个方程的多元线性方程组 O(n^3)
  • 初等行变换：
    • 某行乘以一个非零数
    • 交换两行
    • 某行加上另一行的若干倍
  • 利用初等行变换将方程组化为上三角矩阵
  • 解的情况：
    • 完美阶梯型：唯一解
    • 非完美阶梯型：
      • 0 == 非0：无解
      • 0 == 0：无穷解
  • 步骤:
    • 枚举每一列
      • 找到这一列系数的绝对值最大的一行
      • 将这一行与第一行交换
      • 将改行的第一个数变成一(方程两边同乘某数)
      • 把下面所有行的当前列的系数消成0(某行加上第一行的若干倍)
  • 代码：

    const int N = 110;
    const double esp = 1e-6; //x < esp,则x = 0,否则 x != 0,由于浮点数精度问题,0可能是0.000001
    double a[N][N];
    int n;
    
    int gauss(){
    	int r, c;
    	for(r = 1, c = 1; c <= n; ++ c){//遍历每一列
    		int t = r;
    		for(int i = r; i <= n; ++ i) //找当前列的系数最大的行
    			if(fabs(a[t][c]) < fabs(a[i][c]))
    				t = i; //记录最大行
    		
    		if(fabs(a[t][c]) < esp) continue; //最大行系数为0说明该列系数均为0
    		
    		for(int i = c; i <= n + 1; ++ i) swap(a[r][i], a[t][i]);//否则交换第一行与系数最大行
    		for(int i = n + 1; i >= c; -- i) a[r][i] /= a[r][c]; //将第一行当前列系数变为1
    		
    		for(int i = r + 1; i <= n; ++ i) //将下面所有行的当前列的系数变为0
    			if(fabs(a[i][c]) > esp) //当前列系数非0
    				for(int j = n + 1; j >= c; -- j) 
    					a[i][j] -= a[r][j] * a[i][c];
    
    		r ++ ;
    	}
    	
    	if(r <= n){
    		for(int i = r; i <= n; ++ i)
    			if(fabs(a[i][n + 1]) > esp) return 0; //无解
    		return 1; //无数解
    	}
    	
    	for(int i = n; i >= 1; -- i) //从后往前求解
    		for(int j = i + 1; j <= n; ++ j)
    			a[i][n + 1] -= a[i][j] * a[j][n + 1];
    			
    	return 2; //唯一解
    }
    int main(){
    	cin >> n;
    	
    	for(int i = 1; i <= n; ++ i)
    		for(int j = 1; j <= n + 1; ++ j)
    			cin >> a[i][j];
    	
    	int t = gauss();
    	
    	if(t == 1) cout << "Infinite group solutions";
    	else if(t == 0) cout << "No solution";
    	else if(t == 2){
    		for(int i = 1; i <= n; ++ i) 
    			printf("%.2f\n", a[i][n + 1]);
    	}
    	
    	return 0;
    }
    

• 求组合数

  • \(C_a^b = \frac{a(a-1)(a-2)···(a-b+1)}{b!}=\frac{a!}{b!(a-b)!}\)
  • 一般直接用公式求组合数容易超时
  • 故可以通过预处理的方式，用递推式求出所有可能的组合数
  • 递推式：\(C_a^b = C_{a-1}^{b-1} + C_{a-1}^{b}\) \(O(N^2)\)
  • 当a，b <= 2e3，询问n = 1e4次时
  • 代码：

    //询问10000次
    //a,b <= 2000
    const int N = 2e3 + 10;
    const int mod = 1e9 + 7;
    int c[N][N];
    int n;
    
    void init(){
        for(int i = 0; i < N; ++ i)
            for(int j = 0; j <= i; ++ j)
                if(!j) c[i][j] = 1;
                else c[i][j] = (c[i - 1][j - 1] + c[i - 1][j]) % mod;
    }
    

  • 当a，b <= 1e5，询问n = 1e4次，模数p = 1e9 + 7时
    • 开二维数组会爆内存，于是利用公式，预处理出阶乘以及阶乘的逆 \(O(Nlog(N))\)
    • 代码：

      typedef long long LL;
      const int N = 1e5 + 10;
      const int mod = 1e9 + 7;
      int fact[N], infact[N];
      int n;
      
      //快速幂
      int qmi(int a, int b, int p){
          int res = 1;
          while(b){
              if(b & 1) res = (LL)res * a % p;
              a = (LL)a * a % p;
              b >>= 1;
          }
          return res;
      }
      
      //预处理各阶乘以及阶乘的逆
      void init(){
          fact[0] = infact[0] = 1;
          for(int i = 1; i < N; ++ i){
              fact[i] = (LL)fact[i - 1] * i % mod; //递推式n! = (n - 1)! * n;
              infact[i] = (LL)infact[i - 1] * qmi(i, mod - 2, mod) % mod; //递推式(n!)^-1 = ((n - 1)!)^-1 * n^-1，取模意义下的逆就是逆元
          }
      }
      
      int main(){
          cin >> n;
          init();
          while(n -- ){
              int a, b;
              cin >> a >> b;
              cout << (LL)fact[a] * infact[b] % mod * infact[a - b] % mod << endl;//两个int最大值相乘不会爆longlong，但是三个相乘会爆
          }
          return 0;
      }
      

  • 当a，b <= 1e18，询问n = 2e5次，模数p不定时
    • 两数相乘会爆longlong，利用卢卡斯定理预处理组合数 \(O(plog(N)log(p))\)
    • 定理：\(C_a^b\equiv C_{a \mod p}^{b \mod p}·C_{a/p}^{b/p} \pmod p\)
    • 代码：

      typedef long long LL;
      int p;
      int n;
      
      //快速幂求逆元
      int qmi(int a, int b){
          int res = 1;
          while(b){
              if(b & 1) res = (LL)res * a % p;
              a = (LL)a * a % p;
              b >>= 1;
          }
          return res;
      }
      
      int C(int a, int b){ //求组合数C_a^b
          if(a < b) return 0; //可以不要
          int x = 1, y = 1; //分子分母
          for(int i = 0; i < b; ++ i){
              x = (LL)x * (a - i) % p; //求分子
              y = (LL)y * (i + 1) % p; //求分母
          }
          return (LL)x * qmi(y, p - 2) % p;
      }
      
      //卢卡斯定理
      int lucas(LL a, LL b){
          if(a < p && b < p) return C(a, b); //不必用卢卡斯，可以直接利用公式求出
          return (LL)C(a % p, b % p) * lucas(a / p, b / p) % p; //用卢卡斯定理
      }
      
      int main(){
          cin >> n;
          while(n -- ){
              LL a, b;
              cin >> a >> b >> p;
              cout << lucas(a, b) << endl;
          }
          return 0;
      }
      

  • 当a，b <= 5e5，询问一次，不取模时
    • 此时组合数较大，需要用高精度存储结果
    • 阶乘中素数的个数：
      • \(v_p(n!) = \sum_{i=1}^{\infty} \left\lfloor \frac{n}{p^i} \right\rfloor\)
      • 证明：将 \(n!\) 记为 \(1\times 2\times \cdots \times p\times \cdots \times 2p\times \cdots \times \lfloor n/p\rfloor p\times \cdots \times n\) 那么其中 \(p\) 的倍数有 \(p\times 2p\times \cdots \times \lfloor n/p\rfloor p=p^{\lfloor n/p\rfloor }\lfloor n/p\rfloor !\) 然后在 \(\lfloor n/p\rfloor !\) 中继续寻找 \(p\) 的倍数即可，这是一个递归的过程
    • 代码：

      const int N = 5e3 + 10;
      int primes[N], cnt;
      int st[N];
      int sum[N];
      
      //线性筛，1-n中质数就是n!中的质因子
      void get_primes(int a){
          for(int i = 2; i <= a; ++ i){
              if(!st[i]) primes[ ++ cnt] = i;
              for(int j = 1; primes[j] <= a / i; ++ j){
                  st[primes[j] * i] = 1;
                  if(i % primes[j] == 0) break;
              }
          }
      }
      
      //求a!中质数p的个数
      int get(int a, int p){
          int res = 0;
          while(a){
              res += a / p;
              a /= p;
          }
          return res;
      }
      
      //高精度乘法
      vector<int> mul(vector<int> &A, int b){
          vector<int> c;
          int t = 0;
          for(int i = 0; i < A.size() || t; ++ i){
              if(i < A.size()) t += A[i] * b;
              c.push_back(t % 10);
              t /= 10;
          }
          while(c.size() > 1 && c.back() == 0) c.pop_back();
          return c;
      }
      
      int main(){
      	int a, b;
          cin >> a >> b;
          get_primes(a);
          for(int i = 1; i <= cnt; ++ i){
              int p = primes[i];
              sum[i] = get(a, p) - get(b, p) - get(a - b, p); //a!中质因子个数-b!中质因子个数-(a-b)!中质因子个数
          }
          
          vector<int> res;
          res.push_back(1);
          
          for(int i = 1; i <= cnt; ++ i) //累乘
              for(int j = 1; j <= sum[i]; ++ j)
                  res = mul(res, primes[i]);
          
          for(int i = res.size() - 1; i >= 0; -- i)//输出
              cout << res[i];
          
          return 0;
      }
      

• 卡特兰数

  • \(C_{2n}^n-C_{2n}^{n-1}=\frac{1}{n+1}·C_{2n}^n\)
  • 用于解决序列排列等问题
  • AcWing 889. 满足条件的01序列 - AcWing
  • 代码：

    #include <iostream>
    
    using namespace std;
    
    typedef long long LL;
    const int N = 1e5 + 10;
    const int p = 1e9 + 7;
    int n;
    
    int qmi(int a, int k){
        int res = 1;
        while(k){
            if(k & 1) res = (LL)res * a % p;
            a = (LL)a * a % p;
            k >>= 1;
        }
        return res;
    }
    
    int main(){
        cin >> n;
        
        int res = 1;
        
        int a = 2 * n, b = n;
        int x = 1, y = 1;
        for(int i = 0; i < b; ++ i){
            x = (LL)x * (a - i) % p;
            y = (LL)y * (i + 1) % p;
        }
        
        res = (LL)x * qmi(y, p - 2) % p;
        res = (LL)res * qmi(n + 1, p - 2) % p;
        cout << res;
        
        return 0;
    }