POJ3150【FFT】

转移很好用矩阵表示.然而矩阵乘法复杂度是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;
}

版权声明：本文为博主原创文章，未经博主允许不得转载。