多项式exp

写了一发多项式exp，picks太神辣

#include <bits/stdc++.h>

using namespace std;

typedef long long ll;
const int N = 262144,mod = 998244353,G = 3;

void read(int&x){
    x = 0;char ch = getchar();
    while (ch < '0' || '9' < ch) ch = getchar();
    while ('0' <= ch&& ch <= '9') x = x * 10 + ch - '0',ch = getchar(); 
}
int power(int x,int y,int z){
    int ans = 1;
    while (y){
        if (y&1) ans = (ll)ans*x%z;
        x = (ll)x*x%z;
        y >>= 1;
    }
    return ans;
}
int add(int x,int y){
    x += y;
    while (x >= mod) x -= mod;
    return x;
}


int a[N+10],b[N+10],c[N+10],d[N+10],p[N];
int pg[N+10],pg_inv[N+10],inv[N+10];
int n,m;
void FFT(int*,int,int);
void poly_exp(int*,int*,int*,int*,int);
void poly_inv(int*,int*,int*,int);
void poly_ln(int*,int*,int*,int);
void poly_der(int*,int);
void poly_mot(int*,int); 
int main(){
    read(n);
    for (int i = 0;i < n;i++) read(a[n-i-1]);
    for (m = 1;m < n;m<<=1);
    for(int i = 1;i <= (m<<1);i <<= 1){
        pg[i] = power(G,(mod-1)/i,mod);
        pg_inv[i] = power(pg[i],mod-2,mod);
    }
    inv[1] = 1;
    for (int i = 2;i <= N;i++)inv[i] = mod-(ll)mod/i*inv[mod%i]%mod;
    /* 
    poly_inv(a,b,c,m); 
//    for (int i = 0;i < m<<1;i++)cout << a[i]<<' ';cout<<endl;
//    for (int i = 0;i < m<<1;i++)cout << b[i]<<' ';cout<<endl;
//    cout <<(ll)a[0]*b[0]%mod<<endl;
    FFT(a,1,m<<1);
//    for (int i = 0;i < m<<1;i++)cout <<a[i]<<'*';cout<<endl;
    FFT(b,1,m<<1);
    for(int i = 0;i <m<<1;i++) a[i] =(ll)a[i]*b[i]%mod;
    FFT(a,-1,m<<1);
//    for (int i = 0;i < m;i++)cout << a[i]<<' ';cout<<endl;
    bool flag = 1;
    for(int i = 1;i < m;i++)if (a[i] != 0) flag = 0;
    if(a[0] != 1) flag = 0;
    puts(flag ? "YES" : "NO");
    */ 
//    /*
    poly_exp(a,b,c,d,m);
//    for (int i = 0;i < m<<1;i++) cout << a[i]<<' ';cout << endl;
//    for (int i = 0;i < m<<1;i++)cout << b[i] <<' ';cout<<endl;
    /*
//    b[0] = 1;b[1] = 1;b[2] = inv[2];b[3] = inv[6];
    poly_ln(b,c,d,m);bool flag = 1;
//    for (int i = 0;i < m;i++) cout << c[i]<<' ';cout<<endl;
    for (int i = 0;i <m;i++)if(c[i] != a[i]) flag = 0;
    puts(flag ? "YES" : "NO");
    */
//    */
    return 0;
}
void FFT(int *a,int d,int deg){
    for (int i = 1;i < deg;i <<= 1)
        for (int j = 0;j < i;j++)
            p[i+j] = p[j]+deg/i/2;
    for (int i = 0;i <deg;i++) if (p[i] > i) swap(a[i],a[p[i]]);
    for (int n = 2;n <= deg;n <<= 1){
        int wn = (d == 1 ? pg[n] : pg_inv[n]);
        for (int i = 0;i <deg;i += n){
            int w = 1;
            for (int j = i;j < i + n/2;j++){
                int x = a[j],y = (ll)a[j+n/2]*w%mod;
                a[j] = add(x,y);
                a[j+n/2] = add(x,mod-y);
                w = (ll)w*wn%mod;
            }
        }
    }
    if (d == -1){
        for (int i = 0;i < deg;i++)
            a[i] = (ll)a[i]*inv[deg]%mod;
    }
}
void poly_der(int*a,int deg){
    for (int i = 0;i < deg;i++)
        a[i] = (ll)a[i+1]*(i+1)%mod;
}
void poly_mot(int*a,int deg){
    for (int i = deg;i > 0;i--)
        a[i] = (ll)a[i-1]*inv[i]%mod;
    a[0] = 0;
}
void poly_inv(int*a,int*b,int*c,int deg){
    if (deg == 1){
        if (a[0] <= N) b[0] = inv[a[0]];else 
        b[0] = power(a[0],mod-2,mod);
        return;
    }
    poly_inv(a,b,c,deg/2);
    int deg2 = deg<<1;
    //B = B(2-BA)
    for (int i = 0;i < deg;i++) c[i] = a[i];
    for (int i = deg;i < deg2;i++) c[i] = 0;
    FFT(c,1,deg2);
    FFT(b,1,deg2);
    for (int i = 0;i < deg2;i++) b[i] = (ll)b[i]*add(2,mod-(ll)b[i]*c[i]%mod)%mod;
    FFT(b,-1,deg2);
    for (int i = deg;i < deg2;i++) b[i] = 0;
}
void poly_ln(int*a,int*b,int*c,int deg){
    /*
        b = lna
        b' = a'/a
    */
    int deg2 = deg<<1;
    for (int i = 0;i < deg2;i++)c[i] = 0;
    poly_inv(a,c,b,deg);
    for(int i = 0;i <= deg;i++) b[i] = a[i];
    for(int i = deg+1;i < deg2;i++) b[i] = 0;
    poly_der(b,deg);
    FFT(b,1,deg2);
    FFT(c,1,deg2);
    for (int i = 0;i < deg2;i++) b[i] = (ll)b[i]*c[i]%mod;
    FFT(b,-1,deg2);
    for (int i = deg;i < deg2;i++) b[i] = 0;
    poly_mot(b,deg);
    for(int i = deg;i < deg2;i++)b[i] = 0;
}
void poly_exp(int *a,int *b,int *c,int*d,int deg){
    //B = B(1-lnB+A)
    if(deg == 1){
        //if a[0] = 0
        b[0] = 1;
        return;
    }
    int deg2=deg<<1;    
    poly_exp(a,b,c,d,deg/2);
    for (int i = deg;i < deg2;i++) b[i] = 0;
    poly_ln(b,c,d,deg);

    for (int i = 0;i < deg;i++) c[i] = add(a[i]+(i==0),mod-c[i]);
//    cout <<"1-lnB+A ";for (int i = 0;i <deg2;i++) cout << c[i]<<' ';cout<<endl;    
    for (int i = deg;i < deg2;i++) c[i] = 0;
    FFT(c,1,deg2);
    FFT(b,1,deg2);
    for(int i = 0;i < deg2;i++) b[i] = (ll)b[i]*c[i]%mod;
    FFT(b,-1,deg2);
    for (int i = deg;i < deg2;i++) b[i] = 0;
//    for (int i = 0;i < deg2;i++)cout << b[i]<<' ';cout<<"deg:"<<deg<<"\n";
}

因为只会求e^0的值，所以只能解决常数项为0的多项式的exp。。。。

多项式exp其实求的是e^A的泰勒展开的前若干项，

用牛顿迭代法，求出解的递推式

最后一行在多项式的意义下在这个情况中证明了这个方法至少是平方收敛（多项式下的牛顿迭代应该都是可以类似证精度的）（数意义下的不会证orz）

然后。。我试图用原始的初等的推多项式求逆的那一套方法推这个式子，成功GG。。。。

然后写的时候发现一些问题，就是这个求原函数啊，求模x^n意义下的原函数需要模x^(n+1)意义下的函数，这题因为需要求原函数的函数都是没取模的所以不影响，不然还有注意一下细节。。。