「UOJ500」任意基DFT

三月 27, 2020 · OI 题解

给定 $n$ 次多项式
$f(x) = \sum_{i=0}^n a_i x^i$
$Q$ 次询问，第 $i$ 次询问 $f(q_i)$ 对 $998244353$ 取模的值。

其中 $q_i$ 是一个一阶线性递推，给定 $q_0, x, y$ ，满足
$q_n = x q_{n-1} + y$
$1 \leq n \leq 2.5 \times 10^5, \ 1 \leq Q \leq 10^6, \ 2 \leq x < 998244353, \ 0 \leq q_0, y < 998244353$ 。

Solution Part1

首先这玩意儿肯定没法两个 $\log$ 多点求值，除非你是钱哥哥。

Solution Part2

虽然这东西好像小学生都会，但我自己推的时候怎么莫名搞到特征多项式那个方向去了，越学越傻逼了。。。

$\left\{ \begin{aligned} g_{n+1} = x g_n + y \\ g_n = x g_{n-1} + y \\ \end{aligned} \right.$

两式相减得到

$(g_{n+1}-g_n) = x (g_n - g_{n-1})$

不妨设

$h_n = g_{n+1} - g_n$

可以得到 $h$ 的通项公式

$h_n = x^n h_0 = x^n ((x-1) g_0 + y)$

从而得到 $g$ 的通项公式

$g_n = g_0 + \sum_{i=0}^{n-1} h_i = g_0 + \frac {x^n-1} {x-1} ((x-1)g_0 + y)$

不妨表示成 $ax^n + b$ 的形式，其中

$\begin{cases} a = \dfrac {(x-1) g_0 + y} {x-1} \\ b = \dfrac {y} {x-1} \\ \end{cases}$

Solution Part3

也就是说我们要对于所有 $i \in [1, Q]$ 求出 $f(a x^i + b)$ 的值，不妨设多项式 $g(x) = f(ax + b)$ ，也就是一个多项式平移。

$\begin{aligned} g(x) &= f(ax + b) \\ &= \sum_{i=0}^n a_i (ax + b)^i \\ &= \sum_{i=0}^n a_i \sum_{j=0}^i \frac {i!} {j! (i-j)!} (ax)^j b^{i-j} \\ &= \sum_{i=0}^n a_i i! \sum_{j=0}^i \frac {(ax)^j} {j!} \cdot \frac {b^{i-j}} {(i-j)!} \\ \end{aligned}$

容易发现可以卷积维护，需要反转一个多项式。

Solution Part4

回到之前的式子，现在我们只需要对于所以 $i \in [1, Q]$ 求出 $g(x^i)$ 的值，相当于求

$\begin{aligned} G(z) &= \sum_{i=0}^Q z^i \sum_{j=0}^n a_j x^{ij} \\ &= \sum_{i=0}^Q z^i \sum_{j=0}^n a_j x^{\binom {i+j} 2 - \binom i 2 - \binom j 2} \\ &= \sum_{i=0}^Q z^i x^{-\binom i2} \sum_{j=0}^n a_j x^{-\binom j2} \cdot x^{\binom {i+j} 2} \\ \end{aligned}$

容易发现可以卷积维护，需要反转一个多项式。

题外话：

一开始想过把 $ij$ 拆成 $\binom i 2 + \binom j 2 - \binom {i-j} 2$ ，结果 $i-j$ 的值域有正有负，反而难处理，尽量还是避免这种情况。

Code

#include<bits/stdc++.h>
#define log(...) (void(0))
// #define log(...) fprintf(stderr,__VA_ARGS__)
const int S=1<<21; char ibuf[S],*iS,*iT;
#define getchar() (iS==iT?(iT=(iS=ibuf)+fread(ibuf,1,S,stdin),(iS==iT?EOF:*iS++)):*iS++)
template<class T> inline void read(T &x){
    x=0; register char c=getchar(); register bool f=0;
    while(!isdigit(c))f^=c=='-',c=getchar();
    while(isdigit(c))x=x*10+c-'0',c=getchar(); if(f)x=-x;
}
#ifdef local
const int N=1e3+10,mod=998244353;
#else
const int N=2.1e6+10,mod=998244353;
#endif
int n,q,rev[N];
struct z {
    int x;
    z(int x=0):x(x){}
    friend inline z operator*(z a,z b){return (long long)a.x*b.x%mod;}
    friend inline z operator-(z a,z b){return (a.x-=b.x)<0?a.x+mod:a.x;}
    friend inline z operator+(z a,z b){return (a.x+=b.x)>=mod?a.x-mod:a.x;}
}q0,qx,qy/*orz qy*/,qc,w[N],f[N],g[N],t1[N],t2[N],t3[N],t4[N],fac[N],ifac[N];
inline z fpow(z a,int b){z s=1;for(;b;b>>=1,a=a*a)if(b&1)s=s*a;return s;}
void initFac(int n){
    fac[0]=ifac[0]=ifac[1]=1;
    for(int i=1;i<=n;i++)fac[i]=fac[i-1]*i;
    for(int i=2;i<=n;i++)ifac[i]=(mod-mod/i)*ifac[mod%i];
    for(int i=2;i<=n;i++)ifac[i]=ifac[i-1]*ifac[i];
}
void initPow(z c,int n,z *a,z *b){
    a[0]=b[0]=1;
    for(int i=1;i<n;i++)a[i]=a[i-1]*c;
    for(int i=1;i<n;i++)b[i]=b[i-1]*a[i-1];
}
int init(int n){
    static int wk=1; int k=-1,l=1; while(l<n)l<<=1,++k;
    for(int i=0;i<l;i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<k);
    for(int &k=wk;k<l;k<<=1){
        z c=fpow(3,(mod-1)/k/2); w[k]=1;
        for(int i=1;i<k;i++)w[i+k]=w[i+k-1]*c;
    }return l;
}
void ntt(z *a,z *b,int l){
    for(int i=0;i<l;i++)b[rev[i]]=a[i];
    for(int k=1;k<l;k<<=1)
        for(int i=0;i<l;i+=(k<<1))
            for(int j=0;j<k;j++){
                z x=b[i+j],y=b[i+j+k]*w[j+k];
                b[i+j]=x+y,b[i+j+k]=x-y;
            }
}
void mul(z *a,z *b,z *c,int n,int m){
    static z t[N];
    int l=init(n+m-1); z v=fpow(l,mod-2);
    ntt(a,c,l),ntt(b,t,l);
    for(int i=0;i<l;i++)t[i]=t[i]*c[i];
    std::reverse(t+1,t+l),ntt(t,c,l);
    for(int i=0;i<l;i++)c[i]=c[i]*v;
}
int main(){
#ifdef local
    freopen("1.in","r",stdin);
#endif
    read(n),read(q),initFac(std::max(n,q));
    for(int i=0;i<=n;i++)read((int&)f[i]);
    read((int&)q0),read((int&)qx),read((int&)qy),qc=(q0*(qx-1)+qy)*fpow(qx-1,mod-2);
    for(int i=0;i<=n;i++)log("%d%c",f[i]," \n"[i==n]);
    const auto move=[&](z *a,z c,z k,int n){
        log("\e[33mMOVE %d %d %d\e[0m\n",c,k,n);
        int i; z t;
        for(i=0,t=1;i<n;i++)t1[i]=a[i]*fac[i],t2[n-i-1]=ifac[i]*t,t=t*c;
        mul(t1,t2,t3,n,n);
        for(i=0,t=1;i<n;i++)a[i]=t3[n-1+i]*ifac[i]*t,t=t*k;
    };
    memset(t1,0,sizeof(t1)),memset(t2,0,sizeof(t2));
    move(f,q0-qc,qc,n+1);
    memset(t1,0,sizeof(t1)),memset(t2,0,sizeof(t2));
    for(int i=0;i<=n;i++)log("%d%c",f[i]," \n"[i==n]);
    const auto czt=[&](z *a,z *b,z c,int n,int q){
        log("\e[33mCZT %d %d %d\e[0m\n",c,n,q);
        initPow(c,n+q,t3,t1);
        std::reverse(t1,t1+n+q);
        initPow(fpow(c,mod-2),std::max(n,q),t3,t4);
        for(int i=0;i<n;i++)t2[i]=a[i]*t4[i];
        mul(t1,t2,t3,n+q,n);
        for(int i=0;i<q;i++)b[i]=t3[n+q-1-i]*t4[i];
    };
    czt(f,g,qx,n+1,q+1);
    for(int i=0;i<=q;i++)log("%d%c",g[i]," \n"[i==q]);
    printf("%d\n",std::accumulate(g+1,g+q+1,0,[](int a,z b){return a^b.x;}));
}