HDU3932 再谈佩尔方程

来源：互联网发布：若宇网络编辑：程序博客网时间：2024/06/05 00:54

首先基本的佩尔方程是啥？

一个不定方程形如：x^2-dy^2=1（d>1且d不是完全平方数），为什么不能是完全平方数？

反证法：假设是完全平方数则上式有（x+sqrt(d)*y）*(x-sqrt(d)*y)=1

那么若为1 则有：

x+sqrt(d)*y=1; x+sqrt(y)=-1;

x-sqrt(d)*y=1; x-sqrt(y)=-1;这两种情况去解这个方程显然是无法满足上面说的条件的，因此无解。

那么来看满足条件的d的情况，咱们不妨设(x1,y1)是满足该方程的最小解,则有：

xn^2-yn^2*d=1;---------①

(xn-1)^2-(yn-1)^2*d=1; ------②

x1^2-y1^2*d=1;---------③

把②式做个等量代换有：

（xn-1）^2*(x1^2-y1^2*d)-d*(yn-1)^2*(x1^2-d*y1^2)=1;

展开有：

(xn-1*x1)^2-d*y1^2*(xn-1)^2-d*(yn-1)^2*x1^2+d^2*y1^2*(yn-1)^2=1;

转化成完全平方的形式有：

(xn-1*x1+dyn-1*y1)^2-2*d*(yn-1*y1*xn-1*x1)-d*(y1^2*(xn-1)^2+x1^2*(yn-1)^2)=1;

提取公因子d有：

(xn-1*x1+dyn-1*y1)^2-d(xn-1*y1+yn-1*x1)^2=1

不难发现我们可以得到一个递推的公式：

xn=xn-1*x1+dyn-1*y1

yn=xn-1*y1+yn-1*x1

这样一来，当知道最小特解的时候，我们就能通过递推式构造矩阵来求解满足条件的第k组解。

也即：

|x1 dy1 |^(k-1) | x1 | | xn |

* =

|y1 x1 | | y1 | | yn |

对于本题就是求x^2-N*y^2=1的第k组解。

所以咱们得先暴力的找出x1 y1这组特解，然后再来构造矩阵，当然如果N是个完全平方数，就直接无解了咯，到此这题便很好解决了！

#include<cstdio>#include<cstring>#include<cmath>#include<iostream>#include<algorithm>using namespace std;const long long mod=8191;long long n,k;struct juzhen{long long m[2][2];};juzhen mut(juzhen a,juzhen b){    juzhen ans;    for(int i=0;i<2;i++)    {        for(int j=0;j<2;j++)        {            ans.m[i][j]=0;            for(int k=0;k<2;k++)            {                ans.m[i][j]+=(a.m[i][k]*b.m[k][j])%mod;                ans.m[i][j]%=mod;            }        }    }    return ans;}juzhen power(juzhen a,long long p){    juzhen ans;    memset(ans.m,0,sizeof(ans.m));    for(int i=0;i<2;i++)    {        ans.m[i][i]=1;    }    while(p)    {        if(p&1)            ans=mut(ans,a);        a=mut(a,a);        p>>=1;    }    return ans;}long long x,y;void sousuo(){    y=1;    while(1)    {        x=sqrt(y*y*n+1);        if(x*x-n*y*y==1)            break;        y++;    }}//暴力找第一组特解的过程int main(){    while(cin>>n>>k)    {        long long t=sqrt(n+0.0);        if(t*t==n)        {            cout<<"No answers can meet such conditions"<<endl;        }        else        {            juzhen ans;            sousuo();            ans.m[0][0]=ans.m[1][1]=x;            ans.m[0][1]=n*y;            ans.m[1][0]=y;            ans=power(ans,k-1);            cout<<(ans.m[0][0]*x+ans.m[0][1]*y)%mod<<endl;        }    }    return 0;}

0 0