NOD 1147 还原连分数

0
0
0
1. 云栖社区>
2. 博客>
3. 正文

## NOD 1147 还原连分数

prime7 2013-07-29 19:07:00 浏览699

```#include <iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#define M 1000000007
using namespace std;
const double eps=1e-8;
const int MAX=2;
typedef struct
{
long long m[MAX][MAX];
} Matrix;
Matrix I= {1,0,
0,1
};
Matrix yl= {0,1,
1,0
};
Matrix data,ans,P,Q;
Matrix matrixmul(Matrix a,Matrix b) //矩阵乘法
{
int i,j,k;
Matrix c;
for (i = 0 ; i < MAX; i++)
for (j = 0; j < MAX; j++)
{
c.m[i][j] = 0;
for (k=0; k<MAX; k++)
c.m[i][j]+=((a.m[i][k]%(M))*(b.m[k][j]%(M)))%(M);
c.m[i][j] %=M;
}
return c;
}
Matrix quickpow(long long n)
{
Matrix m = P, b = I;
while (n >= 1)
{
if (n & 1)
b = matrixmul(b,m);
n = n >> 1;
m = matrixmul(m,m);
}
return b;
}
long long gcd(long long a,long long b)
{
return b?gcd(b,a%b):a;
}
int getxun(long long n,long long *a)
{
double k=sqrt(n*1.0);
int cnt=0;
long long q=1,p=(long long)k,tmp=1;
double first=k-p;
while(1)
{
tmp=q;
q=n-p*p;
long long G=gcd(tmp,q);
q/=G;
tmp/=G;
k=(sqrt(1.0*n)+p)*tmp/q;
a[cnt++]=(long long)k;
p=abs(p-a[cnt-1]*q);
if(fabs(k-a[cnt-1]-first)<eps) break;
}
return cnt;
}
int main()
{
int n,k,f;
long long s[10000];
long long p,q;
while(~scanf("%d%d",&n,&k))
{
P=I,ans=I,Q=I;
f=getxun(n,s);
int m=k/f,yu=k%f;
if(m)
{
for(int i=f-1; i>=0; i--)
{
data=yl;
data.m[0][0]=s[i];
P=matrixmul(P,data);
}
Q=quickpow(m);
}
for(int i=yu-1; i>=0; i--)
{
data=yl;
data.m[0][0]=s[i];
ans=matrixmul(ans,data);
}
ans=matrixmul(ans,Q);
p=ans.m[0][0];
q=ans.m[0][1];
long long ans1=(((long long)sqrt(n)*p)%M+q)%M;
cout<<ans1<<"/"<<p<<endl;
}
return 0;
}
```

prime7
+ 关注