多项式全家桶

NTT

太多了不写了,直接贴个板子算了

inline void NTT(int *y , int len , int opt){
    int *rev = new int[len];
    rev[0] = 0;
    for(int i = 1;i < len;i ++) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) * (len >> 1));
    for(int i = 0;i < len;i ++) if(rev[i] > i) swap(y[rev[i]] , y[i]);
    for(int i = 1;i < len;i <<= 1){
        int G1 = quick_pow(3 , (mod - 1) / (i << 1));
        for(int j = 0;j < len;j += (i << 1)){
            for(int k = 0 , g = 1;k < i;k ++ , g = 1LL * g * G1 % mod){
                int res = 1LL * y[i + j + k] * g % mod;
                y[i + j + k] = ((y[j + k] - res) % mod + mod) % mod;
                y[j + k] = (y[j + k] + res) % mod;
            }
        }
    }
    if(opt == -1){
        reverse(y + 1 , y + len);
        for(int i = 0 , inv = quick_pow(len , mod - 2);i < len;i ++) y[i] = 1LL * y[i] * inv % mod;
    }
    delete []rev; rev = 0x0;
}

多项式求逆

(A(x) imes B(x) equiv1pmod {x^{n}})

则有(A(x) imes B(x) equiv1pmod {x^{ leftlceildfrac{n}{2} ight ceil}})

又设(A(x) imes C(x) equiv1pmod {x^{ leftlceildfrac{n}{2} ight ceil}})

那么有(B(x)-C(x)equiv0pmod { x^{leftlceildfrac{n}{2} ight ceil}})

平方可得(B(x)^2-2 imes B(x) imes C(x)-C(x)^2equiv 0pmod{x^n})

那么有(A(x) imes B(x)^2-2 imes A(x) imes B(x) imes C(x) +A(x) imes C(x)^2equiv 0pmod{x^n})

(B(x)-2 imes C(x)+A(x) imes C(x)^2equiv 0pmod{x^n})

(B(x)=2C(x)-A(x) imes C(x)^2pmod{x^n})

inline void poly_inv(int *a , int len){
    if(len == 1){ return void(a[0] = quick_pow(a[0] , mod - 2)); }
    int len1 = (len + 1) / 2;
    int *g = new int[len * 2];
    for(int i = 0;i < len1;i ++) g[i] = a[i];
    for(int i = len1;i < 2 * len;i ++) g[i] = 0;
    poly_inv(g , len1);
    for(int i = len1;i < 2 * len;i ++) g[i] = 0;
    NTT(g , 2 * len , 1); NTT(a , 2 * len , 1);
    for(int i = 0;i < 2 * len;i ++) a[i] = ((((1LL * 2 * g[i]) % mod) - 1LL * a[i] * g[i] % mod * g[i] % mod) % mod + mod) % mod;
    NTT(a , 2 * len , -1);
    for(int i = len;i < 2 * len;i ++) a[i] = 0;
    delete []g; g = 0x0;
}

多项式开根

(C(x)^2equiv A(x)pmod{x^n})

(C(x)^2-A(x)equiv 0pmod{x^n})

(C(x)^4-2 imes C(x)^2 imes A(x)+A(x)^2equiv 0pmod{x^{2n}})

两边同时加上(4 imes C(x)^2 imes A(x))可得

(C(x)^4+2 imes C(x)^2 imes A(x)+A(x)^2equiv 4 imes C(x)^2 imes A(x)pmod{x^{2n}})

(frac{C(x)^4+2 imes C(x)^2 imes A(x)+A(x)^2}{4C(x)^2}equiv A(x)pmod{x^{2n}})

那么最后答案为(B(x)=frac{C(x)^2+ A(x)}{2 imes C(x)})

inline void poly_sqrt(int *a , int len){
    if(len == 1) return;
    int len1 = len / 2;
    int *f1 = new int[len * 2];
    for(int i = 0;i < len1;i ++) f1[i] = a[i];
    for(int i = len1;i < 2 * len;i ++) f1[i] = 0;
    poly_sqrt(f1 , len1);
    for(int i = len1;i < 2 * len;i ++) f1[i] = 0;
    int *f2 = new int[len * 2];
    for(int i = 0;i < 2 * len;i ++) f2[i] = 1LL * f1[i] * 2 % mod;
    poly_inv(f2 , len);
    for(int i = len;i < 2 * len;i ++) f2[i] = 0;
    NTT(f1 , 2 * len , 1);
    for(int i = 0;i < 2 * len;i ++) f1[i] = 1LL * f1[i] * f1[i] % mod;
    NTT(f1 , 2 * len , -1);
    for(int i = 0;i < len;i ++) a[i] = (a[i] + f1[i]) % mod;
    NTT(a , 2 * len , 1); NTT(f2 , 2 * len , 1);
    for(int i = 0;i < 2 * len;i ++) a[i] = 1LL * a[i] * f2[i] % mod;
    NTT(a , 2 * len , -1);
    for(int i = len;i < 2 * len;i ++) a[i] = 0;
    delete []f1; f1 = 0x0;
    delete []f2; f2 = 0x0;
}

多项式求(mathcal{ln})

考虑一个多项式的导数(g(x)^{'}=frac{f(x)^{'}}{f(x)})

(g(x))积分回去可得(int{g(x)}=mathcal{ln}f(x))


inline void poly_dev(int *a , int len){
    for(int i = 1;i < len;i ++) a[i - 1] = 1LL * a[i] * i % mod;
    a[len - 1] = 0;
}

inline void poly_inv_dev(int *a , int len){
    for(int i = len + 1 ; i ; i --) a[i] = 1LL * a[i - 1] * quick_pow(i , mod - 2) % mod;
    a[0] = 0;
}

inline void poly_ln(int *a , int len){
    int *f = new int[2 * len];
    for(int i = 0;i < len;i ++) f[i] = a[i];
    for(int i = len;i < len * 2;i ++) f[i] = 0;
    poly_dev(f , len); poly_inv(a , len);
    NTT(f , 2 * len , 1); NTT(a , 2 * len , 1);
    for(int i = 0;i < 2 * len;i ++) a[i] = 1LL * f[i] * a[i] % mod;
    NTT(a , 2 * len , -1);
    poly_inv_dev(a , len);
    for(int i = len;i < 2 * len;i ++) a[i] = 0;
    delete []f; f = 0x0;
}

多项式求(mathcal{exp})

(g(x)equiv e^{f(x)}pmod{x^n})

同时取(mathcal{ln})可得(lng(x)-f(x)=0)

则有(h(g(x))=ln(g(x))-f(x))

(h(g(x))^{'}=frac{1}{g(x)})

(h(g(x))=frac{h(g_0(x))}{0!}+frac{h(g_0(x)^{'})}{1!}+......)

(h(g_0(x))+h(g_0(x))^{'} imes (g(x)-g(x_0))equiv 0pmod{x^{2n}})

(g(x)equiv g_0(x) imes(f(x)-lng_0(x))+g_0(x)pmod{x^{2n}})

inline void poly_exp(int *a , int len){
    if(len == 1) return void(a[0] ++);
    int *f = new int[len * 2]; int len1 = len / 2;
    for(int i = 0;i < len1;i ++) f[i] = a[i];
    for(int i = len1;i < len;i ++) f[i] = 0;
    poly_exp(f , len1);
    for(int i = len1;i < len * 2;i ++) f[i] = 0;
    int *g = new int[len * 2];
    for(int i = 0;i < len * 2;i ++) g[i] = f[i];
    poly_ln(g , len); a[0] ++;
    for(int i = 0;i < len;i ++) a[i] = ((a[i] - g[i]) % mod + mod) % mod;
    NTT(a , len * 2 , 1); NTT(f , len * 2 , 1);
    for(int i = 0;i < len * 2;i ++) a[i] = 1LL * a[i] * f[i] % mod;
    NTT(a , len * 2 , -1);
    for(int i = len;i < len * 2;i ++) a[i] = 0;
    delete []g; g = 0x0; 
    delete []f; f = 0x0;
}

多项式求三角函数

根据欧拉公式(e^{ix}=cos x+isin x)

我们同理有(e^{-ix}=cos x-isin x)

于是我们相加或者相减即可得到(egin{aligned} sin x&=frac{e^{ix}-e^{-ix}}{2i}end{aligned}) (egin{aligned}cos x&=frac{e^{ix}+e^{-ix}}{2} end{aligned})

直接上(exp)即可,(i)为在(w_4),在膜(998244353)下为(g^{frac{p-1}{4}})

inline void poly_sin(int *a , int len){
    static int f[kato] , g[kato];
    for(int i = 0;i < len;i ++) f[i] = 1LL * a[i] * img % mod;
    poly_exp(f , len);
    for(int i = 0;i < len;i ++) g[i] = f[i];
    poly_inv(g , len);
    for(int i = 0 , inv = quick_pow(2 * img , mod - 2);i < len;i ++) a[i] = 1LL * inv * (((f[i] - g[i]) % mod + mod) % mod) % mod;
}

inline void poly_cos(int *a , int len){
    static int f[kato] , g[kato];
    for(int i = 0;i < len;i ++) f[i] = 1LL * a[i] * img % mod;
    poly_exp(f , len);
    for(int i = 0;i < len;i ++) g[i] = f[i];
    poly_inv(g , len);
    for(int i = 0;i < len;i ++) a[i] = 1LL * inv2 * ((f[i] + g[i]) % mod) % mod;
}

多项式求反三角函数

对于(arcsin)我们有(G(x)equivarcsin F(x)pmod{x^n})

两边求导可得(G'(x)equivfrac{F'(x)}{sqrt{1-F(x)^2}}pmod{x^n})

然后再积回去(G(x)equivintfrac{F'(x)}{sqrt{1-F(x)^2}}mathrm{d}xpmod{x^n})

对于(arccos)同理计算即可最后为(G(x)equiv-intfrac{F'(x)}{sqrt{1-F(x)^2}}mathrm{d}xpmod{x^n})

对于(arctan)也一样最后为(G(x)equivintfrac{F'(x)}{1+F(x)^2}mathrm{d}xpmod{x^n})

求导+开根+求逆+积分一套理论

inline void poly_asin(int *a , int len){
    int *f = new int[kato];
    for(int i = 0;i < len;i ++) f[i] = a[i];
    for(int i = len;i < 2 * len;i ++) f[i] = 0;
    NTT(f , 2 * len , 1);
    for(int i = 0;i < 2 * len;i ++) f[i] = 1LL * f[i] * f[i] % mod;
    NTT(f , 2 * len , -1);
    for(int i = len;i < 2 * len;i ++) f[i] = 0;
    for(int i = 0;i < len;i ++) f[i] = (mod - f[i]) % mod;
    f[0] = 1;
    poly_sqrt(f , len); poly_inv(f , len); poly_dev(a , len);
    NTT(a , 2 * len , 1); NTT(f , 2 * len , 1);
    for(int i = 0;i < 2 * len;i ++) a[i] = 1LL * a[i] * f[i] % mod;
    NTT(a , 2 * len , -1);
    for(int i = len;i < 2 * len;i ++) a[i] = 0;
    poly_dev_inv(a , len);
    delete []f; f = 0x0;
}

inline void poly_atan(int *a , int len){    
    int *f = new int[kato];
    for(int i = 0;i < len;i ++) f[i] = a[i];
    for(int i = len;i < 2 * len;i ++) f[i] = 0;
    NTT(f , 2 * len , 1);
    for(int i = 0;i < 2 * len;i ++) f[i] = 1LL * f[i] * f[i] % mod;
    NTT(f , 2 * len , -1);
    for(int i = len;i < 2 * len;i ++) f[i] = 0;
    f[0] = 1;
    poly_inv(f , len); poly_dev(a , len);
    NTT(a , 2 * len , 1); NTT(f , 2 * len , 1);
    for(int i = 0;i < 2 * len;i ++) a[i] = 1LL * a[i] * f[i] % mod;
    NTT(a , 2 * len , -1);
    for(int i = len;i < 2 * len;i ++) a[i] = 0;
    poly_dev_inv(a , len);
    delete []f; f = 0x0;
}

下降幂多项式乘法

考虑构造关于(A)(B)(EGF)

先来点(EGF)的性质

[sum_{i=0}^{infty}frac{i^{underline{n}}}{i!}x^i=sum_{i=0}^{infty}frac{1}{(i-n)!}x^i=sum_{i=n}^{infty}frac{1}{i!}x^{i-n}=x^nsum_{i=0}^{infty}frac{1}{i!}x^i=x^ne^x ]

考虑(A)的下降幂为(F(n)=sum_{i=0}^{infty}F[i]n^{underline{i}})

那么考虑构造其(EGF)

[sum_{i=0}^{infty}frac{x^i}{i!}sum_{j=0}^{infty}F[j]i^{underline{j}} ]

[=sum_{j=0}^{infty}F[j]sum_{i=0}^{infty}frac{i^{underline{j}}}{i!}x^i ]

[=sum_{j=0}^{infty}F[j]x^je^x ]

卷积然后卷上(e^{-x})即可,最后记得求出来的是(EGF),乘上阶乘

for(int i = 0;i <= yuni;i ++){
    A[i] = fac_inv[i];
    B[i] = (i & 1) ? mod - fac_inv[i] : fac_inv[i];
}
for(len = 1;len <= 2 * yuni;len <<= 1);
NTT(a , len , 1); NTT(A , len , 1);
for(int i = 0;i < len;i ++) a[i] = 1LL * a[i] * A[i] % mod;
NTT(a , len , -1);
NTT(b , len , 1); NTT(B , len , 1);
for(int i = 0;i < len;i ++) b[i] = 1LL * b[i] * A[i] % mod;
NTT(b , len , -1);
for(int i = 0;i <= yuni;i ++) ans[i] = 1LL * a[i] * b[i] % mod * fac[i] % mod;
NTT(ans , len , 1);
for(int i = 0;i < len;i ++) ans[i] = 1LL * ans[i] * B[i] % mod;
NTT(ans , len , -1);
for(int i = 0;i <= n + m;i ++) printf("%d " , ans[i]);

任意模数多项式乘法

假设这一位的答案是 (x),三个模数分别为 (A,B,C),那么:

[egin{aligned}xequiv x_1pmod{A} \ xequiv x_2pmod{B} \ xequiv x_3pmod{C}end{aligned} ]

先对前两个用(CRT)合并我们有

[egin{aligned}x_1+k_1A=x_2+k_2B\x_1+k_1Aequiv x_2pmod{B}\k_1equiv frac{x_2-x_1}Apmod{B}end{aligned} ]

于是求出了 (k_1),也就求出了 (xequiv x_1+k_1Apmod{AB}),记 (x_4=x_1+k_1A)

[egin{aligned}x_4+k_4AB=x_3+k_3C\x_4+k_4ABequiv x_3pmod{C}\k_4equiv dfrac{x_3-x_4}{AB}pmod{C}end{aligned} ]

求出了 (k_4)(xequiv x_4+k_4ABpmod{ABC}),因为 (xlt ABC),所以 (x=x_4+k_4AB)

const int mod1 = 998244353 , mod2 = 1004535809 , mod3 = 469762049;

const LL mod_1_2 = 1LL * mod1 * mod2;

const int inv1 = quick_pow(mod1 , mod2 - 2 , mod2) , inv2 = quick_pow(mod_1_2 % mod3 , mod3 - 2 , mod3);

struct Int{
    int A , B , C;

    Int(int A = 0 , int B = 0 , int C = 0): A(A) , B(B) , C(C){}

    inline void init(int _n){ A = _n; B = _n; C = _n; }

    inline friend Int operator+(const Int &a , const Int &b){
        return Int(add(a.A , b.A , mod1) , add(a.B , b.B , mod2) , add(a.C , b.C , mod3));
    }

    inline friend Int operator-(const Int &a , const Int &b){
        return Int(del(a.A , b.A , mod1) , del(a.B , b.B , mod2) , del(a.C , b.C , mod3));
    }

    inline friend Int operator*(const Int &a , const Int &b){
        return Int(1LL * a.A * b.A % mod1 , 1LL * a.B * b.B % mod2 , 1LL * a.C * b.C % mod3);
    }
    
    inline int get_ans(){
        LL x = static_cast<LL>(B - A + mod2) % mod2 * inv1 % mod2 * mod1 + A;
        return (static_cast<LL>(C - x % mod3 + mod3) % mod3 * inv2 % mod3 * (mod_1_2 % mod) % mod + x) % mod;  
    }
};

二次剩余多项式开根

struct Complex{
    int x , y;
    Complex(int x = 0 , int y = 0): x(x) , y(y) { }
    friend bool operator ==(const Complex &a , const Complex &b){
        return a.x == b.x && a.y == b.y;
    }
    friend Complex operator *(const Complex &a , const Complex &b){
        return Complex((static_cast<LL>(a.x) * b.x % mod + static_cast<LL>(a.y) * b.y % mod * qaq % mod) % mod , (static_cast<LL>(a.x) * b.y % mod + static_cast<LL>(a.y) * b.x % mod) % mod);
    }
};

inline Complex quick_pow(Complex a , int b){
    Complex res = Complex(1 , 0);
    for(; b ; b >>= 1 , a = a * a){
        if(b & 1){
            res = res * a;
        }
    }
    return res;
}

inline bool judge(int x){
    return quick_pow(x , (mod - 1) >> 1) == 1;
}

inline bool Cipolla(int n , int p , int &ans0 , int &ans1){
    if(quick_pow(n , (mod - 1) / 2) == mod - 1) return 0;
    int x = rand() % mod;
    while(!x || judge((static_cast<LL>(x) * x % mod + mod - n) % mod)) x = rand() % mod;
    qaq = (static_cast<LL>(x) * x % mod + mod - n) % mod;
    ans0 = quick_pow(Complex(x , 1) , (mod + 1) >> 1).x;
    ans1 = mod - ans0;
    if(ans0 > ans1) swap(ans0 , ans1);
    return 1;
}

inline void poly_sqrt(int *a , int len){
    if(len == 1){ Cipolla(a[0] , mod , ans0 , ans1) , a[0] = min(ans0 , ans1); return; }
    int len1 = len / 2;
    int *f1 = new int[len * 2];
    for(int i = 0;i < len1;i ++) f1[i] = a[i];
    for(int i = len1;i < 2 * len;i ++) f1[i] = 0;
    poly_sqrt(f1 , len1);
    for(int i = len1;i < 2 * len;i ++) f1[i] = 0;
    int *f2 = new int[len * 2];
    for(int i = 0;i < 2 * len;i ++) f2[i] = 1LL * f1[i] * 2 % mod;
    poly_inv(f2 , len);
    for(int i = len;i < 2 * len;i ++) f2[i] = 0;
    NTT(f1 , 2 * len , 1);
    for(int i = 0;i < 2 * len;i ++) f1[i] = 1LL * f1[i] * f1[i] % mod;
    NTT(f1 , 2 * len , -1);
    for(int i = 0;i < len;i ++) a[i] = (a[i] + f1[i]) % mod;
    NTT(a , 2 * len , 1); NTT(f2 , 2 * len , 1);
    for(int i = 0;i < 2 * len;i ++) a[i] = 1LL * a[i] * f2[i] % mod;
    NTT(a , 2 * len , -1);
    for(int i = len;i < 2 * len;i ++) a[i] = 0;
    delete []f1; f1 = 0x0;
    delete []f2; f2 = 0x0;
}
呐,这份感情,什么时候可以传达呢
原文地址:https://www.cnblogs.com/Ame-sora/p/14311389.html