您的位置:首页 > 其它

POJ3150【FFT】

2015-07-10 22:41 405 查看
转移很好用矩阵表示.然而矩阵乘法复杂度是O(n^3)的.

很容易发现转移矩阵是【循环矩阵】.而且有一个美妙的性质:【循环矩阵 * 循环矩阵 = 循环矩阵】.

所以我们计算矩阵乘法的时候可以只计算第一行.剩下的可以由第一行递推得出.

一次乘法的复杂度降到了O(n^2).这是可以接受的.

#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <ctime>
#include <algorithm>
#include <iostream>
#include <fstream>
#include <vector>
#include <queue>
#include <deque>
#include <map>
#include <set>
#include <string>
#define make make_pair
#define fi first
#define se second

using namespace std;

typedef long long ll;
typedef unsigned long long ull;
typedef pair<int,int> pii;

const int maxn = 1000010;
const int maxm = 1010;
const int maxs = 30;
const int inf = 0x3f3f3f3f;
const int P = 1000000007;
const double error = 1e-9;
const double Pi = 3.141592653589793238462643383279;

inline int read()
{
int x = 0, f = 1; char ch = getchar();
while (ch <= 47 || ch >= 58)
f = (ch == 45 ? -1 : 1), ch = getchar();
while (ch >= 48 && ch <= 57)
x = x * 10 + ch - 48, ch = getchar();
return x * f;
}

struct matrix
{
ll num[501][501];
} f, e;

int n, m, d, k; ll a[501], b[501];

matrix operator * (matrix a, matrix b)
{
matrix c;
for (int i = 1; i <= n; i++)
for (int j = 1; j <= n; j++)
c.num[i][j] = 0;
for (int i = 1; i <= n; i++)
for (int j = 1; j <= n; j++)
(c.num[1][i] += a.num[1][j] * b.num[j][i]) %= m;
for (int i = 2; i <= n; i++)
for (int j = 1; j <= n; j++)
c.num[i][j] = c.num[i - 1][j == 1 ? n : j - 1];
return c;
}

void init()
{
for (int i = 1; i <= n; i++)
e.num[i][i] = 1, a[i] = read();
for (int i = 1; i <= n; i++)
for (int j = 1; j <= n; j++) {
int a = min(i, j), b = max(i, j);
int dis = min(b - a, n + a - b);
if (dis <= d) f.num[i][j] = 1;
}
}

int main()
{
n = read(), m = read();
d = read(), k = read();

init();

for (int i = k; i != 0; f = f * f, i /= 2)
if (i & 1) e = e * f;

for (int i = 1; i <= n; i++)
for (int j = 1; j <= n; j++)
(b[i] += a[j] * e.num[j][i]) %= m;

for (int i = 1; i <= n; i++)
printf("%lld%c", b[i], i == n ? 10 : 32);

return 0;
我们继续考虑.既然循环矩阵可以由第一行递推得出.那么我们可以只保存矩阵的第一行和第一列.

计算乘法的时候.发现很容易变成卷积形式.那么就可以用FFT优化了.

在本题中.转移阵不仅是【循环矩阵】.而且关于主对角线对称.因此第一行和第一列相同.只用保存一个.

这道题用FFT可以做到O(n*logn*logk).在POJ上排名rank1.Excited!



/* I will wait for you */

#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <ctime>
#include <algorithm>
#include <iostream>
#include <fstream>
#include <vector>
#include <queue>
#include <deque>
#include <set>
#include <map>
#include <string>
#define make make_pair
#define fi first
#define se second

using namespace std;

typedef long long ll;
typedef unsigned long long ull;
typedef pair<int,int> pii;

const int maxn = 1000010;
const int maxm = 501;
const int maxs = 30;
const int inf = 0x3f3f3f3f;
const int P = 1000000007;
const double error = 1e-9;
const double Pi = 3.141592653589793238462643383279;

inline int read()
{
int x = 0, f = 1; char ch = getchar();
while (ch <= 47 || ch >= 58)
f = (ch == 45 ? -1 : 1), ch = getchar();
while (ch >= 48 && ch <= 57)
x = x * 10 + ch - 48, ch = getchar();
return x * f;
}

struct complex
{
double re, im;
} _x[maxn], _y[maxn], w[2][maxn];

complex operator + (complex a, complex b)
{
complex c;
c.re = a.re + b.re;
c.im = a.im + b.im;
return c;
}

complex operator - (complex a, complex b)
{
complex c;
c.re = a.re - b.re;
c.im = a.im - b.im;
return c;
}

complex operator * (complex a, complex b)
{
complex c;
c.re = a.re * b.re - a.im * b.im;
c.im = a.im * b.re + a.re * b.im;
return c;
}

struct matrix
{
ll num[maxm];
} a, f;

int n, m, d, k, _n = 1, rev[maxn];

void _init()
{
for (int i = 0; i < _n; i++)
for (int j = i, k = 1; k < _n; k <<= 1, j >>= 1)
(rev[i] <<= 1) |= (j & 1);
for (int i = 0; i < _n; i++) {
w[0][i].re = cos(2 * Pi * i / _n);
w[0][i].im = sin(2 * Pi * i / _n);
w[1][i].re = cos(2 * Pi * i / _n);
w[1][i].im = -sin(2 * Pi * i / _n);
}
}

void FFT(complex *a, int f)
{
for (int i = 0; i < _n; i++)
if (rev[i] > i) swap(a[i], a[rev[i]]);
for (int i = 1; i < _n; i <<= 1)
for (int j = 0, l = _n / (i << 1); j < _n; j += (i << 1))
for (int k = 0, t = 0; k < i; k += 1, t += l) {
complex x = a[j + k], y = w[f][t] * a[i + j + k];
a[j + k] = x + y, a[i + j + k] = x - y;
}
for (int i = 0; f && i < _n; i++) a[i].re /= _n;
}

matrix operator * (matrix a, matrix b)
{
matrix c;
for (int i = 0; i < _n; i++) {
_x[i].re = _x[i].im = 0;
_y[i].re = _y[i].im = 0;
}
for (int i = 0; i < n; i++) {
_x[i].re = a.num[n - i];
_y[i].re = b.num[i + 1];
_y[n + i].re = _y[i].re;
}
FFT(_x, 0), FFT(_y, 0);
for (int i = 0; i < _n; i++)
_x[i] = _x[i] * _y[i];
FFT(_x, 1);
for (int i = 1; i <= n; i++)
c.num[i] = (ll) (_x[2 * n - i].re + 0.5) % m;
return c;
}

void init()
{
for (int i = 1; i <= n; i++)
a.num[i] = read();
for (int i = 1; i <= n; i++) {
int dis = min(i - 1, n + 1 - i);
f.num[i] = (dis <= d);
}
while (_n < n << 2) _n <<= 1;
}

int main()
{
n = read(), m = read();
d = read(), k = read();
init(), _init();

for (int i = k; i != 0; i /= 2, f = f * f)
if (i & 1) a = a * f;
for (int i = 1; i <= n; i++)
printf("%lld ", (ll) (a.num[i] + 0.5));

return 0;
}
内容来自用户分享和网络整理,不保证内容的准确性,如有侵权内容,可联系管理员处理 点击这里给我发消息
标签: