Gauss Fibonacci

Time Limit: 1000/1000 MS (Java/Others)    Memory Limit: 32768/32768 K (Java/Others)
Total Submission(s): 3149    Accepted Submission(s): 1323

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

这个题要利用的东西挺多的。刚开始连斐波拉契的矩阵都弄错了，看来还是要加强训练！

把斐波那契数列转化为矩阵:
A={1 1}
  {1,0};

{f[n+1],f[n]}  
{f[n],f[n-1]} = A^n ;最后输出M.v[1][0]这就是构造斐波拉契数列的矩阵了。

这个题将式子化简之后我们可以得到

sum(f(g(n))) = A^(b)+A^(b+k)+...+A^(b+(n-1)*k) = A^b * (1+A^k+A^2k+...A^(n-1)k)

利用递归函数求解 A^k+A^2k+...A^(n-1)k 然后与 单位矩阵相加，最后还要判断 b是否为 0

很巧妙。

这里给出一个等比数列的求和公式.

LL cal(int p,int n){  ///这里是递归求解等比数列模板 1+p+p^2...+p^n
    if(n==0) return 1;
    if(n&1){//(1+p+p^2+....+p^(n/2))*(1+p^(n/2+1));
         return (1+pow_mod(p,n/2+1))*cal(p,n/2)%mod;
    }
    else { //(1+p+p^2+....+p^(n/2-1))*(1+p^(n/2+1))+p^(n/2);
         return (pow_mod(p,n/2)+(1+pow_mod(p,n/2+1))*cal(p,n/2-1))%mod;
    }
}

#include<stdio.h>
#include<iostream>
#include<string.h>
#include <stdlib.h>
#include<math.h>
#include<algorithm>
using namespace std;
typedef long long LL;
struct Martix{
    LL v[2][2];
}res;
LL k,b,n,M;

Martix mult(Martix a,Martix b){
    Martix temp;
    for(int i=0;i<2;i++){
        for(int j=0;j<2;j++){
            temp.v[i][j] = 0;
            for(int k=0;k<2;k++){
                temp.v[i][j] = (temp.v[i][j]+a.v[i][k]*b.v[k][j])%M;
            }
        }
    }
    return temp;
}

Martix add(Martix a,Martix b){
    for(int i=0;i<2;i++){
        for(int j=0;j<2;j++){
            a.v[i][j]=(a.v[i][j]+b.v[i][j])%M;
        }
    }
    return a;
}
Martix pow_mod(Martix a,LL k){
    Martix ans;
    ans.v[0][0]=ans.v[1][1] = 1;
    ans.v[0][1]= ans.v[1][0]=0;
    while(k){
        if(k&1) ans = mult(ans,a);
        a=mult(a,a);
        k>>=1;
    }
    return ans;
}

Martix cal(Martix p,LL k) ///用二分法求矩阵和,递归实现 计算 a^1+a^2.....+a^p
{
    if(k==1)
        return p;
    else if(k&1)
        return add(cal(p,k-1),pow_mod(p,k));
    else
        return mult(cal(p,k>>1),add(pow_mod(p,k>>1),res));
}

int main(){

    Martix a,t;
    a.v[0][0] = a.v[0][1] = a.v[1][0] = 1; ///斐波拉契数列的特征值矩阵[1 1 1 0]
    a.v[1][1] = 0;
    t.v[0][0] = t.v[1][1] = 1; ///单位矩阵
    t.v[0][1] = t.v[1][0] = 0;
    while(scanf("%lld%lld%lld%lld",&k,&b,&n,&M)!=EOF){
        Martix M1 = pow_mod(a,k);
        res = t;
        res = add(t,cal(M1,n-1));
        if(b!=0){
            Martix M3 = pow_mod(a,b);
            res = mult(res,M3);
        }
        printf("%d\n",res.v[0][1]%M);
    }
    return 0;
}