题解 P5126 【鬼故事】

NaCly_Fish · 2019-08-14 19:23:15 · 题解

update：由于这题标算被爆了多次，就详细详细讲一下各做法吧（合订本既视感）。

出题人做法：

构造一个矩阵来做快速幂，时间复杂度 \Theta(k^3 \log n)。

本人第一次做法：

既然这是线性递推，那我们可以用 BM 算法找出这个递推式，利用优化的线性递推算法，时间复杂度为 \Theta(k^2 \log n) 或 \Theta(k \log k \log n)。

提交记录

能不能给力一点啊？

第二次看这题的时候，已经学习了这题 EI 的特征根性质优化做法（orz EI！），下面做详细讲解。

设

A=\frac{1}{\sqrt 5} \ , \ \alpha=\frac{1+\sqrt 5}{2}

则

\prod_{i=n}^{n+k-1}f_i=\prod_{i=0}^{k-1}(A\alpha^{n+i}-A(-\alpha)^{-n-i})

=A^k\prod_{i=0}^{k-1}(\alpha^n\alpha^i-(-\alpha)^{-n}(-\alpha)^{-i})

=A^k \alpha^{k(k-1)/2} \prod_{i=0}^{k-1}(\alpha^n-(-\alpha)^{-n}(-\alpha^{-2})^i)

分类讨论一下 n 的奇偶性，这个式子就能看成是关于 \alpha^n 的一个 k 次多项式。

再用 x 来代换就有两个多项式：

\frac{A^k\alpha^{k(k-1)/2}}{x^k}\prod_{i=0}^{k-1}(x^2-(-\alpha^{-2})^i) \ \ (n \bmod 2 = 0)

\frac{A^k\alpha^{k(k-1)/2}}{x^k}\prod_{i=0}^{k-1}(x^2+(-\alpha^{-2})^i) \ \ (n \bmod 2 = 1)

设这两个多项式的 n 次项系数分别为 p_n,q_n，答案就是：

\left(\sum_{i=m}^n[2|i]\sum_{j=-k}^kp_j \alpha^{ij}\right)+\left(\sum_{i=m}^n[2|(i-1)]\sum_{j=-k}^kq_j\alpha^{ij} \right)

=\frac 12\sum_{j=-k}^k\left((p_j+q_j)\sum_{i=m}^n\alpha^{ji}\right)+\left( (p_j-q_j)\sum_{i=m}^n(-\alpha^j)^i\right)

还有一点常数优化，就是 p_n=(-1)^{k-n}q_n，求出 q 后就可以直接求出 p。

现在考虑如何求出那两个多项式的系数，只考虑后面那个乘积，是这样一个形式：

f(x)=\prod_{i=0}^{k-1}(x+\beta^i)

（虽然原式中是 x^2，但是没有关系，算出这个后第 n 项对应原式 2n 项。）

这个式子可以倍增处理，以 \Theta(k \log k) 的时间复杂度求出系数 —— 不过这是不必要的，它满足典型的 q-整式递推，我们这样操作：

f(\beta x)= \beta^k \prod_{i=0}^{k-1}(x+\beta^{i-1})=\beta^k \frac{x+\beta^{-1}}{x+\beta^{k-1}}f(x)

(x+\beta^{k-1})f(\beta x)=\beta^k(x+\beta^{-1})f(x)

提取系数即得

\beta^{i-1}f_{i-1}+\beta^{k+i-1}f_i=\beta^kf_{i-1}+\beta^{k-1}f_i

f_{i-1}= \frac{\beta^{k-1}-\beta^{k+i-1}}{\beta^{i-1}-\beta^k}f_i

从上边界 f_k=1 往下递推即可，线性处理逆元即可 \Theta(k) 求出 f(x) 的系数。统计答案的部分涉及等比数列求和，可以只计算常数次快速幂，其余都能线性处理。

当然，你会发现计算快速幂时只有两次指数特别大，而

a_n =\alpha^n \bmod p

显然是有循环节的，所以可以直接把指数对 4p(p+1) 取模再计算。

总时间复杂度为 \Theta(k +\log n)。

参考代码（稍微改改参数就能直接通过加强版）：

#include<cstdio>
#include<iostream>
#include<algorithm>
#include<cstring>
#include<cmath>
#define N 20003
#define ll long long
#define p 1000000007
using namespace std;

inline int add(const int& x,const int& y){ return x+y>=p?x+y-p:x+y; }
inline int dec(const int& x,const int& y){ return x<y?x-y+p:x-y; }

inline int intpow(int a,int t){
    int res = 1;
    while(t){
        if(t&1) res = (ll)res*a%p;
        a = (ll)a*a%p;
        t >>= 1;
    }
    return res;
}

struct complex{
    int r,i;
    inline complex(int _r=0,int _i=0):r(_r),i(_i){}
};

inline complex operator + (const complex& lhs,const complex& rhs){ return complex(add(lhs.r,rhs.r),add(lhs.i,rhs.i)); }
inline complex operator - (const complex& lhs,const complex& rhs){ return complex(dec(lhs.r,rhs.r),dec(lhs.i,rhs.i)); }
inline complex operator * (const complex& lhs,const complex& rhs){ return complex(((ll)lhs.r*rhs.r+5ll*lhs.i*rhs.i)%p,((ll)lhs.r*rhs.i+(ll)lhs.i*rhs.r)%p); }
inline complex operator - (const complex& x){ return complex(dec(0,x.r),dec(0,x.i)); }
inline complex operator * (const complex& lhs,const int& rhs){ return complex((ll)lhs.r*rhs%p,(ll)lhs.i*rhs%p); } 
const complex one = complex(1,0);

inline complex power(complex a,ll t){
    complex res = complex(1,0);
    while(t){
        if(t&1) res = res*a;
        a = a*a;
        t >>= 1;
    }
    return res;
}

inline complex inverse(complex x){
    int d = ((ll)x.r*x.r-5ll*x.i*x.i%p+p)%p;
    return complex(x.r,p-x.i)*intpow(d,p-2);
}

complex pre[N>>1];

void product(int k,complex q,complex *t){
    t[k] = 1;
    static complex pw[N];
    pw[0] = pre[0] = one;
    for(int i=1;i<=k;++i) pw[i] = pw[i-1]*q;
    for(int i=1;i<=k;++i) pre[i] = pre[i-1]*(pw[k]-pw[i-1]);
    complex Inv = inverse(pre[k]),suf = one;
    for(int i=k;i;--i){
        t[i-1] = pw[k-1]*(pw[i]-one)*(Inv*pre[i-1]*suf)*t[i];
        suf = suf*(pw[k]-pw[i-1]);
    }
}

char L[N],R[N];
int k;
complex f[N],g[N],inv1[N>>1],inv2[N>>1],suf[N>>1];
const int ninv2 = p-(p-1)/2;
const ll md = 4ll*p*(p+1);
const __int128_t ten = 10;

int main(){
    scanf("%d",&k);
    scanf("%s%s",L,R);
    int lenl = strlen(L),lenr = strlen(R),lmdp = 0,rmdp = 0;
    ll pwl = 0,pwr = 0;
    for(int i=0;i<lenl;++i) lmdp = (lmdp*10ll+L[i]-'0')%p,pwl = (pwl*ten+L[i]-'0')%md;
    for(int i=0;i<lenr;++i) rmdp = (rmdp*10ll+R[i]-'0')%p,pwr = (pwr*ten+R[i]-'0')%md;
    reverse(L,L+lenl),reverse(R,R+lenr);
    int inv5 = intpow(5,p-2);
    complex alpha = complex(ninv2,ninv2),A = complex(0,inv5);
    complex iva = inverse(alpha),q,beta,ans = 0,a2 = alpha*alpha;
    q = power(iva,k),beta = -iva*iva;
    product(k,beta,g);
    for(int i=0;i<=k;++i) f[i] = (k-i)&1?-g[i]:g[i];
    complex pal = power(alpha,pwl),par = power(alpha,pwr+1);
    complex pa2l = pal*pal,pa2r = par*par,pql = power(inverse(pal),k),pqr = power(inverse(par),k);
    complex pali = one,pari = one,qs,nqs,npql = (L[0]&1)?-pql:pql,npqr = (R[0]&1)?pqr:-pqr;
    inv1[0] = q,inv2[0] = q+one;
    for(int i=1;i<=k;++i){
        inv1[i] = inv1[i-1]*a2;
        inv2[i] = inv1[i]+one;
    }
    for(int i=0;i<=k;++i) inv1[i] = inv1[i]-one;
    if(!(k&1)) inv1[k>>1] = one;
    pre[0] = suf[k+1] = one;
    inv1[0] = inverse(inv1[0]);
    for(int i=1;i<=k;++i) pre[i] = pre[i-1]*inv1[i];
    for(int i=k;i;--i) suf[i] = suf[i+1]*inv1[i];
    complex mul = inverse(pre[k]);
    for(int i=1;i<=k;++i) inv1[i] = pre[i-1]*suf[i+1]*mul;
    inv2[0] = inverse(inv2[0]);
    for(int i=1;i<=k;++i) pre[i] = pre[i-1]*inv2[i];
    for(int i=k;i;--i) suf[i] = suf[i+1]*inv2[i];
    mul = inverse(pre[k]);
    for(int i=1;i<=k;++i) inv2[i] = pre[i-1]*suf[i+1]*mul;
    for(int i=0;i<=k;++i){
        if(q.r==1&&q.i==0) qs = dec(rmdp+1,lmdp);
        else qs = (pqr*pari - pql*pali)*inv1[i];
        nqs = (npqr*pari - npql*pali)*inv2[i];
        ans = ans + (f[i]+g[i])*qs - (f[i]-g[i])*nqs;
        pali = pali*pa2l,pari = pari*pa2r;
        q = q*a2;
    }
    ans = ans*power(A,k)*power(alpha,k*(k-1)>>1)*ninv2;
    printf("%d",ans.r);
    return 0;   
}

关于「能不能再给力一点啊？」这个问题，我对 q-整式递推还不甚了解，就不妄议了。如果哪天复杂度再有优化，还会回来补的（