二次项展开（矩阵加速）

#include <cmath>
#include <stack>
#include <queue>
#include <vector>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <iostream>
#include <algorithm>
#define LL long long
#define ULL unsigned long long
using namespace std;

template <typename T> int read (T &x) {x = 0; T f = 1;char tem = getchar ();while (tem < '0' || tem > '9') {if (tem == '-') f = -1;tem = getchar ();}while (tem >= '0' && tem <= '9') {x = (x << 1) + (x << 3) + tem - '0';tem = getchar ();}x *= f; return 1;}
template <typename T> void write (T x) {if (x < 0) {x = -x;putchar ('-');}if (x > 9) write (x / 10);putchar (x % 10 + '0');}
template <typename T> T Max (T x, T y) { return x > y ? x : y; }
template <typename T> T Min (T x, T y) { return x < y ? x : y; }
template <typename T> T Abs (T x) { return x > 0 ? x : -x; }

const int Maxn = 40;
const LL Mod = 1e9 + 7;

struct Matrix {
	int n, m;
	LL c[Maxn + 5][Maxn + 5];
	
	Matrix () { memset (c, 0, sizeof c); }
	
	void init () {
		for (int i = 0; i <= n; i++)
			c[i][i] = 1;
	}
	
	void cpy (LL a[][Maxn + 5]) {
		for (int i = 0; i < Maxn + 5; i++) {
			for (int j = 0; j < Maxn + 5; j++) {
				c[i][j] = a[i][j];
			}
		}
	}
	
	void read () {
		for (int i = 1; i <= n; i++) {
			for (int j = 1; j <= m; j++) {
				scanf ("%lld", &c[i][j]);
				c[i][j] %= Mod;
				if (c[i][j] < 0) c[i][j] += Mod;
			}
		}
	}
	
	void write () {
		for (int i = 1; i <= n; i++) {
			for (int j = 1; j <= m; j++) {
				printf ("%lld ", c[i][j]);
			}
			puts ("");
		}
	}
	
	Matrix operator * (const Matrix x) const {
		Matrix ans; ans.n = n; ans.m = x.m;
		for (int i = 1; i <= ans.n; i++) {
			for (int k = 1; k <= m; k++) {
				for (int j = 1; j <= ans.m; j++) {
					ans.c[i][j] += c[i][k] * x.c[k][j];
					ans.c[i][j] %= Mod;
					if (ans.c[i][j] < 0) ans.c[i][j] += Mod; 
				}
			}
		}
		return ans;
	}
	
	Matrix operator + (const Matrix x) const {
		Matrix ans; ans.n = n; ans.m = m;
		for (int i = 1; i <= n; i++)  {
			for (int j = 1; j <= m; j++) {
				ans.c[i][j] = c[i][j] + x.c[i][j];
			}
		}
		return ans;
	}
	
	Matrix operator - (const Matrix x) const {
		Matrix ans; ans.n = n; ans.m = m;
		for (int i = 1; i <= n; i++) {
			for (int j = 1; j <= m; j++) {
				ans.c[i][j] = c[i][j] + x.c[i][j];
			}
		}
		return ans;
	}
};

LL fac[Maxn + 5], inv_fac[Maxn + 5];
void exgcd (LL a, LL b, LL &x, LL &y) {
    if (b == 0) {
        x = 1; y = 0;
        return;
    }
    exgcd (b, a % b, y, x);
    y -= a / b * x;
}
LL inv (LL a) {
    LL x, y;
    exgcd (a, Mod, x, y);
    x = ((x % Mod) + Mod) % Mod;
    return x;
}
void Init () {
	fac[0] = 1;
	for (int i = 1; i <= Maxn; i++) {
		fac[i] = fac[i - 1] * i;
		fac[i] %= Mod;
	}
	inv_fac[Maxn] = inv (fac[Maxn]);
	for (int i = Maxn - 1; i >= 0; i--) {
		inv_fac[i] = inv_fac[i + 1] * (i + 1) % Mod;
//		printf ("i = %d, fac[%d] = %lld, inv_fac[%d] = %lld\n", i, i, fac[i], i, inv_fac[i]);
	}
}
LL C (LL a, LL b) {
    return fac[b] * inv_fac[a] % Mod * inv_fac[b - a] % Mod;
}

int n, k;

int main () {
	read (n); read (k);
	Init ();
	Matrix a, b;
	a.n = 1; a.m = b.n = b.m = k + 3;
	for (int i = 1; i <= a.m; i++) {
		a.c[1][i] = 1;
	}
	b.c[2][1] = 1;
	b.c[1][2] = 1;
	b.c[2][2] = 1;
	b.c[3][3] = 1;
	for (int j = 4; j <= b.m; j++)
		for (int i = 3; i <= j; i++)
			b.c[i][j] = C (i - 3, j - 3);
	
	printf ("A:\n");
	a.write ();
	printf ("B:\n");
	b.write ();
	for (int i = 1; i <= 10; i++) {
		printf ("step = %d\n", i);
		a.write ();
		a = a * b;
	}
	return 0;
}