Gauss Fibonacci

Problem Description
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
2 1 4 100
2 0 4 100
Sample Output
21

12

题目中的公式：

f(0)=0 f(1)=1
f(n)=f(n-1)+f(n-2) (n>=2)

g(i)=k*i+b;

对于Fib序列:

g(i)=k*i+b 是等差数列

所以F(g(i)) = F(b) + F(b+k)+F(b+2k)+....+F(b+nk)

           = 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

构造矩阵

<pre name="code" class="cpp">#include<stdio.h>#include<string.h>__int64 k,b,n,m;struct Matrix{__int64 map[2][2];};Matrix A,B;struct Matrix_4{__int64 map[4][4];};Matrix Add(const Matrix &c, const Matrix &d) { Matrix res; int i,j;for(i=0;i<2;i++)for(j=0;j<2;j++)res.map[i][j] = (c.map[i][j]+d.map[i][j]) % m; return res; } Matrix_4 Add_4(const Matrix_4 &c, const Matrix_4 &d) { Matrix_4 res; int i,j;for(i=0;i<4;i++)for(j=0;j<4;j++)res.map[i][j] = (c.map[i][j]+d.map[i][j]) % m; return res; } Matrix Mul(Matrix c,Matrix d)//矩阵的乘法{Matrix tmp;memset(tmp.map,0,sizeof(tmp.map));int i,j,k;for(i=0;i<2;i++)for(j=0;j<2;j++)for(k=0;k<2;k++)tmp.map[i][j]=(tmp.map[i][j]+c.map[i][k]*d.map[k][j]%m);return tmp;}Matrix_4 Mul_4(Matrix_4 c,Matrix_4 d)//矩阵的乘法{Matrix_4 tmp;memset(tmp.map,0,sizeof(tmp.map));int i,j,k;for(i=0;i<4;i++)for(j=0;j<4;j++)for(k=0;k<4;k++)tmp.map[i][j]=(tmp.map[i][j]+c.map[i][k]*d.map[k][j]%m);return tmp;}Matrix pow_mod(Matrix c,__int64 cnt)//矩阵的快速幂{int i;Matrix tmp;memset(tmp.map,0,sizeof(tmp.map));for(i=0;i<2;i++)tmp.map[i][i]=1;if(cnt==0) return tmp;if(cnt==1) return c;tmp=pow_mod(c,cnt/2);tmp=Mul(tmp,tmp);if(cnt&1) tmp=Mul(tmp,c);return tmp;}Matrix_4 pow_mod_4(Matrix_4 c,__int64 cnt)//矩阵的快速幂{int i;Matrix_4 tmp;memset(tmp.map,0,sizeof(tmp.map));for(i=0;i<4;i++)tmp.map[i][i]=1;if(cnt==0) return tmp;if(cnt==1) return c;tmp=pow_mod_4(c,cnt/2);tmp=Mul_4(tmp,tmp);if(cnt&1) tmp=Mul_4(tmp,c);return tmp;}int main(){//freopen("b.txt","r",stdin);while(scanf("%I64d %I64d %I64d %I64d",&k,&b,&n,&m)==4){A.map[0][0]=1,A.map[0][1]=1,A.map[1][0]=1,A.map[1][1]=0;if(n==0) {printf("%d\n",0);continue;}B=pow_mod(A,k);Matrix_4 EB;memset(EB.map,0,sizeof(EB.map));EB.map[0][0]=B.map[0][0];EB.map[0][1]=B.map[0][1];EB.map[1][0]=B.map[1][0];EB.map[1][1]=B.map[1][1];EB.map[0][2]=EB.map[1][3]=1;EB.map[2][2]=EB.map[3][3]=1;EB=pow_mod_4(EB,n);B.map[0][0]=EB.map[0][2];B.map[0][1]=EB.map[0][3];B.map[1][0]=EB.map[1][2];B.map[1][1]=EB.map[1][3];A=Mul(pow_mod(A,b),B);printf("%d\n",A.map[1][0]%m);}return 0;}

#include<stdio.h>#include<string.h>__int64 k,b,n,m;struct Matrix{__int64 map[2][2];};Matrix ma,f;Matrix Add(const Matrix &c, const Matrix &d)  {      Matrix res;  int i,j;for(i=0;i<2;i++)for(j=0;j<2;j++)        res.map[i][j] = (c.map[i][j]+d.map[i][j]) % m;      return res;  }  Matrix Mul(Matrix c,Matrix d)//矩阵的乘法{Matrix tmp;memset(tmp.map,0,sizeof(tmp.map));int i,j,k;for(i=0;i<2;i++)for(j=0;j<2;j++)for(k=0;k<2;k++)tmp.map[i][j]=(tmp.map[i][j]+c.map[i][k]*d.map[k][j]%m);return tmp;}Matrix pow_mod(Matrix c,__int64 cnt)//矩阵的快速幂{int i;Matrix tmp;memset(tmp.map,0,sizeof(tmp.map));for(i=0;i<2;i++)tmp.map[i][i]=1;if(cnt==0) return tmp;if(cnt==1) return c;tmp=pow_mod(c,cnt/2);tmp=Mul(tmp,tmp);if(cnt&1) tmp=Mul(tmp,c);return tmp;}Matrix cal(Matrix c, __int64 cnt)    //分治求（a^0+a^1+a^2+...+a^n）%m  {      if (cnt == 0) return pow_mod(c, 0);      __int64 ans = (cnt+1)/2;      Matrix o = cal(c, ans-1);      Matrix res =Add(o, Mul(pow_mod(c, ans), o));      if (cnt % 2 == 0) res =Add(res, pow_mod(c, cnt));      return res;  }  int main(){freopen("b.txt","r",stdin);while(scanf("%I64d %I64d %I64d %I64d",&k,&b,&n,&m)==4){ma.map[0][0]=1,ma.map[0][1]=1,ma.map[1][0]=1,ma.map[1][1]=0;if(n==0) {printf("%d\n",0);continue;}ma = Mul(pow_mod(ma, b), cal(pow_mod(ma, k), n-1));printf("%I64d\n",ma.map[0][1]%m);}return 0;}

#include<stdio.h>#include<string.h>__int64 k,b,n,m;struct Matrix{    __int64 map[2][2];};Matrix ma,f;Matrix Add(const Matrix &c, const Matrix &d)  {      Matrix res;      int i,j;    for(i=0;i<2;i++)        for(j=0;j<2;j++)        res.map[i][j] = (c.map[i][j]+d.map[i][j]) % m;      return res;  }  Matrix Mul(Matrix c,Matrix d)//矩阵的乘法{    Matrix tmp;    memset(tmp.map,0,sizeof(tmp.map));    int i,j,k;    for(i=0;i<2;i++)        for(j=0;j<2;j++)            for(k=0;k<2;k++)                tmp.map[i][j]=(tmp.map[i][j]+c.map[i][k]*d.map[k][j]%m);    return tmp;}Matrix pow_mod(Matrix c,__int64 cnt)//矩阵的快速幂{    int i;    Matrix tmp;    memset(tmp.map,0,sizeof(tmp.map));    for(i=0;i<2;i++)        tmp.map[i][i]=1;    if(cnt==0) return tmp;    if(cnt==1) return c;    tmp=pow_mod(c,cnt/2);    tmp=Mul(tmp,tmp);    if(cnt&1) tmp=Mul(tmp,c);    return tmp;}Matrix cal(Matrix c, __int64 cnt)    //分治求（a^0+a^1+a^2+...+a^n）%m  {      if (cnt == 0) return pow_mod(c, 0);      if(cnt==1) return c;    Matrix res=cal(c,cnt/2);    if(cnt%2)    {        Matrix tmp=pow_mod(c,cnt/2+1);        res=Add(res,Mul(res,tmp));        res=Add(res,tmp);    }    else         res=Add(res,Mul(res,pow_mod(c,cnt/2)));    return res;  }  int main(){    //freopen("b.txt","r",stdin);    while(scanf("%I64d %I64d %I64d %I64d",&k,&b,&n,&m)==4)    {        ma.map[0][0]=1,ma.map[0][1]=1,ma.map[1][0]=1,ma.map[1][1]=0;        if(n==0) {printf("%d\n",0);continue;}        ma = Mul(pow_mod(ma, b), Add(pow_mod(ma,0),cal(pow_mod(ma, k), n-1)));        printf("%I64d\n",ma.map[0][1]%m);    }    return 0;}