[BZOJ 5221] [Lydsy2017省队十连测]偏题

BZOJ传送门

题目描述

斐波那契数列是一个经典递推数列，即 $Fi{b}_n=Fi{b}_{n-1}+Fi{b}_{n-2}$

在这个问题中，定义一个新数列，对于$n\ge 2$有

$F_n=F_{n-1}+F_{n-2}+2\sqrt{3+F_{n-1}{F}_{n-2}}$

对于给出的四个整数$F_0,F_1,M,n$，求$F_{n}modM$。

数据保证$\sqrt{3+F_{n-1}{F}_{n-2}}$总是整数。

输入输出格式

输入格式

一行四个整数$F_0,F_1,M,n$。

输出格式

一行一个整数，表示$F_n$对$M$取模后的值。

输入输出样例

输入样例#1：

1 1 10 5

输出样例#1：

4

【样例解释1】

$F_0=1$

$F_1=1$

$F_2=1+1+2\sqrt{3+1\times 1}=6$

$F_3=6+1+2\sqrt{3+6\times 1}=13$

$F_4=13+6+2\sqrt{3+13\times 6}=37$

$F_5=37+13+2\sqrt{3+37\times 13}=94$

输入样例#2

2 3 100 6

输出样例#2

82

数据规模

对于30%的数据，$n\le 20$。对于60%的数据，$F_0=F_1=1$。对于80%的数据，$n\le 10^5$。对于100%的数据，$0\le n\le 10_9，1\le M\le 10_9，1\le F_0\le F_1\le 10_6$。

解题分析

我们可以注意到：
$$\begin{gathered} F_n=F_{n-1}+F_{n-2}+2\sqrt{3+F_{n-1}{F}_{n-2}} \\ 3+F_n{F}_{n-1}=F_{n-1}^2+F_{n-2}{F}_{n-1}+2F_{n-1}\sqrt{3+F_{n-1}{F}_{n-2}}+3 \\ \sqrt{3+F_n{F}_{n-1}}=\sqrt{3+F_{n-1}{F}_{n-2}}+F_{n-1} \end{gathered}$$

所以我们可以用一个$3\times 3$的矩阵， 分别代表$\sqrt{3+F_n{F}_{n-1}},F_n,F_n-1$， 就可以矩阵加速递推了。

总复杂度$O(log(n))$。

代码如下：

#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cmath>
#include <algorithm>
#include <cctype>
#define R register
#define IN inline
#define W while
#define gc getchar()
#define ll long long
int mod, n, a, b;
struct Matrix {int mat[3][3];};
IN Matrix operator * (const Matrix &x, const Matrix &y)
{
	Matrix ret; std::memset(ret.mat, 0, sizeof(ret.mat));
	for (R int i = 0; i < 3; ++i)
	for (R int j = 0; j < 3; ++j)
	for (R int k = 0; k < 3; ++k)
	ret.mat[i][j] = (ret.mat[i][j] + 1ll * x.mat[i][k] * y.mat[k][j] % mod) % mod;
	return ret;
}
Matrix unit, res, base, ans;
IN Matrix qpow(R int tim)
{
	W (tim)
	{
		if(tim & 1) base = base * res;
		res = res * res, tim >>= 1;
	}
	return base;
}
int main(void)
{
	scanf("%d%d%d%d", &a, &b, &mod, &n);
	if(n == 0) return printf("%d", a), 0;
	if(n == 1) return printf("%d", b), 0;
	n--;
	unit.mat[0][0] = b, unit.mat[1][0] = a, unit.mat[2][0] = (ll)std::sqrt(3 + 1ll * a * b) % mod;
	res.mat[0][0] = 1, res.mat[0][1] = 1, res.mat[0][2] = 2;
	res.mat[1][0] = 1;
	res.mat[2][0] = 1, res.mat[2][2] = 1;
	base.mat[0][0] = base.mat[1][1] = base.mat[2][2] = 1;
	ans = qpow(n);
	ans = ans * unit;
	printf("%d", ans.mat[0][0]);
}