高维标记

ll md(ll x)
{
    if(0<=x)
    {
        if(x<mod)return x;
        if(x<(mod<<1))return x-mod;
        return x%mod;
    }
    if(x<=-mod)x%=mod;
    return x+mod;
}

struct Matrix
{
    static const int n=3;
    int a[n+1][n+1];
    int* operator[](int x){return a[x];}

    Matrix& operator=(const Matrix &b)//赋值
    {
        for(int i=1;i<=n;++i)
            for(int j=1;j<=n;++j)
                a[i][j]=md(b.a[i][j]);
        return *this;
    }

    Matrix operator*(const Matrix &b)const//矩阵乘法
    {
        Matrix c;c.clear();
        for(int i=1;i<=n;++i)
            for(int k=1;k<=n;++k)//Cache
                for(int j=1;j<=n;++j)
                    c.a[i][j]=md(c.a[i][j]+md(1ll*a[i][k]*b.a[k][j]));
        return c;
    }

    Matrix operator+(const Matrix &b)const//矩阵加法
    {
        Matrix c;c.clear();
        for(int i=1;i<=n;++i)
            for(int j=1;j<=n;++j)
                c.a[i][j]=md(a[i][j]+b.a[i][j]);
        return c;
    }

    void print()
    {
        for(int i=1;i<=n;++i)
            for(int j=1;j<=n;++j)
                pf("%d%c",a[i][j]," \n"[j==n]);
    }

    void clear()
    {
        for(int i=1;i<=n;++i)
            for(int j=1;j<=n;++j)
                a[i][j]=0;
    }
}m[__],sm[__];
//m[i]=pow(M,i),sm[i]=sm[i-1]+m[i];

struct Vector
{
    static const int n=Matrix::n;
    int a[n+1];

    int& operator[](int x){return a[x];}

    Vector& operator=(const Vector &b)//赋值
    {
        for(int i=1;i<=n;++i)
            a[i]=md(b.a[i]);
        return *this;
    }

    Vector operator*(const Matrix &b)const//向量乘矩阵
    {
        Vector c;c.clear();
        for(int j=1;j<=n;++j)
            for(int i=1;i<=n;++i)//Cache
                c.a[i]=md(c.a[i]+md(1ll*a[j]*b.a[j][i]));
        return c;
    }

    Vector& operator+=(const Vector &b)//向量加法
    {
        for(int i=1;i<=n;++i)
            a[i]=md(a[i]+b.a[i]);
        return *this;
    }

    bool empty()
    {
        for(int i=1;i<=n;++i)
            if(a[i])return false;
        return true;
    }

    void clear() {for(int i=1;i<=n;++i)a[i]=0;}
};

//f[i]=a*f[i-1]+b*f[i-2];       (a 1)
//(f[i],f[i-1])=(f[i-1],f[i-2])*(b 0)
void init(int ___,Vector x)//系数向量{0,f[i-1],f[i-2],……f[i-n]}
{
    const int n=Matrix::n;
    Matrix M;//系数阵
    for(int i=1;i<=n;++i)
    {
        M[i][1]=md(x[i]);
        for(int j=2;j<=n;++j)
            if(j==i+1)M[i][j]=1;
            else M[i][j]=0;
    }

    //单位阵
    for(int i=1;i<=n;++i)
        for(int j=1;j<=n;++j)
            if(i==j)m[0][i][j]=sm[0][i][j]=1;
            else m[0][i][j]=sm[0][i][j]=0;

    for(int i=1;i<=___;++i)
    {
        m[i]=m[i-1]*M;
        sm[i]=sm[i-1]+m[i];
    }
}

int a[__];

struct SegmentTree
{
    Vector val;int n,ql,qr;
    struct node
    {
        int val;Vector ad;

        void putadd(int tl,int tr,Vector v)
        {
            ad+=v;
            Vector x=v*sm[tr-tl];
            val=md(val+x[Matrix::n]);
        }

        void clear(){val=0;ad.clear();}
    }t[__<<2];

    void pushup(int x)
    {
        t[x].val=md(t[x<<1].val+t[x<<1|1].val);
    }

    void build(int _n){n=_n;build(1,1,n);}

    void build(int x,int tl,int tr)
    {
        t[x].clear();
        if(tl==tr){t[x].val=a[tl];return;}
        int tm=(tl+tr)>>1;
        build(x<<1,tl,tm);
        build(x<<1|1,tm+1,tr);
        pushup(x);
    }

    void pushdown(int x,int tl,int tm,int tr)
    {
        if(!t[x].ad.empty())
        {
            t[x<<1].putadd(tl,tm,t[x].ad);
            t[x<<1|1].putadd(tm+1,tr,t[x].ad*m[tm+1-tl]);
            t[x].ad.clear();
        }
    }

    void add(int _ql,int _qr,Vector _val)//{f[n],f[n-1],……,f[1]}
    {
        ql=_ql,qr=_qr,val=_val;
        _add(1,1,n);
    }

    void _add(int x,int tl,int tr)
    {
        if(ql>tr || qr<tl)return;
        if(ql<=tl && tr<=qr)
        {
            t[x].putadd(tl,tr,val*m[tl-ql]);
            return;
        }
        int tm=(tl+tr)>>1;
        pushdown(x,tl,tm,tr);
        _add(x<<1,tl,tm);
        _add(x<<1|1,tm+1,tr);
        pushup(x);
    }

    int get_val(int _ql,int _qr)
    {
        ql=_ql,qr=_qr;
        return _get_val(1,1,n);
    }

    int _get_val(int x,int tl,int tr)
    {
        if(ql>tr || qr<tl)return 0;
        if(ql<=tl && tr<=qr)return t[x].val;
        int res=0,tm=(tl+tr)>>1;
        pushdown(x,tl,tm,tr);
        res=md(_get_val(x<<1,tl,tm)+_get_val(x<<1|1,tm+1,tr));
        pushup(x);
        return res;
    }
}T;