https://szkopul.edu.pl/problemset/problem/yGC3v6-xidN3ZW6QNlBAFZSU/site/?key=statement
L(i+1)^x * L(i)^yを求める際に必要な情報を考えてみます。これを(x,y)と表記します。
L(i+1)^x = (L(i) + L(i - 1) + 1) ^ x = Σ(a+b+c <= x)(x!/a!b!c!)L(i)^a * L(i)^b より
L(i+1)^x * L(i)^y = Σ(a+b+c <= x)(k!/a!b!c!)L(i)^(a+y) * L(i)^bとなります。
整理すると、(x,y)を求めるためには(a+y,b)(a+b <= x)が必要となり、a+y=X,b=Yとおくと、XY平面上でY<=x && y<=X<=x+y-Yを満たす格子点に相当する点について情報があればいいです。
上のXYの条件は図示すると(y,x),(y,0),(x+y,0)を頂点とする二等辺三角形の領域と対応し、よってX+Y<=kを満たすような(X,Y)(すなわちL(i)^X*L(i-1)^Y)について情報を持っておけば、(x,y)(すなわちL(i+1)^x*L(i)^Y)についてもx+y<=kをみたす領域について求めることができます。なので行列は(x,y)の情報((k+1)*(k+2)/2)個と、I(n)=Σ(i<=n)L(i)についての情報1個の計M=((k+1)*(k+2)/2+1)個持って累乗すればI(n)を求められます。k<=13のときM=106となるのでO(M^3 * logN)でなんとか間に合います。
M*2とかだとTLEします。A+A^2+A^3...を求める関数が用意してあったので、最初脳死でI(n)をこれを使って求めましたが、行列のサイズがM*2となるので見事にTLEしました。数列の和を求める際は普通に1つだけ要素を増やすだけでいいですね。行列の形にできないと思ったら1つ要素を追加することを考えるのが良さそうです。
void show(ll a) { vi vec; rep(i, 0, 9) { vec.pb(a % 10); a /= 10; } rer(i, 9, 0) cout << vec[i]; cout << endl; } ull N; int K; ll fac[20]; int id[200][200]; void solve() { cin >> N >> K; fac[0] = fac[1] = 1; rep(i, 2, K + 1) fac[i] = fac[i - 1] * i; int M = (K + 1) * (K + 2) / 2 + 1; mat A(M, vl(M, 0)); int cnt = 0; rep(x, 0, K + 1) { rep(y, 0, K + 1) { if(x + y > K) continue; id[x][y] = cnt; cnt++; } } rep(x, 0, K + 1) { rep(y, 0, K + 1) { if(x + y > K) continue; int pid = id[x][y]; rep(a, 0, x + 1) { rep(b, 0, x + 1) { int c = x - a - b; if(c < 0) continue; ll mul = fac[x] / (fac[a] * fac[b] * fac[c]); int tid = id[a + y][b]; A[pid][tid] = mul; } } } } int x = K, y = 0; rep(a, 0, x + 1) { rep(b, 0, x + 1) { int c = x - a - b; if(c < 0) continue; ll mul = fac[x] / (fac[a] * fac[b] * fac[c]); int tid = id[a + y][b]; A[M - 1][tid] = mul; A[M - 1][M - 1] = 1; } } mat B = mat_pow(A, N - 1); ll res = 0; rep(i, 0, M) { if(i != M - 1) ADD(res, B[M - 1][i]); else ADD(res, B[M - 1][i] * 2 % mod); } show(res); }
uLLにしないといけないのもさりげない罠。