HDU 3306 Another kind of Fibonacci （矩阵快速幂）

题目：http://acm.hdu.edu.cn/showproblem.php?pid=3306

题意：A(N) = X * A(N - 1) + Y * A(N - 2) (N >= 2). S(N) = A(0)² +A(1)²+……+A(n)^{2. 给出N，X，Y求S(N)}

^思路：

首先有这些关系：

S(n) = S(n -1) + A(n)^2

A(n) = xA(n - 1) + yA(n - 2)

A(n)^2 = xxA(n - 1)A(n-1) + 2xyA(n - 1)A(n - 2) + yyA(n - 2)A(n - 2)

那么如果我们有一个矩阵1

|S(n - 1)  A(n)^2      A(n)A(n - 1)  A(n - 1)^2|

我们想得到矩阵2

| S(n)   A(n + 1)^2   A(n + 1)A(n)     A(n)^2 |

那么，根据矩阵的乘法，我们应该让矩阵1乘以这样一个矩阵3

  

利用了关系：

S(n) = S(n - 1) + A(n)^2

A(n + 1)^2 = xxA(n)^2 + 2xyA(n)A(n - 1) + yyA(n - 1)^2

      A(n + 1)A(n) = xA(n)^2 + yA(n)A(n - 1)

所以

| S(n)   A(n + 1)^2   A(n + 1)A(n)     A(n)^2 | = |S(0) A(1)^2 A(1)A(0) A(0)^2| *（上述矩阵3的n次方）

#include <stdio.h>
#include <string.h>
#include <algorithm>
#include <math.h>
#include <stdlib.h>
#include <vector>
#define INF 0x7fffffff
#define MOD 10007
using namespace std;
typedef long long ll;
typedef vector<int> vec;
typedef vector<vec> mat;
int size;

mat mul(mat &A, mat &B)
{
    mat C(size, vec(size));
    for(int i = 0; i < size; i++)
    {
        for(int k =  0; k < size; k++)
        {
            for(int j = 0; j < size; j++)
            {
                C[i][j] = (C[i][j] + A[i][k] * B[k][j]) % MOD;
            }
        }
    }
    return C;
}

mat pow(mat A, int n)
{
    mat B(size, vec(size));
    for(int i = 0; i < size; i++)
        B[i][i] = 1;
    while(n > 0)
    {
        if(n & 1) B = mul(B, A);
        A = mul(A, A);
        n >>= 1;
    }
    return B;
}

int main()
{
    #ifdef LOCAL
    freopen("data.in", "r", stdin);
    #endif

    size = 4;
    int n, x, y;
    while(scanf("%d%d%d", &n, &x, &y) != EOF)
    {
        x %= MOD; y %= MOD;
        if(n == 0)
        {
            printf("1\n");
            continue;
        }
        mat A(4, vec(4));
        A[0][0] = 1; A[0][1] = 0; A[0][2] = 0; A[0][3] = 0;
        A[1][0] = 1; A[1][1] = x * x % MOD; A[1][2] = x; A[1][3] = 1;
        A[2][0] = 0; A[2][1] = 2 * x * y % MOD; A[2][2] = y; A[2][3] = 0;
        A[3][0] = 0; A[3][1] = y * y % MOD; A[3][2] = 0; A[3][3] = 0;
        A = pow(A, n);
        /*for(int i = 0; i < 4; i++)
        {
            for(int j = 0; j < 4; j++)
                printf("%d ", A[i][j]);
            printf("\n");
        }*/
        int ans = 0;
        for(int i = 0; i < 4; i++)
            ans += A[i][0];
        printf("%d\n", ans % MOD);
    }

    return 0;
}