Hdu6395 Sequence （分块矩阵快速幂）

Sequence

Time Limit: 4000/2000 MS (Java/Others)    Memory Limit: 262144/262144 K (Java/Others)
Total Submission(s): 1162    Accepted Submission(s): 426

 

Problem Description

Let us define a sequence as below

  Your job is simple, for each task, you should output Fn module 109+7.

 

Input

The first line has only one integer T, indicates the number of tasks.

Then, for the next T lines, each line consists of 6 integers, A , B, C, D, P, n.

1≤T≤200≤A,B,C,D≤1091≤P,n≤109

 

Sample Input

2

3 3 2 1 3 5

3 2 2 2 1 4

 

Sample Output

36

24

 

Source

2018 Multi-University Training Contest 7

 题目大意：

给出Fn的公式，给出n，求Fn。

解题思路：

看上去很像一个类似斐波那契数列的公式，当然就想到矩阵快速幂了，但是公式中又加了p/n，因为n会变，所以导致每个fn都在变，所以不能直接用矩阵快速幂，当然p/n也可能有多个n使得p/n是相同的，所以我们可以将其分块的进行矩阵快速幂，举个例子：p=8，n=6的时候，当i=3时，p/i=2，而i=4的时候，p/i也等于2，所以3，和4是一块的，所以我们进行矩阵快速幂的时候，直接将工具矩阵自乘两次，times=j-p/(p/i)+1,但是也要注意n的限制，当n比p/(p/i)小的时候就要取n

long long j=min(n,P/(P/i)); //进行分块 
matrix out=ans;
out.a[0][2]=P/i;
out=pow(out,j-i+1);

接下来最重要的是怎么建立工具矩阵：我们先来观察两个Fn，假设在i=3和4时P/n相等，设p=P/n，F3=C*A+D*B+p，F4=B*C+D*F3+p=B*C+C*D*A+D*D*B+D*p+p=(C+D*D)*B+(C*D)*A+(D+1)*p

则有F3=C*A+D*B+p   F4=(C+D*D)*B+(C*D)*A+(D+1)*p 那我们就要想办法将F3和F4的(A和B)的系数出现在同一个矩阵，首先我们先把A和B的初始系数确定好，设矩阵为a[3][3],a[0][0]=D(B的系数),a[0][1]=C（A的系数），因为要将D变成D*D+C，那么就要自乘一个D，再加上一个C，又因为B的系数是在第一行第一列，所以对第一行*第一列要得到D*D+C，那么a[1][0]=1,a[2][0]=0,而对于A的系数，则是第一行乘以第二列，要将C变成C*D，那么a[1][1]=a[2][1]=0,而对于p我们不妨设a[0][2]=p,而将P变成p*（D+1），那就是将第1行乘以第三列，那么a[1][2]=0,a[2][2]=1，到此矩阵就建立好了，而且不难发现在自乘之后，该矩阵的第二行就等于原矩阵的第一行，所以我们就能保证Fn-1和Fn-2的系数在同一个矩阵了。

ans.a[0][0]=D; //建立工具矩阵 
ans.a[0][1]=C;
ans.a[1][0]=1;
ans.a[2][2]=1;
long long a=(out.a[1][0]*B+out.a[1][1]*A+out.a[1][2])%mod;
long long b=(out.a[0][0]*B+out.a[0][1]*A+out.a[0][2])%mod;
A=a,B=b;

代码：

#include<bits/stdc++.h>
using namespace std;
const long long mod =1e9+7;
struct matrix{ //设立矩阵 
	long long a[3][3];
	matrix(){
		memset(a,0,sizeof(a));
	}
};
matrix mul_matrix(matrix a,matrix b){ //矩阵的相乘 
	matrix c;
	for(int i=0;i<3;i++){
		for(int j=0;j<3;j++){
			for(int k=0;k<3;k++){
				c.a[i][j]+=a.a[i][k]*b.a[k][j];
				c.a[i][j]%=mod;
			}
		}
	}
	return c;
}
matrix pow(matrix x,long long a){ //矩阵快速幂 
	matrix ans;
	ans.a[0][0]=ans.a[1][1]=ans.a[2][2]=1; //初始化为单位矩阵 
	while(a){
		if(a&1){
			ans=mul_matrix(ans,x);
		}
		x=mul_matrix(x,x);
		a/=2;
	}
	return ans;
}
int main()
{
	int t;
	scanf("%d",&t);
	while(t--){
		long long A,B,C,D,P,n;
		scanf("%lld %lld %lld %lld %lld %lld",&A,&B,&C,&D,&P,&n);
		if(n==1){ //如果是1和2就直接输出 
			printf("%lld\n",A);
			continue;
		}
		else if(n==2){
			printf("%lld\n",B);
			continue;
		}
		matrix ans;
		ans.a[0][0]=D; //建立工具矩阵 
		ans.a[0][1]=C;
		ans.a[1][0]=1;
		ans.a[2][2]=1;
		for(int i=3;i<=n;){
			if(P/i==0){ //当P小于i的时候直接可以进行矩阵快速幂 
				matrix out=ans;
				out=pow(out,n-i+1);
				printf("%lld\n",(A*out.a[0][1]+B*out.a[0][0]+out.a[0][2])%mod);
				goto here;
			}
			long long j=min(n,P/(P/i)); //进行分块 
			matrix out=ans;
			out.a[0][2]=P/i;
			out=pow(out,j-i+1);
			long long a=(out.a[1][0]*B+out.a[1][1]*A+out.a[1][2])%mod;
			long long b=(out.a[0][0]*B+out.a[0][1]*A+out.a[0][2])%mod;
			A=a,B=b;
			i=j+1;
		}
		printf("%lld\n",B);
		here:;
	}
}