CF1106F Lunar New Year and a Recursive Sequence 矩阵加速+BSGS+原根+扩展欧几里得

最新推荐文章于 2021-07-14 21:19:11 发布

PI_PJW

最新推荐文章于 2021-07-14 21:19:11 发布

阅读量254

点赞数

分类专栏：数学文章标签：矩阵加速 BSGS 原根扩展欧几里得

本文链接：https://blog.csdn.net/weixin_43908980/article/details/108932054

版权

数学专栏收录该内容

26 篇文章 2 订阅

订阅专栏

文章目录

一.题目
二.Solution
三.Code
Thanks!

一.题目

有一串 $n(n\leqslant10^9)$ 个数的数列，给你 $b_1\sim b_k(k\leqslant100)$ 当 $i > k i > k$ 时： $f_i=(\prod\limits_{j=1}^kf_{i-j}^{b_j})\bmod{998244353}$ 已知 $f_1=f_2=\cdots=f_{k-1}=1,f_n=m$ ，问最小的正整数 $f_k$ 可能是多少?

二.Solution

这是一道反正我考场上推不出来的的数学题。
首先因为 $f_1\sim f_{k-1}=1$ 并且递推式是累乘，所以 $f_n$ 必定可以表示成 $f_n=f_k^p$ 的形式，然后 $p$ 又满足递推式： $g[i]=\prod^{k}_{j=1}b[j]*g[i-j]$ 所以可以用矩阵加速来计算 $p$ 的大小，构造这样一个递推矩阵 $B$ ：
$\begin{matrix} 0 & 0 & ... & 0 & b[k] \\ 1 & 0 & ... & 0 & b[k - 1] \\ 0 & 1 & ... & 0 & b[k - 2] \\ ... &... &... &...\\ 0 & 0 & ... & 1 & b[1] \end{matrix}$
这样一个答案矩阵 $A$ ：
$\begin{matrix} 0 & 0 & ... & 1 \end{matrix}$
最后快速幂完成后 $p$ 就是 $A . j z [1] [k]$
那么问题就来到 $f_{k}^{p}\equiv m (mod\,998244353)$ 求 $f_k$ 。
这里可以用原根做，我们知道 $998244353$ 的原根是 $3$ ，那么设 $f_k\equiv 3^s(mod\,998244353)$ ，设 $m\equiv 3^t(mod\,998244353)$ ，因为 $m$ 已知，所以可以用BSGS求出 $t$ ，最后又化成： $3^{sp}\equiv 3^t (mod\,998244353)$ 又变成 $sp\equiv t(mod\,998244352)$
这里用扩展欧几里得定理求出 $s$ 即可，答案就是 $3^{sp}$

三.Code

#include <cstdio>
#include <cstring>
#include <iostream>
#include <map>
#include <cmath>
using namespace std;

#define M 105
#define LL long long
const LL mod = 998244353;

int k, n;
LL b[M], m;
map <LL, LL> mp;
struct matrix{
	LL jz[M][M];
	matrix clear (int flag){
		if (! flag){
			memset (jz, 0, sizeof jz);
		}
		else{
			for (int i = 1; i <= k; i ++)
				for (int j = 1; j <= k; j ++)
					jz[i][j] = (i == j) ? 1 : 0;
		}
	}
	matrix operator * (const matrix& rhs) const{
		matrix fina;
		fina.clear(0);
		for (int i = 1; i <= k; i ++)
			for (int j = 1; j <= k; j ++)
				for (int z = 1; z <= k; z ++)
					fina.jz[i][j] = (fina.jz[i][j] + jz[i][z] * rhs.jz[z][j] % (mod - 1)) % (mod - 1);
		return fina;
	}
}A, B;

matrix qkpow1 (matrix x, LL y){
	matrix fina;
	fina.clear(1);
	while (y){
		if (y & 1)
			fina = fina * x;
		x = x * x;
		y >>= 1;
	}
	return fina;
}
LL qkpow2 (LL x, LL y){
	LL res = 1;
	while (y){
		if (y & 1)
			res = res * x % mod;
		x = x * x % mod;
		y >>= 1;
	}
	return res;
}
LL bsgs (LL a, LL b){
	mp.clear();
	mp[1] = 0;
	LL sqt = ceil (sqrt (mod - 1)), t1 = 1, t2 = 1;
	for (int i = 1; i < sqt; i ++){
		t1 = t1 * a % mod;
		mp[t1 * b % mod] = 1ll * i;
	}
	t1 = t1 * a % mod;
	for (int i = 1; i <= sqt; i ++){
		t2 = t1 * t2 % mod;
		if (mp[t2])
			return (1ll * i * sqt - mp[t2]);
	}
}
LL get_gcd (LL x, LL y){
	if (! y)
		return x;
	return get_gcd (y, x % y);
}
void exgcd (LL &x, LL &y, LL a, LL b){
	if (! b){
		x = 1, y = 0;
		return ;
	}
	exgcd (y, x, b, a % b);
	y -= a / b * x;
}
int main (){
	scanf ("%d", &k);
	for (int i = 1; i <= k; i ++)
		scanf ("%lld", &b[i]);
	scanf ("%d %lld", &n, &m);
	A.jz[1][k] = 1;
	for (int i = 1; i <= k; i ++)
		B.jz[i][k] = b[k - i + 1];
	for (int i = 1; i < k; i ++)
		B.jz[i + 1][i] = 1;
	A = A * qkpow1 (B, n - k);
	LL p = A.jz[1][k], t = bsgs (3, m);
	if (t % get_gcd (p, mod - 1))
		printf ("-1\n");
	else{
		LL gcd = get_gcd (p, mod - 1);
		LL x, y, b = mod - 1;
		p /= gcd, t /= gcd, b /= gcd;
		exgcd (x, y, p, b);
		x = (x % (mod - 1) + mod - 1) % (mod - 1);
		x = x * t % (mod - 1);
		printf ("%lld\n", qkpow2 (3, x));
	}
	return 0;
}