LOJ#556. 「Antileaf's Round」咱们去烧菜吧（牛顿迭代）

快速傅里叶变换与多项式乘法

最新推荐文章于 2023-11-21 12:28:31 发布

原创最新推荐文章于 2023-11-21 12:28:31 发布 · 624 阅读

0 ·

CC 4.0 BY-SA版权

牛顿迭代专栏收录该内容

6 篇文章

订阅专栏

本文深入探讨了快速傅里叶变换（FFT）在多项式乘法中的应用，通过高效的算法实现，解决了大规模多项式运算的时间复杂度问题。文章详细介绍了FFT的原理，包括多项式的表示、卷积定理以及逆变换过程，并提供了具体的代码实现，展示了如何利用FFT进行高效的大规模多项式乘法。

传送门

题解：
考前练练手感。

G = \prod F i

$G = \prod F_i$

F i = {1 - x a i ( b i + 1 ) 1 - x a i 1 1 - x a i b i > 0 b i = 0

$F_i = \begin{cases} \frac{1-x^{a_i(b_i+1)}}{1-x^{a_i}} & b_i>0\\ \frac{1}{1-x^{a_i}} & b_i = 0\end{cases}$

$\ln G = \sum \ln\frac{1-x^{a_i(b_i+1)}}{1-x^{a_i}}$
求 $\ln$ 然后相加泰勒展开后是一个调和级数，之后把 $G$ 再 $\exp$ 出来，用牛顿迭代法，时间复杂度 $O(n \log n)$ .

#include <bits/stdc++.h>
using namespace std;
typedef long long LL;

const int RLEN=1<<18|1;
inline char nc() {
    static char ibuf[RLEN],*ib,*ob;
    (ib==ob) && (ob=(ib=ibuf)+fread(ibuf,1,RLEN,stdin));
    return (ib==ob) ? -1 : *ib++;
}
inline int rd() {
    char ch=nc(); int i=0,f=1;
    while(!isdigit(ch)) {if(ch=='-')f=-1; ch=nc();}
    while(isdigit(ch)) {i=(i<<1)+(i<<3)+ch-'0'; ch=nc();}
    return i*f;
}
inline void W(int x) {
    static int buf[50];
    if(!x) {putchar('0'); return;}
    if(x<0) {putchar('-'); x=-x;}
    while(x) {buf[++buf[0]]=x%10; x/=10;}
    while(buf[0]) {putchar(buf[buf[0]--]+'0');}
}

const int N=1e5+50, mod=998244353, G=3;
inline int add(int x,int y) {return (x+y>=mod) ? (x+y-mod) : (x+y);}
inline int dec(int x,int y) {return (x-y<0) ? (x-y+mod) : (x-y);}
inline int mul(int x,int y) {return (LL)x*y%mod;}
inline int power(int a,int b,int rs=1) {for(;b;b>>=1,a=mul(a,a)) if(b&1) rs=mul(rs,a); return rs;}
inline int cinv(int a) {return power(a,mod-2);}

int n,m,inv[N*8],cnt1[N],cnt2[N];
int k,w[N*8],pos[N*8],tp[N*8],tp2[N*8];
inline void init(int nn) {
    for(k=1;k<=nn;k<<=1);
    for(int i=1;i<k;i++) pos[i]=(i&1) ? ((pos[i>>1]>>1)^(k>>1)) : (pos[i>>1]>>1);
}
inline void dft(int *a,int o=1) {
    for(int i=1;i<k;i++) if(pos[i]>i) swap(a[pos[i]],a[i]);
    for(int bl=1;bl<k;bl<<=1) {
        int tl=bl<<1, wn=power(G,(mod-1)/tl);
        w[0]=1; for(int i=1;i<bl;i++) w[i]=mul(w[i-1],wn);
        for(int bg=0;bg<k;bg+=tl)
            for(int j=0;j<bl;j++) {
                int &t1=a[bg+j],&t2=a[bg+j+bl],t=mul(t2,w[j]);
                t2=dec(t1,t); t1=add(t1,t);
            }
    }
    if(!~o) {
        const int inv=cinv(k); reverse(a+1,a+k);
        for(int i=0;i<k;i++) a[i]=mul(a[i],inv); 
    }
}
struct poly {
    vector <int> a; int deg;
    poly(int d=0) {deg=d; a.resize(d+1);} 
    inline poly extend(int d) {poly c=*this; c.deg=d-1; c.a.resize(d); return c;}
    inline poly ig() {
        poly c(this->deg+1);        
        for(int i=c.deg;i>=1;i--) c.a[i]=mul(this->a[i-1],inv[i]);
        return c;
    }
    inline poly dg() {
        poly c(this->deg-1);
        for(int i=0;i<=c.deg;i++) c.a[i]=mul(this->a[i+1],i+1);
        return c;
    }
    inline poly inc(int t) {
        poly c=*this;
        for(int i=0;i<=deg;i++) c.a[i]=mul(a[i],t);
        return c;
    }
    friend inline poly operator *(const poly &a,const poly &b) {
        poly c(a.deg+b.deg); init(c.deg);
        memset(tp,0,sizeof(int)*k); memset(tp2,0,sizeof(int)*k);
        memcpy(tp,a.a.data(),sizeof(int)*(a.deg+1));
        memcpy(tp2,b.a.data(),sizeof(int)*(b.deg+1));
        dft(tp); dft(tp2); 
        for(int i=0;i<k;i++) tp[i]=mul(tp[i],tp2[i]);
        dft(tp,-1); memcpy(c.a.data(),tp,sizeof(int)*(c.deg+1));
        return c;
    }
    friend inline poly operator +(const poly &a,const poly &b) {
        poly c(max(a.deg,b.deg));
        for(int i=0;i<=a.deg;i++) c.a[i]=add(c.a[i],a.a[i]);
        for(int i=0;i<=b.deg;i++) c.a[i]=add(c.a[i],b.a[i]);
        return c;
    }
    friend inline poly operator -(const poly &a,const poly &b) {
        poly c(max(a.deg,b.deg));
        for(int i=0;i<=a.deg;i++) c.a[i]=add(c.a[i],a.a[i]);
        for(int i=0;i<=b.deg;i++) c.a[i]=dec(c.a[i],b.a[i]);
        return c;
    }
};
inline poly calc_inv(poly f,int len) {
    if(len==1) {poly t(0); t.a[0]=cinv(f.a[0]); return t;}
    poly f0=calc_inv(f.extend(len>>1),len>>1);
    poly t=f0.inc(2)-(f*f0).extend(len)*f0;
    return t.extend(len);
}
inline poly calc_ln(poly f,int len) {
    poly t=f.dg(), t2=calc_inv(f,len);
    return ((t*t2).ig()).extend(len);
}
inline poly calc_exp(poly f,int len) {
    if(len==1) {poly t(0); t.a[0]=1; return t;}
    poly f0=calc_exp(f.extend(len>>1),len>>1);
    poly t1=f-calc_ln(f0,len); t1.a[0]++;
    return (f0*t1).extend(len); 
}
int main() {
    n=rd(), m=rd(); inv[0]=inv[1]=1;
    for(int i=2;i<=n*8;i++) inv[i]=mod-mul(mod/i,inv[mod%i]);
    for(int i=1;i<=m;i++) {
        int x=rd(),y=rd();
        if(x<=n) ++cnt2[x];
        if(y && (1ll*x*(y+1))<=n) ++cnt1[x*(y+1)];  
    } poly g(n);
    for(int i=1;i<=n;i++) if(cnt2[i]) for(int j=i;j<=n;j+=i) g.a[j]=add(g.a[j],mul(cnt2[i],inv[j/i]));
    for(int i=1;i<=n;i++) if(cnt1[i]) for(int j=i;j<=n;j+=i) g.a[j]=dec(g.a[j],mul(cnt1[i],inv[j/i]));
    init(n); g=calc_exp(g,k);
    for(int i=1;i<=n;i++) W(g.a[i]), putchar('\n');
}