首页 > 代码库 > Bzoj3527 [Zjoi2014]力
Bzoj3527 [Zjoi2014]力
Submit: 1841 Solved: 1095
Description
给出n个数qi,给出Fj的定义如下:
令Ei=Fi/qi,求Ei.
Input
第一行一个整数n。
接下来n行每行输入一个数,第i行表示qi。
n≤100000,0<qi<1000000000
Output
n行,第i行输出Ei。与标准答案误差不超过1e-2即可。
Sample Input
5
4006373.885184
15375036.435759
1717456.469144
8514941.004912
1410681.345880
4006373.885184
15375036.435759
1717456.469144
8514941.004912
1410681.345880
Sample Output
-16838672.693
3439.793
7509018.566
4595686.886
10903040.872
3439.793
7509018.566
4595686.886
10903040.872
HINT
Source
数学 FFT
式子里的q[j]可以直接约掉。
前后两半式子看上去很相似,所以拆开考虑。
对于前半边:
分母上的(j-i)是等差递减的,离j越近就越小,分子上的qi下标是等差递增的,离j越近就越大。
↑多么像一个优美的卷积
然后构造出一个新多项式B,第i项的系数为1/(i+1)^2,和q[]数组做卷积,结果的第i位就是E[i+1]的前一半 (i下标从0开始);接着把q数组倒序,多项式B整体取负,做出来的结果E[i]就是E[n-(i+1)-1]的后一半。
↑数学多么优美
写完之后的半个小时里各种调试,就是过不了样例……然后想起来【第i位就是E[i+1]的前一半】,啊,需要把答案整体右移一位。最左边的E[1]当然没有前一半可以加,最右边的E[n]当前没有后一半可以减,所以这两位置0
然后就妥妥地过了
之后我在想,这个右移多麻烦啊。
确实呢,多项式B中,只要把第i项的系数设为1/i^2(下标从0开始,所以第0项为0),结果的第i位就是E[i]的结果咯
1 /*by SilverN*/ 2 #include<algorithm> 3 #include<iostream> 4 #include<cstring> 5 #include<cstdio> 6 #include<cmath> 7 #include<vector> 8 using namespace std; 9 const double pi=acos(-1.0);10 const int mxn=800010;11 struct com{12 double x,y;13 com operator + (com b){return (com){x+b.x,y+b.y};}14 com operator - (com b){return (com){x-b.x,y-b.y};}15 com operator * (com b){return (com){x*b.x-y*b.y,x*b.y+y*b.x};}16 }a[mxn],b[mxn],c[mxn];17 int n,l;18 int rev[mxn];19 void FFT(com *a,int flag){20 int i,j,k;21 for(i=0;i<n;i++)22 if(rev[i]>i)swap(a[rev[i]],a[i]);23 for(i=1;i<n;i<<=1){24 com wn=(com){cos(pi/i),flag*sin(pi/i)};25 for(j=0;j<n;j+=(i<<1)){26 com w=(com){1,0};27 for(k=0;k<i;k++,w=w*wn){28 com x=a[j+k],y=w*a[i+j+k];29 a[j+k]=x+y;30 a[i+j+k]=x-y;31 }32 }33 }34 if(flag==-1)for(i=0;i<n;i++) a[i].x/=n;35 return;36 }37 int len;38 double q[mxn];39 int main(){40 int i,j;41 scanf("%d",&len);42 for(i=0;i<len;i++)scanf("%lf",&q[i]);43 //44 int m=len<<1;45 for(n=1;n<m;n<<=1)l++;46 for(i=0;i<n;i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-1));47 //48 for(i=0;i<len;i++)a[i].x=q[i];49 for(i=1;i<len;i++)b[i].x=1.0/(double)(i)/(double)(i);50 //51 FFT(a,1);FFT(b,1);52 for(i=0;i<n;i++)c[i]=a[i]*b[i];53 //54 FFT(c,-1);55 56 57 memset(a,0,sizeof a);58 memset(b,0,sizeof b);59 for(i=0;i<len;i++)a[i].x=q[len-i-1];60 for(i=1;i<len;i++)b[i].x=-1.0/(double)(i)/(double)(i);61 62 FFT(a,1);FFT(b,1);63 for(i=0;i<n;i++)a[i]=a[i]*b[i];64 // FFT(c,-1);65 FFT(a,-1);66 /* 67 for(i=n-1;i;i--){68 c[i]=c[i-1];69 a[i]=a[i-1];70 }71 c[0]=a[0]=(com){0,0};72 */73 for(i=0;i<len;i++)c[i]=c[i]+a[len-i-1];74 for(i=0;i<len;i++){75 printf("%.3f\n",c[i].x);76 }77 return 0;78 }
Bzoj3527 [Zjoi2014]力
声明:以上内容来自用户投稿及互联网公开渠道收集整理发布,本网站不拥有所有权,未作人工编辑处理,也不承担相关法律责任,若内容有误或涉及侵权可进行投诉: 投诉/举报 工作人员会在5个工作日内联系你,一经查实,本站将立刻删除涉嫌侵权内容。