[精准]圆周率
/* PI Calculator Author: stdafx */ # define fastcall __attribute__((optimize("-O3"))) # define IL __inline__ __attribute__((always_inline)) # define rep(_i,_s,_t) for(register int (_i)=(_s);(_i)<=(_t);++(_i)) # define dwn(_i,_s,_t) for(register int (_i)=(_s);(_i)>=(_t);--(_i)) # include <time.h> # include <math.h> # include <string> # include <stdio.h> # include <stdlib.h> # include <assert.h> # include <iostream> # include <string.h> # include <algorithm> # define int_log 31ll # define max_used_len 2097152ll fastcall IL int intSign(long long x) { return x<0; } fastcall IL double global_time() { return (double)clock() / CLOCKS_PER_SEC; } fastcall IL int rand32() { return (int)( (rand() & 0xf) << 0 | (rand() & 0xf) << 8 | (rand() & 0xf) << 16 | (rand() & 0xf) << 24 ) % 1000000000; } fastcall IL std::string number_to_string(long long number) { char c_str[30]; sprintf(c_str,"%lld",number); std::string str = c_str; return str; } typedef double ld; int bin[int_log+1]; struct pdd { int a,b,c; }; struct cpx { ld r, i; fastcall IL cpx conj() { return (cpx){r,-i}; } }; cpx *tf1,*tf2; fastcall cpx operator + (const cpx &a, const cpx &b) { return (cpx) { a.r + b.r, a.i + b.i}; } fastcall cpx operator - (const cpx &a, const cpx &b) { return (cpx) { a.r - b.r, a.i - b.i}; } fastcall cpx operator * (const cpx &a, const cpx &b) { return (cpx) { a.r*b.r - a.i*b.i, a.r*b.i + a.i*b.r}; } fastcall cpx rand_cpx() { return (cpx){ (double)(rand32()%1000), (double)(rand32()%1000) }; } const ld eps=1e-9; void write_in_file(const char *name ,const std::string &str){ FILE *file = fopen(name, "wb"); fwrite(str.c_str(), 1, str.size(), file); fclose(file); return; } namespace FFT { # define max_fft_log 23ll # define max_fft_len 8388608ll # define M_PI 3.14159265358979323846 int *rev[max_fft_log+1]; cpx *prw[max_fft_log+1],*prwi[max_fft_log+1]; bool prepared_Rev[max_fft_log+1],prepared_W[max_fft_log+1]; void prepareRev(int lg) { register int *tmp,len=bin[lg]; rev[lg] = new int [len]; tmp=rev[lg]; tmp[0]=0; for(register int i=0;i<len;i++) tmp[i]=tmp[i>>1]>>1|((i&1)<<(lg-1)); prepared_Rev[lg]=true; return; } void prepareWn(int lg,int p) { register cpx *pw,*pwi; prw [lg] = new cpx [p]; pw = prw [lg]; prwi [lg] = new cpx [p]; pwi = prwi [lg]; register ld omega = M_PI/p,s,c; for(register int i=0;i<p;i++) { pw [i] = (cpx){ c = cos(omega*i) , s = sin(omega*i) }; pwi [i] = (cpx){ c , - s }; } prepared_W[lg] = true; return; } void fft_normal(cpx *a,int n, int flag) { // to compare for (int i = 0, j = 0; i < n; i++) { if (i > j) std::swap(a[i], a[j]); for (int t = n >> 1; (j ^= t) < t; t >>= 1); } int m = 2; cpx wn, w, t, u; for (; m <= n; m <<= 1) { wn = (cpx) {cos(2.0 * M_PI / m), flag * sin(2.0 * M_PI / m)}, w = (cpx) {1, 0}; for (int i = 0; i < n; i += m, w = (cpx) {1, 0}) for (int k = 0, p = m >> 1; k < p; ++k, w = w * wn) t = w * a[i + k + p], u = a[i + k], a[i + k] = u + t, a[i + k + p] = u - t; } if (flag == -1) for (int i = 0; i < n; i++) a[i].r /= n; return; } void fft(register cpx *f,register int lg) { register int n= 1 << lg ; if(!prepared_Rev[lg]) prepareRev(lg); register int *tRev=rev[lg],i,m,log_m,k,p; register cpx *w,*f1,*f2,t; for(i=0;i<n;++i,++tRev) if(i<*tRev) std::swap(f[i],f[*tRev]); for(m=2,log_m=1;m<=n;m<<=1,++log_m) { p=m>>1; if(!prepared_W[log_m]) prepareWn(log_m,p); for(i=0;i<n;i+=m) { w=prw[log_m]; f1=f+i,f2=f+i+p; for(k=0;k<p;++k) { t=(*w)*(*f2); (*f2)=(*f1)-t; f1->r+=t.r; f1->i+=t.i; ++f1,++f2,++w; } } } return; } void ifft(register cpx *f,register int lg) { register int n= 1 << lg ; register int *tRev=rev[lg],i,m,log_m,k,p; register cpx *w,*f1,*f2,t; for(i=0;i<n;++i,++tRev) if(i<*tRev) std::swap(f[i],f[*tRev]); for(m=2,log_m=1;m<=n;m<<=1,++log_m) { p=m>>1; for(i=0;i<n;i+=m) { w=prwi[log_m]; f1=f+i,f2=f+i+p; for(k=0;k<p;++k) { t=(*w)*(*f2); (*f2)=(*f1)-t; f1->r+=t.r; f1->i+=t.i; ++f1,++f2,++w; } } } register ld img = 1.0 / n; for(i=0;i<n;i++) f[i].r*=img; return; } void Test() { double time_tot=0; register cpx *array = new cpx [max_fft_len]; register ld *rp = new double [max_fft_len]; for(int i=0;i<=max_fft_log;i++) { int len= 1 << i ; for(int j=0;j<len;j++) array[j]=(cpx) { rp[j] = (double)(rand32()%1000), 0 }; double st=global_time(); fft(array,i); ifft(array,i); printf("time used = %.3lf \n",global_time()-st); time_tot+=global_time()-st; for(int j=0;j<len;j++) assert( fabs(rp[j]-array[j].r) < eps); printf("accepted on len %d \n\n",len); } delete []array; printf("total time used = %.3lf \n",time_tot); return; } } # define EXTRA_PR 2 struct __float256 { int *array; // digits , 9 digits in one position int exp,L; // exp , total digits | if number = 0 , L = 0 bool sign; /* if sign = true , number >= 0 sign = false , number < 0 */ inline __float256() { // Initialization , 0 array = NULL; exp=L=0; sign=true; return; } fastcall IL void apply_memory(int len) { assert(array==NULL); array = new int [len]; memset(array,0,sizeof(int)*len); return; } fastcall IL void destroy() { if(array!=NULL) delete []array; array = NULL; exp=L=0; sign=true; return; } fastcall IL void set32(int number) { // 0 -> 10^9 - 1 destroy(); if(!number) return; L = 1 , apply_memory(1); if(number<0) sign=false,number=-number; else sign=true; array[0]=number; return; } fastcall IL void negate() { if(!L) return; sign^=1; return; } fastcall IL int at(int mag) { // word at (10 ^ 9) ^ mag if(mag<exp) return 0; if(mag>=L+exp) return 0; return array[ mag - exp ]; } int to_string_trimmed(int digits , std::string &str) { if(!L) { str="0"; return 0; } int exponent = exp; int length = L , *ptr = array ; if(!digits) { // all digits = length * 9; } else { int words = (digits + 17) / 9; if (words < length){ int chop = length - words; exponent += chop; length = words; ptr += chop; } } exponent *= 9; char buffer[10]; str.clear(); int c=length; while(c-- > 0) { register int word = ptr[c]; for(register int i=8;~i;--i) { buffer[i]=word%10+'0'; word/=10; } buffer[9]='\0'; str+=buffer; } int ledz=0; while(str[ledz]=='0') ++ledz; digits+=ledz; if(digits<str.size()) { exponent+=str.size()-digits; str.resize(digits); } return exponent; } std::string to_sci(int digits = 0) { if(!L) return "0"; std::string str; int exponent = to_string_trimmed(digits,str); int ledz=0; while(str[ledz]=='0') ++ledz; str=&str[ledz]; exponent+=str.size()-1; str=str.substr(0,1)+"."+(&str[1]); if(exponent!=0) { str+=" * 10 ^"; str+=number_to_string(exponent); } if(!sign) str="-"+str; return str; } std::string to_string(int digits) { if(!L) return "0"; int mag=exp+L; if(mag>1||mag<0) return to_sci(); // mag = 0 (-1) || mag = 1 ( 0 ) std::string str; int exponent = to_string_trimmed(digits,str); if(mag==0) { if(sign) return "0."+str; else return "-0."+str; } std::string before = number_to_string(array[L-1]); if(exponent>=0) { if(sign) return before; else return "-"+before; } std::string after=str.substr(str.size()+exponent,-exponent); if(sign) return before+"."+after; else return "-"+before+"."+after; } void print(){ std::cout<<to_string(30)<<std::endl; return; } }; fastcall IL pdd getDiv(int x) { register int a,b,c; a=x%1000,x/=1000; b=x%1000,x/=1000; c=x%1000; return (pdd){a,b,c}; } fastcall IL int mergePdd(register int a,register int b,register int c) { return a+b*1000+c*1000000; } int cmp_unsigned(__float256 &x,__float256 &y) { // ignoring sign // 1 = (x<y) , 0 = (x=y) , -1 = (x>y) register int I1=x.exp+x.L,I2=x.exp+x.L; // max_d if(I1!=I2) return I1<I2?1:-1; register int mag=I1,minExp=std::min(x.exp,y.exp); for(register int i=mag;i>=minExp;--i) { I1=x.at(i), I2=y.at(i); if(I1!=I2) return I1<I2?1:-1; } return 0; } __float256 mul(__float256 &x,__float256 &y,int precision) { __float256 z; if(!x.L||!y.L) return z; if(!precision) precision=x.L+y.L; else precision+=EXTRA_PR; register int Aexp=x.exp,Bexp=y.exp,AL=x.L,BL=y.L; register int *pA=x.array,*pB=y.array; if(AL>precision) { int chop=AL-precision; AL=precision; Aexp+=chop; pA+=chop; } if(BL>precision) { int chop=BL-precision; BL=precision; Bexp+=chop; pB+=chop; } z.sign=(x.sign==y.sign); z.exp=Aexp+Bexp; z.L=AL+BL; z.apply_memory(z.L); register pdd dx,dy; int fft_len=1,lg=0,tot=0; for(;fft_len<z.L*3;fft_len<<=1,++lg); tf1=new cpx [fft_len]; for(register int i=0;i<z.L;i++) { dx=getDiv(i<AL?pA[i]:0); dy=getDiv(i<BL?pB[i]:0); tf1[tot++]=(cpx){(ld)dx.a,(ld)dy.a}; tf1[tot++]=(cpx){(ld)dx.b,(ld)dy.b}; tf1[tot++]=(cpx){(ld)dx.c,(ld)dy.c}; } while(tot!=fft_len) tf1[tot++]=(cpx){0,0}; FFT::fft(tf1,lg); tf2=new cpx [fft_len]; for(register int i=0,j;i<fft_len;i++) { j=(fft_len-1)&(fft_len-i); tf2[i]=(tf1[i]*tf1[i]-(tf1[j]*tf1[j]).conj())*(cpx){0,-0.25}; } delete []tf1; FFT::ifft(tf2,lg); register int len; int *temporary = new int[ len = z.L*3 ]; register long long inc=0; for(int i=0;i<len;i++) { temporary[i]=((long long)(tf2[i].r+0.5)+inc)%1000ll; inc=((long long)(tf2[i].r+0.5)+inc)/1000ll; } delete []tf2; while((!temporary[len-1])&&(!temporary[len-2])&&(!temporary[len-3])) len-=3,z.L--; for(register int i=0;i<z.L;i++) z.array[i] = mergePdd(temporary[i*3],temporary[i*3+1],temporary[i*3+2]); delete []temporary; return z; } void mul_to(__float256 &x,__float256 &y,int precision) { __float256 z=mul(x,y,precision); x.destroy(); x=z; return; } __float256 mul(__float256 &x,unsigned int y) { __float256 z; if(!y||!x.L) return z; register long long inc=0; z.L=x.L; z.sign=x.sign; z.apply_memory(z.L+1); z.exp=x.exp; for(int i=0;i<z.L;i++) { z.array[i]=(inc+1ll*x.array[i]*y)%1000000000ll; inc=(inc+1ll*x.array[i]*y)/1000000000ll; } if(inc!=0) z.array[z.L++]=(int)(inc); return z; } void mul_to(__float256 &x,long long y) { __float256 z=mul(x,y); x.destroy(); x=z; return; } __float256 add_unsigned(__float256 &x,__float256 &y,int precision) { __float256 z; int magA=x.exp+x.L,magB=y.exp+y.L; int top = std::max(magA,magB); int bot = std::min(x.exp,y.exp); int TL=top-bot; if(!precision) precision=TL; else precision+=EXTRA_PR; if(TL>precision) bot=top-precision,TL=precision; z.sign=true; z.exp=bot; z.L=TL; z.apply_memory(z.L+1); register long long inc=0; register int I1,I2; for(register int i=0;bot<top;++bot,++i) { I1=x.at(bot),I2=y.at(bot); z.array[i]=(I1+I2+inc)%1000000000; inc=(I1+I2+inc)/1000000000; } if(inc) z.array[z.L++]=inc; return z; } __float256 sub_unsigned(__float256 &x,__float256 &y,int precision) { __float256 z; int magA=x.exp+x.L,magB=y.exp+y.L; int top = std::max(magA,magB); int bot = std::min(x.exp,y.exp); int TL=top-bot; if(!precision) precision=TL; else precision+=EXTRA_PR; if(TL>precision) bot=top-precision,TL=precision; z.sign=true; z.exp=bot; z.L=TL; z.apply_memory(z.L+1); register long long inc=0; register int I1,I2; for(register int i=0;bot<top;++bot,++i) { I1=x.at(bot),I2=y.at(bot); z.array[i]=I1-I2+inc; inc=0; while(z.array[i]<0) z.array[i]+=1000000000,--inc; } while(z.L>0&&z.array[z.L-1]==0) z.L--; if(!z.L) { z.destroy(); } return z; } __float256 add(__float256 &x,__float256 &y,int precision) { __float256 z; if(x.sign==y.sign) { z=add_unsigned(x,y,precision); z.sign=x.sign; } else { if(cmp_unsigned(x,y)<=0) { // |x| > =|y| z=sub_unsigned(x,y,precision); z.sign=x.sign; } else { z=sub_unsigned(y,x,precision); // |x| < |y| z.sign=y.sign; } } return z; } void add_to(__float256 &x,__float256 &y,int precision) { __float256 z; z=add(x,y,precision); x.destroy(); x=z; return; } __float256 sub(__float256 &x,__float256 &y,int precision) { __float256 z; if(x.sign!=y.sign) { // printf("?[]\n"); // x.print(); // y.print(); z=add_unsigned(x,y,precision); z.sign=x.sign; } else { if(cmp_unsigned(x,y)<=0) { // |x| > =|y| z=sub_unsigned(x,y,precision); z.sign=x.sign; } else { z=sub_unsigned(y,x,precision); // |x| < |y| z.sign=x.sign^1; } } return z; } void sub_to(__float256 &x,__float256 &y,int precision) { __float256 z; z=sub(x,y,precision); x.destroy(); x=z; return; } __float256 rcp(__float256 &x,int precision) { assert(x.L); if(!precision) { precision = 3; int Aexp=x.exp,AL=x.L,*pA=x.array; if(AL>precision) { int chop=AL-precision; AL=precision; Aexp+=chop; pA+=chop; } double v=pA[0]; if(AL>=2) v+=pA[1]*1000000000.; if(AL>=3) v+=pA[2]*1000000000000000000.; v=1./v,Aexp=-Aexp; while(v<1000000000.) { v*=1000000000.; --Aexp; } long long v64=(long long)v; __float256 z; z.sign=x.sign; z.apply_memory(2); z.array[0]=v64%1000000000; z.array[1]=v64/1000000000; z.L=2; z.exp=Aexp; return z; } int s = (precision >> 1) + 1; if (precision == 1) s = 0; if (precision == 2) s = 1; __float256 z; z = rcp(x,s); __float256 one ; one.set32(1); __float256 e; e = mul(z,x,precision); sub_to(e,one,precision); mul_to(e,z,precision); sub_to(z,e,precision); e.destroy(); one.destroy(); return z; } __float256 div(__float256 &x,__float256 &y,int precision) { __float256 z; z=rcp(y,precision); mul_to(z,x,precision); return z; } void div_to(__float256 &x,__float256 &y,int precision) { __float256 z; z = div(x,y,precision); x.destroy(); x=z; return ; } void div_to2(__float256 &q,__float256 &p,int precision) { __float256 z; z = div(q,p,precision); p.destroy(); p=z; return ; } __float256 invsqrt(int x,int precision) { // ( r0^2 * x - 1 ) // r1 = r0 - (----------------) * r0 // ( 2 ) assert(x>0); if(!precision) { double v=1.0/sqrt((double)(x)); int exponent=0; while(v<1000000000.) { v*=1000000000.; --exponent; } long long v64=(long long)(v); __float256 z; z.sign=true; z.apply_memory(2); z.array[0]=v64%1000000000ll; z.array[1]=v64/1000000000ll; z.L=2; z.exp=exponent; return z; } int s = (precision >> 1) + 1; if (precision == 1) s = 0; if (precision == 2) s = 1; __float256 z; z = invsqrt(x,s); __float256 sp; sp = mul(z,z,precision); __float256 fig; fig.set32(1); mul_to(sp,x); sub_to(sp,fig,precision); fig.destroy(); mul_to(sp,500000000); --sp.exp; mul_to(sp,z,precision); sub_to(z,sp,precision); sp.destroy(); return z; } void prepareUtility() { bin[0]=1; for(int i=1;i<=int_log;i++) bin[i]=1<<i; return; } void pi_bsr(__float256 &P,__float256 &Q,__float256 &R,int a,int b,int precision,int iter) { assert(iter<20); assert(iter>=0); assert(a<b); assert(a>=0); assert(b>=0); if(b-a==1) { // P = (13591409 + 545140134 b)(2b-1)(6b-5)(6b-1) (-1)^b P.set32(b); mul_to(P,545140134ul); __float256 fig; fig.set32(13591409ul); add_to(P,fig,precision); fig.destroy(); mul_to(P,2*b-1); mul_to(P,6*b-5); mul_to(P,6*b-1); if ( b & 1 ) P.negate(); // Q = 10939058860032000 * b^3 Q.set32(b); __float256 Q2=mul(Q,Q,precision); mul_to(Q,Q2,precision); Q2.destroy(); mul_to(Q,26726400ul); mul_to(Q,409297880ul); // R = (2b-1)(6b-5)(6b-1) R.set32(2*b-1); mul_to(R,6*b-5); mul_to(R,6*b-1); return; } int mid=(a+b)>>1; __float256 P0, Q0, R0, P1, Q1, R1; pi_bsr(P0, Q0, R0, a, mid, precision,iter+1); pi_bsr(P1, Q1, R1, mid, b, precision,iter+1); __float256 I1=mul(P0,Q1,precision); __float256 I2=mul(P1,R0,precision); P0.destroy(); P1.destroy(); P = add(I1,I2,precision); I1.destroy(); I2.destroy(); Q = mul(Q0,Q1,precision); Q0.destroy(); Q1.destroy(); R = mul(R0,R1,precision); R0.destroy(); R1.destroy(); return; } int main() { prepareUtility(); int digits; scanf("%d",&digits); digits++; int precision = (digits + 8) / 9; double terms = (precision * 0.6346230241342037371474889163921741077188431452678) + 1; __float256 P,Q,R; pi_bsr(P,Q,R,0,(int)terms,precision,0); R.destroy(); R=mul(Q,13591409); add_to(P,R,precision); mul_to(Q,4270934400ul); div_to2(Q,P,precision); Q.destroy(); Q=invsqrt(10005,precision); mul_to(P,Q,precision); write_in_file("B",P.to_string(digits)); return 0; }