首页 > 代码库 > BZOJ3527: [Zjoi2014]力

BZOJ3527: [Zjoi2014]力

给出n个数qi,给出Fj的定义如下: 
 
令Ei=Fi/qi,求Ei.

题解:

为什么要乘qi再除以qi。。。

我们发现i-j+j=i为定值,也就是这是个卷积。可以分开算,前半部分

令 a[i]=1/(i*i),b[i]=q[i],然后它们的卷积c[i]就是前半部分的值。

后半部分我们直接把b翻转一下再来一次就可以了。

注意不要给a添加无用的和不存在的项。

代码:比原来的pascal快了10s 233

 

  1 #include<cstdio>  2   3 #include<cstdlib>  4   5 #include<cmath>  6   7 #include<cstring>  8   9 #include<algorithm> 10  11 #include<iostream> 12  13 #include<vector> 14  15 #include<map> 16  17 #include<set> 18  19 #include<queue> 20  21 #include<string> 22  23 #define inf 1000000000 24  25 #define maxn 500000+5 26  27 #define maxm 20000000+5 28  29 #define eps 1e-10 30  31 #define ll long long 32  33 #define pa pair<int,int> 34  35 #define for0(i,n) for(int i=0;i<=(n);i++) 36  37 #define for1(i,n) for(int i=1;i<=(n);i++) 38  39 #define for2(i,x,y) for(int i=(x);i<=(y);i++) 40  41 #define for3(i,x,y) for(int i=(x);i>=(y);i--) 42  43 #define mod 1000000007 44  45 using namespace std; 46  47 inline int read() 48  49 { 50  51     int x=0,f=1;char ch=getchar(); 52  53     while(ch<0||ch>9){if(ch==-)f=-1;ch=getchar();} 54  55     while(ch>=0&&ch<=9){x=10*x+ch-0;ch=getchar();} 56  57     return x*f; 58  59 } 60 struct cp 61 { 62     double x,y; 63     cp operator +(cp b){return (cp){x+b.x,y+b.y};} 64     cp operator -(cp b){return (cp){x-b.x,y-b.y};} 65     cp operator *(cp b){return (cp){x*b.x-y*b.y,x*b.y+y*b.x};} 66 }; 67 double d[maxn],ans[maxn]; 68 //const double PI=3.1415926535898; 69 const double PI=acos(-1.0); 70 cp a[maxn],b[maxn],c[maxn],y[maxn]; 71 int n,nn,m,len,rev[maxn]; 72 void fft(cp *x,int n,int flag) 73 { 74     for0(i,n-1)y[rev[i]]=x[i]; 75     for0(i,n-1)x[i]=y[i]; 76     for(int m=2;m<=n;m<<=1) 77     { 78         cp wn=(cp){cos(2.0*PI/m*flag),sin(2.0*PI/m*flag)}; 79         for(int i=0;i<n;i+=m) 80         { 81             cp w=(cp){1,0};int mid=m>>1; 82             for0(j,mid-1) 83             { 84                 cp u=x[i+j],v=x[i+j+mid]*w; 85                 x[i+j]=u+v;x[i+j+mid]=u-v; 86                 w=w*wn; 87             } 88         } 89     } 90     if(flag==-1)for0(i,n-1)x[i].x/=n; 91 } 92  93 int main() 94  95 { 96  97     freopen("input.txt","r",stdin); 98  99     freopen("output.txt","w",stdout);100 101     nn=n=read();102     n=2*n-1;m=1;103     while(m<=n)m<<=1,len++;n=m;104     for0(i,nn-1)scanf("%lf",&d[i]);105     for0(i,n-1)a[i]=(cp){d[i],0};106     for0(i,nn-1)b[i]=i?(cp){1.0/((double)i*(double)i),0}:(cp){0,0};107     for2(i,nn,n-1)b[i]=(cp){0,0};108     for0(i,n-1)109     {110         int x=i,y=0;111         for1(j,len)y<<=1,y|=x&1,x>>=1;112         rev[i]=y;113     }114     fft(a,n,1);fft(b,n,1);115     for0(i,n-1)c[i]=a[i]*b[i];116     fft(c,n,-1);117     for0(i,nn-1)ans[i]=c[i].x;118     for0(i,nn-1)a[i]=(cp){d[nn-1-i],0};119     for2(i,nn,n-1)a[i]=(cp){0,0};120     for0(i,nn-1)b[i]=i?(cp){1.0/((double)i*(double)i),0}:(cp){0,0};121     for2(i,nn,n-1)b[i]=(cp){0,0};122     fft(a,n,1);fft(b,n,1);123     for0(i,n-1)c[i]=a[i]*b[i];124     fft(c,n,-1);125     for0(i,nn-1)ans[i]-=c[nn-1-i].x;126     for0(i,nn-1)printf("%.5f\n",ans[i]);127 128     return 0;129 130 }  
View Code

 

 

 

BZOJ3527: [Zjoi2014]力