佩尔方程

讲解佩尔方程几个不错的博客地址：http://blog.csdn.net/acdreamers/article/details/8529686

http://hi.baidu.com/aekdycoin/item/a45f7c37850e5b9db80c03d1

http://hi.baidu.com/shouzhewei/item/5ff25ee1624c3419585dd832

ax^2 - by^2 = c       Pell 方程一般解法

以下涉及到定理的证明的部分全部略过.

开题自然少不了介绍,以上的公式就是Pell方程的一般形态.
显然如果告诉你a,b,c,一开始想到的只可能是暴力,可是接下来介绍的纯数学的方法可以很快速的求解几乎大部分解.
1.首先构造一个系数矩阵,显然为了构造这个矩阵，我们需要先得到下面方程的一个最小特解(x,y>0)

至于如何得到，可以使用暴力(当某些情况下暴力几乎求不到最小解)或者使用连分数的方法来求
假设我们得到了以上方程的特解: x0 y0 (x0,y0>0,并是最小的满足条件的解)
2.继续求

的一个最小特解.假设是x1,y1(x1,y1>0)
3.
假设你要求第k个解,那么有

例子:
1.求 x^2 - 3y^2 = 1的解
由于这里a=1,b=3,而c=1,所以我们可以知道x0=x1,y0=y1;
不难解得一个最小特解(2,1)，于是有

假设现在要知道第2个解
那么套用上面的公式得到
x2=7
y2=4
即
49-48=1
其他解类似.

2.求 x^2 - 3y^2 = 13的解
显然x^2 - 3y^2 =1的最小特解在上面已经求出来了
即
x0=2;
y0=1;
现在我们需要知道的是x^2 - 3y^2 = 13的最小特解，显然应该是
(4,1)
于是如果继续套用上面的解，可以得到:

那么得到
(x1,y1) = (4,1)
(x2,y2) = (11,6)

在这儿，特殊的佩尔方程x^2-n*y^2=1的求解

最小特解是(x1,y1),那么就有迭代公式

    xn=x（n-1）*x1+d*y(n-1)*y1;

    yn=x1*y(n-1)+x(n-1)*y1

在这儿 所有的解为就如上面介绍的一样

关于连分数求解佩尔方程

Pell方程连分数定理：D为非平方数，记sqrt(D)的连分数形式为[非循环部分，循环部分] 即[非循环部分，b1,b2,b3…bm]. 

记渐进分数p/q=[非循环部分，b1,b2,b3…bm-1].

则最小整数解为：

m为偶数：x=p,     y=q

m为奇数：x=p*p+D*q*q,     y=2*p*q

1.第一类问题为给D求一对最小解（x,y）

问题1：如何将一个无理数变换为连分数形式

由欧拉前辈的结论：任何w(n)可表示成 g(n)+sqrt(n)/h(n) 和递推关系w(n)=1/(w(n-1)-a(n)) ,再利用数学归纳导出结论：

g(-1)=0   ,   h(-1)=1

g(n)=-g(n-1)+a(n)*h(n-1) ,h(n)= (n-g(n)*g(n)) /h(n-1)

由结论我们可知a(n)=ceil(w(n-1))=> (g(n-1)+a(0))/h(n-1) 成功解决误差问题

问题2：如何判断循环部分（不像分数除法可以判重余数，因为剩余的是小数部分）

只能每算出一个渐进分数p/q 都需要验证一下时候满足x*x-D*y*y 等于1或﹣1.

先给出算法的伪代码:

第一类问题：最小整数解我们有枚举容易得知，问我们第n个整数解

要求精度的话因为是线性齐次递推关系在n比较大的情况下由矩阵乘法完成

参考资料：2009集训队论文，金斌

对于求解最小特解，直接暴力

void find(){    y=1;    while(1)    {        x=sqrt((D*y*y+1)*1.0);        if(x*x-D*y*y==1)            break;        y++;    }}

对于连分数求解佩尔方程的最小特解

#include<cstdio>#include<iostream>#include<cmath>#include<cstdlib>#include<cstring>using namespace std;void Pell(int n){    __int64 p1=1,p2=0,p;    __int64 q1=0,q2=1,q;    __int64 a1=(int)sqrt(n*1.0),a0=a1,a2;    __int64 g1=0,g;    __int64 h1=1,h;    if(a0*a0==n)    {        printf("No solution!\n");        return;    }    while(1)    {        g=-g1+a1*h1;g1=g;        h=(n-g*g)/h1;h1=h;        a2=(g+a0)/h;        p=a1*p1+p2;p2=p1;p1=p;        q=a1*q1+q2;q2=q1;q1=q;        a1=a2;        if(p*p-n*q*q==1) break;    }    printf("%I64d %I64d\n",p,q);}int main(){    int n;    while(scanf("%d",&n)!=EOF)    {        Pell(n);    }    return 0;}

POJ1320
题目连接：http://poj.org/problem?id=1320

/****************** 题意：求出前10个满足佩尔方程的x，y* 思路：直接枚举就可以******************/#include<iostream>#include<algorithm>#include<cmath>#include<cstdio>#include<ctime>#include<climits>using namespace std;int main(){    int k=10;    int x,y;    int vx=3,vy=1;    int vx1=3,vy1=1;    int n=8;    while(k--)    {        x=vx1*vx+n*vy1*vy;        y=vx*vy1+vy*vx1;        vx1=x;        vy1=y;        printf("%10d%10d\n",y,(x-1)/2);    }    return 520;}

HDU3292

题目连接：http://acm.hdu.edu.cn/showproblem.php?pid=3292

/************************ 题意：给定一个n和k，求解第k个满足x^2-n*y^2=1的解x和y* 思路：通过求解最小特解x1和y1，再通过矩阵来求解第k个就可以了，在这儿使用暴力求解最小特解，然后通过矩阵*乘法和幂来求解就可以了***********************/#include<iostream>#include<cstdio>#include<cstring>#include<algorithm>#include<cstring>#include<cstring>#include<cmath>#define maxn 8191using namespace std;struct matrix{    int m[2][2];} per,d;int x,y,D;void find(){    y=1;    while(1)    {        x=sqrt((D*y*y+1)*1.0);        if(x*x-D*y*y==1)            break;        y++;    }}void init(){    d.m[0][0]=x%maxn;    d.m[0][1]=D*y%maxn;    d.m[1][0]=y%maxn;    d.m[1][1]=x%maxn;    for(int i=0; i<2; i++)    {        for(int j=0; j<2; j++)            per.m[i][j]=(i==j);    }}matrix multi(matrix a,matrix b){    matrix c;    for(int i=0; i<2; i++)    {        for(int j=0; j<2; j++)        {            c.m[i][j]=0;            for(int k=0; k<2; k++)            {                c.m[i][j]+=a.m[i][k]*b.m[k][j];            }            c.m[i][j]%=maxn;        }    }    return c;}matrix power(int k){    matrix p=d;    matrix ans=per;    while(k)    {        if(k&1)        {            ans=multi(ans,p);            k--;        }        k/=2;        p=multi(p,p);    }    return ans;}int main(){    int K;    while(cin>>D>>K)    {        int t=sqrt(D*1.0);        if(t*t==D)            cout<<"No answers can meet such conditions"<<endl;        else        {            find();            init();            d=power(K-1);            x=(d.m[0][0]*x%maxn+d.m[0][1]*y%maxn)%maxn;            printf("%d\n",x);        }    }    return 0;}

POJ2427

题目连接：http://poj.org/problem?id=2427

/**************** 因为是高精度的问题，所以用java的BigInteger类来求解****************/import java.math.BigInteger;  import java.util.Scanner;  public class Main{      public static void solve(int n)     {          BigInteger N, p1, p2, q1, q2, a0, a1, a2, g1, g2, h1, h2,p,q;          g1 = q2 = p1 = BigInteger.ZERO;          h1 = q1 = p2 = BigInteger.ONE;         N = BigInteger.valueOf(n);           a0 = a1 = BigInteger.valueOf((int)Math.sqrt(1.0*n));        if (a1.pow(2).compareTo(BigInteger.valueOf(n)) == 0)        {        System.out.println("No solution!");        return;        }        while (true)         {              g2 = a1.multiply(h1).subtract(g1); g1=g2;                h2 = N.subtract(g2.pow(2)).divide(h1); h1=h2;            a2 = g2.add(a0).divide(h2);                      p = a1.multiply(p2).add(p1); p1=p2;p2=p;                      q = a1.multiply(q2).add(q1);q1=q2;q2=q;a1=a2;                      if (p.pow(2).subtract(N.multiply(q.pow(2))).compareTo(BigInteger.ONE) == 0) break;        }        System.out.println(p+" "+q);    }      public static void main(String[] args)     {          Scanner cin = new Scanner(System.in);          while(cin.hasNextInt())        {            solve(cin.nextInt());           }    }  }

HDU3005

题目连接：http://acm.hdu.edu.cn/showproblem.php?pid=3005

这个题目目前还没有AC,还是太菜！！！自己在慢慢想哈吧