P5126 鬼故事

P5126 鬼故事

题目背景

给大家讲个鬼故事

有一天晚上，下着大暴雨。小 K 正在他那小小的书房里做着老师给他布置的、数也数不完的信息题目。关着门，关着窗难免有些闷热。小 K 起身，将书桌前的窗子打开一个小小的细缝，小到没有雨点能透过缝隙飘进来。

今天正是农历七月十五日，中元节，俗称鬼节。小 K 从来都没有在这种日子这么晚睡过，因为小 K 迷信，害怕午夜之后，便有鬼怪出没。然而今天，小 K 无可奈何。

小 K 看了看时间：$23:54$。看到 $4$，小 K 皱了皱眉眉头。$4$，谐音是“死”，特别不吉利。在这种日子看到这样的字眼，往往都是不祥之兆。

小 K 的眼皮在打架。他从来都不会做毒瘤题。他索性趴到了书桌上面，两只眼睛渐渐朦胧了起来。

“那儿有一个本子。”他想着。不知道何时，他的书桌靠窗的一角上，静静地躺着一个湿漉漉的本子，好像是刚刚淋过雨。“它是怎么进来的？”小 K 喃喃道。他下意识地翻开那本本子，看到里面有写了一些字。不知道为什么，那些字在发黄的纸页上看起来也那么红。

左边的那一页写着：

$$4^{4-4}\le M\le N\le 4^{4^{(4+4-\frac{4}{4})}},\sqrt{4}\le K\le 4^{\sqrt{4}}\times (4-\frac{4}{4})+\sqrt{4}$$

右边的那一页写着的似乎比左边的要长：

$$a_{\frac{4}{4}}=a_{\sqrt{4}}=\frac{4}{4},a_n=\frac{\sqrt{4}}{\sqrt{4}}\times a_{n-\frac{4}{4}}+4^{4-4}\times a_{n-\sqrt{4}}(n\ge \sqrt{4}+\frac{4}{4})$$

$$b_n=\prod _{i=n}^{n+K-\frac{4}{4}}{a}_i$$

求 $\sum _{i=m}^n{b}_i$

角落里还有一行小字：_ 不要翻到最后一页，不然会有可怕的事情发生 _。但是在这个时候，小 K 早已经闭上了双眼，鼾声和远处的雷声混成一片。

一阵微风吹来，轻轻地，谁也没有意识到。本子的一角被风扬起，滑过一个优美的弧线，落在了本子的另一边。风一阵一阵的吹来，拂过本子发黄的书页。渐渐地，右边的书页少了，左边的书页多了。风停了，本子的倒数第二页停在半空中。在一刹那，似乎一切都静止了。然后，它轻轻地落在了其它书页的最上面。

最后一页上，赫然用鲜红色的歪歪扭扭的大字写着：

这道题你已经拖了一个月了！限明天之前做完！

这时候，你夜观天象，预测到了小 K 的这场劫难。时间已是 $23:59:59:400$，如果在这 $1000-400=600$ 毫秒内没有做完，小 K 的检讨将在劫难逃。身为小 K 的好朋友，你能帮他解决这个问题吗？

题目描述

给定 $k,m,n$，求：

$$\sum _{i=m}^n\prod _{j=i}^{i+k-1}{a}_j$$

答案对 $10^9+7$ 取模。
其中 { a } \{ a\} {a} 为 fibonacci 数列。

输入格式

三个正整数，分别表示 $k,m,n$。

输出格式

输出一行一个整数表示答案。

输入输出样例 #1

输入 #1

4 1 3

输出 #1

276

输入输出样例 #2

输入 #2

3 2 3

输出 #2

36

说明/提示

$a_1=1,a_2=1,a_3=2,a_4=3,a_5=5,a_6=8$。

对于样例1：

$$K=4$$
$$b_1=1\times 1\times 2\times 3=6,b_2=1\times 2\times 3\times 5=30,b_3=2\times 3\times 5\times 8=240$$
$$\sum _{i=1}^3{b}_i=276$$

对于样例2：

$$K=3$$
$$b_2=1\times 2\times 3=6,b_3=2\times 3\times 5=30$$
$$\sum _{i=2}^3{b}_i=36$$

本题共有 $20$ 个数据点，每个数据点的分数均为 $5$ 分，总分为 $100$ 分。每个数据点的性质如下：

(出题人不想再用 $4$ 表示任何数了！真香)

编号	$K,M,N$范围	特殊性质
$1$	$1\le m\le n\le 10^6,k=4$	无
$2$	$1\le m\le n\le 10^{18},k=4$	$n-m\le 10^6$
$3\sim 4$	$1\le m\le n\le 10^{18},k=4$	无
$5\sim 6$	$1\le m\le n\le 4^{4^4},k=4$	$n-m\le 10^6$
$7\sim 10$	$1\le m\le n\le 4^{4^7},k=4$	无
$11\sim 12$	$1\le m\le n\le 4^{6000},2\le k\le 10$	无
$13\sim 14$	$1\le m\le n\le 10^{41},2\le k\le 10$	无
$15\sim 20$	$1\le m\le n\le 10^{41},2\le k\le 50$	无

（注意，题面中的数据范围只是大致描述，请以以上具体范围为准）

$a^{b^c}=a^{(b^c)}$

答案

#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;   
}