hdu1588 Gauss Fibonacci （矩阵快速幂）

P.S. 感谢大神，借鉴于此：http://blog.csdn.net/shiyuankongbu/article/details/8458459  讲的真心很明白，，， 
     

题目大意：给定k, n, b, M , 求斐波那契数列第b项到第 k*i+b 项（ 0<=i<n）的和sum(n) % M的值。

解题思路：对于F(b) + F(k+b) + F(2*k+b) + … + F( (n-1) * k+b )   我们用矩阵快速幂处理 可得：T = B^b + B^(k+b) + B^(2*k+b) +…+ B^((n-1)*k +b)  = B^b * [ E + B^k + B^(2*k) +…+B^((n-1)*k) ] (其中E为单位阵)；然后中括号里面用二分解决，，，很麻烦，，，

代码如下：

#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cmath>
#include <algorithm>
using namespace std;
typedef struct
{
long long m[2][2];
}Matrix;

Matrix B, E;
long long n, k, b, M;
void init()
{
B.m[0][0] = B.m[0][1] = B.m[1][0] = 1;
B.m[1][1] = 0;
E.m[0][0] = E.m[1][1] = 1;
E.m[0][1] = E.m[1][0] = 0;
}

Matrix multi(Matrix a, Matrix b)
{
Matrix c;
for(int i = 0; i < 2; i ++)
{
for(int j = 0; j < 2; j ++)
{
c.m[i][j] = 0;
for(int k = 0; k < 2; k ++)
{
c.m[i][j] = c.m[i][j] + a.m[i][k] * b.m[k][j] % M;
}
c.m[i][j] = c.m[i][j] % M;
}
}
return c;
}

Matrix Power(Matrix A, int k)
{
Matrix ans = E;
while(k)
{
if(k & 1) ans = multi(ans, A);
k = k / 2;
A = multi(A, A);
}
return ans;
}

Matrix Add(Matrix a, Matrix b)
{
Matrix c;
for(int i = 0; i < 2; i ++)
{
for(int j = 0; j < 2; j ++)
{
c.m[i][j] = (a.m[i][j] + b.m[i][j])% M;
}
}
return c;
}

Matrix Sum(Matrix a, long long n)
{
Matrix ans;
if(n == 0)
{
ans.m[0][0] = ans.m[0][1] = ans.m[1][0] = ans.m[1][1] = 0;
return ans;
}
if(n == 1)
return a;
else
{
// cnt=A^k
// sum(n) = A^k+A^2k+...A^(n-1)k
// = cnt+cnt^2+cnt^3+...cnt^(n-1)
// t=sum(n/2)
// sum(n) = sum(n/2) + cnt^(n/2+1)+cnt^(n/2+2)+...cnt^(n-1)
// sum(n) = t + cnt^(n/2+1)+cnt^(n/2+2)+...cnt^(n-1)
// sum(n) = t + t * cnt^(n/2)
// 递归到最后的时候，n为奇数时直接算，相当于 n=39时
// sum(39) = cnt^1+cnt^2+...+cnt^39 t=19
// sum(39) = sum(19)+sum(19)*cnt^19+ cnt^39 (最后一项是因为n为奇数所余留的一项)
Matrix res, t;
t = Sum(a, n/2);
res = Add(t, multi(Power(a, n/2), t));
if(n & 1)
res = Add(res, Power(a, n));
return res;
}
}
int main()
{
while(~scanf("%lld %lld %lld %lld", &k, &b, &n, &M) && k + b + n + M)
{
init();
Matrix M1 = Power(B, k);
Matrix M2 = Power(B, b);
Matrix res = multi(M2, Add(E, Sum(M1, n-1)));
printf("%lld\n", res.m[0][1]);
}
return 0;
}