矩阵相乘，快速算法 HDOJ 4291 A Short problem

对于g（x） % mod；很明显，x一定有一个循环结；

证明：

只要存在g（x）= g（k）， g（x - 1） =  g（k - 1），那么就会出现循环，而g（x）的值是在有限域（1 - mod）内变化，那么g(x) ，g(x - 1)的排列也是有限的，而x是无限的，所以一定会出现g（x）= g（k）， g（x - 1） =  g（k - 1）；

 

那么可以求出每层函数中的mod，依次为

mod1 = 1000000007;

mod2 = 222222224；
mod3 = 183120;
mod4 = 240;

所以i层函数中余modi，中间值余mod（i-1）;

 

而g（n） = 3 * g(n - 1) + g(n - 2);那么有

g[1]                 1                  g[n]                          3          1                            g[n - 1]

             =                                                   =                                      *             

g[0]                 0                  g[n - 1]                     1           0                          g[n - 2]

 

可以得到矩阵的幂，类似与快速幂计算方法，二分计算

 

#include<stdio.h>

const __int64 mod1 = 1000000007;
const __int64 mod3 = 183120;
const __int64 mod2 = 222222224, mod4 = 240;

__int64 result[241];

int copy(__int64 b[][2], __int64 a[][2])
{
	b[0][0] = a[0][0];
	b[0][1] = a[0][1];
	b[1][0] = a[1][0];
	b[1][1] = a[1][1];
	return 0;
}

__int64 f(__int64 a[][2], __int64 x, __int64 b[][2], __int64 mod)
{
	__int64 temp[2][2];
	if(x == 1)
	{
		copy(b, a);
		return 0;
	}
	f(a, x / 2, temp, mod);
	b[0][0] = (temp[0][0] * temp[0][0] + temp[0][1] * temp[1][0]) % mod;
	b[0][1] = (temp[0][0] * temp[0][1] + temp[0][1] * temp[1][1]) % mod;
	b[1][0] = (temp[1][0] * temp[0][0] + temp[1][1] * temp[1][0]) % mod;
	b[1][1] = (temp[1][0] * temp[0][1] + temp[1][1] * temp[1][1]) % mod;
	if(x % 2 == 0)
		return 0;
	copy(temp, b);
	b[0][0] = (temp[0][0] * a[0][0] + temp[0][1] * a[1][0]) % mod;
	b[0][1] = (temp[0][0] * a[0][1] + temp[0][1] * a[1][1]) % mod;
	b[1][0] = (temp[1][0] * a[0][0] + temp[1][1] * a[1][0]) % mod;
	b[1][1] = (temp[1][0] * a[0][1] + temp[1][1] * a[1][1]) % mod;
	return 0;
}

__int64 g(__int64 x, __int64 mod)
{
	__int64 a[2][2], b[2][2];
	if(x < 2)
		return x;
	a[0][0] = 3;
	a[0][1] = 1;
	a[1][0] = 1;
	a[1][1] = 0;
	f(a, x - 1, b, mod);
	return b[0][0];
}

int init()
{
	int i;
	result[0] = 0;
	result[1] = 1;
	for(i = 2; i <= 240; i ++)
	{
		result[i] = g(g(g(i % mod4, mod3) % mod3, mod2) % mod2, mod1) % mod1;
	}
	return 0;
}

int main()
{
	__int64 n;
	init();
	while(scanf("%I64d", &n) != EOF)
	{
		printf("%I64d\n", result[n % mod4]);
	}
	return 0;
}