[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;
}
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值