2 1 4 100
2 0 4 100 

Sample Output

 21
12 

题意：g(i)=k*i+b , 求所有的f(g(i))在n的范围里

思路：借鉴：构造矩阵：
       |1  1|      | f(2)    f(1)|
    A= |1  0|   =  | f(1)    f(0)|
         
           
     |1  1| ^b   | f(b+1)    f(b)|             
A^b =|1  0|  =   | f(b)    f(b-1)|
                 
    f(b) = matrix[0][1]=matrix[1][0];
                       
    首项是：A^b
    公比是：A^k
    项数是：N
    可以把问题进一步简化
    因为矩阵的加法对乘法也符合分配律，我们提出一个A^b来，形成这样的式子：
    A^b*( I + A^k + (A^k)^2 + .... + (A^k)^(N-1) )
    A^b 和 A^k 显然都可以用我们之前说过的方法计算出来，这剩下一部分累加怎么解决呢
    设A^k=B
    要求 G(N)=I + ... + B^(N-1)，
    i=N/2
    若N为偶数，G(N)=G(i)+G(i)*B^i = G(i) *( I+B^(i));
    若N为奇数，G(N)=I+ G(i)*B + G(i) * (B^(i+1)) = G(N-1)+B^N; (前一个等式可能要快点，但是后面更简练)
                                   
    我们来设置这样一个矩阵
    B I
    O I
    其中O是零矩阵，I是单位矩阵
    将它乘方，得到
    B^2 I+B
    O   I
    乘三方，得到
    B^3 I+B+B^2
    O   I
    乘四方，得到
    B^4  I+B+B^2+B^3
    O    I                                 
既然已经转换成矩阵的幂了，继续用我们的二分或者二进制法，直接求出幂就可以了

#include <iostream>
#include <cstring>
#include <cstdio>
#include <algorithm>
#include <cmath>
typedef long long ll;
using namespace std;
const int maxn = 2;

int m;
struct Matrix {
	ll v[maxn][maxn];
	Matrix() {}
	Matrix(int x) {
		init();
		for (int i = 0; i < maxn; i++) 
			v[i][i] = x;
	}
	void init() {
		memset(v, 0, sizeof(v));
	}
	Matrix operator *(Matrix const &b) const {
		Matrix c;
		c.init();
		for (int i = 0; i < maxn; i++)
			for (int j = 0; j < maxn; j++)
				for (int k = 0; k < maxn; k++)
					c.v[i][j] = (c.v[i][j] + (v[i][k]*b.v[k][j])%m) % m;
		return c;
	}
	Matrix operator ^(int b) {
		Matrix a = *this, res(1);
		while (b) {
			if (b & 1)
				res = res * a;
			a = a * a;
			b >>= 1;
		}
		return res;
	}
} u, em;

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

Matrix BinarySum(Matrix a, int n) {
	if (n == 1)
		return a;
	if (n & 1)
		return Add(BinarySum(a, n-1), a^n);
	else return BinarySum(a, n>>1) * Add(u, a^(n>>1));
}

int main() {
	int k, b, n;
	u.init(), em.init();
	u.v[0][0] = 1, u.v[0][1] = 0, u.v[1][0] = 0, u.v[1][1] = 1;
	em.v[0][0] = 1, em.v[0][1] = 1, em.v[1][0] = 1, em.v[1][1] = 0;
	while (scanf("%d%d%d%d", &k, &b, &n, &m) != EOF) {
		Matrix t1, t2, ans;
		t1 = em^b;
		t2 = em^k;
		ans = Add(u, BinarySum(t2, n-1)) * t1;
		cout << ans.v[0][1] << endl;
	}
	return 0;
}