洛谷P5050 【模板】多项式多点求值

https://www.luogu.org/problemnew/show/P5050

给定多项式A(x)，求$A(x_l)$,$A(x_{l+1})$,..,$A(x_r)$

分治：（如果r-l+1=1，直接O(deg(A))暴力求出即可）

首先设$mid=lfloorfrac{l+r}{2} floor$，$P^{[0]}(x)=prod_{i=l}^{mid}(x-x_i)$，$P^{[1]}(x)=prod_{i=mid+1}^{r}(x-x_i)$

以[l,mid]的求值为例：设$A^{[0]}(x)=A(x)\,mod\,P^{[0]}(x)$

即$A(x)=P^{[0]}(x)B^{[0]}(x)+A^{[0]}(x)$（$B^{[0]}$为某个多项式）

可以发现，将$x_l$,$x_{l+1}$,..,$x_{mid}$带入$P^{[0]}(x)$，值都为0

因此对于$l<=i<=mid$，$A(x_i)=A^{[0]}(x_i)$，递归下去算就行；[mid+1,r]的求值同理

这个P可以在分治过程中处理出来

时间复杂度大概是$O(n\,log^2\,n)$（未区分n=r-l+1,m=deg(A)）

版本1：基于版本1，加了小范围暴力，预处理了P方便快速插值

注意：这种分治FFT的题，NTT里面wn需要预处理否则可能慢很多（听说多一个log，没有仔细分析）！

#prag
ma GCC optimize(2)
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<vector>
#include<cmath>
using namespace std;
#define fi first
#define se second
#define mp make_pair
#define pb push_back
typedef long long ll;
typedef unsigned long long ull;
const int md=998244353;
const int N=131072;
#define delto(a,b) ((a)-=(b),((a)<0)&&((a)+=md))
inline int del(int a,int b)
{
    a-=b;
    return a<0?a+md:a;
}
int rev[N];
void init(int len)
{
    int bit=0,i;
    while((1<<(bit+1))<=len)    ++bit;
    for(i=1;i<len;++i)
        rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
}
ull poww(ull a,ull b)
{
    ull ans=1;
    for(;b;b>>=1,a=a*a%md)
        if(b&1)
            ans=ans*a%md;
    return ans;
}
int inv[300011],pw_1[300011],pw_2[300011];
void dft(int *a,int len,int idx)//要求len为2的幂
{
    int i,j,k,t1,t2;ull wn,wnk;
    for(i=0;i<len;++i)
        if(i<rev[i])
            swap(a[i],a[rev[i]]);
    for(i=1;i<len;i<<=1)
    {
        wn=idx==1?pw_1[i]:pw_2[i];
        //wn=poww(idx==1?3:332748118,(md-1)/(i<<1));
        for(j=0;j<len;j+=(i<<1))
        {
            wnk=1;
            for(k=j;k<j+i;++k,wnk=wnk*wn%md)
            {
                t1=a[k];t2=a[k+i]*wnk%md;
                a[k]+=t2;
                (a[k]>=md)&&(a[k]-=md);
                a[k+i]=t1-t2;
                (a[k+i]<0)&&(a[k+i]+=md);
            }
        }
    }
    if(idx==-1)
    {
        ull ilen=inv[len];
        for(i=0;i<len;++i)
            a[i]=a[i]*ilen%md;
    }
}
void p_inv(int *f,int *g,int len)//g=f^(-1);f,g数组的长度不小于2len(需要足够长用于临时存放元素);要求len是2的幂
{
    static int t1[N],t2[N];
    g[0]=poww(f[0],md-2);
    for(int i=2,j;i<=len;i<<=1)
    {
        memcpy(t1,f,sizeof(int)*i);
        memcpy(t2,g,sizeof(int)*(i>>1));
        memset(t2+(i>>1),0,sizeof(int)*(i>>1));
        init(i);
        dft(t1,i,1);dft(t2,i,1);
        for(j=0;j<i;++j)
            t1[j]=ull(t1[j])*t2[j]%md;
        dft(t1,i,-1);
        for(j=0;j<(i>>1);++j)
            t1[j]=t1[j+(i>>1)];
        memset(t1+(i>>1),0,sizeof(int)*(i>>1));
        dft(t1,i,1);
        for(j=0;j<i;++j)
            t1[j]=ull(t1[j])*t2[j]%md;
        dft(t1,i,-1);
        for(j=i>>1;j<i;++j)
            g[j]=md-t1[j-(i>>1)];
    }
}
inline void p_de(int *f,int len)//derivative求导;f=f'
{
    for(int i=0;i<len-1;++i)
        f[i]=ull(i+1)*f[i+1]%md;
    f[len-1]=0;
}
inline void p_in(int *f,int len)//integral积分;f=?f
{
    for(int i=len-1;i>=1;--i)
        f[i]=ull(f[i-1])*inv[i]%md;
    f[0]=0;
}
void p_ln(int *f,int len)//要求len为2的幂，f[0]=1
{
    static int t3[N];
    p_inv(f,t3,len);p_de(f,len);
    init(len<<1);
    dft(f,len<<1,1);dft(t3,len<<1,1);
    for(int i=0;i<(len<<1);++i)
        f[i]=ull(f[i])*t3[i]%md;
    dft(f,len<<1,-1);p_in(f,len);
}
void p_exp(int *f,int *g,int len)//要求len为2的幂，f[0]=0
{
    static int t1[N],t2[N];
    g[0]=1;
    for(int i=2,j;i<=len;i<<=1)
    {
        memcpy(t1,g,sizeof(int)*(i>>1));
        memset(t1+(i>>1),0,sizeof(int)*(i>>1));
        p_ln(t1,i);
        for(j=0;j<(i>>1);++j)
            t1[j]=del(f[j+(i>>1)],t1[j+(i>>1)]);
        memset(t1+(i>>1),0,sizeof(int)*(i>>1));
        init(i);
        dft(t1,i,1);
        memcpy(t2,g,sizeof(int)*(i>>1));
        memset(t2+(i>>1),0,sizeof(int)*(i>>1));
        dft(t2,i,1);
        for(j=0;j<i;++j)
            t1[j]=ull(t1[j])*t2[j]%md;
        dft(t1,i,-1);
        for(j=i>>1;j<i;++j)
            g[j]=t1[j-(i>>1)];
    }
}
void p_div(int *a,int *b,int *c,int n,int m)//c=a/b;deg(a)=n,deg(b)=m,deg(c)=n-m;a,b无前导0;n>=m
{
    reverse(a,a+n+1);reverse(b,b+m+1);
    int x=n-m+1,t=1;
    for(;t<x;t<<=1);
    memset(b+m+1,0,sizeof(int)*max(t-m-1,0));
    p_inv(b,c,t);
    memset(c+x,0,sizeof(int)*((t<<1)-x));
    memset(a+x,0,sizeof(int)*((t<<1)-x));
    init(t<<1);
    dft(a,t<<1,1);dft(c,t<<1,1);
    for(int i=0;i<(t<<1);++i)
        c[i]=ull(c[i])*a[i]%md;
    dft(c,t<<1,-1);
    memset(c+(n-m+1),0,sizeof(int)*((t<<1)-n+m-1));
    reverse(c,c+x);
}
void p_divmod(int *a,int *b,int *c,int *d,int n,int m)//c=a/b,d=a%b,deg(d)=(<=)m-1;其余同上
{
    static int t1[N];
    memcpy(d,a,sizeof(int)*(m+1));
    int x=n+1,t=1;
    for(;t<x;t<<=1);
    memcpy(t1,b,sizeof(int)*(m+1));
    memset(t1+m+1,0,sizeof(int)*max(t-m-1,0));
    p_div(a,b,c,n,m);
    memcpy(a,c,sizeof(int)*(n-m+1));
    memset(a+n-m+1,0,sizeof(int)*(t-n+m-1));
    init(t);
    dft(a,t,1);dft(t1,t,1);
    for(int i=0;i<t;++i)
        t1[i]=ull(t1[i])*a[i]%md;
    dft(t1,t,-1);
    for(int i=0;i<=m;++i)
        delto(d[i],t1[i]);
}
namespace P_me
{
    int *ta[N];//用线段树的方法给递归的每一层一个编号，ta[i]表示编号为i的层的P函数的各项系数
    int data[N*40],*tp;//内存池
    int *a,*x,*y;
#define LC (u<<1)
#define RC (u<<1|1)
    int mt1[N];
    const int T=200;//小范围暴力阀值
    void _p_me1(int l,int r,int u)//计算(x-x_l)(x-x_{l+1})..(x-x_r)并存下来
    {
        if(r-l<=T)
        {
            int i,j;
            tp[0]=1;
            for(i=l;i<=r;++i)
            {
                tp[i-l+1]=tp[i-l];
                for(j=i-l;j>=1;--j)
                {
                    tp[j]=(ull(tp[j])*(md-x[i])+tp[j-1])%md;
                }
                tp[0]=ull(tp[0])*(md-x[i])%md;
            }
            ta[u]=tp;tp+=r-l+2;
            return;
        }
        int mid=(l+r)>>1;
        _p_me1(l,mid,LC);_p_me1(mid+1,r,RC);
        int x=r-l+2,t=1;//x=(mid-l+1)+(r-mid)+1
        for(;t<x;t<<=1);
        memcpy(mt1,ta[LC],sizeof(int)*(mid-l+2));
        memset(mt1+mid-l+2,0,sizeof(int)*(t-mid+l-2));
        memcpy(tp,ta[RC],sizeof(int)*(r-mid+1));
        memset(tp+r-mid+1,0,sizeof(int)*(t-r+mid-1));
        init(t);
        dft(mt1,t,1);dft(tp,t,1);
        for(int i=0;i<t;++i)
            tp[i]=ull(tp[i])*mt1[i]%md;
        dft(tp,t,-1);
        ta[u]=tp;tp+=r-l+2;
    }
    int mt2[N],mt3[N];
    void _p_me2(int *a,int n,int l,int r,int u)//a是A的系数,deg(A)<=n；求A(x_l)到A(x_r),放入y_l到y_r
    {
        if(r-l<=T)
        {
            int t,i,j;
            for(i=l;i<=r;++i)
            {
                t=a[n];
                for(j=n-1;j>=0;--j)
                    t=(ull(t)*x[i]+a[j])%md;
                y[i]=t;
            }
            return;
        }
        int x=(n+1)<<1,t=1;
        for(;t<x;t<<=1);
        int mt4[t];//根据需要改成new?
        int mid=(l+r)>>1,n1;
        memcpy(mt1,a,sizeof(int)*(n+1));
        for(n1=n;n1>=0 && mt1[n1]==0;)    --n1;
        if(n1<0)
        {
            memset(y+l,0,sizeof(int)*(r-l+1));
            return;
        }
        memcpy(mt2,ta[LC],sizeof(int)*(mid-l+2));
        if(n1<mid-l+1)
        {
            memcpy(mt4,mt1,sizeof(int)*(n1+1));
            _p_me2(mt4,n1,l,mid,LC);
        }
        else
        {
            p_divmod(mt1,mt2,mt3,mt4,n1,mid-l+1);
            _p_me2(mt4,mid-l,l,mid,LC);    
        }
        memcpy(mt1,a,sizeof(int)*(n+1));
        for(n1=n;n1>=0 && mt1[n1]==0;)    --n1;
        memcpy(mt2,ta[RC],sizeof(int)*(r-mid+1));
        if(n1<r-mid)
        {
            memcpy(mt4,mt1,sizeof(int)*(n1+1));
            _p_me2(mt4,n1,mid+1,r,RC);
        }
        else
        {
            p_divmod(mt1,mt2,mt3,mt4,n1,r-mid);
            _p_me2(mt4,r-mid-1,mid+1,r,RC);
        }
    }
    void p_multieval(int *a0,int *x0,int *y0,int n,int m)//deg(a)=n,x有m个数
    {
        tp=data;
        a=a0;x=x0;y=y0;
        _p_me1(0,m-1,1);
        _p_me2(a,n,0,m-1,1);
    }
}
using P_me::p_multieval;
int a[N],x[N],y[N];
int n,m;
int main()
{
    int i;
    inv[1]=1;
    for(i=2;i<=300000;++i)
        inv[i]=ull(md-md/i)*inv[md%i]%md;
    for(i=1;i<300000;i<<=1)
    {
        pw_1[i]=poww(3,(md-1)/(i<<1));
        pw_2[i]=poww(332748118,(md-1)/(i<<1));
    }
    //n=100000;m=100000;
    scanf("%d%d",&n,&m);
    for(i=0;i<=n;++i)
        //a[i]=rand()%md;
        scanf("%d",a+i);
    for(i=0;i<m;++i)
        //x[i]=rand()%md;
        scanf("%d",x+i);
    p_multieval(a,x,y,n,m);
    for(i=0;i<m;++i)
        printf("%d
",y[i]);
    return 0;
}