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
",c[i].x);
76     }
77     return 0;
78 }
原文地址:https://www.cnblogs.com/SilverNebula/p/6516134.html