首页 > 代码库 > 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 }
BZOJ3527: [Zjoi2014]力
声明:以上内容来自用户投稿及互联网公开渠道收集整理发布,本网站不拥有所有权,未作人工编辑处理,也不承担相关法律责任,若内容有误或涉及侵权可进行投诉: 投诉/举报 工作人员会在5个工作日内联系你,一经查实,本站将立刻删除涉嫌侵权内容。