[ACNOI2022]数学压轴

题目

题目背景
“卷爷，去看下 $T3$ 吧。 ”

“哼，我叫你看 $T3$ 的时候你怎么不去看？”

$\mathsf{OneInDark}$ 啼笑皆非。“因为你是最强的啊，你要有最强者的担当！”

$\mathsf{OL}\mathsf{IYE}$ 睁开了双眼，印着直巴的纹，蕴着血红的底。获得这双眼睛的人，都曾经历过杀戮。

“我的想法很简单：谁敢润我，我就把他卷飞——这样，你们应该能够稍微体会我的孤独与愤怒了吧！”

题目描述
递推数列 $f_n$ 由如下方程给出：$f_0=f_1=1,f_n=9f_{n-1}+12f_{n-2}(n⩾2)$ 。

现有 $T$ 次询问，每次给定 $a,b,c,d,P$，请求出 $\gcd (a{f}_n+b{f}_{n+1},c{f}_n+d{f}_{n+1})\mathrm{mod}P$ 。

数据范围与约定
数据组数 $T⩽10^5$，而 $n⩽10^9$ 。还保证 $P⩽10^9+7$ 以及 $a,b,c,d\in [0,10^3]$ 。

思路

这 $100\%$ 的数学题啊。只怪我数学太差。

首先我们至少需要会解决 $\gcd (f_n,f_{n+1})$ 吧。然而就在这一步，我便停滞不前了！更可悲可叹的是，强如 $\mathsf{OL}\mathsf{IYE}$ 者早已发现，其有固定结论
$${\rho }_n=\gcd (f_n,f_{n+1})=3^{\lfloor \frac{n}{2}\rfloor }$$

证明则可使用归纳法。在此之前，先说明

引理：$f_n$ 中 $3$ 的次数为 ${\kappa }_n=\lfloor \frac{n}{2}\rfloor$，即 $3^{{\kappa }_n}\mid f_n$ 但 $3^{{\kappa }_n+1}\nmid f_n$ 。

证明：归纳法。首先 $3^{{\kappa }_n}\mid f_n$ 显然，因为系数 $9,12$ 都是 $3$ 的倍数。而 $f_n\equiv 12f_{n-2}(\mathrm{mod}{3}^{{\kappa }_n+1})$，因为 ${\kappa }_{n-1}+2⩾{\kappa }_n+1$ 。由归纳法 $12f_{n-2}$ 中 $3$ 的次数为 ${\kappa }_{n-2}+1={\kappa }_n$，所以 $3^{{\kappa }_n+1}\nmid f_n$ 得证。

而 $2\nmid f_n$ 是显然的。于是 ${\rho }_n=\gcd (f_n,12f_{n-1})=\gcd (f_n,3f_{n-1})=3^{[{\kappa }_n>{\kappa }_{n-1}]}\cdot {\rho }_{n-1}=3^{{\kappa }_n}$，得证。

总算搞定了这个超级弱化版问题。考虑原问题，通过辗转相减总可以使得 $d=0$ 。所以我们实际上只需解决 $\gcd (a{f}_n+b{f}_{n+1},c{f}_n)$ 。

考虑到这个答案可能很大，想要求出，必然是半道上取模。也就是说将该 $\gcd$ 表为若干乘积。可以看到 $c$ 孤零零的，那么我们可以把 $c$ 提出来考虑。

要么先求 $\gcd (\dots ,c)$，然后 $\dots$ 部分除以该值，再与 $f_n$ 取 $\gcd$ 。要么先与 $f_n$ 取 $\gcd$，剩余部分直接与 $c$ 求 $\gcd$ 。前者对问题的简化的帮助不大；选择后者，先求出
$$g=\gcd (a{f}_n+b{f}_{n+1},f_n)=\gcd (b{f}_{n+1},f_n)={\rho }_n\gcd \left(b,\frac{f_n}{{\rho }_n}\right)$$

最后一步的思想，跟目前求 $\gcd (\cdots ,c{f}_n)$ 一样：先计算 $f$，就会剩下一个含常数的 $\gcd$ 式。算 $\frac{f_n}{{\rho }_n}\mathrm{mod}b$ 可以用矩阵加速，过程中将 $3$ 除掉即可，因为转移矩阵的平方满足每一项都是 $3$ 的倍数。

然后再求
$$\gcd \left(\frac{a{f}_n+b{f}_{n+1}}{g},c\right)$$

这个看上去没法搞了。巧妙的是，${\rho }_n\mid g$，于是左式就等于
$$\left(a\frac{f_n}{{\rho }_n}+b\frac{f_{n+1}}{{\rho }_n}\right){\left[\gcd \left(b,\frac{f_n}{{\rho }_n}\right)\right]}^{-1}$$

然后对 $c$ 取模，分子就只需要对 $c\gcd (b,\frac{f_n}{{\rho }_n})$ 取模。这是可接受的了。计算方法不变，时间复杂度 $\mathcal{O}(T\log n)$ 。

代码

#include <cstdio> // megalomaniac JZM yydJUNK!!!
#include <iostream> // Almighty XJX yyds!!!
#include <algorithm> // decent XYX yydLONELY!!!
#include <cstring> // Pet DDG yydOLDGOD!!!
#include <cctype> // oracle: ZXY yydBUSlut!!!
typedef long long llong;
# define rep(i,a,b) for(int i=(a); i<=(b); ++i)
# define drep(i,a,b) for(int i=(a); i>=(b); --i)
# define rep0(i,a,b) for(int i=(a); i!=(b); ++i)
inline int readint(){
	int a = 0, c = getchar(), f = 1;
	for(; !isdigit(c); c=getchar()) if(c == '-') f = -f;
	for(; isdigit(c); c=getchar()) a = a*10+(c^48);
	return a*f;
}
inline int getGcd(int a, int b){
	for(; b; a%=b,a^=b^=a^=b){}; return a;
}

struct Matrix{ int a[2][2]; };
Matrix mat_mul(const Matrix &a, const Matrix &b, const int &mod){
	Matrix c; rep(i,0,1) rep(j,0,1) // use for loop
	c.a[i][j] = int((llong(a.a[i][0])*b.a[0][j]
		+llong(a.a[i][1])*b.a[1][j])%mod);
	return c; // hope that it's fast enough
}
Matrix qkpow(Matrix b, int q, const int &mod){
	if(!q) return Matrix{1,0,0,1};
	Matrix a = b; -- q; // avoid unit Matrix
	for(; q; q>>=1,b=mat_mul(b,b,mod))
		if(q&1) a = mat_mul(a,b,mod);
	return a;
}
inline llong qkpow(llong b, int q, const int &mod){
	llong a = 1;
	for(; q; q>>=1,b=b*b%mod) if(q&1) a = a*b%mod;
	return a;
}

Matrix norm; int coe[2][2];
int main(){
	norm.a[0][1] = 1, norm.a[1][0] = 12, norm.a[1][1] = 9;
	Matrix trans = mat_mul(norm,norm,1e9+7);
	rep(i,0,1) rep(j,0,1) trans.a[i][j] /= 3;
	for(int T=readint(); T; --T){
		int n = readint(), mod = readint();
		rep(i,0,1) rep(j,0,1) coe[i][j] = readint();
		for(int f; coe[1][1]; std::swap(coe[0],coe[1])){
			f = coe[0][1]/coe[1][1]; // ratios
			rep(i,0,1) coe[0][i] -= f*coe[1][i];
		}
		if(!coe[1][0]){ // just the other one
			Matrix ans = qkpow(norm,n,mod);
			int fn = ans.a[0][0]+ans.a[0][1];
			int fn_1 = ans.a[1][0]+ans.a[1][1];
			printf("%lld\n",(llong(fn_1)*coe[0][1]
				+llong(fn)*coe[0][0])%mod); continue;
		}
		if(coe[1][0] < 0) coe[1][0] = -coe[1][0];
		coe[0][0] = coe[0][0]%coe[1][0]+coe[1][0];
		if(coe[0][1] == 0){
			Matrix ans = qkpow(norm,n,mod);
			llong fn = ans.a[0][0]+ans.a[0][1]; // get this
			printf("%lld\n",getGcd(coe[0][0],
				coe[1][0])*fn%mod); continue;
		}
		Matrix tmp = qkpow(trans,n>>1,coe[0][1]); // [2m,2m+1]
		int bgcd = tmp.a[n&1][0]+tmp.a[n&1][1];
		bgcd = getGcd(bgcd,coe[0][1]); // part of g
		int cmod = coe[1][0]*bgcd; // modules
		tmp = qkpow(trans,n>>1,cmod); // [2m,2m+1]
		int frhon = tmp.a[n&1][0]+tmp.a[n&1][1], frhon_1;
		if(n&1) frhon_1 = int((llong(9)*frhon+
			llong(12)*(tmp.a[0][0]+tmp.a[0][1]))%cmod);
		else frhon_1 = tmp.a[1][0]+tmp.a[1][1];
		int cgcd = int((llong(coe[0][0])*frhon
			+llong(coe[0][1])*frhon_1)%cmod);
		cgcd = getGcd(cgcd/bgcd,coe[1][0]);
		printf("%lld\n",qkpow(3,n>>1,mod)*bgcd%mod*cgcd%mod);
	}
	return 0;
}

后记

我本来想构建模 $p^k$ 的域，但我失败了。形式幂级数？形式洛朗级数？好像都不行。拜托，如果这种东西真的存在，那肯定早就被提出与使用了；用得着我吗！

我去问数竞同学 $\mathsf{YJ}$，他表示他们不会研究这种东西，只要 “能在有限步内结束” 就 $\mathrm{OK}$ 了。也就是说他也没看出来三的次幂的结论。