poj2821 循环卷积 bluestein算法

(6 mins to read)

已知多项式A*B=C,其中*为循环卷积,给出B和C求A对B和C进行任意长度的DFT(bluestein,狭义的CZT),然后点值相除再变换回去即可$ij = \binom{i+j}{2} - \binom{i}{2} - \binom{j}{2}$

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
#include <cstdio>
#include <cmath>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;
const int N = (1<<19) + 5;
#define db double
struct cp
{
db x, y;
cp() {}
cp(db _x, db _y): x(_x), y(_y) {}
cp operator + (cp &oth) { return cp(x+oth.x, y+oth.y); }
cp operator - (cp &oth) { return cp(x-oth.x, y-oth.y); }
cp operator * (cp &oth) { return cp(x*oth.x-y*oth.y, x*oth.y+y*oth.x); }
cp operator / (cp &oth) { return cp((x*oth.x+y*oth.y)/(oth.x*oth.x+oth.y*oth.y), (y*oth.x-x*oth.y)/(oth.x*oth.x+oth.y*oth.y)); }
};
int up, rev[N];
cp w[N], pw[N];
const db pi = acos(-1.0);
void init(int n)
{
up = 1; int l = 0;
while(up<=n) up <<= 1, l++;
for(int i=0; i<up; i++) rev[i] = (rev[i>>1]>>1)|((i&1)<<(l-1));
for(int i=1; i<up; i<<=1)
for(int j=0; j<i; j++)
w[i+j] = cp(cos(pi*j/i), sin(pi*j/i));
}
void DFT(cp *a, int l)
{
static cp t;
for(int i=0; i<l; i++) if(i<rev[i]) swap(a[i], a[rev[i]]);
for(int i=1; i^l; i<<=1)
for(int j=0, d=i<<1; j^l; j+=d)
for(int k=0; k<i; k++)
t = a[i|j|k]*w[i|k], a[i|j|k] = a[j|k]-t, a[j|k] = a[j|k]+t;
}
void IDFT(cp *a, int l)
{
reverse(a+1, a+l); DFT(a, l);
for(int i=0; i<l; i++) a[i].x /= l, a[i].y /= l;
}
int n;
cp a[N], b[N], c[N];
int comb(int x)
{
return 1ll*x*(x-1)/2%n;
}
void bluestein(cp *a, int n)
{
static cp x[N], y[N];
for(int i=0; i<up; i++) x[i].x = x[i].y = y[i].x = y[i].y = 0.0;
for(int i=0; i<n; i++) x[i] = pw[n-comb(i)]*a[i];
for(int i=0; i<2*n-1; i++) y[i] = pw[comb(i)];
reverse(x, x+n);
DFT(x, up); DFT(y, up);
for(int i=0; i<up; i++) x[i] = x[i]*y[i];
IDFT(x, up);
for(int i=n-1; i<=2*n-2; i++) a[i-n+1] = x[i]*pw[n-comb(i-n+1)];
}
void Ibluestein(cp *a, int n)
{
reverse(a+1, a+n); bluestein(a, n);
for(int i=0; i<n; i++) a[i].x /= n, a[i].y /= n;
}
void div(cp *a, cp * b, cp * c, int n)
{
bluestein(a, n);
bluestein(b, n);
for(int i=0; i<n; i++) c[i] = a[i]/b[i];
Ibluestein(c, n);
}
int main()
{
scanf("%d", &n);
for(int i=0; i<n; i++) scanf("%lf", &b[i].x);
for(int i=0; i<n; i++) scanf("%lf", &c[i].x);
init(3*n-3);
pw[0] = cp(1, 0); pw[1] = cp(cos(2*pi/n), sin(2*pi/n));
for(int i=2; i<=n; i++) pw[i] = pw[i-1]*pw[1];
div(c, b, a, n);
for(int i=0; i<n; i++) printf("%.4f\n", a[i].x);
return 0;
}