P5050 【模板】多项式多点求值 题解

以下题解仅供学习参考使用。

抄袭、复制题解,以达到刷AC率/AC数量或其他目的的行为,在洛谷是严格禁止的。

洛谷非常重视学术诚信。此类行为将会导致您成为作弊者。具体细则请查看洛谷社区规则

评论

  • GNAQ
    什么毒瘤玩意
  • 子曰子悦
    哈哈哈哈哈哈这题解有毒
  • sermoon
    仿佛预感到了插值的惨烈
  • 142857cs
    插值都出到1e5了,这个应该也可以?
  • 142857cs
    插值应该没那么惨吧,O(n^2)的插值有多少人会?
  • ButterflyDew
    哈哈哈哈哈额
  • rEdWhitE_uMbrElla
    毒瘤。。。
  • 南城忆潇湘
    原谅我笑了出来
  • olinr
    卡常巨佬啊%%%%
  • olinr
    666666
作者: 玫葵之蝶 更新时间: 2018-12-17 23:21  在Ta的博客查看 举报    64  

脑子想想就知道这种两个log还常数大的题可以暴力艹啊

显然我们可以有一个O(nm)的秦九韶暴力,这个不会你就可以退役了

开始卡常:

1.首先加register,能加的都加,出了存系数的那个不加

2.写快读快写,这个也不必非写fread,我就没写

3.循环展开,实测4层最优

4.最后开个unsigned long long,你循环展开的过程量就不用取模了

5.写个unordered_map,如果已经算过了,就直接输出(这个貌似没卵用)

然后找个夜深人静的好时候,提交就好了,你就可以AC了

#pragma GCC optimize("Ofast")
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<cmath>
#include<iostream>
#include<algorithm>
#include<unordered_map>
#define LL unsigned long long
#define R register
using namespace std;
inline void read(int &x){
    x=0;int f=1;char ch=getchar();
    while(ch<'0' || ch>'9'){if(ch=='-')f=-1;ch=getchar();}
    while(ch>='0' && ch<='9')x=(x<<3)+(x<<1)+(ch^48),ch=getchar();
    if(f==-1)x=-x;
}
char s[20];int cnt;
inline void write(int x){
    cnt=0;
    while(x){
        s[++cnt]=x%10+'0';
        x/=10;
    }
    for(R int i=cnt;i;--i)putchar(s[i]);
    putchar('\n');
}
const int N=64005,mod=998244353;
unordered_map<int,int> mp;
int a[N];
int main(){
    R int n,m,i,j,x;
    R LL b[17],c1,c2,c3,c4,now;
    read(n);read(m);
    for(i=0;i<=n;++i)read(a[i]);
    for(i=1;i<=m;++i){
        read(x);
        if(mp[x]){write(mp[x]);continue;}
        b[0]=1;
        for(j=1;j<=16;++j)b[j]=b[j-1]*x%mod;
        now=a[n];
        for(j=n-1;j-15>=0;j-=16){
            c1=now*b[16]+a[j]*b[15]+a[j-1]*b[14]+a[j-2]*b[13];
            c2=a[j-3]*b[12]+a[j-4]*b[11]+a[j-5]*b[10]+a[j-6]*b[9];
            c3=a[j-7]*b[8]+a[j-8]*b[7]+a[j-9]*b[6]+a[j-10]*b[5];
            c4=a[j-11]*b[4]+a[j-12]*b[3]+a[j-13]*b[2]+a[j-14]*b[1];
            now=(c1+c2+c3+c4+a[j-15])%mod;
        }
        for(j=n%16-1;~j;--j)now=(now*x+a[j])%mod;
        write(mp[x]=now);
    }
    return 0;
}

评论

  • 皎月半洒花
    一开始构造的多项式是m次的吧
作者: mrsrz  珂愛 更新时间: 2018-12-26 18:44  在Ta的博客查看 举报    19  

我们将要求值的点均分成两份,构造多项式$P_0(x)=\prod\limits_{i=1}^{\lfloor\frac n 2\rfloor}(x-x_i)$,$P_1(x)=\prod\limits_{i=\lfloor\frac n 2\rfloor+1}^{n}(x-x_i)$。

显然,对于$\forall i\in[1,\lfloor\frac n 2 \rfloor]$,有$P_0(x_i)=0$。$P_1$同理。

我们假设多项式$D(x),R(x)$满足:$D(x)$是一个$n-\lfloor\frac n 2\rfloor$次多项式,$R(x)$是一个次数不超过$\lfloor\frac n 2\rfloor-1$的多项式,且$A(x)=P_0(x)D(x)+R(x)$。

那么对于$\forall i\in[1,\lfloor\frac n 2 \rfloor]$,有$A(x_i)=R(x_i)$。$P_1$同理可得。

由于$R(x)$的次数是$A(x)$的次数的一半,所以我们把一个规模为$n$的问题,用$O(n\log n)$的复杂度(多项式取模,详见多项式除法),转化为两个规模为$\frac n 2$的问题。

分治计算即可,由主定理得时间复杂度$O(n\log^2 n)$。

求每次的$P_0(x),P_1(x)$,可以开始用一次分治NTT预处理。此处时间复杂度也为$O(n\log^2 n)$。

常数极大就是了QAQ。

Code:

#include<cstdio>
#include<algorithm>
typedef long long LL;
const int md=998244353,N=262145;
//Poly begin
int rev[N],lim,g[20][N],M,vv;
inline void upd(int&a){a+=a>>31&md;}
inline int pow(int a,int b){
    int ret=1;
    for(;b;b>>=1,a=(LL)a*a%md)if(b&1)ret=(LL)ret*a%md;return ret;
}
inline void init(int n){
    int l=-1;
    for(lim=1;lim<n;lim<<=1)++l;M=l+1;
    for(int i=1;i<lim;++i)rev[i]=(rev[i>>1]>>1)|((i&1)<<l);vv=pow(lim,md-2);
}
void NTT(int*a,int f){
    for(int i=1;i<lim;++i)if(i<rev[i])std::swap(a[i],a[rev[i]]);
    for(int i=0;i<M;++i){
        const int*G=g[i],c=1<<i;
        for(int j=0;j<lim;j+=c<<1)
        for(int k=0;k<c;++k){
            const int x=a[j+k],y=a[j+k+c]*(LL)G[k]%md;
            upd(a[j+k]+=y-md),upd(a[j+k+c]=x-y);
        }
    }
    if(!f){
        for(int i=0;i<lim;++i)a[i]=(LL)a[i]*vv%md;
        std::reverse(a+1,a+lim);
    }
}
void INV(const int*a,int*B,int n){
    if(n==1){
        *B=pow(*a,md-2);
        return;
    }
    INV(a,B,n+1>>1);
    init(n<<1);
    static int A[N];
    for(int i=0;i<n;++i)A[i]=a[i];
    for(int i=n;i<lim;++i)A[i]=0;
    NTT(A,1),NTT(B,1);
    for(int i=0;i<lim;++i)B[i]=B[i]*(2-(LL)A[i]*B[i]%md+md)%md;
    NTT(B,0);
    for(int i=n;i<lim;++i)B[i]=0;
}
void REV(int*A,int n){std::reverse(A,A+n+1);}
void MOD(const int*a,const int*b,int*R,int n,int m){
    static int A[N],B[N],D[N];
    for(int i=0;i<n<<1;++i)D[i]=0;
    for(int i=0;i<=m;++i)B[i]=b[i];
    REV(B,m);
    for(int i=n-m+1;i<n<<1;++i)B[i]=0;
    INV(B,D,n-m+1);
    init(n-m+1<<1);
    for(int i=0;i<=n-m+1;++i)A[i]=a[n-i];
    for(int i=n-m+2;i<lim;++i)A[i]=0;
    NTT(A,1),NTT(D,1);
    for(int i=0;i<lim;++i)D[i]=(LL)D[i]*A[i]%md;
    NTT(D,0);
    REV(D,n-m);
    init(n+1<<1);
    for(int i=n-m+1;i<lim;++i)D[i]=0;
    for(int i=0;i<lim;++i)A[i]=B[i]=0;
    for(int i=0;i<=m;++i)B[i]=b[i];
    for(int i=0;i<=n;++i)A[i]=a[i];
    NTT(B,1),NTT(D,1);
    for(int i=0;i<lim;++i)B[i]=(LL)B[i]*D[i]%md;
    NTT(B,0);
    for(int i=0;i<m;++i)upd(R[i]=A[i]-B[i]);
}
//Poly end
int n,m,A[N],a[N],*P[N],len[N];
void solve(int l,int r,int o){
    if(l==r){
        len[o]=1;
        P[o]=new int[2];
        upd(P[o][0]=-a[l]),P[o][1]=1;
        return;
    }
    const int mid=l+r>>1,L=o<<1,R=L|1;
    solve(l,mid,L),solve(mid+1,r,R);
    len[o]=len[L]+len[R];
    P[o]=new int[len[o]+1];
    init(len[o]+1<<1);
    static int A[N],B[N];
    for(int i=len[L];~i;--i)A[i]=P[L][i];
    for(int i=len[L]+1;i<lim;++i)A[i]=0;
    for(int i=len[R];~i;--i)B[i]=P[R][i];
    for(int i=len[R]+1;i<lim;++i)B[i]=0;
    NTT(A,1),NTT(B,1);
    for(int i=0;i<lim;++i)A[i]=(LL)A[i]*B[i]%md;
    NTT(A,0);
    for(int i=len[o];~i;--i)P[o][i]=A[i];
}
void solve(int l,int r,int o,const int*A){
    if(l==r){printf("%d\n",*A);return;}
    const int mid=l+r>>1,L=o<<1,R=L|1;
    int B[len[o]+2<<1];
    MOD(A,P[L],B,len[o]-1,len[L]);
    solve(l,mid,L,B);
    MOD(A,P[R],B,len[o]-1,len[R]);
    solve(mid+1,r,R,B);
}
int main(){
    for(int i=0;i<19;++i){
        int*G=g[i];
        G[0]=1;
        const int gi=G[1]=pow(3,(md-1)/(1<<i+1));
        for(int j=2;j<1<<i;++j)G[j]=(LL)G[j-1]*gi%md;
    }
    scanf("%d%d",&n,&m);if(!m)return 0;
    for(int i=0;i<=n;++i)scanf("%d",A+i);
    for(int i=1;i<=m;++i)scanf("%d",a+i);
    solve(1,m,1);
    if(n>=m)MOD(A,P[1],A,n,m);
    solve(1,m,1,A);
    return 0;
}

评论

  • Jμdge
    还人傻常熟大,每次用时都是我的 $\frac{1}{n}$ ,$n >> 1$
作者: bztMinamoto 更新时间: 2019-02-13 19:52  在Ta的博客查看 举报    7  

传送门

人傻常数大.jpg

因为求逆的时候没清零结果调了几个小时……

前置芝士

多项式除法,多项式求逆

什么?你不会?左转你谷模板区,包教包会

题解

首先我们要知道一个结论$$f(x_0)\equiv f(x)\pmod{(x-x_0)}$$

其中$x_0$为一个常量,$f(x_0)$也为一个常量

证明如下,设$f(x)=g(x)(x-x_0)+A$,也就是说$A$是$f(x)$对$(x-x_0)$这个多项式取模之后的结果

因为$(x-x_0)$的最高次项为$1$,所以$A$的最高次项为$0$,也就是说$A$是一个常数,即$f(x)\equiv A\pmod{(x-x_0)}$

我们把$x_0$代入上式,得$f(x_0)=g(x_0)(x_0-x_0)+A$,同理可得$f(x_0)\equiv A\pmod{(x-x_0)}$

于是我们知道上式成立

这有毛用啊$O(n\log n)$多项式取模还没我暴力快

乍一看的确没啥卵用,但是考虑取模的过程是否能优化呢?

答案是可以的,我们考虑分治。设当前分治区间为$[l,r]$,令$P_0=\prod_{i=l}^{mid}(x-x_i)$,$P_1=\prod_{i=mid+1}^r (x-x_i)$,当前已经算出了$A\equiv f(x)\pmod{\prod_{i=l}^r(x-x_i)}$,那么只要分别用$A$对$P_0$和$P_1$取模,然后继续递归下去就行了。取模之后$A(x)$的最高次项的次数变为原来的一半,问题规模也就变为原来的一半。继续递归下去就行了

时间复杂度为$O(n\log^2n)$

upd:改了改代码,常数应该会小一点,比方说分治到某个时候暴力秦九韶展开

//minamoto
#include<bits/stdc++.h>
#define R register
#define fp(i,a,b) for(R int i=(a),I=(b)+1;i<I;++i)
#define fd(i,a,b) for(R int i=(a),I=(b)-1;i>I;--i)
#define go(u) for(int i=head[u],v=e[i].v;i;i=e[i].nx,v=e[i].v)
using namespace std;
char buf[1<<21],*p1=buf,*p2=buf;
inline char getc(){return p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<21,stdin),p1==p2)?EOF:*p1++;}
int read(){
    R int res,f=1;R char ch;
    while((ch=getc())>'9'||ch<'0')(ch=='-')&&(f=-1);
    for(res=ch-'0';(ch=getc())>='0'&&ch<='9';res=res*10+ch-'0');
    return res*f;
}
char sr[1<<21],z[20];int C=-1,Z=0;
inline void Ot(){fwrite(sr,1,C+1,stdout),C=-1;}
void print(R int x){
    if(C>1<<20)Ot();if(x<0)sr[++C]='-',x=-x;
    while(z[++Z]=x%10+48,x/=10);
    while(sr[++C]=z[Z],--Z);sr[++C]='\n';
}
const int N=(1<<17)+5,P=998244353;
inline int add(R int x,R int y){return x+y>=P?x+y-P:x+y;}
inline int dec(R int x,R int y){return x-y<0?x-y+P:x-y;}
inline int mul(R int x,R int y){return 1ll*x*y-1ll*x*y/P*P;}
int ksm(R int x,R int y){
    R int res=1;
    for(;y;y>>=1,x=mul(x,x))(y&1)?res=mul(res,x):0;
    return res;
}
int r[19][N],w[2][N],lg[N],inv[19];
void Pre(){
    fp(d,1,17){
        fp(i,1,(1<<d)-1)r[d][i]=(r[d][i>>1]>>1)|((i&1)<<(d-1));
        lg[1<<d]=d,inv[d]=ksm(1<<d,P-2);
    }
    for(R int t=(P-1)>>1,i=1,x,y;i<131072;i<<=1,t>>=1){
        x=ksm(3,t),y=ksm(332748118,t),w[0][i]=w[1][i]=1;
        fp(k,1,i-1)
            w[1][k+i]=mul(w[1][k+i-1],x),
            w[0][k+i]=mul(w[0][k+i-1],y);
    }
}
int lim,d,n,m;
inline void init(R int len){lim=1,d=0;while(lim<len)lim<<=1,++d;}
void NTT(int *A,int ty){
    fp(i,0,lim-1)if(i<r[d][i])swap(A[i],A[r[d][i]]);
    for(R int mid=1;mid<lim;mid<<=1)
        for(R int j=0,t;j<lim;j+=(mid<<1))
            fp(k,0,mid-1)
                A[j+k+mid]=dec(A[j+k],t=mul(w[ty][mid+k],A[j+k+mid])),
                A[j+k]=add(A[j+k],t);
    if(!ty)fp(i,0,lim-1)A[i]=mul(A[i],inv[d]);
}
void Inv(int *a,int *b,int len){
    if(len==1)return b[0]=ksm(a[0],P-2),void();
    Inv(a,b,len>>1),lim=(len<<1),d=lg[lim];
    static int A[N],B[N];
    fp(i,0,len-1)A[i]=a[i],B[i]=b[i];fp(i,len,lim-1)A[i]=B[i]=0;
    NTT(A,1),NTT(B,1);
    fp(i,0,lim-1)A[i]=mul(A[i],mul(B[i],B[i]));
    NTT(A,0);
    fp(i,0,len-1)b[i]=dec(add(b[i],b[i]),A[i]);
    fp(i,len,lim-1)b[i]=0;
}
struct node{
    node *lc,*rc;vector<int>vec;int deg;
    void Mod(const int *a,int *r,int n){
        static int A[N],B[N],D[N];
        int len=1;while(len<=n-deg)len<<=1;
        fp(i,0,n)A[i]=a[n-i];fp(i,0,deg)B[i]=vec[deg-i];
        fp(i,n-deg+1,len-1)B[i]=0;
        Inv(B,D,len);
        lim=(len<<1),d=lg[lim];
        fp(i,n-deg+1,lim-1)A[i]=D[i]=0;
        NTT(A,1),NTT(D,1);
        fp(i,0,lim-1)A[i]=mul(A[i],D[i]);
        NTT(A,0);
        reverse(A,A+n-deg+1);
        init(n+1);
        fp(i,n-deg+1,lim-1)A[i]=0;
        fp(i,0,deg)B[i]=vec[i];fp(i,deg+1,lim-1)B[i]=0;
        NTT(A,1),NTT(B,1);
        fp(i,0,lim-1)A[i]=mul(A[i],B[i]);
        NTT(A,0);
        fp(i,0,deg-1)r[i]=dec(a[i],A[i]);
    }
    void Mul(){
        static int A[N],B[N];deg=lc->deg+rc->deg,vec.resize(deg+1),init(deg+1);
        fp(i,0,lc->deg)A[i]=lc->vec[i];fp(i,lc->deg+1,lim-1)A[i]=0;
        fp(i,0,rc->deg)B[i]=rc->vec[i];fp(i,rc->deg+1,lim-1)B[i]=0;
        NTT(A,1),NTT(B,1);
        fp(i,0,lim-1)A[i]=mul(A[i],B[i]);
        NTT(A,0);
        fp(i,0,deg)vec[i]=A[i];
    }
}pool[N],*rt;
int A[N],a[N],tot;
inline node* newnode(){return &pool[tot++];}
void solve(node* &p,int l,int r){
    p=newnode();
    if(l==r)return p->deg=1,p->vec.resize(2),p->vec[0]=P-a[l],p->vec[1]=1,void();
    int mid=(l+r)>>1;
    solve(p->lc,l,mid),solve(p->rc,mid+1,r);
    p->Mul();
}
int b[25];
void calc(node* p,int l,int r,const int *A){
    if(r-l<=512){
        fp(i,l,r){
            int x=a[i],c1,c2,c3,c4,now=A[r-l];
            b[0]=1;fp(j,1,16)b[j]=mul(b[j-1],x);
            for(R int j=r-l-1;j-15>=0;j-=16){
                c1=(1ll*now*b[16]+1ll*A[j]*b[15]+1ll*A[j-1]*b[14]+1ll*A[j-2]*b[13])%P,
                c2=(1ll*A[j-3]*b[12]+1ll*A[j-4]*b[11]+1ll*A[j-5]*b[10]+1ll*A[j-6]*b[9])%P,
                c3=(1ll*A[j-7]*b[8]+1ll*A[j-8]*b[7]+1ll*A[j-9]*b[6]+1ll*A[j-10]*b[5])%P,
                c4=(1ll*A[j-11]*b[4]+1ll*A[j-12]*b[3]+1ll*A[j-13]*b[2]+1ll*A[j-14]*b[1])%P,
                now=(0ll+c1+c2+c3+c4+A[j-15])%P;
            }
            fd(j,(r-l)%16-1,0)now=(1ll*now*x+A[j])%P;
            print(now);
        }
        return;
    }
    int mid=(l+r)>>1,b[p->deg+1];
    p->lc->Mod(A,b,p->deg-1),calc(p->lc,l,mid,b);
    p->rc->Mod(A,b,p->deg-1),calc(p->rc,mid+1,r,b);
}
int main(){
//  freopen("testdata.in","r",stdin);
    n=read(),m=read();if(!m)return 0;
    Pre();
    fp(i,0,n)A[i]=read();
    fp(i,1,m)a[i]=read();
    solve(rt,1,m);
    if(n>=m)rt->Mod(A,A,n);
    calc(rt,1,m,A);
    return Ot(),0;
}

评论

  • p_b_p_b
    orz
  • p_b_p_b
    看这惨烈的代码,仿佛您把能错的地方都错了一遍
  • AThousandSuns
    @p_b_p_b 的确是这样……加了注释的地方都错过……
  • David_Sktain_Li
    至少您没有在读入的时候把m写成n/绝望
  • JohnVictor
    @David_Sktain_Li 至少您没有写过scanf("&d",&m)然后就RE了
作者: AThousandSuns  更新时间: 2019-01-30 12:03  在Ta的博客查看 举报    4  

调了一个月终于调出来了……

这里就加点详细的注释表明哪里容易写错好了。毕竟没题解调这种题会死人的……


我们发现如果构造一个多项式 $P_{l,r}(x)=\prod\limits^r_{i=l}(x-a_i)$,那么对于所有的 $l\le i\le r$,都有 $P_{l,r}(a_i)=0$。

存在一个 $n-(r-l+1)$ 次多项式 $Q$,和 $r-l$ 次多项式 $R$ 满足 $F(a_i)=P_{l,r}(a_i)Q(a_i)+R(a_i)$。因为 $P_{l,r}(a_i)=0$,所以 $F(a_i)=R(a_i)$。我们便成功把一个 $n$ 次多项式变成了 $r-l$ 次多项式。

求 $R(a_i)$?发现它实际上是 $F$ 模 $P_{l,r}$,多项式除法/取模!

现在考虑分治。令 $F_{l,r}(x)$ 为当前区间的多项式,那么分成 $[l,mid]$ 和 $[mid+1,r]$ 计算。以 $[l,mid]$ 为例,$F_{l,r}(a_i)=P_{l,mid}(a_i)Q(a_i)+F_{l,mid}(a_i)$。

边界:$l=r$ 时 $F$ 次数为 $0$,那么 $F(a_i)=[x^0]F_{l,r}(x)$。

但是 $P_{l,r}(x)$ 直接求也不可承受。实际上我们可以另外再来一次分治,$P_{l,r}(x)=P_{l,mid}(x)P_{mid+1,r}(x)$。边界也是 $l=r$。

时间复杂度:$T(n)=2T(\frac{n}{2})+O(n\log n)=O(n\log^2n)$。

因为自带大常数,我就模仿了@officeyutong 的做法,在 $r-l\le100$ 的时候直接秦九韶爆算。实际上会跑得快很多。

下面是我的代码+注释:(我用vector写的,可能会比较丑)

#include<bits/stdc++.h>
using namespace std;
const int mod=998244353,g=3;
#define lson o<<1,l,mid
#define rson o<<1|1,mid+1,r
#define FOR(i,a,b) for(int i=(a);i<=(b);i++)
#define ROF(i,a,b) for(int i=(a);i>=(b);i--)
inline int read(){
    int x=0,f=0;char ch=getchar();
    while(ch<'0' || ch>'9') f|=ch=='-',ch=getchar();
    while(ch>='0' && ch<='9') x=x*10+ch-'0',ch=getchar();
    return f?-x:x;
}
int n,m,b[64444],ans[64444];
int limit,l,rev[266666];
vector<int> poly[266666],a;
inline void init(int upr){
    for(limit=1,l=0;limit<=upr;limit<<=1,l++);
    FOR(i,0,limit-1) rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-1));
}
inline int add(int x,int y){return x+y>=mod?x+y-mod:x+y;}
inline int sub(int x,int y){return x<y?x-y+mod:x-y;}
inline int qpow(int a,int b){
    int ans=1;
    for(;b;b>>=1,a=1ll*a*a%mod) if(b&1) ans=1ll*ans*a%mod;
    return ans;
}
void NTT(vector<int> &A,int tp){
    while(A.size()<limit) A.push_back(0);   //加上这话!不然会越界!
    FOR(i,0,limit-1) if(i<rev[i]) swap(A[i],A[rev[i]]);
    for(int i=1;i<limit;i<<=1)
        for(int j=0,r=i<<1,Wn=qpow(g,mod-1+tp*(mod-1)/r);j<limit;j+=r)
            for(int k=0,w=1;k<i;k++,w=1ll*w*Wn%mod){
                int x=A[j+k],y=1ll*A[i+j+k]*w%mod;
                A[j+k]=add(x,y);A[i+j+k]=sub(x,y);
            }
    if(tp==-1){
        int linv=qpow(limit,mod-2);
        FOR(i,0,limit-1) A[i]=1ll*A[i]*linv%mod;
    }
}
void NTT(vector<int> A,vector<int> B,vector<int> &C){
    C.clear();  //清空!
    NTT(A,1);NTT(B,1);
    FOR(i,0,limit-1) C.push_back(1ll*A[i]*B[i]%mod);    //不能直接用下标!
    NTT(C,-1);
}
void inv(vector<int> A,vector<int> &B,int deg){ //逆元(必须保证调用时B是空的)
    if(deg==1) return B.push_back(qpow(A[0],mod-2));
    inv(A,B,(deg+1)>>1);
    init(deg<<1);
    while(A.size()<limit) A.push_back(0);   //虽然NTT内有补全到limit的话,但这也要写
    FOR(i,deg,limit-1) A[i]=0;  //因为这里需要把deg后的设为0
    NTT(A,1);NTT(B,1);
    FOR(i,0,limit-1) B[i]=1ll*sub(2,1ll*A[i]*B[i]%mod)*B[i]%mod;
    NTT(B,-1);
    FOR(i,deg,limit-1) B[i]=0;
}
void modulo(vector<int> A,vector<int> B,vector<int> &D,int n,int m){
    //取模(须保证B高位为空,A随便)
    while(A.size()<=n) A.push_back(0);
    while(B.size()<=m) B.push_back(0);  //把空位补全!(具体原因下文会说)
    if(n<m) return void(D=A);   //如果除不了,那么答案就是被除数!
    vector<int> Arev,Brev,Brevinv,Crev,C;D.clear();
    FOR(i,0,n) Arev.push_back(A[n-i]);
    FOR(i,0,m) Brev.push_back(B[m-i]);
    FOR(i,n-m+2,Arev.size()-1) Arev[i]=0;
    FOR(i,n-m+1,Brev.size()-1) Brev[i]=0;   //记得清零
    inv(Brev,Brevinv,n-m+1);
    init((n-m+1)<<1);
    NTT(Arev,Brevinv,Crev);
    FOR(i,0,n-m) C.push_back(Crev[n-m-i]);
    init(n<<1);
    NTT(B,C,D);
    FOR(i,0,m-1) D[i]=sub(A[i],D[i]);
    FOR(i,m,limit-1) D[i]=0;
}
void get_poly(int o,int l,int r){
    if(l==r) return void((poly[o].push_back(b[l]?mod-b[l]:0),poly[o].push_back(1)));
    //l=r时就是x-a_l
    int mid=(l+r)>>1;
    get_poly(lson);get_poly(rson);  //分治两边
    init(r-l+1);
    NTT(poly[o<<1],poly[o<<1|1],poly[o]);   //乘到自己这
}
void solve(vector<int> A,int o,int l,int r){
    if(r-l<=100){   //常数(?)优化
        FOR(i,l,r){
            int s=0;
            ROF(j,A.size()-1,0) s=add(1ll*s*b[i]%mod,A[j]); //直接秦九韶
            ans[i]=s;
        }
        return;
    }
    vector<int> B;int mid=(l+r)>>1;
    modulo(A,poly[o<<1],B,r-l,mid-l+1); //取模到左边
    //此时A的次数应该<=r-l
    solve(B,lson);
    modulo(A,poly[o<<1|1],B,r-l,r-mid); //取模到右边
    //因为A一开始的次数(就是n)可能巨小无比,甚至<r-l,所以在取模中要把高位补齐,不然可能越界
    solve(B,rson);
}
int main(){
    n=read();m=read();
    FOR(i,0,n) a.push_back(read());
    FOR(i,1,m) b[i]=read();
    get_poly(1,1,m);
    modulo(a,poly[1],a,n,m);    //提前取模(注意!就是这里可能n<m,所以在取模中就要特判)
    solve(a,1,1,m);
    FOR(i,1,m) printf("%d\n",ans[i]);
}

评论

  • 还没有评论
作者: Fuyuki  更新时间: 2020-02-04 22:14  在Ta的博客查看 举报    3  

记待求值的多项式为 $F(x)$,需要求出其在 $x_0$ 处的点值 $F(x_0)$。

如果 $F(x)=Q(x)G(x)+R(x)$,则有 $F(x_0)=Q(x_0)G(x_0)+R(x_0)$。 如果构造 $G(x)$ 使得 $G(x_0)=0$,那么就有 $F(x_0)=R(x_0)$。

那么构造 $G(x)=x-x_0$,则 $F(x)$ 对 $G(x)$ 做多项式取模后得到的常数项就是 $F(x_0)$(因为取模后只有这一项了)。

多项式取模的过程可以分治进行,将 $F(x)$ 对 $\prod_{i=1}^{m/2} (x-x_i)$ 和 $\prod_{i=m/2+1}^{m} (x-x_i)$ 分别取模,再将取模后的结果向两边递归下去。每次 $O(mlogm)$ 的多项式取模都可以使得原问题的规模减半,就可以在 $O(nlogn+mlog^2m)$ 的时间内求出所有的点值(最开始还需要将 $F(x)$ 对 $\prod_{i=1}^{m}(x-x_i)$ 取模)。

这是大部分人的做法,但是鉴于多项式取模的大常数,这个做法的效率并不高。

重新分析一下,$R(x)$ 只有常数项,如果能算出 $Q(x)$ 的常数项,那么就可以计算出 $R(x)$。

考虑一下多项式除法是进行了一个怎么样的过程,假设 $F(x),G(x)$ 的最高次项次数分别为 $n,m$:

$$F(x)=Q(x)G(x)+R(x)$$

$$x^nF(\frac{1}{x})=x^{n-m}Q(\frac{1}{x})x^mG(\frac{1}{x})+x^nR(\frac{1}{x})$$

$$F_r(x)=Q_r(x)G_r(x)+x^{n-m+1}R_r(x)$$

$$Q_r(x)=F_r(x)G_r^{-1}(x)\pmod{x^{n-m+1}}$$

在这里,$F_r(x)$ 表示的是 $F_r(x)$ 的系数翻转后的结果。

设 $G(x)=\prod_{i=1}^{m}(x-x_i),G_0(x)=\prod_{i=1}^{m/2}(x-x_i),G_1(x)=\prod_{i=m/2+1}^{m}(x-x_i)$ 。

可以想到先求出 $Q_r(x)=F_r(x)G_r^{-1}(x)$ 然后用利用 $G_{0r}^{-1}(x)=G_r^{-1}(x)G_1(x)$ 来向下递归。

最底层的 $Q(x)$ 有 $n$ 项,也就是说为了保证底层计算结果是正确的,所有的计算都必须至少在 $\pmod{x^n}$ 的意义下进行,这样的复杂度显然是无法接受的。

但是注意到根本不需要计算出 $Q(x)$,只需要计算出 $Q(x)$ 的常数项,或者说 $[x^{n-1}]Q_r(x)$。换言之,只有某一项可能对最高项产生贡献才需要保留,否则可以直接舍弃。

如果接下来乘上的若干 $G(x)$ 之积的次数界为 $k$,那么就只需要保留 $Q(x)$ 的最后 $k$ 项,因为前面的若干项不可能对最高项产生贡献。

即,递归到区间 $[l,r]$ 的 $Q(x)$ 只需要保留最高的 $r-l+1$ 项。

这样,每次只需要一次多项式乘法而不是多项式取模就可以使得问题规模减半,效率得到了显著的提升(对我而言,五倍)。

此外,还需要注意:最底层需要第 $n$ 项,但是递归计算时当成只有 $m$ 项来处理,这样会导致 $n>m$ 时答案会错。

这有两个方法解决,第一种就是将 $F(x)$ 对 $G(x)$ 取模。第二种更简单,令 $m=\max(m,n)$,多出来的一些询问不会对答案造成影响。

(一千个人就有一千个多项式板子,代码仅供参考)

#include<bits/stdc++.h>
using namespace std;
#define I inline int
#define V inline void
#define ll long long int
#define isnum(ch) ('0'<=ch&&ch<='9')
#define FOR(i,a,b) for(int i=a;i<=b;i++)
#define ROF(i,a,b) for(int i=a;i>=b;i--)
#define gc (_op==_ed&&(_ed=(_op=_buf)+fread(_buf,1,100000,stdin),_op==_ed)?EOF:*_op++)
char _buf[100000],*_op(_buf),*_ed(_buf);
const int N=1<<18|1,mod=998244353;
V check(int&x){x-=mod,x+=x>>31&mod;}
I getint(){
    int _s=0;char _ch=gc;
    while(!isnum(_ch))_ch=gc;
    while(isnum(_ch))_s=_s*10+_ch-48,_ch=gc;
    return _s;
}
I Pow(ll t,int x){
    ll s=1;
    for(;x;x>>=1,t=t*t%mod)if(x&1)s=s*t%mod;
    return s;
}
namespace poly{
    int lmt,w[N],r[N];
    V init(int n){
        int l=-1,wn;
        for(lmt=1;lmt<=n;)lmt<<=1,l++;
        FOR(i,0,lmt-1)r[i]=(r[i>>1]>>1)|((i&1)<<l);
        wn=Pow(3,mod>>++l),w[lmt>>1]=1;
        FOR(i,(lmt>>1)+1,lmt-1)w[i]=1ll*w[i-1]*wn%mod;
        ROF(i,(lmt>>1)-1,1)w[i]=w[i<<1];
    }
    V cl(int*a,int n){memset(a,0,n<<2);}
    I getLen(int n){return 1<<32-__builtin_clz(n);}
    V mul(int*a,int x,int n,int*b){while(n--)*b++=1ll**a++*x%mod;}
    V dot(int*a,int*b,int n,int*c){while(n--)*c++=1ll**a++**b++%mod;}
    V DFT(int*a,int l){
        static unsigned ll tmp[N];
        int u=__builtin_ctz(lmt/l),t;
        FOR(i,0,l-1)tmp[i]=a[r[i]>>u];
        for(int i=1;i^l;i<<=1)for(int j=0,d=i<<1;j^l;j+=d)FOR(k,0,i-1)
            t=tmp[i|j|k]*w[i|k]%mod,tmp[i|j|k]=tmp[j|k]+mod-t,tmp[j|k]+=t;
        FOR(i,0,l-1)a[i]=tmp[i]%mod;
    }
    V IDFT(int*a,int l){reverse(a+1,a+l),DFT(a,l),mul(a,mod-mod/l,l,a);}
    V Inv(const int*a,int n,int*b){
        static int A[N],B[N],tmp[N],d,l;
        tmp[0]=Pow(a[0],mod-2),cl(A,d),cl(B,d);
        for(d=1,l=2;d<n;d<<=1,l<<=1){
            copy(a,a+min(l,n),A),copy(tmp,tmp+d,B);
            DFT(A,l),DFT(B,l),dot(A,B,l,A),IDFT(A,l);
            cl(A,d),DFT(A,l),dot(A,B,l,A),IDFT(A,l);
            copy(A+d,A+l,tmp+d),mul(tmp+d,mod-1,d,tmp+d);
        }
        copy(tmp,tmp+n,b);
    }
    int*f[N],*g[N],bin[N<<5],*np(bin);
    V Mul(int*a,int*b,int n,int m,int*c){
        static int A[N],B[N],l;
        l=getLen(n+m-1),copy(a,a+n,A),copy(b,b+m,B);
        DFT(A,l),DFT(B,l),dot(A,B,l,A),IDFT(A,l);
        copy(A,A+n+m-1,c),cl(A,l),cl(B,l);
    }
    V eva_init(int p,int l,int r,int*a){
        g[p]=np,np+=r-l+2,f[p]=np,np+=r-l+2;
        if(l==r)return g[p][0]=1,check(g[p][1]=mod-a[l]);
        int lc=p<<1,rc=lc|1,mid=l+r>>1,len1=mid-l+2,len2=r-mid+1;
        eva_init(lc,l,mid,a),eva_init(rc,mid+1,r,a);
        Mul(g[lc],g[rc],len1,len2,g[p]);
    }
    V Mult(int*a,int*b,int n,int m,int*c){
        static int A[N],B[N],l;
        l=getLen(n),copy(a,a+n,A),reverse_copy(b,b+m,B);
        DFT(A,l),DFT(B,l),dot(A,B,l,A),IDFT(A,l);
        copy(A+m-1,A+n,c);
        cl(A,l),cl(B,l);
    }
    V eva_work(int p,int l,int r,int*a){
        if(l==r)return void(a[l]=f[p][0]);
        int lc=p<<1,rc=lc|1,mid=l+r>>1,len1=mid-l+2,len2=r-mid+1;
        Mult(f[p],g[rc],r-l+1,len2,f[lc]);
        eva_work(lc,l,mid,a);
        Mult(f[p],g[lc],r-l+1,len1,f[rc]);
        eva_work(rc,mid+1,r,a);
    }
    V eva(int*a,int*b,int n,int m,int*c){
        static int X[N],Y[N],l;
        eva_init(1,1,m,b),Inv(g[1],m+1,X);
        reverse(X,X+m+1),Mul(a,X,n,m+1,Y);
        copy(Y+n,Y+n+m,f[1]),eva_work(1,1,m,c);
        FOR(i,1,m)check(c[i]=1ll*c[i]*b[i]%mod+a[0]);
    }
}
int n,m,a[N],b[N],c[N];
int main(){
    n=getint()+1,m=getint(),poly::init(max(n,m+1)<<1);
    FOR(i,0,n-1)a[i]=getint();
    FOR(i,1,m)b[i]=getint();
    n=max(n,m+1),poly::eva(a,b,n,max(n-1,m),c);
    FOR(i,1,m)cout<<c[i]<<'\n';
    return 0;
}

P.S. 代码里还使用了一个优化,就是从上往下计算 $Q$ 的时候,会有长度为 $l$ 的多项式乘上长度为 $\frac{l}{2}$ 的多项式然后只取后 $\frac{l}{2}$ 项的操作,这部分可以利用 $FFT$ 循环卷积的性质只进行长度为 $l$ 的 $DFT$。

 
反馈
如果你认为某个题解有问题,欢迎向洛谷反馈,以帮助更多的同学。



请具体说明理由,以增加反馈的可信度。