矩阵相乘，快速算法 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
                          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;
}