关闭

【Bluestein's Algorithm】[POJ2821]TN's Kingdom III - Assassination

标签: fftcztpoj卷积bluestein
610人阅读 评论(0) 收藏 举报
分类:

题目大意

TN要暗杀Dzx,为了保密,他想到了这样一种方式:首先,把信息编码为N个实数,组成序列α,之后再随便搞一个长度为N的实数序列β。然后按照下面的步骤计算序列γ:
0、做一个空序列γ。
1、把β倒过来。
2、把β向右平移一个元素。最右侧的元素补到左边。
3、计算此时α和β对应元素的积的和。将其加到γ的末尾。
4、如果γ还不足N个元素,重复步骤2和3。

虽然这种加密方法是很弱的,可是那些笨蛋刺客们却没法破解。然后,这些东西就被Dzx拿到了,于是他就躲过了暗杀……现在作为历史学家的你得到了β和γ,你需要解出α来获得研究TN的材料。

分析

说得更直白一点,γ就是αβ的循环卷积。

运用FFT所进行的卷积本身就是循环卷积

我们假设cab两个序列的卷积,即ck=i=0kak×bki
令n为大于等于c的长度而且是2的整数次幂的数,那么卷积也可以写成这个形式ck=i,jakbj[i+j0(modn)]
FFT就是根据第二个式子算出来的,当a,b长度相等而且n等于他们的长度的时候,根据这个式子算出来的就是循环卷积。平时计算的时候由于nc的长度,所以和线性卷积的式子没有区别。
很显然,我们对γβDFT,然后相除,再IDFT即可,但是用FFT进行DFT要求长度必须是2的整数次幂,如果原长不是2的整数次幂,我们强行弄成2的整数次幂的话,在做γDFT和最后IDFT时就会出现问题,因为你根本不知道你这样做完了是一个什么东西。
所以,我们需要一个可以做任意长度卷积的东西,Bluesteins  Algorithm或者混合基FFT

Bluestein’s Algorithm

我们考虑DFT的式子

(11)Ak=j=0N1ajωNjk(12)=j=0N1aje2πiNjk(13)jk=(jk)2+j2+k22(14)Ak=j=0N1ajeπiN((jk)2+j2+k2)(15)=eπiNk2ajeπiNj2eπiN((kj)2)

Xj=ajeπiNj2,Yj=eπiNj2

Ak=eπiNk2j=0NXjYkj

这是一个卷积的形式,我们发现我们可以通过一次卷积来做一次DFT.
ps:其实是CZT,fft是dft的快速算法,其实就是N点dft算法,就是计算量小一点。N点dft的本质是z变换后,在z域单位圆上等间距N点连续采样。
IDFT的时候,可以类比用FFT进行IDFT时做出的改变即可。
我们发现kj可能小于0,我们只需要将Y(x)右移N位即可。

代码

#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
using namespace std;
const double pi=acos(-1);
#define MAXN (1<<17)
int n;
struct cpx{
    double r,i;
    inline cpx(){
    }
    inline cpx(double r,double i=0):r(r),i(i){
    }
    inline cpx operator+(const cpx &b)const{
        return cpx(r+b.r,i+b.i);
    }
    inline cpx operator-(const cpx &b)const{
        return cpx(r-b.r,i-b.i);
    }
    inline cpx operator*(const cpx &b)const{
        return cpx(r*b.r-i*b.i,r*b.i+i*b.r);
    }
    inline cpx operator/(const cpx &b)const{
        return cpx((r*b.r+i*b.i)/(b.r*b.r+b.i*b.i),(i*b.r-r*b.i)/(b.r*b.r+b.i*b.i));
    }
    inline cpx& operator*=(const cpx &b){
        return *this=*this*b;
    }
    inline cpx &operator/=(const cpx &b){
        return *this=*this/b;
    }
}beta [MAXN+10],gamma[MAXN+10],alpha[MAXN+10],A[MAXN*4+10],B[MAXN*4+10];
void fft(cpx *a,int N,int f){
    int i,j,t,k;
    for(i=1,j=0;i<N-1;i++){
        for(t=N;j^=(t>>=1),~j&t;);
        if(i<j)
            swap(a[i],a[j]);
    }
    for(i=1;i<N;i<<=1){
        cpx wn(cos(pi/i),f*sin(pi/i));
        t=i<<1;
        for(j=0;j<N;j+=t){
            cpx w(1,0);
            for(k=0;k<i;k++,w*=wn){
                cpx x(a[j+k]),y(w*a[j+i+k]);
                a[j+k]=x+y;
                a[j+i+k]=x-y;
            }
        }
    }
    if(f==-1){
        for(i=0;i<N;i++)
            a[i]/=N;
    }
}
void bluestein(cpx *a,int n,int f){
    int N,i;
    memset(A,0,sizeof A);
    memset(B,0,sizeof B);
    for(i=0;i<n;i++)
        A[i]=cpx(cos(pi*i*i/n),f*sin(pi*i*i/n))*a[i];
    for(i=0;i<(n<<1);i++)
        B[i]=cpx(cos(pi*(i-n)*(i-n)/n),-f*sin(pi*(i-n)*(i-n)/n));
    for(N=1;N<(n<<2);N<<=1);
    fft(A,N,1);
    fft(B,N,1);
    for(i=0;i<N;i++)
        A[i]*=B[i];
    fft(A,N,-1);
    for(i=0;i<n;i++){
        a[i]=A[i+n]*cpx(cos(pi*i*i/n),f*sin(pi*i*i/n));
        if(f==-1)
            a[i]/=n;
    }
}
void read(){
    scanf("%d",&n);
    int i;
    for(i=0;i<n;i++)
        scanf("%lf",&beta[i].r);
    for(i=0;i<n;i++)
        scanf("%lf",&gamma[i].r);
}
void solve(){
    bluestein(beta,n,1);
    bluestein(gamma,n,1);
    for(int i=0;i<n;i++)
        alpha[i]=gamma[i]/beta[i];
    bluestein(alpha,n,-1);
}
void print(){
    int i;
    for(i=0;i<n;i++)
        printf("%.4f\n",alpha[i].r);
}
int main()
{
    read();
    solve();
    print();
}
1
0
查看评论

[多维FFT Bluestein′s Algorithm] Codechef October Challenge 2017 .Chef and Horcrux

题目里那个 其实不重要,只要能算出 pip_i 就行了pip_i 的话发现就是个多维FFT的转移形式到Hillan大佬博客里拷了个代码改一改就好了…Bluestein算法里的FFT可以用暴力卷积代替,因为数据范围小FFT常数大#include<cstdio> #include<i...
  • Coldef
  • Coldef
  • 2017-10-17 07:47
  • 245

理解快速傅里叶变换(FFT)算法

http://blog.jobbole.com/58246/
  • danteLiujie
  • danteLiujie
  • 2016-06-09 00:03
  • 3641

[FWT][BlueStein Algorithm]CodeChef XORTREEH

SolutionSolution裸的多维FFT,用BlueSteinBlueStein算法。。把里面的复数域的变换,改成膜意义下的变换就好了。 因为这道题维数真的不高。。而且FNT复杂度比较大,所以直接暴力卷积了。。 所以就要打一个FWT一样的东西。 改了Hillan的模板~#include ...
  • Vectorxj
  • Vectorxj
  • 2017-12-21 14:27
  • 91

【bzoj 2821】 作诗 分块

2821: 作诗(Poetize) Time Limit: 50 Sec Memory Limit: 128 MB Submit: 2365 Solved: 670 [Submit][Status][Discuss] Description神犇SJY虐完HEOI之后给傻×LYD...
  • ALPS233
  • ALPS233
  • 2016-03-17 22:04
  • 408

HDU1025:Constructing Roads In JGShining's Kingdom(LIS)

Problem Description JGShining's kingdom consists of 2n(n is no more than 500,000) small cities which are located in two parallel lines. Half o...
  • libin56842
  • libin56842
  • 2013-07-30 17:18
  • 1672

hdu 1025 Constructing Roads In JGShining's Kingdom(LIS)题意较难转换成LIS

1、http://acm.hdu.edu.cn/showproblem.php?pid=1025 题目好难懂,懂了之后也没想到这是一个LIS的题目,注意输出的细节,wrong了n遍 输出,注意一条多条的不同输出,样例之间打印空行 2、题目大意 有2n个城市,其中有n个富有的城市,n个贫穷的城...
  • sdjzping
  • sdjzping
  • 2013-04-05 21:17
  • 2551

【程序18】求s=a+aa+aaa+aaaa+aa...a的值

求s=a+aa+aaa+aaaa+aa…a的值,其中a是一个数字。例如2+22+222+2222+22222(此时共有5个数相加),几个数相加由键盘控制。Tn = 0 Sn = [] n = int(input('n = ')) a = int(input('a = '...
  • fanren224
  • fanren224
  • 2017-03-23 16:54
  • 378

Kadane's algorithm(Kadane算法)

在刷LeetCode的时候遇到了。。查了一下维百,Kadane是卡内基梅隆大学的教授,这个算法是为了解决最大子序列的和(maximum subarray)提出的。 以下资料全部来自维基百科: 1、什么是maximum subarray problem? In computer...
  • lishichengyan
  • lishichengyan
  • 2017-08-13 21:09
  • 1163

杭电OJ——1025 Constructing Roads In JGShining's Kingdom(比较有趣的一道题目,思路详解)

Constructing Roads In JGShining's Kingdom Problem Description JGShining's kingdom consists of 2n(n is no more than 500,000) small cities wh...
  • lishuhuakai
  • lishuhuakai
  • 2012-11-10 01:48
  • 3294

Majority Element问题

问题: 数组中有一个数字出现的次数超过了数组changs
  • a130737
  • a130737
  • 2014-10-04 20:12
  • 2094
    个人资料
    • 访问:99107次
    • 积分:2711
    • 等级:
    • 排名:第15765名
    • 原创:171篇
    • 转载:1篇
    • 译文:0篇
    • 评论:52条
    Contact me
    欢迎各位神犇对我博客中的错误进行指正,或和我进行交流
    最新评论