FOJ 1692 Key problem

       更新：2014-02-07    清晨

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

       来源：http://acm.fzu.edu.cn/problem.php?pid=1692

       概述：有个大小为n的初值"向量"A[i]，顺时针围成一个圈，然后进行m轮操作。每轮每个A[i]进行操作：A[i] = A[i] + L*A[(i-1+n)%n] + R*A[(i+1)%n]。最后的A[i]对M取余。

       勘误：原题L和R上下文没有匹配，一种修正方法是：Sample Input的输入顺序是"n m R L M"，其它不变。

       补充：

       （1）A[i] = A[i] + L*A[(i-1+n)%n] + R*A[(i+1)%n]解释：这是围成圈的数值，常用的索引方法，要非常熟悉。如果元素按0~n-1编号，顺时针围成一个圈，那么A[(i-1+n)%n] 就是A[i]左边的值，A[(i+1)%n]就是A[i]右边的值。

       （2）为了更好的理解题意，我制作了图1，粗斜体数代表输入值。

       （3）建议先了解一道更简单的题：http://blog.csdn.net/code4101/article/details/18893319

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

理论分析

        

       

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

编程分析

       （1）题目最大对1e6取余，所以元素乘积最大为1e12，会超出int的范围，保险起见，我们对所有的向量、矩阵都用__int64类型。

       （2）最后的Y向量，可以不用数组完成。

       （3）POW函数中，我的快速幂模板，变量名都是固定有规律的，T数组作为中间存储变量，P存储幂运算前面的系数（在数值里P初始为1，在矩阵里P初始为单位阵），A存储幂运算的底数，b存储幂运算的指数。

       （4）效率优化最主要的是“理论分析”中的结论。其次，因为100个乘积累加和最大为1e14，所以可以求和后，再一次性取MOD，而不是每次运算都取MOD。实验证明，这样效率能提高47ms。

#include <cstdio>
#include <cstring>
using namespace std;

typedef __int64 i64;

#define N 100
const int SUM = sizeof(i64)*N*N;
int n, m, L, R, MOD;

// 根据该题特殊矩阵制定的算法
void MUL(i64 a[][N], i64 b[][N], i64 c[][N])
{
	int i, j, k;

	// 只算第一行 
	for (j = 0; j < n; j++)
	{
		for (c[0][j] = k = 0; k < n; k++)
			c[0][j] += a[0][k]*b[k][j];
		c[0][j] %= MOD;
	}
	
	// 利用性质,完成接下来几行 
	for (i = 1; i < n; i++)
	{
		c[i][0] = c[i-1][n-1];
		for (j = 1; j < n; j++)
			c[i][j] = c[i-1][j-1];
	}
}

void POW(int b, i64 P[][N])
{
	int i, j, k;
	i64 A[N][N] = {0}, T[N][N];
	
	// init P为单位阵 
	memset(P, 0, SUM); for (i = 0; i < n; i++) P[i][i] = 1;
	
	// init A
	A[0][0] = 1;   A[0][1] = R;     A[0][n-1] = L;
	A[n-1][0] = R; A[n-1][n-2] = L; A[n-1][n-1] = 1;
	for (i = 1; i < n - 1; i++){A[i][i]	= 1; A[i][i-1] = L; A[i][i+1] = R;}
	
	// 快速幂模 
	while (b)
    {  
        if (b&1)
        {
            MUL(P, A, T);
            memcpy(P, T, SUM);
        }
        b >>= 1;
        MUL(A, A, T);
        memcpy(A, T, SUM);
	}
}

int main()
{
	int i, k, T;
	i64 P[N][N], X[N], Y;
	scanf("%d", &T);
	while (T--)
	{
		scanf("%d%d%d%d%d", &n, &m, &R, &L, &MOD);
		for (i = 0; i < n; i++) scanf("%I64d", &X[i]);
		
		POW(m, P);
		
		for (i = 0; i < n - 1; i++)
		{
			for (Y = k = 0; k < n; k++) Y += P[i][k]*X[k];
			printf("%I64d ", Y%MOD);
		}
		for (Y = k = 0; k < n; k++) Y += P[i][k]*X[k];
		printf("%I64d\n", Y%MOD);
	}
	return 0;
}

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

补充：memset与for的效率比较

       经过多次代码测试证明：在该题，使用memset、memcpy对N*N的整个数组进行操作，与用for进行n*n的局部数组操作，效率不变，都是187ms。

       另外，把变量L,R等和P[N][N],T[N][N]设置为局部变量或全局变量，效率也没有差异，所以以后编程怎么方便就怎么写，不要纠结这种小效率的优化，遇到特别大的数组时，再考虑是否设为全局变量。

       还有，VC++效率也比G++快，用G++慢了31ms。

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

方法进阶

       通过查阅资料，了解到这种矩阵称为循环矩阵。

       请先阅读该篇博客：http://blog.csdn.net/code4101/article/details/18956275，文章对本文的理论分析最后部分给出了证明。

       读完该篇博客，相信读者对循环矩阵也有了更进一步的了解，从C=circ(c[0],c[1],...c[n-1])的形式发现，循环矩阵只需用一维数组进行存储就行了。所以上述代码还能进一步优化，把循环矩阵的二维数据结构改成一维，相应对MUL函数等进行更新。

       我的代码中，用i,j代表行列，用k代表下标。

#include <cstdio>
#include <cstring>
using namespace std;

typedef __int64 i64;

#define N 100
int n, m, L, R, MOD;

// circ matrix multiplicative
void MUL(i64 a[N], i64 b[N], i64 c[N])
{
	int i, j, k;
	for (i = 0; i < n; i++)
	{
		for (c[i] = j = 0, k = i + n; j < n; j++, k--)
			c[i] += a[j]*b[k%n];
		c[i] %= MOD;
	}
}

void POW(int b, i64 P[N])
{
	int i, j, k, SUM = sizeof(i64)*n;
	i64 A[N] = {0}, T[N];
	
	// init P为单位阵 ,单位阵也是特殊的循环矩阵 
	memset(P, 0, SUM); P[0] = 1;
	
	// init A
	A[0] = 1; A[1] = R; A[n-1] = L;
	
	// 快速幂模 
	while (b)
    {  
        if (b&1)
        {
            MUL(P, A, T);
            memcpy(P, T, SUM);
        }
        b >>= 1;
        MUL(A, A, T);
        memcpy(A, T, SUM);
	}
}

int main()
{
	int i, j, k, T;
	i64 P[N], X[N], Y;
	scanf("%d", &T);
	while (T--)
	{
		scanf("%d%d%d%d%d", &n, &m, &R, &L, &MOD);
		for (i = 0; i < n; i++) scanf("%I64d", &X[i]);
		
		POW(m, P);
		
		for (i = 0; i < n - 1; i++)
		{
			for (Y = j = 0, k = n - i; j < n; k++, j++) Y += P[k%n]*X[j];
			printf("%I64d ", Y%MOD);
		}
		for (Y = j = 0, k = n - i; j < n; k++, j++) Y += P[k%n]*X[j];
		printf("%I64d\n", Y%MOD);
	}
	return 0;
}

       进阶方法效率上又更进一步，用时109ms。

       花絮：奇怪的是，换成SUM = sizeof(i64)*N;竟然用时125ms。