BZOJ 1213, 高精度開根
阿新 • • 發佈:2018-12-24
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;
}