首页 > 代码库 > Bzoj3527 [Zjoi2014]力

Bzoj3527 [Zjoi2014]力

Time Limit: 30 Sec  Memory Limit: 256 MBSec  Special Judge
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

Sample Output

-16838672.693
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]力