HDU1588_Gauss Fibonacci_fib有关题目的升级版

这道题目又跟fib有关，不过比一般的fib矩阵难度提高了，关键要有比较强的数学思想，

还要掌握牢固的矩阵公式。

/*
*State: 1588    0MS    312K    4305 B    
*题目大意：
*        有等差数列：g(i)=k*i+b;
*        fib数列:
*            f(0)=0
*            f(1)=1
*            f(n)=f(n-1)+f(n-2) (n>=2)
*        有 k,b,n ,calculate the sum of every f(g(i)) for 0<=i<n
*解题思路：
*        如下。其实挺简单。
*解题感想：
*        题目还有困惑，就是最后怎么是取ans(0, 1)?
*/
/*
把斐波那契数列转化为矩阵:
A={1 1}
  {1,0};
{f[n+1],f[n]}   
{f[n],f[n-1]} = A^n ;最后输出右上角那项
或者用
{f[n+2],f[n+1]}
{f[n+1], f[n] } =  A^(n+1); 最后输出右下角那项
我们用第一个公式
所求即为A^b + A^(k+b) + A^(2*k+b) + ... + A^((n-1)*k+b)
=A^b * ( A^0 + A^k + A^(2*k) + ... + A^((n-1)*k) )
=A^b * ( (A^k)^0 + (A^k)^1 + (A^k)^2 + ...+ (A^k)^(n-1) );
B=A^k;
上式=A^b * ( B^0 + B^1 + B^2 + ... + B^(n-1) );
B^0 + B^1 + B^2 + ... + B^(n-1)用上篇介绍到的递归二分 方法求解
最后输出矩阵的第二项(右上角)即可；

对于求解 B^0 + B^1 + B^2 + ... + B^(n-1) 
我们也可以构造矩阵的方法。
我们来设置这样一个矩阵
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
用快速幂求出n方，让第二项再与A^b相乘即可。
*/

#include <stdio.h>
#include <string.h>
#include <iostream>
#define MAX_DIMENSION 5 

using namespace std;

typedef int MATRIX_TYPE;
typedef __int64 MAX_INT_TYPE; //temporary var to avoid overflowing
typedef MATRIX_TYPE Matrix[MAX_DIMENSION][MAX_DIMENSION];
int ndim=MAX_DIMENSION;
int mod=1;//mod

void m_zero(Matrix  x)
{
    memset(x, 0, sizeof(MATRIX_TYPE)*MAX_DIMENSION*MAX_DIMENSION);
}

void m_one(Matrix  x)
{
    int i;
    m_zero(x);
    for(i=0;i<ndim;++i)x[i][i]=1;
}

void m_copy(Matrix  src,Matrix  dest)
{
    memcpy(dest,src, sizeof(MATRIX_TYPE)*MAX_DIMENSION*MAX_DIMENSION);
}

//z=x+y;
void m_add(Matrix  x,Matrix  y,Matrix  z)
{
    int i,j;
    for(i=0;i<ndim;i++)
        for(j=0;j<ndim;j++)
            if(mod<=1)z[i][j]=x[i][j]+y[i][j];
            else z[i][j]=(x[i][j]+(MAX_INT_TYPE)y[i][j])%mod;//module
}


//c=a*b
void m_multiple(Matrix  a,Matrix b,Matrix c)
{
    int i,j,k;
    MAX_INT_TYPE t;

    for(i=0;i<ndim;i++)
        for(j=0;j<ndim;j++)
        {
            t=0;
            if(mod<=1)
                for(k=0;k<ndim;k++) t+=a[i][k]*b[k][j];//module
            else
                for(k=0;k<ndim;k++){
                    t+=(a[i][k]*(MAX_INT_TYPE)b[k][j])%mod;
                    t%=mod;
                }//module
            c[i][j]=t;
        }
}

void m_pow(Matrix x, unsigned int n, Matrix y)
{
    Matrix temp,temp_x;
    m_one(y);
    m_copy(x,temp_x);
    for(;n;)
    {
        if ((n & 1) != 0)
        {
            m_multiple(temp_x,y,temp);
            m_copy(temp,y);
        }
        if ((n >>= 1))
        {
            m_multiple(temp_x,temp_x,temp);
            m_copy(temp,temp_x);
        }
    }
}


void m_pow_with_saved(Matrix  x_p[],unsigned int n, Matrix y)
{
    Matrix temp;
    m_one(y);
    for(int k=1;n;++k,n>>=1)
    {
        if ((n & 1) != 0)
        {
            m_multiple(x_p[k],y,temp);
            m_copy(temp,y);
        }
    }
}


void m_pow_sum1(Matrix  x_p[],unsigned int n, Matrix y)
{
    m_zero(y);
    if(n==0)return;
    if(n==1) m_copy(x_p[1],y);
    else
    {
        Matrix temp;
        //calculate x^1+...+^(n/2)
        m_pow_sum1(x_p,n>>1,temp);
        //add to y
        m_add(temp,y,y);
        //calculate x^(1+n/2)+...+x^n
        Matrix temp2;
        m_pow_with_saved(x_p,n>>1,temp2);
        Matrix temp3;
        m_multiple(temp,temp2,temp3);
        //add to y
        m_add(temp3,y,y);
        if(n&1)
        {
            //calculate x^(n-1)
            m_multiple(temp2,temp2,temp3);
            //calculate x^n
            m_multiple(temp3,x_p[1],temp2);
            //add x^n
            m_add(temp2,y,y);
        }
    }

}

void m_pow_sum(Matrix x, unsigned int n, Matrix y)
{
    unsigned int i=0,logn,n2=n;
    for(logn=0,n2=n;n2;logn++,n2 >>= 1);
    Matrix *pow_arry=new Matrix[logn==0?2:(logn+1)];
    m_one(pow_arry[0]);
    m_copy(x,pow_arry[1]);
    for(i=1;i<logn;i++)
    {
        m_multiple(pow_arry[i],pow_arry[i],pow_arry[i+1]);
    }

    m_pow_sum1(pow_arry,n,y);
    delete []pow_arry;
}

void fib(int m)
{
    Matrix a={{1,1},{1,0}},y;
    ndim=2;
    mod=1024*1024;

    m_pow(a,m,y);
    //fib(m+1)
    printf("fib(%d)=%d\n",m,y[0][0]);

    m_pow_sum(a,m,y);
    //fib(m+1)
    printf("sum(fib(i))=%d,i=1,..,%d\n",y[0][0],m);

}

void view_arr(Matrix a, int n)
{
    for(int i = 0; i < n; i++)
    {
        for(int j = 0; j < n;j++)
            cout << a[i][j] << " ";
        cout << endl;
    }
}

int main()
{
#ifndef ONLINE_JUDGE
    freopen("in.txt", "r", stdin);
#endif
    Matrix a = {{1, 1}, {1, 0}};
    //view_arr(a, 2);
    ndim = 2;
    int k, b, n, M;
    while(scanf("%d %d %d %d", &k, &b, &n, &M) == 4)
    {
        mod = M;

        Matrix a_k, a_b, a_k_sum, res;
        m_one(res);
        m_pow(a, k, a_k);
        m_pow(a, b, a_b);
        m_pow_sum(a_k, n - 1, a_k_sum);
        Matrix addans, ans;
        m_add(res, a_k_sum, addans);
        m_multiple(addans, a_b, ans);
        printf("%d\n", ans[0][1]);
        //view_arr(ans, 2);
        //cout << endl;
    }
    return 0;
}