POJ3150【FFT】

原创 2015年07月10日 22:41:46

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

HDU 1402 FFT 求 大数乘法

这题的数据量是5w, 也就是传统意义上的n^2算法是不可取的。这里就用到了FFT FFT一般的作用就是使得多项式乘法的复杂度降到nlogn。利用FFT可以快速求出循环卷积。 那么卷积又是什...
  • sdj222555
  • sdj222555
  • 2013年08月06日 10:13
  • 6737

【特征多项式解线性递推】poj2118

叉姐论文:http://www.docin.com/p-724323397.html a[i]=sigma(a[j]*b[k-j]) 求第n项 标准的常系数线性递推,用矩阵乘法可以做到o(k^...
  • huyuncong
  • huyuncong
  • 2014年01月12日 17:06
  • 2342

POJ 难度及题型题目分类

OJ上的一些水题(可用来练手和增加自信) (poj3299,poj2159,poj2739,poj1083,poj2262,poj1503,poj3006,poj2255,poj3094) 初期:...
  • u010871244
  • u010871244
  • 2013年07月16日 10:34
  • 1933

POJ_3150_矩阵快速幂+循环矩阵的乘法

//开始没有主要矩阵实际上是循环的,所以用n^3的算法一直简单优化,比如把中 //最后层,移动到中间 //但,实际上不论E还是自己构造出来的F都是循环的,比如E单位矩阵,实际上就是1的位置发生了移动而...
  • qq_29660153
  • qq_29660153
  • 2017年07月23日 17:42
  • 95

Poj 3150/UVA 1386/UVALive 3704 Cellular Automaton 循环矩阵

#include #include #include using namespace std; #define LL long long LL n,m,d,k,e[505],s[505]; LL...
  • a601025382s
  • a601025382s
  • 2013年08月08日 19:20
  • 1037

poj 3150 Cellular Automaton(矩阵快速幂)

http://poj.org/problem?id=3150 大致题意:给出n个数,问经过K次变换每个位置上的数变为多少。第i位置上的数经过一次变换定义为所有满足 min( abs(i-j),n-...
  • u013081425
  • u013081425
  • 2014年06月12日 21:09
  • 839

POJ 3150 循环矩阵的应用

思路: 首先 先普及一个性质: 循环矩阵*循环矩阵=循环矩阵 由于此题是距离小于d的都加上一个数。 那么 构造矩阵的时候 我们发现 诶呦 这是个循环矩阵 看看数据范围 n...
  • qq_31785871
  • qq_31785871
  • 2016年07月21日 16:36
  • 564

POJ 3150 / Uva 1386 Cellular Automaton 解题报告(循环矩阵)

题目大意:细胞上的数字会随着周围数字的和而变化,求第k轮之后数字的情况。     解题报告:很容易可以想到矩阵乘法。用快速幂的确可以很快的计算出k次方来。     可是n最大有500,矩阵乘法的复杂度...
  • kbdwo
  • kbdwo
  • 2013年08月13日 10:40
  • 865

POJ 3150 Cellular Automaton

矩阵乘法
  • Rlt1296
  • Rlt1296
  • 2017年04月11日 09:16
  • 142

[POJ 3150] Cellular Automaton (矩阵快速幂 + 矩阵乘法优化)

Cellular Automaton Time Limit: 12000MS   Memory Limit: 65536K Total Submissions: 304...
  • SIOFive
  • SIOFive
  • 2014年07月15日 14:14
  • 707
内容举报
返回顶部
收藏助手
不良信息举报
您举报文章:POJ3150【FFT】
举报原因:
原因补充:

(最多只允许输入30个字)