AtCoder Regular Contest 182 C. Sum of Number of Divisors of Product（矩阵快速幂+多项式贡献）

题目

求所有长为k的序列(k∈[1,n](n<=1e18))的score，第i个数可以在[1,m](m<=16)中任取

一个确定了的序列的score，为这个序列所有数的乘积的因子数

思路来源

官方题解 Editorial - AtCoder Regular Contest 182

poj3233 Matrix Power Series(矩阵快速幂 求前1-k次幂矩阵和)_矩阵前k次幂之和-CSDN博客

题解

两个题的拼凑，多项式贡献经典题+如何把矩阵快速幂从只求第n项改成求前缀和

先是第一部分，多项式贡献，

首先[1,16]内的质因子，只有2、3、5、7、11、13这6个

而实际要求的式子是一个(1+x1)*(1+x2)*(1+x3)*(1+x4)*(1+x5)*(1+x6)，

其中xi是第i个质因子出现的个数

我们可以维护2*2*2*2*2*2个多项式中的项，共64项，

比如，000111就可以表示当前(1+x1)*(1+x2)*(1+x3)的值，

那么，矩阵快速幂，每次从dp[i]转移到dp[i+1]上的时候，

枚举第i+1位填什么，j从[1,m]枚举，

对j做质因数分解，看每个质因子出现的个数，会得到一个长为6的向量

这里用长为3的状态000111举例，且枚举的j分解后得到的质因数向量是(a,b,0,0,0,0)，

对于(x1+a)*(x2+b)*(x3)来说，拆解开，发现是等于x1x2x3+a*x2*x3+b*x1x3+a*b*x3

那么(x3)这个状态转移到(x1x2x3)这个状态，系数就是a*b，

同理有，(x1x3)状态转移到(x1x2x3)这个状态，系数是b，

(x2x3)状态转移到(x1x2x3)这个状态，系数是a，

(x1x2x3)状态转移到(x1x2x3)这个状态，系数是1，

在转移矩阵上加上对应的数即可

不难发现，从y状态转移到x状态，一定需要满足y状态是x状态的子集，

且多出来的增量状态（也就是x^y这部分），确定了转移系数，

对于一个固定的j来说，j对应的质因数向量，在增量状态这些位上的乘积，就是转移系数

枚举所有j，求一个加和，得到转移一次的系数矩阵b，

[向量]=[转移矩阵]*[初始向量]，

最开始的时候x1=x2=x3=x4=x5=x6=0，也就是初始向量63这个状态的值为1

对b做矩阵快速幂b^n即可，但是这样只能求得第n项的值，

现在需要求前缀和，就需要用到第二部分

第二部分则是一个套路，这里63位置的值是我们要求的dp[i]，

我们新增一个状态64，表示前缀和，

最后只求前缀和就可以了，所以不用再管63，直接置初始向量64这个状态初始值为1

钦定每次转移的时候63都可以往64转移，也就是sum[i+1]+=dp[i]

64也可以往64转移，也就是sum[i+1]+=sum[i]

这样最后的sum[n+1]就是所求，但是此处会多加dp[0]这一种空方案，减去即可

从dag转移的思路上，之前poj3233的套路就好理解很多了

代码

#include <bits/stdc++.h>
#include<iostream>
#include<cstdio>
#include<vector>
#include<map>
using namespace std;
#define rep(i,a,b) for(int i=(a);i<=(b);++i)
#define per(i,a,b) for(int i=(a);i>=(b);--i)
typedef long long ll;
typedef double db;
typedef pair<int,int> P;
#define fi first
#define se second
#define pb push_back
#define dbg(x) cerr<<(#x)<<":"<<x<<" ";
#define dbg2(x) cerr<<(#x)<<":"<<x<<endl;
#define SZ(a) (int)(a.size())
#define sci(a) scanf("%d",&(a))
#define pt(a) printf("%d",a);
#define pte(a) printf("%d\n",a)
#define ptlle(a) printf("%lld\n",a)
#define debug(...) fprintf(stderr, __VA_ARGS__)
using namespace std;
using namespace std;
typedef long long ll;
const int mod=998244353;
struct mat {
    static const int MAXN=65;
    ll c[MAXN][MAXN];
    int m, n;
    mat(){
    	memset(c, 0, sizeof(c));
    	m=n=MAXN;
    }
    mat(int a, int b) : m(a), n(b) {
        memset(c, 0, sizeof(c));
    }
    void clear(){
		memset(c, 0, sizeof(c));
    }
    mat operator * (const mat& temp) {
        mat ans(m, temp.n);
        for (int i = 0; i < m; i ++)
            for (int j = 0; j < temp.n; j ++){
                for (int k = 0; k < n; k ++){
                    ans.c[i][j] += c[i][k] * temp.c[k][j];
                    ans.c[i][j]%=mod;
                }
            }
        return ans;
    }
    mat operator ^(ll n){
        mat M(*this),ans(M.m, M.m);
    	for (int i = 0; i < M.m; i ++)
        ans.c[i][i] = 1;
   		while (n > 0) {
            if (n & 1) ans = ans * M;
            M = M * M;
            n >>= 1;
    	}
        return ans;
	}
}b;
ll n;
int m,a[10]={2,3,5,7,11,13};
int c[17][10],ans;
int has;
int main(){
    scanf("%lld%d",&n,&m);
    rep(i,0,5){
        rep(j,1,m){
            int p=j;
            while(p%a[i]==0)p/=a[i],c[j][i]++;
        }
    }
    rep(i,1,m){
        rep(x,0,63){
            rep(y,0,63){//y->x
                if((y&x)!=y)continue;
                int xs=1;
                rep(z,0,5){
                    int u=x>>z&1,v=y>>z&1;
                    if(!u)continue;
                    if(!v)xs=1ll*xs*c[i][z]%mod;
                }
                b.c[y][x]+=xs;
            }
        }
    }
    b.c[63][64]=1;
    b.c[64][64]=1;
    b=b^(n+1);
    rep(i,0,63){
        ans=(ans+b.c[i][64])%mod;
    }
    ans=(ans+mod-1)%mod;
    pte(ans);
    return 0;
}