zoukankan      html  css  js  c++  java
  • P5205 【模板】多项式开根

    思路

    按如下式子计算即可

    [B(x)=frac{A(x)+B'^2(x)}{2B'(x)} ]

    代码

    // luogu-judger-enable-o2
    #include <cstdio>
    #include <cstring>
    #include <algorithm>
    #define int long long
    using namespace std;
    const int MAXN = 300000;
    const int G = 3;
    const int invG = 332748118;
    const int MOD = 998244353;
    int n;
    struct Poly{
        int t;//次数界
        int data[MAXN];
        Poly(){}
        Poly(int x,int val[]){
            for(int i=0;i<=x;i++)
                data[i]=val[i];
        }
    };
    int pow(int a,int b){
        int ans=1;
        while(b){
            if(b&1)
                ans=(1LL*ans*a)%MOD;
            a=(1LL*a*a)%MOD;
            b>>=1;
        }
        return ans;
    }
    void rever(Poly &a){
        for(int i=0,j=a.t;i<j;i++,j--){
            swap(a.data[i],a.data[j]);
        }
    }
    void save(Poly &a,int top){
        for(int i=top+1;i<=a.t;i++)
            a.data[i]=0;
        a.t=top;
    }
    void output(Poly a){
        putchar('
    ');
        printf("a.times=%lld
    ",a.t);
        putchar('
    ');
        for(int i=0;i<=a.t;i++)
            printf("%lld ",a.data[i]);
        putchar('
    ');
        putchar('
    ');
    }
    void NTT(Poly &a,int opt,int n){//1 DFT 0 IDFT
        int lim=0;
        while((1<<(lim))<n)
            lim++;
        n=(1<<lim);
        for(int i=0;i<n;i++){
            int t=0;
            for(int j=0;j<lim;j++)
                if((i>>j)&1)
                    t|=(1<<(lim-j-1));
            if(i<t)
                swap(a.data[i],a.data[t]);        
        }
        for(int i=2;i<=n;i<<=1){
            int len=i/2;
            int tmp=pow((opt)?G:invG,(MOD-1)/i);
            for(int j=0;j<n;j+=i){
                int arr=1;
                for(int k=j;k<j+len;k++){
                    int t=(1LL*a.data[k+len]*arr)%MOD;
                    a.data[k+len]=(a.data[k]-t+MOD)%MOD;
                    a.data[k]=(a.data[k]+t)%MOD;
                    arr=(1LL*arr*tmp)%MOD;
                }
            }
        }
        if(!opt){
            int invN = pow(n,MOD-2);
            for(int i=0;i<n;i++){
                a.data[i]=(a.data[i]*invN)%MOD;
            }
        }
    }
    void mul(Poly &a,Poly b){//a=a*b
        int num=(a.t+b.t),lim=0;
        while((1<<(lim))<=((num+2)))
            lim++;
        lim=(1<<lim);
        NTT(a,1,lim);
        NTT(b,1,lim);
        for(int i=0;i<lim;i++)
            a.data[i]=(1LL*a.data[i]*b.data[i])%MOD;
        NTT(a,0,lim);
        a.t=num;
        for(int i=num+1;i<lim;i++)
            a.data[i]=0;
    }
    void Inv(Poly a,Poly &inv,int dep,int &len){//
        if(dep==1){
            inv.data[0]=pow(a.data[0],MOD-2);
            inv.t=dep-1;
            return;
        }
        Inv(a,inv,(dep+1)>>1,len);
        static Poly tmp1,tmp2,tmp;
        while((dep<<1)>len)
            len<<=1;
        for(int i=0;i<dep;i++)
            tmp.data[i]=a.data[i];
        for(int i=dep;i<len;i++)
            tmp.data[i]=0;
        NTT(tmp,1,len);
        NTT(inv,1,len);
        for(int i=0;i<len;i++)
            inv.data[i]=1LL*inv.data[i]*((2-1LL*inv.data[i]*tmp.data[i])%MOD+MOD)%MOD;
        NTT(inv,0,len);
        for(int i=dep;i<len;i++)
            inv.data[i]=0;
        inv.t=dep-1;
    }
    void div(Poly a,Poly b,Poly &D,Poly &R){
        static Poly tmp1,tmp2;
        int Up=a.t-b.t+1,midlen=1;
        tmp1=b;
        rever(tmp1);
        Inv(tmp1,tmp2,Up,midlen);
        tmp1=a;
        rever(tmp1);
        mul(tmp2,tmp1);
        save(tmp2,a.t-b.t);
        rever(tmp2);
        D=tmp2;
        mul(tmp2,b);
        for(int i=0;i<b.t;i++)
            R.data[i]=(a.data[i]-tmp2.data[i]+MOD)%MOD;
        R.t=b.t-1;
    }
    void sqrt(Poly a,Poly &b,int &midlen,int dep){
        // printf("dep=%lld
    ",dep);
        if(dep==1){
            b.data[0]=1;
            b.t=dep-1;
            return;
        }
        sqrt(a,b,midlen,(dep+1)>>1);
        // printf("dep=%lld
    ",dep);
        while((dep<<1)>(midlen))
            midlen<<=1;
        static Poly tmp1,tmp2,tmp3;
        tmp1=b;tmp3=b;
        save(tmp1,dep-1);
        save(tmp3,dep-1);	
        save(tmp2,-1);
        int midlent=1;
        for(int i=0;i<dep;i++)
            tmp1.data[i]=(tmp1.data[i]*2)%MOD;
        Inv(tmp1,tmp2,dep,midlent);
        mul(b,tmp3);
        for(int i=0;i<dep;i++)
            b.data[i]=(b.data[i]+a.data[i])%MOD;
        mul(b,tmp2);
        for(int i=dep;i<midlen;i++)	
            b.data[i]=0;		
        b.t=dep-1;
        // output(b);
    }
    Poly a,b;
    signed main(){
        scanf("%lld",&n);
        for(int i=0;i<n;i++)
            scanf("%lld",&a.data[i]);
        a.t=n-1;
        int midlen=1;
        sqrt(a,b,midlen,n);
        for(int i=0;i<=b.t;i++)
            printf("%lld ",b.data[i]);
        putchar('
    ');
        return 0;
    }
    
  • 相关阅读:
    3D集合图元:最小边界框/包围盒(boundingbox)
    vs2012下 error4996
    将自己的类封装为lib的方法
    3D特征:关于HFM和HBB
    C++的Matlab接口
    BigDataMini导论
    vs2012编译boost_1_54_0
    RGB_D_开发征程(使用Kinect)
    **PCD数据获取:Kinect+OpenNI+PCL对接(代码)
    PCL:全程详解 VS2010+PCL配置
  • 原文地址:https://www.cnblogs.com/dreagonm/p/10651696.html
Copyright © 2011-2022 走看看