洛谷P5110:块速递推(特征根方程、光速幂)
解析
去你的搬磚生成函數,特征根太香了。
一開始我是用生成函數解的,和特征根相比有億點點搬磚…
但是這個東西原理似乎使用一些神奇的等比差分,有些玄學,生成函數較易理解。
背下來背下來!
就以本題為情境講一下特征根方程的解題流程吧。
對于遞推式:
an=233an?1+666an?2a_n=233a_{n-1}+666a_{n-2}an?=233an?1?+666an?2?
其特征方程就是:
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?=2233?56953??
那么通項公式就可以寫成
an=c1p1n+c2p2na_n=c_1p_1^n+c2p_2^nan?=c1?p1n?+c2p2n?
當解得的 p1=p2=pp_1=p_2=pp1?=p2?=p 時,通項會變成:an=(c1+c2n)pna_n=(c_1+c_2n)p^nan?=(c1?+c2?n)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=c1?p1?+c2?p2?
解得
c1=156953c_1=\frac{1}{\sqrt{56953}}c1?=56953?1?
c2=?156953c_2=-\frac{1}{\sqrt{56953}}c2?=?56953?1?
帶回原式:
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?(2233?56953??)n?
然后有一個很重要的技巧:尋找 569535695356953 在對 1e9+71e9+71e9+7 的模意義下是否是二次剩余。
出題人總會給人一條活路,寫個暴力跑一會就能找到:
1883058372≡56953(mod109+7)188305837^2\equiv 56953\pmod{10^9+7}1883058372≡56953(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?(2233?188305837?)n?
化簡得:
an=94153035n?905847205n188305837a_n=\frac{94153035^n-905847205^n}{188305837}an?=18830583794153035n?905847205n?
再把逆元提前算出來:
an=233230706(94153035n?905847205n)a_n=233230706(94153035^n-905847205^n)an?=233230706(94153035n?905847205n)
其中的指數 nnn 可以用歐拉定理縮小到 O(mod)O(mod)O(mod) 范圍。
然后問題就變成了如何快速求這個東西。
可以使用光速冪
原理很好理解,使用類似分塊的思路,設 p=modp=\sqrt{mod}p=mod?,O((mod))O(\sqrt(mod))O((?mod)) 預處理所有 xkp(k≤?modp?)x^{kp}(k\le\lfloor \frac{mod}{p}\rfloor)xkp(k≤?pmod??) 和 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}x?pk???p×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; } /**/總結
以上是生活随笔為你收集整理的洛谷P5110:块速递推(特征根方程、光速幂)的全部內容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: AI怎么做全开海报(ai怎么做全开海报图
- 下一篇: ddos攻击能够防止吗(ddos攻击如果