自动搬运

查看原文

来自洛谷，原作者为

	NaCly_Fish 北海虽赊，扶摇可接。

搬运于2025-08-24 22:07:25，当前版本为作者最后更新于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-整式递推 还不甚了解，就不妄议了。如果哪天复杂度再有优化，还会回来补的（