(Luogu Solution)SP5152:BFALG - Brute-force Algorithm EXTREME

• 【题目链接】

Link:SP5152

• 【前置知识】

  • 矩阵快速幂，矩阵加速

  • 扩展欧拉定理

• 【解题思路】

看到这个数据范围，我们应该能想到，这是一道结论题。

再看一下SPOJ几页的 $0.00s$，我们应该会觉得这题正解 $\Theta (1)$。

那么，我们推一下结论吧。

首先，我们设 $f(n)=1$，那么计算 $F_{n,f}$ 中调用 $f$ 的次数即为 $F_{n,f}$ 的返回值。（这里 $F$ 请用题目翻译中的定义。）

然后，我们又可以知道， $F_{a,F_{b,f}}=F_{a,f}{F}_{b,f}$。因为，显然，$F_{a,F_{b,f}}$ 会把 $F_{b,f}$ 调用 $F_{a,f}$ 次。

看起来也没有什么思路，我最开始用的是暴力展开。然后，有一个问题，就是不可避免的要计算任意 $n$ 对应的 $F_{n,f}$ 的值。于是又绕回了原来的问题。不过，我们可以试着，如果把所有的 $F_{n,f}$ 的值列出来，会怎样呢？

接着，我们就来从 $F_1,f$ 开始写出它对应的值：

$F_{1,f}=a,F_{2,f}=b,F_{3,f}=F_{2,F_{1,f}}=ab,F_{4,f}=F_{3,F_{2,f}}=a{b}^2,F_{5,f}=F_{4,F_{3,f}}=a^2{b}^3\cdots$

发现了吗？$a$ 与 $b$ 的指数是斐波那契数列！

其实这个性质的证明也很显然，因为把前两项指数相加就是这一项的指数，而这正是斐波那契数列的定义。

那么，我们对斐波那契数列的边界条件做点微调：设数列 $g$，它满足以下条件：

$$g_i=\{\begin{array}{ll}1& i=1\\ 0& i=2\\ g_{i-1}+g_{i-2}& i>2\end{array}$$

则答案即为 $a^{g_n}{b}^{g_{n+1}}\mathrm{m}\mathrm{o}\mathrm{d}p$。

而我们发现，$g_n$ 与 $g_{n+1}$ 可以通过矩阵加速 $\Theta ({\log }_{2}n)$ 求出！

那么，接下来就是一个快速幂，没了。复杂度 $\Theta ({\log }_{2}n)$。

没了？真没了？

好好再想一想，我们该对指数怎么处理？如果像平常那样把指数对 $p$ 取模，你会发现样例都过不了。我们总不能用高精算一遍再取模吧？

你可能想，对 $p-1$ 取模不就行了？然后，你会发现，题目中没有说过 $p$ 是质数，所以费马小定理不能用。 而且，最阴间的是，题目同样没有对 $a,b$ 作额外限定。

那么，还有一种方法：扩展欧拉定理。它不需要 $a,b$ 与 $p$ 互质。

到这里，这个问题就被真正的解决了。

也许你可能会说：开始不是说要 $\Theta (1)$ 吗？

是的。但是，这个 $O(\sqrt{p}{\log }_{2}n)$ 算法确实是本题的正解。之所以我会说是 $\Theta (1)$ 的结论题，是因为这个猜想会引导我们往这个方向去想，从而导出真正的正解。

• 【实现细节】

问题虽然被解决了，代码还是要注意的。

扩展欧拉定理怎么实现，就看你自己了。我的做法是在矩阵乘法函数里判定是否大于 $\phi (p)$。

还有一种极为阴间的情况：$p=1$！注意这一点，特判好了，就能过了。

对了，我曾经特判输出 $0$ 的时候忘记输出测试信息了，调了好久……

• 【代码实现】

#include <iostream>
#include <vector>

using namespace std;

typedef unsigned long long lrg_l;
typedef vector<vector<lrg_l>> matrix;

lrg_l moder,prime;
matrix greater_varphi;

long long fast_pow(long long base,long long exp)
{
    long long result;
    for(result=1;exp;exp&1?result=result*base%prime:true,base=base*base%prime,exp>>=1);
    return result;
}

matrix operator*(matrix& src,matrix& param)
{
	const int len=src.size();
	matrix result(len,vector<lrg_l>(len,0));
	for(int k=0;k<len;k++)
		for(int i=0;i<len;i++)
			for(int j=0;j<len;j++)
			{
				result[i][j]=result[i][j]+src[i][k]*param[k][j];
				if(result[i][j]>moder)
					result[i][j]%=moder,greater_varphi[i][j]=true;
			}
	return result;
}

int main(int argc,char *argv[],char *envp[])
{
	ios::sync_with_stdio(false);
	cin.tie(nullptr);
	long long testcnt,power,base_A,base_B,tmp;
	cin>>testcnt;
	matrix result(2,vector<lrg_l>(2)),base(2,vector<lrg_l>(2)),
		initstate={{1,1},{1,0}},resstate={{0,1},{0,0}},greater_varphi_init={{false,false},{0,0}};
	for(int t=1;t<=testcnt;t++)
	{
		cin>>base_A>>base_B>>prime>>power;
		power--;
		if(prime==1)
		{
			cout<<"Case #"<<t<<": 0\n";
			continue;
		}
		base=initstate;
		result=resstate;
		greater_varphi=greater_varphi_init;
		moder=prime;
		tmp=prime;
		for(int i=2;i*i<=prime;i++)
			if(prime%i==0)
			{
				moder=moder/i*(i-1);
				while(prime%i==0)
					prime/=i;
			}
		if(prime>1)
			moder=moder/prime*(prime-1);
		prime=tmp;
		while(power)
		{
			if(power&1)
				result=result*base;
			base=base*base;
			power>>=1;
		}
		cout<<"Case #"<<t<<": "<<fast_pow(base_A,result[0][1]+greater_varphi[0][1]*moder)*fast_pow(base_B,result[0][0]+greater_varphi[0][0]*moder)%prime<<'\n';
	}
 	return 0;
}