11
23
2016
0

51nod 1147 连分数

即ural1814

这里说明了一个无理数的连分数项是会重复的当且仅当它是一个一元二次方程的解。

还说明了$\sqrt{n}$的连分数项的循环节规模大约在$O(\ln(n)\sqrt{n})$虽然并没有说明为什么

还有就是循环节的最后一项一定是$\sqrt{n}$的整数部分的2倍。虽然也没有说明为什么

那么就好做了。

我们可以写一个struct来表示$\frac {k\sqrt{n}+a} {b}$形式的数,然后爆下去就可以找出循环节就行了。

考虑从下往上倒着做的时候,假如当前数是x,当前连分数项为a,那么x会变成$\frac {1+ax} {x}$原谅我语文水平有限。。。

那么假如x是$\frac {u} {v}$的形式,那么$\frac {1+ax} {x}=\frac {u*a+v} {u}$

这里要注意,假如说$u$和$v$互质的话,$u*a+v$和$u$必定互质。而且初始$u=1$所以一定不会出现约分的情况。直接爆爆爆就行了。

可以注意到每次做一个a相当于是乘上一个2*2的矩阵。

所以后面m%cnt的部分直接做,中间的部分快速幂就行了喽。

总复杂度$O(\ln(n)\sqrt{n}+log_{2}(m/(\ln(n)\sqrt{n})))$什么鬼

萎住了不知道什么时候能放代码

#include<cstdio>
#include<cmath>
#define mo 1000000007
using namespace std;
int n,m,cnt,k;
long double X;
int a[30005];
inline int gcd(int x,int y){
	if(y==0)return x;
	return x==0?y:gcd(y%x,x);
}
inline int abs(int k){return k>0?k:-k;}
struct ele{
	int fz,fm,fzn;
	inline int f(){return int((fz+fzn*X)/fm);}
	ele operator +(const int k)const{
		ele ans;
		ans.fz=k*fm+fz;ans.fm=fm;ans.fzn=fzn;
		int tmp=gcd(abs(ans.fz),abs(ans.fm));
		tmp=gcd(tmp,abs(ans.fzn));
		ans.fz/=tmp;ans.fm/=tmp;ans.fzn/=tmp;
		return ans;
	}
	ele operator -(const int k)const{
		ele ans;
		ans.fz=fz-k*fm;ans.fm=fm;ans.fzn=fzn;
		int tmp=gcd(abs(ans.fz),abs(ans.fm));
		tmp=gcd(tmp,abs(ans.fzn));
		ans.fz/=tmp;ans.fm/=tmp;ans.fzn/=tmp;
		return ans;
	}
	void ud(){
		int t=fm;
		fm=fzn*fzn*n-fz*fz;
		fzn=t*fzn;fz=t*-fz;
		if(fm<0)fm*=-1,fzn*=-1,fz*=-1;
		int tmp=gcd(abs(fz),abs(fm));
		tmp=gcd(tmp,abs(fzn));
		fz/=tmp;fm/=tmp;fzn/=tmp;
	}
}x;
struct matrix{
	int x[2][2];
	matrix operator +(const matrix k)const{
		matrix ans;
		for(int i=0;i<2;i++)
			for(int j=0;j<2;j++)ans.x[i][j]=(x[i][j]+k.x[i][j])%mo;
		return ans;
	}
	matrix operator *(const matrix k)const{
		matrix ans;
		for(int i=0;i<2;i++)
			for(int j=0;j<2;j++)ans.x[i][j]=(1LL*x[i][0]*k.x[0][j]+1LL*x[i][1]*k.x[1][j])%mo;
		return ans;		
	}
	matrix pow(long long k){
		matrix ans,X;
		for(int i=0;i<2;i++)
			for(int j=0;j<2;j++)X.x[i][j]=x[i][j];
		ans.x[0][0]=ans.x[1][1]=1;ans.x[0][1]=ans.x[1][0]=0;
		for(;k;k>>=1,X=X*X)if(k&1)ans=ans*X;
		return ans;
	}
}p,now,i;
int ans[2];
void updata(matrix p){
	int tmp[2];
	tmp[0]=ans[0];tmp[1]=ans[1];
	ans[0]=(1LL*tmp[0]*p.x[0][0]+1LL*tmp[1]*p.x[1][0])%mo;
	ans[1]=(1LL*tmp[0]*p.x[0][1]+1LL*tmp[1]*p.x[1][1])%mo;
}
int main(){
	scanf("%d%d",&n,&m);
	a[0]=int(sqrt(n));
	X=sqrt(n);
	x.fm=1;x.fzn=1;
	x=x-x.f();x.ud();
	while(true){
		a[++cnt]=x.f();
		if(a[cnt]==2*a[0])break;
		x=x-x.f();x.ud();
	}
	m--;ans[0]=a[m%cnt+1];ans[1]=1;
	p.x[0][1]=1;p.x[1][0]=1;
	for(int i=m%cnt;i;i--){
		p.x[0][0]=a[i];
		updata(p);
	}
	now.x[1][1]=1;now.x[0][0]=1;
	for(int i=cnt;i;i--){
		p.x[0][0]=a[i];
		now=now*p;
	}
	k=m/cnt;
	for(;k;k>>=1,now=now*now)if(k&1)updata(now);
	ans[1]=(ans[1]+1LL*ans[0]*a[0])%mo;
	printf("%d/%d",ans[1],ans[0]);
}

 

Category: 51nod | Tags: 数论 矩阵乘法 51nod | Read Count: 380

登录 *


loading captcha image...
(输入验证码)
or Ctrl+Enter

Host by is-Programmer.com | Power by Chito 1.3.3 beta | Theme: Aeros 2.0 by TheBuckmaker.com