[NOI2011]兔农 矩阵乘法

Description
你知道普通的Fibonacci数列吗，他的递推式为：
f[i]=f[i-1]+f[i-2]
但这题有点不一样。
这题给出一个k，如果当前f[i]%k==1，f[i]–。
其他还是照样递推，让你输出f[n]%mod。

---

Sample Input
6 7 100

---

Sample Output
7

---

这题很容易发现他是有循环节的，若无循环节那么减1的f数量一定很有限，剩下的做矩阵乘法快速幂即可。
如何找出循环节呢，我使用的是exgcd求逆元的方式，因为假设你遇到一个i将他减1变成0，后面的第j项其实就是：fib[j]*f[i-1]。（其中f[i]表示在即成数列中的权值，fib[i]表示在原Fibonacci数列中的第i项的值，那么其实你就是要找一个fib[j]*f[i-1]%k==1，你先预处理出每个值的位置，解一个exgcd即可。
然后中间是要用矩阵乘法加速的，矩阵有两种：
一种是普通Fibonacci矩阵的递推矩阵，
还有一种是减一的矩阵，自己推吧。。。
又臭又长的代码

---

#include <cstdio>
#include <cstring>

using namespace std;
typedef long long LL;

LL mod;
struct matrix {
	LL a[3][3];
	matrix() {memset(a, 0, sizeof(a));}
	friend matrix operator * (matrix a, matrix b) {
		matrix c;
		for(int i = 0; i < 3; i++) {
			for(int j = 0; j < 3; j++) {
				for(int k = 0; k < 3; k++) {
					(c.a[i][j] += a.a[i][k] * b.a[k][j] % mod) %= mod;
				}
			}
		} return c;
	}
} ans, A, h1, h2;
int p[1100000];
LL pos[1100000], f[3100000];
int v[1100000];

LL exgcd(LL a, LL b, LL &x, LL &y) {
	if(b == 0) {x = 1, y = 0; return a;}
	else {
		LL tx, ty, d = exgcd(b, a % b, tx, ty);
		x = ty, y = tx - ty * (a / b);
		return d;
	}
}

int main() {
	LL n; scanf("%lld", &n);
	LL k; scanf("%lld%lld", &k, &mod);
	LL now = 1, hh = 0; bool bk = 0;
	f[1] = f[2] = 1;
	for(int i = 3; ; i++) {
		f[i] = (f[i - 1] + f[i - 2]) % k;
		if(!pos[f[i]]) pos[f[i]] = i;
		if(f[i] == 1 && f[i - 1] == 1) break;
	}
	h1.a[0][1] = 1, h1.a[0][0] = h1.a[0][2] = 0;
	h1.a[1][0] = h1.a[1][1] = 1, h1.a[1][2] = 0;
	h1.a[2][0] = h1.a[2][1] = 0, h1.a[2][2] = 1;
	h2.a[0][0] = 1, h2.a[0][1] = h2.a[0][2] = 0;
	h2.a[1][0] = h2.a[1][2] = 0, h2.a[1][1] = 1;
	h2.a[2][0] = h2.a[2][1] = -1, h2.a[2][2] = 1;
	A.a[0][0] = A.a[1][1] = A.a[2][2] = 1;
	LL uy, ux = 1;
	for(int i = 1; ; i++) {
		LL x, y, d = exgcd(now, k, x, y);
		if(d > 1) {bk = 1; uy = i - 1; break;}
		x = (x % k + k) % k;
		if(v[x]) {ux = v[x], uy = i - 1; break;}
		hh += pos[x]; v[x] = i;
		p[i] = hh;
		now = now * f[pos[x] - 1] % k;
	} ans.a[0][0] = 0; ans.a[0][1] = ans.a[0][2] = 1;
	for(int i = ux; i <= uy; i++) {
		int u = p[i] - p[i - 1];
		matrix o = h1;
		while(u) {
			if(u & 1) A = A * o;
			o = o * o; u /= 2;
		} A = A * h2;
	}
	if(bk) {
		if(n < hh) {
			LL f = n; LL u;
			for(int i = 1; i <= uy; i++) {
				if(p[i] > f) {
					u = f - p[i - 1];
					matrix o = h1;
					while(u) {
						if(u & 1) ans = ans * o;
						o = o * o; u /= 2;
					} break;
				}
				u = p[i] - p[i - 1];
				matrix o = h1;
				while(u) {
					if(u & 1) ans = ans * o;
					o = o * o; u /= 2;
				} ans = ans * h2;
				if(p[i] == f) break;
			}
		} else {
			n -= hh;
			ans = ans * A;
			while(n) {
				if(n & 1) ans = ans * h1;
				h1 = h1 * h1; n /= 2;
			}
		} printf("%lld\n", (ans.a[0][0] + mod) % mod);
	} else {
		if(n <= p[uy]) {
			LL f = n, u;
			for(int i = 1; i <= uy; i++) {
				if(p[i] > f) {
					u = f - p[i - 1];
					matrix o = h1;
					while(u) {
						if(u & 1) ans = ans * o;
						o = o * o; u /= 2;
					} break;
				}
				u = p[i] - p[i - 1];
				matrix o = h1;
				while(u) {
					if(u & 1) ans = ans * o;
					o = o * o; u /= 2;
				} ans = ans * h2;
				if(p[i] == f) break;
			} printf("%lld\n", (ans.a[0][0] + mod) % mod);
		} else {
			for(int i = 1; i < ux; i++) {
				int u = p[i] - p[i - 1];
				matrix o = h1;
				while(u) {
					if(u & 1) ans = ans * o;
					o = o * o; u /= 2;
				} ans = ans * h2;
			} LL gg = p[ux - 1]; n -= gg;
			for(int i = ux; i <= uy; i++) p[i - ux + 1] = p[i] - gg;
			uy = uy - ux + 1;
			LL f = n % p[uy];
			n /= p[uy];
			while(n) {
				if(n & 1) ans = ans * A;
				A = A * A; n /= 2;
			} int u;
			for(int i = 1; i <= uy; i++) {
				if(p[i] > f) {
					u = f - p[i - 1];
					matrix o = h1;
					while(u) {
						if(u & 1) ans = ans * o;
						o = o * o; u /= 2;
					} break;
				}
				u = p[i] - p[i - 1];
				matrix o = h1;
				while(u) {
					if(u & 1) ans = ans * o;
					o = o * o; u /= 2;
				} ans = ans * h2;
				if(p[i] == f) break;
			} printf("%lld\n", (ans.a[0][0] + mod) % mod);
		}
	}
	return 0;
}