[矩阵快速幂][等比矩阵求和]Gauss Fibonacci HDU1588

Without expecting, Angel replied quickly.She says: "I'v heard that you'r a very clever boy. So if you wanna me be your GF, you should solve the problem called GF~. "
How good an opportunity that Gardon can not give up! The "Problem GF" told by Angel is actually "Gauss Fibonacci".
As we know ,Gauss is the famous mathematician who worked out the sum from 1 to 100 very quickly, and Fibonacci is the crazy man who invented some numbers.

Arithmetic progression:
g(i)=k*i+b;
We assume k and b are both non-nagetive integers.

Fibonacci Numbers:
f(0)=0
f(1)=1
f(n)=f(n-1)+f(n-2) (n>=2)

The Gauss Fibonacci problem is described as follows:
Given k,b,n ,calculate the sum of every f(g(i)) for 0<=i<n
The answer may be very large, so you should divide this answer by M and just output the remainder instead.

Input

The input contains serveral lines. For each line there are four non-nagetive integers: k,b,n,M
Each of them will not exceed 1,000,000,000.

Output

For each line input, out the value described above.

Sample

Input	Output
2 1 4 100 2 0 4 100	21 12

题意： 给出b，k，n，mod，设f(i)表示Fibonacci第i项，求f(b)+f(b+k)+......+f(b+(n-1)*k)的值余mod后的结果。

分析： 首先先利用矩阵化简一下表达式，设矩阵a为[1, 1; 1, 0]，也就是求Fibonacci第n项用到的那个矩阵，根据之前的知识，我们可以将f(i)转换为a^(i-1)*[f(1); f(0)]，设矩阵[f(1); f(0)]为m，同样的道理我们可以将表达式转化为a^(b-1)*m+a^(b+k-1)*m+......+a^(b+(n-1)*k-1)*m，将a^(b-1)提出来得到a^(b-1)*(E+a^k+(a^k)^2+(a^k)^3+......+(a^k)^(n-1))*m，括号里面是一个等比矩阵求和，这部分计算方法如下：设M为某矩阵，若n为偶数，M^1+M^2+......+M^n可以分成两半，一部分是M^1+M^2+......+M^(n/2)，另一部分是M^(n/2+1)+......+M^(n-1)+M^n，可以从后面那部分提取出M^(n/2)，然后两部分都可以提出M^1+M^2+......+M^(n/2)，最终式子转换为(E+M^(n/2))*(M^1+M^2+......+M^(n/2))，(E+M^(n/2))就是答案的一部分，而(M^1+M^2+......+M^(n/2))则是接下来要继续求解的部分，可以发现这会将n次矩阵加法转化为logn次矩阵加法，通过这种方法就可求解等比矩阵的和了，然后将结果代入一开始的表达式就可以得到答案了。要注意给出的参数从0~1e9，对于边界情况一定要仔细斟酌！

具体代码如下： 

#include <iostream>
#include <cstdio>
#include <algorithm>
#include <cstring>
#include <cmath>
#include <string>
#define int long long
using namespace std;

const int N = 3;
int E[N][N] = {0, 0, 0, 0, 1, 0, 0, 0, 1}, a[N][N] = {0, 0, 0, 0, 1, 1, 0, 1, 0};
int mod;

void cpy(int a[N][N], int b[N][N]){
	for(int i = 1; i < N; i++)
		for(int j = 1; j < N; j++)
			a[i][j] = b[i][j];
}

void mul(int a[N][N], int b[N][N]){
	int temp[N][N];
	memset(temp, 0, sizeof temp);
	for(int i = 1; i < N; i++)
		for(int j = 1; j < N; j++)
			for(int k = 1; k < N; k++)
				temp[i][j] = (temp[i][j]+a[i][k]*b[k][j])%mod;
	cpy(a, temp);
}

void add(int a[N][N], int b[N][N]){
	for(int i = 1; i < N; i++)
		for(int j = 1; j < N; j++)
			a[i][j] = (a[i][j]+b[i][j])%mod;
}

void ksm(int base[N][N], int power){
	int ans[N][N];
	ans[1][1] = 1%mod, ans[1][2] = 0, ans[2][1] = 0, ans[2][2] = 1%mod;
	while(power){
		if(power&1) mul(ans, base);
		mul(base, base);
		power>>=1;
	}
	cpy(base, ans);
}


void sum(int a[N][N], int x){//a^1+a^2+......+a^x 
	int ans[N][N], sum[N][N], temp[N][N];
	memset(sum, 0, sizeof sum);
	if(x == 0){
		cpy(a, sum);
		return;
	}
	ans[1][1] = 1, ans[1][2] = 0, ans[2][1] = 0, ans[2][2] = 1;
	while(x){
		if(x == 1){
			mul(ans, a);
			break;
		}
		if(x&1){
			cpy(temp, a);
			ksm(temp, x);
			int t[N][N];
			cpy(t, ans);
			mul(t, temp);
			add(sum, t);
		}
		if(x > 1){
			cpy(temp, a);
			ksm(temp, x/2);
			add(temp, E);
			mul(ans, temp);
		}
		x >>= 1;
	}
	add(ans, sum);
	cpy(a, ans);
} 

signed main()
{
	int k, b, n;
	while(~scanf("%lld%lld%lld%lld", &k, &b, &n, &mod)){
		if(b == 0) b = k, n--;
		if(b == 0){
			puts("0");
			continue;
		}
		if(n < 1){
			puts("0");
			continue;
		}
		int temp1[N][N], temp2[N][N];
		cpy(temp1, a);
		cpy(temp2, a);
		ksm(temp1, b-1);
		ksm(temp2, k);
		sum(temp2, n-1);
		add(temp2, E);
		mul(temp1, temp2);
		int ans = 0;
		ans = (ans+temp1[1][1])%mod;
		printf("%lld\n", ans);
	}
	
    return 0;
}