• bzoj3240 [NOI2013]矩阵游戏


    bzoj3240 [NOI2013]矩阵游戏

    有一个递推式 (f_{i, j}) 满足:

    [egin{cases}f_{1, 1}=1\f_{i, j}=a imes f_{i, j-1}+b&(j eq1)\f_{i, 1}=c imes f_{i-1, m}+d&(i eq1)end{cases} ]

    (f_{n, m} mod (10^9+7))

    (n, mleq10^{10^6}, a, b, c, din[1, 10^9])

    矩阵加速,数论


    易知上式是有循环节的,循环节长度为 (P-1) ,因此 (f_{n, m}=f_{nmod (P-[a eq1]), mmod (P-[c eq1])}) (需特判 (a=1)(c=1) 的情况 )

    对于每一行,很容易就能构造出 (egin{bmatrix}f_{i, j}&bend{bmatrix}egin{bmatrix}a&0\1&1end{bmatrix}=egin{bmatrix}f_{i, j+1}&bend{bmatrix}) ,即 (egin{bmatrix}f_{i, 1}&bend{bmatrix}egin{bmatrix}a&0\1&1end{bmatrix}^{m-1}=egin{bmatrix}f_{i, m}&bend{bmatrix})

    但是相邻行的转移貌 (egin{bmatrix}f_{i,m}&dend{bmatrix}egin{bmatrix}c&0\1&1end{bmatrix}=egin{bmatrix}f_{i+1,1}&dend{bmatrix}) 似并不能与上式结合……

    (egin{bmatrix}f_{i+1,1}&dend{bmatrix}) 看着很不爽,我们将它化成 (egin{bmatrix}f_{i, m}&bend{bmatrix}egin{bmatrix}c&0\frac{d}{b}&1end{bmatrix}=egin{bmatrix}f_{i, j+1}&bend{bmatrix})

    (M=egin{bmatrix}a&0\1&1end{bmatrix}^{m-1}, N=egin{bmatrix}c&0\frac{d}{b}&1end{bmatrix})

    可以发现最终答案 (A=egin{bmatrix}f_{1, 1}&bend{bmatrix} imes M imes N imes M imes Ncdots imes M) (即从 (f_{i, 1}) 转移到 (f_{i, m}) 再转移到 (f_{i+1, 1}) ,且第 (n) 行不用乘以 (N)

    因此将 (M imes N) 结合一下,原式 (A=egin{bmatrix}f_{1, 1}&bend{bmatrix} imes (M imes N)^{n-1} imes M)

    时间复杂度 (O(log n+log m))

    代码

    #include <bits/stdc++.h>
    using namespace std;
    
    #define rep(i) for (int i = 0; i < 2; i++)
    const int maxn = 1e6 + 10, P = 1e9 + 7;
    int n, m, a, b, c, d;
    char str1[maxn], str2[maxn];
    
    struct matrix {
      int arr[2][2];
    
      matrix() { memset(arr, 0, sizeof arr); }
    
      int* operator [] (const int &pos) {
        return arr[pos];
      }
    } E;
    
    inline int inc(int x, int y) {
      return x + y < P ? x + y : x + y - P;
    }
    
    matrix operator + (matrix a, matrix b) {
      rep(i) rep(j) a[i][j] = inc(a[i][j], b[i][j]);
      return a;
    }
    
    matrix operator * (matrix a, matrix b) {
      matrix res;
      rep(k) rep(i) rep(j) {
        res[i][j] = inc(res[i][j], 1ll * a[i][k] * b[k][j] % P);
      }
      return res;
    }
    
    matrix qp(matrix a, int k) {
      matrix res = E;
      for (; k; k >>= 1, a = a * a) {
        if (k & 1) res = res * a;
      }
      return res;
    }
    
    int qp(int a, int k) {
      int res = 1;
      for (; k; k >>= 1, a = 1ll * a * a % P) {
        if (k & 1) res = 1ll * res * a % P;
      }
      return res;
    }
    
    int main() {
      E[0][0] = E[1][1] = 1;
      scanf("%s %s %d %d %d %d", str1 + 1, str2 + 1, &a, &b, &c, &d);
      int nP = P - (a != 1), mP = P - (c != 1);
      int len1 = strlen(str1 + 1), len2 = strlen(str2 + 1);
      for (int i = 1; i <= len1; i++) {
        n = (10ll * n + str1[i] - 48) % nP;
      }
      for (int i = 1; i <= len2; i++) {
        m = (10ll * m + str2[i] - 48) % mP;
      }
      if (!n) n = nP;
      if (!m) m = mP;
      matrix A, M, N;
      M[1][0] = M[1][1] = N[1][1] = 1;
      N[1][0] = 1ll * d * qp(b, P - 2) % P;
      A[0][0] = 1, A[0][1] = b, M[0][0] = a, N[0][0] = c;
      M = qp(M, m - 1);
      A = A * qp(M * N, n - 1) * M;
      printf("%d", A[0][0]);
      return 0;
    }
    
  • 相关阅读:
    k近邻 KNN
    聚类之k-means
    支持向量机SVM、优化问题、核函数
    [THUSC 2016] 补退选 (Trie树)
    [CQOI2016] 手机号码 (数位dp)
    [CQOI2012] 交换棋子 (费用流)
    [SCOI2016] 背单词 (Trie树)
    [JSOI2009] 球队收益 (费用流)
    [BZOJ1878][SDOI2009] HH的项链 (树状数组)
    [BZOJ2151]种树
  • 原文地址:https://www.cnblogs.com/Juanzhang/p/10788636.html
Copyright © 2020-2023  润新知