hdu 4291 A Short problem 矩阵快速幂,找循环节

A Short problem

Time Limit: 2000/1000 MS (Java/Others)    Memory Limit: 32768/32768 K (Java/Others)
Total Submission(s): 1560    Accepted Submission(s): 568

Problem Description

　　According to a research, VIM users tend to have shorter fingers, compared with Emacs users.
　　Hence they prefer problems short, too. Here is a short one:
　　Given n (1 <= n <= 10 ¹⁸), You should solve for 
g(g(g(n))) mod 10 ⁹ + 7
　　where
g(n) = 3g(n - 1) + g(n - 2)
g(1) = 1
g(0) = 0

 

Input

　　There are several test cases. For each test case there is an integer n in a single line.
　　Please process until EOF (End Of File).

 

Output

　　For each test case, please print a single line with a integer, the corresponding answer to this case.

 

Sample Input

  
  

   
   0
1
2
  
  

 

Sample Output

  
  

   
   0
1
42837
  
  

这题和以往的矩阵快速幂的题目不同在于g(g(g(n))) mod 10⁹ + 7;   这个取模 是把最外层的g( ) 函数的 因变量改了. g()函数是比斐波那契函数的增长速度还要快的. 所以最外层g( ) 的自变量 g(g(n)) 是非常大的.  所以考虑找循环节.

最外层的g()的因变量是要%1000000007的,而且g(0)=0,g(1)=1,所以可以通过下面的程序找到 g(g(n))的循环节是222222224.

所以第二层的g()的因变量是要%222222224;而且它同样有g(0)=0,g(1)=1;所以可以把下面的程序的①号处的mod 改为mod1; 

从而能算出g(n)的循环节是183120;

所以n带入后 用矩阵快速幂 算出的答案 首先要mod183120; 然后把结果再带入g(), 算出 g(g(n))mod222222224 的值;  在带入g(), 算出g(g(g(n)))mod1000000007的值.

#include <stdio.h>
#define mod 1000000007
#define mod1 222222224
#define mod2 183120
int main()
{
__int64 a=0,b=1,c,i;
for(i=2;;i++){//求循环节 
c=(b*3+a)%mod;//①
a=b;
b=c;
if(a == 0 && b == 1)
{
printf("%I64d\n",i-1);
break;
}
}
return 0;
}

#include<stdio.h>
#include<string.h>
#define Matr 3 //矩阵大小,注意能小就小   矩阵从1开始   所以Matr 要+1 
#define ll __int64
struct mat//矩阵结构体，a表示内容，size大小 矩阵从1开始   但size不用加一
{
    ll a[Matr][Matr];
    mat()//构造函数
    {
        memset(a,0,sizeof(a));
    }
};
int Size=2;
ll mod;

mat multi(mat m1,mat m2)//两个相等矩阵的乘法，对于稀疏矩阵，有0处不用运算的优化 
{
    mat ans=mat(); 
    for(int i=1;i<=Size;i++)
        for(int j=1;j<=Size;j++)
            if(m1.a[i][j])//稀疏矩阵优化 
                for(int k=1;k<=Size;k++)
                    ans.a[i][k]=(ans.a[i][k]+m1.a[i][j]*m2.a[j][k])%mod; //i行k列第j项
    return ans;
}

mat quickmulti(mat m,ll n)//二分快速幂 
{
    mat ans=mat();
    int i;
    for(i=1;i<=Size;i++)ans.a[i][i]=1;
    while(n)
    {
        if(n&1) ans=multi(m,ans);//奇乘偶子乘 挺好记的.
        m=multi(m,m);
        n>>=1;
    }
    return ans;
}
mat gouzao=mat(),chu=mat();
ll solve(ll n)
{
	if(n==0)
	{
		return 0;
	}
	if(n==1)
	{
		return 1;
	}
	mat ans;
	ans=multi(chu,quickmulti(gouzao,n-1));
	return ans.a[1][2]%=mod;
}
int main()
{
	ll n;
	ll tem;
	ll kmod[3]={183120,222222224,1000000007};
	chu.a[1][1]=0;
	chu.a[1][2]=1;
	gouzao.a[1][2]=gouzao.a[2][1]=1;
	gouzao.a[2][2]=3; 
	while(~scanf("%I64d",&n))
	{
		tem=n;
		for(int i=0;i<3;i++)
		{
			mod=kmod[i];
			if(tem>=2)
			tem=solve(tem);
		}
		printf("%I64d\n",tem);
	}
	return 0;
}