zoukankan      html  css  js  c++  java
  • 多项式全家桶

    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;
    }
    
    呐,这份感情,什么时候可以传达呢
  • 相关阅读:
    oracle(Xe)数据库远程连接需修改配置参数
    oracl 权限循环查询
    控件网站
    java常用类(1)
    关于webdriver和谷歌浏览器的那些事
    2020年第27周,24.75h,完成计算智能/物联网/数据挖掘大作业
    2020年第26周,24.75h,计算智能的大小作业
    2020年第25周,25.5h,随机过程考试、report和计算智能作业
    2020年24周,11.75h,以完成作业和考试为主,看了一点点论文
    2020年第23周,11h,努力完成课程作业
  • 原文地址:https://www.cnblogs.com/Ame-sora/p/14311389.html
Copyright © 2011-2022 走看看