hdu 2256 Problem of Precision -矩阵快速幂

Input

　　The first line of input gives the number of cases, T. T test cases follow, each on a separate line. Each test case contains one positive integer n. (1 <= n <= 10^9) 
Output

　　For each input case, you should output the answer in one line. 
Sample Input

3
1
2
5

Sample Output

9
97
841

---

　　首先化简一下式子

　　接着我们设

　　为了证明这个展开式一定是成立(n >= 1)，现在进行用归纳法证明

　　①当n = 1时，结论显然成立

　　②当n = k时，假设当n = k - 1时成立，那么有

　　显然等号右边符合这个结构，由此结论成立。

　　同时，我们可以得出x和y的递推关系

　　于是我们可以构造出矩阵，用来快速幂

 

　　现在唯一的问题就是，怎么就这个式子的向下取整的值。正解一个很诡异的解法：

　　用类似的方法可以得出

　　所以有

　　又因为

　　所以会小于1，

　　所以

Code

/**
 * hdu
 * Problem#2256
 * Accepted
 * Time:15ms
 * Memory:2920k
 */
#include<iostream>
#include<cstdio>
#include<cctype>
#include<cstring>
#include<cstdlib>
#include<fstream>
#include<sstream>
#include<algorithm>
#include<map>
#include<set>
#include<queue>
#include<vector>
#include<stack>
using namespace std;
typedef bool boolean;
#define INF 0xfffffff
#define smin(a, b) a = min(a, b)
#define smax(a, b) a = max(a, b)
template<typename T>
inline void readInteger(T& u){
    char x;
    int aFlag = 1;
    while(!isdigit((x = getchar())) && x != '-');
    if(x == '-'){
        x = getchar();
        aFlag = -1;
    }
    for(u = x - '0'; isdigit((x = getchar())); u = (u << 1) + (u << 3) + x - '0');
    ungetc(x, stdin);
    u *= aFlag;
}

#define moder 1024

typedef class Matrix {
    public:
        int* p;
        int lines, cols;
        Matrix():p(NULL), lines(0), cols(0)    {    } 
        Matrix(int lines, int cols):lines(lines), cols(cols) {
            p = new int[(lines * cols)];
        }
        
        int* operator [](int pos) {
            return p + pos * cols;
        }
        
        Matrix operator * (Matrix b) {
            if(cols != b.lines)    return Matrix();
            Matrix res(lines, b.cols);
            for(int i = 0; i < lines; i++) {
                for(int j = 0; j < b.cols; j++) {
                    res[i][j] = 0;
                    for(int k = 0; k < cols; k++) {
                        (res[i][j] += (*this)[i][k] * b[k][j] % moder) %= moder; 
                    } 
                }
            }
            return res;
        }
}Matrix;

int n;
Matrix I2;
Matrix org;
Matrix unit;

template<typename T>
T pow(T a, int pos) {
    if(pos == 0)    return I2;
    if(pos == 1)    return a;
    T temp = pow(a, pos / 2);
    if(pos & 1)    return temp * temp * a;
    return temp * temp; 
}

inline void init_matrix() {
    I2 = Matrix(2, 2);
    I2[0][0] = I2[1][1] = 1;
    I2[1][0] = I2[0][1] = 0;
    org = Matrix(2, 1);
    org[0][0] = 5, org[1][0] = 2;
    unit = Matrix(2, 2);
    unit[0][0] = 5, unit[0][1] = 12;
    unit[1][0] = 2, unit[1][1] = 5;
}

inline void init() {
    readInteger(n);
}

inline void solve() {
    Matrix res = pow(unit, n - 1) * org;
    printf("%d
", (res[0][0] * 2 - 1) % moder);
}

int T;
int main() {
    readInteger(T);
    init_matrix();
    while(T--) {
        init();
        solve();
    }
    return 0;
}