3527: [Zjoi2014]力

此题bzoj未提供题意--

把定义式直接除以qi,右边的东西再弄成逆序,就是两个卷积,直接套FFT



第一次写FFT把不理解的地方就补在这里了

感觉Rader算法求倒位序很麻烦。。。网上很多人都是用预处理rev数组,个人感觉简单暴力

然后FFT中使用的多项式要补齐到2^k项,这样才能保证蝴蝶操作的正确性

此处2^k取一个刚好>=2*n的东西

以及代码中的h数组不能多赋值,,该是n项就是n项,不然点值表达式是错的

#include<iostream>  
#include<cstdio>  
#include<cstring>  
#include<algorithm>  
#include<cmath>  
using namespace std;  
   
const int maxn = 4E5 + 10;  
typedef double DB;  
const DB PI = acos(-1.0);  
   
struct Virt{  
    DB r,i;  
    Virt(){}  
    Virt(DB r,DB i): r(r),i(i){}  
    Virt operator + (const Virt &b) {return Virt(r + b.r,i + b.i);}  
    Virt operator - (const Virt &b) {return Virt(r - b.r,i - b.i);}  
    Virt operator * (const Virt &b) {return Virt(r*b.r - i*b.i,r*b.i + i*b.r);}  
    Virt operator * (const DB &t) {return Virt(r*t,i*t);}  
}q[maxn],p[maxn],h[maxn],e1[maxn],e2[maxn],a[maxn];  
   
int n,m,N,L,rev[maxn],dig[maxn];  
DB Now[maxn];  
   
void FFT(Virt *A,int on)  
{  
    for (int i = 0; i < N; i++) a[i] = A[rev[i]];  
    for (int i = 0; i < N; i++) A[i] = a[i];  
    DB t1 = on*2.00;  
    for (int k = 2; k <= N; k <<= 1)  {  
        DB t2 = k;  
        Virt wn = Virt(cos(PI*t1/t2),sin(PI*t1/t2));  
        for (int i = 0; i < N; i += k) {  
            Virt w = Virt(1.00,0.00);  
            for (int j = i; j < i + (k>>1); j++) {  
                Virt u = A[j];  
                Virt t = w*A[j+(k>>1)];  
                A[j] = u + t;  
                A[j+(k>>1)] = u - t;  
                w = w*wn;  
            }  
        }  
    }  
    if (on == -1)  
        for (int i = 0; i < N; i++)  
            A[i].r /= (DB)(N);  
}  
   
int main()  
{  
    #ifdef DMC  
        freopen("DMC.txt","r",stdin);  
    #endif  
       
    cin >> n;  
    for (N = 1; N < n; N <<= 1,++L);  
    N <<= 1; ++L;  
    for (int i = 0; i < N; i++) {  
        for (int x = i,len = 0; x; x >>= 1)  
            dig[len++] = x&1;  
        for (int j = 0; j < L; j++)  
            rev[i] = rev[i]*2 + dig[j];  
    }  
    for (int i = 0; i < n; i++) scanf("%lf",&Now[i]);  
    for (int i = 0; i < N; i++) 
        p[N-i-1] = q[i] = Virt(Now[i],0);    
    for (int i = 1; i < n; i++) 
        h[i] = Virt(1.00/(DB)(i)/(DB)(i),0);
    FFT(h,1); FFT(q,1); FFT(p,1); 
    for (int i = 0; i < N; i++)  
        e1[i] = q[i]*h[i];  
    for (int i = 0; i < N; i++)  
        e2[i] = p[i]*h[i];  
    FFT(e1,-1); FFT(e2,-1);  
    for (int i = 0; i < n; i++)  
        printf("%.8lf\n",e1[i].r - e2[N-i-1].r);  
    return 0;  
}  

  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值