hoj2060-Fibonacci Problem Again（矩阵快速幂模板）

Fibonacci Problem Again

My Tags	(Edit)

	Source : HCPC 2005 FALL
	Time limit : 1 sec		Memory limit : 32 M

Submitted : 380, Accepted : 137

As we know , the Fibonacci numbers are defined as follows:

 """"

Given two numbers a and b , calculate . """"

Input

The input contains several test cases. Each test case consists of two non-negative integer numbers a and b (0 ≤ a ≤ b ≤1,000,000,000). Input is terminated by a = b = 0.

Output

For each test case, output S mod 1,000,000,000, since S may be quite large.

Sample Input

1 1
3 5
10 1000
0 0

Sample Output

1
16
496035733

#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <vector>
#include <queue>
#include <algorithm>
#define LL long long
#define M 1000000000
using namespace std;
struct Matrix
{
    LL A[2][2];
    int size;
    Matrix() { memset(this, 0, sizeof(*this)); }
};

LL mymod(LL x)
{
    return (x % M + M) % M;
}
Matrix operator+(Matrix m1, Matrix m2) {
    Matrix ret;
    ret.size = m1.size;
    for(int i = 0; i < ret.size; i++)
        for(int j = 0; j < ret.size; j++)
            ret.A[i][j] = mymod(m1.A[i][j] + m2.A[i][j]);
    return ret;
}

Matrix operator*(Matrix m1, Matrix m2) {
    Matrix ret;
    ret.size = m1.size;
    for(int i = 0; i < ret.size; i++)
        for(int j = 0; j < ret.size; j++) {
            ret.A[i][j] = 0;
            for(int k = 0; k < ret.size; k++){
                ret.A[i][j] += m1.A[i][k] * m2.A[k][j];
                ret.A[i][j] = mymod(ret.A[i][j]);
            }
           }
    return ret;
}

Matrix mypower(Matrix m, LL n) {
    Matrix ret, tmp;
    ret.size = m.size;
    if(n == 0) {
        for(int i = 0; i < ret.size; i++)
            ret.A[i][i] = 1 % M;
        return ret;
    }
    tmp = mypower(m, n / 2);
    if(n & 1)
        return tmp * tmp * m;
    else return tmp * tmp;
}

Matrix sumpower(Matrix m, LL n)
{
    Matrix tmp;
    if(n == 1) return m;
    tmp = sumpower(m, n / 2);
    if(n & 1)
        return mypower(m, n / 2) * tmp + tmp + mypower(m, n);
    return mypower(m, n / 2) * tmp + tmp;
}

int main()
{
    LL a,b,f1,f2;
    Matrix m;
    while(scanf("%lld%lld",&a,&b)==2&&(a||b)){
        m.A[0][0]=1;m.A[0][1]=1;m.A[1][0]=1;m.A[1][1]=0;m.size=2;
        if(a==0)f2=0;
        else{
            Matrix ans2=sumpower(m,a);
            f2=ans2.A[1][0];
        }
        Matrix ans1=sumpower(m,b+1);
        f1=ans1.A[1][0];
        if(f1<f2)printf("%lld\n",f1-f2+1000000000);
        else printf("%lld\n",f1-f2);
    }
    return 0;
}

这道题可用作矩阵快速幂的模板，mypower是求快速幂，sumpower是求和，即s[n]

构造结构体：

struct Matrix 
{
    long long A[N][N];
    int size;
    Matrix() { memset(this, 0, sizeof(*this)); }
};

基本矩阵操作：

+：

Matrix operator+(Matrix m1, Matrix m2) {
    Matrix ret;
    ret.size = m1.size;
    for(int i = 0; i < ret.size; ++ i)
        for(int j = 0; j < ret.size; ++ j)
            ret.A[i][j] = mymod(m1.A[i][j] + m2.A[i][j]);
    return ret;
}

-:

Matrix operator-(Matrix m1, Matrix m2) {
    Matrix ret;
    ret.size = m1.size;
    for(int i = 0; i < ret.size; ++ i)
        for(int j = 0; j < ret.size; ++ j)
            ret.A[i][j] = mymod(m1.A[i][j] - m2.A[i][j]);
    return ret;
}

P.S.: 关于可能出现的负数取模问题:

long long mymod(long long x)
 {
    return (x % M + M) % M;
}

*:

Matrix operator*(Matrix m1, Matrix m2) {
    Matrix ret;
    ret.size = m1.size;
    for(int i = 0; i < ret.size; ++ i)
        for(int j = 0; j < ret.size; ++ j) {
            ret.A[i][j] = 0;
            for(int k = 0; k < ret.size; ++ k){
                ret.A[i][j] += m1.A[i][k] * m2.A[k][j];
                ret.A[i][j] = mymod(ret.A[i][j]);
            }
       }
    return ret;
}

快速幂：

Matrix mypower(Matrix m, int n) {
    Matrix ret, tmp;
    ret.size = m.size;
    if(n == 0) {
        for(int i = 0; i < ret.size; ++ i)
            ret.A[i][i] = 1 % M;
        return ret;
    }
    tmp = mypower(m, n / 2);
    if(n & 1)
        return tmp * tmp * m;
    else return tmp * tmp;
}

求和快速幂：

Matrix sumpower(Matrix m, int n) 
{
    Matrix tmp;
    if(n == 1) return m;
    tmp = sumpower(m, n / 2);
    if(n & 1)
        return mypower(m, n / 2) * tmp + tmp + mypower(m, n);
    return mypower(m, n / 2) * tmp + tmp;
}