HDU 1588 Gauss Fibonacci（矩阵高速幂+二分等比序列求和）
ACM
题意：
g(i)=k*i+b;i为变量。
给出k,b,n,M，问( f(g(0)) + f(g(1)) + ... + f(g(n)) ) % M的值。
分析：
把斐波那契的矩阵带进去，会发现这个是个等比序列。
推倒：
S(g(i))= F(b) + F(b+k) + F(b+2k) + .... + F(b+nk)// 设 A = {1,1,0,1}, （花括号表示矩阵...）// 也就是fib数的变化矩阵，F(x) = (A^x) * {1,0}= F(b) + (A^k)F(b) + (A^2k)F(b)+….+(A^nk)F(b)// 提取公因式 F(b)= F(b) [ E +A^k + A^2k + ….+ A^nk] // (E表示的是单位矩阵)// 令 K = A^k 后E +A^k + A^2k + ….+ A^nk 变成 K^0+K^1+K^2+…+K^n
然后等比数列是能够二分求和的：数论_等比数列二分求和
代码：
/*
* Author: illuz <iilluzen[at]gmail.com>
* Blog: http://blog.csdn.net/hcbbt
* File: 1588.cpp
* Create Date: 2014-08-04 16:13:51
* Descripton: Matrix
*/ #include <cstdio>
#include <cstring>
#include <iostream>
#include <map>
#include <algorithm>
using namespace std;
#define repf(i,a,b) for(int i=(a);i<=(b);i++)
typedef long long ll; const int N = 20;
const int SIZE = 2; // max size of the matrix ll MOD;
ll k, b, n; struct Mat{
int n;
ll v[SIZE][SIZE]; // value of matrix Mat(int _n = SIZE) {
n = _n;
} void init(ll _v = 0) {
memset(v, 0, sizeof(v));
if (_v)
repf (i, 0, n - 1)
v[i][i] = _v;
} void output() {
repf (i, 0, n - 1) {
repf (j, 0, n - 1)
printf("%lld ", v[i][j]);
puts("");
}
puts("");
}
} a, B, C; Mat operator * (Mat a, Mat b) {
Mat c(a.n);
repf (i, 0, a.n - 1) {
repf (j, 0, a.n - 1) {
c.v[i][j] = 0;
repf (k, 0, a.n - 1) {
c.v[i][j] += (a.v[i][k] * b.v[k][j]) % MOD;
c.v[i][j] %= MOD;
}
}
}
return c;
} Mat operator ^ (Mat a, ll k) {
Mat c(a.n);
c.init(1);
while (k) {
if (k&1) c = a * c;
a = a * a;
k >>= 1;
}
return c;
} Mat operator + (Mat a, Mat b) {
Mat c(a.n);
repf (i, 0, a.n - 1)
repf (j, 0, a.n - 1)
c.v[i][j] = (b.v[i][j] + a.v[i][j]) % MOD;
return c;
} Mat operator + (Mat a, ll b) {
Mat c = a;
repf (i, 0, a.n - 1)
c.v[i][i] = (a.v[i][i] + b) % MOD;
return c;
} // 二分求和1..n
Mat calc(Mat a, int n) {
if (n == 1)
return a;
if (n&1)
return (a^n) + calc(a, n - 1);
else
return calc(a, n/2) * ((a^(n/2)) + 1);
} int main() {
a.init();
a.v[0][0] = a.v[0][1] = a.v[1][0] = 1;
while (~scanf("%lld%lld%lld%lld", &k, &b, &n, &MOD)) {
B = (a^k);
C = calc(B, n - 1) + (a^0);
C = (a^b) * C;
printf("%lld\n", C.v[0][1]);
}
return 0;
}