佩尔方程

讲解佩尔方程几个不错的博客地址：/article/2371901.html

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 8191
using 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,还是太菜！！！自己在慢慢想哈吧