[精准]圆周率

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