[机房测试]NYG的序列拆分

Description

问长为 (n) ，每个数值域在 ([l,r]) 的所有 ((r-l+1)^n) 个序列的每个序列的拆分的方式之和。拆分合法当且仅当，将序列分为若干段，每一段都是一个等比数列。(1leq lleq rleq 10^7)，(nleq 10^{18})。

Solution

容易发现问题可以拆解为求不同长度的等比数列个数和将这些数列组合起来。

先抛开公比是一的情况。那么对于长为 (n) ，公比为 (frac{x}{y}) 的数列，一定长成

[ty^{n-1}，ty^{n-2}x，dots，tx^{n-1} ]

的样子。(n=1) 和 (n=2) 比较特殊。单独考虑，容易得到 (cnt_{1}=r-l+1) 和 (cnt_2=(r-l+1) imes(r-l))。对于其它情况，(ngeq 3)，那么有 (max(x^2,y^2)leq r)，所以只用在 (sqrt r) 下枚举。容易发现这样的数列长度最大为 (min(log_x n,log_y n))。那么如果枚举长度，就可以快速统计该种长度的数列个数 (max(0,lfloor frac{r}{x^{n-1}} floor -lfloor frac{l-1}{y^{n-1}} floor))。

如果题目保证公比不为一，那么就可以直接矩阵快速算答案了。因为有转移

[dp_i=sum_{j=1}^{lfloor log r floor} dp_{i-j} imes f_j ]

(f) 表示等比数列个数。

现在加上公比是一的情况，有

[dp_i=sum_{j=1}^{lfloor log r floor} dp_{i-j} imes f_j+(r-l+1) imessum_{j=2}^{i} dp_{i-j} ]

一晃眼，以为 (i) 的转移和所有的 (j<i) 都有关。但实际上需要的只是一个前缀和，所以直接考虑在矩阵里多维护一个 (S_{i-2}) 即可。

枚举 (x,y) 的时候要保证 ((x,y)=1)，所以有复杂度 (O(rlog r+log nlog^3 r))

#include<stdio.h>
#include<algorithm>
#include<string.h>
#include<math.h>
using namespace std;

typedef long long ll;

inline ll read(){
    ll x=0,flag=1; char c=getchar();
    while(c<'0'||c>'9'){if(c=='-')flag=0;c=getchar();}
    while(c>='0'&&c<='9'){x=(x<<1)+(x<<3)+c-48;c=getchar();}
    return flag? x:-x;
}

const int N=24;
const int Mod=1e9+7;

struct Matrix{
    ll a[N+1][N+1];
    Matrix(){memset(a,0,sizeof(a));}
    Matrix operator *(const Matrix &X) const{
        Matrix C;
        for(int k=0;k<=N;k++)
            for(int i=0;i<=N;i++)
                for(int j=0;j<=N;j++)
                    C.a[i][j]=(C.a[i][j]+a[i][k]*X.a[k][j]%Mod)%Mod;
        return C;
    }
    void print(){
        for(int i=0;i<=N;i++){
            for(int j=0;j<=N;j++)
                printf("%lld ",a[i][j]);
            putchar('
');
        }
    }
};

ll n,l,r,cnt[N];

ll gcd(ll x,ll y){return y? gcd(y,x%y):x;}

int main(){
    freopen("excellent.in","r",stdin);
    freopen("excellent.out","w",stdout);
    int T=read();
    while(T--){
        n=read(),l=read(),r=read();
        for(int i=1;i<N;i++) cnt[i]=0;
        int rg=(int)(sqrt(r)+0.5);
        cnt[1]=r-l+1,cnt[2]=(r-l+1)*(r-l)%Mod;
        for(int i=1;i<=rg;i++)
            for(int j=1;j<i;j++){
                if(gcd(i,j)!=1) continue;
                for(int pw=3,x=i*i,y=j*j;x<=r;x*=i,y*=j,pw++)
                    cnt[pw]=(cnt[pw]+max(0ll,r/x-(l-1)/y))%Mod;
            }
        for(int i=3;i<N;i++) cnt[i]=(cnt[i]<<1)%Mod;
        Matrix A,ret;
        for(int i=1;i<=N;i++) A.a[i-1][0]=cnt[i];
        for(int i=1;i<N;i++) A.a[i-1][i]=1;
        A.a[N][N]=A.a[0][N]=1,A.a[N][0]=r-l+1;
    //    A.print();
        for(int i=0;i<=N;i++) ret.a[i][i]=1;
        while(n){ if(n&1) ret=A*ret; A=A*A,n>>=1; }
        printf("%lld
",ret.a[0][0]);
    }
}