HDU 3292 No more tricks, Mr Nanguo(佩尔方程)

最新推荐文章于 2019-11-15 00:02:04 发布

w20810

最新推荐文章于 2019-11-15 00:02:04 发布

阅读量482

点赞数

分类专栏： ACM-数学文章标签：佩尔方程

本文链接：https://blog.youkuaiyun.com/w20810/article/details/43563801

版权

ACM-数学专栏收录该内容

73 篇文章

订阅专栏

题意：输入N，K，求方程x^2+N*(y^2)=1的第K大解，结果MOD8191

分析：利用公式，然后直接矩阵快速幂。

http://blog.youkuaiyun.com/w20810/article/details/43527357

代码：

#include <iostream>
#include <cstdio>
#include <cmath>

using namespace std;

#define M 8191

typedef long long int64;

int D,K;

struct Matrax
{
	int mat[4][4];
}U,F;

void Search(int d,int &x,int &y)
{
	y=1;
	while(1)
	{
		x=sqrt(1.0+d*y*y);
		if(x*x-d*y*y==1)
			break;
		y++;
	}
}

Matrax multi(Matrax a,Matrax b)
{
	Matrax ans;
	for(int i=0;i<2;i++)
		for(int j=0;j<2;j++)
		{
			ans.mat[i][j]=0;
			for(int k=0;k<2;k++)
				ans.mat[i][j]+=a.mat[i][k]*b.mat[k][j];
			ans.mat[i][j]%=M;
		}
	return ans;
}

Matrax mkpow(int n)
{
	Matrax ans=U,p=F;
	while(n)
	{
		if(n&1)
			ans=multi(ans,p);
		n>>=1;
		p=multi(p,p);
	}
	return ans;
}

int main()
{
	int x1,y1;
	while(scanf("%d%d",&D,&K)!=EOF)
	{
		int t=sqrt(1.0*D);
		if(D==t*t)
		{
			printf("No answers can meet such conditions\n");
			continue;
		}
		Search(D,x1,y1);
		
		
		U.mat[0][0]=1;
		U.mat[0][1]=0;
		U.mat[1][0]=0;
		U.mat[1][1]=1;
		
		F.mat[0][0]=x1%M;
		F.mat[0][1]=D*y1%M;
		F.mat[1][0]=y1%M;
		F.mat[1][1]=x1%M;
		
		Matrax ans=mkpow(K-1);
		int X=ans.mat[0][0]*x1+ans.mat[0][1]*y1;
		//int Y=ans.mat[1][0]*x1+ans.mat[1][1]*y1;
		printf("%d\n",X%M);
	}
	return 0;
}