高斯消元法 矩阵快速幂

高斯消元

高斯消元是解 \(n\) 个 \(n\) 次方程的解, 可以判断是否无解, 是否有无数个解, 或者输出那唯一组解。

我们解二元一次方程组运用的是加减消元法, 在高斯消元中同样适用。

我们要把它解成

• 只有第一个方程的第一项的系数非0。
• 只有第二个方程的第二项的系数非0。
  \(\dots\)
• 只有第 \(n\) 个方程的第 \(n\) 项的系数非 0。

对于每一项, 我们可以枚举哪一个在这一项的系数非0, 然后用这个方程对其他方程做加减法, 把其他系数消掉。

对于无解和有无数解的情况, 如果他的某一项经过前面的某些操作后没有一个不为 0, 有两种情况 : 无解, 有无数解。 如果你这是就判断这个, 你会被以下数据 hack

3
0 1 1
0 0 0

这个数据有无数个解, 应为你可以调换他们的顺序。

处理方法 : 如果遇到:经过前面的某些操作后没有一个不为 0, 先不管他, 后面有可能某一项会用到他。

例题 : ，洛谷P2455

#include<bits/stdc++.h>

using namespace std;

const int N = 105;
const double EPS = 1e-12;

double op;
vector<double>e[N];
int n, flag;

int main(){
  cin >> n;
  for(int i = 1; i <= n; ++i){
    e[i].resize(n + 3);
    for(int j = 1; j <= n + 1; j++){
      cin >> e[i][j];
    }
  }
  for(int i = 1; i <= n; ++i){
    for(int j = 1; j <= n; ++j){
      flag = 1;
      for(int k = 1; k < i; ++k){
        if(fabs(e[j][k]) >= EPS){
          flag = 0;
        }
      }
      if(fabs(e[j][i]) >= EPS && flag){
        swap(e[i], e[j]);
        break;
      }
    }
    if(fabs(e[i][i]) < EPS){
      continue;
    }
    for(int j = 1; j <= n; ++j){
      if(i != j){
        op = e[j][i] / e[i][i];
        for(int k = 1; k <= n + 1; ++k){
          e[j][k] -= e[i][k] * op;
        }
      }
    }
  }
  for(int i = 1; i <= n; ++i){
    if(fabs(e[i][i]) < EPS && fabs(e[i][n + 1]) >= EPS){
      cout << -1;
      return 0;
    }
  }
  for(int i = 1; i <= n; ++i){
    if(fabs(e[i][i]) < EPS && fabs(e[i][n + 1]) < EPS){
      cout << 0;
      return 0;
    }
  }
  for(int i = 1; i <= n; ++i){
    cout << "x" << i << "=" << fixed << setprecision(2) << e[i][n + 1] / e[i][i] << '\n';
  }
  return 0;
}

矩阵快速幂

可以通过式子不断的迭代, 算出经过 \(2^x\) 次操作后的通用公式。

例题 : luoguP1962

推式子

\(F(n), F(n - 1) \to F(n + 1)， F(n)\)

\(F(n + 1) = F(n) + F(n - 1)\)

\(F(n) = F(n)\)

\(F(n + 1), F(n) \to F(n + 2), F(n + 1)\)

\(F(n + 2) = F(n + 1) + F(n) = 2F(n) + F(n - 1)\)

\(F(n + 1) = F(n + 1) = F(n) + F(n - 1)\)

\(F(n + 1), F(n + 1) \to F(n + 4), F(n + 3)\)

\(F(n + 4) = 2F(n + 2) + F(n + 1) = 5F(n) + 3F(n - 1)\)

\(F(n + 3) = F(n + 2) + F(n + 1) = 3F(n) + 2F(n - 1)\)

\(ax_1 + bx_2 = y_1\)

\(cx_1 + dx_2 = y_2\)

\(ey_1 + fy_2 = z_1\)

\(gy_1 + hy_2 = z_2\)

\((ea + fc)x_1 + (eb + fd)x_2 = z_1\)

\((ga + hc)x_1 + (gb + hd)x_2 = z_2\)

不断的迭代就可以了。

#include<bits/stdc++.h>

using namespace std;

const int mod = 1e9 + 7;

int mp[3][3], a[3][3], now0, now1, f[3];

long long n;

int main(){
  cin >> n;
  if(n <= 2){
    cout << 1;
    return 0;
  }
  n -= 2;
  f[0] = f[1] = 1;
  mp[0][0] = mp[0][1] = mp[1][0] = 1;
  for(; n; n >>= 1){
    if(n & 1){
      now0 = (1ll * f[0] * mp[0][0] % mod + 1ll * f[1] * mp[0][1] % mod) % mod;
      now1 = (1ll * f[0] * mp[1][0] % mod + 1ll * f[1] * mp[1][1] % mod) % mod;
      f[0] = now0, f[1] = now1;
    }
    for(int i = 0; i <= 1; ++i){
      for(int j = 0; j <= 1; ++j){
        a[i][j] = 0;
        for(int k = 0; k <= 1; ++k){
          a[i][j] = (a[i][j] + 1ll * mp[i][k] * mp[k][j] % mod) % mod;
        }
      }
    }
    swap(a, mp);
  }
  cout << f[0];
  return 0;
}