洛谷P5110:块速递推(特征根方程、光速幂)

解析

去你的搬砖生成函数,特征根太香了。
一开始我是用生成函数解的,和特征根相比有亿点点搬砖…
但是这个东西原理似乎使用一些神奇的等比差分,有些玄学,生成函数较易理解。

背下来背下来!

就以本题为情境讲一下特征根方程的解题流程吧。

对于递推式:
an=233an−1+666an−2a_n=233a_{n-1}+666a_{n-2}an=233an1+666an2
其特征方程就是:
x2=233x+666x^2=233x+666x2=233x+666
解得两个根:
p1=233+569532p_1=\frac{233+\sqrt{56953}}{2}p1=2233+56953
p2=233−569532p_2=\frac{233-\sqrt{56953}}{2}p2=223356953
那么通项公式就可以写成
an=c1p1n+c2p2na_n=c_1p_1^n+c2p_2^nan=c1p1n+c2p2n

当解得的 p1=p2=pp_1=p_2=pp1=p2=p 时,通项会变成:an=(c1+c2n)pna_n=(c_1+c_2n)p^nan=(c1+c2n)pn

然后带入 a0,a1a_0,a_1a0,a1 解出 c0,c1c_0,c_1c0,c1
0=c1+c20=c_1+c_20=c1+c2
1=c1p1+c2p21=c_1p_1+c_2p_21=c1p1+c2p2
解得
c1=156953c_1=\frac{1}{\sqrt{56953}}c1=569531
c2=−156953c_2=-\frac{1}{\sqrt{56953}}c2=569531
带回原式:
an=(233+569532)n−(233−569532)n56953a_n=\frac{(\dfrac{233+\sqrt{56953}}{2})^n-(\dfrac{233-\sqrt{56953}}{2})^n}{\sqrt{56953}}an=56953(2233+56953)n(223356953)n
然后有一个很重要的技巧:寻找 569535695356953 在对 1e9+71e9+71e9+7 的模意义下是否是二次剩余
出题人总会给人一条活路,写个暴力跑一会就能找到:
1883058372≡56953(mod109+7)188305837^2\equiv 56953\pmod{10^9+7}188305837256953(mod109+7)
所以式子就变成了:
an=(233+1883058372)n−(233−1883058372)n188305837a_n=\frac{(\dfrac{233+188305837}{2})^n-(\dfrac{233-188305837}{2})^n}{188305837}an=188305837(2233+188305837)n(2233188305837)n
化简得:
an=94153035n−905847205n188305837a_n=\frac{94153035^n-905847205^n}{188305837}an=18830583794153035n905847205n
再把逆元提前算出来:
an=233230706(94153035n−905847205n)a_n=233230706(94153035^n-905847205^n)an=233230706(94153035n905847205n)
其中的指数 nnn 可以用欧拉定理缩小到 O(mod)O(mod)O(mod) 范围。
然后问题就变成了如何快速求这个东西。

可以使用光速幂
原理很好理解,使用类似分块的思路,设 p=modp=\sqrt{mod}p=modO((mod))O(\sqrt(mod))O((mod)) 预处理所有 xkp(k≤⌊modp⌋)x^{kp}(k\le\lfloor \frac{mod}{p}\rfloor)xkp(kpmod)xk(k<p)x^k(k<p)xk(k<p),然后任意幂次都可以转化成 x⌊kp⌋∗p×xk%px^{\lfloor\frac{k}{p}\rfloor*p}\times x^{k\%p}xpkp×xk%p,然后实现 O(1)O(1)O(1) 计算。

代码

#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define ull unsigned long long
#define debug(...) fprintf(stderr,__VA_ARGS__)
inline ll read() {ll x(0),f(1);char c=getchar();while(!isdigit(c)) {if(c=='-')f=-1;c=getchar();}while(isdigit(c)) {x=(x<<1)+(x<<3)+c-'0';c=getchar();}return x*f;
}
const int N=4e5+100;
const int mod=1e9+7;
int n,m,k;
inline ll ksm(ll x,ll k){ll res=1;while(k){if(k&1) res=res*x%mod;x=x*x%mod;k>>=1;}return res;
}
int niv2=ksm(2,mod-2);
int r[N];
void init(int n,int &lim){lim=1;int L=0;while(lim<n) lim<<=1,L++;for(int i=1;i<lim;i++) r[i]=(r[i>>1]>>1)|((i&1)<<(L-1));
}
void NTT(ll *x,int lim,int op){for(int i=0;i<lim;i++) if(i<r[i]) swap(x[i],x[r[i]]);for(int l=1;l<lim;l<<=1){ll w=ksm(3,(mod-1)/(l<<1));if(op==-1) w=ksm(w,mod-2);for(int st=0;st<lim;st+=(l<<1)){for(ll i=0,now=1;i<l;i++,now=now*w%mod){ll u=x[st+i],v=now*x[st+i+l]%mod;x[st+i]=u+v>=mod?u+v-mod:u+v;x[st+i+l]=u-v<0?u-v+mod:u-v;}}}if(op==-1){ll ni=ksm(lim,mod-2);for(int i=0;i<lim;i++) x[i]=x[i]*ni%mod;}return;
}
void copy(ll *a,ll *b,int n,int lim){assert(n<=lim);memcpy(a,b,sizeof(ll)*n);fill(a+n,a+lim,0);return;
}
void mul(ll *a,ll *b,ll *c,int n,int m){static ll u[N],v[N];static int lim;//for(int i=0;i<n;i++) printf("%lld ",a[i]);putchar('\n');//for(int i=0;i<m;i++) printf("%lld ",b[i]);putchar('\n');init(n+m-1,lim);copy(u,a,n,lim);copy(v,b,m,lim);NTT(u,lim,1);NTT(v,lim,1);for(int i=0;i<lim;i++) c[i]=u[i]*v[i]%mod;NTT(c,lim,-1);//for(int i=0;i<n+m-1;i++) printf("%lld ",c[i]);putchar('\n');//putchar('\n');return;
}
void inv(ll *h,ll *f,int n){static ll t1[N],t2[N];static int lim;if(n==1){f[0]=ksm(h[0],mod-2);return;}inv(h,f,(n+1)>>1);init(n<<1,lim);fill(f+((n+1)>>1),f+lim,0);copy(t1,f,n,lim);copy(t2,h,n,lim);NTT(t1,lim,1);NTT(t2,lim,1);for(int i=0;i<lim;i++) t1[i]=(2*t1[i]-t1[i]*t1[i]%mod*t2[i]%mod+mod)%mod;NTT(t1,lim,-1);memcpy(f,t1,sizeof(ll)*n);return;
}
//499122177
void Sqrt(ll *h,ll *f,int n){static ll t1[N],t2[N];static int lim;if(n==1){f[0]=1;return;}Sqrt(h,f,(n+1)>>1);init(n<<1,lim);fill(f+((n+1)>>1),f+lim,0);inv(f,t1,n);fill(t1+n,t1+lim,0);mul(h,t1,t1,n,n);copy(t2,f,n,lim);NTT(t1,lim,1);NTT(t2,lim,1);for(int i=0;i<lim;i++) t1[i]=(t1[i]+t2[i])*niv2%mod;NTT(t1,lim,-1);memcpy(f,t1,sizeof(ll)*n);return;
}
void dao(ll *h,ll *f,int n){static ll t[N];static int lim;init(n<<1,lim);copy(t,h,n,lim);f[n-1]=0;for(int i=0;i<n-1;i++) f[i]=t[i+1]*(i+1)%mod;fill(f+n,f+lim,0);return;
}
void jifen(ll *h,ll *f,int n){static ll t[N];static int lim;init(n<<1,lim);copy(t,h,n,lim);f[0]=0;for(int i=1;i<n;i++) f[i]=h[i-1]*ksm(i,mod-2)%mod;fill(f+n,f+lim,0);
}
void Ln(ll *h,ll *f,int n){static ll t1[N],t2[N];static int lim;init(n<<1,lim);inv(h,t1,n);fill(t1+n,t1+lim,0);dao(h,t2,n);mul(t1,t2,t1,n,n);jifen(t1,f,n);return;
}
void Exp(ll *h,ll *f,int n){static ll t1[N],t2[N],t3[N];static int lim;if(n==1){f[0]=1;return;}Exp(h,f,(n+1)>>1);init(n<<1,lim);fill(f+((n+1)>>1),f+lim,0);copy(t1,f,n,lim);copy(t2,h,n,lim);Ln(f,t3,n);fill(t3+n,t3+lim,0);NTT(t1,lim,1);NTT(t2,lim,1);NTT(t3,lim,1);for(int i=0;i<lim;i++) t1[i]=t1[i]*(1+t2[i]-t3[i]+mod)%mod;NTT(t1,lim,-1);memcpy(f,t1,sizeof(ll)*n);return;
}
struct KSM{int w,x;ll u[N],v[N];void init(int base){x=base;w=sqrt(mod);v[0]=1;for(int i=1;i<w;i++) v[i]=v[i-1]*x%mod;ll o=v[w-1]*x%mod;u[0]=1;for(int i=1;i<=mod/w;i++) u[i]=u[i-1]*o%mod;		return;}inline ll calc(int k){int a=k/w,b=k%w;return u[a]*v[b]%mod;}
}a,b;
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;
}
signed main() {
#ifndef ONLINE_JUDGE//freopen("a.in","r",stdin);//freopen("a.out","w",stdout);
#endifint T=read();init();a.init(94153035);b.init(905847205);ll ans(0);while(T--){ll o=Rand()%(mod-1);ans^=(233230706ll*(a.calc(o)+mod-b.calc(o))%mod);}printf("%lld\n",ans);return 0;
}
/**/

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.mzph.cn/news/317144.shtml

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

Product 1 Modulo N CodeForces - 1514C

Product 1 Modulo N CodeForces - 1514C 题意&#xff1a; 在[1,n-1]中选x个数&#xff0c;使得乘积mod n 1&#xff0c;求x的最大值&#xff0c;并输出所选的数 题解&#xff1a; 我们设S为所选x个数的乘积 S%n 1说明gcd(S,n)1,即所选的x个数均与n互质&#xff0c;如果不…

k8s使用helm打包chart并上传到腾讯云TencentHub

本文只涉及Helm的Chart操作&#xff0c;不会对其他知识进行过多描述。至于安装这块&#xff0c;麻烦自行百度吧&#xff0c;一大堆呢。在容器化的时代&#xff0c;我们很多应用都可以部署在docker&#xff0c;很方便&#xff0c;而再进一步&#xff0c;我们还有工具可以对docke…

数据结构之基环树——骑士,Island,旅行加强版,Number of Simple Paths,Traffic Network in Numazu,Card Game

文章目录[ZJOI2008]骑士[IOI2008] Island[NOIP2018 提高组] 旅行 加强版CF1454E Number of Simple PathsTraffic Network in NumazuCard Game基环树的常见解法若干个基环树互相独立断环为链&#xff08;随便断一条&#xff09;环外树和环外树之间的树形DP环变链后整体可以用数据…

1090. 绿色通道

1090. 绿色通道 题意&#xff1a; n个题&#xff0c;每个题所花时间为ai&#xff0c;最多只用不超过t分钟做这个&#xff0c;肯定会有题目做不完&#xff0c;下标连续的一些空题称为一个空题段&#xff0c;问在t分钟内最长的空题段至少有多长&#xff1f; 题解&#xff1a; …

模板:BSGS(数论)

所谓 BSGS&#xff0c;就是北上广深。 &#xff08;逃&#xff09; BSGS 给出 a,b,pa,b,pa,b,p&#xff0c;请处出满足 ax≡b(modp)a^x\equiv b\pmod pax≡b(modp) 的最小非负正数解或者报告无解。 a,b,p≤109,gcd⁡(a,p)1a,b,p\le 10^9,\gcd(a,p)1a,b,p≤109,gcd(a,p)1 由于 …

如何在ASP.NET Core中自定义Azure Storage File Provider

主题&#xff1a;如何在ASP.NET Core中自定义Azure Storage File Provider作者&#xff1a; Lamond Lu地址: https://www.cnblogs.com/lwqlun/p/10406566.html项目源代码&#xff1a; https://github.com/lamondlu/AzureFileProvider背景ASP.NET Core是一个扩展性非常高的框架…

splay/fhq-treap 问卷调查反馈—— [JSOI2008]火星人prefix(splay),Strange Queries(fhq-treap)

文章目录[JSOI2008]火星人prefixStrange Queries[JSOI2008]火星人prefix BZOJ1014 思路很好想&#xff0c;哈希字符串即可 只是平衡树的码量大 注意因为splay加入哨兵的原因&#xff0c;每个点在平衡树内的排名比真实排名大111&#xff08;有极小值的占位&#xff09; 考虑…

AcWing 1091. 理想的正方形

AcWing 1091. 理想的正方形 题意&#xff1a; 有一个 ab 的整数组成的矩阵&#xff0c;现请你从中找出一个 nn 的正方形区域&#xff0c;使得该区域所有数中的最大值和最小值的差最小。 题解&#xff1a; 前置知识&#xff1a;已经学会了一维的单调队列优化dp 在本题中要求…

美好生活从撸好代码开始

楔子 昨天晚上做了个梦&#xff0c;梦到老板对我说了一番道理&#xff0c;他说对家庭要用爱心&#xff0c;做人对社会要有包容心&#xff0c;对工作要有责任心&#xff0c;对老板要有同理心。 我深以为然。现在的老板确实太不容易了&#xff0c;尤其是作为一家承载梦想&#xf…

模板:莫比乌斯反演(数论)

文章目录前言整除分块代码积性函数线性筛狄利克雷卷积莫比乌斯反演trick所谓莫比乌斯反演&#xff0c;就是莫比乌斯进行的反演 &#xff08;逃&#xff09; 前言 在一些需要整除的式子和 gcd⁡,lcm⁡\gcd,\operatorname{lcm}gcd,lcm 等问题中发挥作用。 整除分块 整除分块是…

[TJOI2013]拯救小矮人(反悔贪心证明),「ICPC World Finals 2019」Hobson 的火车(基环树,差分)

2021-09-07 test[TJOI2013]拯救小矮人「ICPC World Finals 2019」Hobson 的火车[TJOI2013]拯救小矮人 luogu4823 考试题目的数据加强为2e5&#xff0c;所以此题做法应为O(nlog⁡n)O(n\log n)O(nlogn)的反悔贪心 这种有多元属性&#xff0c;选择最优的问题 如果发现简单的贪心…

Dotnet全平台下APM-Trace探索

随着支撑的内部业务系统越来越多&#xff0c;向着服务化架构进化&#xff0c;在整个迭代过程中&#xff0c;会逐渐暴露出以下问题。传统依赖于应用服务器日志等手段的排除故障原因的复杂度越来越高&#xff0c;传统的监控服务已经无法满足需求。终端--> Nginx --> IIS --…

生成函数全家桶

文章目录有用的式子1.&#xff08;牛顿二项式定理&#xff09;2.普通生成函数&#xff08;OGF&#xff09;常见封闭形式&#xff1a;1.2.3.4.指数生成函数&#xff08;EGF&#xff09;排列与圆排列有用的式子 1.&#xff08;牛顿二项式定理&#xff09; 我们把组合数的定义推…

2020年牛客多校第五场C题-easy(纯组合计数不要生成函数的做法)

文章目录descriptionsolutioncodedescription 有TTT组测试数据 对于两个长度为KKK的数列{a}\{a\}{a}和{b}\{b\}{b}&#xff0c;满足∑i1KaiN,∑i1KbiM\sum_{i1}^Ka_iN,\sum_{i1}^Kb_iM∑i1K​ai​N,∑i1K​bi​M 对于这两个数列&#xff0c;定义权值为P∏i1Kmin⁡(ai,bi)P\p…

部署Chart应用并使用.net core读取Kubernetes中的configMap

上一篇文章讲了 k8s使用helm打包chart并上传到腾讯云TencentHub&#xff0c;今天就讲一下使用Helm部署应用并使用configMap代替asp.net core 中的appsettings.json文件。把Chart上传到TencentHub之后&#xff0c;我们就可以通过腾讯云的容器服务&#xff0c;直接部署Helm应用了…

Vases and Flowers HDU - 4614

Vases and Flowers HDU - 4614 题意: 一排空瓶子放花&#xff0c;操作1&#xff1a;从第x个瓶子开始放花&#xff0c;放y朵花&#xff0c;每个瓶子就一朵花&#xff0c;如果碰到已经有花的瓶子跳过这个瓶子&#xff0c;看下一个&#xff0c;当花没了&#xff0c;或者瓶子不够…

洛谷P3327:[SDOI2015]约数个数和(莫比乌斯反演)

枚举倍数的一种灵活的变形&#xff1a;g(d)∑d∣inf(i)∑i1⌊nd⌋f(i⋅d)g(d)\sum_{d|i}^nf(i)\sum_{i1}^{\lfloor\frac{n}{d}\rfloor}f(i\cdot d)g(d)∑d∣in​f(i)∑i1⌊dn​⌋​f(i⋅d) 很显然&#xff0c;但有时能发挥大作用。 其实本质还是要理解西格玛究竟是在算什么 解析…

EFCore Lazy Loading + Inheritance = 干净的数据表 (一)

前言α角 与 β角关于α角 与 β角的介绍&#xff0c;请见上文 如何用EFCore Lazy Loading实现Entity Split。本篇会继续有关于β角的彩蛋在等着大家去发掘。/斜眼笑其他本篇的程序&#xff0c;可以在 https://github.com/kentliu2007/EFCoreDemo/tree/master/InheritanceWithE…

专题突破之反悔贪心——建筑抢修,Cow Coupons G, Voting (Hard Version),Cardboard Box

文章目录[JSOI2007]建筑抢修[USACO12FEB]Cow Coupons GCF1251E2 Voting (Hard Version)CF436E Cardboard Box[JSOI2007]建筑抢修 luogu4053 将建筑按照结束时间从小到大排序 然后记录一下已经修理的建筑总共的花费时间 如果花费时间加上现在这个建筑的修建时间超过了这个建…

Max Sum Plus Plus HDU - 1024

Max Sum Plus Plus HDU - 1024 题意&#xff1a; 给你n个数&#xff0c;选m个子段&#xff0c;各个子段连续且不相交&#xff0c;长度可以为1&#xff0c;设maxn为各个子区间的和&#xff0c;求最大的maxn。 题解&#xff1a; 设dp[i][j]表示前j个数分成i段的最大值 对于第…