HDU 2256 Problem of Precision （矩阵快速幂）（推算）

Problem of Precision

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

Problem Description

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

【分析】

这个题目算是矩阵快速幂的比较难推的一个题目。题目要求 （sqrt(2)+sqrt(3)）的 2^n并%1024，要求出值来并不难，构造矩阵即可，但是要mod1024就有问题了，小数不能直接mod，但是如果你取整之后再mod，结果绝逼出问题，因为浮点数的精度问题。

所以从斌牛的博客上看到如此推算，推算第一块不难，而且很容易求出Xn 和 Yn，但是问题又出来了，要是求出来后，直接用（int）（Xn+Yn*sqrt（6））%1024，又会出问题，还是浮点数取整问题，我一开始就这么算的，导致结果奇葩。看来在mod的时候有浮点数要格外注意，直接处理的话，不管怎么取整，都会出问题。

所以分割线下面的推算就避开了这个问题，这个确实好难想到，通过变换一下，得到最终的结果必定是2Xn-（0.101...）^n,因为最终mod是用不大于浮点数的最大整数在mod，所以最终结果就是2Xn-1.第二条确实好难想到！

题解转载于 http://www.cnblogs.com/kkrisen/p/3437710.html；

#include <iostream>
#include <cstring>
#include <cstdio>
#include <algorithm>
#include <cmath>
#include <string>
#include <stack>
#include <queue>
#include <vector>
#define inf 0x3f3f3f3f
#define met(a,b) memset(a,b,sizeof a)
#define pb push_back
using namespace std;
typedef long long ll;
const ll N = 2;

ll f1,f2,k;
ll mod = 1024;
ll n;

struct Fast_Matrax {
    ll a[N][N];
    Fast_Matrax() {
        memset(a,0,sizeof(a));
    }
    void init() {
        for(int i=0; i<N; i++)
            for(int j=0; j<N; j++)
                a[i][j]=(i==j);
    }
    Fast_Matrax operator * (const Fast_Matrax &B)const {
        Fast_Matrax C;
        for(int i=0; i<N; i++)
            for(int k=0; k<N; k++)
                for(int j=0; j<N; j++)
                    C.a[i][j]=(C.a[i][j]+1LL*a[i][k]*B.a[k][j]%mod+mod)%mod;
        return C;
    }
    Fast_Matrax operator ^ (const ll &t)const {
        Fast_Matrax A=(*this),res;
        res.init();
        ll p=t;
        while(p) {
            if(p&1)res=res*A;
            A=A*A;
            p>>=1;
        }
        return res;
    }
} ans,tmp,x;
int main() {
    x.a[0][0]=5;x.a[1][0]=2;
    int T;
    scanf("%d",&T);
    while(T--){
        scanf("%lld",&n);
        if(n<=1){
            puts("9");
        }
        else {
            tmp.a[0][0]=5;tmp.a[0][1]=12;
            tmp.a[1][0]=2;tmp.a[1][1]=5;
            ans=(tmp^(n-1))*x;
            printf("%lld
",(2*ans.a[0][0]-1)%mod);
        }
    }
    return 0;
}