HDU 1588 Gauss Fibonacci(分块矩阵优化)

Gauss Fibonacci

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

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

题意

给你$k,b,n,M$,求解${\sum }_{i=0}^{n-1}f(g(i))$,其中$g(i)=ki+b$,$f(i)$为斐波那契数列,$f(0)=0,f(1)=1,f(n)=f(n-1)+f(n-2)$

思路

先把斐波那契数列写成矩阵形式
$$\left[\begin{array}{c}f(n+1)\\ f(n)\end{array}\right]={\left[\begin{array}{cc}1& 1\\ 1& 0\end{array}\right]}^n\left[\begin{array}{c}f(1)\\ f(0)\end{array}\right]={\left[\begin{array}{cc}1& 1\\ 1& 0\end{array}\right]}^n\left[\begin{array}{c}1\\ 0\end{array}\right]$$
我们记
$$A=\left[\begin{array}{cc}1& 1\\ 1& 0\end{array}\right]$$
那么
$$f(n)=A^n.m[1][0]$$
那么我们要求
$$\sum _{i=0}^{n-1}f(g(i))=f(b)+f(k+b)+f(2k+b)+\cdots +f((n-1)k+b)$$

$$=(A^b+A^{k+b}+A^{2k+b}+\cdots +A^{(n-1)k+b}).m[1][0]$$
$$=(A^b(I+A^k+A^{2k}+\cdots +A^{(n-1)k})).m[1][0]$$
我们设$R=A^k$,那么就有
$$I+A^k+A^{2k}+\cdots +A^{(n-1)k}=I+R+R^2+\cdots +R^{n-1}$$
那么现在的问题就转化为求解上面的$R$的式子了
我们可以构造这样一个矩阵
$$\left[\begin{array}{cc}R& I\\ 0& I\end{array}\right]$$
这个矩阵有如下的性质即
$${\left[\begin{array}{cc}R& I\\ 0& I\end{array}\right]}^2=\left[\begin{array}{cc}{R}^2& RI+I\\ 0& I\end{array}\right]$$
$${\left[\begin{array}{cc}R& I\\ 0& I\end{array}\right]}^3=\left[\begin{array}{cc}{R}^3& R^{2}I+RI+I\\ 0& I\end{array}\right]$$
$$⋮$$
$${\left[\begin{array}{cc}R& I\\ 0& I\end{array}\right]}^n=\left[\begin{array}{cc}{R}^n& R^{n}I+\cdots +R^{2}I+RI+I\\ 0& I\end{array}\right]$$
$I$为单位矩阵
所以我们可以用矩阵快速幂得到我们要求的$I+R+R^2+\cdots +R^{n-1}$
然后再和$A^b$做矩阵乘法即可

#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
using namespace std;
struct Matrix
{
    long long m[2][2];
};
struct Matrix2
{
    long long m[4][4];
};
long long M;
Matrix P1= {1,1,
            1,0};
Matrix I1= {1,0,
            0,1};
Matrix2 P2 = {0,1,1,0,
            0,0,0,1,
            0,0,1,0,
            0,0,0,1
           };

Matrix2 I2 ={1,0,0,0,
            0,1,0,0,
            0,0,1,0,
            0,0,0,1
           };

Matrix Matrixmul(Matrix a,Matrix b)
{
    int i,j,k;
    Matrix c;
    for(i=0; i<2; i++)
        for(j=0; j<2; j++)
        {
            c.m[i][j]=0;
            for(k=0; k<2; k++)
            {
                c.m[i][j]+=((a.m[i][k]%M)*(b.m[k][j]%M))%M;
            }
            c.m[i][j]%=M;
        }
    return c;
}
Matrix quickpow(long long n)
{
    Matrix m=P1,b=I1;
    while(n>0)
    {
        if(n%2==1)
            b=Matrixmul(b,m);
        n=n/2;
        m=Matrixmul(m,m);
    }
    return b;
}
Matrix2 Matrixmul2(Matrix2 a,Matrix2 b)
{
    int i,j,k;
    Matrix2 c;
    for(i=0; i<4; i++)
        for(j=0; j<4; j++)
        {
            c.m[i][j]=0;
            for(k=0; k<4; k++)
            {
                c.m[i][j]+=((a.m[i][k]%M)*(b.m[k][j]%M))%M;
            }
            c.m[i][j]%=M;
        }
    return c;
}
Matrix2 quickpow2(long long n)
{
    Matrix2 m=P2,b=I2;
    while(n>0)
    {
        if(n%2==1)
            b=Matrixmul2(b,m);
        n=n/2;
        m=Matrixmul2(m,m);
    }
    return b;
}
int main()
{
    long long k,b,n;
    while(scanf("%lld%lld%lld%lld",&k,&b,&n,&M)!=EOF)
    {
        Matrix B=quickpow(b);
        Matrix R=quickpow(k);
        P2.m[0][0]=R.m[0][0];
        P2.m[0][1]=R.m[0][1];
        P2.m[1][0]=R.m[1][0];
        P2.m[1][1]=R.m[1][1];
        Matrix2 A=quickpow2(n);
        R.m[0][0]=A.m[0][2];
        R.m[0][1]=A.m[0][3];
        R.m[1][0]=A.m[1][2];
        R.m[1][1]=A.m[1][3];
        R=Matrixmul(R,B);
        printf("%lld\n",R.m[1][0]);
    }
    return 0;
}