AGC182 C - Sum of Number of Divisors of Product

题目大意：求长度为1,2,...N,的好序列因数个数。好序列满足每个元素\(1\leq x \leq M\)

\(N \leq 1e18, M \leq 16\)

很容易想到维护所有好序列质因数的指数+1的乘积。
\(\prod b_i, 1 \leq i \leq 6\).

考虑每个数对这个乘积的影响。假设我们只有2, 3, 5这三个质数，序列末端添加15,
乘积从\(b_1 b_2 b_3\) 变成 \(b_1(b_2+1)(b_3+1)=b_1 b_2 b_3 + b_1 b_2 + b_1 b_3 + b_1\)
这样\(b_1 b_2 b_3\) 任意子集的乘积也需要计算。

所有不包含这些质因子的数不会对乘积产生影响，但是会产生方案数, 比如末尾添加7， 乘积保持不变。

每个数对这些乘积会产生什么样的影响呢？用二进制表示是否包含某个因子，我们有111代表\(b_1 b_2 b_3\)乘积，101代表\(b_1 b_3\)乘积 以此类推，用S表示质因子的集合
可以发现如果x包含\(b_2\)(3)这个质因子，任何包含3的集合S都会从去掉3的集合S^(1<<1)转移,系数为质因子的指数。比如\(x=9\), \(b_1(b_2+2)b_3=b_1 b_2 b_3 + 2 b_1 b_3\)，
\(x=14\), \((b_1 + 1)b_2 b_3=b_1 b_2 b_3 + b_2 b_3\)，对于每一个x,我们都可以得到一个转移矩阵。把所有转移矩阵相加，再利用矩阵快速幂，就可以求得最终的答案。

这样转移系数就满足加法原理，比如末尾的两个选择，9, 15. \(b_1 b_2 b_3\)乘积的变化就是
$b_1(b_2+2)b_3 + b_1 (b_2 + 1) (b_3 + 1) =(b_1 b_2 b_3 + 2 b_1 b_3) + (b_1 b_2 b_3 + b_1 b_3 + b_1 b_2 + b_1) = 2 b_1 b_2 b_3 + 3 b_1 b_3 + b_1 b_2 + b_1 $，

即，枚举1到m，每个数包含3的因子的个数和就是 每个包含3的集合S都会从去掉3的集合S^(1<<1)转移的系数。
对于上述类似\(b_1 b_2 b_3\)从\(b_1\)转移的系数，我们要指数的乘积, 比如末端是\(12=2^2 * 3^1\)，\(b_1\)的系数是\(2 * 1\)
这样求出末端添加x的转移系数后， 枚举S的子集得到dp的转移矩阵系数。

空集代表的意义是有多少个序列。比如末端填3, \((b_2 + 1) = b_2 + 1\). 这里对上一次操作的所有序列都要+1.
矩阵快速幂求前缀和，再添加一行\((0, 0, 0, ..., 1, 1)\)。注意算n+1次才能得到前n个的和。还要消除第0次求和的影响。

#include<bits/stdc++.h>
using namespace std;
#define LL long long
const int mod = 998244353;
int prime[7] = {2, 3, 5, 7, 11, 13};
const int NUM_P = 6;
LL mat[65][65], base[65][65], res[65][65], lala[65][65];
LL coef[65];
void mat_mul(LL A[][65], LL B[][65], LL C[][65]){
  for(int i = 0;i <= 64; i++){
    for(int j = 0;j <= 64; j++){
      LL tmp = 0;
      for(int k = 0;k <= 64; k++){
        (tmp += B[i][k] * C[k][j] % mod) %= mod;
      }
      lala[i][j] = tmp;
    }
  }
  for(int i = 0;i <= 64; i++){
    for(int j = 0;j <= 64; j++){
      A[i][j] = lala[i][j];
    }
  }
}
void ksm(LL x){
  for(int i = 0;i <= 1 << NUM_P; i++)base[i][i] = res[i][i] = 1;
  mat_mul(base, base, mat);
  while(x){
    if(x & 1){
      mat_mul(res, res, base);
    }
    mat_mul(base, base, base);
    x >>= 1;
  }
  for(int i = 0;i <= 64; i++){
    for(int j = 0;j <= 64; j++)
      mat[i][j] = res[i][j];
  }
}
int main(){
  LL n;
  int m;
  cin>>n>>m;
  for(int i = 0;i < 1 << NUM_P; i++){
    for(int j = 1;j <= m; j++){
      int add = 1;
      for(int k = 0;k < NUM_P; k++){
        if(!((i >> k) & 1))continue;
        int tmp = j;
        int cnt = 0;
        while(tmp % prime[k] == 0){
          tmp /= prime[k];
          cnt++;
        }
        add *= cnt;
      }
      coef[i] += add;
    }
  }
  for(int i = 0;i < 1 << NUM_P; i++){
    for(int s = i;s > 0;s = (s - 1) & i){
      mat[i][i^s] += coef[s];
    }
    mat[i][i] += coef[0];
  }
  mat[1<<NUM_P][1<<NUM_P]++;
  mat[1<<NUM_P][(1<<NUM_P)-1]++;
  ksm(n + 1);
  int res = 0;
  for(int i = 0;i < 1 << NUM_P; i++){
    (res += mat[1<<NUM_P][i]) %= mod;
  }
  cout<<(res - 1 + mod) % mod<<endl;
  return 0;
}