URAL 1814 Continued Fraction 数学 矩阵乘法

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

题目大意：

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



由于连分数可能是无穷无尽的，因此指定输出到某一位即可。

给2个数x和k，求x√用连分数表示到ak的结果。

由于分子分母结果巨大，只需要输出分子和分母分别mod 1000000007的结果即可。

2≤x≤106,1≤k≤109

题解：

设pkqk=[a1,a2,a3,……,ak]，我们可以发现：

pk=⎧⎩⎨a1,a1a2+1,akpk−1+pk−2,k=1k=2k≥3

qk=⎧⎩⎨1,a2,akqk−1+qk−2,k=1k=2k≥3

再将其转化为矩阵乘法的模式：

[pkpk−1]=[ak110][pk−1pk−2]

[qkqk−1]=[ak110][qk−1qk−2]

也就是说，我们只要能求出所有ak就能算出pk和qk

那么如何找到ak？

最初的式子的形式是这样的：

n√=a0+1x0

显然a0就是n√的整数部分

如果我们想继续求a1,那么a1就是⌊x0⌋

我们可以将每次的式子化成同一个形式：

x=n√+pkq

其中p0=a0，pk=ak∗(n−pk−12)−pk−1，q=n−p2kn−pk−12（可以证明q为整数）

这样每次通过ak=⌊xk−1⌋就可以求得所有ak了

最后我们又通过打表发现将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;
}