HDU3932 再谈佩尔方程

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

一个不定方程形如：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;
}