3527: [Zjoi2014]力

最新推荐文章于 2018-11-26 22:27:00 发布

转载最新推荐文章于 2018-11-26 22:27:00 发布 · 146 阅读

CC 4.0 BY-SA版权

原文链接：http://www.cnblogs.com/htwx/articles/5588849.html

本文深入解析了快速傅立叶变换（FFT）算法在处理特定数学问题中的应用，通过将复杂的数学公式转换为卷积操作，展示了如何利用FFT解决一系列涉及大量数值计算的问题。文章提供了详细的代码实现，包括数据预处理、FFT核心算法及结果后处理等关键步骤。

给出n个数q_i，给出F_j的定义如下：

令E_i=F_i/q_i，求E_i.

化简一下，

F_J=Σ((q_j)/(i-j)²)-Σ((q_j)/(i-j)²)

另t_i=1/i² ,那么t_i-j=1/(i-j)²。。。

那么Σ((q_j)/(i-j)²)就变成了一个卷积，后面的把数组翻转一下就好了。。。

 1 #include<iostream>
 2 #include<cstdlib>
 3 #include<cmath>
 4 #include<cstring>
 5 #include<cstdio>
 6 #include<algorithm>
 7 #include<string>
 8 #include<map>
 9 #include<queue>
10 #include<vector>
11 #include<set>
12 #define inf 1000000000
13 #define maxn 500000+5
14 #define maxm 10000+5
15 #define eps 1e-10
16 #define ll long long
17 #define for0(i,n) for(int i=0;i<=(n);i++)
18 #define for1(i,n) for(int i=1;i<=(n);i++)
19 #define for2(i,x,y) for(int i=(x);i<=(y);i++)
20 #define for3(i,x,y) for(int i=(x);i>=(y);i--)
21 #define for4(i,x) for(int i=head[x],y=e[i].go;i;i=e[i].next,y=e[i].go)
22 using namespace std;
23 int read(){
24     int x=0,f=1;char ch=getchar();
25     while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
26     while(ch>='0'&&ch<='9'){x=10*x+ch-'0';ch=getchar();}
27     return x*f;
28 }
29 struct cp{
30     double x,y;
31     cp operator +(cp b){return (cp){x+b.x,y+b.y};}
32     cp operator -(cp b){return (cp){x-b.x,y-b.y};}
33     cp operator *(cp b){return (cp){x*b.x-y*b.y,x*b.y+y*b.x};}
34 };
35 double d[maxn],ans[maxn];
36 const double PI=acos(-1.0);
37 cp a[maxn],b[maxn],c[maxn],y[maxn];
38 int n,nn,m,len,rev[maxn];
39 void fft(cp *x,int n,int flag){
40     for0(i,n-1)y[rev[i]]=x[i];
41     for0(i,n-1)x[i]=y[i];
42     for(int m=2;m<=n;m<<=1){
43         cp wn=(cp){cos(2.0*PI/m*flag),sin(2.0*PI/m*flag)};
44         for(int i=0;i<n;i+=m){
45             cp w=(cp){1.0,0};int mid=m>>1;
46             for0(j,mid-1){
47                 cp u=x[i+j],v=x[i+j+mid]*w;
48                 x[i+j]=u+v;x[i+j+mid]=u-v;
49                 w=w*wn;
50             }
51         }
52     }
53     if(flag==-1)for0(i,n-1)x[i].x/=n;
54 }
55 int main(){
56     //freopen("input.txt","r",stdin);
57     //freopen("output.txt","w",stdout);
58      nn=n=read();
59     n=2*n-1;m=1;
60     while(m<=n)m<<=1,len++;n=m;
61     for0(i,nn-1)scanf("%lf",&d[i]);
62     for0(i,n-1)a[i]=(cp){d[i],0};
63     for0(i,nn-1)b[i]=i?(cp){1.0/((double)i*(double)i),0}:(cp){0,0};
64     for2(i,nn,n-1)b[i]=(cp){0,0};
65     for0(i,n-1){
66         int x=i,y=0;
67         for1(j,len)y<<=1,y|=x&1,x>>=1;
68         rev[i]=y;
69     }
70     fft(a,n,1);fft(b,n,1);
71     for0(i,n-1)c[i]=a[i]*b[i];
72     fft(c,n,-1);
73     for0(i,nn-1)ans[i]=c[i].x;
74     for0(i,nn-1)a[i]=(cp){d[nn-1-i],0};
75     for2(i,nn,n-1)a[i]=(cp){0,0};
76     for0(i,nn-1)b[i]=i?(cp){1.0/((double)i*(double)i),0}:(cp){0,0};
77     for2(i,nn,n-1)b[i]=(cp){0,0};
78     fft(a,n,1);fft(b,n,1);
79     for0(i,n-1)c[i]=a[i]*b[i];
80     fft(c,n,-1);
81     for0(i,nn-1)ans[i]-=c[nn-1-i].x;
82     for0(i,nn-1)printf("%.5f\n",ans[i]);
83     return 0;
84 
85 }

View Code

转载于:https://www.cnblogs.com/htwx/articles/5588849.html