URAL 1814 Continued Fraction 数学 矩阵乘法

题目链接：http://acm.timus.ru/problem.aspx?space=1&num=1814

题目大意：
任何整数的平方根都可以表示为连分数，就是下面的形式。

由于连分数可能是无穷无尽的，因此指定输出到某一位即可。
给$2$个数$x$和$k$，求$\sqrt{x}$用连分数表示到$a_k$的结果。
由于分子分母结果巨大，只需要输出分子和分母分别$mod$ $1000000007$的结果即可。
$2\le{x}\le10^6,1\le{k}\le10^9$

题解：
设${{p_k}\over{q_k}}=[a_1,a_2,a_3,……,a_k]$，我们可以发现：
${p_k= \begin{cases} a_1,&k=1\\ a_1a_2+1,&k=2\\ a_kp_{k-1}+p_{k-2},&k\ge3\\ \end{cases}}$

${q_k= \begin{cases} 1,&k=1\\ a_2,&k=2\\ a_kq_{k-1}+q_{k-2},&k\ge3\\ \end{cases}}$
再将其转化为矩阵乘法的模式：
${\begin{bmatrix}p_k\\p_{k-1}\end{bmatrix}}={\begin{bmatrix}a_k&1\\[1.4ex]1&0\end{bmatrix}}{\begin{bmatrix}p_{k-1}\\p_{k-2}\end{bmatrix}}$

${\begin{bmatrix}q_k\\q_{k-1}\end{bmatrix}}={\begin{bmatrix}a_k&1\\[1.4ex]1&0\end{bmatrix}}{\begin{bmatrix}q_{k-1}\\q_{k-2}\end{bmatrix}}$

也就是说，我们只要能求出所有$a_k$就能算出$p_k$和$q_k$
那么如何找到$a_k$？
最初的式子的形式是这样的：

$$\sqrt{n}=a_0+{{1}\over{x_0}}$$
显然 $a_0$就是 $\sqrt{n}$的整数部分
如果我们想继续求 $a_1$,那么 $a_1$就是 $\lfloor{x_0}\rfloor$
我们可以将每次的式子化成同一个形式：
$${x}={{\sqrt{n}+p_k}\over{q}}$$
其中 $p_0=a_0$， $p_k=a_k*(n-{p_{k-1}}^2)-p_{k-1}$， $q={{n-p_k^2}\over{n-{p_{k-1}}^2}}$（可以证明 $q$为整数）
这样每次通过$a_k=\lfloor{x_{k-1}}\rfloor$就可以求得所有 $a_k$了
最后我们又通过打表发现将 $\sqrt{x}$表示为连分数的话会存在一定的循环节，
将循环节中的结果预处理出来，再做一次快速幂，剩下的部分直接乘上去就可以了。
注意矩阵乘法没有交换律，所以要明确 $a*b$与 $b*a$的区别。

代码：

#include<stdio.h>
#include<string.h>
#include<math.h>
#include<stdlib.h>
#define MOD 1000000007
#define EPS 1e-9
typedef long long ll;
int n;
int k;
struct matrix
{
    ll x[2][2];
    matrix()
    {
        memset(x,0,sizeof(x));
    }
    friend matrix operator*(matrix a,matrix b)
    {
        matrix re;
        for(int i=0;i<2;i++)
        {
            for(int j=0;j<2;j++)
            {
                for(int k=0;k<2;k++)
                {
                    re.x[i][j]=(re.x[i][j]+a.x[i][k]*b.x[k][j])%MOD;
                }
            }
        }
        return re;
    }
};
matrix quickpower(matrix a,ll b)
{
    matrix re;
    re.x[0][0]=1;
    re.x[1][1]=1;
    while(b)
    {
        if(b&1)
        {
            re=re*a;
        }
        a=a*a;
        b>>=1;
    }
    return re;
}
ll a[50005];
int cnt;
ll gcd(ll a,ll b)
{
    return b?gcd(b,a%b):a;
}
void pre()
{
    double k=sqrt((double)n);
    ll q=1,p=(ll)k;
    a[cnt]=p;
    double st=k-p;
    do{
        q=(n-p*p)/q;
        k=(sqrt((double)n)+p)/q;
        a[++cnt]=(ll)k;
        p=a[cnt]*q-p;
    }while(fabs(k-a[cnt]-st)>EPS);
}
int main()
{
    scanf("%d%d",&n,&k);
    pre();
    matrix ans;
    ans.x[0][0]=1;
    ans.x[1][1]=1;
    matrix t;
    t.x[0][1]=1;
    t.x[1][0]=1;
    for(int i=1;i<=cnt;i++)
    {
        t.x[0][0]=a[i];
        ans=ans*t;
    }
    ans=quickpower(ans,k/cnt);
    for(int i=1;i<=k%cnt;i++)
    {
        t.x[0][0]=a[i];
        ans=ans*t;
    }
    t.x[0][0]=a[0];
    ans=t*ans;
    printf("%lld/%lld\n",ans.x[0][0],ans.x[1][0]);
    return 0;
}