【codejam2008_Round1A_C】Numbers（矩阵快速幂）

题目大意：
计算$(3+\sqrt 5)^n$整数部分末三位。

题解：
题解原网址https://code.google.com/codejam/contest/32016/dashboard#s=a&a=2

设$\alpha=3+\sqrt 5$，$\beta=3-\sqrt 5$，$X_n=\alpha^n+\beta^n$
可证明$X_n$为整数（二项式定理）

$$X_n=\alpha^n+\beta^n=\sum^n_{i=0}\left(^n_i\right)3^{n-i}{(-\sqrt 5)}^i+\sum^n_{i=0}\left(^n_i\right)3^{n-i}\sqrt 5^i=2\left(\sum^{\lfloor n/2\rfloor}_{i=0}\left(^n_{2i}\right)3^{n-2i}5^i\right)$$
$∵0<\beta^n<1$
∴ $\alpha^n$的整数部分一定等于 $X^n-1$

Solution1：（共轭）

设$\alpha^n=a_n+b_n\sqrt 5$，$\beta^n=a_n-b_n\sqrt 5$
$∴X_n=2a_n$
又$∵\alpha^{n+1}=(3+\sqrt 5)(a_n+b_n\sqrt 5)=(3a_n+5b_n)+(a_n+3b_n)\sqrt 5$
$∴a_{n+1}=3a_n+5b_n$，$b_{n+1}=a_n+3b_n$

$$\left( \begin{matrix} a_n\\ b_n \end{matrix}\right)= A\left( \begin{matrix} a_{n-1}\\ b_{n-1} \end{matrix}\right)= A^n\left( \begin{matrix} a_0\\ b_0 \end{matrix}\right) , A= \left( \begin{matrix} 3&5\\1&3 \end{matrix} \right)$$
初始状态： $\alpha^0=1$，
$$\left( \begin{matrix} a_0\\ b_0 \end{matrix}\right)= \left( \begin{matrix} 1\\ 0 \end{matrix}\right)$$
可使用矩阵快速幂，时间复杂度 $O(log n)$

方法2：

$∵\alpha+\beta=6$，$\alpha\beta=4$
所以α和β是$x_2-6x+4=0$的两个根
$∴\alpha^2=6\alpha-4$，$\beta^2=6\beta-4$
$∴\alpha^{n+2}=6\alpha^{n+1}-4\alpha^n$，$\beta^{n+2}=6\beta^{n+1}-4\beta^n$
两式相加，得：

$$X_{n+2}=6X_{n+1}-4X_n$$
所以可使用矩阵乘法：
$$\left( \begin{matrix} X_{n+1}\\ X_n \end{matrix}\right)= B\left( \begin{matrix} X_n\\ X_{n-1} \end{matrix}\right)= B^n\left( \begin{matrix} X_1\\ X_0 \end{matrix}\right) , B= \left( \begin{matrix} 6&{-4}\\1&0 \end{matrix} \right)$$

方法3：（方法2存在循环，可优化）

我们只需要最后结果的最后3位，经过法2的多次计算后发现，结果从$n=4$至$n=103$是一个长度为100的循环，然后继续下去。
如果n很大，输出$(n-3)mod 100+3$

代码：

#include<cstdio>
#include<vector>
using namespace std;
typedef vector<int> Vector;
typedef vector<Vector> Matrix;
Matrix operator * (Matrix a,Matrix b)
{
    int as1=a.size(),as2=a[0].size(),bs1=b.size(),bs2=b[0].size();
    if(as2!=bs1)
        return Matrix(0,Vector(0));
    Matrix ret(as2,Vector(bs1,0));
    for(int i=0;i<as1;i++)
        for(int j=0;j<bs2;j++)
            for(int k=0;k<as2;k++)
                ret[i][j]=(ret[i][j]+(a[i][k]*b[k][j])%1000)%1000;
    return ret;
}
Matrix pow(Matrix a,int b)
{
    Matrix c(a.size(),Vector(a.size(),0));
    for(int i=0;i<(int)a.size();i++)
        c[i][i]=1;      
    for(;b;b>>=1,a=a*a)
        if(b&1)
            c=c*a;
    return c;
}
namespace SolutionA
{
    int solve(int n)
    {
        Matrix a(2,Vector(2,0));
        a[0][0]=3;a[0][1]=5;
        a[1][0]=1;a[1][1]=3;
        a=pow(a,n);
        return (a[0][0]*2+999)%1000;
    }
}
namespace SolutionB
{
    int solve(int n)
    {
        Matrix b(2,Vector(2));
        b[0][0]=6;b[0][1]=-4;
        b[1][0]=1;b[1][1]=0;
        b=pow(b,n);
        int ret=b[1][0]*6+b[1][1]*2+999;
        while(ret<0)
            ret+=1000;
        return ret%1000;
    }
}
namespace SolutionC
{
    int solve(int n)
    {
        if(n>103)
            n=(n-3)%100+3;
        return SolutionB::solve(n);
    }
}
int main()
{
    int T,n;
    scanf("%d",&T);
    for(int i=1;i<=T;i++)
    {
        scanf("%d",&n);
        //printf("Case #%d: %03d\n",i,SolutionA::solve(n));
        //printf("Case #%d: %03d\n",i,SolutionB::solve(n));
        printf("Case #%d: %03d\n",i,SolutionC::solve(n));
    }
    return 0;
}