特殊的不定方程——佩尔方程

暴力法求特解：由x^2-dy^2=11,得x^2=dy^2+1,即x=根号(dy^2+1),那么令y=1,求得x并检验此时是否满足x^2-dy^2=1;如果满足，那么(x,y)即为所求，否则y+1；

代码如下：

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

例题一：Street Numbers(pku 1320)

题目可转述为，求解两个不相等的正整数n,m(n<m)，使得

1+2+3......+n=(n+1)+.........+m

输入：无

输出：输出前十组满足条件的(n,m)(从小到大)，每组占一行，每个数占十格，向右对齐

输出样例：

         6         8
        35        49

分析：1+2+……+n = (n+1) + ……+m,那么 n(n+1)/2 = (m-n)(m+n+1)/2，化简 (2m+1)^2 - 8*n^2 = 1,令x = 2*m+1,y=n,那么方程可以转化为x^2 - 8*y^2 = 1,看到这里的话 这个是典型的佩尔方程，比较偏，有迭代公式可以利用 令x1 = 3,y1 = 1;

由迭代公式可得：

Xn=3Xn-1+8Yn-1

Yn=Xn-1+3Yn;

#include<cstdio>
#include<iostream>
#include<cstring>
#include<cmath>
using namespace std;
int main()
{
    int x,y;
    int x1,y1;///最小特解
    int px,py;///相当于Xn-1，Yn-1
    int d;
    x1=3;
    y1=1;
    px=3;///刚开始时为x1,y1的值
    py=1;
    d=8;
    for(int i=1;i<=10;i++)
    {
        x=px*x1+d*py*y1;
        y=px*y1+py*x1;
        printf("%10d%10d\n",y,(x-1)/2);
        px=x;
        py=y;
    }
    return 0;
}

例题二：No more tricks, Mr Nanguo(hdu 3292)

现在Sailormoon的女孩想告诉你一个古老的成语故事“滥竽充数”。

在战国时期（公元前475-221），有一个国家叫齐。齐王非常喜欢玉风乐器，他每天下午都有许多音乐家为他演奏。音乐家的数量是一个平方数。每一列都有X的音乐家。

国王不知道乐队的成员Nan Guo。事实上，Nan Guo对竽一无所知。但他不知何故设法把自己作为一个竽的球员坐在后面，假装玩乐器。国王一点也不聪明。齐王死后，国王的儿子继位。新国王不像他父亲，他决定把乐队的音乐家分成若干小部分。他还希望每个部分的数目是平方数。当然，Nan Guo很快就意识到自己的愚蠢会暴露出来，他发现自己再也没有一个乐队可以躲起来了，所以他很快就跑掉了。

他离开后，乐队的数量令人满意。因为现在乐队的数量将分成若干等份，每个部分的数目也都是一个平方数，每一行每一列都有Y个音乐家。

输入：输入数据有多组，每组占一行包含两个正整数N(2<=N<29),K(K<10^9);

输出：输出第K大的满足条件的X值，如果没有这样的X存在，则输出"No answers can meet such conditions".

Sample Input

2 999888
3 1000001
4 8373

Sample Output

7181
600 
No answers can meet such conditions

分析：对题目易得出X^2-NY^2=1;就是求佩尔方程的解

首先，若N是完全平方数，则肯定无解；否则，用暴力法求出方程的解，然后由

结合快速幂就可以快速求出方程的第K大的解

#include<cstdio>
#include<iostream>
#include<cstring>
#include<cmath>
#define MAXN 40
#define M 8191///为啥要对8191取模
using namespace std;
typedef struct
{
    int m[MAXN][MAXN];
}Matrix;
Matrix per,d;
int n;
int x,y,D;
void search()
{
    y=1;
    while(1)
    {
        x=(long long)sqrt(D*y*y+1);
        if(x*x-D*y*y==1)
            break;
        y++;
    }
}
void init()
{
    d.m[0][0]=x%M;
    d.m[0][1]=D*y%M;
    d.m[1][0]=y%M;
    d.m[1][1]=x%M;
    for(int i=0;i<n;i++)
        for(int j=0;j<n;j++)
        per.m[i][j]=(i==j);
}
Matrix multi(Matrix a,Matrix b)///矩阵乘法
{
    Matrix c;
    int i,j,k;
    for(i=0;i<n;i++)
        for(j=0;j<n;j++)
    {
        c.m[i][j]=0;
        for(k=0;k<n;k++)
        {
            c.m[i][j]+=a.m[i][k]*b.m[k][j];
        }
        c.m[i][j]%=M;
    }
    return c;
}
Matrix power(int k)///矩阵快速幂
{
    Matrix p,ans=per;
    p=d;
    while(k)
    {
        if(k&1)
        {
            ans=multi(ans,p);
            k--;
        }
        k>>=1;
        p=multi(p,p);
    }
    return ans;
}
int main()
{
    int K;
    while(scanf("%d%d",&D,&K)!=EOF)
    {
        int ad=sqrt(D+0.0);
        if(ad*ad==D)
        {
            printf("No answers can meet such conditions\n");
            continue;
        }
        search();
        n=2;///矩阵的边界
        init();
        d=power(K-1);///求出前面矩阵的
        x=(d.m[0][0]*x%M+d.m[0][1]*y%M)%M;
        printf("%d\n",x);
    }
    return 0;
}