BZOJ 1213, 高精度开根

Problem

传送门

Mean

输入一个正整数m（1≤m≤50）和一个整数n（0≤n≤10^10000），求开根后取整的结果。

Analysis

二分+高精度+NTT
倍增地确定二分边界L，R；
由于n过大，高精度乘法不够快，所以需要用NTT来加速乘法（FFT比较慢，而且精度不够）。
（但是我的NTT貌似打萎了，m过小的时候就会出错……不知什么……好在m小的时候直接上高精乘也是可以的……于是…………大牛要是看出我的NTT哪里错了请不吝赐教QAQ）
（其实是第一次打NTT……最后的感想是人生苦短，请用Python……雾……）

FFT/NTT参考资料

Code

#include<cstdio>
#include<string>
#include<cstring>
#include<iostream>
using namespace std;
typedef long long ll;
const int P=998244353,G=3,N=16384,K=13;
int m,i,inv[N+1],n,g[K+1],ng[K+1],x[N+5],y[N+5];
string s;
struct BigInt{
    static const int BASE=10000;
    static const int WIDTH=4;
    int s[10010];
    BigInt(int num=0){*this=num;}
    BigInt(string num){*this=num;}
    BigInt operator = (int& num){
        s[0]=0;
        while(num){
            s[++s[0]]=num%BASE;
            num/=BASE;
        }
        return *this;
    }
    BigInt operator = (const string& str){
        s[0]=0;
        int len=(str.length()-1)/WIDTH+1;
        for(int i=0;i<len;i++){
            int end=str.length()-i*WIDTH,start=end>WIDTH?end-WIDTH:0;
            sscanf(str.substr(start,end-start).c_str(),"%d",&s[++s[0]]);
        }
        return *this;
    }
    void NTT(int a[],int n,int t){
        for(int i=1,j=0;i<n-1;i++){
            for(int s=n;j^=s>>=1,~j&s;);
            if(i<j){int tmp=a[i];a[i]=a[j];a[j]=tmp;}
        }
        for(int d=0;(1<<d)<n;d++){
            int m=1<<d,m2=m<<1,_w=t>0?g[d]:ng[d];
            for(i=0;i<n;i+=m2)
              for(int w=1,j=0;j<m;j++){
                int &A=a[i+j+m],&B=a[i+j],t=(ll)w*A%P;
                A=B-t;if(A<0) A+=P;
                B+=t;if(B>=P) B-=P;
                w=(ll)w*_w%P;
              }
        }
        if(t<0) for(i=0;i<n;i++) a[i]=(ll)a[i]*inv[n]%P;
    }
    BigInt operator * (BigInt& b) {
        BigInt c;
        memset(c.s,0,sizeof(c.s));
        if(m>20){
            int n=1,Max=s[0]>b.s[0]?s[0]:b.s[0];
            while(n<(Max<<1)) n<<=1;
            for(i=0;i<s[0];i++) x[i]=s[i+1];
            for(i=0;i<b.s[0];i++) y[i]=b.s[i+1];
            NTT(x,n,1),NTT(y,n,1);
            for(int i=0;i<n;i++) x[i]=(ll)x[i]*y[i]%P;
            NTT(x,n,-1);
            for(int i=0;i<n;i++){
                c.s[++c.s[0]]+=x[i];
                c.s[c.s[0]+1]+=c.s[c.s[0]]/BASE;
                c.s[c.s[0]]%=BASE;
            }
            for(i=0;i<n;i++) x[i]=y[i]=0;
        }else{
            c.s[0]=s[0]+b.s[0];
            for(i=1;i<=s[0];i++)
            for(int j=1,g=0;;j++){
                if(j>b.s[0] && !g) break;
                int u=g;
                if(j<=b.s[0]) u+=s[i]*b.s[j];
                c.s[i+j-1]+=u%BASE;
                g=u/BASE;
                c.s[i+j]+=c.s[i+j-1]/BASE;
                c.s[i+j-1]%=BASE;
            }
        }
        while(!c.s[c.s[0]] && c.s[0]) c.s[0]--;
        return c;
    }
    BigInt operator *= (BigInt& b){return *this=*this*b;}
    BigInt operator + (const BigInt& b) const {
        BigInt c;
        c.s[0]=0;
        for(int i=1,g=0;;i++){
            if(i>s[0] && i>b.s[0] && !g) break;
            int x=g;
            if(i<=s[0]) x+=s[i];
            if(i<=b.s[0]) x+=b.s[i];
            c.s[++c.s[0]]=x%BASE;
            g=x/BASE;
        }
        return c;
    }
    BigInt operator - (const BigInt& b){
        BigInt c;
        c.s[0]=s[0];
        for(int i=1;i<=s[0];i++){
            if(i<=b.s[0]){
                if(s[i]<b.s[i]) s[i]+=BASE,s[i+1]--;
                c.s[i]=s[i]-b.s[i];
            }else{
                if(s[i]<0) s[i]+=BASE,s[i+1]--;
                c.s[i]=s[i];
            }
        }
        while(!c.s[c.s[0]] && c.s[0]) c.s[0]--;
        return c;
    }
    BigInt operator / (const int& b) const {
        BigInt c;
        c.s[0]=s[0];
        for(int i=s[0],x=0;i>0;i--){
            x=x*BASE+s[i];
            c.s[i]=x/b;
            x%=b;
        }
        while(!c.s[c.s[0]] && c.s[0]) c.s[0]--;
        return c;
    }
};
int cmp(const BigInt& a,const BigInt& b){
    if(a.s[0]!=b.s[0]) return a.s[0]<b.s[0];
    for(int i=a.s[0];i;i--) if(a.s[i]!=b.s[i]) return a.s[i]<b.s[i];
    return -1;
}
int pow(int a,int n){
    int ans=1;
    for(;n;n>>=1,a=(ll)a*a%P) if(n&1) ans=(ll)ans*a%P;
    return ans;
}
BigInt pow(BigInt a,int n){
    BigInt ans=1;
    bool p=1;
    for(;n;n>>=1,a*=a){
        if(n&1){
            if(p) ans=a,p=0;
            else ans*=a;
        } 
        if(n==1) break;
    }
    return ans;
}
int main(){
    for(g[K]=pow(G,(P-1)/N),ng[K]=pow(g[K],P-2),i=K-1;~i;i--)
        g[i]=(ll)g[i+1]*g[i+1]%P,ng[i]=(ll)ng[i+1]*ng[i+1]%P;
    for(inv[1]=1,i=2;i<=N;i++) inv[i]=(ll)(P-inv[P%i])*(P/i)%P;
    scanf("%d",&m);
    cin>>s;
    if(m==1){cout<<s;return 0;}
    if(s=="0"){printf("0");return 0;}
    BigInt v=s,Ans,l=1,L=1,r;
    for(;;){
        L.s[L.s[0]]=0;
        L.s[L.s[0]+=m]=1;
        if(cmp(L,v)==1) l.s[l.s[0]]=0,l.s[++l.s[0]]=1;
        else break;
    }
    r.s[0]=l.s[0]+1,r.s[r.s[0]]=1;
    for(int u=cmp(l,r);u==-1 || u==1;u=cmp(l,r)){
        BigInt mid=(l+r)/2,ans=1;
        ans=pow(mid,m);
        int p=cmp(ans,v);
        Ans=mid;
        if(p==-1) break;
        else if(p==1){
            l=mid+1;
            if(!cmp(pow(l,m),v)) break;
        }else{
            BigInt tmp=mid;
            r=tmp-1;
        }
    }
    printf("%d",Ans.s[Ans.s[0]]);
    for(i=Ans.s[0]-1;i>0;i--) printf("%04d",Ans.s[i]);
    return 0;
}