给出n个数qi,给出Fj的定义如下:
令Ei=Fi/qi,求Ei.
标签:color mem lag 构造 memset 递增 技术分享 hint log
n行,第i行输出Ei。与标准答案误差不超过1e-2即可。
数学 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 }
标签:color mem lag 构造 memset 递增 技术分享 hint log
原文地址:http://www.cnblogs.com/SilverNebula/p/6516134.html