[BZOJ 2326] [HNOI2011] 数学作业 【矩阵乘法】

#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <algorithm>
#include <cmath>

using namespace std;

typedef long long LL;

LL n, Mod;

struct Matrix
{
	int x, y;
	LL A[5][5];
	void Clear() {
		memset(A, 0, sizeof(A));
	}
	void SetXY(int a, int b) {
		x = a; y = b;
	}
} M0, M_Ans, M_t;

LL MulNum(LL a, LL b) {
	LL f = a, ret = 0;
	while (b) {
		if (b & 1) {
			ret += f;
			if (ret > Mod) ret %= Mod;
		}
		b >>= 1;
		f <<= 1;
		if (f > Mod) f %= Mod;
	}
	return ret;
}

Matrix Mul(Matrix Ma, Matrix Mb) {
	Matrix ret;
	ret.SetXY(Ma.x, Mb.y);
	ret.Clear();
	for (int i = 1; i <= ret.x; ++i) {
		for (int j = 1; j <= ret.y; ++j) {
			for (int k = 1; k <= Ma.y; ++k) {
				ret.A[i][j] += MulNum(Ma.A[i][k], Mb.A[k][j]);
				ret.A[i][j] %= Mod;
			}
		}
	}
	return ret;
}

Matrix Pow(Matrix Ma, LL b) {
	Matrix f, ret;
	f = Ma;
	ret.SetXY(Ma.x, Ma.y);
	ret.Clear();
	for (int i = 1; i <= ret.x; ++i) ret.A[i][i] = 1;
	while (b) {
		if (b & 1) ret = Mul(ret, f);
		b >>= 1;
		f = Mul(f, f);
	}
	return ret;
}

int main() 
{
	scanf("%lld%lld", &n, &Mod);
	LL Temp, Ans;
	M0.SetXY(1, 3);
	M0.Clear();
	M0.A[1][1] = 0; M0.A[1][2] = 1; M0.A[1][3] = 1;
	M_t.SetXY(3, 3);
	Temp = 1;
	for (int i = 1; i <= 18; ++i) {
		Temp *= 10ll;
		if (Temp > n) break;
		M_t.Clear();
		M_t.A[1][1] = Temp;
		M_t.A[2][1] = M_t.A[2][2] = M_t.A[3][2] = M_t.A[3][3] = 1;
		M_t = Pow(M_t, Temp - Temp / 10);
		M0 = Mul(M0, M_t);	
	}
	M_t.Clear();
	M_t.A[1][1] = Temp;
	M_t.A[2][1] = M_t.A[2][2] = M_t.A[3][2] = M_t.A[3][3] = 1;
	M_t = Pow(M_t, n - Temp / 10 + 1);
	M0 = Mul(M0, M_t);
	Ans = M0.A[1][1];
	printf("%lld\n", Ans);
	return 0;
}