luogu P3338 [ZJOI2014]力

最新推荐文章于 2024-02-28 22:19:43 发布

weixin_30379531

最新推荐文章于 2024-02-28 22:19:43 发布

阅读量45

点赞数

CC 4.0 BY-SA版权

原文链接：http://www.cnblogs.com/smyjr/p/10073912.html

本文介绍了一种新颖的方法来计算电荷系统的势能，通过将电荷分布转换为多项式并利用快速傅立叶变换（FFT）进行高效计算。这种方法避免了直接求和的复杂度，提供了更快速的解决方案。

传送门

首先化简原式\[F_j=\sum_{i<j}\frac{q_iq_j}{(i-j)^2}-\sum_{i>j}\frac{q_iq_j}{(i-j)^2},E_j=F_j/q_j\]

把所有\(q_j\)提出来,则显然\[E_j=\sum_{i<j}\frac{q_i}{(i-j)^2}-\sum_{i>j}\frac{q_i}{(i-j)^2}\]\[E_j=...-\frac{q_{j-2}}{2^2}-\frac{q_{j-1}}{1^2}+0+\frac{q_{j+1}}{1^2}+\frac{q_{j+2}}{2^2}...\]

然后设多项式\(A,B\),其中\(A[i]=q_{i+1}(0\le i<n)\)\(,B[i]=\frac{1}{(i-(n-1))|i-(n-1)|}(0\le i<2n-1)(B[n-1]=0)\),然后把\(A,B\)乘起来,答案要求的\(E_i\)也就是多项式的第\(n-2+i\)项的系数

#include<bits/stdc++.h>
#define LL long long
#define db double
#define il inline
#define re register
  
using namespace std;
const int N=100000+10,M=550000+10;
const db pi=acos(-1);
il int rd()
{
  int x=0,w=1;char ch=0;
  while(ch<'0'||ch>'9') {if(ch=='-') w=-1;ch=getchar();}
  while(ch>='0'&&ch<='9') {x=(x<<3)+(x<<1)+(ch^48);ch=getchar();}
  return x*w;
}
struct comp
{
  db r,i;
  comp(){r=i=0;}
  comp(db nr,db ni){r=nr,i=ni;}
  il comp operator + (const comp &bb) const {return comp(r+bb.r,i+bb.i);}
  il comp operator - (const comp &bb) const {return comp(r-bb.r,i-bb.i);}
  il comp operator * (const comp &bb) const {return comp(r*bb.r-i*bb.i,r*bb.i+i*bb.r);}
}a[M],b[M];
int n,m,nn,l,rdr[M];
void fft(comp *a,int op)
{
  comp W,w,x,y;
  for(int i=0;i<nn;++i) if(i<rdr[i]) swap(a[i],a[rdr[i]]);
  for(int i=1;i<nn;i<<=1)
    {
      W=comp(cos(pi/i),sin(pi/i)*op);
      for(int j=0;j<nn;j+=i<<1)
        {
          w=comp(1,0);
          for(int k=0;k<i;++k,w=w*W)
            {
              x=a[j+k],y=w*a[j+k+i];
              a[j+k]=x+y,a[j+k+i]=x-y;
            }
        }
    }
    
}
  
int main()
{
  n=rd();
  for(int i=0;i<n;++i) scanf("%lf",&a[i].r);
  m=-1;
  for(int i=-n+1;i<=n-1;++i) b[++m].r=1/(db)i/fabs(i);
  b[n-1].r=0;
  m+=n-1;
  for(nn=1;nn<=m;nn<<=1) ++l;
  for(int i=0;i<nn;++i) rdr[i]=(rdr[i>>1]>>1)|((i&1)<<(l-1));
  fft(a,1),fft(b,1);
  for(int i=0;i<nn;++i) a[i]=a[i]*b[i];
  fft(a,-1);
  for(int i=n-1;i<n+n-1;i++) printf("%.3lf\n",a[i].r/nn);
  return 0;
}

转载于:https://www.cnblogs.com/smyjr/p/10073912.html