[数论第三节]高斯消元法/求组合数/卡特兰数
  snDytuRCG8Kr 2023年11月01日 70 0
  • 高斯消元

    • 求解含有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;
      }
      
【版权声明】本文内容来自摩杜云社区用户原创、第三方投稿、转载,内容版权归原作者所有。本网站的目的在于传递更多信息,不拥有版权,亦不承担相应法律责任。如果您发现本社区中有涉嫌抄袭的内容,欢迎发送邮件进行举报,并提供相关证据,一经查实,本社区将立刻删除涉嫌侵权内容,举报邮箱: cloudbbs@moduyun.com

  1. 分享:
最后一次编辑于 2023年11月08日 0

暂无评论

推荐阅读
  jTMfQq5cr55P   2024年05月17日   42   0   0 算法与数据结构
  jTMfQq5cr55P   2024年05月17日   38   0   0 算法与数据结构
snDytuRCG8Kr