• [Zjoi2014]力(FFT,卷积)


    [Zjoi2014]力(FFT,卷积)

    题意:给定\(n\)个点电荷,排在单位数轴上,求每个点的场强

    考虑每个\(i\)对于每个\(j\)的贡献,分析式子

    \(E=\cfrac{q_i}{(j-i)^2}\)

    \(f(x)=\sum q_ix^i\)

    \(g(x)=\sum a_ix^i,a_i=i<0?-\frac{1}{i^2}:\frac{1}{i^2}\)

    \(g(x)\)每一项\(x\)的指数其实是\(j-i\)的值

    \(f(x)\cdot g(x)\)即可,注意负数系数的话偏移一下

    是一个简单的作差卷积,适合作为\(FFT\)入门题

    #include<bits/stdc++.h>
    using namespace std;
    
    #define double long double
    
    #define reg register
    typedef long long ll;
    #define rep(i,a,b) for(reg int i=a,i##end=b;i<=i##end;++i)
    #define drep(i,a,b) for(reg int i=a,i##end=b;i>=i##end;--i)
    
    template <class T> inline void cmin(T &a,T b){ ((a>b)&&(a=b)); } 
    template <class T> inline void cmax(T &a,T b){ ((a<b)&&(a=b)); } 
    #define double long double
    
    char IO;
    int rd(){
    	int s=0,f=0;
    	while(!isdigit(IO=getchar())) if(IO=='-') f=1;
    	do s=(s<<1)+(s<<3)+(IO^'0');
    	while(isdigit(IO=getchar()));
    	return f?-s:s;
    }
    
    const double PI=acos(-1);
    
    const int N=(1<<19)+4,P=998244353;
    const int g=3;
    
    bool be;
    int n,m;
    struct Cp{
    	double x,y;
    	Cp(){}
    	Cp(double _x,double _y){ x=_x,y=_y; }
    	Cp operator + (const Cp t){ return Cp(x+t.x,y+t.y); }
    	Cp operator - (const Cp t){ return Cp(x-t.x,y-t.y); }
    	Cp operator * (const Cp t){ return Cp(x*t.x-y*t.y,x*t.y+y*t.x); }
    } a[N],b[N];
    int rev[N];
    
    void FFT(int n,Cp *a,int f) {
    	rep(i,0,n-1) if(rev[i]>i) swap(a[i],a[rev[i]]);
    	for(reg int i=1;i<n;i<<=1) {
    		Cp w(cos(PI/i),f*sin(PI/i));
    		for(reg int l=0;l<n;l+=i*2) {
    			Cp e(1,0);
    			for(reg int j=l;j<l+i;++j,e=e*w) {
    				Cp t=a[j+i]*e;
    				a[j+i]=a[j]-t;
    				a[j]=a[j]+t;
    			}
    		}
    	}
    	if(f==-1) rep(i,0,n-1) a[i].x/=n;
    }
    
    
    bool ed;
    int main(){
    	n=rd();
    	int R=1,c=-1;
    	while(R<=n*3) R<<=1,c++;
    	rep(i,1,R) rev[i]=(rev[i>>1]>>1)|((i&1)<<c);
    	rep(i,0,n-1) scanf("%Lf",&a[i].x);
    	rep(i,-n+1,n-1) if(i) {
    		b[i+n].x=1.0/(1.0*i*i);
    		if(i<0) b[i+n].x*=-1;
    	}
    	FFT(R,a,1),FFT(R,b,1);
    	rep(i,0,R) a[i]=a[i]*b[i];
    	FFT(R,a,-1);
    	rep(i,0,n-1) printf("%.3Lf\n",a[i+n].x);
    }
    
    
    
    
    
    
  • 相关阅读:
    JBPM工作流(四)——管理流程定义
    JBPM工作流(三)——ProcessEngine与Service API
    JBPM工作流(二)——数据库表说明
    JBPM工作流(一)——实现一个简单的工作流例子
    jbpm与spring hibernate struts整合
    SpringMVC12拦截器
    SpringMVC11文件上传
    阅读代码的方法
    关于linux系统的资料
    关于图灵机的介绍(相见恨晚,太赞了)
  • 原文地址:https://www.cnblogs.com/chasedeath/p/12092758.html
Copyright © 2020-2023  润新知