其他分享
首页 > 其他分享> > Luogu P5110 块速递推

Luogu P5110 块速递推

作者:互联网

传送门

又是一道常系数线性递推的题,一下子便能想到矩阵快速幂,但一看数据规模

\[1≤T≤5×10^7 \]

我们可以知道,普通的矩阵快速幂肯定是过不了的

显然,本题只能接受\(O(T)\)时间复杂度的方法(LYT:?)

考虑光速幂

对于\(x^y=(x^b)^{\lfloor\frac{y}{b}\rfloor}x^{y\%b}\),,我们预处理出\(x^n(n<b)\),再预处理出\((x^b)^n\),令\(b=a^k(k\in (0,1))\),则预处理时空复杂度为\(O(b+\frac{y}{b})\),单次查询\(O(1)\)。

更一般的,\(x^y=x^{a_0}(x^b)^{a_1}(x^{b^2})^{a_2}\cdots\Rightarrow y=a_0+a_1b+a_2b^2+a_3b^3\cdots\),就是将\(y\)进行\(b\)进制分解,对每一位进行预处理,令\(b=a^k(k\in (0,1))\) 预处理空间复杂度为\(O(b=a^k)\),查询时间复杂度为\(O(\frac{1}{k})\).

1.矩阵光速幂

但我们可以发现,因为\(a_i< 2^{64}\),所以\(k=0.25\),对于\(5*10^7\),哪怕处理了光速幂,\(3\)次矩阵乘法(约\(30\)的常数)也是难以通过的

2.特征多项式?

我不会,这里推荐一个大佬的博客:【学习笔记】特征多项式

但\(3\)次多项式乘法(约\(15\)的常数)依然难以通过

3.使用生成函数解出该式通解

先推导一些生成函数对应的函数

\[\frac{1}{1-ax}=g\\ 1=g-axg\\ axg+1=g\\ g_0=1,g_{n+1}=ag_n\\ \]

\[\frac{1}{a-x}=g\\ 1=ag-gx\\ g=\frac{1+gx}{a}\\ g_0=\frac{1}{a},g_{n+1}=\frac{g_n}{a} \]

进入正题

\[a_n=233a_{n-1}+666a_{n-2},a_0=0,a_1=1\\ 我们设函数g为序列a的生成函数\\ 则有:\\ \begin{aligned} &233xg+666x^2g+x=g\\ \Leftrightarrow &g=\frac{x}{1-233x-666x^2}\\ \Leftrightarrow &g=-x\left(\frac{1}{666x^2+233x-1}\right)\\ \Leftrightarrow &-\frac{g}{x}=\frac{1}{666(x-\frac{-233+\sqrt{56953}}{1332})(x-\frac{-233-\sqrt{56953}}{1332})}\\ \Leftrightarrow &-\frac{666g}{x}=\frac{1}{(x-\frac{-233+\sqrt{56953}}{1332})(x-\frac{-233-\sqrt{56953}}{1332})}\\ &-\frac{666g}{x}=\frac{\frac{1}{x-\frac{-233+\sqrt{56953}}{1332}}-\frac{1}{x-\frac{-233-\sqrt{56953}}{1332}}}{\frac{\sqrt{56953}}{666}}\\ &\frac{g}{x}=\frac{\frac{1}{\frac{-233+\sqrt{56953}}{1332}-x}-\frac{1}{\frac{-233-\sqrt{56953}}{1332}-x}}{\sqrt{56953}}\\ &a_n=\frac{\left(\frac{233+\sqrt{56953}}{2}\right)^n+\left(\frac{233-\sqrt{56953}}{2}\right)^n}{\sqrt{56953}} \end{aligned} \]

这里补充一下最后一步的推导

\[令a=\frac{-233+\sqrt{56953}}{1332}\\ \Rightarrow \frac{1}{a}=\frac{233+\sqrt{56953}}{2}\\ \frac{1}{\frac{-233+\sqrt{56953}}{1332}-x}\\ =\sum_{i=0}^{+\infty}\frac{x^i}{a^{i+1}} \]

到了这一步,我们就可以通过\(2\)种方式来解决

Ⅰ.扩域

可以看出,上式中出现了无理数,但在模意义下无理数是无意义的,因为模运算是定义在\(L=Z\cap [0,mod=1000000007)\)上的。(如无特殊说明,下文中四则运算皆为模意义下)

我们不妨定义域\(K=(\{a+b\sqrt{56953}|a,b\in L\},+,*,-)\)

证明\(K\)为\(L\)的扩域

\[(a+b\sqrt{56953})(c+d\sqrt{56953})=(ac+56953bd)+(ad+bc)\sqrt{56953}\in K \]

之后就可以像复数运算一样操作了

#include<bits/stdc++.h>
using namespace std;
# define ll long long
# define read read1<ll>()
# define Type template<typename T>
# define fre(k) freopen(k".in","r",stdin);freopen(k".out","w",stdout)
Type T read1(){
	T t=0;
	char k;
	bool vis=0;
	do (k=getchar())=='-'&&(vis=1);while('0'>k||k>'9');
	while('0'<=k&&k<='9')t=(t<<3)+(t<<1)+(k^'0'),k=getchar();
	return vis?-t:t;
}
namespace Mker
{
	unsigned long long SA,SB,SC;
	void init(){scanf("%llu%llu%llu",&SA,&SB,&SC);}
	unsigned long long rand()
	{
	    SA^=SA<<32,SA^=SA>>13,SA^=SA<<1;
	    unsigned long long t=SA;
		SA=SB,SB=SC,SC^=t^SA;return SC;
	}
}
# define mod 1000000007
# define M 32767
# define N 32768
# define inv2 500000004ll
struct A{
	int x,y;
	A& operator *=(const A &b){
		int t=y;
		y=(1ll*x*b.y+1ll*y*b.x)%mod;
		x=(1ll*x*b.x+56953ll*t%mod*b.y)%mod;
		return *this;
	}
	A operator *(const A &b){A t=*this;return t*=b;}
}f[2][N|1],f1[2][N|1],o,p;
int main(){
	int T=read,ans=0;
	Mker::init();o.x=233*inv2%mod,o.y=inv2;p.x=o.x,p.y=mod-inv2;
	f[0][0].x=f1[0][0].x=1;f[1][0].x=f1[1][0].x=1;
	for(int i=1;i<=N;++i)f[0][i]=f[0][i-1]*o,f[1][i]=f[1][i-1]*p;
	for(int i=1;i<=N;++i)f1[0][i]=f1[0][i-1]*f[0][N],f1[1][i]=f1[1][i-1]*f[1][N];
	while(T--){
		unsigned ll n=Mker::rand()%(mod-1);
		int w=n>>15&M,t=(1ll*f[0][n&M].x*f1[0][w].y+1ll*f[0][n&M].y*f1[0][w].x-1ll*f[1][n&M].x*f1[1][w].y-1ll*f[1][n&M].y*f1[1][w].x)%mod;
		if(t<0)t+=mod;ans^=t;
	}cout<<ans;
	return 0;
}

Ⅱ.二次剩余

我们经过二次剩余可以知道,对于素数\(mod\),我们可以找到唯一的\(x\)满足\(x^2\equiv 56953 \mod mod\).

经过打表可知,对于本题,\(188305837 \equiv 56953 \mod mod\)

(可能你想知道

Mker::rand()%(mod-1)

是为什么

通过\(扩域\)我们可以知道域\(K\)是满足域\(L\)的运算规则的,自然而然满足费马定理,所以在模意义下\(a_n与a_{n+k(mod-1)}\)相等

标签:frac,Luogu,速递,sqrt,56953,P5110,mod,1332,233
来源: https://www.cnblogs.com/SYDevil/p/14180205.html