2018 Multi-University Training Contest 7--HDU 6395 Sequence（矩阵快速幂+分块）

题意：

给出一个公式，让你求出第n个数值是多少。

题解：

看到公式第一个就想起了矩阵乘法，因为很像斐波那契数列（可以用矩阵乘法得到），但是后面有个尾巴，所以要把尾巴考虑进去。

推出公式如下：

这样我们只考虑左矩阵的次方，再乘上右矩阵就是我们要找到答案，可以用矩阵快速幂来缩短这段时间。

struct matrix{
	ll a[3][3];
	
	matrix(){
		mem(a);
	}
	matrix operator* (matrix &temp){
		matrix ans;
		for(int i = 0; i < 3; i++){
			for(int j = 0; j < 3; j++){
				for(int k = 0; k < 3; k++){
					ans.a[i][j] = ans.a[i][j] + temp.a[k][j]*a[i][k];
				}
				ans.a[i][j] %= mod;
			}
		}
		return ans;
	}
}m1;

matrix qpow(matrix a, int b){
	matrix sum = m1;
	while(b){
		if(b&1){
			sum = sum*a;
		}
		a = a*a;
		b >>= 1;
	}
	return sum;
}

但是 p/n 是一个不定因素，所以不能直接就用矩阵快速幂。

假设 p = 40, 我们要求的是 n = 40:

p/3 = 13,     p/4 = 10,      p/5 = 8,     p/6 = 6,     p/8 = 5,     p/10 = 4…

这样我们可以发现[3，3]的时候 p/n 都为13，[4，4]的时候 p/n 都为10，[5，5]的时候 p/n 都为8，[6，6]的时候 p/n 都为6，[7，8]的时候 p/n 都为5，[9，10]的时候 p/n 都为4...

不难看出，我们可以对这些区间进行矩阵快速幂，这样就快了很多。

比赛的时候，分区分错了（从 p/3 =  13开始递减到 0，没有考虑到会跳的，太菜了，下次一定要数据开大点看！！！），导致tle，石乐志233333。

#include <algorithm>
#include  <iostream>
#include   <cstdlib>
#include   <cstring>
#include    <cstdio>
#include    <string>
#include    <vector>
#include    <bitset>
#include     <stack>
#include     <cmath>
#include     <deque>
#include     <queue>
#include      <list>
#include       <set>
#include       <map>
#define mem(a) memset(a, 0, sizeof(a))
#define pi acos(-1)
using namespace std;
typedef long long ll;

const int mod = 1e9+7;
struct matrix{
	ll a[3][3];
	
	matrix(){
		mem(a);
	}
	matrix operator* (matrix &temp){
		matrix ans;
		for(int i = 0; i < 3; i++){
			for(int j = 0; j < 3; j++){
				for(int k = 0; k < 3; k++){
					ans.a[i][j] = ans.a[i][j] + temp.a[k][j]*a[i][k];
				}
				ans.a[i][j] %= mod;
			}
		}
		return ans;
	}
}m1;

matrix qpow(matrix a, int b){
	matrix sum = m1;
	while(b){
		if(b&1){
			sum = sum*a;
		}
		a = a*a;
		b >>= 1;
	}
	return sum;
}

int main(){
	m1.a[0][0] = m1.a[1][1] = m1.a[2][2] = 1;
	int t;
	scanf("%d", &t);
	while(t--){
		int a, b, c, d, p, n;
		scanf("%d %d %d %d %d %d", &a, &b, &c, &d, &p, &n);
		if(n == 2){
			printf("%d\n", b);
			continue;
		}
		if(n == 1){
			printf("%d\n", a);
			continue;
		}
		ll ans[3] = {b, a, 0};
		matrix temp;
		temp.a[0][0] = d;
		temp.a[0][1] = c;
		temp.a[0][2] = 1;
		temp.a[1][0] = 1;
		temp.a[2][2] = 1;
		int s = 3;
		for(int kl = p/3.0; kl >= 0;){
			int l = s,r;
			if(kl == 0){
				r = n;
			}
			else{
				r = min(p/kl, n);
			}
			cout<<kl<<" "<<l<<" "<<r<<endl;
			int len = r-l+1;
			ans[2] = kl;
			s = r+1;
			matrix qmatrix = qpow(temp, len);
			ll ans1, ans2;
			ans1 = (qmatrix.a[0][0]*ans[0] + qmatrix.a[0][1]*ans[1] + qmatrix.a[0][2]*ans[2]) % mod;
			ans2 = (qmatrix.a[1][0]*ans[0] + qmatrix.a[1][1]*ans[1] + qmatrix.a[1][2]*ans[2]) % mod;
			ans[0] = ans1;
			ans[1] = ans2;
			if(r == n){
				break;
			}
			kl = p/s;
		}
		printf("%d\n", ans[0]);
	}
}
//1000000000 1000000000 1000000000 1000000000 1000000000 1000000000