1. 程式人生 > >BZOJ 1213, 高精度開根

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……霧……)

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; }