ICode9

精准搜索请尝试: 精确搜索
首页 > 其他分享> 文章详细

Luogu P5110 块速递推

2020-12-23 18:32:19  阅读:237  来源: 互联网

标签:frac Luogu 速递 sqrt 56953 P5110 mod 1332 233


传送门

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

\[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

本站声明: 1. iCode9 技术分享网(下文简称本站)提供的所有内容,仅供技术学习、探讨和分享;
2. 关于本站的所有留言、评论、转载及引用,纯属内容发起人的个人观点,与本站观点和立场无关;
3. 关于本站的所有言论和文字,纯属内容发起人的个人观点,与本站观点和立场无关;
4. 本站文章均是网友提供,不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属;如您发现该文章侵犯了您的权益,可联系我们第一时间进行删除;
5. 本站为非盈利性的个人网站,所有内容不会用来进行牟利,也不会利用任何形式的广告来间接获益,纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。

专注分享技术,共同学习,共同进步。侵权联系[81616952@qq.com]

Copyright (C)ICode9.com, All Rights Reserved.

ICode9版权所有