【HDU 1588】 Gauss Fibonacci

最新推荐文章于 2020-07-11 22:22:21 发布

even_bao

最新推荐文章于 2020-07-11 22:22:21 发布

阅读量254

点赞数 1

分类专栏：矩阵乘法

本文链接：https://blog.csdn.net/even_bao/article/details/80456679

版权

矩阵乘法专栏收录该内容

7 篇文章 0 订阅

订阅专栏

【题目链接】

点击打开链接

【算法】

要求

f(g(0)) + f(g(1)) + f(g(2)) + ... + f(g(n-1))

因为g(i) = k * i + b

所以原式 = f(b) + f(k+b) + f(2k+b) + .... + f((n-1)k+b)

令矩阵A = {1,1,0,1}（求斐波那契数的矩阵）

那么，式子就可以写成A^b + A^(k + b) + A ^ (2k + b) + .... + A ^ ((n - 1)k + b)

因为矩阵符合乘法分配律，所以可以将A^b提出，式子被写成 :

A ^ b( E + A ^ k + A ^ 2k + ... + A ^ (n - 1)k ) (其中E为2阶单位阵）

令矩阵S = A ^ k

那么式子就被进一步化简为： A^b( S^0 + S^1 + S^2 + .. + S^(n-1) )

A^b可以通过矩阵乘法快速幂求出

而后面的S^0 + S^1 + S ^ 2 + ... S^(n-1)则可以通过二分求解（也就是POJ 3233的方法）

【代码】

#include<bits/stdc++.h>
using namespace std;

int n,b,k,m;
struct Matrix
{
		long long mat[3][3];		
} A,E,ans,s;

inline Matrix mul(Matrix a,Matrix b)
{
		int i,j,k;
		Matrix ans;
		memset(ans.mat,0,sizeof(ans.mat));
		for (i = 1; i <= 2; i++)
		{
				for (j = 1; j <= 2; j++)
				{
						for (k = 1; k <= 2; k++)
						{
								ans.mat[i][j] = (ans.mat[i][j] + a.mat[i][k] * b.mat[k][j]) % m;
						}
				}
		}
		return ans;
}
inline Matrix add(Matrix a,Matrix b)
{
		int i,j;
		Matrix ans;
		memset(ans.mat,0,sizeof(ans.mat));
		for (i = 1; i <= 2; i++)
		{
				for (j = 1; j <= 2; j++)
				{
						ans.mat[i][j] = (a.mat[i][j] + b.mat[i][j]) % m;
				}
		}
		return ans;
}
inline Matrix power(Matrix a,int n)
{
		int i,j;
		Matrix ans,p = a;
		for (i = 1; i <= 2; i++)
		{
				for (j = 1; j <= 2; j++)
				{
						ans.mat[i][j] = (i == j);	
				}	
		}		
		while (n > 0)
		{
				if (n & 1) ans = mul(ans,p);
				p = mul(p,p);
				n >>= 1;
		}
		return ans;
}
inline Matrix solve(Matrix a,int n)
{
		Matrix tmp;
		if (n == 1) return a;
		if (n % 2 == 1) return add(solve(a,n-1),power(a,n));	
		else
		{
				tmp = solve(a,n/2);
				return add(tmp,mul(power(a,n/2),tmp));
		}
}

int main() {
		
		E.mat[1][1] = 1; E.mat[2][2] = 1;
		E.mat[1][2] = E.mat[2][1] = 0;
		while (scanf("%d%d%d%d",&k,&b,&n,&m) != EOF)
		{
				A.mat[1][1] = A.mat[1][2] = A.mat[2][1] = 1;
				A.mat[2][2] = 0;
				s = power(A,k);
				ans = mul(power(A,b),add(E,solve(s,n-1)));
				printf("%lld\n",ans.mat[1][2]);
		}
		
		return 0;
	
}