ABC381G

· · 题解

Solution

为方便描述，我们让 a 从 0 开始（即令 a_0=x,a_1=y，后同）。

考虑 a 的通项。具体地，我们设 A(z)=\sum_{n\ge 0}a_nz^n，那么有：

A(z)=zA(z)+z^2A(z)-xz+x+yz\\ A(z)=\frac{x+(y-x)z}{1-z-z^2}

使用具体数学 7.3 中的方法（可以在我的 博客 中查看），我们可以得到 a 的通项：记 c_1=\frac{1}{\sqrt 5}(\frac{\sqrt 5-1}{2}\cdot x+y),c_2=\frac{1}{\sqrt 5}(\frac{\sqrt 5+1}{2}\cdot x-y)，则有

a_n=c_1\Big(\frac{1+\sqrt 5}{2}\Big)^n+c_2\Big(\frac{1-\sqrt 5}{2}\Big)^n

考虑 a 的在 {}\bmod 998244353 意义下的循环节。打表发现，a_{1996488708}=x（即 a_{2\times (998244353+1)}=x）。这意味着我们可以将 n 对 1996488708 取模，只需处理 n<1996488708 的情况。

那么现在我们要求的即为：

\prod_{i=0}^{n-1}(c_1A^i+c_2B^i)

其中 A=\frac{1+\sqrt 5}{2},B=\frac{1-\sqrt 5}{2}。

下面我们来计算这个式子。令 D=\frac{A}{B}，则原式即为

\prod_{i=0}^{n-1}B^i(c_1D^i+c_2)=B^{\frac{n(n-1)}{2}}\prod_{i=0}^{n-1}(c_1D^i+c_2)

考虑根号分治计算后面那个式子。设块长为 m。下文假设 m|n，若不是直接 O(m) 调整即可。则原式即为

\prod_{i=0}^{\frac{n}{m}-1}\prod_{j=0}^{m-1}(c_1D^{im+j}+c_2)=\prod_{i=0}^{\frac{n}{m}-1}\prod_{j=0}^{m-1}(c_1D^jX^i+c_2)

其中 X=D^m。

设 F(z)=\prod_{j=0}^{m-1}(c_1D^jz+c_2)，则我们可以 O(m\log^2 m) 分治 NTT 或 O(m\log m) 倍增展开 F(z)。

关于倍增展开的说明：设我们已经计算出了 F(z)\bmod z^{m}，现在计算 F(z)\bmod z^{2m}。记前者为 F_0(z)，后者为 F_1(z)，那么有：

F_1(z)&=\Big(\prod_{i=0}^{m-1}(c_1D^iz+c_2)\Big)\times\Big(\prod_{i=m}^{2m-1}(c_1D^iz+c_2)\Big)\\ &=F_0(z)\times \prod_{i=0}^{m-1}(c_1D^i\times D^mz+c_2)\\ &=F_0(z)\times F_0(D^mz) \end{aligned}

也就是说我们将 F_0(z) 的第 i 项乘上 D^{im} 就能得到右半边的系数表示。

现在原问题可以表述为

\prod_{i=0}^{\frac{n}{m}-1}F(X^i)

我们可以通过 O(\frac{n}{m}\log^2 \frac nm) 的多项式多点求值或 O(\frac nm\log \frac nm) 的 Chirp Z Transform 求出每个 F(X^i)。

取 m=\sqrt n，有最优时间复杂度即为 O(\sqrt n\log n)，可以通过。

这里讲一下 Chirp Z Transform：

设我们要求的是 F(1),F(X),F(X^2),\cdots,F(X^p)。考虑这样一个事实：

ki={k+i\choose 2}-{i\choose 2}-{k\choose 2}

组合意义不难证明。应用到我们的式子中：

F(X^k)&=\sum_{i=0}^{\deg F}[z^i]F(z)\times X^{ki}\\ &=\sum_{i=0}^{\deg F}[z^i]F(z)\times X^{{k+i\choose 2}-{i\choose 2}-{k\choose 2}}\\ &=X^{-{k\choose 2}}\sum_{i=0}^{\deg F}[z^i]F(z)\times X^{{k+i\choose 2}-{i\choose 2}} \end{aligned}

令 F_0(z)=\sum_{i=0}^{\deg F}[z^i]F(z)\times X^{-{i\choose 2}}，F_1(z)=\sum_{i=0}^{p+\deg F}X^{{^{i\choose 2}}}，则不难通过一次差卷积计算出所有 F(X^k)。

另外，由于 5 不是模 998244353 意义下的二次剩余，我们需要扩域。单位根沿用 3 即可，毕竟 NTT 的过程中并没有用到单位根的幂必须取遍整个域这个性质。

Code

#include<bits/stdc++.h>
using namespace std;
using ui=unsigned int;
using ll=long long;
using ull=unsigned long long;
using i128=__int128;
using u128=__uint128_t;
using pii=pair<int,int>;
#define fi first
#define se second
constexpr int N=262145,mod=998244353,loop=(mod+1)<<1,g=3,inv2=(mod+1)>>1;
ll n,m,x,y;
inline ll qpow(ll a,ll b){
    ll res=1;
    for(;b;b>>=1,a=a*a%mod)
        if(b&1)res=res*a%mod;
    return res;
}
inline ll add(ll x,ll y){return (x+=y)>=mod&&(x-=mod),x;}
inline ll Add(ll &x,ll y){return x=add(x,y);}
inline ll sub(ll x,ll y){return (x-=y)<0&&(x+=mod),x;}
inline ll Sub(ll &x,ll y){return x=sub(x,y);}
struct node{ // a + b * sqrt(5)
    ll a,b;
    node(){a=b=0;}
    node(ll _a,ll _b){a=_a,b=_b;}
    inline node operator +(ll _)const{return node(add(a,_),b);}
    inline node operator -(ll _)const{return node(sub(a,_),b);}
    inline node operator *(ll _)const{return node(a*_%mod,b*_%mod);}
    inline node operator +(const node &_)const{return node(add(a,_.a),add(b,_.b));}
    inline node operator -(const node &_)const{return node(sub(a,_.a),sub(b,_.b));}
    inline node operator *(const node &_)const{return node((a*_.a+b*_.b*5)%mod,(a*_.b+b*_.a)%mod);}
    inline node inv()const{
        static ll tmp;tmp=qpow((a*a-b*b*5%mod+mod)%mod,mod-2);
        return node(a*tmp%mod,(mod-b)*tmp%mod);
    }
    inline friend node operator /(const node &a,const node &b){return a*b.inv();}
};
inline ostream& operator <<(ostream& ouf,const node &x){
    return ouf<<'('<<x.a<<", "<<x.b<<')';
}
const node A(inv2,inv2),B(inv2,mod-inv2),D=A/B,I(0,1);
node c1,c2;
inline node qpow(node a,ll b){
    node res(1,0);
    for(;b;b>>=1,a=a*a)
        if(b&1)res=res*a;
    return res;
}
int lim,ilim,rev[N],w[N];
void initNTT(int n){
    lim=1;
    while(lim<=n)lim<<=1;
    ilim=qpow(lim,mod-2);
    for(int i=1;i<lim;i++)
        rev[i]=(rev[i>>1]>>1)|((i&1)*(lim>>1));
    for(int i=1;i<lim;i<<=1){
        ll wn=qpow(g,(mod-1)/(i<<1)),cur=1;
        for(int j=0;j<i;j++,cur=cur*wn%mod)
            w[i|j]=cur;
    }
}
using poly=vector<node>;
void NTT(node *F,int type){
    for(int i=0;i<lim;i++)
        if(i<rev[i])
            swap(F[i],F[rev[i]]);
    node x,y;
    for(int i=1;i<lim;i<<=1)
        for(int j=0;j<lim;j+=i<<1)
            for(int k=0;k<i;k++){
                x=F[j|k],y=F[i|j|k]*w[i|k];
                F[j|k]=x+y,F[i|j|k]=x-y;
            }
    if(type==1)return;
    reverse(F+1,F+lim);
    for(int i=0;i<lim;i++)F[i]=F[i]*ilim;
}
void mul(const poly &F,const poly &G,poly &H){
    static node A[N],B[N],C[N];
    int n=(int)F.size()-1,m=(int)G.size()-1;
    for(int i=0;i<=n;i++)A[i]=F[i];
    for(int i=0;i<=m;i++)B[i]=G[i];
    if(1ll*n*m<=1000){
        for(int i=0;i<=n+m;i++)C[i]=node(0,0);
        for(int i=0;i<=n;i++)
            for(int j=0;j<=m;j++)
                C[i+j]=C[i+j]+A[i]*B[j];
    }
    else{
        initNTT(n+m);
        for(int i=n+1;i<lim;i++)A[i]=node(0,0);
        for(int i=m+1;i<lim;i++)B[i]=node(0,0);
        NTT(A,1),NTT(B,1);
        for(int i=0;i<lim;i++)C[i]=A[i]*B[i];
        NTT(C,-1);
    }
    H.resize(n+m+1);
    for(int i=0;i<=n+m;i++)H[i]=C[i];
}
void solve1(poly &F,int n){
    if(n==1){
        F.resize(2);
        F[0]=c2,F[1]=c1;
        return;
    }
    static poly G;
    int m=n>>1;
    solve1(F,m);
    node cur(1,0),X=qpow(D,m);G.resize(m+1);
    for(int i=0;i<=m;i++,cur=cur*X)G[i]=cur*F[i];
    mul(F,G,F);
    if(n&1){
        F.emplace_back(0,0);
        node tmp=c1*qpow(D,n-1);
        for(int i=n;i>=0;i--)F[i]=(i>0?tmp*F[i-1]:node(0,0))+c2*F[i];
    }
}
inline ll C(int n){return n*(n-1ll)>>1;}
void solve2(const poly &F,const node &X,poly &G,int n){
    if(!n){G.resize(0);return;}
    static node iX;static poly F0,F1;
    iX=X.inv();
    F0.resize(F.size());
    F1.resize((int)F.size()+n-1);
    for(int i=0;i<(int)F0.size();i++)F0[i]=F[i]*qpow(iX,C(i));
    for(int i=0;i<(int)F1.size();i++)F1[i]=qpow(X,C((int)F1.size()-i-1));
    mul(F0,F1,F0);
    G.resize(n);
    for(int i=0;i<n;i++)G[i]=qpow(iX,C(i))*F0[(int)F1.size()-i-1];
}
poly F;
node calc(int n){
    if(!n)return node(1,0);
    m=sqrt(n);
    solve1(F,m);
    solve2(F,qpow(D,m),F,n/m);
    node ans(1,0);
    for(int i=0;i<(int)F.size();i++)ans=ans*F[i];
    node cur=qpow(D,n/m*m);
    for(int i=n/m*m;i<n;i++,cur=cur*D)ans=ans*(c1*cur+c2);
    return ans*qpow(B,C(n));
}
void solve(){
    cin>>n>>x>>y;
    c1=(node((mod-1)>>1,inv2)*x+y)/I;
    c2=(node(inv2,inv2)*x-y)/I;
    cout<<(qpow(calc(loop),n/loop)*calc(n%loop)).a<<'\n';
}
int main(){
    ios::sync_with_stdio(false);
    cin.tie(0),cout.tie(0);
    int _Test;cin>>_Test;
    while(_Test--)solve();
    return 0;
}

所以有没有人会证这个循环节的结论啊 /kel