Íŕçŕä ę ëó÷řčě đĺřĺíč˙ě Status: AC, problem ZSQRT2, contest ZEL07. By levlam (Levin Aleksey), 2007-02-15 10:23:55.
#include <stdio.h>
#include <math.h>
#include <memory.h>
#define base 10000
#define N 2000010
#define MaxN N/4
#define S2 0.7071067811865475244008443621
#define doFHT(x, y)\
 x=new double[newl];\
 memset(x, 0, newl<<3);\
 for (i=0; i<l; i++) x[i]=y[i];\
 FHT(x, newl);
#define FHTB(i1,i2,C,S)\
 t1=L[i1]-R[i1], t2=L[i2]-R[i2];\
 L[i1]+=R[i1];\
 L[i2]+=R[i2];\
 R[i1]=t1*C+t2*S;\
 R[i2]=t1*S-t2*C;
 
int a[4*MaxN+10], b[2*MaxN+10], c[2*MaxN+10], d[4*MaxN+10], z[4*MaxN+10];
double st[80];
 
void Sum(int *x, int *y, int *z, int l1, int l2)
{
int *p=y, *p1=y+l1, *p2=z, *p3=x, *p4=z+l2, um=0;
for (; p<p1 && p2<p4; p++, p2++, p3++)
{
 *p3+=um+*p+*p2;
 for (um=0; *p3>=base; um++) *p3-=base;
}
while (p<p1)
{
 *p3+=um+*p;
 um=0;
 if (*p3>=base) um=1, *p3-=base;
 p++, p3++;
}
while (p2<p4)
{
 *p3+=um+*p2;
 um=0;
 if (*p3>=base) um=1, *p3-=base;
 p2++, p3++;
}
while (um>0)
{
 *p3+=um;
 um=0;
 if (*p3>=base) um=1, *p3-=base;
 p3++;
}
}
 
void Umen(int *x, int *y, int l)
{
int *p=x, *p1=y, *p2=x+l, t;
for (; p<=p2; p1++)
{
 t=*p-*p1;
 if (t>=0) *p++=t;
 else
  *p=base+t, (*(++p))--;
}
while ((*p)<0)
 *p+=base, (*(++p))--;
}
 
void FHT(double *a, int l)
{
if (l==4)
{
 double d0=a[0], d1=a[1], d2=a[2], d3=a[3];
 a[0]=d0+d1+d2+d3;
 a[1]=d0-d1+d2-d3;
 a[2]=d0+d1-d2-d3;
 a[3]=d0-d1-d2+d3;
 return;
}
l>>=1;
int x, l2=l>>1, l4=l>>2;
for (x=-1; l; x++) l>>=1;
l=l2<<1;
double *L=a, *R=a+l, t1, t2, t, S0=st[x], C0=st[x+1], S, C;
C0=-2.0*C0*C0;
S=S0, C=1.0+C0;
t1=L[0], t2=R[0];
L[0]=t1+t2, R[0]=t1-t2;
t1=L[l2], t2=R[l2];
L[l2]=t1+t2, R[l2]=t1-t2;
for (x=1;x<l4;x++)
{
 FHTB(x,l-x,C,S);
 FHTB(l2-x,l2+x,S,C);
 t=C;
 C=C*C0-S*S0+C;
 S=S*C0+t*S0+S;
}
FHTB(l4,l-l4,S2,S2);
FHT(L,l);
FHT(R,l);
}
 
void FHTReOrder(double *a, int l)
{
double t;
int i, tt, r=0;
if (l<=2) return;
for (i=1; i<l; i++)
{
 tt=l;
 do
 {
  l>>=1;
  r^=l;
 }while ((r&l)==0);
 l=tt;
 if (r>i) t=a[i], a[i]=a[r], a[r]=t;
}
}
 
void Umn(int *x, int *y, int *z, int l)
{
int newl=8, i, j, k;
double *a, *b, ad, sd, t, um=0;
while (newl<l+l) newl<<=1;
doFHT(a,y);
doFHT(b,z);
a[0]*=b[0]+b[0];
a[1]*=b[1]+b[1];
for (k=2; k<newl; k*=2)
 for (i=k, j=k+k-1; j>=k; i+=2, j-=2)
 {
  ad=a[i]+a[j];
  sd=a[i]-a[j];
  a[i]=b[i]*ad+b[j]*sd;
  a[j]=b[j]*ad-b[i]*sd;
 }
FHTReOrder(a,newl);
FHT(a,newl);
FHTReOrder(a,newl);
for (j=0; j<l+l; j++)
{
 t=floor(a[j]*0.5/newl+um+x[j]+0.5);
 um=floor(t/base);
 x[j]=t-um*base;
}
while (um!=0)
{
 t=floor(um+x[j]+0.5);
 um=floor(t/base);
 x[j++]=t-um*base;
}
delete[] a;
delete[] b;
return;
}
 
int main(void)
{
int i, k, l=1, deg, dd, t;
long long z0;
a[0]=3;
b[0]=2;
for (i=0, k=1; k<=(N<<3); k*=2)
 st[i++]=sin(3.1415926535897932384626433832/k);
while (2*l<MaxN)
{
 Umn(d,a,b,l);
 Sum(d,d,NULL,(l<<1)+2,0);
 Umn(c,b,b,l);
 Sum(c,c,NULL,(l<<1)+2,0);
 Umn(c,a,a,l);
 l=(l<<1)+2;
 while (c[l-1]==0 && d[l-1]==0) l--;
 memcpy(a, c, l<<2);
 memcpy(b, d, l<<2);
 memset(c, 0, l<<2);
 memset(d, 0, l<<2);
}
for (deg=1, dd=0; deg<MaxN; deg<<=1, dd++) ;
deg+=3;
for (i=deg-1; i-deg+l>=0 ; i--)
 b[i]=b[i-deg+l];
memset(b, 0, (deg-l)<<2);
while (b[deg-1]==0) deg--;
z0=1.0*base*base*base*base*base/((1.0*b[deg-1]*base+b[deg-2])*base+b[deg-3]);
z[2]=(z0/base)/base;
z[1]=(z0/base)%base;
z[0]=z0%base;
for (k=0, t=1; k<dd; k++)
{
 Umn(c,z,z,t+2);
 if (k<dd-1)
  Umn(d,c,b+(deg-(2*t+3)),2*t+4);
 else
  Umn(d+2*(2*t+2-MaxN),c+(2*t+2-MaxN),b+(deg-MaxN-1), MaxN+2);
 Sum(z+(t+t+t+4),z,z,t+2,t+2);
 memset(z, 0, (t+2)<<2);
 Umen(z, d, (t<<2)+7);
 memcpy(z, z+(t+t+4), (t+t+2)<<2);
 memset(z+(t+t+2), 0, (t+t+5)<<2);
 c[0]=1;
 Sum(z,c,NULL,1,0);
 memset(c, 0, (t+t+4)<<2);
 memset(d, 0, (4*t+6)<<2);
 t<<=1;
}
Umn(c, a, z+(t+2-MaxN), MaxN+1);
deg=MaxN+MaxN+2;
while (c[deg]==0) deg--;
if (c[deg]>=1000) printf("1.414");
else if (c[deg]>=100) printf("1.41");
else if (c[deg]>=10) printf("1.4");
else printf("1.");
for (i=1; i<MaxN-2; i++)
 printf("%04d", c[deg-i]);
printf("76");
return 0;
}